3

1 - Introduzione

L'elaborato richiede la simulazione del funzionamento di un combustore premiscelato

allo scopo di valutarne la massima potenza termica. Nel nostro caso-studio

l'approccio usato quello LES (Large Eddy Simulation) che consiste nel mediare le

equazioni nello spazio.

I file utilizzati sono stati attinti dal tutorial PitzDaily presente in OpenFoam, in

particolare prelevando le cartelle 0, Costant e System utili a definire le caratteristiche

di partenza. In seguito, in figura 1, riportata la struttura del combustore.

Figura 1: Struttura del combustore con mesh predefinito.

Il primo passaggio effettuato stato quello relativo alla modifica di alcuni parametri

come richiesto dalla traccia.

In particolare, nella cartella 0 sono stati modificati i seguenti parametri presenti:

file "p" : la pressione di esercizio stata portata da 1 bar (1e05 Pa) a 10 bar

(10e5 Pa);

file "k" : l'energia cinetica turbolenta stata portata da 2e-05 a 4e-05 [J/kg];

file "U" : mantenendo costante il vettore velocit (13.3 0 0), sono state

modificate le fluttuazioni in ingresso portandole da (0.04 0.02 0.02) a

(0.04 0.04 0.02).

Nel file "p" la modifica stata apportata uniformemente sia nel campo InternalField

che in Outlet. Questo accorgimento fornisce alla miscela il tempo necessario a far

avvenire la reazione di combustione.

4

Le modifiche apportate ai files "k" e "U" riguardano strettamente il problema della

turbolenza. Difatti, la turbolenza caratterizzata da due contributi principali :

un contributo diretto ed uno indiretto.

Per quanto riguarda il contributo diretto, esso viene espresso tramite le fluttuazioni

della velocit in ingresso. Il contributo indiretto, invece, fa riferimento alla

dissipazione dell'energia cinetica in turbolenza. Detto ci, lecito assumere l'energia

cinetica turbolenta proprio pari all'intensit della turbolenza non risolta.

L'ultima parte dell'elaborato richiede di valutare la stabilit ragionando sul numero

CFL (Courant-Friedrichs-Lewy) definito come:

=

dove:

U la velocit in ingresso della miscela di combustibile;

t il time-step;

x il passo di discretizzazione del mesh.

Al fine di valutare la stabilit del metodo numerico utilizzato, al variare del CFL,

sono stati modificati alcuni parametri presenti nel file controlDict, in modo

particolare il deltaT e il maxCo.

La condizione di stabilit necessaria (ma non sufficiente poich inerisce in essa

anche il concetto di consistenza) per garantire la condizione di convergenza della

soluzione. Si vedr come essa, in condizioni di divergenza, mostrer valori che

oscillano con ampiezza sempre maggiore, ossia il tipico andamento che si manifesta

in caso di instabilit numerica.

L'analisi della stabilit stata svolta con due diverse metodologie. La prima consiste

nel fissare un maxCo ed impostare il comando adjustTimeStep sulla condizione "yes".

Questo equivale a dire che il calcolatore modeller il time-step al fine di garantire che

il CFL sia sempre minore del maxCo predefinito. Se tale parametro viene impostato a

valori "troppo" alti, la soluzione, essendo instabile, non andr a convergenza dato che

la macchina non riuscir pi ad espletare i calcoli. La seconda metodologia consiste

nell'impostare il comando adjustTimeStep sulla condizione "no", in questo modo tutti

i calcoli verranno eseguiti lasciando invariato il time-step predefinito. In questo caso,

come si vedr, non sar possibile "giocare" eccessivamente sul time-step che dovr

essere mantenuto sufficientemente limitato rispetto al valore utilizzato nella

metodologia sopra citata, a causa di condizioni di divergenza troppo "rapide".

5

Si riportano e si discutono, di seguito, i risultati relativi alle simulazioni pi notevoli

e per ognuna di esse verr valutata la massima potenza termica sviluppata dal

combustore.

2 - Analisi delle simulazioni

2.1 - Simulazione con maxCo pari a 0.5, deltaT pari a 5e-0.6 ed

adjustTimeStep impostato su yes

La prima simulazione stata lanciata impostando il numero di Courant massimo pari

a 0.5 (valore predefinito nel file controlDict), il deltaT pari a 5e-06 e adjustTimeStep

su yes, salvo conservando tutte le modifiche gi citate.

Si noti che con questi parametri la macchina ha eseguito i calcoli senza manifestare

errori dato che il maxCo stato mantenuto al di sotto del valore di soglia prestabilito,

inoltre, il tempo di calcolo non stato eccessivamente elevato e si mantenuto

nell'ordine dei 30 minuti.

Al termine della simulazione stato utilizzato il software di post-processing

(visualizzazione grafica dei risultati) ParaView per analizzare il comportamento delle

principali variabili di esercizio del combustore (U, T, k e b).

La velocit, nei primi istanti, data la struttura del combustore e le fluttuazioni pre-

impostate, mostra un andamento vorticoso soprattutto nella zona di discontinuit

all'ingresso (Fig. 2).

Figura 2 : Andamento della velocit all'istante di tempo 0.007.

6

Si nota, invece, che l'andamento della velocit rimane pressoch invariato negli

istanti di tempo successivi, come si pu notare in Fig. 3:

Figura 3 : Andamento della velocit all'istante di tempo 0.122.

Dalla figura, inoltre, si pu notare grazie alla scala graduata che la velocit aumenta

in prossimit dell'uscita a causa del restringimento della sezione di passaggio.

Lo stesso profilo turbolento si pu notare dall'andamento della temperatura nei primi

istanti di tempo. Si riporta, a titolo di esempio, la Fig.4:

Figura 4 : Profilo di temperatura all'istante di tempo 0.007.

7

Dato che lo scopo dell'elaborato quello di analizzare il comportamento della

potenza termica del combustore, si riportano in seguito i profili di temperatura relativi

ad alcuni istanti di tempo caratteristici ( Fig. da 5 a 8), che sono: 0.05, 0.1, 0.15 e 0.2.

Figura 5 : Profilo di temperatura all'istante 0.05.

Figura 6 : Profilo di temperatura all'istante 0.1.

Figura 7 : Profilo di temperatura all'istante 0.15.

8

Figura 8 : Profilo di temperatura all'istante 0.2.

In piena attinenza allo scopo del problema, ossia il calcolo della potenza termica, si

riporta in Tab. 1 il profilo di temperatura in forma tabellare:

Tabella 1 : Temperature agli istanti caratteristici a maxCo=0.5.

T(t50) [K] T(t100) [K] T(t150) [K] T(t200) [K]

1602.55 1576.54 1592.92 1618.09

1646.72 1639.37 1628.25 1640.15

1648.32 1646.81 1635.75 1646.44

1650.63 1648.71 1644.61 1648.08

1645.84 1648.1 1645.05 1289.16

1377.61 1200.31 1292.03 680.78

315.94 418.43 308.84 294.51

293.98 293.25 293.05 293.13

293.93 293.86 293.05 293.09

293.87 293.88 293.07 293.08

Bisogna notare che in paraView queste temperature sono state ottenute impostando il

cursore di sezione su quella di uscita e con una risoluzione dei punti pari a 10.

Si procede successivamente al calcolo della potenza termica, per il quale stata

utilizzata la seguente formula:

=

in cui:

Q la potenza termica;

la portata massica della miscela in ingresso al combustore;

9

cp il calore specifico a pressione costante;

T la variazione di temperatura tra ingresso ed uscita.

Al fine di procedere con il calcolo, sono state assunte delle ipotesi di partenza tra le

quali: il propano (gas combustibile utilizzato) ha un comportamento da "gas perfetto",

il cp calcolato in funzione della variazione di temperatura tramite la formula

= + + 2 + 31, la Tin pari a 293 K e la U (velocit d'ingresso del gas)

impostata a 13.3 m/s.

Per quanto riguarda il calcolo della portata massica, quest'ultimo stato ricavato a

partire da:

=

nella quale:

la densit del propano pari a 493 kg/m3;

U la gi citata velocit d'ingresso del gas;

= (rispettivamente le dimensioni lungo y e lungo z del combustore)

la sezione d'ingresso del combustore pari a 9.7e-04 m2, il cui valore stato

possibile ottenerlo grazie alla conoscenza del file blockMeshDict.

Il valore della restituito da tale procedura di 6.36 kg/s. Grazie alla conoscenza di

tutte queste variabili, si pu procedere, infine, alla valutazione della potenza termica

riportata in Tab. 2:

Tabella 2 : Potenza termica ai diversi istanti di tempo a maxCo=0.5.

P[W] a t50 P[W] a t100 P[W] a t150 P[W] a t200

195461.91 186502.72 192111.62 200951.50

211337.68 208637.30 204596.59 208922.75

211928.66 211370.90 207315.99 211234.38

212783.86 212072.88 210560.05 211839.94

211013.13 211847.33 210722.05 105121.97

126921.83 85853.90 105787.29 16578.95

389.92 2940.20 262.51 23.74

15.38 3.91 0.78 2.03

14.59 13.49 0.78 1.41

13.65 13.80 1.09 1.25

1 Per il propano (C3H8) i valori utili al calcolo del cp(T) sono: a=9.580e-02, b=6.946e-03, c=3.598e-06,

d=7.290e-10.

10

Per dare significato a tali valori, si evince che la potenza massima, erogata dal

combustore, si aggira intorno ai 210 kW.

2.2 - Simulazione con maxCo pari a 1.5, deltaT pari a 5e-0.6 ed

adjustTimeStep impostato su yes

Al fine di velocizzare i calcoli, stato aumentato il massimo numero di Courant

(CFL) lasciando attiva l'opzione adjustTimeStep. Lo scopo stato quello di ricercare

il massimo CFL che rendesse possibile la convergenza della soluzione. Il primo

tentativo stato fatto per un maxCo pari a 1.5. Ci ha velocizzato i calcoli senza

indurre divergenza. Dato che i profili di temperatura sono risultati essere molto

simili, si riportano solo le tabelle relative alla temperatura e alla potenza termica

(Tab. 3 e 4). Tabella 3 : Temperature agli istanti caratteristici a maxCo=1.5.

T(t50) [K] T(t100) [K] T(t150) [K] T(t200) [K]

1551.51 1473.89 1458.95 1526.91

1694.91 1683.64 1666.69 1677.91

1688.17 1679.63 1664.60 1671.50

1681.31 1669.94 1658.27 1665.08

1673.20 1659.98 1654.64 1659.67

1435.36 1473.27 1624.60 1049.14

835.90 460.15 423.72 611.22

406.63 351.32 305.78 306.98

300.15 300.68 295.27 295.46

296.38 295.48 294.44 294.68

295.82 294.73 294.15 294.41

Tabella 4 : Potenza termica ai diversi istanti di tempo a maxCo=1.5.

P[W] a t50 P[W] a t100 P[W] a t150 P[W] a t200

178147.43 153846.53 149440.46 170185.53

229635.16 225262.92 218794.61 223061.90

227013.45 223721.04 218005.93 220617.11

224366.13 220024.90 215629.05 218186.89

221263.70 216269.39 214274.00 216153.21

142656.49 153661.97 203281.96 58601.32

30530.61 4385.17 3109.72 11843.93

2575.93 1117.14 209.47 230.20

114.79 123.49 35.76 38.85

53.47 39.17 22.56 26.46

11

In questo caso la potenza termica massima valutata di circa 225 kW. Questo

aumento riconducibile ad una perdita di accuratezza dovuta ad utilizzo di un

maggiore passo di tempo.

2.3 - Simulazione con maxCo pari a 1.91, deltaT pari a 5e-0.6 ed

adjustTimeStep impostato su yes

L'ultimo numero di Courant, ottenuto grazie ad una iterazione di calcolo, che

garantisce la convergenza pari a 1.91. Con valori maggiori la simulazione si

bloccata ad istanti di tempo minori in relazione a CFL maggiori.

L'errore riportato nella figura seguente:

Figura 9 : Errore riscontrato.

Esso afferisce all'impossibilit del calcolatore di effettuare il calcolo in virgola

mobile in condizioni di divergenza.

Si riportano i valori di temperatura e potenza termica nelle tabelle 5 e 6:

Tabella 5 : Temperature agli istanti caratteristici a maxCo=1.91.

T(t50) [K] T(t100) [K] T(t150) [K] T(t200) [K]

1355.05 1416.50 1303.51 1316.14

1696.52 1734.95 1669.92 1722.63

1688.57 1715.84 1663.80 1735.70

1677.12 1712.89 1656.20 1714.28

1640.84 1466.64 1652.63 1687.39

1595.13 1224.77 1653.25 1377.23

1075.99 890.02 1649.56 703.20

318.98 587.03 435.71 463.63

296.52 410.36 302.40 342.90

295.78 324.13 295.85 301.95

295.50 302.75 295.19 298.17

12

Tabella 6 : Potenza termica ai diversi istanti di tempo a maxCo=1.91.

P[W] a t50 P[W] a t100 P[W] a t150 P[W] a t200

121102.85 137383.04 108476.21 111485.94

230264.46 245638.00 220017.31 240635.35

227168.47 237908.25 217704.55 245944.82

222759.59 236730.05 214855.62 237284.70

209175.49 151697.67 213526.19 226711.38

192877.02 90906.19 213756.66 126822.32

62963.55 36612.01 212387.43 18296.30

446.34 10398.94 3508.66 4517.01

55.76 2688.96 152.18 929.91

43.92 544.49 45.04 144.66

Anche in questa simulazione, in virt del sopra citato ragionamento, vi un aumento

della potenza termica massima che si aggira intorno ai 240 kW.

2.4 - Simulazione con adjustTimeStep impostato su no e deltaT pari a

5e-0.7

In questa simulazione il numero CFL non viene mantenuto al di sotto della soglia

dato che il calcolatore non adatta ad ogni istante di tempo la relativa discretizzazione

spaziale. Ci imputabile all'impostazione della funzione adjustTimeStep su "no".

Prima di giungere alla scelta del deltaT pari a 5e-07, sono stati provati diversi time-

step, maggiori e uguali di quello di default (5e-06), per, in tutti questi casi il

calcolatore ha presentato l'errore gi citato (cfr. paragrafo 2.3). Data l'impossibilit di

procedere in ogni altro modo, l'unica scelta possibile stata quella di diminuire

ulteriormente il deltaT. In questo caso il tempo impiegato dal calcolatore per la

simulazione aumentato drasticamente (10 ore).

In questo caso, data la maggiore accuratezza della simulazione, si riportano i profili

di temperatura agli stessi istanti di tempo, vedendo come si passa da un profilo

lievemente increspato ad uno marcatamente corrugato.

Seguono le immagini (Fig. da 10 a 13).

13

Figura 10 : Profilo di temperatura all'istante 0.05 a deltaT=5e-07.

Figura 11 : Profilo di temperatura all'istante 0.1 a deltaT=5e-07.

Figura 12 : Profilo di temperatura all'istante 0.15 a deltaT=5e-07.

14

Figura 13 : Profilo di temperatura all'istante 0.2 a deltaT=5e-07.

Sulla base di quanto gi effettuato (cfr. paragrafo 2.1), si inseriscono le tabelle

relative alle temperature e alle potenze termiche (Tab. 7 e 8).

Tabella 7 : Temperature agli istanti caratteristici a deltaT=5e-07.

T(t50) [K] T(t100) [K] T(t150) [K] T(t200) [K]

1413.96 1299.24 1470.34 1418.26

1623.74 1594.54 1643.97 1557.46

1090.29 1614.79 1643.92 1604.41

845.83 1633.98 1642.93 1522.57

829.29 1634.48 1644.10 1581.31

847.92 1539.09 1618.15 1537.46

1012.93 1048.26 1056.22 1421.98

987.90 373.43 321.71 423.12

538.69 294.12 293.63 296.71

303.34 293.24 293.30 293.96

15

Tabella 8 : Potenza termica ai diversi istanti di tempo a deltaT=5e-07.

P[W] a t50 P[W] a t100 P[W] a t150 P[W] a t200

136682.92 107470.89 152791.80 137869.57

202973.05 192672.48 210324.57 180110.19

65368.41 199777.12 210306.18 196113.54

31596.38 206672.01 209942.29 168806.24

29832.60 206853.79 210372.39 188124.48

31823.59 174096.97 200972.89 173570.04

53027.72 58461.69 59732.38 138901.67

49377.15 1651.97 498.04 3090.35

7806.29 17.58 9.87 58.81

167.98 3.68 4.62 15.06

Da queste tabelle si evince come si riscontrano delle maggiori fluttuazioni dei valori

di temperatura e potenza termica, con quest'ultima che nel suo valore massimo si

aggira a ca. 209 kW.

3 - Considerazioni e conclusioni

Da tutte le simulazioni effettuate si evince che un aumento del numero CFL comporta

una diminuzione del tempo di calcolo a scapito di una perdita di accuratezza. Questo

ragionamento non pu essere spinto eccessivamente verso CFL troppo grandi,

poich, in tali circostanze, incombono condizioni di divergenza.

Per meglio comprendere la relazione che sussiste tra il caso base e quello a CFL

massimo, si procede ad una comparazione grafica riportata in Fig. 14:

Figura 14 : Grafico comparativo all'istante t100.

0,00

50000,00

100000,00

150000,00

200000,00

250000,00

300000,00

1 2 3 4 5 6 7 8 9 10

0.5

1.91

16

In questo grafico si vede come il caso a CFL massimo sia caratterizzato da potenze

termiche sensibilmente differenti rispetto a quello base: si nota come il combustore

viene sovradimensionato di ca. 40 kW rispetto al caso base. Ci riconducibile ad

una perdita di accuratezza poich la macchina, nel calcolo del CFL massimo, ha

utilizzato dei time-step maggiori che, come ben noto dalle propriet del metodo di

Eulero esplicito, comportano un aumento dell'errore globale.

Per meglio chiarificare il vincolo che lega il CFL all'accuratezza della soluzione,

stata lanciata la simulazione a deltaT pari a 5e-07. A tal fine riportato in Fig. 15 un

grafico significativo sul confronto delle potenze.

Figura 15 : 2 grafico comparativo a t100.

Si evince come le due curve possano quasi essere sovrapposte presentando livelli

assimilabili di accuratezza e, quindi, di potenze termiche in uscita: in questo caso si

percepisce un lieve sovradimensionamento (ca. 5 kW) nel caso base che risulta essere

il meno accurato tra i due.

0,00

50000,00

100000,00

150000,00

200000,00

250000,00

1 2 3 4 5 6 7 8 9 10

0.5

e-07