Transcript Visualizza/apri - POLITesi
POLITECNICO DI MILANO
Scuola di Ingegneria Industriale Corso di Laurea in Ingegneria Energetica Dipartimento di Energia
MODELLAZIONE DEL PROCESSO DI COMBUSTIONE IN MOTORI DIESEL MEDIANTE MODELLI DI COMBUSTIONE BASATI SU CINETICA CHIMICA TABULATA Relatore
: Prof. Tommaso LUCCHINI
Tesi di Laurea di:
Maurizio MASTRAPASQUA Matr. 820200 Anno Accademico 2015 - 2016
Indice
Elenco delle Figure Elenco delle tabelle Sommario Introduzione 1 La combustione nei motori Diesel
1.1
Ritardo d’accensione . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2
1.3
Combustione diffusiva . . . . . . . . . . . . . . . . . . . . . . . . .
Modalità avanzate di combustione . . . . . . . . . . . . . . . . . . .
1.3.1
1.3.2
1.3.3
HCCI PCCI RCCI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4
1.5
Cinetica di reazione . . . . . . . . . . . . . . . . . . . . . . . . . . .
Simulazione della combustione Diesel . . . . . . . . . . . . . . . . .
1.5.1
1.5.2
1.5.3
Equazioni di conservazione . . . . . . . . . . . . . . . . . . .
Turbolenza . . . . . . . . . . . . . . . . . . . . . . . . . . .
Determinazione del reaction rate . . . . . . . . . . . . . . .
3
2 Tabulazione della cinetica chimica
2.1
Inquadramento del problema . . . . . . . . . . . . . . . . . . . . . .
2.2
2.3
2.4
ISAT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ILDM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
FGM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3 Generazione della tabella
3.1
3.2
3.3
Modello PSR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Progress variable . . . . . . . . . . . . . . . . . . . . . . . . . . . .
PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
iii
vii
xi
xv
1
iv INDICE
3.4
3.5
Mixing Line . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Costruzione della tabella . . . . . . . . . . . . . . . . . . . . . . . .
3.5.1
Correzione della tabulazione . . . . . . . . . . . . . . . . . .
3.5.2
3.5.3
Virtual Species . . . . . . . . . . . . . . . . . . . . . . . . .
Convergenza automatica . . . . . . . . . . . . . . . . . . . .
4 Il solver FlameletDieselFoam
4.1
Equazioni del modello . . . . . . . . . . . . . . . . . . . . . . . . .
4.2
Flamelet combustion library . . . . . . . . . . . . . . . . . . . . . .
4.2.1
Tecniche di interpolazione . . . . . . . . . . . . . . . . . . .
41
5 Verifica della consistenza del modello
5.1
Prove a pressione costante . . . . . . . . . . . . . . . . . . . . . . .
5.1.1
5.1.2
Effetto della Z . . . . . . . . . . . . . . . . . . . . . . . . . .
Effetto della T . . . . . . . . . . . . . . . . . . . . . . . . .
5.2
5.3
5.1.3
Effetto della p . . . . . . . . . . . . . . . . . . . . . . . . . .
Prove a volume costante . . . . . . . . . . . . . . . . . . . . . . . .
5.2.1
5.2.2
Effetto della Z . . . . . . . . . . . . . . . . . . . . . . . . . .
Effetto della T . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.3
Effetto della p . . . . . . . . . . . . . . . . . . . . . . . . . .
Calcolo degli ignition delay con l’n-dodecano . . . . . . . . . . . . .
47
6 Validazione
6.1
Simulazione della combustione
6.1.1
Condizioni termodinamiche . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . .
6.1.2
6.1.3
6.1.4
Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Modelli utilizzati . . . . . . . . . . . . . . . . . . . . . . . .
Analisi dei risultati . . . . . . . . . . . . . . . . . . . . . . .
6.2
Campagna sperimentale Spray A 6.2.1
Condizioni termodinamiche . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . .
6.2.2
6.2.3
Mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Iniezione di combustibile . . . . . . . . . . . . . . . . . . . .
6.2.4
6.2.5
Modelli utilizzati . . . . . . . . . . . . . . . . . . . . . . . .
Analisi dei risultati . . . . . . . . . . . . . . . . . . . . . . .
61
Conclusioni
85
Acronimi
1 Homogeneous Charge Compression Ignition
87
INDICE v Bibliografia
Riferimenti citati nel testo . . . . . . . . . . . . . . . . . . . . . . . . . .
Pubblicazioni e Manuali . . . . . . . . . . . . . . . . . . . . . . . .
Materiale Online . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Ulteriore materiale consultato . . . . . . . . . . . . . . . . . . . . . . . .
Pubblicazioni e Manuali . . . . . . . . . . . . . . . . . . . . . . . .
Materiale Online . . . . . . . . . . . . . . . . . . . . . . . . . . . .
89
Elenco delle figure
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
Andamento della pressione in funzione dell’angolo di manovella, in assenza di combustione (linea rossa tratteggiata) e con regolare iniezione di combustibile (linea nera continua) nel cilindro di un motore Diesel sovralimentato.
Insieme alla curva della frazione di massa bruciata x b e della legge di rilascio di energia d Q b /dϑ permette di distinguere le quattro fasi principali del processo di ,
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Schematizzazione del processo di combustione di una goccia di com bustibile, circondata da un fronte di fiamma, quando si trova ancora
allo stato liquido [5]. . . . . . . . . . . . . . . . . . . . . . . . . . .
Modello concettuale della struttura di una fiamma durante la fase di
combustione mixing-controlled [19].
. . . . . . . . . . . . . . . . .
Andamento del rapporto di equivalenza locale in funzione della
temperatura di fiamma con differenti processi di combustione [9].
.
Schema cinetico semplificato per una fiamma premiscelata metano-
aria stechiometrica, a p=100kPa e T=298 K [28]. . . . . . . . . . .
Schema cinetico semplificato per una fiamma premiscelata metano-
aria ricca, a p=100kPa, T=298 K [28]. . . . . . . . . . . . . . . . .
Dipendenza dalla temperatura del tasso di produzione della potenza termica prodotta Q ˙ prod e persa Q ˙ lost
[24]. . . . . . . . . . . . . . . .
Diagramma di accensione per una miscela di idrocarburi in aria [27].
2.1
2.2
2.3
Andamento di una traiettoria di reazione dalle condizioni iniziali φ 0
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Andamento di due traiettorie di reazione, quella tabulata, che parte da φ 0 , e quella richiesta che parte da φ q
. . . . . . . . . . . .
Esempio di traiettorie nello spazio delle frazioni molari di CO 2 − H per una miscela di CO/H 2
. . . . . . . . . . . . . . . .
vii
viii ELENCO DELLE FIGURE
2.4
Rappresentazione schematica di una fiamma premiscelata con la curva x ( s )
che attraversa il fronte di fiamma [20]. . . . . . . . . . .
3.1
3.2
3.3
Andamento della temperatura per dei PSR tabulati al variare della T u , a 60 bar e Z 0.044 . . . . . . . . . . . . . . . . . . . . . . . . . .
Andamento della funzione β -PDF per diversi valori di α e β
Andamento della temperatura (curva rossa) e del rapporto Y CO /Y CO 2 (curva verde) per un
di condizioni iniziali: T 0 = 800 K, P = 50 bar , Z = 0 .
0602965 (a) e T 0 = 800 K, P = 50 bar , Z = 0 .
0333968
5.1
5.2
5.3
5.4
5.5
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 800 K e Z = 0.024
. . .
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 800 K e Z = 0.026
. . .
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 800 K e Z = 0.018
. . .
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 812.5 K e Z = 0.024
. .
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 1025 K e Z = 0.024 . . .
5.6
5.7
5.8
5.9
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 50 bar, con temperatura di 800 K e Z = 0.024
. . .
Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 80 bar, con temperatura di 800 K e Z = 0.024
. . .
Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura 800 K e Z = 0.024. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura 800 K e Z = 0.026 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.10 Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura di 800 K e Z = 0.018 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.11 Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura di 812.5
K e Z = 0.024 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 Perfectly-Stirred Reactor
ELENCO DELLE FIGURE ix
5.12 Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura di 1025 K e Z = 0.024 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.13 Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 50 bar, temperatura di 800 K e Z = 0.024 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.14 Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 80 bar, temperatura di 800 K e Z = 0.024 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.15 Valori di mixture fraction in funzione del rapporto di equivalenza. .
5.16 Punti operativi scelti per il calcolo dell’ignition delay. . . . . . . . .
5.17 Calcolo degli ignition delay dei punti operativi mostrati in figura
5.16, per un EGR del 35 % . . . . . . . . . . . . . . . . . . . . . . .
5.18 Calcolo degli ignition delay dei punti operativi mostrati in figura 5.16,
senza EGR (colonna di sinistra) e con un EGR del 25% (colonna di destra) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1
6.2
6.3
6.4
6.5
6.6
6.7
6.8
Mesh 2-D del motore in esame in vista laterale (a) e dall’alto (b). .
Andamento di pressione,
e temperatura per la simulazione HCCI con rapporto di equivalenza φ = 0 .
14 . . . . . . . . . . . . .
Andamtento di pressione,
e temperatura per la simulazione HCCI con rapporto di equivalenza φ = 0 .
20 . . . . . . . . . . . . .
Andamtento di pressione,
e temperatura per la simulazione HCCI con rapporto di equivalenza φ = 0 .
28 . . . . . . . . . . . . .
Distribuzione dei valori di temperatura all’interno del cilindro nel caso con φ = 0 .
20 . Nella colonna sinistra sono presenti i valori calcolati
col modello MZWM, a destra con la tabulazione della cinetica chimica. 66
Distribuzione dei valori di monossido di carbonio all’interno del cilindro nel caso con φ = 0 .
20 . Nella colonna sinistra sono presenti i valori calcolati col modello MZWM, a destra con la tabulazione della cinetica chimica.
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Fotografia dell’esterno (a) e rappresentazione dell’interno (b) della camera di prova Sandia. . . . . . . . . . . . . . . . . . . . . . . . .
Mesh 2-D della camera di combustione in esame,vista laterale (a) e dall’alto (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3 Rate of Heat Release
x ELENCO DELLE FIGURE
6.9
(a)Legge di iniezione calcolata fornita dall’ ECN 4
(b) confronto tra la penetrazione,calcolata e sperimentale, di liquido e vapore. . . . . . .
6.10 Andamento delle curve di rilascio del calore sperimentale e calcolate con diversi scemi cinetici . . . . . . . . . . . . . . . . . . . . . . . .
6.11 Andamento del LOL sperimentale e calcolato con diversi scemi cinetici 72
6.12 Andamento delle curve di rilascio del calore sperimentale e calcolate con e senza la segregation della mixture fraction. . . . . . . . . . . .
6.13 Andamento delle curve di rilascio del calore sperimentale e calcolate
con e senza la segregation della mixture fraction. . . . . . . . . . . .
6.14 Andamento del LOL sperimentale e calcolato con e senza l’introdu zione della segregation. . . . . . . . . . . . . . . . . . . . . . . . . .
6.15 Chemiluminescenza rilevata a 0.8 ms, 1.2ms e 2.5 ms. La linea blu
indica un livello di luminosità pari al 50% del massimo rilevabile.
6.16 Confronto tra la struttura di fiamma sperimentale e calcolata, con (b) e senza (a) l’uso della β − P DF a 0.8 ms
.
identifica i punti a miscela stechiometrica, quella nera i punti a 1400 K. 75
6.17 Confronto tra la struttura di fiamma sperimentale e calcolata, con (b) e senza (a) l’uso della β − P DF a 1.2 ms
identifica i punti a miscela stechiometrica, quella nera i punti a 1400 K. 76
6.18 Confronto tra la struttura di fiamma sperimentale e calcolata, con (b) e senza (a) l’uso della β − P DF a 2.5 ms
identifica i punti a miscela stechiometrica, quella nera i punti a 1400 K. 77
6.19 Confronto tra il campo di formaldeide sperimentale e calcolato, con (b) e senza (a) l’uso della β − P DF a 0.3 ms
. . . . . . . . .
6.20 Confronto tra il campo di formaldeide sperimentale e calcolato, con (b) e senza (a) l’uso della β − P DF a 0.4 ms
. . . . . . . . .
6.21 Confronto tra valori calcolati e sperimentali di ritardo di accensione al variare della temperatura del vessel. . . . . . . . . . . . . . . . .
6.22 Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della temperatura del vessel. . . . . . . . . . . . . . . . .
6.23 Confronto tra valori calcolati e sperimentali di ritardo di accensione, al variare della % di O 2 nel vessel. . . . . . . . . . . . . . . . . . . .
6.24 Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della % di O 2 nel vessel. . . . . . . . . . . . . . . . . . . .
6.25 Confronto tra valori calcolati e sperimentali di ritardo di accensione, al variare della pressione di iniezione del combustibile. . . . . . . . .
4 Engine Combustion Network 5 After Start Of Injection
ELENCO DELLE FIGURE xi
6.26 Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della pressione di iniezione del combustibile. . . . . . . . .
6.27 Confronto tra valori calcolati e sperimentali di ritardo di accensione, al variare della densità dell’aria nel vessel. . . . . . . . . . . . . . .
6.28 Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della densità dell’aria nel vessel. . . . . . . . . . . . . . .
Elenco delle tabelle
1.1
Costanti di chiusura Standard k − ε . . . . . . . . . . . . . . . . . .
3.1
3.2
Suddivisione del dominio di integrazione per c e Z col metodo dei local vectors.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Confronto tra i tempi di di solo calcolo (cpu time) e dell’intera simulazione (clock time) con e senza meccanismi di velocizzazione. .
4.1
Valori delle costanti di chiusura utilizzate nelle equazioni di 00 g 2 e c 00 f 2
5.1
5.2
5.3
Set-up della tabella per l’iso-ottano usata per le simulazioni PSR.
La temperatura ha un intervallo di tabulazione che non è costante, oltre i 100 K si utilizza un ∆ T di 50 K. . . . . . . . . . . . . . . . .
Punti operativi analizzati per le simulazioni
a pressione e volume costante. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Set-up della tabella per l’n-dodecano usata per le simulazioni PSR a volume costante.
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1
6.2
6.3
Condizioni operative per le simulazioni
HCCI . . . . . . . . . . . . . .
Condizioni operative standard per le simulazioni Spray A.
. . . . .
Condizioni operative analizzate nel corso della validazione parame trica del modello in configurazione Spray. . . . . . . . . . . . . . . .
xiii
Sommario
In questo lavoro di Tesi, è stato studiato il processo di combustione, con partico lare riferimento ai motori Diesel. Con questo intento è stato sviluppato e validato un nuovo solutore all’interno della libreria lib-ICE, sviluppata dal gruppo di motori del Politecnico di Milano, basata sul softwere open source OpenFOAM. Questo modello sfrutta la tabulazione della cinetica chimica, che permette l’utilizzo di schemi cinetici dettagliati ad un costo computazionale molto ridotto. In questo lavoro, la tabella è stata generata a partire da simulazioni di reattori perfettamente miscelati a pressione costante, con la possibilità di utilizzare l’approccio β − P DF , per tenere in considerazione gli effetti della turbulence-chemistry interaction. La validazione del modello è stata eseguita attraverso diverse tipologie di simulazioni.
Per prima cosa è stata verificata la corretta costruzione della tabella attraverso delle simulazioni
PSR . Successivamente ne è stata testata l’accuratezza in caso di
combustione premiscelata all’interno di un motore HCCI. Il modello è risultato in accordo con i dati sperimentali e si è dimostrato, a parità di accuratezza, molto più veloce del modello well-mixed, che utilizza l’integrazione diretta della chimica. Infi ne, utilizzando i dati sperimentali resi disponibili dall’Engine Combustion Network, si è potuto verificare il comportamento del modello con una fiamma turbolenta non premiscelata, in varie condizioni operative. I risultati ottenuti hanno mostrato una buona accuratezza nel rappresentare i trend di variazione dell’ignition delay e del lift-off, oltre che la struttura di fiamma diffusiva. Si è riscontrata però una sottostima del valore assoluto del lift-off.
Parole chiave:
CFD, combustione, Diesel, chimica tabulata, PSR
xv
Abstract
This thesis is focused on the study of the combustion process, with particular reference to the Diesel engines. With this intent a new solver has been developed, and validated, in the lib-ICE library, developed by the ICE group of Politecnico di Milano, based on open source code OpenFOAM. This model uses the tabulation of chemical kinetics, which allows the use of detailed kinetic schemes with a reduced computational cost. In this work, the tables were generated from simulations of perfectly stirred reactor at constant pressure, with the possibility of using the β − P DF approach, to take into account the effects of the turbulence-chemistry interaction. The validation of the model was performed through different types of simulations. First, table generation has been tested through
simulations.
Subsequently the accuracy of the model has been tested in premixed combustion in a HCCI engine. The model meets the experimental data and it shows higly speed-up-factor while maintaining simulation results comparable to solving on-the fly chemical kinetics. Finally, using the experimental data available from Engine Combustion Network, it was possible to verify the behavior of the model with a turbulent non-premixed flame in various operating conditions. The results obtained shows a good accuracy in representing the variation trend of ignition delay and lift-off, as well as the diffusive flame structure. However, it was observed an underestimation of the absolute value of the lift-off.
Keywords:
CFD, combustion, Diesel, tabulated chemistry, PSR
xvi
Introduzione
I motori a combustione interna, in particolare quelli Diesel, grazie ai loro relativamente alti rapporti di compressione e alle basse perdite di rendimento al variare del carico, sono tra le macchine termiche a più alta efficienza. Ciò li ha portati ad una vasta diffusione sia per utilizzo industriale, che per i trasporti, rendendoli uno dei fattori chiave per lo sviluppo della società odierna. Costituiscono tuttavia una delle maggiori fonti di emissioni inquinanti a livello urbano. Negli ultimi anni sono state svolte numerose ricerche, con l’obiettivo di ridurre le emissione di N O x e particolato, molecole molto nocive per la salute. Aumentare l’efficienza dei motori e contemporaneamente ridurne le emissioni è un obiettivo difficile da perseguire. Si stanno quindi studiando, con questo intento, nuove modalità di combustione, nuovi sistemi di trattamento dei gas combusti e nuovi combustibili.
Nello sviluppo di queste tecnologie, la modellazione
svolge un ruolo chiave.
Permette infatti di ottenere, con un notevole risparmio di risorse, tutte le informazio ni necessarie relative ad un determinato processo. Il continuo aumento di precisione richiesto, soprattutto nella previsione degli inquinanti e nella comprensione dell’au toaccensione, ha determinato la necessità di utilizzare simulazioni multidimensionali con schemi cinetici dettagliati. Nella simulazione di un processo complesso come la combustione, ciò ha portato ad un aumento rilevante dell’onere computazionale. Per cercare di ridurre i tempi di calcolo sono state sviluppate diverse strategie. Alcune riguardano la formulazione di nuovi schemi cinetici, in modo da ridurre il numero di equazioni di conservazione da risolvere. Altre intervengono sull’integrazione delle equazioni per la chimica.
La tabulazione della cinetica chimica è una possibile soluzione a questo problema.
Permette infatti di utilizzare schemi cinetici complessi, con una notevole riduzione dell’onere computazionale. Questa tecnica consiste nella creazione preventiva di una tabella, da cui durante la simulazione vengono ricavate le informazioni necessarie alla risoluzione della chimica del fenomeno che si sta studiando.
In questo progetto di tesi è stato analizzato il processo di combustione, in 6 Computational Fluid Dynamics
1
2 Introduzione
particolare all’interno dei motori Diesel, con l’obiettivo di sviluppare e validare un modello che sfrutti la tabulazione della cinetica chimica. Questo approccio è stato implementato all’interno della libreria Lib-ICE, sviluppata dal gruppo di motori del Politecnico di Milano, basata sul softwere open source
Per la validazione del modello sono state simulate diverse configurazioni di fiamme. I risultati così ottenuti sono stati confrontati con i dati sperimentali presenti in letteratura, in particolare per le simulazioni degli spray a volume costante, con quelli resi disponibili dall’Engine Combustion Network.
Il testo della Tesi è così strutturato: • nel primo capitolo è descritta la fenomenologia della combustione, con par ticolare riferimento ai motori Diesel, per poi passare ad una presentazione preliminare delle equazioni di conservazione e dei modelli di turbolenza.
• nel secondo capitolo vengono descritti alcuni modelli per la combustione che prevedono la tabulazione, sia on-line che off-line, della cinetica chimica; • nel terzo capitolo è presentata la tecnica di tabulazione utilizzata e alcune metodologie per la riduzione dei tempi di costruzione; • nel quarto capitolo sono riportate le equazioni di conservazione utilizzate per la risoluzione del dominio
e le funzioni necessarie all’estrapolazione delle informazioni dalla tabella; • nel quinto capitolo vengono riportati i risultati preliminari, ottenuti attraverso simulazioni
a pressione e volume costante; • nel sesto capitolo sono presentati i risultati delle simulazioni in configura zione
e SprayA con particolare attenzione alla capacità del modello di prevedere correttamente il ritardo di accensione e la struttura della fiamma.
7 Open source Field Operation And Manipulation
Capitolo 1 La combustione nei motori Diesel
La combustione all’interno dei motori Diesel è un processo controllato dalle condizioni operative, dalla natura fisico-chimica del combustibile e dalle condizioni di miscelamento della carica.
Figura 1.1: Andamento della pressione in funzione dell’angolo di manovella, in assenza di combustione (linea rossa tratteggiata) e con regolare iniezione di combustibile (linea nera continua) nel cilindro di un motore Diesel sovralimentato. Insieme alla curva della frazione di massa bruciata x b e della legge di rilascio di energia d Q b /dϑ , permette di distinguere le quattro fasi principali del processo di
Il combustibile viene iniettato sotto forma di un fine spray all’interno della camera di combustione, dove vaporizza e si miscela con i gas ad alta pressione e temperatura (aria e gas combusti ricircolati). In corrispondenza delle zone dove il rapporto di miscela è prossimo allo stechiometrico, la miscela aria-combustibile si
3
4 Capitolo 1. La combustione nei motori Diesel
accende spontaneamente dopo un ritardo dell’ordine del millisecondo (Fase AB). Il conseguente aumento di pressione nel cilindro accelera le reazioni di preossidazione della miscela aria-combustibile formatasi; questa brucia rapidamente innalzando ulteriormente temperatura e pressione con conseguente evaporazione del combusti bile liquido residuo (Fase BC). L’iniezione continua fino a quando la quantità di combustibile richiesta è stata introdotta, con una combustione regolata dal processo di diffusione (Fase CD). La quantità di energia liberata può essere controllata modulando la legge di iniezione, eventualmente dividendola in più parti. Durante la fase di espansione il miscelamento dell’aria rimasta nel cilindro con i gas combusti o parzialmente ossidati porta al completamento delle reazioni.
1.1
Ritardo d’accensione
Grande importanza assume lo studio del ritardo di accensione, in quanto respon sabile della successiva fase di combustione a volume quasi costante. Dalla quantità di carica premiscelata che si viene a formare dipende infatti la velocità di incremento della pressione e il suo massimo valore, elementi che influenzano la rumorosità, le vibrazioni e le sollecitazioni meccaniche e termiche dei vari componenti. È definito come l’intervallo temporale compreso tra l’apertura dell’iniettore
e l’inizio della combustione
. In letteratura sono presenti varie strategie per l’identificazione
del
• l’istante in cui l’incremento di pressione dovuto alla combustione fa staccare la curva della pressione da quella di semplice compressione oltre una certa soglia (1%); • l’istante in cui abbiamo il massimo della sua derivata; • l’istante in cui appare la prima fiamma nella camera di combustione rilevata con un sensore di luminosità.
Il processo che porta all’accensione di un combustibile liquido, iniettato sotto forma di spray in una massa di aria calda, segue leggi simili a quelle dell’accensione di una miscela omogenea. La differenza principale sta nel fatto che il ritardo complessivo è la somma di due tipi di processi: uno di natura fisica, legato al tempo necessario affinché si venga a formare la miscela di aria combustibile; l’altro di natura chimica, dovuto alle basse velocità di reazione nelle prime fasi del processo di combustione.
Scendendo maggiormente nel dettaglio, attraverso i processi fisici otteniamo una 1 Start of Injection 2 Start of Combustion
1.2. Combustione diffusiva 5
progressiva disintegrazione del getto di combustibile in fini gocce che, riscaldate dall’aria, evaporano e diffondono all’interno dell’aria. Su questi processi influiscono soprattutto la polverizzazione dello spray, le condizioni di moto dell’aria, nonché la sua pressione e temperatura. I processi chimici invece, sono responsabili, per effetto termico, della decomposizione degli idrocarburi a più alta massa molecolare in composti più leggeri, del loro attaccamento da parte dell’ossigeno con formazione
di composti intermedi e dell’avvio delle reazioni a catena (1.4). La natura del
combustibile è molto importante, soprattutto per i primi stadi, dove la grandezza e la forma delle molecole lo rendono più o meno attaccabile dall’ossigeno.
1.2
Combustione diffusiva
Avvenuta l’accensione e avviata la combustione della carica premiscelata, la fiamma si propaga verso il centro dello spray dove trova gocce di combustibile di dimensioni maggiori. Queste ricevono calore dal fronte di fiamma, evaporano velocemente e i vapori diffondono nella carica circostante. Molte gocce, specialmente nella zona centrale e in condizioni di pieno carico, vengono circondate dal fronte di fiamma mentre sono ancora allo stato liquido. Nella figura sottostante è rappre sentata una visione semplificata in cui abbiamo una sola goccia di combustibile, normalmente si hanno nuvole di gocce con fronte di reazione fortemente influenzato dal moto della carica.
Figura 1.2: Schematizzazione del processo di combustione di una goccia di combustibile, circondata da un fronte di fiamma, quando si trova ancora allo stato liquido
6 Capitolo 1. La combustione nei motori Diesel
I gas caldi che circondano la goccia forniscono il calore necessario alla sua graduale evaporazione, i vapori così generati sulla superficie si diffondono nell’aria circostante formando una miscela aria-combustibile che, in corrispondenza di com posizioni prossime al valore stechiometrico, è in grado di bruciare. Questo processo è indipendente dal rapporto aria-combustibile globale, per cui, un motore Diesel ai bassi carichi può funzionare con miscele molto magre. Alle alte temperature in cui si opera durante la parte centrale della combustione, i tempi caratteristici della cinetica chimica sono più piccoli di quelli legati all’evaporazione e alla diffu sione. La velocità di reazione durante questa fase, detta mixing-controlled, sarà quindi limitata dai fenomeni fisici, in particolare dal miscelamento dei vapori di combustibile con l’aria legato al moto turbolento della carica.
La fiamma diffusiva non si estende fino all’uscita dell’iniettore, ma rimane ad una certa distanza, detta Lift-off lenght. Dato che la fiamma è sollevata, aria e combustibile si miscelano a monte della fiamma, con la formazione di una zona di combustione premiscelata vicino al
Figura 1.3: Modello concettuale della struttura di una fiamma durante la fase di
combustione mixing-controlled [19].
Dopo il transitorio iniziale, durante il quale si ha la penetrazione dello spray con la conseguente accensione e combustione premiscelata, la struttura della fiamma rimane pressoché invariata fino al termine dell’iniezione.
3 Lift-off Lenght
1.3. Modalità avanzate di combustione
1.3
Modalità avanzate di combustione
7
Il processo fino a questo punto descritto, e in generale il motore Diesel, è caratterizzato da buoni consumi e basse emissioni di anidride carbonica. Più rilevante è invece la produzione di ossidi d’azoto ( N O x
quali risultano soggetti a controlli sempre più stringenti. Per ridurre le emissioni si può operare in tre direzioni: combustibili alternativi, processi di combustione innovativi e ottimizzazione delle tecniche di after-treatment dei gas di scarico. Passi in avanti sono stati fatti in tutti e tre gli ambiti, ma in un’ottica di riduzione della complessità e dei costi del motore, la ricerca si è maggiormente applicata sul controllo della produzione delle specie inquinanti attraverso l’applicazione di
tecniche di combustione a bassa temperatura ( LTC 5
1.4, questi processi operano a temperature inferiori rispetto a quelle caratteristiche
della formazione degli N O x e con un rapporto di equivalenza locale più basso rispetto a quello in cui si ha la produzione del
Figura 1.4: Andamento del rapporto di equivalenza locale in funzione della temperatura
di fiamma con differenti processi di combustione [9].
Possiamo distinguere due modalità principali di
quale la combustione è governata dalla cinetica chimica e scollegata dall’iniezione; l’altra, detta
, nella quale la combustione è notevolmente influenzata dalla
legge d’iniezione. Recentemente è stata proposta una nuova metodologia detta
, che ha le potenzialità di superare i limiti delle precedenti.
4 Particulate Matter 5 Low Temperature Combustion 6 Premixed Charge Compression Ignition 7 Reactivity Controlled Compression Ignition
8 Capitolo 1. La combustione nei motori Diesel
1.3.1
HCCI
La modalità di combustione
è derivata dal tentativo di combinare il grande vantaggio dei motori
di bruciare una carica omogenea, con i vantaggi di quelli
. Il combustibile viene iniettato molto prima dell’inizio della combustione, in
modo che possa verificarsi il miscelamento aria-combustibile. L’accensione della miscela così ottenuta avverrà simultaneamente in più punti mostrando due picchi di rilascio del calore: il primo legato alle fiamme fredde, il secondo relativo alla combustione completa. L’autoaccensione è controllata dalla cinetica chimica a bassa temperatura e il secondo picco di rilascio dall’ossidazione del CO. Il grande vantaggio della combustione
è la riduzione delle emissioni di
e N O x e una maggior efficienza. é però molto difficile controllare l’istante di autoaccensione. Per questo sono state studiate numerose strategie, sia per l’iniezione (multiple-injection, late e early direct injection, port fuel injection, narrow-angle injection) che per il
interno o esterno, valvole a fasatura variabile, temperatura d’ingresso variabile). Un’ultima strategia è quella di combinare combustibili con diversi punti di autoaccensione, in modo che, variandone la quantità, si possa spostare l’istante di ignizione. In un’ipotetica applicazione, queste tecniche non verrebbero usate singolarmente, ma combinate tra loro, con ottime possibilità in termini di controllo.
1.3.2
PCCI
Questa modalità di combustione deriva dall’ HCCI
con l’intento di avere un maggior controllo sull’istante del
e sul rumore dovuto all’accensione simultanea in più punti con elevato picco di rilascio del calore. Viene utilizzata una carica non perfettamente omogenea, ottenendo l’istante di accensione desiderato attraverso un’elevata turbolenza in camera, ridotto rapporto di compressione, maggior pressione di iniezione e massiccio uso di
EGR . Il combustibile è iniettato in camera attraverso
tre strategie di iniezione possibili: port-fuel injection e early e late fuel injection. In relazione alle prime due tecniche, può accadere che parte del combustibile si depositi sulle pareti del cilindro, con conseguente aumento degli idrocarburi incombusti e del
CO. Tale problema è risolvibile mediante un aumento dell’ EGR
e con la variazione dell’angolo di iniezione. La terza strategia, la late fuel injection, è invece importante per poter tornare, agli alti carichi, alla modalità di combustione convenzionale.
8 Spark Ignition 9 Compression Ignition 10 Exhaust Gas Recirculation
1.4. Cinetica di reazione 9
Lo spazio di miglioramento è legato all’aumento del carico a cui si può applicare questo stile di combustione. Un nuovo approccio alla combustione PCCI prevede un’early injection, che determina la formazione della miscela aria-combustibile, seguita da una late injection che determina l’inizio della combustione. Questa miscela aria-combustibile variabile, evita che si verifichi un’accensione simultanea, ottenendo un miglior controllo sulla combustione e sul rilascio del calore.
1.3.3
RCCI
La più recente applicazione dell’iniezione di diversi combustibili per regolare l’istante del
SOC , prende il nome di RCCI . Tale metodo prevede l’iniezione di
un combustibile con reattività relativamente bassa (port injection), che avviene con molto anticipo rispetto alla combustione in modo da ottenere una miscela aria-combustibile perfettamente omogenea. Successivamente, direttamente dentro al cilindro, viene iniettato un combustibile molto più reattivo, creando così zone più ricche e reattive che determineranno l’inizio della combustione.
Questa metodologia permette di ridurre l’ EGR
agli alti carichi, riducendo il picco di rilascio di calore e controllando il tempismo dell’accensione. Per ottenere il massimo però, sarà necessario poter utilizzare un ampio spettro di miscele di combustibili, dato che la reattività del sistema è influenzata dalla condizione operativa a cui stiamo lavorando.
1.4
Cinetica di reazione
Una generica reazione chimica può essere sempre scritta nella forma: ν a A a + ν b A b + . . .
ν c A c + ν d A d + . . .
(1.1) dove ν i rappresenta il coefficiente stechiometrico della specie i-esima. Tale reazione, su scala microscopica, procede generalmente in entrambi i versi, mentre su scala globale è diretta sempre verso l’equilibrio, dato dal bilancio tra reazione diretta e inversa. Le informazioni relative a quanto velocemente si arrivi alle condizioni
di equilibrio sono ottenuti dall’analisi della cinetica chimica. Per la reazione (1.1)
la variazione temporale della concentrazione di una specie generica c , può essere espressa come: d [ A c ] dt = ν c k f [ A a ] ν a [ A b ] ν b | {z f orward } − k r [ A c ] ν c [ A d ] ν d | {z reverse } (1.2)
10 Capitolo 1. La combustione nei motori Diesel
Dove k f e k r sono rispettivamente i coefficienti dei reaction rate della reazione diretta e inversa. Questi coefficienti possono essere determinati sperimentalmente, o ricavati da formule empiriche come ad esempio quelle di Arrhenius, in cui è evidente la forte dipendenza dalla temperatura: k = AT b exp E act RT (1.3) Le costanti A e b , come anche l’energia di attivazione E act sono tabulate a partire da dati sperimentali.
La combustione, anche di un semplice idrocarburo come il metano ( CH 4 ) è descritta da un meccanismo di reazione molto complesso, che coinvolge un gran numero di specie e reazioni differenti. Per idrocarburi più complessi, come quelli tipicamente usati nei motori Diesel, i meccanismi di reazione diventano molto più complicati. Per l’n-eptano, tipicamente usato come surrogato del Diesel a causa delle simili proprietà di accensione, i meccanismi dettagliati includono centinaia di specie e migliaia di reazioni. La soluzione diretta di questi complessi sistemi di reazioni non è ancora compatibile con le necessità industriali a causa degli elevati tempi di calcolo. Di conseguenza, è utile ricorrere a degli schemi ridotti, con un numero molto più basso di specie e reazioni, ma sempre in grado di descrivere accuratamente le principali caratteristiche del processo di combustione. Per un idrocarburo la reazione generale di combustione è: C x H y + x + y 4 O 2 → xCO 2 + y 2 H 2 O (1.4) e il reaction rate globale può essere approssimato nella seguente formula, come
d [ C x H y ] dt = − A exp − E act RT [ C x H y ] m [ O 2 ] n (1.5) in cui A , m , n e E act sono parametri tabulati, ricavati sperimentalmente. Questi meccanismi globali permettono di ottenere stime approssimative sul rilascio del calore e sulle variazioni di composizione, ma non sono in grado di fornire informazioni accurate sui processi di ossidazione delle specie intermedie, fondamentali per la formazioni degli inquinanti. Per questa ragione, sono necessari meccanismi a più reazioni, in modo da controllare l’evolvere delle specie intermedie che più interessano. In letteratura sono presenti numerose tecniche per la riduzione de meccanismi cinetici, generalmente sono basati sull’osservazione che le reazioni procedono a velocità diverse, anche di ordini di grandezza. Questo permette di ignorare quelle più veloci concentrandosi su quelle che controllano maggiormente la
1.4. Cinetica di reazione
reazione globale.
11
Figura 1.5: Schema cinetico semplificato per una fiamma premiscelata metano-aria
stechiometrica, a p=100kPa e T=298 K [28].
Figura 1.6: Schema cinetico semplificato per una fiamma premiscelata metano-aria ricca,
Nelle figure 1.5 e 1.6, sono mostrati due esempi di schemi ridotti, si può notare
però come questi siano fortemente dipendenti dalle condizioni al contorno. Le ipotesi semplificative che li hanno generati sono valide solo per determinate condizioni operative (temperatura, pressione, rapporto di equivalenza) e configurazioni di fiamma (diffusiva, premiscelata). Una volta che sono state determinate le reazioni più importanti, possono essere fatte ulteriori semplificazioni assumendo che alcune specie siano in equilibrio parziale tra loro, dato che le reazioni che le coinvolgono
12 Capitolo 1. La combustione nei motori Diesel
sono molto veloci. Si può anche assumere che altre specie siano in stato quasi stazionario. Questo descrive la condizione in cui le specie appena formate vengono immediatamente decomposte, in modo che la variazione della loro concentrazione sia globalmente nulla.
I processi di accensione possono essere divisi in due categorie: esplosioni termiche e reazioni chimiche a catena. La prima tipologia può essere spiegata attraverso
la teoria di Semenov [24] , secondo cui l’andamento della temperatura all’interno
di un sistema è dato dalla differenza tra l’energia termica prodotta dalle reazioni chimiche e quella persa a causa dello scambio termico con l’esterno.
ρV c p ∂T ∂t = Q ˙ prod − Q ˙ loss (1.6) Il termine di produzione segue una legge di tipo Arrhenius, mentre lo scambio con l’esterno è stimato usando l’ipotesi di scambio convettivo.
Q ˙ prod ∼ [ f uel ] ν F A exp − E act RT (1.7) Q ˙ lost = h A w ( T − T w ) (1.8) dove h è il coefficiente convettivo e A w è l’area di scambio. Vediamo come il Q ˙ prod cresca esponenzialmente con la temperatura, mentre Q ˙ loss linearmente.
Figura 1.7: Dipendenza dalla temperatura del tasso di produzione della potenza termica prodotta Q ˙ prod e persa Q ˙ lost
Le due curve, come mostrato nella figura 1.7, possono intersecarsi in due punti.
S 1 è detto stabile in quanto per temperature inferiori a T S 1 la produzione è maggiore delle perdite e tende ad avvicinare il sistema ad S 1 per temperature T S 1 < T < T S 2 le perdite sono superiori alla produzione, quindi il sistema tende a raffreddarsi tornando sempre ad S 1 . Il secondo punto di intersezione S 2 è detto instabile in quanto se T > T S 2 le perdite sono sempre inferiori alla produzione con un
1.4. Cinetica di reazione 13
conseguente aumento incontrollato della temperatura. In figura 1.7 è mostrata
anche una curva di produzione che non ha punti di intersezione con quella delle perdite. In questo caso non esistono punti stabili, ma per qualsiasi temperatura del sistema si va in contro ad un’esplosione.
In contrasto con le esplosioni termiche le reazioni a catena non sono direttamente funzione della temperatura del sistema, ma piuttosto del procedere delle reazioni.
L’evolversi delle reazioni e il corrispondente rilascio di energia viene infatti accelerato quando il numero di radicali presenti nel sistema aumenta. Ad esempio, gli step fondamentali nella combustione dell’idrogeno in ossigeno sono: (i) H 2 + O 2 (ii) (iii) + H 2 = H 2 O + O 2 = (iv) (v) (vi) + H 2 = = 1 2 H 2 + O 2 + M = HO 2 + M La reazione (i) è quella che da inizio alle reazioni a catena (chain initiation), dato che crea dei radicali a partire da molecole stabili. Attraverso la reazione (ii) i radicali si trasformano senza però che ne venga modificato il numero (chain propagation), mentre con le reazioni (iii) e (iv) si ottiene l’incremento di radicali necessario all’esplosione delle reazioni a catena (chain branching). Infine i radicali tornano a formare molecole stabili (chain termination) attraverso le reazioni (v) e (vi), ricombinandosi rispettivamente a parete e in fase gas.
Il processo di accensione di un generico idrocarburo può essere ricondotto al modello sopra descritto, ma risulta molto più complesso in termini di aumento del numero di reazioni e di specie coinvolte. Nella prima fase del processo di accensione, le reazioni procedono lentamente dato che sono coinvolte molecole stabili. Di conseguenza la temperatura del sistema cresce lentamente. Solo dopo il superamento di una certa soglia si ha l’inizio dei fenomeni di propagazione e ramificazione dei radicali, che portano in breve tempo all’esplosione. Il tempo necessario alla fase iniziale è fortemente influenzato dalla temperatura, e può essere ridotto significativamente con l’apporto di energia dall’esterno. Molto importante è conoscere il rapporto di equivalenza oltre il quale una miscela è infiammabile. Tale limite è funzione sia della temperatura che della pressione.
14 Capitolo 1. La combustione nei motori Diesel
Figura 1.8:
Diagramma di accensione per una miscela di idrocarburi in aria [27].
Per una data temperatura e una pressione molto bassa, la miscela non si accende.
Ciò avviene perché tutti i radicali che si sono formati diffondono verso le pareti dove ritornano ad essere molecole stabili. Il trasporto diffusivo è inversamente proporzionale alla pressione, quindi molto veloce in queste condizioni. Aumentando la pressione, lasciando costante la temperatura, la diffusione rallenta fino a che la velocità di formazione di radicali supera quella con cui si ricombinano, arrivando ad avviare le reazioni a catena. Incrementando ulteriormente la pressione, la miscela torna ad essere stabile. La reazione di ricombinazione nella fase gas diviene infatti dominante confrontata con quelle di ramificazione e propagazione. Il terzo limite di esplosione si ha aumentando ancora la pressione e identifica l’esplosione termica.
Questo limite dipende infatti dall’equilibrio tra la produzione di energia legata alle reazioni chimiche e lo scambio termico con l’esterno. Il tasso di produzione di energia aumenta con l’aumentare della pressione, ottenendo nuovamente una miscela infiammabile. Per gli idrocarburi possiamo osservare zone di accensione più complesse, specialmente in corrispondenza di basse temperature e alte pressioni. In queste condizioni possiamo avere ignizioni multiple e fiamme fredde. I particolare quest’ultimo fenomeno è osservato per temperature comprese fra 600 e 800 K e prevede una combustione lenta con un incremento di temperatura di circa 100 K.
1.5. Simulazione della combustione Diesel 15
1.5
Simulazione della combustione Diesel
Per lo studio e l’ottimizzazione dei processi di combustione all’interno dei motori per utilizzi industriali, è necessario ricorrere, oltre che a campagne sperimentali, anche a simulazione numeriche. Queste permettono, con un notevole risparmio di costi, di testare condizioni scarsamente riproducibili e ottenere dati difficilmente reperibili con prove sperimentali. Dalla descrizione fenomenologica proposta nei paragrafi precedenti, risulta evidente come una simulazione adeguata richieda modelli dettagliati per tutti i sottro processi con cui abbiamo a che fare. I modelli multidimensionali del fenomeno sono basati sulle equazioni fondamentali della termo fluidodinamica le quali, nel corso degli anni, sono state sviluppate ed affiancate a modelli semplificati in grado di ridurre l’onere computazionale delle simulazioni senza perdere di accuratezza.
1.5.1
Equazioni di conservazione
Come detto in precedenza, i modelli multidimensionali si basano su un set di equazioni che descrivono la fisica del problema. Il sistema di base che viene usato è composto dalle equazioni di Navier-Stokes. Di seguito vengono presentate nel caso di flusso reagente tridimensionale.
• Conservazione della massa ∂ρ ∂t + ∇ · ( ρ u ) = ˙ (1.9) Con questa equazione andiamo ad imporre la conservazione della grandezza scalare ρ , che rappresenta la densità della fase gas. Il termine di sorgente S ˙ rappresenta la porzione di combustibile, che evaporando, va ad aggiungersi alla miscela gassosa.
Dato che le specie chimiche sono continuamente create e distrutte, col pro gredire delle reazioni, è necessario affiancare all’equazione di conservazione globale l’equazione di conservazione della massa di ognuna di esse.
∂ρY k ∂t + ∇ · ( ρY k u ) = ∇ · ( ρD i ∇ Y k ) + ˙ k + ˙ k (1.10) Dove ∇ · ( ρD k ∇ ( Y k )) descrive il trasporto diffusivo della specie k − esima , mentre il termine sorgente ω ˙ k è legato alle reazioni chimiche e S ˙ k all’evapora zione della fase liquida.
16 Capitolo 1. La combustione nei motori Diesel
• Conservazione della quantità di moto ∂ ( ρ u ) + ∇ · ( ρ uu ) = ρ g − ∇ p + ∇ · τ ∂t (1.11) Deriva dalla seconda legge di Newton, secondo la quale la variazione della quantità di moto nel tempo è uguale alla risultante delle forze agenti sul fluido nelle tre direzioni cartesiane. Nel secondo membro possiamo notare la presenza del tensore degli sforzi τ , le nove componenti che lo costituiscono sono incognite, ma per un fluido newtoniano possono essere espresse in funzione del tasso di deformazione.
• Conservazione dell’energia τ ij = 2 µs ij − 2 3 δ ij ∇ · u (1.12) Esistono diverse formulazioni di questa equazione, che dipendono dalla gran dezza fisica che andiamo a conservare. Nelle applicazioni motoristiche molto diffusa è la conservazione dell’entalpia specifica: ∂ ( ρh ) ∂t − ∂p ∂t + ∇ · ( ρh u + j q ) = ∇ · ( p u ) − p : ∇ u + q r + q s (1.13) Dove q r che rappresenta l’apporto dello scambio termico radiativo e q s quello dello scambio termico con lo spray.
1.5.2
Turbolenza
Durante la combustione Diesel, come nella maggior parte dei processi di com bustione, il moto del fluido è di tipo turbolento.
Possiamo descriverlo come instazionario, tridimensionale e caratterizzato da vortici che si sviluppano, a partire dal flusso medio, su una vasta scala dimensionale con un ampio spettro di frequenze.
L’origine del fenomeno risiede nel fatto che ad elevati numeri di Reynolds, le forze viscose non sono più in grado di controbilanciare le forze inerziali, generando così fluttuazioni su larga scala di tutte le grandezze caratterizzanti il flusso.
La simulazione diretta ( DNS 11
) di un flusso turbolento risulta tutt’ora molto
complessa e onerosa dal punto di vista computazionale. Ciò è dovuto al fatto che la dimensione delle piccole scale richiede una risoluzione della mesh molto alta; allo stesso modo, l’alta frequenza dei piccoli vortici richiede un passo temporale molto 11 Direct Numerical Simulation
1.5. Simulazione della combustione Diesel 17
piccolo. Ciò deriva dalla teoria di Kolgomorov che riuscì anche a stimare l’ordine di grandezza dei rapporti, temporale e spaziale, tra microscala e macroscala: l = Re − 2 4 (1.14) t = Re − 1 2 (1.15) Analizzando le scale turbolente, vediamo come quelle più piccole diventino sempre più indipendenti dalla geometria del sistema. Possiamo allora identificare una dimensione da cui il flusso può essere considerato isotropo. Questo porta a differenti approcci per modellare i flussi turbolenti, che differiscono per la scelta della dimen sione minima della scala turbolenta che risolvono. Le simulazioni
arrivano fino alla scala di Kolgomorov, richiedendo quindi un elevatissimo onere computazionale, mentre nelle
poniamo il limite a metà tra la microscala e la macroscala. Per la maggior parte degli utilizzi industriali è sufficiente fermarsi alle scale più grandi,
risolvendo le equazioni per le grandezze medie ( RANS 13
La logica
prevede di considerare ogni grandezza come la somma di due termini: il valore medio all’interno di un certo intervallo temporale e la componente fluttuante legata alla turbolenza.
q ( r , t ) = q ( r , t ) + q 0 ( r , t ) (1.16) I processi di combustione sono però caratterizzati da forti variazioni della densità, è quindi utile introdurre un altro tipo di media, detta di Favre, pesata rispetto alla densità. Per una generica grandezza q risulta: = ρq ρ (1.17)
come nell’equazione (1.16). Possiamo quindi separare la componente media da
quella fluttuante.
q ( − ) = ˜ ( → ) + q 00 ( r , t ) (1.18) Andando a mediare le equazioni di bilancio del sistema, otteniamo dei termini incogniti aggiuntivi legati alle componenti fluttuanti. Di particolare importanza sono gli sforzi di Reynolds ρ ] i j , presenti nell’equazione mediata della quantità di moto: 12 Large-Eddy simulation 13 Reynolds-Averaged Navier-Stokes
18 Capitolo 1. La combustione nei motori Diesel
∂ ( ρ ˜ i ) ∂t + ∂ ( ρ ˜ i ˜ j ) ∂x i = ∂ρ − ∂x j + ∂ ∂x i ( τ ij − ρ ] i j ) (1.19) Tali sforzi sono determinati attraverso modelli di turbolenza, tra questi quello più utilizzato nelle simulazioni dei processi di combustione è il k − ε . Sfruttando l’ipotesi di Boussinesq, gli sforzi di Reynolds vengono considerati proporzionali alle rispettive componenti del tensore degli sforzi τ ij : ρ ] i j = − µ t ∂ ˜ i ∂x j + ∂ ˜ j ∂x i − 2 3 δ ij ∂ ˜ k ∂x k + 2 3 ρk (1.20) Da qui rimane da determinare µ t , la viscosità dinamica turbolenta. Storicamente,a seconda del numero di equazioni di trasporto che vanno risolte per determinarla, sono stati proposti diversi approcci. I modelli a zero equazioni la determinano attraverso un’unica equazione algebrica, ad esempio quella della lunghezza di miscelazione di Prandtl lega µ t al gradiente della velocità: µ t = ρl 2 m | ˜ | (1.21) dove rappresenta il tensore medio degli sforzi: ˜ ij = 1 2 ∂ ˜ i ∂x j + ∂ ˜ j ∂x i (1.22) e l m è la lunghezza di miscelazione data da relazioni empiriche che sono fortemente dipendenti dalla geometria del flusso.
Formulazioni più generali richiedono la soluzione di un’equazione di trasporto, ad esempio dell’energia cinetica turbolenta k come nel modello Prandtl-Kolmogorov.
µ t = ρC µ l pk √ k (1.23) Con la costante C µ = 0 .
09 e l pk lunghezza caratteristica determinata empiricamente.
I modelli più diffusi sono quelli a due equazioni, in particolare per la modellazione della combustione il modello k − ε . In questo approccio la viscosità turbolenta è espressa secondo la relazione: scalar dissipation rate ε : µ t = ρC µ k 2 ε (1.24) Le due equazioni necessarie sono quelle dell’energia cinetica turbolenta k e dello ∂ ( ρk ) ∂t + ∇ · ( ρ ˜ k ) − ∇ 2 ( D ε ef f ε ) = C 1 G ε k − C 3 ε ρ ∇ · (˜ ) − C 2 ρ ε k ε (1.25)
1.5. Simulazione della combustione Diesel 19
∂ ( ρε ) + ∇ · ( ρ ˜ ε ) − ∇ 2 ( D k ef f ∂t k ) = G − ρ ε k k (1.26) Questo modello necessita di 6 costanti di chiusura, e in caso di flusso comprimibile, del numero di Prandtl turbolento.
Tabella 1.1: Costanti di chiusura Standard k − ε Costante Valore C µ C 1 C 2 C 3 σ k σ ε P r t 0 0 0 .
.
.
09 1 .
5 1 .
92 33 1 .
0 1 .
4 7 Questo modello è molto utilizzato anche per la sua semplicità nonché per il basso rapporto costo-beneficio; ne esistono tuttavia altre versioni che cercano di migliorarne alcuni aspetti. Ad esempio il modello k − ε
utilizza delle equazioni derivate da quelle di Navier-Stokes con metodi statistici rigorosi. La formulazione così ottenuta presenta un termine di sorgente aggiuntivo nell’equazione di ε e una diversa formula per il P r t che migliora le previsioni dei flussi di transizione, con grande swirl, e scambio termico e di massa a parete.
1.5.3
Determinazione del reaction rate
Uno dei termini da modellare che derivano dalle equazioni mediate è la velocità di reazione media ˙ k , presente nell’equazione di conservazione per una generica specie k : ∂ ( ρ Y ˜ k ) + ∂t ∂ ( ρ ˜ i Y ˜ k ) ∂x i = − ∂ ∂x i ( V k,i Y k ij − ρ ] i k 00 ) + ˙ k (1.27) La descrizione completa delle reazioni chimiche che costituiscono la combustione coinvolge centinaia di specie e migliaia di reazioni, gestire un processo così complesso è tutt’ora molto oneroso: • è necessaria un’equazione di bilancio per ogni specie presente nel sistema; 14 Renormalization Group Method
20 Capitolo 1. La combustione nei motori Diesel
• i reaction rates delle varie reazioni sono funzioni complesse, nel caso di una semplice reazione tra un combustibile (F) e un ossidante (O), F + sO → (1 + s ) P (1.28)
sono espressi secondo una legge di tipo Arrhenius (eqn. 1.2), fortemente non
lineare. Di conseguenza non è facile derivarne il valore medio come funzione delle frazioni massiche medie del combustibile ( Y ˜ F ), dell’ossidante ( Y ˜ O ), della temperatura media ( ˜ ) e della densità media ( ρ ); • un ulteriore difficoltà è data dall’accoppiamento tra turbolenza e combustione.
I due fenomeni avvengono infatti su scale temporali differenti, con la necessità di utilizzare time step molto piccoli per poterne descrivere correttamente l’interazione.
Per cercare di dare risposta a questi problemi la chimica e la fluidodinamica vengono risolte separatamente. Per la determinazione dei reaction rates si sono sviluppate due strategie: l’integrazione diretta delle equazioni di conservazione e la tabulazione della cinetica chimica. Nel primo caso, i solutori ottengono i reaction rates integrando le equazioni, in accordo con i meccanismi cinetici, utilizzando intervalli temporali dell’ordine di 10 − 7 − 10 − 9 . Un approccio di questo tipo offre grande flessibilità sia per il tipo di combustibile che per le condizioni operative, ma risulta molto oneroso dal punto di vista computazionale. La tabulazione della chimica invece, prevede l’utilizzo di tabelle in cui sono contenuti i reaction rates, la composizione chimica e altre grandezze utili, come temperatura e pressione, in funzione di una quantità detta progress variable. Nota quest’ultima grandezza, dall’integrazione dell’equazione di conservazione corrispondente, siamo in grado di ricostruire la composizione e i reaction rates in ogni cella. Nel successivo capitolo verranno presi in considerazione più nel dettaglio alcuni modelli che sfruttano la chimica tabulata, in particolare quelli off-line.
Capitolo 2 Tabulazione della cinetica chimica
L’utilizzo di schemi cinetici dettagliati delle reazioni sta diventando un pre requisito fondamentale per poter ottenere delle simulazioni accurate inerenti alla combustione nei motori a combustione interna. Grazie a questi schemi è infat ti possibile descrivere il processo in un ampio intervallo di condizioni operative, ottenendo informazioni specifiche sulla formazione degli inquinanti (soot, N O x e incombusti) e sulla legge di rilascio del calore anche in condizioni di funzionamento
2.1
Inquadramento del problema
Per un generico flusso reagente, lo stato termodinamico è individuato dalle frazioni massiche Y k ( k = 1 , 2 , . . . , n s ) delle n s specie, dall’entalpia h e dalla pressione p . Lo stato è definito da: = { φ ˆ 1 , φ ˆ 2 , . . . , ˆ n +2 } = { Y 1 , Y 2 , . . . , Y n s , h, p } .
(2.1) Le componenti di non sono però tutte linearmente indipendenti, in quanto la somma delle frazioni massiche deve essere pari a uno, e andrà rispettata la conservazione degli atomi e dell’entalpia.
Se vi sono n l dipendenti, i gradi di libertà del sistema termochimico sono: variabili linearmente D = n s + 2 − n l 1 In Situ Adaptive Tabulation 2 Intrinsic Low-Dimensional Manifold 3 Flamelet-Generated Manifolds
21
(2.2)
22 Capitolo 2. Tabulazione della cinetica chimica
Possiamo quindi definire la composizione φ ( t ) = { φ 1 , φ 2 , . . . , φ D } , come un vettore o come un punto dello spazio delle composizioni.Si può anche notare come non tutti i punti dello spazio della composizione siano fisicamente realizzabili, in quanto vanno rispettate le seguenti condizioni: T ( h, p, φ ) p > 0 > 0 0 ≤ Y k ≤ 1 (2.3) (2.4) (2.5) Nel dominio fluidodinamico, la composizione evolve seguendo un set di equazioni che può essere riscritto nella forma: dφ ( t ) = S [ φ ( t )] + M ( t ) dt (2.6) S [ φ ( t )] è il termine di sorgente legato alle reazioni chimiche, mentre M ( t ) rappresenta i fenomeni di trasporto che causano il miscelamento. É importante valutare le
diverse scale temporali caratteristiche dei fenomeni descritti dall’equazione 2.6, che
per la combustione variano tra 10 − 9 s a più di 1 secondo. I processi di miscelamento non sono mai più rapidi di 10 − 3 − 10 − 4 s e possiamo quindi separarli dalle reazioni chimiche risolvendo due equazioni separate: dφ ( t ) = M ( t ) dt (2.7) dφ ( t ) = S [ φ ( t )] dt (2.8)
La soluzione dell’equazione 2.8, partendo da una composizione chimica iniziale
φ 0 , corrisponde ad una traiettoria nel dominio delle composizioni, che porta ad una condizione di equilibrio chimico. All’istante t 0 + ∆ t , con ∆ t sufficientemente Figura 2.1: Andamento di una traiettoria di reazione dalle condizioni iniziali φ 0
piccolo, la composizione sarà φ ( t 0 + ∆ t ) = R ( φ 0 ) con R ( φ 0 ) , detto reaction mapping,
2.2. ISAT
funzione solamente della composizione iniziale.
23
2.2
ISAT
La regione da tabulare (accessed region) è l’insieme delle composizioni che si presentano nel flusso o nella simulazione, gruppo ben più piccolo rispetto alla totalità delle possibili composizioni. Queste dipendono da numerosi fattori, quali cinetica, fenomeni di trasporto e condizioni al contorno: non possono quindi essere conosciute a priori. Per questa ragione la tabella non è costruita in una fase di pre-processing, ma durante la simulazione del flusso reagente. Ogni nuova composizione che si ottiene dai calcoli verrà aggiunta.
Le varie voci della tabella contengono la composizione termochimica φ 0 , il reaction mapping R ( φ 0 ) e il mapping gradient A ( φ 0 ) = ∂R i ( φ ) ∂φ j che rappresenta la sensibilità del reaction mapping alla variazione della composizione necessaria per l’interpolazione lineare.
R ( φ q ) ≈ R l ( φ q ) ≡ R ( φ 0 ) + δR l (2.9) con δR l ≡ Aδφ = δR + o ( | δφ | 2 ) (2.10) Figura 2.2: Andamento di due traiettorie di reazione, quella tabulata, che parte da e quella richiesta che parte da φ q
φ 0 , È necessaria anche la dimensione della regione di accuratezza, per cui il reaction mapping viene determinato attraverso un’interpolazione lineare in modo che l’errore sia inferiore ad una certa tolleranza ε max prestabilita. Si assume che tale regione
24 Capitolo 2. Tabulazione della cinetica chimica
sia adeguatamente approssimata dal un iper-ellissoide rappresentato dalla matrice unitaria Q e dalla matrice diagonale Λ .
Quando è necessario un certo reaction mapping R ( φ q ) viene determinata la più vicina condizione R ( φ 0 ) seguendo una logica ad albero binario: • se siamo nella regione di accuratezza di un R ( φ 0 ) esistente, viene attuata una
semplice interpolazione lineare (Eqn. 2.9);
• se siamo fuori dalla regione di accuratezza calcoliamo R ( φ q ) attraverso l’in-
tegrazione dell’equazione 2.8 e confrontando la soluzione con il risultato
dell’interpolazione lineare determiniamo ε . Vanno poi considerati due casi:
–
se ε < ε max ingrandita; la regione di accuratezza della composizione φ 0 viene
–
se ε > ε max viene generata una nuova voce nella tabella corrispondente a φ q .
è un meccanismo efficiente per la riduzione dell’onere computazionale delle simulazioni della cinetica chimica con modelli dettagliati, ma rimane superiore a quello richiesto per schemi ridotti. Per questo sono stati proposti nuovi modelli che combinano la tabulazione con degli schemi che vadano a ridurre i meccanismi cinetici.
Un esempio è il
, che unisce al meccanismo di tabulazione on-line, il
il quale riduce il meccanismo cinetico in funzione delle condizioni termodinamiche locali.
2.3
ILDM
Questo metodo prevede la costruzione di uno schema cinetico semplificato partendo dall’ipotesi che inizialmente l’evoluzione temporale di un sistema chimico dipenda dalle condizioni iniziali, ma dopo un certo intervallo temporale (di solito abbastanza piccolo), questo dipenda da un numero ridotto di variabili. La curva su cui collidono le varie traiettorie è il cosiddetto Low-dimensional manifold. Il formalismo matematico che sta alla base di questa semplificazione è molto solido, si basa su un’analisi agli autovalori e autovettori dello Jacobiano del sistema di equazioni di conservazione che rivela come ci siano n s scale temporali caratteristiche (associate agli autovalori) e n s direzioni preferenziali (associate agli autovettori).
Come già detto nei paragrafi precedenti, le reazioni hanno scale temporali che 4 Tabulation of Dynamic Adaptive Chemistry 5 Dynamic Adaptive Chemistry
2.3. ILDM 25
Figura 2.3: Esempio di traiettorie nello spazio delle frazioni molari di CO 2 miscela di CO/H 2
− H per una coprono un range molto vasto. Assumendo che quelle più veloci raggiungano prima le condizioni di equilibrio, possiamo disaccoppiarle da quelle più lente che regoleranno gli sviluppi successivi del sistema. Otteniamo quindi un sistema caratterizzato da n g = n s − n f gradi di libertà, detto low-dimensional manifold, con le proprietà che ci permettono di confinare lungo lo stesso i movimenti delle reazioni più lente.
Le assunzioni di stato stazionario ed equilibrio parziale sono usate frequentemente per ridurre gli schemi cinetici. La prima prevede che alcune specie, ad esempio alcuni radicali, abbiano raggiunto lo stato di equilibrio. Questo significa che la loro composizione non varierà più nel tempo, e che i loro reaction rates sono trascurabili.
La seconda assume che alcune reazioni siano arrivate all’equilibrio: questo implica che il reaction rate della reazione in un verso è uguale a quello nell’altro. Queste due ipotesi hanno però il grosso svantaggio di essere imposte a priori. É necessaria un’analisi approfondita dei sistemi di reazioni, affiancata da considerazioni sullo stato termochimico del sistema da simulare. Temperature, pressioni e composizioni della miscela influenzano molto l’applicazione delle ipotesi di equilibrio parziale e
stato stazionario. L’ ILDM , attraverso l’analisi agli autovalori, è in grado di attuare
le migliori semplificazioni, dato che queste non sono imposte a priori ma derivano dallo studio del sistema nelle condizioni istantanee. Andranno considerate solo le reazioni che corrispondono agli autovalori positivi, in quanto quelli nulli o negativi rappresentano le scale temporali trascurabili.
Ottenuto il Low-dimensional manifold, dobbiamo scegliere i parametri da usare per descriverlo. Le frazioni massiche e i reaction rate possono essere determinati come funzione di un numero ridotto di variabili, tra cui una "progress variable"
26 Capitolo 2. Tabulazione della cinetica chimica
(definita ad esempio come combinazione lineare delle frazioni massiche di alcune specie),che da una misura del progresso delle reazioni, e una grandezza che descriva il miscelamento, detta mixture fraction, come una frazione massica di qualche specie. Ottenuta la descrizione del low-dimensional manifold, possiamo costruire una tabella in cui ogni voce corrisponde ad un valore delle variabili da cui dipende.
2.4
FGM
L’ILDM è un ottimo metodo per descrivere reazioni ad alta temperatura e vicine all’equilibrio, mentre perde significativamente di accuratezza diminuendo la temperatura, dato che le condizioni non sono tabulate, ma ottenute mediante interpolazione lineare. Questo non è un problema se dobbiamo descrivere fiamme stazionarie, ma lo diventa se vogliamo simulare correttamente fiamme diffusive e fenomeni instazionari come l’autoaccensione, aspetto chiave delle moderne tecniche
di combustione per i motori Diesel ( HCCI
e
per cui una fiamma di qualsiasi struttura può essere rappresentata come un insieme di fiammelle monodimensionali, con un’implementazione che prevede la costruzione di un manifold all’interno dello spazio delle composizioni. Andremo ad usare un certo numero di parametri per descriverlo a seconda dell’accuratezza che vogliamo ottenere.
Il sistema differenziale composto dalle equazioni di conservazione, presentate nel
paragrafo 1.5.1, è numericamente oneroso da risolvere a causa del grande numero di
specie chimiche coinvolte. Inoltre l’ampio spettro di scale temporali, caratteristiche delle reazioni chimiche, porta il sistema ad essere matematicamente rigido (stiff). La semplificazione del problema, proposta con i flamelets, considera come la maggior parte dei processi di diffusione e convezione siano presenti anche nelle fiamme 1D, e che quindi una combinazione di fiamme monodimensionali possa rappresentare adeguatamente strutture di fiamme più complesse. Per analizzare la struttura interna della fiamma andiamo a definire un sistema di riferimento curvilineo, nel quale la coordinata s è localmente perpendicolare alle iso-superfici corrispondenti ad un generico scalare ϕ . Questo potrebbe essere una grandezza qualsiasi che risponda alla condizione ∇ ϕ = 0 , ma normalmente se ne sceglie una che sia monotona crescente o decrescente con s . Una superficie di fiamma è una iso-superfice di ϕ che si muove in accordo con l’equazione Dϕ ≡ Dt ∂ϕ + u f ∂t · ∇ ϕ = 0 , (2.11)
2.4. FGM 27
dove u f è la velocità superficiale di fiamma. Un nuovo set di equazioni di conser vazione può quindi essere derivato, al fine di descrivere la struttura interna delle fiamme in funzione del nuovo sistema di coordinate. In questo modo abbiamo ottenuto un sistema di pari dimensioni rispetto a quello di partenza, ma dato che il fronte di fiamma può essere considerato sottile, il trasporto lungo la direzione tangenziale alla iso-superficie di ϕ è trascurabile rispetto a quello nella direzione normale. Il raggio di curvatura del fronte di fiamma è molto più grande del suo spessore.
In questo modo possiamo andare a risolvere le equazioni per i flamelets separa tamente e andare a tabularne i risultati in funzione di un certo numero di variabili.
Prima di sceglierle è fondamentale definire adeguatamente la tipologia di flamelet che andremo a risolvere. Dovremo usare un modello di combustione che sia rap presentativa del fenomeno che andremo a studiare nelle successive simulazioni. Ad esempio, per ottenere un’accurata rappresentazione della combustione con tecniche
LTC , caratterizzate dall’autoaccensione di miscele omogenee, scegliamo un flamelet
strutturato come un reattore omogeneo. Invece, se fossimo interessati a modellare metodi di combustione convenzionali, potremmo usare flamelet diffusivi. Anche la definizione delle variabili caratterizzanti il manifold dipende dal fenomeno che vogliamo studiare, nella combustione Diesel convenzionale o a bassa temperatura fondamentali sono il miscelamento delle specie e il procedere delle razioni. Dovremo quindi definire almeno due variabili indipendenti: la prima, • mixture fraction Z , rappresenta lo stato locale di miscelamento, della quale vi sono diverse espressioni,tutte caratterizzate dal fatto che questa vari tra 0 (ossidante puro) e 1 (combustibile puro).
• progress variable c , descrive il procedere delle reazioni, ed è definita attraverso Figura 2.4: Rappresentazione schematica di una fiamma premiscelata con la curva x ( s )
che attraversa il fronte di fiamma [20].
28 Capitolo 2. Tabulazione della cinetica chimica
una combinazione lineare delle frazioni massiche di alcune specie, in modo che il suo andamento sia monotono nel passaggi o da una miscela incombusta ad una completamente ossidata.
Capitolo 3 Generazione della tabella
Il modello che è stato studiato e validato in questo progetto di tesi, è implemen tato all’interno della libreria Lib-ICE. Questa è stata interamente sviluppata dal gruppo di Motori a Combustione Interna del dipartimento di Energia del Politecnico di Milano. É basata sul softwere
open source
e contiene numerosi solver e utilities per la simulazione di tutti i processi che avvengono nei motori a combustione interna. Possiamo distinguere due parti principali del modello: la prima permette la costruzione della tabella per la cinetica chimica, la seconda è il so lutore CFD, che lega il dominio fluidodinamico a quello chimico. In questo capitolo
verrà analizzata la tecnica di tabulazione, mentre nel prossimo (4) il solutore.
3.1
Modello PSR
La metodologia di tabulazione, come già visto per il modello
vede la soluzione di flamelets laminari, schematizzati come reattori perfettamente
miscelati ( PSR ) a pressione costante. Ogni PSR
è identificato da un set di variabili che ne costituiscono le condizioni iniziali: temperatura della miscela incombusta T u , pressione e mixture fraction Z.
Ad ogni time-step vengono ricavati i valori dei reaction rates e, risolvendo le equazioni di conservazione delle singole specie e dell’energia, otteniamo il nuovo stato della miscela.
∂ ( ρY i ) ∂t = ˙ i i = 1 , . . . , N s (3.1) ∂h = ˙ c ∂t (3.2) L’assenza di termini avvettivi e diffusivi è giustificata dal fatto che la geometria del sistema è composta da un’unica cella a composizione omogenea. Le variazioni di composizione sono dovute esclusivamente alle reazioni chimiche, rappresentate
29
30 Capitolo 3. Generazione della tabella
dai termini sorgente ω ˙ i e S ˙ c . Nella figura sottostante è rappresentato il classico andamento della curva di temperatura di una simulazione
valori di T u relativamente bassi, può essere distinta una prima fase di reazione a bassa T seguita dall’accensione principale, mentre per valori più alti il processo
avviene in un unica fase. (Fig. 3.1).
Figura 3.1: Andamento della temperatura per dei PSR tabulati al variare della T u , a 60 bar e Z 0.044 .
3.2
Progress variable
Un parametro fondamentale per tutte le metodologie di tabulazione della cinetica chimica è la progress variable. Questa grandezza costituisce infatti il collegamento tra la fluidodinamica e il progredire delle reazioni che trasformano la miscela fresca in gas combusti. La definizione di tale variabile non è però univoca ed errori nella sua formulazione possono portare a risultati privi di significato fisico. Una buona scelta deve rispettare alcuni principi: • deve essere una grandezza trasportabile nel dominio fluidodinamico; • gli scalari usati nella sua formulazione devono evolvere con scale temporali comparabili; • ogni valore che assume nel tempo deve caratterizzare un singolo stato termo dinamico.
3.2. Progress variable 31
Nel presente codice c’è la possibilità di scegliere tra 3 tipologie di progress variable. Possiamo usare: • una combinazione lineare delle frazioni massiche delle specie rilevanti ai fini della combustione: c ( T u , p, Z, t ) = N s X α i Y i ( T u , p, Z, t ) i (3.3) con α i generico coefficiente moltiplicativo per la specie i-esima. Le specie da usare ( speciesN ames ) e i rispettivi coefficienti ( coef f icients ) vengono impostati dall’utente nell’apposita subDictionary multiSpeciesCoef f s ; • l’entalpia sensibile della miscela: c ( T u , p, Z, t ) = h s ( T, t ) • l’entalpia di formazione della miscela a 298 K: c ( T u , p, Z, t ) = h 298 (298 K, t ) (3.4) (3.5) La formulazione della progress variable basata sulle specie, è condizionata dalle dif ficoltà inerenti alla scelta delle specie e alla definizione dei coefficienti moltiplicativi.
Il vincolo principale a cui si deve sottostare è che tale variabile sia monotona col procedere delle reazioni, in modo che, passando da miscela fresca a completamente bruciata, il suo valore normalizzato passi da 0 a 1. In letteratura sono presenti
alcuni metodi per l’ottimizzazione della sua formulazione, come ad esempio [8].
Durante lo sviluppo di questo solver è stata fatta un analisi di sensitività che ha portato alla scelta di una formulazione basata sull’entalpia. L’utilizzo, in particolare
dell’entalpia chimica di formazione, usata in [11] e [10], porta migliori risultati in
termini di predizione del lift-off e simulazione di condizioni con grande
32 Capitolo 3. Generazione della tabella
3.3
Per poter calcolare i termini mediati delle equazioni del trasporto, tenendo conto delle fluttuazioni turbolente di scala inferiore alla griglia di calcolo, viene introdotto un metodo statistico basato sull’impiego di una
φ , la funzione di probabilità è definita come F φ (Ψ) = P rob [Ψ ≤ φ ≤ Ψ + d Ψ] , mentre la probability density function : P φ (Ψ) = dF φ (Ψ) d Ψ Dalla definizione ricaviamo le due proprietà che la caratterizzano: (3.6) • P φ (Ψ) ≥ 0 • R 1 0 P φ (Ψ) d Ψ = 1 Possiamo derivare una formulazione anche per il calcolo di qualsiasi grandezza mediata alla Favre ( f ˜ ), introducendo una
pesata sulla massa( ˜ ): ρ f ˜ = ρf = Z 0 1 ρ (Ψ) f (Ψ) P φ (Ψ) d Ψ = ρ Z 0 1 f (Ψ) ˜ φ (Ψ) d Ψ (3.7) L’utilizzo delle
ha il vantaggio di poter descrivere tutti i flussi reagenti insta zionari determinandone tutte le informazioni necessarie. Il problema è calcolare il valore della
PDF , dato che cambia in ogni punto del dominio. Sono state sviluppate
due metodologie : la risoluzione di equazioni al trasporto per P φ (Ψ) e l’assunzione della tipologia di funzione (presumed PDF) Nel solver descritto in questo capitolo, è stato scelto l’approccio presumed PDF, assumendo come modello la β -PDF (in questo caso riferita alla mixture fraction).
˜ ( Z ) = Z α − 1 (1 − Z ) β − 1 (3.8) Questa soluzione permette di avere una funzione determinata da soli due parame tri α e β , che possono essere calcolati a partire dal valore medio ( ˜ ) e dalla varianza ( 00 g 2 ) di Z. L’utilizzo di questa funzione è giustificato dal fatto che, al variare dei due parametri, può cambiare forma in maniera continua. Questo modello è stato inserito nel codice in modo che ogni
venga calcolato con e senza l’approccio
permettendo la creazione di due tabelle indipendenti. Vi è la possibilità di calcolare la
per la mixture fraction Z (utile nelle simulazioni degli spray), per la progress variable c (utile per le simulazioni HCCI), e per entrambe contemporaneamente. La 1 Probability Density Function
3.3. PDF 33
Figura 3.2: Andamento della funzione β -PDF per diversi valori di α e β
scelta della grandezza a cui fare riferimento è implicita al momento in cui forniamo i valori da tabulare. I valori che andiamo a tabulare non sono però i valori della varianza di c o Z, ma i loro valori segregati: seg c 00 f i 2 = c 00 f i 2 c i (1 − c i ) (3.9) seg g i 00 2 = 00 g i 2 Z i (1 − Z i ) (3.10) Questa scelta è stata fatta per facilitare la costruzione della tabella in quanto La segregazione della varianza porta ad una sua normalizzazione, ottenendo quindi delle variabili comprese tra 0 e 1, i cui valori hanno significato più generale. Per ogni combinazione di c i e c 00 f i 2 o Z i e 00 g i 2 , possiamo calcolare i coefficienti α della
PDF . Per la mixture fraction sono semplicemente dati da:
e β α ( Z i , 00 g i 2 ) = Z i Z i (1 − 00 g i 2 Z i ) − 1 !
β ( Z i , 00 g i 2 ) = (1 − Z i ) Z i (1 − 00 g i 2 Z i ) − 1 !
(3.11) Per la progress variable, sono calcolati con le stesse formule: α ( c i , c 00 f i 2 ) = c i c i (1 − c i ) c 00 f i 2 − 1 !
β ( c i , c 00 f i 2 ) = (1 − c i ) c i (1 − c 00 f i 2 c i ) − 1 !
(3.12)
34 Capitolo 3. Generazione della tabella
Va tenuto conto della possibilità che i coefficienti α e β possano raggiungere valori molto alti. Assunto come valore massimo 500, nel caso di superamento di tale soglia: α β = 500 = 500 β = α = ( α − 1 .
0) − f M ax ( α − 2 .
0) 1 .
0 + f M ax f M ax ( β − 2 .
0) 1 − f M ax (3.13) (3.14) f M ax = 1 + 1 β − 1 α − 1 (3.15) Se, dopo questa operazione, il loro valore è ancora superiore a 500 si impone nulla la segregation della variabile corrispondente: seg g i 00 2 = 0 se α Z,i > 500 o β Z,i > 500 seg c 00 f i 2 = 0 se α c,i > 500 o β c,i > 500 Per l’integrazione della
si usa il metodo dei trapezi e vengono definiti dei local vectors che forniscono i valori su cui effettuare l’integrazione. I punti così definiti, non sono distribuiti uniformemente nell’intervallo di integrazione [0 , 1] , ma avvicinandosi agli estremi il passo di integrazione diminuisce fino ad un minimo imposto dall’utente (in automatico imposto a 10 − 7 ). La creazione dei vettori prevede la divisione del dominio di integrazione in 11 intervalli, ognuno dei quali è a sua volta suddiviso in un determinato numero di parti. Il numero di sotto intervalli Nvec, j 8 9 10 11 1 2 3 4 5 6 7 ∆ [10 − 7 , [10 − 6 , [10 − 5 , [10 − 4 , [1 [10 − 3 , [10 − 2 , 1 − 10 − 2 , [1 [1 [1 [1 − − − − 10 − 3 , 10 − 4 , 10 − 5 , 10 − 6 , 10 − 6 ] 10 − 5 ] 10 − 4 ] 10 − 3 ] 10 − 2 ] − 10 − 2 ] 1 − 10 − 3 ] 1 1 1 1 − − − − 10 − 4 ] 10 − 5 ] 10 − 6 ] 10 − 7 ] N Int 20 20 20 20 20 200 20 20 20 20 20 Tabella 3.1: Suddivisione del dominio di integrazione per c e Z col metodo dei local vectors.
( N int ) può essere modificato andando a definire con dei valori opportuni M P DF , per
3.3. PDF 35
l’intervallo centrale ( j = 6 ), e N P DF per i restanti. I singoli valori di c e Z su cui fare l’integrazione vengono successivamente calcolati, per ogni j-esimo intervallo, come: c i = c i − 1 + δ j (3.16) Z i = Z i − 1 + δ j (3.17) in cui δ j = ∆ j /N j Int . Ottenuti gli andamenti di c e Z possiamo procedere con l’integrazione della
PDF , ottenendo i valori di tutte le grandezze (
φ ) per la nuova tabella: • se segC 2 = 0 e segZ 2 = 0 segφ ( φ, segZ 2) = P DF Z = N c X 1 2 " P DF φ P DF Z ˜ ( Z i +1 P ( Z i ) # Z i +1 − Z i i =0 P DF φ = N c X 2 i =0 1 " ˜ ( Z i +1 ) φ i +1 + ˜ ( Z i ) φi # Z i +1 − Z i (3.18) (3.19) (3.20) • se segZ 2 = 0 e segC 2 = 0 , P DF c = N c X " 1 P DF φ segφ ( φ, segC 2) = P DF c ˜ ( c i +1 P ( c i ) # i =0 2 c i +1 − c i ˜ ( c i +1 ) φ i +1 + ˜ ( c i ) φ ( k ) i # P DF φ = N Z X " 1 i =0 2 c i +1 − c i (3.21) (3.22) (3.23) • se segZ 2 = 0 e segC 2 = 0 , segφ ( φ, segC 2 , segZ 2) = sumP DF φ sumP DF (3.24) Cambia il metodo di integrazione dato che dobbiamo andare a discretizzare un integrale doppio. Per la sola funzione
risulta: P DF = A T r ( P DF Hm 0 + P DF Hm 1 ) P DF Hm 0 = P DF Hm 1 = 1 h ˜ ( c i ) ˜ ( Z i ) + ˜ ( c i ) ˜ ( Z i +1 ) + ˜ ( c i +1 P ( Z i +1 ) i 3 1 h ˜ ( c i ) ˜ ( Z i ) + ˜ ( c i +1 ) ˜ ( Z i +1 ) + ˜ ( c i +1 P ( Z i ) i 3 (3.25) (3.26) (3.27)
36 Capitolo 3. Generazione della tabella
mentre per una generica grandezza da tabulare φ : P DF φ = A T r ( Hm 0 + Hm 1) (3.28) H m 0 = 1 3 h ˜ ( c i ) ˜ ( Z i ) φ ( c i , Z i ) + ˜ ( c i P ( Z i +1 ) φ ( c i , Z i +1 )+ ˜ ( c i +1 P ( Z i +1 ) φ ( c i +1 , Z i +1 ) i (3.29) H m 1 = 1 h ˜ ( c i ) ˜ ( Z i ) φ ( c i , Z i ) + ˜ ( c i +1 P ( Z i +1 ) φ ( c i +1 , Z i +1 )+ 3 ˜ ( c i +1 P ( Z i ) φ ( c i +1 , Z i ) i (3.30) Dove, in entrambi i casi, A T r = 1 2 ( c i +1 − c i )( Z i +1 − Z i ) .
3.4
Mixing Line
Nella modellazione della combustione Diesel un importante fattore da considerare è che il combustibile è iniettato allo stato liquido e ad una temperatura più bassa rispetto all’aria presente nel cilindro. Il suo riscaldamento e la sua evaporazione sono possibili grazie al raffreddamento dei gas che lo circondano. La metodologia di tabulazione che abbiamo fin qui descritto, classifica i
in base alla temperatura della miscela fresca ottenuta combinando combustibile e aria alla stessa temperatura.
Questa approssimazione è lecita se non consideriamo la possibilità che le grandezze
possano variare all’interno di una singola cella (3.3). Il fatto che la mixture fraction
abbia una certa distribuzione, all’interno di ogni cella, implica che la temperatura vari allo stesso modo. Questo perchè, cambiando la quantità di combustibile con cui l’aria si è miscelata, cambia la temperatura della miscela che si è venuta a formare.
Per tenere in considerazione queste variazioni abbiamo inserito la possibilità di calcolare gli effetti della miscelazione e dell’evaporazione al momento di generare la tabella. In questo modo, per ogni
PSR , viene calcolata la nuova temperatura
T P SR = T mix a partire dall’entalpia sensibile della miscela aria-combustibile: h s,mix = (1 − Z ) h s,air ( p, T air ) + Zh s,f uel ( p, T f uel ) (3.31) Dove la temperatura del combustibile viene imposta dall’utente. L’influenza dell’eva porazione del combustibile è inserita separatamente e modifica il calcolo dell’entalpia
3.5. Costruzione della tabella 37
specifica della miscela, andando a togliere il calore latente di evaporazione: h s,mix = (1 − Z ) h s,air ( p, T air ) + Zh s,f uel ( p, T f uel ) − ∆ h evap (3.32) Con il ∆ h evap calcolato come somma dei contributi di tutte le n f s specie da cui è composto, a 1 bar e T = T f uel .
∆ h evap = n f s X x i ∆ h evap,i i =1 (3.33)
3.5
Costruzione della tabella
Ogni volta che viene terminata una simulazione
PSR , il codice passa alla creazione
della tabella. La sua dimensione finale sarà pari al prodotto del numero di valori che andiamo ad inserire nel set-up iniziale: N P SR table = n cV alues · n zV alues · n tV alues · n pV alues (3.34) L’aggiunta della segregation porterà ad ingrandire la tabella in base a quanti valori di seg g i 00 2 e seg c 00 f i 2 andremo a tabulare: N seg table = N P SR table · n z 2 V alues · n c 2 V alues (3.35) Gli andamenti delle grandezze d’interesse sono ottenuti dai calcoli come funzione del tempo; prima di costruire la tabella ognuno di essi andrà trasformato in una funzione della progress variable, in modo da poterla usare come coordinata di riferimento al momento di ricavare i dati utili dalla tabella. Non vengono utilizzati i valori assoluti di c, ma si è preferito tabularne i valori normalizzati c k norm = c k c M ax .
Per ogni combinazione di valori in input, la tabella conterrà: • c ˙ , valore del termine sorgente della progress variable. Per l’iterazione k − 1 viene calcolato come: c ˙ k − 1 = c k − c k − 1 t k − t k − 1 (3.36) Vi sono condizioni per cui è imposto nullo, come ad esempio se ci troviamo al termine della combustione, forzata o naturale, o se ci troviamo agli estremi dei valori di Z ( 0 e 1 ); • c M ax e c M in , definiti rispettivamente come il valore massimo e minimo della progress variable;
38 Capitolo 3. Generazione della tabella
• Virtual Species, frazioni massiche di tutte le specie virtuali che abbiamo
• Output Species, frazioni massiche di specie reali utili per post processing come ad esempio OH, CH 2 O e CO . Vengono scelte dall’utente al momento del set up della generazione della tabella; • t c , corrisponde agli istanti temporali per ogni valore della progress variable.
3.5.1
Correzione della tabulazione
Dato che la tabella che andremo a generare sarà utilizzata per la simulazione della combustione Diesel, le condizioni iniziali dei
corrisponderanno a punti operativi tipici dei motori
CI . Questo comporta un range di condizioni operative molto vario,
soprattutto in termini di Z, passando da miscele molto magre a molto ricche. Vi saranno quindi dei
che, in quanto al di sotto del limite di infiammabilità o molto lenti ad accendersi, porterebbero dei problemi alla costruzione della tabella.
Per questa ragione è stata inserita la possibilità di forzare l’accensione dei reattori dopo un certo intervallo di tempo, tipicamente fissato a metà del tempo massimo per il calcolo del singolo
PSR . L’entalpia dei prodotti è calcolata come funzione
di composizione, pressione e T ai corrispondente alla temperatura calcolata dopo l’accensione forzata: T ai = T q = T q − 1 + ( t q − t q − 1 ) ( T eq − T ∆ t ign q − 1 ) (3.37) T q è la temperatura alla q-esima iterazione, mentre T eq è la temperatura della miscela all’equilibrio, ossia la temperatura adiabatica di fiamma ottenuta bruciando tutto il combustibile o consumando tutto l’ O 2 presente nel
un controllo per evitare di forzare dei reattori che, seppur lentamente, si stanno accendendo. Viene infatti monitorato il rapporto: R T T q − T 0 = T eq − T 0 (3.38) in modo che in caso di superamento di una determinata soglia, impostata preventi vamente, si escluda la possibilità di forzare l’accensione.
Un altro problema, riscontrato in alcuni casi, è che la progress variable, una volta raggiunto il valore massimo, diminuisca. Ciò crea gravi problemi al momento della scrittura della tabella, dato che allo stesso valore di c corrispondono più valori delle altre grandezze. Per risolvere questo problema è stata introdotta la possibilità
3.5. Costruzione della tabella 39
di bloccare il valore della progress variable a quello corrispondente all’istante in cui si ha la massima temperatura.
3.5.2
Virtual Species
Il numero delle specie coinvolte nelle reazioni di combustione è molto elevato, per ridurre il numero delle equazioni di trasporto da risolvere e in questo caso anche per ridurre la dimensione della tabella, si è deciso di utilizzare le virtual species.
Questo approccio, proposto per l’applicazione nel modello
scelta di un set di specie che siano rappresentative dell’intero sistema, imponendo che vengano conservate alcune proprietà fondamentali della miscela. Nel nostro caso è stato scelto un sistema a 5 specie: CO 2 , H 2 O, N 2 , O 2 e f uel ; per determinarne le frazioni molari andremo a risolvere le equazioni di conservazione degli atomi e dell’entalpia assoluta della miscela: P n vs P i =1 n vs i =1 P n vs i =1 P n vs P i =1 n vs i =1 x v,i n C,i x v,i n H,i x v,i n O,i x v,i n N,i = = = = P n i =1 P n s P s i =1 n i =1 P n s s i =1 x v,i w i H i ( T, p ) = x i n C,i x x x i i i n n n P H,i O,i N,i n s i =1 x i w i H i ( T, p ) (3.39) In cui n j,i indica il numero di atomi di tipo j contenuti nella specie i − esima e w i H i ( T, p ) l’entalpia assoluta della specie i − esima espressa in J/kmol . Nota la composizione virtuale molare, le frazioni massiche verranno calcolate come: Y v,i = x v,i w i (3.40)
3.5.3
Convergenza automatica
A causa della grande variabilità delle condizioni operative dei
una grande diversità di durata dei processi di combustione. Per essere sicuri che tutti i reattori arrivino a combustione completa, sarebbe necessario imporre un end time pari a quello necessario al
più lento per arrivare all’equilibrio. Ciò implicherebbe che la maggior parte delle simulazioni abbiano una durata superiore a quella necessaria, con conseguente incremento del tempo utile alla generazione della tabella. Per velocizzare il processo di calcolo è stata sviluppata una procedura, che interrompa la simulazione una volta soddisfatto un determinato vincolo. Dall’analisi 2 Representative Interactive Flamelet
40 Capitolo 3. Generazione della tabella
dell’evoluzione delle specie, si è visto come il rapporto R conv = Y CO Y CO 2 abbia un andamento tale da permettere di identificare il raggiungimento delle condizioni di equilibrio.
(a) (b) Figura 3.3: Andamento della temperatura (curva rossa) e del rapporto Y CO /Y CO 2 verde) per un 0 .
0602965 (a) e
T 0 di condizioni iniziali: = 800 K, P = 50 bar , Z T 0 = 800 K, P = 0 .
0333968 (b) = 50 bar , (curva Z = Arrivati al termine delle reazioni il valore del rapporto in esame tende a rimanere costante, a meno di una tolleranza fissata a priori, di default pari a 10 − 6 . In questo modo è stato possibile diminuire significativamente i tempi di calcolo, per una tabella di 118422 elementi: Convergenza automatica CPU time [ s ] Clock time [ s ] ON OF F speed-up factor 700 .
9 2823 .
09 4 .
03 703 2829 4 .
02 Tabella 3.2: Confronto tra i tempi di di solo calcolo (cpu time) e dell’intera simulazione (clock time) con e senza meccanismi di velocizzazione.
Capitolo 4 Il solver FlameletDieselFoam
All’interno della libreria Lib-ICE, parallelamente al generatore della tabella per
la cinetica chimica (Cap. 3), è stato implementato il solver in grado di risolvere
il dominio termo-chimico accoppiando la risoluzione delle equazioni di conserva zione con l’estrapolazione dei risultati dalla tabella precedentemente generata. Di seguito verrà descritta la formulazione delle equazioni e successivamente la libreria contenente tutte le funzioni utilizzate per estrapolare i dati necessari dalla tabella.
4.1
Equazioni del modello
Per prima cosa il codice richiama la libreria nella quale è implementato il solver lagrangiano che, attraverso vari modelli di atomizzazione, breakup ed evaporazione, determina la struttura dello spray. Da questa vengono quindi ricavate informazioni fondamentali per la risoluzione delle equazioni di conservazione. Il sistema di equazioni utilizzato dal solutore in esame è costituito da una rielaborazione di
quello presentate nella sezione 1.5.1.
Per imporre la conservazione della massa, viene risolta l’equazione della densità: ∂ρ + ∇ · ( ρ u ) = ˙ evap ∂t (4.1) dove S ˙ evap è il termine sorgente legato all’evaporazione del combustibile, ottenuto dal modello dello spray, espresso in Kg/m 3 s .
All’interno dell’algoritmo PIMPLE viene risolta l’Equazione del momentum: ∂ρ u + ∇ · ( ρ uu ) + ˙ u ,turb ∂t = ρ g + ˙ u ,spray (4.2) S ˙ u ,spray è il termine sorgente legato allo scambio di quantità di moto con lo spray,
41
42 Capitolo 4. Il solver FlameletDieselFoam
mentre S ˙ u ,turb è legato allo scambio di quantità di moto dovuto alla turbolenza. In particolare nel caso di utilizzo del modello standard k ε , quest’ultimo viene definito come: S ˙ u ,turb = −∇ 2 ( µ ef f u ) − ∇ · ( µ ef f ) T ∇ u − 2 3 I tr ( T ∇ u ) (4.3) Per ottenere tutte le grandezze necessarie ad accedere alla tabella, vanno aggiunte le equazioni di trasporto per la progress variable e mixture fraction. A seconda della definizione che ne abbiamo dato, cambia la formulazione dell’equazione di c e in particolare il calcolo del termine sorgente legato allo spray: • Specie, S ˙ c,spray calcolato come: è legato all’evaporazione dello spray di combustibile, sarà S ˙ c, spray = N s, pv X α i S ˙ i, spray i =0 (4.4) dove N s, pv è il numero di specie di cui composta la progress variable e α i il coefficiente dell’i-esima specie. L’equazione di c risulta quindi: ∂ρc + ∇ · ( ρ u c ) − ∇ 2 ∂t · [( µ + µ t Sc t ) c S c + ˙ c, spray (4.5) Dove S ˙ c dipende dalla chimica e viene estrapolato dalla tabella generata in precedenza.
• Entalpia sensibile, andremo a considerare sia lo scambio termico con lo spray, che l’evaporazione del combustibile: S ˙ c, spray = ˙ ht, spray (4.6) L’equazione della progress variable cambia leggermente rispetto al caso pre cedente dato che il termine diffusivo non sarà più legato alla massa, ma all’energia termica: ∂ρc + ∇ · ( ρ u c ) − ∇ 2 · ( µ c c ) = ˙ c + ˙ c, spray ∂t (4.7) con µ c = α ef f = µ t /P r , diffusività termica efficace, estratta dal modello di turbolenza.
• Entalpia chimica, viene considerato lo scambio di massa con lo spray: S ˙ c, spray = N s X S ˙ i, spray h 298 , i i =0 (4.8)
4.1. Equazioni del modello 43
L’equazione della progress variable ha quindi la stessa formulazione presentata nel caso di c basata sulle specie: ∂ρc + ∇ · ( ρ u c ) − ∇ 2 ∂t · [( µ + µ t Sc t ) c S c + ˙ c, spray (4.9) In tutte le formulazioni precedenti dell’equazione della progress variable è presente un termine sorgente, S ˙ c , legato alle reazioni chimiche. Il loro valore non viene calcolato, ma estratto dalla tabella precedentemente costruita, secondo la
metodologie presentata nel paragrafo 4.2.
S ˙ c = f (∆ t, Z, seg ” g 2 , c norm , seg c ” f 2 , χ st , T u , p, S ˙ h ) (4.10) L’equazione per la mixture fraction risulta: ∂ρZ ∂t + ∇ · ( ρ u Z ) − ∇ 2 · [( µ + µ t Sc t ) Z S Z (4.11) con S ˙ Z termine sorgente legato all’evaporazione del combustibile.
Se andiamo a considerare anche l’effetto dell’interazione della turbolenza con la chimica sono necessarie due ulteriori equazioni di trasporto: • L’equazione per c 00 f 2 , ∂ρ c 00 f 2 ∂t + ∇ · ( ρ u c 00 f 2 ) −∇ 2 2 µ · Sc [( ef f 00 µ g 2 + |∇ µ t Sc 00 g 2 ) 00 g 00 2 g + c norm | 2 − C χ ρ ε k ∇ c 00 f 2 (4.12) • L’equazione da risolvere per la 00 g 2 , ∂ρ 00 g 2 + ∇ · ∂t ( ρ u µ ef f 2 Sc 00 g 2 00 g |∇ 2 ) Z | − ∇ 2 2 − · C χ ρ µ t µ + Sc 00 g 2 ε k ∇ Z 00 2 00 g 2 = (4.13) Nelle equazioni fin qui descritte, sono presenti alcune costanti legate all’interazio ne tra la chimica e la turbolenza. Queste vengono inizializzate , per le simulazioni
presentate nei capitoli successivi, come mostrato nella tabella successiva (Tab. 4.1).
Una grandezza fondamentale nella modellazione dell’interazione tra la chimica e la turbolenza, è lo scalar dissipation rate χ . É definito come l’inverso di un tempo [ 1 s ] e svolge un ruolo analogo alla diffusività. All’interno del codice viene modellato,
44 Capitolo 4. Il solver FlameletDieselFoam
Costante Valore C χ Sc Sc t Sc 00 g 2 0 0 0 .
.
.
2 7 7 7 Tabella 4.1: Valori delle costanti di chiusura utilizzate nelle equazioni di 00 g 2 e 00 2 g .
attraverso un’equazione algebrica per la mixture fraction: χ Z = 2 µ ef f |∇ Z | 2 (4.14) Analogamente, anche per la progress variable possiamo calcolare: χ C = C χ ε k ∇ c 00 f 2 (4.15) Da questi valori possiamo ricavare il valore dello scalar dissipation rate in corri spondenza della mixture fraction stechiometrica χ st necessario per il calcolo del S ˙ c
Sono infine risolte due equazioni di conservazione dell’energia, la prima è espressa in funzione dell’entalpia totale: ∂ρh ∂h + ∇ · ( ρ u h ) − ∇ 2 · ( α ef f h ) = Dp + ˙ h, spray Dt (4.16) in cui S ˙ h, spray è il termine sorgente legato allo scambio termico con lo spray. La seconda in funzione dell’entalpia della frazione di gas che non è ancora reagito: ∂ρh u ∂h u + ∇ · ( ρ u h u ) − ∇ 2 · ( α ef f h u ) = ρ ρ u Dp Dt + ρ ρ u S ˙ h u , spray (4.17) Nel caso in cui la tabella sia stata generata assumendo una mixing line, quest’ultima equazione cambia leggermente: ∂ρh u ∂h u + ∇ · ( ρ u h u ) − ∇ 2 · ( α ef f h u ) = ρ ρ u Dp + ˙ evap h u Dt (4.18) A partire dall’ h e dall’ h u vengono ricavate la temperatura della fase gas e la temperatura della miscela incombusta T u parametro necessario all’estrapolazione dei dati dalla tabella.
4.2. Flamelet combustion library 45
4.2
Flamelet combustion library
La flamelet combustion library contiene tutte le funzioni necessarie ad accoppiare la generazione della tabella alla risoluzione del dominio fluidodinamico. É possibile scegliere tra due alternative, definite allo stesso modo, con l’unica differenza che venga incluso o no il modello β
-PDF, introdotto nel paragrafo 3.3.
La funzione fondamentale, contenuta in questa libreria, è quella che permette di calcolare il termine sorgente per la progress variable. Entrando più nel dettaglio, considerando l’implementazione nella libreria che include il modello β -PDF, vi sono due approcci per il calcolo del c ˙ : • basato su t c
(Sez. 3.5), viene calcolato come:
c ˙ i = c + − c − t + c − t − c (4.19) in cui, dato il valore di c i e delle altre variabili di input, attraverso una funzioni appositamente implementate, vengono determinati i valori di progress varivable precedente ( c − ) e successiva ( c + ) e i corrispondenti valori di t − c e t + c , • estrapolato direttamente dai valori di cDot
Per migliorare l’accuratezza con cui viene calcolato il c ˙ è stata inserita l’opportunità di suddividere l’intervallo temporale su cui si effettua il calcolo. Questo risulta utile soprattutto in corrispondenza dei transitori, accensione e completamento della combustione, in cui il passo di integrazione potrebbe risultare troppo grande per valutare correttamente l’andamento delle grandezze d’interesse.
Attraverso altre funzioni della libreria in esame è possibile ricavare altre infor mazioni utili dalla tabella. Una volta calcolato il c ˙ , vengono aggiornati i valori delle frazioni massiche delle virtual species e delle output species considerate. Op pure vi sono funzioni specifiche per ricavare i valori di c M ax e c M in , utili per la normalizzazione della progress variable.
4.2.1
Tecniche di interpolazione
Accanto alle funzioni fin qui descritte, richiamate direttamente dal codice, ne sono presenti altre che vengono usate all’interno della libreria, per la ricerca e l’interpolazione dei dati nella tabella. Ne esistono di vari tipi, applicabili sia solo ad uno scalare, che anche un campo di scalari. L’implementazione di queste funzioni è analoga, cambia solamente che in un caso è possibile ottenere in output più di un valore interpolato contemporaneamente. Entrando maggiormente nel dettaglio,
46 Capitolo 4. Il solver FlameletDieselFoam
vediamo che il valore ricercato di una generica variabile tabulata φ , scalare, è espresso come: φ interp = f ( i c 2+ , i c 2 − , i Z 2+ , i Z 2 − ,i c + , i c − , i Z + , i Z − , i T + , i T − , i p + , i p − , c 2 i , Z i 2 , c i , Z i , T i , p i ) (4.20) in cui le variabili con pedici + e -, sono gli indici delle coordinate tabulate più vicine al valore i-esimo. Il risultato finale è dato dalla media pesata sulla distanza dai valori tabulati, secondo la formula: φ interp = P N index k ( φ k w k ) V tot Il volume totale è calcolato come: V tot = N coord Y | φ j + − φ j − | j (4.21) (4.22) mentre i singoli pesi sulle coordinate, w k = N coord Y | φ i − φ j,oppIndex | j (4.23)
Capitolo 5 Verifica della consistenza del modello
Per una prima verifica del corretto funzionamento del modello, descritto nei capitoli precedenti, si è scelto di svolgere delle prove di autoaccensione a pressione e volume costante. Per le simulazioni in esame vengono utilizzati due solutori, P SRF oam e tableF oam sviluppati appositamente.
Entrambi per prima cosa generano on-the-fly una mesh formata da un unica cella e successivamente risolvono le equazioni di conservazione. I due solver si distinguono per l’approccio alla cinetica chimica in quanto il P SRF oam integra direttamente le equazioni per le specie
tableF oam estrapola l’ ω ˙ c dalla tabella, come mostrato per il f lameletDieselF oam
Per la generazione della tabella e per il solver P SRF oam , l’integrazione del si stema di equazioni differenziali è stato effettuata utilizzando il solver ODE
che garantisce stabilità e accuratezza anche nel caso di sistemi stiff. La risoluzione diretta di schemi cinetici complessi comporta tuttavia oneri computazionali molto elevati a causa dell’elevato numero di sotto-iterazioni necessarie a soddisfare le tolleranze imposte e al calcolo di jacobiani di grandi dimensioni. Per questa ragione viene usato un metodo di riduzione degli schemi cinetici, il
in quanto garantisce maggiore accuratezza e stabilità [12] e [13].
Le tabelle, qui di seguito utilizzate, sono state generate utilizzando come com-
bustibile l’iso-ottano, seguendo lo schema cinetico di Pitsch [22], composto da 109
specie e 393 reazioni. Il set-up della tabella è descritto nella tabella 5.1.
1 Stiff linearly implicit EULer EXtrapolation 2 Direct Relation Graph
47
48 Capitolo 5. Verifica della consistenza del modello
Grandezza T p Z Range 700 − 1200 K 10 − 100 bar 0 − 0 .
032 , 1 Intervallo 25 K 15 bar 0 .
04 Tabella 5.1: Set-up della tabella per l’iso-ottano usata per le simulazioni PSR. La temperatura ha un intervallo di tabulazione che non è costante, oltre i 100 K si utilizza un ∆ T di 50 K.
La formulazione usata per la progress variable è quella basata sull’entalpia chimica, dati i migliori risultati ottenuti, soprattutto nelle simulazioni Spray, che
verranno mostrate successivamente (Cap. 6). Molto importante è anche la scelta
della quantità di punti da tabulare per la progress variable. É stato scelto di usare un intervallo ∆ c = 0 .
01 , con una raffinazione molto marcata per i valori vicini a 0 e 1. Operando in questo modo si è in grado di avere una maggiore accuratezza nelle fasi di accensione e completamento della combustione.
5.1
Prove a pressione costante
Con le prove a pressione costante è possibile verificare il corretto funzionamento delle sole tecniche ti tabulazione e interpolazione dei dati, in quanto la tabella stessa è costruita a partire da simulazioni PSR a pressione fissata.
Punto Temperatura [K] Pressione [bar] Mixture fraction [-] T ab Z 1 Z 2 T 1 T 2 p 1 p 2 800 800 800 812 .
5 1025 800 800 40 40 40 40 40 50 80 0 0 0 .
.
.
026 018 024 0 .
026 0 .
026 0 .
026 0 .
026 Tabella 5.2: Punti operativi analizzati per le simulazioni
costante.
a pressione e volume Come prima verifica sono stati simulati i punti operativi tabulati, verranno quindi analizzati separatamente gli effetti della variazione delle singole grandezze tabulate. In questo modo è possibile valutare la qualità dei risultati nelle condizioni di massimo errore, corrispondenti ai valori intermedi tra quelli tabulati. Con questa
5.1. Prove a pressione costante 49
metodologia è possibile effettuare delle prove di sensitività per verificare che le scelte dei set-up della tabella siano corrette.
Si può vedere, dalla figura (5.1), che i risultati ricavati dalla tabulazione delle
chimica sono perfettamente sovrapponibili a quelli che ne prevedono un’integrazione diretta. Ciò dimostra una corretta implementazione delle tecniche di costruzione della tabella.
Figura 5.1: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 800 K e Z = 0.024
5.1.1
Effetto della Z
La prima grandezza analizzata è la mixture fraction. É stato scelto, per rendere Figura 5.2: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 800 K e Z = 0.026
50 Capitolo 5. Verifica della consistenza del modello
più significativo lo studio, di andare ad analizzare punti operativi simili a quelli in cui si trova a funzionare un motore
Figura 5.3: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 800 K e Z = 0.018
Dai risultati ottenuti possiamo vedere che, in entrambi i casi, le due fasi del processo di accensione vengono descritte correttamente. Solo nel caso più magro
(Fig. 5.3) si ha un leggero divario tra le due curve nella fase di completamento della
combustione.
5.1.2
Effetto della T
La temperatura ha invece un effetto più marcato, vediamo infatti che a bassa temperatura l’errore massimo commesso è rilevante. Definendo l’ignition delay come l’istante in cui la derivata della temperatura è massima, t ID = dT dt max (5.1) la tabulazione risulta in anticipo t ID,tab = 0 .
03863 s rispetto all’integrazione diretta t ID,int dir = 0 .
043 s
. Osservando i risultati ad alta temperatura (Fig. 5.5) si nota
come si sia annullato l’anticipo della tabulazione. Viene quindi giustificata la scelta di incrementare l’intervallo di tabulazione per la temperatura superati i 1000 K
5.1. Prove a pressione costante 51
Figura 5.4: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 812.5 K e Z = 0.024
Figura 5.5: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 40 bar, con temperatura di 1025 K e Z = 0.024
5.1.3
Effetto della p
L’ultima grandezza analizzata è la pressione e si vede che, in linea con quanto osservato per la Temperatura, l’anticipo riscontrato nei risultati della tabulazione è inversamente proporzionale ad essa. In questo caso l’intervallo con cui è tabulata la pressione è piuttosto piccolo e anche l’anticipo misurato risulta più contenuto di
quello ottenuto precedentemente (Sez. 5.1.2). Nel caso della simulazione mostrata
nella figura sottostante (Fig. 5.6), il ritardo di accensione per la tabulazione è
t ID,tab = 0 .
028964 s mentre per l’integrazione diretta t ID,int dir = 0 .
030311 s
52 Capitolo 5. Verifica della consistenza del modello
Figura 5.6: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 50 bar, con temperatura di 800 K e Z = 0.024
Figura 5.7: Confronto tra chimica dettagliata e tabulata, in un PSR a pressione costante pari a 80 bar, con temperatura di 800 K e Z = 0.024
5.2
Prove a volume costante
Con le prove a volume costante, andiamo a valutare il comportamento del modello in condizioni di autoaccensione più simili a quelle di un motore.
Rispetto alle simulazioni condotte a pressione costante, vediamo come, per le medesime condizioni iniziali, l’anticipo ottenuto con la cinetica tabulata tenda ad aumentare. Questo è dovuto al fatto che la tabella è stata generata a partire da dei PSR a pressione costante. La diversa configurazione della simulazione porta ad un aumento inevitabile degli errori di interpolazione. Nonostante ciò i valori massimi di temperatura e pressione risultano correttamente predetti.
Allo stesso modo di come abbiamo operato nella sezione precedente andremo a
5.2. Prove a volume costante 53
(a) (b) Figura 5.8: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura 800 K e Z = 0.024.
valutare gli effetti delle singole grandezze tabulate.
5.2.1
Effetto della Z
In entrambi i casi presentati, variando la frazione di combustibile, il primo stadio di accensione è sempre correttamente rappresentato. L’accensione principale riscontra invece un leggero aumento di anticipo andando a smagrire la miscela.
Figura 5.9: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura 800 K e Z = 0.026
54 Capitolo 5. Verifica della consistenza del modello
Figura 5.10: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura di 800 K e Z = 0.018
5.2.2
Effetto della T
Per quanto riguarda la temperatura, l’effetto della variazione della configurazione dei
risulta meno importante. Ad alti valori (Fig. 5.12) l’anticipo della cinetica
tabulata è trascurabile. Scendendo con la temperatura (Fig. 5.11), lo sfasamento
tra chimica dettagliata e tabulata aumenta, pur rimanendo entro valori contenuti.
Nel caso di una temperatura iniziale di 812.5 K i ritardi di accensione ottenuti sono: t ID,tab = 0 .
028569 s e t ID,int dir = 0 .
0324573 s .
Figura 5.11: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura di 812.5 K e Z = 0.024
5.2. Prove a volume costante 55
Figura 5.12: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 40 bar, temperatura di 1025 K e Z = 0.024
5.2.3
Effetto della p
Come riscontrato nelle simulazioni precedenti (Sez. 5.1 ), l’effetto della pressione
è analogo a quello della temperatura. Notiamo infatti come, all’aumentare di tale grandezza, l’anticipo della cinetica tabulata vada riducendosi fino ad annullarsi,
come nel caso mostrato in figura 5.14.
Figura 5.13: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 50 bar, temperatura di 800 K e Z = 0.024
56 Capitolo 5. Verifica della consistenza del modello
Figura 5.14: Confronto tra chimica dettagliata e tabulata, in un PSR a volume costante di condizioni iniziali: pressione 80 bar, temperatura di 800 K e Z = 0.024
5.3
Calcolo degli ignition delay con l’n-dodecano
Nello studio della combustione Diesel il ritardo di accensione svolge un ruolo
fondamentale, come già spiegato nel paragrafo 1.1. In questa sezione andremo
quindi a verificare che il modello sia in grado di predire con precisione tale grandezza.
Come combustibile, per la generazione delle tabelle, è stato usato l’n-dodecano,
seguendo lo schema cinetico di Pei [21], composto da 58 specie e 272 reazioni. La
formulazione della progress variable è sempre quella basata sull’entalpia chimica di
formazione e il set-up è schematizzato nella tabella 5.3.
Grandezza T p Z Range 600 900 − 900 − 1200 K K 20 − 240 bar 0 − 0 .
2 , 1 Intervallo 25 K 50 K 20 bar ∼ 0 .
04 Tabella 5.3: Set-up della tabella per l’n-dodecano usata per le simulazioni PSR a volume costante.
Per valutare il comportamento del modello al variare delle condizioni ambiente, sono state generate tre tabelle con lo stesso set-up, ma diversa composizione dell’aria. Le condizioni scelte corrispondono ad una percentuale di EGR pari allo 0%, 25% e 35%. Diversamente dalle tabelle generate in precedenza per l’iso-ottano, l’intervallo scelto per la tabulazione della Z non è costante. Si è deciso infatti di partire da un egual frazionamento del rapporto d’equivalenza, dal quale sono stati
ricavati i valori di mixture fraction (Fig. 5.15). Questa scelta è legata alla volontà
5.3. Calcolo degli ignition delay con l’n-dodecano 57
di mantenere una distribuzione costante dei valori intorno allo stechiometrico al variare della composizione dell’aria. Il ritardo di accensione di una miscela ha infatti un andamento che presenta un minimo in corrispondenza delle condizioni stechiometriche. Risulta quindi fondamentale rappresentare adeguatamente questa zona per calcolare correttamente l’ignition delay.
Figura 5.15: Valori di mixture fraction in funzione del rapporto di equivalenza.
Per le simulazioni sono stati scelti punti operativi caratteristici di diverse condizioni di accensione: PCCI, iniezione pilota e iniezione principale; e diversa quantità di combustibile in miscela, stechiometrico ( φ = 1 ), magro ( φ = 0 .
5 ) e ricco ( φ = 1 .
5 ).
Figura 5.16: Punti operativi scelti per il calcolo dell’ignition delay.
58 Capitolo 5. Verifica della consistenza del modello
(a) Miscela magra, φ =0.5
.
(b) Miscela stechiometrica, φ = 1 .
(c) Miscela ricca, φ = 1.5
.
Figura 5.17:
Calcolo degli ignition delay dei punti operativi mostrati in figura 5.16, per
un EGR del 35 % Dai risultati ottenuti si vede come il modello sia in grado di riprodurre con buona accuratezza gli istanti di accensione per un ampio range di condizioni operative. Si può notare anche che il suo comportamento, osservato per le prove con l’iso-ottano, sia confermato. Vediamo infatti che la tabulazione, in caso di errore, è sempre in anticipo rispetto all’integrazione diretta. In accordo con i risultati precedenti è anche la riduzione dell’anticipo all’aumentare di pressione e temperatura.
5.3. Calcolo degli ignition delay con l’n-dodecano 59
(a) Miscela magra, φ =0.5 senza EGR .
(b) Miscela magra, φ =0.5 e EGR del 25% .
(c) Miscela stechiometrica, φ = 1 senza EGR .
(d) Miscela stechiometrica, 25% .
φ = 1 e EGR del (e) Miscela ricca, φ = 1.5 senza EGR .
(f ) Miscela ricca, φ = 1.5 e EGR del 25% .
Figura 5.18:
Calcolo degli ignition delay dei punti operativi mostrati in figura 5.16,
senza EGR (colonna di sinistra) e con un EGR del 25% (colonna di destra)
Capitolo 6 Validazione
Si procede ora, in questo capitolo, alla validazione del modello con due configu razioni di fiamma più complesse. Per prima cosa si andrà a simulare la combustione all’interno di un motore in assetto
HCCI . Successivamente verrà svolta un’ampia
casistica di simulazioni spray A.
6.1
Simulazione della combustione
HCCI
6.1.1
Condizioni termodinamiche
Il motore che si andrà a simulare nel seguito è un classico esempio di Diesel
medium-duty, descritto esaustivamente in [7]. Per tutti i casi analizzati la velocità
di rotazione è 1200 rpm e la carica, di aria e iso-ottano, è perfettamente miscelata ottenendo una combustione di tipo
e i rapporti di equivalenza possibili variano da φ = 0 .
08 a φ = 0 .
028 . Per i risultati descritti nelle prossime pagine sono stati scelti tre punti di funzionamento: φ = 0 .
14 , φ = 0 .
20 e φ = 0 .
28 .
φ Z T cil p cil 0 .
14 0 .
00884 476 .
30 K 1 .
361 bar 0 .
20 0 .
01265 474 .
10 K 1 .
378 bar 0 .
28 0 .
01774 468 .
30 K 1 .
357 bar Tabella 6.1: Condizioni operative per le simulazioni
61
62 Capitolo 6. Validazione
6.1.2
Mesh
Vista la natura assialsimmetrica del problema, la geometria del cilindro è riprodotta nel dettaglio all’interno di una mesh bidimensionale avente 3207 celle al
e 815 al
(a) (b) Figura 6.1: Mesh 2-D del motore in esame in vista laterale (a) e dall’alto (b).
6.1.3
Modelli utilizzati
Le simulazioni in esame sono state condotte utilizzando due diversi modelli
per la combustione. Il primo è il Well-mixed con multi-zone chemistry ( MZWM 3
proposto da [1]. Questo modello prevede la creazione di insiemi di celle con simili
temperatura e rapporto di equivalenza, a meno di una tolleranza imposta dall’utente.
Le equazioni per la cinetica chimica vengono quindi risolte una sola volta per gruppo e i reaction rates così calcolati vengono poi utilizzati nel dominio
questo modello è presente all’interno della Lib-ICE e ampiamente validato in [17]. Il
secondo, che prevede la tabulazione della cinetica chimica, è quello di cui è oggetto questa tesi. I set-up per la generazione della tabella e lo schema cinetico usato, sono quelli già presentati per le simulazioni
a pressione e volume costante (Sez.
Per la turbolenza è stato usato lo standard k − ε , con funzione a parete di
Angelberger [18] per lo scambio termico.
1 Bottom Dead Centre 2 Top Dead Centre 3 Multi-Zone Well-Mixed
6.1. Simulazione della combustione
63
6.1.4
Analisi dei risultati
Si procede ora col mostrare gli andamenti delle grandezze di interesse per un range di angoli di manovella che comprende interamente il processo di combustione ( − 30 ≤ ϑ ≤ 30 ). Il confronto verrà fatto tra i due modelli di combustione e, quando presenti, con i dati sperimentali. Analizzando il primo punto operativo, φ = 0 .
14 , (a) (b) (c) Figura 6.2: Andamento di pressione,
e temperatura per la simulazione HCCI con rapporto di equivalenza φ = 0 .
14 possiamo vedere come i risultati ottenuti dalla tabulazione della cinetica chimica siano perfettamente consistenti con quelli del modello Well-mixed. Si nota anche un ottimo accordo con i dati sperimentali, sia in termini di istante di accensione
che di valore massimo di pressione raggiunto (Fig. 6.2a). Dall’andamento del
vediamo che non è presente alcun anticipo da parte della tabulazione. Anche la
64 Capitolo 6. Validazione
curva di pressione è correttamente calcolata, con una leggera sovrastima del suo valore massimo, rispetto all’integrazione diretta della chimica.
(a) (b) (c) Figura 6.3: Andamtento di pressione,
e temperatura per la simulazione HCCI con rapporto di equivalenza φ = 0 .
20 Anche per il punto operativo successivo, φ = 0 .
20 , notiamo un buon accordo tra i due modelli di combustione. Vengono predetti gli stessi valori massimi raggiunti di temperatura e pressione, che lo stesso istante di accensione. Rispetto alla curva
sperimentale (Fig. 6.3a) notiamo invece un leggero anticipo da parte di entrambi i
modelli. Per ridurre questo tipo di errore potrebbe essere opportuno migliorare il set-up CFD della simulazione o migliorare la qualità della mesh procedendo con una raffinazione. Anche per il caso con φ = 0 .
28 vengono ottenuti risultati simili a quelli dei casi precedenti. Si può notare un leggero anticipo della chimica tabulata rispetto al modello Well-Mixed
6.1. Simulazione della combustione
65
(a) (b) (c) Figura 6.4: Andamtento di pressione,
e temperatura per la simulazione HCCI con rapporto di equivalenza φ = 0 .
28 Si va ora analizzare più nel dettaglio il caso φ = 0 .
020 . Per prima cosa viene mostrata la distribuzione della temperatura all’interno del cilindro per alcuni istanti
nell’intorno del processo di combustione. Dalle immagini riportate (Fig 6.5) si nota
come la tabulazione della cinetica chimica sia in grado di rappresentare nel dettaglio il processo di combustione, in accordo col modello well mixed. Si vede come la parte centrale del cilindro rimane sempre la zona più calda, mentre le regioni periferiche sono significativamente più fredde. Molto importante, a questo proposito, è la rappresentazione della regione compresa tra le pareti del cilindro e il pistone. Queste zone rimangono più fredde a causa del più alto scambio termico verso l’esterno, le reazioni avvengono quindi più lentamente, contribuendo in maniera significativa al totale delle emissioni del motore.
66 Capitolo 6. Validazione
(a) Temperatura K (b) -10 CAD (c) -5 CAD (d) 0 CAD (e) 5 CAD (f ) 10 CAD Figura 6.5: Distribuzione dei valori di temperatura all’interno del cilindro nel caso con φ = 0 .
20 . Nella colonna sinistra sono presenti i valori calcolati col modello MZWM, a destra con la tabulazione della cinetica chimica.
A tal proposito risulta significativo a guardare ad esempio la distribuzione del
monossido di carbonio (Fig. 6.6). Si vede che, come riscontrato in [7], il CO si
formi inizialmente nella regione centrale, per poi spostarsi verso le fessure laterali, dove vi arriva intorno 5 ◦ dopo il
6.1. Simulazione della combustione
(a) Y CO
67
(b) -10 CAD (c) -5 CAD (d) 0 CAD (e) 5 CAD (f ) 10 CAD Figura 6.6: Distribuzione dei valori di monossido di carbonio all’interno del cilindro nel caso con φ = 0 .
20 . Nella colonna sinistra sono presenti i valori calcolati col modello MZWM, a destra con la tabulazione della cinetica chimica.
68 Capitolo 6. Validazione
6.2
Campagna sperimentale Spray A
Quest’ultima validazione è stata effettuata sfruttando i dati sperimentali messi a disposizione dall’Engine Combustione Network , costituito da varie istituzioni accademiche con l’intento di condividere e coordinare numerose sperimentazioni nell’ambito dei motori a combustione interna. Le prove sperimentali sono condotte in reattori a volume costante e usando combustibili simili, per proprietà fisiche
e chimiche, al Diesel stradale. La camera utilizzata (Fig. 6.7) è un volume di
forma cubica, di lato 108 mm. Su ogni lato sono presenti delle aperture circolari di diametro 105 mm da cui è possibile ottenere le informazioni necessarie attraverso rivelazioni ottiche, laser o schlieren. Le condizioni operative all’interno del vessel (a) (b) Figura 6.7: Fotografia dell’esterno (a) e rappresentazione dell’interno (b) della camera di prova Sandia.
possono essere regolate in modo da essere rappresentative di ogni punto operativo in cui un motore Diesel può trovarsi ad operare. I laboratori Sandia hanno a disposizione ad esempio una camera in cui la temperatura può variare da 450 K a 1400 K, la densità da 1 kg/m 3 a 60 kg/m 3 e la pressione fino a 350 bar. Il combustibile è alimentato usando degli iniettori common-rail con pressioni che variano da 40 a 200 MPa e diametri da 0.05 a 0.5 mm.
6.2.1
Condizioni termodinamiche
ha identificato alcune condizioni sperimentali che sono di particolare interesse, come ad esempio quelle definite dallo standard "spray A". Queste condizioni sono caratterizzate da basse temperature, tipiche di motori che usano un moderato
. A partire dalla configurazione di base (Tab. 6.3), sono
6.2. Campagna sperimentale Spray A 69
disponibili i risultati sperimentali ottenuti con la variazione di alcune grandezze, utili per la successiva fase di validazione parametrica del modello di combustione.
grandezza x O 2 T [ K ] ρ [ kg/m 3 ] p inj [ bar ] f uel valore 0 .
15 900 22 .
8 n 1500 − dodecano Tabella 6.2: Condizioni operative standard per le simulazioni Spray A.
6.2.2
Mesh
Vista la geometria della camera di combustione utilizzata, la mesh è anche in questo caso bidimensionale e assialsimmetrica. La griglia, rappresentata in figura
6.8, è composta da 23328 (108x216) celle, con un pesante infittimento nella zona
dell’iniettore, caratterizzato da un fattore di scalatura f scal = ∆ max ∆ min = 10 .
(a) (b) Figura 6.8: Mesh 2-D della camera di combustione in esame,vista laterale (a) e dall’alto (b).
6.2.3
Iniezione di combustibile
Lo spray di combustibile è stato risolto utilizzando un modello lagrangiano
implementato all’interno della Lib-ICE e ampiamente validato in [14]. Lo spray è
70 Capitolo 6. Validazione
quindi modellato attraverso un numero finito di particelle, ognuna rappresentante un gruppo di gocce di combustibile accomunate dalle stesse proprietà fisiche. La soluzione delle equazioni lagrangiane permette di studiare la traiettoria di tali particelle fino alla loro completa evaporazione. Come mostrato nelle equazioni del
paragrafo 4.1, l’accoppiamento col dominio euleriano, non si limita unicamente allo
scambio di massa, da fase liquida a vapore, ma anche a quello di quantità di moto ed energia. Le equazioni lagrangiane non sono però in grado di descrivere tutti i fenomeni che si susseguono durante la formazione, penetrazione ed evaporazione dello spray. Sono quindi necessari dei sottomodelli fenomenologici, come ad esempio il
per il breackup secondario, utilizzato nelle successive simulazioni.
fornisce, parallelamente ai risultati sperimentali, le specifiche tecniche degli iniettori utilizzati: data la legge di iniezione è possibile quindi andare a simulare la struttura dello spray. I casi che saranno simulati in seguito hanno una
legge di iniezione, mostrata in figura 6.9a, caratterizzata da rampe di inizio e fine
iniezione molto ripide. Questo genera una forte interazione tra l’aria, inizialmente in quiete, e il combustibile. Nonostante ciò il modello per lo spray risulta accurato
nella determinazione della penetrazione del liquido e del vapore (Fig. 6.9b).
(a) (b) Figura 6.9:
(a)Legge di iniezione calcolata fornita dall’ ECN
(b) confronto tra la penetrazione,calcolata e sperimentale, di liquido e vapore.
6.2.4
Modelli utilizzati
Per le successive simulazioni sono state utilizzate delle tabelle con lo stesso
set-up presentato per le simulazioni a volume costante (Tab. 6.1). Verrà fatta
però un’analisi sull’influenza del meccanismo cinetico. Per la combustione dell’n 4 Kelvin-Helmotz-Rayleigh-Taylor
6.2. Campagna sperimentale Spray A 71
dodecano sono disponibili diversi schemi ridotti, tra questi sono stati messi a
Per la turbolenza è stato usato il modello standard k − epsilon , con costante di chiusura C 1
6.2.5
Analisi dei risultati
Effetto del meccanismo cinetico
L’analisi di sensitività al meccanismo cinetico è stata effettuata usando le condi zioni operative standard per lo spray A, solo la temperatura è stata modificata, portandola ad 800 K, in modo che la fase di accensione sia più lenta e maggiormente comprensibile. L’influenza dello schema di reazione è particolarmente evidente sulla curva di rilascio del calore e sul ritardo di accensione, funzioni del grado di avanzamento delle reazioni chimiche e quindi della reattività dei meccanismi cinetici.
Figura 6.10: Andamento delle curve di rilascio del calore sperimentale e calcolate con diversi scemi cinetici
In figura 6.10 è mostrata la curva di rilascio del calore sperimentale e dei tre
schemi cinetici considerati . Si può vedere che la curva di rilascio del calore in fase diffusiva risulta ben calcolata da tutti gli schemi mentre, per quanto riguarda le reazioni a bassa temperatura, lo schema di Pei è quello che maggiormente si avvicina all’ignition delay sperimentale, t ID,exp = 0 .
94 ms . Per tutti e tre gli schemi è però presente una sovrastima del picco di rilascio di calore e una sottostima del
72 Capitolo 6. Validazione
LOL . La sottostima del LOL è particolarmente marcata nel caso dello schema di Pei
e Cai dato che la loro elevata reattività sposta la zona di combustione più vicina all’iniettore. Visti i risultati ottenuti, per le simulazioni successive, si è quindi scelto di utilizzare lo schema di Pei.
Figura 6.11: Andamento del LOL sperimentale e calcolato con diversi scemi cinetici
Effetto della segregation della Z
Il problema della stima dell’ID potrebbe essere legato ad una mancanza della modellazione dell’interazione tra la chimica e la turbolenza. Per cercare di migliorare questo aspetto sono state svolte simulazioni utilizzando una tabella calcolata con l’approccio β
-PDF riferita alla mixture fraction, descritto nel paragrafo 3.3.
Figura 6.12: Andamento delle curve di rilascio del calore sperimentale e calcolate con e senza la segregation della mixture fraction.
6.2. Campagna sperimentale Spray A 73
Per la costruzione della tabella è stato usato lo stesso set-up delle tabelle precedenti con l’aggiunta dei valori della segregation di 00 2 g . Andando sempre ad
osservare la curva di rilascio del calore (Fig. 6.14), possiamo osservare chiaramente
l’effetto dell’introduzione dell’interazione chimica-turbolenza.
Si può notare come l’istante di accensione venga spostato più a destra, andando a sovrastimare leggermente il valore sperimentale. Questo effetto è legato alla turbolenza che, date le alte velocità con cui il combustibile esce dall’iniettore, porta ad allungare il ritardo di accensione.
Figura 6.13: Andamento delle curve di rilascio del calore sperimentale e calcolate con e senza la segregation della mixture fraction.
Figura 6.14: Andamento del LOL sperimentale e calcolato con e senza l’introduzione della segregation.
Il secondo cambiamento che si nota è l’aumento del valore massimo di rilascio del calore. Non è un effetto della turbulence-chemistry interaction, ma una mancanza
74 Capitolo 6. Validazione
di accuratezza del modello, risolta con l’introduzione della mixing line (Sez. 3.4).
Possiamo infatti vedere come il picco massimo di rilascio del calore sia di nuovo pari al caso senza la segregation, mantenendo l’aumento dell’ignition delay legato all’effetto della turbolenza sulla chimica. Utile è anche guardare l’andamento della
lunghezza del lift-off (Fig. 6.14). Si può notare come l’effetto della turbulence-
chemistry interaction sia maggiormente evidente nei primi istanti di combustione.
La zona di accensione viene infatti spostata più lontana dall’iniettore rispetto al caso di base. Tuttavia in regime di combustione diffusiva, il LOL tende a ridursi, stabilizzandosi a valori lontani dai dati sperimentali, riscontrati anche senza l’uso della β − P DF .
Passiamo ora all’analisi della struttura di fiamma calcolata, in riferimento alle
rilevazioni ottiche effettuate in bomba, per le condizioni spray A (Tab. 6.3). Il primo
confronto effettuato è quello con le rilevazioni di chemiluminescenza, effettuate
nell’intorno dell’istante di accensione. Dai fotogrammi, riportati in figura 6.15,
Figura 6.15: Chemiluminescenza rilevata a 0.8 ms, 1.2ms e 2.5 ms. La linea blu indica un livello di luminosità pari al 50% del massimo rilevabile.
possiamo dividere la fiamma in tre zone. Partendo dall’uscita dell’iniettore, fino alla lunghezza di lift-off vediamo che non è presente alcuna luminosità. Questo è dovuto al fatto che prima che si scaturisca la combustione vera e propria, devono
6.2. Campagna sperimentale Spray A 75
avvenire delle reazioni preliminari (Sez. 1.4), inoltre la turbolenza generata dal-
l’iniezione di combustibile ad alta pressione nel vessel, tende a spostare lontano dall’iniettore la zona di accensione. Successivamente è presente una zona a bassa luminosità, caratterizzata da una miscela piuttosto ricca in cui avvengono reazioni a bassa temperatura. Infine troviamo la zona di combustione ad alta temperatura, caratterizzata da una forte luminosità.
(a) (b) Figura 6.16: Confronto tra la struttura di fiamma sperimentale e calcolata, con (b) e senza (a) l’uso della β − P DF a 0.8 ms
ASOI . La linea bianca identifica i
punti a miscela stechiometrica, quella nera i punti a 1400 K.
Osservando le strutture di fiamma calcolate, vediamo come non rispecchino
esattamente quelle sperimentali. Come già visto nel grafico in figura 6.14, il lift-off
76 Capitolo 6. Validazione
risulta molto sottostimato, portando l’inizio della fiamma molto vicino all’iniettore.
La regione di alta temperatura è invece sempre correttamente posizionata, in corrispondenza di una miscela stechiometrica.
(a) (b) Figura 6.17: Confronto tra la struttura di fiamma sperimentale e calcolata, con (b) e senza (a) l’uso della β − P DF a 1.2 ms
ASOI . La linea bianca identifica i
punti a miscela stechiometrica, quella nera i punti a 1400 K.
Andiamo quindi a confrontare le strutture di fiamma calcolate con e senza l’utilizzo dell β − P DF . Si può vedere che l’effetto della turbulence-chemistry interaction introdotto non è sufficiente a rappresentare correttamente la zona di estinzione turbolenta, lasciando il
inalterato rispetto al caso senza segragation.
Si osserva invece un effetto nella zona di alta temperatura, dove si ha una maggiore
6.2. Campagna sperimentale Spray A 77
dimensione della zona di reazione. Questo effetto è particolarmente visibile dal campo di OH, che presenta una maggior estensione radiale e una distribuzione più uniforme.
(a) (b) Figura 6.18: Confronto tra la struttura di fiamma sperimentale e calcolata, con (b) e senza (a) l’uso della β − P DF a 2.5 ms
ASOI . La linea bianca identifica i
punti a miscela stechiometrica, quella nera i punti a 1400 K.
Si passa ora ad analizzare la struttura di fiamma ricavata dalle rilevazioni simultanee ottenute mediante tecnica Schlieren e
tecnica di imaging Schlieren si basa sul principio che gradienti di temperatura e composizione all’interno di un gas, possono essere visualizzati come variazioni 5 Planar Laser-Induced Fluorescence
78 Capitolo 6. Validazione
dell’indice di rifrazione. Ad esempio può essere rilevato il confine tra lo spray in evaporazione e l’aria ambiente dato che il combustibile, evaporando, raffredda l’aria che lo circonda creando un gradiente di temperatura all’interfaccia tra aria e combustibile.
(a) (b) Figura 6.19: Confronto tra il campo di formaldeide sperimentale e calcolato, con (b) e senza (a) l’uso della β − P DF a 0.3 ms
In maniera analoga possono anche essere rilevati i gradienti caratteristici delle fasi di accensione sia ad alta che a bassa temperatura. In alcuni casi i gradienti sono meno evidenti a causa della presenza di alcune specie intermedie che si formano in
6.2. Campagna sperimentale Spray A 79
fase vapore. Quindi questa tecnica risulta poco chiara. Si utilizza allora la rilevazione della formaldeide che è uno dei principali prodotti intermedi della combustione. La corretta rilevazione del CH2O fornisce quindi un indice della corretta riproduzione dell’accensione a bassa temperatura.
(a) (b) Figura 6.20: Confronto tra il campo di formaldeide sperimentale e calcolato, con (b) e senza (a) l’uso della β − P DF a 0.4 ms
Nel primo istante temporale scelto, corrispondente a circa 0.3 ms
(Fig.
6.19), possiamo vedere dal campo di temperatura che le prime reazioni a bassa
temperatura sono già cominciate. Guardando la distribuzione della formaldeide
80 Capitolo 6. Validazione
si nota come la simulazione che considera l’interazione tra chimica e turbolenza si avvicini maggiormente al dato sperimentale. Le reazioni a bassa temperatura nel caso senza segregation risultano infatti più sviluppate e la presenza di formal deide si estende maggiormente verso l’iniettore e verso il centro dello spray. In corrispondenza di circa 0.4 ms
ASOI , cominciano le reazioni ad alta temperatura,
come indicato dalla marcata dilatazione radiale ed assiale della testa dello spray.
La zona ad alta concentrazione di formaldeide è correttamente rappresentata da entrambi, al centro dello spray, ma per il modello che non utilizza la β − P DF vi è un eccessivo sviluppo assiale della distribuzione di CH2O. Entrambi i modelli non sono in grado di rappresentare la regione a bassa concentrazione di formaldeide
vicino all’iniettore, molto importante per la stabilizzazione della fiamma [25].
Validazione parametrica
Come ultima analisi sono state svolte delle simulazioni al variare di alcune grandezze, partendo dalle condizioni spray A, per investigare la resa del modello al variare delle condizioni operative in termini di
e ignition delay.
grandezza T [ K ] x O 2 ρ [ kg/m 3 ] p inj [ bar ] valori 800 , 850 , 900 , 1000 , 1100 , 1200 0 .
15 , 0 .
13 , 0 .
21 22 .
8 , 1500 , 15 1000 , .
2 500 Tabella 6.3: Condizioni operative analizzate nel corso della validazione parametrica del modello in configurazione Spray.
Figura 6.21: Confronto tra valori calcolati e sperimentali di ritardo di accensione al variare della temperatura del vessel.
6.2. Campagna sperimentale Spray A 81
Nei risultati presentati, l’ignition delay è stato calcolato in maniera consistente
rispetto ai valori forniti dall’ ECN . Vieni quindi definito come l’istante in cui si
riscontra un incremento della pressione del vessel, rispetto al valore iniziale, pari a 3 kPa.
La prima grandezza che si è scelto di analizzare è la temperatura dell’aria
all’interno del vessel. Dalla figura 6.21 si può vedere come l’ignition delay sia
sempre predetto con una buona accuratezza. Il calcolo del
è invece affetto da un errore maggiore, soprattutto per temperature più basse, in cui si ha una pesante
sottostima del valore sperimentale (Fig. 6.22).
Figura 6.22: Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della temperatura del vessel.
Andando a guardare i risultati al variare della percentuale di ossigeno nell’aria, si osserva come ancora una volta i ritardi di accensione siano correttamente calcolati
continua ad essere sottostimato (Fig. 6.24).
Figura 6.23: Confronto tra valori calcolati e sperimentali di ritardo di accensione, al variare della % di O 2 nel vessel.
82 Capitolo 6. Validazione
Figura 6.24: Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della % di O 2 nel vessel.
Il comportamento del modello è confermato anche al variare della pressione di
trend di variazione che rispecchiano quelli sperimentali, al contratio della lunghezza di lift-off che mantiene un certo errore rispetto al dato misurato.
Figura 6.25: Confronto tra valori calcolati e sperimentali di ritardo di accensione, al variare della pressione di iniezione del combustibile.
6.2. Campagna sperimentale Spray A 83
Figura 6.26: Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della pressione di iniezione del combustibile.
Figura 6.27: Confronto tra valori calcolati e sperimentali di ritardo di accensione, al variare della densità dell’aria nel vessel.
Figura 6.28: Confronto tra valori calcolati e sperimentali della lunghezza di lift-off, al variare della densità dell’aria nel vessel.
Conclusioni
In questo lavoro di tesi è stato analizzato il processo di combustione, in particolare all’interno dei motori Diesel. A questo proposito è stato sviluppato un modello basato sulla tabulazione della cinetica chimica. Per prima cosa sono state svolte delle simulazioni
per poi passare alla modellazione di problemi più complessi come la combustione in un motore HCCI e di uno spray di combustibile a volume costante. Confrontando i risultati ottenuti con le rilevazioni sperimentali presenti in letteratura si è potuto quindi investigare il comportamento del modello nelle diverse condizioni.
Dalle simulazioni
PSR , abbiamo dimostrato il corretto funzionamento delle
tecniche di creazione ed utilizzo della tabella, al variare di temperatura, pressione, mixture fraction, combustibile e frazione di ossigeno nell’aria.
Attraverso le simulazioni HCCI, abbiamo potuto osservare che la tabulazione dei
porta ad ottimi risultati in caso di combustione premiscelata. É stato riscontrato infatti un ottimo accordo sia con i dati sperimentali che con il modello well-mixed, dimezzando il tempo di calcolo necessario.
Il modello ha invece dimostrato una scarsa accuratezza nella simulazione di fiamme turbolente non premiscelate. Vengono calcolati correttamente, e in tempi molto rapidi, le tendenze in termini di ignition delay e lift-off, al variare di diffe renti condizioni operative. Anche la struttura della fiamma in regime diffusivo, è riprodotta in maniera abbastanza accurata. Si è riscontrato però una sottostima costante del valore assoluto del
β − P DF , per la mo dellazione della turbulence-chemistry interaction, non è stata in grado di risolvere tutti i problemi. Porta ad una correzione dell’ignition delay, e ad una fiamma con un’estensione radiale maggiore. Tuttavia per temperature basse, i valori di lift-off a cui si stabilizza sono ancora molto diversi da quelli rilevati sperimentalmente.
Questi problemi richiedono quindi un’analisi più approfondita del processo di accensione e di come questo sia influenzato dalla diffusione della progress variable all’interno del dominio di calcolo. A questo proposito, ulteriori sviluppi per il modello di tabulazione potrebbero prevedere l’introduzione diverse strutture di fiamme da tabulare, come ad esempio fiamme diffusive. Potrebbe essere utile anche
85
86 Conclusioni
andare a studiare nuove formulazioni di progress variable, più efficaci nel caso di combustione non premiscelata. Infine, il modello, sia con l’utilizzo della β − P DF che senza, è risultato sensibile al meccanismo cinetico utilizzato. Lo sviluppo di nuovi schemi cinetici più accurati potrebbe quindi portare risultati più precisi in termini di lift-off, ritardo di accensione e picco di rilascio del calore.
LES RANS RNG ISAT ILDM DAC FGM SI CI EGR SOI SOC LOL PM LTC HCCI PCCI RCCI DNS
Acronimi
Spark Ignition Compression Ignition Exhaust Gas Recirculation Start of Injection Start of Combustion Lift-off Lenght Particulate Matter Low Temperature Combustion Homogeneous Charge Compression Ignition Premixed Charge Compression Ignition Reactivity Controlled Compression Ignition Direct Numerical Simulation Large-Eddy simulation Reynolds-Averaged Navier-Stokes Renormalization Group Method In Situ Adaptive Tabulation Intrinsic Low-Dimensional Manifold Dynamic Adaptive Chemistry Flamelet-Generated Manifolds
87
88
TDAC PSR PDF DRG RIF BDC TDC EGR MZWM RoHR ECN KHRT SEULEX PLIF ASOI CFD OpenFOAM Tabulation of Dynamic Adaptive Chemistry Perfectly-Stirred Reactor Probability Density Function Direct Relation Graph Representative Interactive Flamelet Bottom Dead Centre Top Dead Centre Exhaust Gas Recirculation Multi-Zone Well-Mixed Rate of Heat Release Engine Combustion Network Kelvin-Helmotz-Rayleigh-Taylor Stiff linearly implicit EULer EXtrapolation Planar Laser-Induced Fluorescence After Start Of Injection Computational Fluid Dynamics Open source Field Operation And Manipulation
Acronimi
Bibliografia Riferimenti citati nel testo Pubblicazioni e Manuali
[1] Salvador M Aceves et al.
A multi-zone model for prediction of HCCI combu stion and emissions
. Rapp. tecn. SAE Technical paper, 2000 (cit. a p. 62).
[2] Liming Cai et al. «Optimized reaction mechanism rate rules for ignition of normal alkanes». In: Combustion and Flame
[3] Francesco Contino et al. «Coupling of in situ adaptive tabulation and dy namic adaptive chemistry: An effective method for solving combustion in engine simulations». In: Proceedings of the Combustion Institute 33.2 (2011),
[4] Gianluca D’Errico et al. «CFD modelling of combustion in Heavy-Duty Diesel Engines». In: THIESEL, Conference on Thermo- and Fluid Dynamic Processes in Direct Injection Engines
[5] Giancarlo Ferrari.
Motori a combustione interna . Il capitello, 2008 (cit. alle
[6] Alessio Frassoldati et al. «Reduced kinetic mechanisms of diesel fuel surrogate for engine CFD simulations». In: Combustion and Flame 162.10 (2015),
[7] Randy P Hessel et al.
Modeling iso-octane HCCI using CFD with multi-zone detailed chemistry; comparison to detailed speciation data over a range of lean equivalence ratios
. Rapp. tecn. SAE Technical Paper, 2008 (cit. alle pp. 61,
[8] Matthias Ihme, Lee Shunn e Jian Zhang. «Regularization of reaction progress variable for application to flamelet-based combustion models». In: Journal of Computational Physics
231.23 (2012), pp. 7715–7721 (cit. a p. 31).
89
90 Bibliografia
[9] S Imtenan et al. «Impact of low temperature combustion attaining strategies on diesel engine emissions for diesel and biodiesels: a review». In: Energy Conversion and Management
80 (2014), pp. 329–356 (cit. a p. 7).
[10] Harry Lehtiniemi et al.
Efficient 3-D CFD combustion modeling with transient flamelet models
. Rapp. tecn. SAE technical paper, 2008 (cit. a p. 31).
[11] Harry Lehtiniemi et al. «Modeling diesel spray ignition using detailed chemi stry with a progress variable approach». In: Combustion science and technology
178.10-11 (2006), pp. 1977–1997 (cit. a p. 31).
[12] Tianfeng Lu e Chung K Law. «A directed relation graph method for me chanism reduction». In: Proceedings of the Combustion Institute 30.1 (2005),
[13] Tianfeng Lu e Chung K Law. «Linear time reduction of large kinetic mechani sms with directed relation graph: n-Heptane and iso-octane». In: Combustion and Flame
144.1 (2006), pp. 24–36 (cit. a p. 47).
[14] Tommaso Lucchini, Gianluca D’Errico e Daniele Ettorre. «Numerical investiga tion of the spray–mesh–turbulence interactions for high-pressure, evaporating sprays at engine conditions». In: International Journal of Heat and Fluid Flow
32.1 (2011), pp. 285–297 (cit. a p. 69).
[15] Ulrich Maas e Stephen B Pope. «Implementation of simplified chemical kinetics based on intrinsic low-dimensional manifolds». In: Symposium (International) on Combustion
. Vol. 24. 1. Elsevier. 1992, pp. 103–112 (cit. a p. 25).
[16] Ulrich Maas e Stephen B Pope. «Simplifying chemical kinetics: intrinsic low dimensional manifolds in composition space». In: Combustion and flame 88.3
(1992), pp. 239–264 (cit. a p. 21).
[17] Amin Maghbouli et al.
Parametric Comparison of Well-Mixed and Flamelet n-dodecane Spray Combustion with Engine Experiments at Well Controlled Boundary Conditions
. Rapp. tecn. SAE Technical Paper, 2016 (cit. a p. 62).
[18] E Mastorakos, TA Baritaud e TJ Poinsot. «Numerical simulations of autoi gnition in turbulent mixing flows». In: Combustion and Flame 109.1 (1997),
[19] Mark PB Musculus, Paul C Miles e Lyle M Pickett. «Conceptual models for partially premixed low-temperature diesel combustion». In: Progress in Energy and Combustion Science
39.2 (2013), pp. 246–283 (cit. a p. 6).
Bibliografia 91
[20] JA van Oijen e LPH De Goey. «Modelling of premixed laminar flames using flamelet-generated manifolds». In: Combustion Science and Technology 161.1
(2000), pp. 113–137 (cit. alle pp. 21, 27).
[21] Yuanjiang Pei et al. «Modelling n-dodecane spray and combustion with the transported probability density function method». In: Combustion and Flame
162.5 (2015), pp. 2006–2019 (cit. alle pp. 56, 71).
[22] Perrine Pepiot-Desjardins e Heinz Pitsch. «An automatic chemical lumping me thod for the reduction of large chemical kinetic mechanisms». In: Combustion Theory and Modelling
12.6 (2008), pp. 1089–1108 (cit. a p. 47).
[23] Stephen B Pope. «Computationally efficient implementation of combustion
chemistry using in situ adaptive tabulation». In: (1997) (cit. alle pp. 21–23).
[24] Nikolai N Semenov.
Chemical Kinetics and Chain Reactions . Oxford Uniersity
Press, London, UK, 1935 (cit. a p. 12).
[25] Scott A Skeen, Julien Manin e Lyle M Pickett. «Simultaneous formaldehyde PLIF and high-speed schlieren imaging for ignition visualization in high pressure spray flames». In: Proceedings of the Combustion Institute 35.3
(2015), pp. 3167–3174 (cit. alle pp. 77, 80).
[26] Henk Kaarle Versteeg e Weeratunge Malalasekera.
An introduction to compu tational fluid dynamics: the finite volume method . Pearson Education, 2007
[27] J Warnatz, U Maas e RW Dibble.
Combustion: physical and chemical funda mentals, modeling and simulation, experiments, pollutant formation.
Springer,
[28] Jürgen Warnatz. «Rate coefficients in the C/H/O system». In: Combustion chemistry
. Springer, 1984, pp. 197–360 (cit. a p. 11).
[29] Charles K Westbrook e Frederick L Dryer. «Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames». In: Combustion science and technology
27.1-2 (1981), pp. 31–43 (cit. a p. 10).
Ulteriore materiale consultato Pubblicazioni e Manuali
[30] A Babajimopoulos et al. «A fully coupled computational fluid dynamics and multi-zone model with detailed chemical kinetics for the simulation of
92 Bibliografia
premixed charge compression ignition engines». In: International journal of engine research 6.5 (2005), pp. 497–512.
[31] Olivier Colin, António Pires da Cruz e Stéphane Jay. «Detailed chemistry based auto-ignition model including low temperature phenomena applied to 3-D engine calculations». In: Proceedings of the Combustion Institute 30.2
(2005), pp. 2649–2656.
[32] Roberto Finesso e Ezio Spessa. «Ignition delay prediction of multiple injections in diesel engines». In: Fuel 119 (2014), pp. 170–190.
[33] Patrick F Flynn et al. «Diesel combustion: an integrated view combining laser diagnostics, chemical kinetics, and empirical validation». In: (1999).
[34] N Peters. «Laminar flamelet concepts in turbulent combustion». In: Sym posium (International) on Combustion . Vol. 21. 1. Elsevier. 1988, pp. 1231– 1250.
[35] Thierry Poinsot e Denis Veynante.
Theoretical and numerical combustion . RT Edwards, Inc., 2005.
[36] MA Singer e SB Pope. «Exploiting ISAT to solve the reaction-diffusion equation». In: Combustion theory and modelling 8.2 (2004), pp. 361–383.
[37] MA Singer, SB Pope e HN Najm. «Operator-splitting with ISAT to model reacting flow with detailed chemistry». In: Combustion Theory and Modelling 10.2 (2006), pp. 199–217.
[38] Michael A Singer, Stephen B Pope e Habib N Najm. «Modeling unsteady reacting flow with operator splitting and ISAT». In: Combustion and Flame 147.1 (2006), pp. 150–162.
[39] Gunnar Stiesch.
Modeling Engine Spray and Combustion Processes . Springer Science & Business Media, 2003.
[40] Mingfa Yao, Zhaolei Zheng e Haifeng Liu. «Progress and recent trends in homogeneous charge compression ignition (HCCI) engines». In: Progress in Energy and Combustion Science 35.5 (2009), pp. 398–437.
Materiale Online
[41] Engine Combustion Network.
url : http://www.sandia.gov/ecn/ .