Filippo Zonta

Una tecnologia che trova un crescente impiego - soprattutto nei Paesi del nord - nel mercato dei sistemi di riscaldamento e condizionamento, è costituita dalle pompe di calore accoppiate a terreno. Si tratta di dispositivi che utilizzano il terreno come mezzo per realizzare grandi accumuli di calore, consentendo lo stoccaggio a basso costo di calore (e/o freddo) per periodi anche di 6 mesi. É così possibile sfruttare al meglio le risorse energetiche disponibili, accumulando calore estivo gratuito (o calore di scarto industriale) per estrarlo poi per il riscaldamento invernale nella stagione successiva. Viceversa, è possibile smaltire il calore nel sottosuolo in un ciclo di raffrescamento estivo senza ricorrere ai costosi condizionatori tradizionali.

Figura 1: Esempio di campo termico generato da un accumulo a terreno (soluzione di Laplace Ñ ² T=0 in regime stazionario).
I cerchi nella zona centrale indicano le sezioni dei tubi verticali.

Per realizzare accumuli a terreno di grandi dimensioni, minimizzando nel contempo l’ingombro richiesto dall’impianto in superficie, in genere si ricorre a particolari scambiatori a tubi verticali installati nel sottosuolo a profondità comprese tra i 20 e i 100 m. In Fig. è riportata una sezione trasversale di un accumulo a tubi verticali.

I metodi attuali di progettazione prevedono l’uso intensivo della simulazione numerica al calcolatore per giudicare l’idoneità funzionale di un progetto e ottenere stime analitiche dei costi operativi di funzionamento degli impianti. I modelli da impiegare sono strettamente legati alla situazione idrogeologica del terreno dell’accumulo.

La presente ricerca è mirata a quantificare il diverso comportamento di uno scambiatore verticale in due situazioni idrogeologiche estreme:

• suolo roccioso impermeabile puramente conduttivo;
• terreno poroso saturo d’acqua con meccanismi di trasmissione del calore influenzati dal moto dell’acqua (convezione termica).

Il modello qui adottato considera un singolo scambiatore a tubi ad U oppure a due tubi concentrici. Lo scambio termico fluido-terreno viene semplificato ricorrendo al concetto di resistenza termica effettiva fluido-terreno (Hellstrom, 1991). L’acquifero, oltre una certa distanza dalla parete dello scambiatore, è perfettamente stagnante e l’unico elemento di perturbazione dello stato di quiete della falda è dovuto alle spinte ascensionali subite dal fluido riscaldato in prossimità dello scambiatore interrato.

Quando un fluido scorre attraverso una matrice solida porosa, per esempio attraverso un terreno sabbioso o ghiaioso, la soluzione delle equazioni del moto diventa intrattabile, a causa della tortuosità degli interstizi tra le particelle solide della matrice.

Equazioni approssimate per descrivere i profili di velocità e temperatura sono ottenibili mediando le grandezze d’interesse in un opportuno volume di controllo, maggiore del volume delle particelle solide e degli interstizi del mezzo poroso (Burmeister,1993):

Equazione di continuità (conservazione della massa).

La porosità del mezzo (saturo), indicata con f , rappresenta la frazione del volume di controllo che è occupata dal fluido di densità r e capacità termica C_p, mentre la restante parte (1-f ) è occupata dal solido di densità r _s e capacità termica C_ps.

La capacità termica equivalente per unità di volume si esprime in termini di media pesata delle corrispondenti grandezze, utilizzando come pesi le frazioni volumetriche del fluido e del solido:

In modo analogo si definisce anche il termine relativo alla generazione interna di calore e il valore effettivo della conducibilità termica.

Il modello assume che solido e fluido siano sempre in equilibrio termico fra loro, ossia si trovino entrambi alla stessa temperatura T.

Il vettore q indica la portata specifica (m/s) nel punto considerato, ovvero il volume d’acqua fluente nell’unità di tempo attraverso una sezione di area unitaria normale alla direzione del flusso.

K è la conduttività idraulica (m/s) del mezzo poroso, una grandezza che misura la facilità con cui il fluido è trasportato attraverso la matrice solida.

Il modello generale (1-3) va applicato al particolare processo di scambio termico d’interesse. A tal fine occorre considerare le seguenti ipotesi semplificative:

• presenza nel sottosuolo di un acquifero stagnante, privo di qualsiasi componente di velocità su larga scala;
• il moto dell’acqua, governato esclusivamente dalla convezione naturale, è assialsimmetrico (q=q_r i_r+q_z i_z) come il campo termico che lo genera (T=T(r,z));
• la densità dell’acqua viene considerata indipendente dalla pressione, ma variabile con la temperatura T secondo una legge lineare nell’intorno di una temperatura di riferimento T₀:

• i fenomeni termici sono considerati dinamici, mentre il campo di moto nella falda, avendo tempi di assestamento attorno al pozzo di pochi minuti, viene considerato stazionario;
• i termini inerziali dell’equazione (2), assieme alle estensioni di Brinkman e Forchheimer, vengono trascurati poiché i valori delle velocità in gioco sono modesti;
• tutte le propietà fisiche del mezzo poroso (eccetto la densità dell’acqua) vengono considerate costanti nello spazio (omogeneità) ed invarianti rispetto alla direzione (isotropia);

Alla luce delle ipotesi fatte, le equazioni (1-3) divengono:

Per eliminare la pressione dall’equazione del moto, è sufficiente applicare l’operatore vettoriale rotore (curl) ad entrambi i membri dell’equazione (7). La (7) diviene allora:

In conclusione, scrivendo in coordinate cilindriche anche l’equazione di continuità e l’equazione termica, si ottiene:

La formulazione standard del metodo ai residui pesati è impiegata per ottenere una discretizzazione spaziale delle equazioni precedenti.

Opportune funzioni polinomiali N(x), note con il nome di funzioni di forma (shape functions), sono introdotte per interpolare, nello spazio dei punti x del dominio di integrazione W Ì R², i valori nodali della soluzione incognita. Gli n nodi costituiscono i punti, scelti arbitrariamente in W , grazie ai quali viene ottenuta una partizione (o tassellatura) dello spazio in elementi finiti.

In seguito, si farà riferimento esclusivamente agli elementi triangolari lineari, per i quali le funzioni di forma divengono:

dove A^(e) indica l’area del triangolo (e), mentre a_i, b_i, c_i, sono dei coefficienti esprimibili in funzione delle coordinate x_i e y_i dei tre nodi di (e). In notazione vettoriale, i parametri per la costruzione delle N_i divengono:

L’integrazione delle funzioni di forma lineari (e delle loro derivate) è resa particolarmente semplice dalle seguenti formule valide, rispettivamente, per gli integrali di linea lungo uno spigolo di lunghezza L e per quelli doppi all’interno di un triangolo di area A (Pepper, 1992):

La formulazione residuale delle equazioni differenziali (10-12) fornisce le seguenti espressioni:

Rispetto alla (11) si è scelto di orientare l’asse z verso il basso, per cui è stato modificato il segno del gradiente di temperatura che compare nella (17), responsabile dell’effetto ascensionale nelle zone più calde del fluido.

Ponendo dV=2p rdrdz=2p rdW , con W regione rettangolare del piano r-z sede dei fenomeni d’interesse, le (17-19) divengono:

Figura 2: Condizioni al contorno impiegate nel modello: a sinistra le condizioni per il campo termico; a destra quelle relative al campo di velocità. In nero è indicato il pozzo di raggio R0 dotato di tubi ad U multipli o di due tubi coassiali.

Sostituendo alle funzioni incognite q_r , q_z , T le corrispondenti approssimazioni polinomiali (indicate con la soprascritta ^) e utilizzando la formula generalizzata d’integrazione per parti per semplificare il termine diffusivo della (22), si arriva alla cosiddetta formulazione debole (weak statement) del modello:

h_c (W/m²K) indica il coefficiente di convezione sulla superficie del terreno (tratto g _{rÎ ¶ W} ), esposta alla temperatura T_air dell’aria esterna, mentre q_p (W/m²) è il flusso termico scambiato tra fluido e terreno sulla parete del pozzo (tratto g _{zÎ ¶ W} ), orientato dal fluido verso il terreno (cfr. Figura 2).

Esplicitando nelle weak statements il legame tra le incognite e le funzioni di forma, ovvero inserendo nelle equazioni precedenti le espressioni:

si può riformulare il modello in forma matriciale:

Il modello diviene un insieme di n equazioni differenziali ordinarie lineari (a coefficienti variabili con le velocità di falda) scritte per ciascuno degli n nodi del reticolo e accoppiate con altre 2n equazioni lineari algebriche associate al campo di moto del fluido a terreno:

Per eseguire agevolmente gli integrali coinvolti nella definizione delle precedenti matrici, conviene approssimare il raggio r mediante interpolazione lineare tra i valori nodali r con le stesse funzioni di forma N usate per le descrivere l’andamento spaziale delle incognite u,v, T:

r » N(r,z)^T r (32)

Per quanto concerne le funzioni peso W_i, è stata fatta la seguente scelta:

• nelle equazioni del moto (26-27) si pone W_i(r,z)º 1 (W=e), ottenendo così il metodo ai volumi finiti; i vantaggi riguardano soprattutto la (26), che in tal modo diventa una discretizzazione diretta dell’equazione di continuità espressa in forma integrale. Il rispetto del principio di conservazione della massa in un qualsiasi volume di controllo in W , assicura risultati accurati per il campo di velocità.
• nell’equazione termica (28) si pone W_i=N_i, seguendo il metodo standard di Galerkin. Per motivi di stabilità numerica, il termine convettivo rappresentato dalla matrice A(u,v), viene trattato usando il metodo di Petrov-Galerkin (Tezduyar e al.,1990):

Infine, confrontando le (26-28) con le (29-31), si arriva alla determinazione esplicita delle componenti delle matrici del modello:

Per eseguire le integrazioni precedenti sono state usate direttamente le espressioni (15-16). Le matrici M_3x3 ed S_3x3 si ricavano con semplici passaggi algebrici.

Nell’equazione termica (31) compare il flusso scambiato sulla parete del pozzo tra fluido e terreno. Per valutare l’entità di tale flusso, conviene utilizzare un modello semplificato ed esprimere la temperatura d’uscita da un pozzo (contenente tubi ad U multipli o due tubi concentrici) con la seguente relazione (Hellström, 1991):

RH _b (mK/W) è la resistenza effettiva fluido-terreno dello scambiatore, m’ (kg/s) è la portata complessiva che fluisce nello scambiatore, c_pf (J/kgK) è il calore specifico del fluido termovettore. Nell’ipotesi che la temperatura T_p della parete esterna del pozzo sia uniforme, il valore di T_p può essere approssimato con la media della temperatura del terreno lungo il bordo g _z di W (cfr. Figura 2).

La potenza totale scambiata tra fluido e terreno, per unità di superficie, diviene:

Noti i valori di temperatura e portata all’ingresso dello scambiatore, la (48) consente il calcolo del flusso q’_p, assumendo nota la temperatura T_p del terreno a contatto con la parete. Quest’ultima, però, dipende dalla soluzione T del sistema (29-31), che a sua volta è funzione del valore di q’_p. Il metodo più semplice per arrivare alla soluzione, consiste nell’iterare ripetutamente le equazioni del loop algebrico finché una o più variabili di controllo (nel nostro caso T e T_out) convergono entro una prefissata tolleranza.

Per integrare il sistema sparso di equazioni differenziali (31), è stato impiegato lo schema temporale alle differenze finite di Crank-Nicolson:

dove l’apice n indica che le temperature sono valutate all’istante nD t, con D t ampiezza del passo temporale d’integrazione. Sostituendo le (49) nella (31), si ottiene il seguente schema (implicito) di avanzamento temporale:

Nella (50) T_ait dovrebbe essere la media aritmetica dei valori noti al passo n e n+1, ma per semplicità viene mantenuta ad una valore medio costante. Inoltre, si nota che non è l’unico elemento responsabile della formazione di una loop algebrico, perché T dipende (in modo non lineare) anche dalle velocità di falda u e v, le quali a loro volta sono in relazione (lineare) con la temperatura (cfr. (29-30)):

ovvero in forma compatta:

Per risolvere il duplice anello viene impiegato un metodo iterativo diretto che continua a ricalcolare le equazioni accoppiate finché le temperature non convergono entro una specificata tolleranza.

Nel primo esempio proposto è stato preso in considerazione un terreno poroso saturo con livello di falda è H_f=0 e velocità iniziale dell’acquifero nulla in tutti i punti del dominio d’integrazione W.

Gli ingressi del sistema termico sono mantenuti a valori costanti nel tempo:

• portata totale fluente nello scambiatore Q=1 1/s
• temperatura di ingresso del fluido T_in=60°C
• temperatura dell'aria esterna T_ar=20°C

Infine, il passo temporale d’avanzamento è D t=1h con un orizzonte temporale di simulazione di 3 giorni.

La Figura 3 illustra una mappa in falsi colori del campo termico ottenuto dopo 3 giorni di intenso (ed ininterrotto) riscaldamento del suolo. La Figura 4 traccia l’evoluzione del profilo di temperatura all’uscita dallo scambiatore.

Le Figure 5 e 6 mostrano, invece, i risultati relativi al campo di velocità nella falda. L’evidente sviluppo di una corrente ascensionale calda in prossimità del pozzo, è in piena sintonia con il fenomeno atteso di spinta idrostatica verticale (forza di galleggiamento) nelle zone di minore densità (maggiore temperatura) del fluido. Inoltre, nella Figura 6 si notano chiaramente delle celle convettive (moti di circolazione) alimentate dalla corrente ascensionale.

I parametri impiegati in questo secondo test sono gli stessi di quelli appena descritti nell’esempio precedente. L’unica differenza riguarda l’assenza dell’acqua di falda a terreno. Si vuole così valutare l’impatto della falda (e dei connessi effetti convettivi) sulle prestazioni dello scambiatore.

L’ipotesi di pura conduzione a terreno trasforma il problema non lineare dell’esempio 1 (con tre gradi di libertà per nodo: u,v, T) in un problema lineare con un solo grado di libertà per nodo (T): i tempi di calcolo (con un PC486/80) si riducono sensibilmente, passando da 10.2 minuti ad appena 1.25 minuti (generazione del reticolo con 256 nodi inclusa).

La Figura 7 mostra in falsi colori il campo termico puramente conduttivo. La Figura 8 riporta la temperatura in uscita dallo scambiatore in funzione del tempo. Un confronto diretto tra le Figure 3-4 e 7-8 mette in luce alcuni aspetti rilevanti:

• il campo termico conduttivo è caratterizzato da temperature leggermente superiori rispetto a quello convettivo; in quest’ultimo caso si osserva chiaramente come lo strato limite termico (cioè il raggio oltre il quale il terreno è indisturbato) abbia uno spessore più contenuto; la causa è dovuta al moto ascensionale dell’acqua vicino alla parete del pozzo, fenomeno che rende fortemente anisotropo lo scambio termico, privilegiando la direzione verticale piuttosto che quella radiale;
• a regime, la temperatura d’uscita dallo scambiatore convettivo è minore di circa 1°C rispetto al caso conduttivo, confermando l’ipotesi secondo cui un accumulo realizzato in un terreno saturo, a parità di altre condizioni, ha una maggiore capacità di scambio termico. Una differenza unitaria di temperatura, subita da un fluido termovettore che scorre nel terreno con una portata di 1 kg/s si traduce, ricordando il valore della capacità termica dell’acqua (2 kJ/kgK), in un guadagno netto di potenza scambiata con il terreno di 4.2 kW, valutabile attorno al 50% della potenza nominale installata in un impianto residenziale di piccole dimensioni.

In conclusione, si può affermare che l’influenza della falda sulle prestazioni dell’accumulo è decisamente significativa. Dovendo eseguire una stima analitica dei costi operativi di funzionamento di un sistema termico dotato di scambiatori accoppiati a terreno, non sarebbe ammissibile prescindere da un’accurata valutazione degli effetti della convezione termica.

• Burmeister, 1993, "Convective heat transfer", 2° Ed., John Whiley & Sons.
• Deans, Park, Tezduyar, 1990, "Finite elements procedures for time dependent convection-diffusion-reaction systems", Finite Elements in Fluids - Vol. 7, Gallagher, Glowinski, Gresho, Oden, Zienkiewicz ed., McGraw Hill, New York.
• Heinrich J.C., Pepper D.W., 1992, "Finite element method", McGraw-Hill, New York.
• Hellström G., 1991, "Ground heat storage, Vol. I. Thermal analyses of duct storage systems. Theory", Dep. of Mathematical Physics, University of Lund, Sweden.