UN METODO MARKOVIANO DI GENERAZIONE DI SEQUENZE DIINSOLAZIONE PER IL DIMENSIONAMENTO E LA SIMULAZIONE

DI IMPIANTI DI CAPTAZIONE

Giuseppe Oliveti - Natale Arcuri - Salvatore Ruffolo

Dipartimento di Ingegneria Meccanica- Università della Calabria 87030 Rende (Cosenza) - e-mail:giuseppe.oliveti@unical.it

SOMMARIO

Viene presentato un metodo di generazione di sequenze di insolazione che utilizza una libreria di matrici diMarkov e i parametri delle distribuzioni normali degli indici di serenità mensili, idoneo a descrivere statisticamentela variabilità mensile e annuale dell’energia solare incidente. Il metodo di generazione è stato utilizzato per valutaregli effetti della variabilità climatica sugli impianti solari con accumulo interstagionale.

ABSTRACT

An insolation sequence generation method that uses a Markov’s library of transition matrices are presentedtogether with normal distribution parameters of monthly clearness indices capable of statistically describing monthlyand yearly solar incident energy variability. The generation method has been used to evaluate the effects of climaticvariability on solar plants with interseasonal storage.

1. INTRODUZIONE

Il dimensionamento e le prestazioni degli impianti che utilizzano la radiazione solare sono di solito ottenuteimpiegando i dati climatici relativi all’anno storico, e pertanto nelle valutazioni termiche ed economiche, la frazionedel fabbisogno fornito è assunta costante nel corso della vita dell’impianto. Nel funzionamento effettivo la variabilitàclimatica rende mutevole la frazione solare fornita. Più realistiche sono quindi le valutazioni in condizioni climatichevariabili attraverso la determinazione della curva di probabilità cumulata, la quale stima i livelli di probabilitàassociati all’energia fornita dall’impianto. Queste valutazioni richiedono la disponibilità di sequenze annuali dellaradiazione, rappresentative della variabilità climatica della località considerata.

Le matrici di Markov a un passo esprimono il più semplice modello di previsione stocastica delle sequenze diinsolazione. Se si dispone di dati sperimentali di insolazione, ad esempio una successione annuale dell’indice diserenità giornaliero k, rapporto tra l’energia incidente sul piano orizzontale e la corrispondente energia extratmosferica,la costruzione della matrice di Markov è facile, e richiede i seguenti passi:a) discretizzazione dei valori di k in n classi;b) costruzione di una matrice quadrata di dimensioni nxn il cui elemento generico (i,j) contiene il numero di

transizioni dalla classe i alla classe j. Ogni elemento della matrice è in sostanza un contatore che registra ilnumero di volte in cui le due classi si sono succedute;

c) normalizzazione della matrice ottenuta dividendo ogni elemento di una riga per la somma degli elementi presentinella riga considerata.

Si è ottenuta in tal modo una “Matrice di Transizione Autoregressiva” ad un passo, i cui elementi p(i,j) esprimonola probabilità di transizione dell’indice di serenità dalla classe i alla classe j.

La generazione di una sequenza annuale si ottiene operando sulla matrice nel modo seguente. Inizialmente sisceglie casualmente una delle righe e si genera un numero casuale R, utilizzando una distribuzione uniformenell’intervallo tra zero e uno. Si sommano sulla riga individuata gli elementi a partire dal primo, fino a quando laloro somma risulta non minore di R, si individua in tal modo una classe o indice di colonna. All’interno di taleclasse si estrae in maniera casuale un valore dell’indice di serenità che costituisce il primo elemento della sequenzasimulata. La classe prima individuata costituisce il nuovo indice di riga da cui ripartire per generare un nuovo valoredell’indice di serenità.

In questo lavoro sono state utilizzate matrici di transizione di Markov (MTM) messe a punto da Aguiar et al.,[1], attraverso dati sperimentali relativi a diverse località. La libreria è formata da dieci matrici, tutte di dimensioni10x10, ogni matrice è relativa ad un definito intervallo dell’indice di serenità mensile sul piano orizzontale. Per unadeterminata località in genere sono noti i valori degli indici di serenità mensile relativi all’anno storico. Operandosulla matrice corrispondente è possibile generare sequenze mensili dell’indice di serenità giornaliero, statisticamenteindistinguibile dalla sequenza relativa all’anno storico.

2. IL METODO DI GENERAZIONE DELLE SEQUENZE DI INSOLAZIONE

A causa della variabilità climatica l’indice di serenità mensile relativo ad un generico mese, differisce dal valoremedio fornito dall’anno storico. La variabilità è stata simulata facendo l’ipotesi che la distribuzione dell’indice diserenità mensile sia di tipo normale. Se quindi per i singoli mesi dell’anno storico si conosce oltre al valore medio

k anche la deviazione standard σ k , è possibile utilizzare tali distribuzioni per generare casualmente un valore

dell’indice di serenità mensile, e successivamente, attraverso l’impiego delle MTM, disporre di una sequenza mensiledell’indice di serenità giornaliero. Generando casualmente dodici valori di k, uno per ogni mese secondo lacorrispondente distribuzione normale è possibile ottenere una sequenza annuale della radiazione. La proceduraconsente di generare un numero qualsivoglia di sequenze annuali, che si caratterizzano per avere ciascuna valorimensili e annuali dell’energia diversi da quelli relativi all’anno storico. Nella generazione mensile si è posto ilvincolo che la sequenza generata restituisca un valore di energia mensile corrispondente al valore di k generato.Tale procedura di generazione simula l’effettiva variabilità climatica. Come caso particolare è possibile generaresequenze con il vincolo che per ogni mese l’energia sia quella corrispondente all’anno storico.

Per l’utilizzo delle sequenze generate occorre riportare i valori dell’energia giornaliera dal piano orizzontale alpiano di captazione. Tra le varie correlazioni disponibili in letteratura è stata scelta quella di Liu e Jordan [2],perché ha fornito una più accurata verifica tra i valori di radiazione misurati sul piano orizzontale e su un pianoinclinato di 39,3° (latitudine della stazione attinometrica dell’Università della Calabria).

Nella fig.1 è riportato un confronto tra l’energia mensilesul piano di captazione, determinata sperimentalmente, equella ottenuta con le sequenze generate sul pianoorizzontale e la correlazione di Liu e Jordan per il riportosulla superficie inclinata. Per la generazione sul pianoorizzontale sono stati utilizzati i valori dell’indice diserenità mensile ricavati sperimentalmente, con il vincoloche la distribuzione restituisca lo stesso valore mensiledell’energia misurata. I dati di radiazione sul pianoorizzontale e sul piano di captazione sono stati misuratinel periodo gennaio 1993 - dicembre 1996. Con i valorimensili dell’indice di serenità relativi all’anno storico, sonostati generati con le MTM dieci anni successivi, in mododa ottenere in uno stesso mese diverse distribuzionidell’indice di serenità giornaliero. Il confronto mostral’ottimo accordo tra i valori mensili dell’energia misuratae calcolata, ed evidenzia che l’algoritmo di calcolo conducead una sottostima dell’energia sul piano di captazione,imputabile in gran parte alla correlazione, con unadistribuzione dello scarto relativo percentuale, fig.2, cheassegna l’82,3% dei valori alla banda –6% +2%; il 91,7%dei valori tra –10% +2%; e il 94,2% tra –10% +10%. Notal’energia incidente sul piano di captazione, l’energia

200

400

600

800

200 400 600 800

Fig. 1 - Energie mensili misurate e calcolatesul piano di captazione.

Energia misurata(MJ/m2)

Ene

rgia

cal

cola

ta (

MJ/

m2 )

raccolta può essere determinata stimando l’efficienza media giornaliera del collettore. Ciò richiede la conoscenzadelle caratteristiche del collettore installato, le condizioni di alimentazione, e il valore dell’irraggiamento medio,calcolabile una volta nota la durata di insolazione. I miglioririsultati sono stati ottenuti stimando la durata di insolazionesul piano orizzontale con una procedura che operacasualmente tra un valore minimo e il valore massimoextratmosferico. Il valore minimo è stato determinato comerapporto tra l’energia giornaliera generata e l’irraggiamentomedio, relativo a giornata serena, valutato con la proceduraASHRAE [3]. I valori estremi della durata di insolazionecalcolati in tal modo sono quindi relativi a giornate o frazionidi giornate con cielo sereno. Le durate di insolazione sulpiano orizzontale generate casualmente tra tali valori estremisono state poi riportate sul piano di captazione con i consuetimetodi.

Disponendo di dati sperimentali relativi alfunzionamento di un sistema di accumulo interstagionaledell’energia solare, costituito da un campo di captazione di91,2 m2, di un sistema di accumulo di 500 m3, e di unedificio di 1750 m3 [4], la procedura di generazione messaa punto è stata sottoposta a verifica sperimentale. Nella fig.3sono riportati gli andamenti sperimentali della temperaturadi alimentazione del campo di captazione, dell’energiaincidente, dell’energia captata registrati nel periodo giugno 1996, maggio 1997.

L’energia giornaliera captata con le sequenze generate è stata calcolata considerando la caratteristica del campodi captazione a tubi di calore dell’impianto sperimentale [5], identiche condizioni di alimentazione e di temperaturadell’aria esterna, e gli stessi vincoli impiegati nell’impianto sperimentale per il funzionamento automatico. Inparticolare è stata valutata la possibilità di raccolta giornaliera dell’energia mediante un confronto tra l’irraggiamentodi soglia richiesto per ottenere un incremento di temperatura della portata di 1°C, tra l’ingresso e l’uscita dal campodi captazione, e l’irraggiamento medio ottenuto attraverso la procedura di generazione. I risultati del confronto trai valori mensili dell’energia captata nell’impianto sperimentale e quella calcolata con le sequenze generate, sonoriportati nella fig.4. Anche in questo confronto sono stati considerati dieci anni di sequenze generate, con valorimensili e annuali dell’energia sul piano orizzontale identici a quelli sperimentali. Il massimo scarto mensile si èverificato nel mese di dicembre con una sottostima dell’energia sperimentale del 20%. Se si considera l’energiaannuale captata, l’impiego delle sequenze generate conduce ad una sottostima dell’energia che non supera il 2,5%.

0

10

20

30

40

50

60%

-inf

÷-14

-14÷

-10

-10÷

-6

-6÷-

2

-2÷2 2÷6

6÷10

10÷1

4

14÷i

nf

Fig. 2 - Distribuzione dell'errore relativo percentuale tra l'energia misurata e quella calcolata sul piano di captazione.

0

20

40

60

80

°C

0

20

40

60

80

GJ

giug

no

lugl

ioag

osto

sett

embr

eot

tobr

e

nove

mbr

e

dice

mbr

e

genn

aio

febb

raio

mar

zoap

rile

mag

gio

Tf

Ec

Ei

Fig. 3 - Temperatura di alimentazione del campo di captazione, energia incidente ed energia captata nell'impianto sperimentale (giugno 1996 - maggio 1997).

0

10

20

30

40

GJ

Gen

naio

Febb

raio

Mar

zo

Apr

ile

Mag

gio

Giu

gno

Lug

lio

Ago

sto

Sette

mbr

e

Otto

bre

Nov

embr

eD

icem

bre

Energie captateEnergie calcolate

Fig. 4 - Confronto tra le energie mensili captate e quelle calcolate con le sequenze generate.

3.VARIABILITÀ DELL’ENERGIA UTILE NEGLI IMPIANTI CON ACCUMULO INTERSTAGIONALE

In questi ultimi anni sono stati costruiti impianti solari con accumulo interstagionale dell’energia, destinati alriscaldamento invernale e alla produzione di acqua calda sanitaria. Questi impianti sono noti come CSHPSS (CentralSolar Heating Plants with Seasonal Storage), e si prestano alla fornitura di energia a complessi abitativi[6],[7],[8],[9],[10].Nella nostra indagine, il dimensionamento termico è stato ottenuto con una procedura automatica iterativa, cheimpiega tre codici di calcolo disposti in cascata [11], [12] e per la simulazione un codice di calcolo convalidatomediante confronti con dati sperimentali [4]. I dati climatici utilizzati sono quelli relativi all’anno storico.

Queste procedure automatiche di calcolo sono state impiegate per una estesa analisi termica considerando im-pianti di taglia diversa ubicati in tre località climaticamente differenti del territorio italiano (Milano, Pisa e Cosenza).All’analisi termica è stata poi affiancata una analisi di tipo economico con la formulazione di una funzione costoglobale, il cui minimo rispetto all’area dei pannelli e al volume di accumulo, ottenuto tenendo conto dei vincolitermici tra queste due grandezze, consente di determinare l’energia che deve essere erogata in un anno dall’impian-to. Il resto del fabbisogno è fornito dall’impianto di integrazione [13].

In questo articolo per brevità è stata riportata l’analisi della variabilità climatica per gli impianti ubicati a Pisa.Nella Tab. I sono riportati i valori ottenuti con l’analisi economica, dell’area di captazione, del volume di accumuloe dell’energia che l’accumulatore eroga per quattro impianti di taglia diversa. Nella fig.5 è riportato un esempio diandamento annuale della temperatura media dell’accumulatore, dell’energia mensile incidente sui collettori nell’annostorico, dell’energia estratta dall’accumulatore. Nella Tab. II sono indicati i parametri delle distribuzioni mensilidell’indice di serenità [14].

Con la procedura automatica messa a punto, sono stategenerate sequenze annuali di radiazione, e con il codice disimulazione dell’impianto, calcolati i valori dell’energia chel’accumulatore è in grado di erogare. In condizioni divariabilità climatica, per determinare il prelievo di energiadall’accumulatore si è utilizzato il seguente criterio: dovendol’accumulatore funzionare con continuità, l’energia estraibileè stata determinata in modo da ottenere un profilo ditemperatura dell’accumulatore prossimo a quello dell’annostorico. In tal modo le perdite termiche dell’accumulatorediventano quasi indipendenti dalle variazioni dell’energiaintrodotta. Ciò è stato ottenuto confrontando il valore della temperatura media dell’accumulatore con il corrispondentevalore relativo all’anno storico, e distribuendo la maggiore o minore energia disponibile in modo proporzionale aiprelievi corrispondenti all’anno storico.

Per ridurre i tempi di elaborazione, le prestazioni degli impianti al variare delle sequenze annuali di insolazione,sono state valutate con una procedurache sfrutta le risposte del sistema campodi captazione-accumulatore allevariazioni mensili dell’energia incidente.A tal proposito sono stati calcolati conil codice di simulazione dell’impianto ivalori dell’energia annuale fornitadall’accumulatore considerandoincidente sul piano di captazione lasequenza dell’anno storico modificataogni volta in un solo mese. La modificaè stata realizzata generando con le MTMuna prima sequenza mensile con un

valore di k kk

= +3σ , e una seconda

con k kk

= −3σ . I parametri della

distribuzione sono quelli relativi al meseconsiderato.

Ac V Qu,0

m2 m3 GJ

2479 13101 54571252 6484 2729640 3298 1364338 1694 677

Tab. I - Area di captazione, volume di accumulo ed energia utile erogata nelle condizioni di ottimo economico per gli impianti ubicati a Pisa.

k σ kGennaio 0.322 0.060Febbraio 0.356 0.066Marzo 0.375 0.062Aprile 0.420 0.057Maggio 0.452 0.047Giugno 0.466 0.043Luglio 0.513 0.051Agosto 0.490 0.043Settembre 0.459 0.047Ottobre 0.424 0.052Novembre 0.334 0.060Dicembre 0.325 0.053

Tab. II - Parametri delle distribuzioni normali dell'indice di serenità mensile per Pisa.

Se si ipotizza che, relativamente all’energia estratta, il comportamento del sistema campo di captazione-accumulatore nell’intorno del punto di funzionamento corrispondente all’anno storico, sia linearizzabile, allora èpossibile porre l’energia annuale che l’accumulatore può erogare sotto la forma:

Q Q S a E Eu u i i ii

= + ⋅ −=∑, ,( )0 0

1

12

(1)

La relazione precedente consente il calcolo dell’energia utile Qu a partire dell’energia utile che l’impianto eroga

nell’anno storico Qu,o

attraverso la somma delle risposte del sistema alle singole variazioni mensili dell’energiaincidente(Ei-Ei,o), calcolate con riferimento all’anno storico. S è l’area di captazione, a

i sono dei coefficienti che

hanno il significato di fattori di utilizzo delle variazioni mensili dell’energia incidente, e sono calcolabili, nellaformulazione adottata, per ogni mese come rapporto tra la variazione di energia utile annuale conseguente allavariazione mensile dell’energia incidente, e quest’ultima. La (1) consente quindi il calcolo dell’energia chel’accumulatore eroga qualunque sia la sequenza annuale incidente.

Nella fig.6 è riportato per Pisa l’andamento dei valori mensili dei coefficienti di utilizzazione a partire da aprileche è il mese di fine riscaldamento invernale e di inizio della carica dell’accumulatore. La figura mostra chemediamente la quota di energia utile annuale conseguente alle variazioni mensili dell’energia incidente è variabileda 0,31 (aprile) a 0,94 (marzo). Ciò evidenzia che il sistema è particolarmente sensibile alle variazioni di energiaincidente nel periodo di riscaldamento invernale (novembre- marzo), per il fatto che tali variazioni in larga misurarappresentano variazioni di energia utile erogabile dall’accumulatore. Inoltre tali coefficienti presentano unadipendenza non marcata dal volume di accumulo.

Un confronto sistematico tra l’energia utile calcolata con il codice di calcolo e quella valutata con la relazioneprecedente, impiegando sequenze di radiazione generate per un numero di anni pari a 100, ha evidenziato chel’errore relativo percentuale non supera l’1,8% ( vedi fig.7).

La procedura di calcolo messa a punto è stata utilizzata per una estesa campionatura dell’energia utile erogabiledagli impianti considerati. E’ stato considerato un numero elevato di sequenze annuali in modo da ottenere per tuttii mesi dell’anno, e per i diversi anni, una estesa rappresentazione della variabilità climatica. I risultati ottenuti sonostati successivamente sottoposti ad una indagine statistica con la costruzione delle curve di probabilità cumulatadell’energia estratta dall’accumulatore, e ancora meglio, della frazione di energia estratta, calcolata con riferimentoall’energia fornita nell’anno storico. Queste curve forniscono la probabilità che l’impianto eroghi una frazione dienergia maggiore o uguale al valore corrispondente, e sono di grande interesse nella fase di studio e di analisi delleprestazioni, ed essenziali per una corretta progettazione del sistema di produzione dell’energia di integrazione. Inparticolare forniscono la frazione di energia che l’impianto è certamente in grado di erogare ( p = 99%), o di nonerogare (ad esempio p = 10%), e la conferma che ad un valore di energia utile maggiore o uguale a quello relativo

0

0.2

0.4

0.6

0.8

1

Apr

ile

Mag

gio

Giu

gno

Lug

lio

Ago

sto

Set

tem

bre

Ott

obre

Nov

embr

e

Dic

embr

e

Gen

naio

Feb

brai

o

Mar

zo

V = 13004 m3V = 6484 m3

V = 3298 m3

V = 1694 m3

Fig. 6 - Coefficienti di utilizzo delle variazioni mensili dell'energia incidente per gli impianti ubicati a Pisa

1000

2000

3000

4000

5000

6000

Ene

rgia

util

e co

dice

(G

J)

1000 2000 3000 4000 5000 6000

Energia utile (GJ)

V = 13101 m3

V = 6484 m3

V = 3298 m3

V = 1694 m3

Fig. 7 - Confronto tra l'energia utile calcolata con il codice di calcolo e quella valutata con la relazione (1).

all’anno storico è associata una probabilità del 50%.Un esame dei valori della superficie di captazione e del

volume di accumulo richiesti nelle diverse tipologie diimpianti, ha evidenziato che, nelle condizioni di ottimoeconomico considerato, il legame funzionale tra l’area dicaptazione e il volume di accumulo è di tipo lineare eomogeneo, e dello stesso tipo è pure il legame tra l’energiautile fornita e il volume di accumulo. Ciò rende,considerando la relazione (1) e la fig.6, la frazione di energiaQu/Qu,o sostanzialmente indipendente dal volume diaccumulo, e così pure la curve di probabilità cumulata. Nellafig.8 sono state riportate le curve di probabilità cumulataottenute per Pisa. Gli impianti sono in grado di erogare conuna probabilità del 99%, 85% dell’energia erogabilenell’anno storico; con una probabilità del 90%, il 92%; conuna probabilità del 70% il 97%; con una probabilità del50% il 100%; con una probabilità del 30%, il 104%, e infineuna probabilità del 10%, il 109%.

4. CONCLUSIONI

E’ stato messo a punto un metodo di generazione di sequenze annuali della radiazione in grado di descriverestatisticamente la variabilità climatica. La validità della procedura, è stata accertata mediante un confronto tral’energia captata in un impianto sperimentale e quella ottenuta con un modello di un campo di captazione cheutilizza le sequenze generate.

Il metodo di generazione è stato impiegato per studiare la variabilità delle prestazioni di impianti solari conaccumulo interstagionale progettati con un criterio economico, utilizzando le sequenze di radiazione relative all’annostorico. L’analisi ha condotto alla formulazione di un metodo di calcolo dell’energia utile annuale fornitadall’accumulatore in condizioni di variabilità climatica, che impiega l’energia erogata nell’anno storico e i coefficientidi utilizzo delle variazioni mensili dell’energia incidente. Quest’ultimi sono calcolabili attraverso le sequenze gen-erate.

L’indagine, inoltre, ha consentito di ricavare le curve di probabilità cumulata delle frazioni annuali dell’energiautile, calcolate con riferimento all’energia fornita nell’anno storico, e di stimare i corrispondenti livelli di probabilità.In particolare è stata valutata la frazione di energia utile che l’impianto è certamente in grado di erogare (p = 99%),e di non erogare. L’indagine ha evidenziato che le curve di probabilità possono essere ritenute indipendenti dalvolume di accumulo.

Queste analisi completano la conoscenza delle prestazioni degli impianti di captazione consentendo di verificarein termini probabilistici la bontà del dimensionamento e le prestazioni ottenibili, inoltre forniscono prezioseinformazioni per il dimensionamento dei sistemi di integrazione.

BIBLIOGRAFIA

[1] Aguiar R. J., Collares-Pereira M., Conde J.P., Simple Procedure for Generating Sequences of Daily RadiationValues Using a Library of Markov Transition Matrices, Solar Energy, 1988 Vol. 40, pp. 269-279.

[2] Klein S. A., A Design Procedure for Solar Heating System” Solar Energy, 1977, Vol. 19, 325

[3] Ashrae Handbook – Fundamentals – 1989

[4] Oliveti G., Arcuri N., Ruffolo S., First Experimental Results from a Prototype Plant for the InterseasonalStorage of Solar Energy for the Winter Heating of Buildings, Solar Energy, 1998, Vol. 62, No 4, pp. 281-290

[5] Arcuri N., Cucumo M., Oliveti G., Determinazione Sperimentale dell’Efficienza Istantanea di un Campo diCollettori, 1995, CDA Ricerca N. 7, pp. 707-715

0

0.25

0.5

0.75

1

0.8 0.9 1 1.1 1.2Qu/Qu,0

V = 13101 m2V = 6484 m2V = 3298 m2V = 1694 m2

Pro

babi

lità

cum

ulat

a

Fig. 8 - Curve di probabilità cumulata della frazione di energia utile per gli impianti ubicati a Pisa.

[6] Grabrielsson E., Seasonal Storage of Thermal Energy – Swedish Experience,. Trasm. ASME, Journal ofSolar Engineering, 1988, 110, 203.

[7] Bokhoven T. P. and de Geus A.C., Large Solar Heating Syistem with Seasonal Storage for Bylb Drying inLisse the Netherlands,. EuroSun’ (1996) , Freiburg, Germania, Vol. I pp. 77-82.

[8] Schulz M. E. and Fisch M. N., Central Solar Heating Plant with Seasonal Storage in Hamburg-Gramfeld,.EuroSun, Freiburg, Germania, Vol. I pp. 357-362.

[9] Dalenback J. O., Central Solar Heating Plants with Seasonal Storage. Status Report IEA Solar Heating andCooling Programme, Task 7 Technical Report, Swedish Council for Building” Research, Document D14,Stockholm, Sweden, 1990

[10] Fisch M. N. Guigas M. and Dalendak J. O., Large-Scale Solar District Heating – Status and Future inEurope, EuroSun’, Freiburg, Germania, 1966 Vol. I pp. 157-170.

[11] Cucumo M., MarinelliV., Oliveti G. and Sabato A., Long Term Solar Storage througt Underground WaterTanks for the Heating of Housing, Final Report EUR 9774 EN, 1985

[12] Oliveti G. and Arcuri N., Prototype Experimental Plant for the Interseasonal Storage of Solar Energy for theWinter Heating of Buildings: Descirption of Plant and its Functions,. Solar Energ, 1995 Vol. 54 pp. 85-97

[13] Oliveti G., Arcuri N., Ruffolo S., Impianti di Riscaldamento con Accumulo Interstagionale dell’EnergiaSolare. Limiti di Convenienza Economica,. CDA, 1998, No. 5 Maggio, pp. 507-516.

[14] Atlante Europeo della Radiazione Solare (1984) - Commissione delle Comunità Europee – Verlang TuvRheinland.