ComplementiDiIdraulica new 02 WEB - image.unipd.it€¦ · Moto vario in ipotesi anelastiche...

- 1 -

MOTO VARIO NELLE CONDOTTE Alcuni appunti

Andrea Defina

Novembre 2010

Moto vario in ipotesi anelastiche

Prefazione

Sono qui raccolti alcuni appunti relativi agli argomenti trattati nel corso di Complementi di Idraulica, limitatamente alla parte del moto vario nei sistemi in pressione. Si tratta ancora di una BOZZA e, presumibilmente, rimarrà tale in eterno, anche se sarà oggetto di frequenti (sono un ottimista) aggiornamenti. A questo proposito, vi chiedo di segnalarmi gli immancabili svarioni (anche se non siete sicuri che si tratti di uno sbaglio segnalatemelo lo stesso) e di indicarmi le parti che necessitano di maggiore chiarezza (prefazione compresa).

Non vi ringrazio fin d’ora, lo farò di volta in volta.

N.B. Le parti evidenziate in giallo sono buchi che prima o poi intendo riempire

Andrea Defina

c/o Dipartimento di Ingegneria Idraulica, Marittima Ambientale e Geotecnica - Università di Padova. via Loredan, 20 – 35131 PADOVA email: [email protected]


Indice

1 INTRODUZIONE ............................................................................................................. 5

2 MOTO VARIO IN IPOTESI ANELASTICHE ............................................................... 7

2.1 Avviamento di una condotta .............................................................................. 7 2.2 Funzionamento di un serbatoio ....................................................................... 11

2.2.1 Vuotamento di un serbatoio attraverso una lunga condotta ................... 11 2.2.2 Riempimento di un serbatoio ....................................................................... 15

2.3 Oscillazioni di massa ......................................................................................... 16 2.3.1 Oscillazione di massa in un tubo a “U”....................................................... 16 2.3.2 Pozzi piezometrici .......................................................................................... 19 2.3.3 Casse d’aria .................................................................................................... 31 2.3.4 Oscillazioni forzate ........................................................................................ 39

2.4 Colpo d’ariete ....................................................................................................... 49 2.4.1 Riduzione lineare della velocità ................................................................... 50 2.4.2 Chiusura lineare dell’otturatore ................................................................... 51

3 MOTO VARIO IN IPOTESI ELASTICHE .................................................................. 55

3.1 Descrizione qualitativa del fenomeno del colpo d’ariete ......................... 55 3.2 Le equazioni semplificate del moto ................................................................ 63

3.2.1 Forma integrata delle equazioni del moto ................................................. 66 3.2.2 Esempio di applicazione ............................................................................... 67

3.3 Modalità di impiego delle equazioni del moto ............................................. 70 3.3.1 Condizioni al contorno .................................................................................. 71 3.3.2 Intercettazione del flusso in una condotta (colpo d’ariete) ...................... 71 3.3.3 Avviamento del flusso in una condotta (tempo di avviamento) .............. 79 3.3.4 Arresto di una pompa .................................................................................... 81 3.3.5 Arresto di una pompa con cassa d’aria ...................................................... 84

4 RETI DI CONDOTTE .................................................................................................... 89

4.1 Reti di condotte trattate in ipotesi anelastiche ........................................... 89 4.1.1 Condizioni iniziali e al contorno ................................................................... 93 4.1.2 Dispositivi ed organi di regolazione ............................................................ 94

4.2 Reti di condotte trattate in ipotesi elastiche .............................................. 103

Introduzione

- 5 -

1 INTRODUZIONE Il moto che si sviluppa in un sistema di condotte in pressione viene generalmente

studiato mediante l’impiego delle cosiddette equazioni unidimensionali:

jtv

gsE

−∂∂

−=∂∂ 1 (1)

0=∂

∂+

∂∂

tA

sQ ρρ (2)

in cui E è l’energia per unità di peso del fluido, j è la dissipazione di energia per unità di lunghezza misurata lungo il percorso s, v è la velocità, Q è la portata, A è l’area della sezione trasversale, ρ è la densità del fluido, g è l’accelerazione di gravità, s è l’ascissa curvilinea e t è il tempo.

L’equazione dinamica (1) contiene già una semplificazione, comunemente adottata, che consiste nell’aver trascurato, nel termine E, gli effetti legati alle variazioni spaziali della densità (A questo proposito si veda la nota, riportata in calce, della prima pagina del capitolo 3.2).

Nel caso più generale, le variabili dipendenti sono la pressione p (contenuta nell’espressione per energia E), la velocità, la densità e l’area della sezione trasversale. Queste quattro grandezze saranno, in generale, dipendenti dallo spazio (s) e dal tempo (t). Il sistema composto dalle equazioni (1) e (2) non è quindi completo e altre relazioni (equazioni costitutive) vanno aggiunte per rendere determinato il problema.

In particolare si assume che il moto sia barotropico, ovvero che la densità sia funzione della sola pressione: ρ=ρ(p) e che il condotto sia deformabile elasticamente per effetto della pressione per cui A=A(p). Si parla in tal caso di “ipotesi elastiche”.

Spesso, però, è possibile assumere, con buona approssimazione, costante nel tempo e nello spazio la densità del fluido e costante nel tempo l’area della sezione trasversale. In tal caso le variabili del problema si riducono alla pressione e alla velocità (funzioni del tempo e dello spazio) e le equazioni (1) e (2), a cui vanno associate le opportune condizioni al contorno, sono sufficienti a rendere determinato il problema. Si parla in tal caso di “ipotesi anelastiche”.

Non è facile enunciare le regole generali che consentano, a priori, di stabilire se un certo problema possa essere trattato in ipotesi anelastiche, oppure richieda la soluzione delle equazioni complete. Solo l’esperienza e la conoscenza qualitativa di come si sviluppa nella realtà il moto, possono indirizzare verso la scelta corretta.


- 7 -

2 MOTO VARIO IN IPOTESI ANELASTICHE La soluzione dei problemi in ipotesi anelastiche è sensibilmente più agevole

rispetto alla trattazione completa. Pertanto, ogni qual volta il problema da risolvere consente di invocare tali ipotesi conviene, almeno in via preliminare, ricorrere ad esse.

In ipotesi anelastiche le equazioni del moto, assunta la densità costante nel tempo e nello spazio e la generica sezione trasversale della condotta costante nel tempo, diventano

jtv

gsE

−∂∂

−=∂∂ 1 (3)

0=∂∂

sQ (4)

Nel seguito sono illustrati alcuni problemi che possono essere trattati e risolti in ipotesi anelastiche. Si tratta di semplici esempi che consentono, tra l’altro, di introdurre alcuni parametri caratteristici in grado di qualificare sinteticamente il moto vario.

2.1 Avviamento di una condotta Si consideri il sistema schematicamente illustrato in Figura 2.1 costituito da un

serbatoio e da una lunga condotta di scarico munita, al termine, di un rubinetto (R) inizialmente chiuso. Immaginiamo di aprire istantaneamente e completamente il rubinetto. Il moto, all’interno della condotta si avvia e la velocità cresce, a partire dal valore nullo iniziale, raggiungendo un valore di regime v0 costante se il livello nel serbatoio si dovesse mantenere anch’esso costante nel tempo.

Figura 2.1

Nella realtà, durante l’avviamento della condotta, il serbatoio gradualmente si vuota. E’ possibile però immaginare che la superficie del serbatoio sia sufficientemente grande e la durata del fenomeno di avviamento, talmente breve da poter trascurare la variazione di livello nel serbatoio.

Il primo passo, nella soluzione dell’equazione (3), consiste nell’integrazione della stessa lungo il percorso tra i punti 1 (serbatoio) e 2 (sbocco).


∫∫∫ −∂∂

−=∂∂ 2

1

2

1

2

1

1 dsjdstv

gds

sE (5)

Essendo diverse le velocità (e quindi anche le derivate temporali della velocità) nei tratti 1-1’ e 1’-2, conviene spezzare i due integrali a destra come segue

∫∫∫∫ −−∂∂

−∂∂

−=−2

'1

'1

1

2

'1

'1

112

11 dsjdsjdstv

gds

tv

gEE (6)

E’ facile mostrare1 che sia la velocità, sia il termine ∂v/∂t sono piccoli lungo il tratto 1-1’ cosicché i corrispondenti integrali possono essere trascurati. Lungo il tratto 1’-2, inoltre, sia l’accelerazione temporale ∂v/∂t che le dissipazioni di energia per unità di lunghezza j sono costanti2 e possono essere portate fuori dall’integrale. Non dipendendo la velocità dallo spazio ma solo dal tempo, la derivata parziale rispetto al tempo può essere sostituita con la derivata totale, ottenendo

Ljdtdv

gLEE −−=− 12 (7)

in cui L è la lunghezza della condotta e v è la velocità in condotta, variabile nel tempo.

Nel serbatoio è E1=h0, costante per le ipotesi fatte. Allo sbocco, a partire dall’istante in cui viene completamente aperto il rubinetto, è E2=v2/2g. Le dissipazioni di energia, continue e localizzate, possono essere espresse come proporzionali al carico cinetico. Si può porre quindi: jL=r.v2/2g. L’equazione (7) può essere pertanto riscritta come segue

gvr

dtdv

gLh

gv

22

2

0

2

−−=− (8)

ovvero

gvr

dtdv

gLh

2)1(

2

0 ++= (9)

Avendo ipotizzato costante il livello nel serbatoio, la precedente equazione differenziale presenta una sola variabile dipendente dal tempo, la velocità. E’ sufficiente quindi, trattandosi di un’equazione di primo ordine, una sola condizione al contorno. Tale condizione è: per t=0, v=0. 1 Dall’equazione di continuità, essendo ∂Q/∂s=0, sarà anche ∂/∂t[∂Q/∂s]=0. Invertendo l’ordine di derivazione si ha: ∂/∂s[∂Q/∂t]=0. Essendo l’area costante nel tempo si ha: ∂/∂s[A ∂v/∂t]=0. Cioè il prodotto A ∂v/∂t è costante lungo s. Lungo il tratto 1-1’ la sezione trasversale è molto più grande di quella della condotta (tratto 1’-2), di conseguenza ∂v/∂t lungo il tratto 1-1’ sarà trascurabilmente piccolo rispetto allo stesso termine lungo la condotta. Sempre dall’equazione di continuità che sancisce la costanza della portata lungo s, essendo Q=vA, si deduce che nel tratto 1-1’ le velocità sono estremamente basse e trascurabilmente piccole quindi le dissipazioni di energia. 2 Dalla nota precedente si ha: ∂/∂s[A ∂v/∂t]=0. Lungo la condotta, inoltre, l’area A è costante. Si ha quindi: ∂/∂s[∂v/∂t]=0, ovvero ∂v/∂t è costante lungo s. Esprimendo j mediante la formula di Darcy-Weisbach, ad esempio, si ha:j=f/d(v2/2g) in cui diametro d e velocità v (per continuità) sono costanti lungo s. Costante lungo s è anche la funzione di resistenza f che dipende dalla scabrezza e dal numero di Reynolds (cioè da velocità e diametro).


Di un qualche interesse è la soluzione della precedente equazione nella forma semplificata che approssima la soluzione corretta nei primissimi istanti del fenomeno. In un intorno di t=0 le velocità sono ancora molto piccole e risulta quindi trascurabile il termine proporzionale al carico cinetico. In queste ipotesi semplificative l’equazione (9) diventa

dtdv

gLh =0 (10)

L’equazione (10) è a variabili separabili. La soluzione è:

thgvLt cos

0

+= (11)

Con la condizione al contorno prima indicata, la costante di integrazione è nulla. Il tempo necessario all’avviamento della condotta, detto “tempo di avviamento”, Ta, si trova dalla (11) sostituendo alla velocità il suo valore di regime v0.

La velocità risulta pertanto linearmente variabile nel tempo e data dalla seguente relazione:

0

0

0 hgvL

TTt

vv

aa

== (12)

Il tempo di avviamento Ta non è ovviamente il tempo necessario al moto per avviarsi completamente. La soluzione espressa dalla (12), infatti, è stata ottenuta sotto ipotesi semplificative pesanti. Tuttavia il parametro Ta è una misura che dà una precisa idea della durata che caratterizza questo fenomeno.

Consideriamo ora la soluzione dell’equazione (9) completa. Nelle condizioni finali, di regime, la velocità in condotta sarà v=v0, mentre sarà nullo il termine che esprime la variazione temporale di velocità. Nelle condizioni di regime, quindi, l’equazione (9) si riduce alla seguente

gvrh2

)1(20

00 += (13)

E’ da osservare che il coefficiente r che compare nella (9), in generale dipende dal numero di Reynolds e quindi dalla velocità che a sua volta varia nel tempo. Pertanto, il valore r0 di regime sarà diverso dal valore istantaneo durante l’evolversi del fenomeno. Possiamo però assumere, con buona approssimazione, che sia r=r0. In questa ipotesi, esplicitata l’equazione (13) rispetto a (1+r0) e sostituita questa espressione al posto del termine (1+r) che compare nella (9), si trova:

20

2

00 vvh

dtdv

gLh += (14)

La precedente equazione può essere riorganizzata come segue:

dtdv

vhgvL

vv

00

020

2 11 =− (15)

e, ricordando l’espressione per il tempo di avviamento Ta, si può scrivere:


20

0

)/(1)/()/(

vvvvdTtd a −

= (16)

L’equazione differenziale (16), con la precedente condizione al contorno, presenta la seguente soluzione:

⎟⎟⎠

⎞⎜⎜⎝

⎛=

aTt

vv tanh

0

(17)

La Figura 2.2 illustra l’andamento della soluzione espressa dalla (17) e, sovrapposta, la soluzione lineare (12). Si osserva che nei primi istanti, cioè quando la velocità è trascurabilmente piccola, le due soluzioni sono poco dissimili tra loro. Si può osservare inoltre che pur essendo a rigore infinito il tempo necessario all’avviamento della condotta, già per t=3Ta la velocità è praticamente indistinguibile da quella di regime essendo v≈0.995.v0.

Figura 2.2


2.2 Funzionamento di un serbatoio Il funzionamento di un serbatoio è sostanzialmente governato dall’equazione di

continuità nella forma detta “equazione dei serbatoi” espressa dalla relazione (19). In questa equazione la portata entrante (Qentr) e quella uscente (Qusc) sono generalmente variabili nel tempo con andamenti che dipendono dalle modalità di alimentazione e di scarico del serbatoio. In questo capitolo si riportano solo alcuni esempi. Nei due casi che seguono (paragrafi 2.2.1 e 2.2.2) si assume che lo scarico del serbatoio avvenga mediante una lunga condotta. E’ questa una modalità di scarico che si incontra frequentemente nella pratica e che, dal punto di vista didattico, consente di fare qualche interessante considerazione sull’importanza relativa dei diversi termini che compaiono nell’equazione dinamica (3).

2.2.1 Vuotamento di un serbatoio attraverso una lunga condotta

Un problema in qualche misura associato al precedente (avviamento di una condotta) riguarda il fenomeno di vuotamento di un serbatoio. Facciamo riferimento per semplicità alla stessa Figura 2.1. Ripetendo l’integrazione spaziale dell’equazione del moto ritroviamo l’equazione (9) nella quale, questa volta, il livello del serbatoio sarà variabile nel tempo mentre, con buona approssimazione possiamo considerare costante il parametro r. L’equazione dinamica che governa il processo di vuotamento per il sistema considerato è quindi:

gvr

dtdv

gLh

2)1(

2

++= (18)

In questa situazione le variabili dipendenti del problema sono due: la velocità in condotta e il livello del serbatoio. E’ necessario quindi associare alla precedente, una seconda equazione contenente almeno una delle due variabili. L’equazione che si introduce è quella di continuità nella forma detta “equazione dei serbatoi”:

dtdhhQQ uscentr )(Ω=− (19)

nella quale Qentr e Qusc sono le portate entranti ed uscenti dal serbatoio e Ω è l’area della superficie orizzontale del serbatoio al generico livello h.

Nel caso in esame si ha Qentr=0, Qusc=v.A essendo A l’area della sezione trasversale della condotta. L’equazione (19), può quindi essere riscritta, esplicitata rispetto alla velocità, come:

dtdh

Ahv )(Ω

−= (20)

Il sistema composto dalle equazioni (18) e (20) è non lineare e la sua soluzione analitica è di fatto impraticabile. Nel caso di vuotamento di un serbatoio, fatta eccezione per la breve fase iniziale di avviamento del moto, vista nel paragrafo precedente, le variazioni temporali di velocità sono spesso modeste. Il termine di accelerazione temporale può quindi spesso essere trascurato.

In questa ipotesi, l’equazione (18), esplicitata rispetto alla velocità, diventa:


2/1

12 h

rgv ⋅

+= (21)

Uguagliando tra loro le equazioni (20) e (21), per eliminare la variabile velocità, si trova:

dtr

gAdhh

h⋅

+−=

Ω12)(

2/1 (22)

La precedente equazione può agevolmente essere integrata una volta nota la legge Ω(h). Nel caso particolare di serbatoio cilindrico (Ω indipendente dal livello h), L’equazione (22) presenta una semplice soluzione analitica:

ttr

gAh cos12

22/1 +⋅

+Ω−= (23)

La costante di integrazione si determina utilizzando la condizione al contorno: per t=0 è h=h0. Si trova così

2

0 12

2 ⎥⎦

⎤⎢⎣

⎡⋅

+Ω−= t

rgAhh (24)

Il completo vuotamento del serbatoio si determina dalla (24) ponendo h=0. Osservando, dalla (21) che è )1/(2 00 rghv += , il tempo di vuotamento, tv, per questo particolare caso, vale:

0

02vh

Atv

Ω= (25)

Esempio. Consideriamo il sistema illustrato in Figura 2.1. La condotta è lunga 2000m, ha un diametro di 0.5m ed è caratterizzata da un coefficiente di scabrezza secondo Gauckler-Strickler di 80 m1/3/s. Il serbatoio è cilindrico con una superficie orizzontale Ω=1000 m2. Il livello iniziale nel serbatoio è h0=50 m. Immaginando un’apertura istantanea del rubinetto, valutiamo la legge di vuotamento del serbatoio.

Con i dati del problema, il coefficiente r=r0 vale

98.123/420 ==

Hs RKLgr

Utilizzando l’equazione (13) si trova v0=3.146 m/s. Utilizzando infine le equazioni (12) e (25) si trova che il tempo di avviamento vale Ta=12.8 s e il tempo di vuotamento vale approssimativamente tv=45 ore. E’ da sottolineare, con riferimento ai risultati trovati, che la velocità nella condotta passa da zero a v0 in poco più di mezzo minuto mentre, nella successiva fase di vuotamento, passa da v0 a zero in quasi due giorni. E’ evidente la trascurabilità del termine di accelerazione temporale in questa seconda fase.

Se anche proviamo a cambiare i dati del problema, riducendo la lunghezza della condotta a L=50 m e limitando l’area della superficie orizzontale del serbatoio Ω ad


appena 10 m2, si trova Ta=1.7 s e tv=302 s. Il processo di vuotamento risulta essere ancora molto lento.

Soluzione numerica del problema di avviamento e vuotamento. Nel seguito si illustra una semplice soluzione numerica del problema utilizzando un metodo di integrazione di tipo esplicito.

Le equazioni (18) e (20) sono riscritte come segue

⎥⎦

⎤⎢⎣

⎡+−=

gvrh

Lg

dtdv

2)1(

2

)(hAv

dtdh

Ω−=

(26)

L’asse dei tempi è diviso in intervalli Δt tali per cui

nktkt ,,3,2,1 K=Δ=

La generica variabile all’istante t, è indicata con l’apice k. Ad esempio h(t) è indicato con hk. Discretizzando le due derivate che compaiono nelle (26) alle differenze finite, e con riferimento ad un metodo di integrazione esplicito,si trova

kkk

gvrh

Lg

tvv

⎥⎦

⎤⎢⎣

⎡+−=

Δ−+

2)1(

21

kkk

hAv

thh

⎥⎦

⎤⎢⎣

⎡Ω

−=Δ−+

)(

1

(27)

Le precedenti relazioni possono essere riscritte come segue: k

kk

gvrh

Lgtvv ⎥

⎦

⎤⎢⎣

⎡+−Δ+=+

2)1(

21

kkk

hAvthh ⎥

⎦

⎤⎢⎣

⎡Ω

Δ−=+

)(1

(28)

Note le condizioni al contorno: per t=k=0 è v=0 e h=h0, le precedenti relazioni, utilizzate ricorsivamente, consentono di ricostruire la soluzione h(t) e v(t) per passi temporali discreti Δt. Ovviamente deve essere preventivamente scelto il passo di integrazione Δt. Nel caso in esame non vi sono particolari problemi di stabilità e la scelta va effettuata considerando la qualità, ovvero l’accuratezza che si vuole caratterizzi la soluzione. Si potrebbe pensare di assumere un passo temporale pari a 0.1÷0.2 Ta. La soluzione, con i dati iniziali del problema, è illustrata in Figura 2.3.


Figura 2.3 – Caso 1:L=2000m, Ω=1000m2, Δt=2s

Nel grafico dei livelli si può osservare come, durante la fase di avviamento che dura meno di 40 secondi, il livello nel serbatoio scende di pochi centimetri e l’ipotesi di assumere h=h0 risulta più che accettabile. Con riferimento al processo di vuotamento, l’ipotesi di trascurare la fase iniziale di avviamento risulta, sia per i livelli che per le velocità, certamente adeguata.

La soluzione con la seconda seri di dati è illustrata in Figura 2.4.

Figura 2.4 – Caso 2:L=50m, Ω=10m2, Δt=0.2s

In questo caso il processo di vuotamento è molto rapido. Durante la fase di avviamento, che dura approssimativamente 5 secondi, il livello nel serbatoio scende di circa 1.5 m (corrispondente al 3% del carico totale). Ciò nonostante l’andamento delle velocità predetto nell’ipotesi di livello costante risulta abbastanza accurato. Il processo di vuotamento, nell’ipotesi di trascurare la fase iniziale di avviamento è meno accurato, con riferimento ai livelli, rispetto al caso precedente. Tuttavia la soluzione può ancora ritenersi accettabile.


2.2.2 Riempimento di un serbatoio Il problema opposto al precedente è rappresentato dal processo di riempimento di

un serbatoio. Facciamo riferimento per semplicità alle condizioni illustrate nella stessa Figura 2.1 per le quali, però, si assume che vi sia una portata Qentr, non nulla in ingresso al serbatoio. Assumiamo inoltre che il livello iniziale nel serbatoio sia relativamente piccolo (ma superiore alla quota di imbocco della condotta) e tale per cui la portata scaricata mediante la condotta sia inferiore a quella entrante.

Ovviamente, per la soluzione di questo problema valgono le equazioni (18) e (19). In particolare, come per il caso precedente, possiamo trascurare il termine che descrive l’accelerazione temporale nell’equazione (18) e scrivere inoltre Qusc=v.A. Si perviene così al seguente sistema di equazioni

gvrh2

)1(2

+= (29)

dtdhhAvQentr )(Ω=− (30)

Assumiamo Ω indipendente dal livello e consideriamo una portata entrante costante nel tempo. Esplicitata l’equazione (29) rispetto alla velocità e sostituita l’espressione di quest’ultima nell’equazione (30), si trova

rhgAQ

dtdh entr

+Ω−

Ω=

12 (31)

Si tratta di un’equazione differenziale a variabili separabili la cui soluzione è

( ) trhgA

rhg

QAQ

Agrt

entrentr cos

12

1211

2 +⎥⎥⎦

⎤

⎢⎢⎣

⎡

++⎟

⎟⎠

⎞⎜⎜⎝

⎛

+−

Ω+−= ln (32)

La costante di integrazione che compare nella precedente relazione si determina imponendo le condizioni al contorno: per t=0 si ha h=h0. Si trova così

( ) ( )( )

( )⎥⎥⎦

⎤

⎢⎢⎣

⎡−

++

+−

+−Ω+−= hh

rgA

rhgAQrhgAQ

QAgrt

entr

entrentr 0

02 1

21/21/21 ln (33)

Il massimo livello nel serbatoio si raggiunge quando il sistema è a regime e può essere agevolmente determinato utilizzando l’equazione (31) nella quale si ponga dh/dt=0. Risulta

2

21

⎟⎠

⎞⎜⎝

⎛+=

AQ

grh entr

MAX (34)


2.3 Oscillazioni di massa Il fenomeno dell’oscillazione di massa si presenta con notevole frequenza nei

sistemi in pressione. Si tratta di un fenomeno che riveste notevole importanza dal punto di vista ingegneristico e, nella maggior parte dei casi, può essere trattato in ipotesi anelastiche.

Per inquadrare il fenomeno esaminiamo inizialmente un semplice esempio. 2.3.1 Oscillazione di massa in un tubo a “U”

Consideriamo il sistema illustrato in Figura 2.5 costituito da una tubazione a “U” di diametro costante contenente un liquido. In condizioni di equilibrio, il livello nei due rami del tubo si pone alla stessa quota (Figura 2.5-sinistra). Assumiamo questa quota come livello di riferimento.

Figura 2.5

Una qualsiasi perturbazione di questa situazione di equilibrio innesca un’oscillazione di massa che, per effetto delle dissipazioni di energia, si smorza nel tempo fino a ritornare alle condizioni di equilibrio iniziali.

Consideriamo, per semplicità, una particolare perturbazione. Tappiamo l’estremità del ramo di destra del tubo e creiamo molto lentamente una depressione p0 (vedi Figura 2.5-destra). In questo modo il livello nel ramo di destra si alza della quantità ΔH=-p0/2γ, essendo γ il peso specifico del liquido. Ovviamente il livello nel ramo di sinistra si abbassa della stessa quantità. Anche la nuova condizione è di equilibrio statico.

A partire da questa situazione eliminiamo rapidamente il tappo dall’estremo del ramo di destra ripristinando istantaneamente condizioni di pressione atmosferica. La situazione che si viene a determinare è ora squilibrata ed innesca l’oscillazione della colonna di liquido contenuta nel tubo.

Per studiare le caratteristiche di questa oscillazione, facciamo riferimento alle equazioni del moto vario in ipotesi anelastiche.


Effettuando l’integrazione spaziale dell’equazione del moto ritroviamo l’equazione (7) nella quale E2=z+v2/2g, E1=-z+v2/2g e jL=r.v v /2g. L’equazione (7) può essere pertanto riscritta come segue

gvv

rdtdv

gLz

22 −−= (35)

Da osservare che potendo la velocità essere positiva o negativa (il moto cambia verso nel tempo) è necessario utilizzare il modulo della velocità nel termine che esprime le dissipazioni di energia. Come per il problema di vuotamento di un serbatoio, le variabili dipendenti del problema sono due: la velocità v nel tubo e il livello z nel ramo di destra. E’ necessario quindi introdurre una seconda equazione contenente almeno una delle due variabili. L’equazione che si introduce è quella di continuità:

dtdzv = (36)

Per ottenere un’agevole soluzione analitica del sistema composto dalle equazioni (35) e (36), assumiamo di trascurare le dissipazioni di energia. In tal caso, sostituita l’espressione della velocità data dalla (36) nell’equazione (35), quest’ultima può essere riscritta nella forma

022

2

=+ zLg

dtzd (37)

Si tratta di un’equazione differenziale lineare a coefficienti costanti che presenta un’agevole soluzione analitica. Posto ω2=2g/L, l’equazione (37) è riscritta come

022

2

=+ zdt

zd ω (38)

Che presenta la soluzione generale:

tCtsenCz ωω cos21 += (39)

Le costanti di integrazione C1 e C2 vanno determinate imponendo le condizioni al contorno che sono: per t=0, z=ΔH e per t=0, v=0.

Dalla prima condizione si ha:

HCtCtsenCztt

Δ==+=== 20210

cosωω (40)

Dalla seconda, si ha

0cos 10210

0==−==

==

=CtsenCtC

dtdzv

tt

tωωωωω (41)

La soluzione è quindi

tLgHtHz 2coscos Δ=Δ= ω (42)


Utilizzando la (36), si ha inoltre

tLgsen

LgHv 22

Δ−= (43)

Si tratta di una oscillazione di ampiezza ΔH e periodo gLT 2/2π= .

Soluzione numerica del problema di oscillazione in un tubo a “U”. Nel seguito si illustra una semplice soluzione numerica del problema utilizzando un metodo di integrazione di tipo esplicito. E’ da osservare che la soluzione con un metodo esplicito, in assenza di dissipazioni di energia, per questo sistema di equazioni, risulta essere instabile. In presenza di dissipazioni di energia la stabilità può essere assicurata scegliendo un passo di integrazione temporale sufficientemente piccolo.


⎥⎥⎦

⎤

⎢⎢⎣

⎡+−=

gvv

rzLg

dtdv

22

vdtdz

=

(44)


nktkt ,,3,2,1 K=Δ=

La generica variabile all’istante t, è indicata con l’apice k. Ad esempio z(t) è indicato con zk. Discretizzando le due derivate che compaiono nelle (46) alle differenze finite, e con riferimento ad un metodo di integrazione esplicito,si trova

kkk

gvv

rzLg

tvv

⎥⎥⎦

⎤

⎢⎢⎣

⎡+−=

Δ−+

22

1

kkk

vt

zz=

Δ−+1

(45)


kk

gvv

rzLgtvv

⎥⎥⎦

⎤

⎢⎢⎣

⎡+Δ−=+

221

kkk vtzz Δ+=+1

(46)

Note le condizioni al contorno: per t=k=0 è v=0 e z=ΔH, le precedenti relazioni, utilizzate ricorsivamente, consentono di ricostruire la soluzione z(t) e v(t) per passi temporali discreti Δt.

Assumiamo, ad esempio, che sia L=10 m, d=0.1 m, ΔH=0.5 m e che la tubazione sia caratterizzata da un coefficiente di scabrezza secondo Gaukler-Strickler di 50 m1/3/s.


In tal caso, il parametro relativo alle dissipazioni vale r=10.74. Il periodo teorico delle oscillazioni in assenza di dissipazioni vale T=4.5 s. Si può quindi ragionevolmente fissare un passo di integrazione temporale Δt=0.05 s.

La soluzione numerica e quella teorica in assenza di dissipazioni, espressa dalle equazioni (42) e (43), sono confrontate in Figura 2.6.

Figura 2.6 – Linea continua: soluzione analitica, linea tratteggiata:soluzione numerica con dissipazioni.

Si osserva come la soluzione analitica semplificata (ipotesi di dissipazioni nulle) preveda correttamente il periodo dell’oscillazione e, con buona approssimazione, anche il valore massimo dell’ampiezza dell’oscillazione del livello e della velocità. La soluzione numerica evidenzia l’effetto delle dissipazioni che producono un’oscillazione di tipo smorzato.

2.3.2 Pozzi piezometrici

Consideriamo il sistema illustrato in Figura 2.7 che schematizza una comune derivazione idroelettrica in pressione. Il sistema è costituito da un serbatoio di alimentazione, da una condotta o galleria di derivazione, spesso molto lunga, da un pozzo piezometrico e da una condotta forzata che adduce alla centrale.

Figura 2.7


Per effetto della regolazione dell’otturatore, posto al termine della condotta forzata, legata alle esigenze di produzione della centrale, si creano, in corrispondenza dell’otturatore stesso, delle variazioni di pressione che possono essere anche sensibilmente elevate (vedi paragrafo 2.4 e capitolo 1). Queste onde di pressione si propagano rapidamente lungo la condotta forzata e, in assenza del pozzo piezometrico, andrebbero ad interessare l’intera galleria di derivazione, costringendo, per quest’ultima, a dimensionamenti costosi. Per questo motivo, al termine della galleria di derivazione viene posto un pozzo piezometrico che di fatto è un dispositivo di protezione della galleria assorbendo gran parte delle sovrapressioni provenienti dalla condotta forzata.

Le equazioni che governano l’oscillazione di massa nel pozzo piezometrico sono ancora quelle per il moto vario unidimensionale in ipotesi anelastiche. Effettuando l’integrazione spaziale dell’equazione del moto tra le sezioni 1 e 3 indicate in Figura 2.7, si ha

∫∫∫∫∫∫ −−−∂∂

−∂∂

−∂∂

−=−3

2

2

'1

'1

1

3

2

2

'1

'1

113

111 dsjdsjdsjdstv

gds

tv

gds

tv

gEE (47)

nella quale l’integrale dei termini di inerzia e dissipativi è stato spezzato nei tratti 1-1’, 1’-2 e 2-3.

Nella precedente relazione sono certamente trascurabili i contributi relativi al tratto 1-1’. Con riferimento alla configurazione schematizzata in Figura 2.7, inoltre, si può ammettere che la lunghezza del pozzo e le velocità al suo interno siano modeste e risulti pertanto trascurabile il contributo inerziale lungo il tratto 2-3. Assunto il riferimento indicato in Figura 2.7 e trascurato il carico cinetico nel pozzo si ha E3=z. Si può inoltre assumere che l’estensione del serbatoio sia notevole e risultino pertanto trascurabili le variazioni di quota della sua superficie. Si assume cioè E1=0. In queste ipotesi, l’equazione (47) può essere riscritta come segue

32−Δ−−−= EjLdtdv

gLz (48)

in cui v è la velocità in galleria e ΔE2-3 indica il complesso di perdite di energia che si realizzano nel pozzo. Nella pratica, tali perdite non sono trascurabili e vengo spesso appositamente create per smorzare rapidamente le oscillazioni di massa.

All’equazione (48), nelle variabili z e v, va associata l’equazione di continuità che, nel caso in esame, è un equazione di nodo (nodo 2 di Figura 2.7):

ppfg AvQAv += (49)

nella quale Ag e Ap sono le aree delle sezioni della galleria e del pozzo rispettivamente, vp è la velocità nel pozzo e Qf e la portata che fluisce attraverso la condotta forzata. Potendo scrivere vp=dz/dt, e assumendo noto l’andamento nel tempo della portata Qf scaricata attraverso la condotta forzata, il sistema composto dalle equazioni (48) e (49), con le necessarie condizioni al contorno, può essere risolto. A tale proposito, si vedrà più avanti un esempio di soluzione numerica del problema.

Per ottenere una semplice soluzione analitica di questo problema assumiamo di trascurare tutte le dissipazioni di energia. Dalla (49) si ha


fgg

p QAdt

dzAA

v 1+= (50)

Sostituita la (50) nella (48), quest’ultima diventa

dtdQ

AgL

dtzd

AA

gLz f

gg

p −−= 2

2

(51)

Ancora per semplicità ammettiamo che la manovra all’otturatore che innesca l’oscillazione di massa sia una variazione istantanea della portata Qf che si realizzi all’istante t=0, dalle condizioni precedenti Qf=v0Ag (v0 è la velocità in galleria all’istante t=0) a quelle valide per t>0: Qf=Qf0. Con questa semplificazione, per t>0 risulta dQf/dt=0 e l’equazione (51) può essere ulteriormente semplificata.

Posto

p

g

ALAg

=ω (52)

L’equazione differenziale (51), a coefficienti costanti, presenta la soluzione generale

)cos()( tBtsenAz ωω += (53)

in cui le costanti A e B vanno determinate utilizzando le condizioni al contorno. All’istante t=0, il livello nel pozzo è quello indisturbato z=0. Dalla (53) si ha dunque

BBsenA =+= )0cos()0(0 (54)

Ancora per t=0 vale, ovviamente l’equazione di continuità (50).

g

f

tg

p

AQ

dtdz

AA

v 0

00 +=

=

(55)

Derivando la (53) e sostituendo tale derivata nella (55) si trova A=(Q0-Qf0)/ωAp, in cui Q0=v0

.Ag è la portata fluente in galleria all’istante t=0. La soluzione, pertanto è

)(00 tsenAQQz

p

f ωω

−= (56)

Il coefficiente A risulta quindi essere l’ampiezza delle oscillazioni di livello nel pozzo piezometrico e dipende dalle caratteristiche geometriche del sistema e dalla variazione di portata Q0-Qf0 che si determina attraverso la condotta forzata.

E’ possibile in questo modo verificare che il livello massimo che si raggiunge nel pozzo sia contenuto all’interno della vasca e che il livello minimo non determini lo scoprimento della base del pozzo causando fastidiosi ingressi d’aria nel sistema. Nella pratica, le verifiche sono più severe ipotizzando particolari successioni di manovre di apertura e chiusura dell’otturatore tali da produrre valori estremi per i livelli, più penalizzanti rispetto a quelli che si hanno con una manovra singola.


Utilizzando l’equazione di continuità (50) è possibile infine valutare la legge con cui varia nel tempo la velocità (e quindi la portata) nella galleria di derivazione. Si trova

g

f

g

f

AQt

AQQv 000 )cos( +

−= ω (57)

Soluzione numerica del problema di oscillazione in un pozzo piezometrico. Nel seguito si illustra una semplice soluzione numerica del problema utilizzando un metodo di integrazione di tipo esplicito. E’ da osservare, come per i casi precedenti, che la soluzione con un metodo esplicito, in assenza di dissipazioni di energia risulta essere instabile. In presenza di dissipazioni di energia la stabilità può essere assicurata scegliendo un passo di integrazione temporale sufficientemente piccolo.


( )32−Δ++−= ELjzLg

dtdv (58)

p

f

p

g

AQ

AAv

dtdz

−= (59)


nktkt ,,3,2,1 K=Δ=

La generica variabile all’istante t, è indicata con l’apice k. Ad esempio z(t) è indicato con zk. Discretizzando le due derivate che compaiono nelle equazioni (58) e (59) alle differenze finite, e con riferimento ad un metodo di integrazione esplicito,si trova:

[ ]kkk

ELjzLg

tvv

32

1

−

+

Δ++−=Δ

−

p

kfg

kkk

AQAv

tzz −

=Δ

−+1

(60)

le quali possono essere esplicitate rispetto a vk+1 e zk+1. Il sistema è completato dalle condizioni iniziali e al contorno per il particolare

problema che si vuole studiare. Consideriamo ad esempio (Figura 2.7) una galleria di derivazione a sezione

circolare di diametro d=4.0 m e lunga L=2500 m. Assumiamo, come riferimento per le quote, il livello nel serbatoio di alimentazione. Ipotizziamo inoltre, per semplicità, che le dissipazioni di energia lungo la galleria di derivazione possano essere calcolate mediante la formula di Darcy-Weisbach nella quale la funzione di resistenza assuma un valore costante e pari a f=0.02, mentre le dissipazioni alla base del pozzo piezometrico siano espresse dalla relazione pp QQkE 023 =Δ con k0=0.01 s2/m5,


essendo Qp=vpAp la portata che entra nel pozzo e Ap=45 m2 l’area di una generica sezione orizzontale del pozzo.

Consideriamo dapprima il caso in cui a partire da Qf=0, si operi una manovra di apertura istantanea a seguito della quale risulti Qf=10 m3/s costante per t>0. In tal caso le condizioni iniziali sono rappresentate da: z=0, v=vp=0.

La soluzione, ottenuta utilizzando un passo di integrazione temporale t=1.0s è illustrata in Figura 2.8. Nella stessa figura, per confronto, è riportata la soluzione analitica ricavata in precedenza ed espressa dall’equazione (56) nella quale si pone Q0=0 e Qf0=10 m3/s.

Si osservi, in Figura 2.8, che in condizioni di regime si ha: z=-0.40 m (z=-JL), vp=0 e, ovviamente, Q=10.0 m3/s.

Figura 2.8 – Manovra di apertura: andamento nel tempo del livello nel pozzo piezometrico

(considerando e trascurando la dissipazione localizzata alla base del pozzo) e confronto con la soluzione analitica.

A partire da queste condizioni di moto stazionario (condizioni iniziali) è stata simulata una manovra di chiusura. In questo caso le condizioni al contorno, rappresentate dall’andamento nel tempo delle portate verso la condotta forzata, sono descritte dalle seguenti relazioni: Qf=10.0 m3/s per t=0 e Qf=0 per t>0.

I risultati del calcolo effettuato utilizzando un passo di integrazione temporale Δt=1.0s sono illustrati in Figura 2.9. Nella stessa figura, per confronto, è riportata la soluzione analitica ricavata in precedenza ed espressa dall’equazione (56) nella quale si pone Q0=10 e Qf0=0 m3/s.


Figura 2.9 – Manovra di chiusura: andamento nel tempo del livello nel pozzo piezometrico

(considerando e trascurando la dissipazione localizzata alla base del pozzo) e confronto con la soluzione analitica. I simboli sono gli stessi della precedente Figura 2.8.

Si osserva che l’ampiezza delle oscillazioni si riduce molto più lentamente rispetto al caso precedente. Ciò è dovuto al fatto che le portate lungo la galleria di derivazione oscillano attorno ad un valore nullo e pertanto gli effetti dissipativi sono sensibilmente inferiori. E’ da osservare inoltre che le oscillazioni, dopo una prima riduzione, sembrano mantenere un’ampiezza costante nel tempo. Tale comportamento è da attribuirsi alla natura intrinsecamente instabile del metodo per cui la riduzione di ampiezza prodotta dalle (deboli) dissipazioni di energia è bilanciata dalla tendenza all’amplificazione delle oscillazioni a causa dell’instabilità del metodo.

Le precedenti valutazioni sono state effettuate riducendo il passo di integrazione temporale a Δt=0.2 s. Per i due casi (manovra di apertura e manovra di chiusura) i confronti sono illustrati in Figura 2.10 e in Figura 2.11, rispettivamente. Per entrambe le situazioni è stata considerata la dissipazione localizzata alla base del pozzo piezometrico. Modeste sono le differenze riscontrabili per il caso di manovra di apertura, mentre nel caso di chiusura si osserva che l’ampiezza delle oscillazioni si riduce sensibilmente prima che, per il meccanismo numerico prima evidenziato, si stabilisca una condizione di moto periodico non smorzato.


Figura 2.10 – Manovra di apertura: confronto dei livelli nel pozzo piezometrico calcolati con due

diversi passi di integrazione temporale.

Figura 2.11 – Manovra di chiusura: confronto dei livelli nel pozzo piezometrico calcolati con due

diversi passi di integrazione temporale.

Soluzione analitica con dissipazioni linearizzate. Nel seguito si riporta una semplice soluzione analitica del problema ottenuta linearizzando opportunamente le dissipazioni di energia lungo la galleria e trascurando quelle alla base del pozzo.

Immaginiamo di linearizzare le dissipazioni di energia introducendo una velocità caratteristica vM e ponendo

gvv

df

gvv

dfj M

22≅= (61)

In tal caso l’equazione dinamica (48), nella quale si trascura il termine ΔE2-3 che indica il complesso di perdite di energia che si realizzano nel pozzo, si scrive


gvv

dfL

dtzd

AA

gLjL

dtzd

AA

gLz M

g

p

g

p

22

2

2

2

−−=−−= (62)

Potendo esprimere la velocità in galleria ancora mediante la (50), si avrà

⎟⎟⎠

⎞⎜⎜⎝

⎛+−−=

g

f

g

pM

g

p

AQ

dtdz

AA

gv

dfL

dtzd

AA

gLz

22

2

(63)

ovvero

p

fM

p

gM

AQ

dvfz

LAgA

dtdz

dvf

dtzd

222

2

−=++ (64)

Posto

p

g

ALAg

=ω d

vf M

4=ψ

p

fM

AQ

dvfE 22 ω

=Δ (65)

l’equazioni del moto diventa

222

2

2 ωωψ Ezdtdz

dtzd

Δ−=++ (66)

Consideriamo il polinomio caratteristico dell’omogenea associata

02 22 =++ ωψλλ (67)

le cui soluzioni sono

222,1 ωψψλ −±−= (68)

La soluzione dell’omogenea associata sarà pertanto

tt eBeAz 21 λλ += (69)

A questa sommiamo un integrale particolare della (66) ottenendo così la soluzione completa che sarà

EeBeAz tt Δ−+= 21 λλ (70)

nella quale, in generale, i coefficienti A e B che dovranno essere determinati imponendo le condizioni al contorno, come pure gli esponenti λ1 e λ2, sono numeri complessi.

Consideriamo dapprima il caso in cui sia ψ>ω e le radici λ1 e λ2 siano pertanto reali.

All’istante t=0, il livello nel pozzo, z0, è quello relativo al moto stazionario precedente l’istante in cui viene effettuata una manovra. Dalla (70) si ha dunque

EBAz Δ−+=0 (71)


Ancora per t=0 vale, ovviamente la condizione al contorno (55) discendente dall’equazione di continuità. Essendo

tt eBeAdtdz

2121

λλ λλ += (72)

si avrà

( )g

f

g

p

AQ

BAAA

v 0210 ++= λλ (73)

Risolvendo il sistema lineare costituito dalle due equazioni (71) e (73), si trova

⎟⎟⎠

⎞⎜⎜⎝

⎛Δ+−

−−

= )(102

00

21

EzA

QQA

p

f λλλ

(74)

⎟⎟⎠

⎞⎜⎜⎝

⎛Δ+−

−−

−= )(101

00

21

EzA

QQB

p

f λλλ

(75)

Utilizzando l’equazione (50) possiamo inoltre determinare l’andamento nel tempo delle velocità in galleria che risulta espresso dalla seguente relazione

( )g

ftt

g

p

AQ

eBeAAA

v ++= 2121

λλ λλ (76)

Consideriamo ora il caso in cui sia ψ<ω e le radici λ1 e λ2 siano pertanto complesse. Posto

22 ψωω −=D (77)

si avrà

Dii ωψψωψλ ±−=−±−= 22 (78)

e la soluzione completa, comprensiva dell’integrale particolare, risulta

( ) EeBeAez titit DD Δ−+= −− ωωψ (79)

La precedente relazione, nella quale i coefficienti A e B sono complessi (in particolare sono complessi coniugati) può essere scritta nella seguente forma1 1 Ricordando che è

)()cos( ααα senie i += si ha:

[ ] EtsenitBtsenitAez DDDDt Δ−++−+−= − ))()(cos())()(cos( ωωωωψ

[ ] EtsenitBtsenitAez DDDDt Δ−++−= − ))()(cos())()(cos( ωωωωψ

[ ] EtsenBAitBAez DDt Δ−−−+= − )()()cos()( ωωψ

Posto: C1=-i(A-B) e C2=A+B si trova, infine [ ] EtCtsenCez DD

t Δ−+= − )cos()( 21 ωωψ


[ ] EtCtsenCez DDt Δ−+= − )cos()( 21 ωωψ (80)

All’istante t=0, il livello nel pozzo, z0, è quello relativo al moto stazionario precedente l’istante in cui viene effettuata una manovra. Dalla (80) si ha dunque

ECz Δ−= 20 (81)

Ancora per t=0 vale, ovviamente la condizione al contorno (55) discendente dall’equazione di continuità. Essendo

[ ] [ ])cos()()()cos( 1212 tCtsenCetsenCtCedtdz

DDt

DDDt ωωωωωψ ψψ −−+−= −− (82)

si avrà

( )210

00 CCAdtdzAQQ Dp

tpf ψω −==−

=

(83)

Risolvendo il sistema lineare costituito dalle due equazioni (81) e (83), si trova

DpD

f EzAQQC

ωψ

ω)( 000

1Δ+

+−

= (84)

EzC Δ+= 02 (85)

e quindi

EtsenEz

AQQ

tEzez DDpD

fD

t Δ−⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛ Δ++

−+Δ+= − )(

)()cos()( 000

0 ωω

ψω

ωψ (86)

Conviene riscrivere la precedente equazione (86) come segue

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛+

Δ+−

+Δ++Δ−= − )()(

)cos()(0

000 tsen

EzAQQ

teEzEz DDpD

fD

t ωωψ

ωωψ (87)

nella quale è evidente che per t→0, z→z0 (condizione di moto stazionario iniziale) e per t→∞, z→-ΔE, corrispondente alla nuova condizione di moto stazionario.

Utilizzando l’equazione (50) possiamo inoltre determinare l’andamento nel tempo delle velocità in galleria che, con qualche passaggio, risulta espresso dalla seguente relazione

⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛++

−

Δ+−

−+= − )()1(

)()cos( 2

2

00

0000 tsenQQ

AEzte

AQQ

AQ

v DDDf

pDD

t

g

f

g

f ωωψ

ωψω

ωψ (88)

nella quale è ancora evidente che per t→0, v→Q0/Ag (condizione di moto stazionario iniziale) e per t→∞, v→Qf0/Ag, corrispondente alla nuova condizione di moto stazionario.


Per quanto riguarda la stima della velocità caratteristica vM non è possibile stabilire un criterio che consenta di individuare per questo parametro un valore compatibile al tempo stesso con le condizioni di regime iniziale e finale e con le dissipazioni che si determinano durante il processo di moto vario.

Per considerare correttamente il primo aspetto, per manovre di completa chiusura o completa apertura conviene assumere per vM la velocità che si stabilisce nella galleria di derivazione rispettivamente all’istante iniziale (vM=Q0/Ag) ovvero a regime (vM=Qf0/Ag).

Per ottenere una stima verosimile della rapidità con cui le oscillazioni si smorzano nel tempo conviene assegnare a vM il valore massimo della velocità in galleria dedotto, ad esempio, dalla soluzione in assenza di dissipazioni fornita dall’equazione (57).

ESEMPI

Manovre multiple. Nel seguito si illustra brevemente il comportamento del sistema quando ad una prima manovra (effettuata a partire da una condizione di moto stazionario) ne segue una seconda a breve distanza temporale. Per semplicità, nella discussione farà riferimento al caso di dissipazioni di energia trascurabili.

A partire da una condizione stazionaria, in un istante generico t0, il livello nel pozzo piezometrico e la velocità in galleria sono quelli forniti dalle equazioni (56) e (57), rispettivamente, ovvero

)()( 000

0 tsenAQQtz

p

f ωω

−= (89)

g

f

g

f

AQt

AQQtv 0

000

0 )cos()( +−

= ω (90)

Consideriamo ora la soluzione generale del problema in assenza di dissipazioni, prima dell’imposizione delle condizioni al contorno, fornita dalla (53).

Imponiamo le nuove condizioni al contorno: per t=t0, z=z(t0) e, sempre per t=t0, Qf=Qf1.

Dalle (53) e (89) si avrà

)()cos()()( 000

000 tsenAQQtBtsenAtz

p

f ωω

ωω−

=+= (91)

Mentre dalle (50) e (90) si avrà

[ ]

g

f

g

f

g

f

g

p

g

f

ttg

p

AQ

tA

QQAQtsenBtA

AA

AQ

dtdz

AA

tv

00

00

100

10

)cos(

)()cos()(0

+−

=

=+−=+==

ω

ωωω (92)

Dalla soluzione del sistema lineare composto dalle equazioni (91) e (92) si determinano i coefficienti A e B. Sostituite, quindi, le espressioni per questi coefficienti nella soluzione generale (53) si ottiene


)cos()()()cos( 001

00100 ttsen

AQQ

tsentAQQ

AQQ

zp

ff

p

ff

p

f ωωω

ωωωω ⎥

⎥⎦

⎤

⎢⎢⎣

⎡ −+

⎥⎥⎦

⎤

⎢⎢⎣

⎡ −−

−= (93)

valida, evidentemente, per t>t0. Consideriamo, ad esempio, una manovra iniziale di chiusura (Qf0=0) e una

successiva manovra di apertura effettuata quando il livello nel pozzo ha raggiunto il suo valore minimo. Dalla (56), questo istante corrisponde a t=t0=0.75T. Sostituite queste condizione nella (93), si avrà

)cos()( 10 tA

QtsenA

Qz

p

f

p

ωω

ωω ⎥

⎥⎦

⎤

⎢⎢⎣

⎡−+

⎥⎥⎦

⎤

⎢⎢⎣

⎡= (94)

Immaginando che a seguito dell’apertura, la portata verso la condotta forzata, Qf1, corrisponda alla portata verso la condotta forzata che si aveva antecedentemente la prima manovra di chiusura, cioè la portata Q0 fluente in galleria, si avrà

[ ])cos()(0 ttsenA

Qz

p

ωωω

−= (95)

la quale può essere riscritta nella forma

)4/(2 0 πωω

−= tsenA

Qzp

(96)

La precedente relazione mostra che, a partire dall’istante t=t0 in cui si effettua la seconda manovra, l’oscillazione resta caratterizzata da un’ampiezza incrementata del fattore 2 .

La situazione più gravosa si verifica quando la manovra di chiusura avviene all’istante t0 corrispondente al momento in cui il livello nel pozzo, nella sua discesa, raggiunge il valore z=0 e la portata in galleria, che sta andando verso il serbatoio, raggiunge il suo valore massimo. Questo istante corrisponde a t=t0=0.5T. Sostituite queste condizione nella (93), si avrà

)(10 tsenA

QA

Qz

p

f

p

ωωω ⎥

⎥⎦

⎤

⎢⎢⎣

⎡+= (97)

Assunto, come in precedenza, Qf1=Q0, si avrà

)(2 0 tsenA

Qz

p

ωω

= (98)

la quale mostra che, a partire dall’istante t=t0 in cui si effettua la seconda manovra, l’oscillazione resta caratterizzata da un’ampiezza doppia rispetto a quella precedente l’istante t0.

ESEMPI


2.3.3 Casse d’aria Un dispositivo di protezione analogo a quello visto nel paragrafo precedente è

necessario inserire in un impianto di sollevamento immediatamente a valle delle pompe. Un pozzo piezometrico potrebbe assolvere il compito di protezione della condotta di mandata dai fenomeni di colpo d’ariete innescati dall’attacco e dallo stacco delle pompe, ma per la sua collocazione (immediatamente a valle delle pompe) dovrebbe essere caratterizzato da dimensioni verticali spesso eccessive (vedi Figura 2.12), superiori al dislivello geodetico tra serbatoio di recapito e vasca di alimentazione. Per tale motivo il pozzo piezometrico è normalmente sostituito da un piccolo serbatoio in pressione, detto cassa d’aria, che sfrutta la comprimibilità dell’aria per assorbire le rapide variazioni dinamiche imposte dall’attacco o dall’arresto della pompa.

Consideriamo il sistema illustrato in Figura 2.12 che schematizza un impianto di sollevamento. Il sistema è costituito da un serbatoio di alimentazione A, da una breve condotta di aspirazione, da una pompa, da una cassa d’aria posta immediatamente a valle della pompa e da una lunga condotta di mandata fino al serbatoio di recapito B.

Figura 2.12

Per effetto della regolazione della pompa si innesca un fenomeno di oscillazione di massa che interessa la cassa d’aria, la condotta di mandata e il serbatoio di recapito. Il problema è del tutto analogo a quello dell’oscillazione di massa in un pozzo piezometrico illustrato nel paragrafo precedente.

Le equazioni che governano l’oscillazione di massa nella cassa d’aria sono ancora quelle per il moto vario unidimensionale in ipotesi anelastiche. Effettuando l’integrazione spaziale dell’equazione del moto tra le sezioni 1 e 2 indicate in Figura 2.12 e trascurando i termini di inerzia e quelli dissipativi nella cassa e nel serbatoio, si ha

jLdtdv

gLEE −−=− 12 (99)

nella quale L, v e j rappresentano la lunghezza, la velocità e la dissipazione di energia per unità di lunghezza relativi alla condotta di mandata.


Assunto il riferimento indicato in Figura 2.12 e trascurato il carico cinetico nella cassa si ha E1=z+p/γ, essendo p la pressione dell’aria nella cassa. Si può inoltre assumere che l’estensione del serbatoio di recapito sia notevole e siano pertanto trascurabili le variazioni di quota della sua superficie. Si può assumere cioè E2=h0. In queste ipotesi, l’equazione (99) può essere riscritta come segue

jLdtdv

gLpzh −−=+− )(0 γ

(100)

All’equazione (100), nelle variabili z, v e p, va innanzitutto associata l’equazione di continuità che, nel caso in esame, è un equazione di nodo:

ccp AvAvQ += (101)

nella quale A e Ac sono le aree delle sezioni della condotta di mandata e della cassa d’aria rispettivamente, vc è la velocità nella cassa d’aria e Qp e la portata pompata. Potendo scrivere vc=dz/dt, la precedente relazione diventa:

cp AdtdzAvQ += (102)

A queste due relazioni ne va associata necessariamente una terza che contenga almeno una delle variabili del problema. L’equazione che possiamo associare è quella che descrive la trasformazione termodinamica subita dall’aria contenuta nella cassa. Tale trasformazione sarà in generale una politropica del tipo:

kk UpUp *0

*0 = (103)

Nella quale p* e U sono la pressione assoluta e il volume dell’aria nella cassa ad un generico istante, p0

* e U0 sono le stesse quantità all’istante iniziale e k è l’esponente della politropica che per l’aria varia tra k=1 (trasformazione isoterma) e k=1.4 (trasformazione adiabatica. Indicata con zM la quota del cielo della cassa, la precedente relazione può essere scritta come segue

kMcatm

kMcatm zzAppzzApp )()()()( 00 −+=−+ (104)

ovvero k

Matmk

Matm zzppzzpp )()()()( 00 −+=−+ (105)

in cui patm è la pressione atmosferica (patm=101.33 kPa) e z0=0 è il livello iniziale dell’acqua nella cassa.

Assumendo noto l’andamento nel tempo della portata Qp, il sistema composto dalle equazioni (100), (102) e (105) può essere risolto per fornire, in particolare, l’andamento nel tempo dei livelli e delle pressioni nella cassa d’aria. A tale proposito, si vedrà più avanti un esempio di soluzione numerica del problema.

Per ottenere una semplice soluzione analitica di questo problema assumiamo innanzitutto di trascurare le dissipazioni di energia. Dalla (102) si ha


dtdz

AA

AQ

v cp −= (106)

Sostituita la (106) nella (100), quest’ultima diventa

dtdQ

AgL

dtzd

AA

gLpzh pc −=+− 2

2

0 )(γ

(107)

Ancora per semplicità ammettiamo che la manovra alla pompa che innesca l’oscillazione di massa sia una variazione istantanea della portata Qp che si realizzi all’istante t=0, dalle condizioni precedenti Qp=v0A (v0 è la velocità nella condotta di mandata all’istante t=0) a quelle valide per t>0, Qp=Qp0. Con questa semplificazione, per t>0 risulta dQp/dt=0 e l’equazione (107) può essere ulteriormente semplificata.

Una seconda, fastidiosa, non linearità è presente nell’equazione (105). A tale proposito va osservato che le oscillazioni di livello nella cassa d’aria sono relativamente contenute essendo legate alla comprimibilità dell’aria. Possiamo quindi pensare di linearizzare la suddetta equazione mediante uno sviluppo in serie di Taylor in un intorno di z=z0=0. L’espressione per la pressione p, esplicitata dall’equazione (105) si scrive:

atmkM

kMatm p

zzzzppp −

−−+=

)(1)()( 00 (108)

Lo sviluppo in serie di Taylor, arrestato al termine del primo ordine si scrive:

zzppp

zz

00

== ∂

∂+≅ (109)

Sostituita l’espressione per p fornita dall’equazione (108) nella (109), con qualche passaggio si ottiene:

zz

ppkpp

M

atm )( 00

++≅ (110)

Sostituita infine quest’ultima relazione nell’equazione (107), si ha:

2

200

0 ))((dt

zdAA

gLz

zppkp

zh c

M

atm =+

++−γγ

(111)

La precedente è un’equazione differenziale lineare del secondo ordine, a coefficienti costanti, che può essere facilmente risolta.

Si osservi innanzitutto che è z0+p0/γ=p0/γ=h0. La relazione precedente è riscritta come segue:

0))(

1( 02

2

=⎥⎦

⎤⎢⎣

⎡ +++ z

AA

Lg

zppk

dtzd

cM

atm

γ (112)

Posto

cM

atm

AA

Lg

zppk

))(

1( 0

γω

++= (113)


l’equazione differenziale (112), analoga a quella vista per i pozzi piezometrici, presenta la soluzione generale

)cos()( 21 tCtsenCz ωω += (114)

in cui le costanti C1 e C2 vanno determinate utilizzando le condizioni al contorno. All’istante t=0, il livello nella cassa d’aria è z=0. Dalla (114) si ha dunque

221 )0cos()0(0 CCsenC =+= (115)

Ancora per t=0 vale, ovviamente, l’equazione di continuità (106).

0

00

=

−=t

cp

dtdz

AA

AQ

v (116)

Derivando la (114) e sostituendo tale derivata nella (116) si trova C1=(Qp0-Q0)/ωAc, in cui Q0=v0

.A è la portata fluente lungo la condotta di mandata all’istante t=0. La soluzione, pertanto è

)(00 tsenA

QQz

c

p ωω

−= (117)

Il coefficiente C1 risulta quindi essere l’ampiezza delle oscillazioni di livello nella cassa d’aria. Sostituendo questa espressione per z(t) nell’equazione di continuità è possibile determinare l’andamento nel tempo della velocità o della portata lungo la condotta di mandata, mentre dall’equazione (108) si ottiene l’andamento nel tempo delle pressioni. In alternativa, la pressione potrebbe essere calcolata, istante per istante, sostituendo l’espressione per z(t) nell’equazione linearizzata (110). Nell’esempio che segue sono analizzate entrambe queste possibilità.

Soluzione numerica del problema di oscillazione in una cassa d’aria. Nel seguito si illustra una semplice soluzione numerica del problema utilizzando un metodo di integrazione di tipo esplicito. E’ da osservare, come per i casi precedenti, che la soluzione con un metodo esplicito, in assenza di dissipazioni di energia, per questo sistema di equazioni, risulta essere instabile. In presenza di dissipazioni di energia la stabilità può essere assicurata scegliendo un passo di integrazione temporale sufficientemente piccolo.

Le equazioni (100), (102) e (105) sono riscritte come segue

⎥⎦

⎤⎢⎣

⎡−−+= jLhpz

Lg

dtdv

0γ (118)

c

p

AAvQ

dtdz −

= (119)

[ ] atmkM

kMatm p

zzzzppp −

−−+=

)(1)()( 00 (120)



nktkt ,,3,2,1 K=Δ=

La generica variabile all’istante t, è indicata con l’apice k. Ad esempio z(t) è indicato con zk. Discretizzando le due derivate che compaiono nelle equazioni (118) e (119) alle differenze finite, e con riferimento ad un metodo di integrazione esplicito,si trova:

kkk

jLhpzLg

tvv

⎥⎦

⎤⎢⎣

⎡−−+=

Δ−+

0

1

γ

c

kkp

kk

AAvQ

tzz −

=Δ

−+1

(121)


kk jLhpzL

tgvv ⎥⎦

⎤⎢⎣

⎡−−+

Δ+=+

01

γ

tA

AvQzz

c

kkpkk Δ

−+=+1

(122)

Ad ogni passo temporale, poi, mediante l’equazione (120) si calcola il corrispondente valore per la pressione p.

Il sistema è completato dalle condizioni iniziali e al contorno per il particolare problema che si vuole studiare.

Consideriamo ad esempio (Figura 2.12) una condotta con le seguenti caratteristiche: L=1000 m e d=0.2 m. Assumiamo, come riferimento per le quote il livello nella cassa d’aria ad impianto fermo come indicato in Figura 2.12. Ipotizziamo inoltre, per semplicità, che le dissipazioni di energia possano essere calcolate mediante la formula di Darcy-Weisbach nella quale la funzione di resistenza assuma un valore costante e pari a f=0.02. Le caratteristiche della cassa d’aria siano Ac=0.5 m2 e zM=1.0 m. Sia, infine, h0=40.0 m.

Nel caso di attacco di una pompa, le condizioni iniziali sono rappresentate da: z=0, v=0 e p=γh0=392.4 kPa. Le condizioni al contorno sono rappresentate dall’andamento nel tempo delle portate pompate. Si può, ad esempio assumere: Qp=0 per t=0 e Qp=0.05 m3/s per t>0.

La soluzione, ottenuta utilizzando un passo di integrazione temporale Δt=0.02 s è illustrata in Figura 2.13. Nella stessa figura, per confronto, è riportata la soluzione analitica ricavata in precedenza.


Figura 2.13 – Andamento nel tempo del livello z e dell’altezza di pressione p/γ nella cassa d’aria e

della portata Q lungo la condotta di mandata nel caso di attacco della pompa. La soluzione 1 usa l’equazione (110) per la pressione mentre la soluzione 2 usa l’equazione (108).


Nella successiva Figura 2.14 è illustrata la soluzione dello stesso problema

nell’ipotesi di ridurre praticamente a zero le dissipazioni di energia (si assume f=0.0002). In questo caso la soluzione numerica si avvicina sensibilmente a quella analitica. Le differenze che permangono sono da attribuire alla linearizzazione dell’equazione termodinamica necessaria per poter ottenere una soluzione analitica del problema.

Si osservi, in Figura 2.13, che in condizioni di regime si ha: z=0.149 m, p=517.59k Pa (p/γ=52.76 m) e, ovviamente, Q=0.05 m3/s.

A partire da queste condizioni di moto stazionario è stato simulato lo stacco della pompa. In questo caso le condizioni al contorno, rappresentate dall’andamento nel tempo delle portate pompate, sono descritte dalle seguenti relazioni: Qp=0.05 m3/s per t=0 e Qp=0 per t>0.

I risultati della simulazione numerica, ottenuti utilizzando un passo di integrazione temporale Δt=0.02s sono illustrati in Figura 2.15. Anche in questo caso, nella stessa figura è riportata per confronto la soluzione analitica precedentemente ricavata.


della portata Q lungo la condotta di mandata e confronto con la soluzione teorica. I risultati sono stati ottenuti riducendo il coefficiente di resistenza al valore f=0.0002. I simboli sono gli stessi della precedente Figura 2.13.



della portata Q lungo la condotta di mandata nel caso di stacco della pompa. La soluzione 1 usa l’equazione (110) per la pressione mentre la soluzione 2 usa l’equazione (108).


2.3.4 Oscillazioni forzate

Ci sono situazioni in cui l’oscillazione è non indotta da una qualche manovra, più o meno rapida, di un organo di regolazione come negli esempi visti nei paragrafi precedenti, ma è forzata da una variazione continua, spesso periodica o quasi periodica, di una condizione al contorno.

Anche in situazioni di questo tipo la soluzione numerica del problema non lineare è di gran lunga più efficace quando si debbano fare precise valutazioni quantitative. Per inquadrare il problema e studiarne le caratteristiche più rilevanti è comunque utile procedere, in prima approssimazione, attraverso la soluzione analitica delle equazioni che governano il fenomeno, in qualche modo linea rizzate.

Consideriamo l’esempio illustrato in Figura 2.16 nel quale il sistema è costituito da una condotta che collega un pozzetto ad un bacino. Assumiamo che nel bacino il livello oscilli (periodicamente) nel tempo seguendo la legge h=a.sen(ψt) che descrive un’oscillazione attorno al livello medio nel bacino, assunto come riferimento, di ampiezza a e pulsazione ψ e la risposta del sistema a questa forzante.

Figura 2.16

Le equazioni che governano l’oscillazione forzata sono ancora quelle per il moto vario unidimensionale in ipotesi anelastiche. Possiamo integrare nello spazio l’equazione del moto tra le sezioni 1 e 2 indicate in Figura 2.16 e, mantenendo i contributi dei termini dissipativi e delle accelerazioni temporali relativi alla sola condotta, si trova

∫∫ −∂∂

−=−'2

'1

'2

'112

1 dsjdstv

gEE (123)

Assunto il riferimento indicato in Figura 2.16 e trascurato il carico cinetico nel pozzetto si ha E2=z e E1=h=a.sen(ψt). Pertanto l’equazione (47) può essere riscritta come segue

jLdtdv

gLtsenaz −−=− )(ψ (124)

in cui v è la velocità nella condotta. All’equazione (48), nelle variabili z e v, va associata l’equazione di continuità che, come per il caso del pozzo piezometrico, è un equazione di nodo (nodo 2’ di Figura 2.16):


Ω= pvAv (125)

nella quale A e Ω sono le aree delle sezioni della condotta e del pozzetto rispettivamente, vp è la velocità nel pozzetto che si può esprimere in funzione del livello z essendo vp=dz/dt. Il sistema composto dalle equazioni (48) e (49), con le necessarie condizioni al contorno, può essere agevolmente risolto numericamente. Soluzioni analitiche si possono ottenere introducendo alcune ulteriori semplificazioni.

Consideriamo dapprima il caso in cui si trascurino le dissipazioni di energia. Combinando le equazioni (48) e (49) e con la usuale posizione

Ω=

LgA2ω

si trova

tsenazdt

zd ψωω 222

2

=+ (126)

La soluzione generale dell’omogenea associata all’equazione completa (126) è quella già vista nei precedenti problemi di oscillazione di massa. A questa soluzione va aggiunto l’integrale particolare zp

1

)(22

2

tsenazp ψψω

ω−

= (127)

per ottenere, nel complesso

)()(cos)( 22

2

21 tsenatCtsenCz ψψω

ωωω−

++= (128)

Le costanti di integrazione C1 e C2 si determinano utilizzando le opportune condizioni al contorno che dipendono dal particolare problema indagato.

Esempio. Consideriamo il sistema illustrato in Figura 2.16. La condotta è lunga 6 m e ha un diametro di 0.2 m cui corrisponde un’area A=0.0314 m2. Il pozzetto è cilindrico di diametro D=0.6 m cui corrisponde una superficie orizzontale Ω=0.283 m2. Assumiamo inoltre che l’oscillazione del livello nel bacino sia caratterizzata da un’ampiezza a=0.5 m e da un periodo TW=10 s (che potrebbe corrispondere ad un modo ondoso). Assumiamo, inoltre, che fino all’istante t=0 la condotta e il pozzetto siano isolati dal bacino (possiamo immaginare, ad esempio, che vi sia una saracinesca chiusa in corrispondenza della sezione 1’ di Figura 2.16) e il livello nel 1 L’integrale particolare si determina a partire da una sua verosimile struttura algebrica. Nel caso in esame si può assumere zp=B1 sen(yt)+B2 cos(yt) in cui le costanti B1 e B2 vanno determinate in modo tale che l’espressione per zp verifichi l’equazione differenziale completa. Sostituita dunque quest’espressione per zp nell’equazione (126) si ottiene il cosiddetto residuo R il quale deve risultare identicamente nullo. Nel caso in esame si trova

[ ] [ ] )(cos)( 22

22

221

21 tBBtsenaBBR ψωψψωωψ +−+−+−=

Il residuo è identicamente nullo quando sono nulli i coefficienti delle funzioni sen(ψt) e cos(ψt). Imponendo questa condizione si trova

22

2

1 ψωω−

=aB 02 =B


pozzetto coincida con quello medio nel bacino. All’istante t=0 la saracinesca viene rapidamente (istantaneamente) aperta realizzando il collegamento dinamico tra bacino e pozzetto attraverso la condotta.

Nelle ipotesi qui assunte, per la valutazione delle costanti di integrazione C1 e C2

si impongono le seguenti condizioni al contorno

⎪⎩

⎪⎨⎧

=Ω

=

==

0

00 0 Av

dtdzz

t (129)

la seconda delle quali è l’equazione di continuità (125) esplicitata rispetto a dz/dt e calcolata per t=0.

Con queste condizioni, si trova

221 ψωωψ−

−=aC 02 =C

L’oscillazione del livello nel pozzetto sarà quindi descritta dalla seguente relazione

⎥⎦⎤

⎢⎣⎡ −

−= )()(22

2

tsentsenaz ωωψψ

ψωω (130)

nella quale, con i dati del problema, si ha ψ=2π/TW=0.628 s-1 e ω=0.426 s-1. L’oscillazione del livello nel pozzo, descritta analiticamente dall’equazione (130), è illustrata in Figura 2.17.

Si osservi che con le caratteristiche geometriche ipotizzate e la particolare forzante assunta, le oscillazioni di livello nel pozzetto sono caratterizzate da un’ampiezza superiore (praticamente doppia) a quella della forzante. Ciò è dovuto al fatto che differenza tra la frequenza della forzante ψ /2π e quella propria del sistema ω /2π è relativamente piccola e siamo dunque vicini alla condizione di risonanza.

Dall’equazione di continuità (125) si determina l’andamento delle velocità in condotta che sono quindi espresse dalla seguente relazione (vedi Figura 2.17)

[ ])(cos)(cos22

2

ttaA

v ωψψωψω

−−

Ω= (131)

Consideriamo ora il problema in presenza di dissipazioni di energia. In tal caso il

sistema composto dalle equazioni (124) e (125), unitamente alle condizioni al contorno, può essere agevolmente risolto per via numerica, come si vedrà più avanti. Per ottenere una soluzione analitica è necessario linearizzare il termine che descrive le dissipazioni di energia assumendo che queste siano proporzionali alla velocità in condotta.

Il modo di procedere è del tutto analogo a quello sviluppato nel paragrafo 2.3.2. In particolare si introduce una velocità caratteristica vM e si pone

gvv

df

gvv

dfj M

22≅=


Figura 2.17

In questo modo l’equazione dinamica (124) si scrive

dtdz

gv

AdfL

dtzd

gL

Atsenaz M

2)( 2

2 Ω−

Ω−=− ψ

ovvero

)(22

2

tsenLgAaz

LgA

dtdz

dvf

dtzd M ψ

Ω=

Ω++

Posto

Ω=

LAgω

dvf M

4=χ (132)

l’equazioni del moto diventa

)(2 222

2

tsenazdtdz

dtzd ψωωχ =++ (133)


La soluzione dell’omogenea associata (zo) è già stata illustrata in dettaglio nel paragrafo 2.3.2. In particolare, se χ>ω risulta

tto eCeCz 21

21λλ += con 22

2,1 ωχχλ −±−= (134)

Nel caso in cui sia χ<ω, posto

22 χωω −=D (135)

si avrà

[ ])cos()( 21 tCtsenCez DDt

o ωωχ += − (136)

Le costanti C1 e C2 che compaiono nelle equazioni (134) e (136) andranno determinate utilizzando le condizioni al contorno. Alla soluzione generale rappresentata da queste equazioni si deve poi sommare l’integrale particolare.

Per la valutazione dell’integrale particolare si procede come già descritto nella soluzione in assenza di dissipazioni. Si trova

[ ])cos(2)()()()2(

222222

2

ttsenazp ψχψψψωψωχψ

ω−−

−+= (137)

La soluzione completa sarà dunque z=zo+zp.

Esempio. Consideriamo l’esempio appena illustrato ed assumiamo per le dissipazioni di energia continue un valore costante per la funzione di resistenza f=0.02.

Assumiamo, per la velocità caratteristica un valore verosimile, vM=2.0 m/s (si veda

a questo proposito il grafico di Figura 2.17. Con i dati del problema risulta χ=0.05 s-1. Essendo χ<ω, la soluzione del problema sarà del tipo

[ ] pDDt ztCtsenCez ++= − )cos()( 21 ωωχ (138)

in cui zp è espresso dalla (137) e ωD=0.423 s-1, poco inferiore a ω. Per determinare le costanti C1 e C2 si utilizzano le stesse condizioni al contorno

viste in precedenza ed espresse dalla (129). Dopo qualche (laborioso) passaggio si trova

2222

2222

1 )()2(2

ψωχψψωχ

ωψω

−++−

=D

aC

2222

2

2 )()2(2

ψωχψωχψ

−+= aC

(139)

Con i dati del problema si trova C1=0.596 m e C2=0.116 m. L’oscillazione del livello nel pozzo, descritta analiticamente dall’equazione (130), è

illustrata in Figura 2.18. Dall’equazione di continuità (125) si determina quindi l’andamento delle velocità in

condotta che sono espresse dalla seguente relazione (vedi Figura 2.18)


[ ][ ])(2)(cos)(

)()()/()cos()()()2(

22

22222222

2

tsent

tsenteaA

v DDDt

ψχψψψω

ωψωωχωψωψωχψ

ψω χ

+−

+++−−−+

Ω= −

(140)

Figura 2.18

Con riferimento a questa soluzione vi sono alcuni aspetti che è opportuno sottolineare.

Nella soluzione generale che si ottiene sommando l’integrale dell’omogenea associata (136) all’integrale particolare (137) le costanti di integrazione C1 e C2 che dipendono dalle condizioni iniziali del sistema (condizioni al contorno nel tempo) interessano soltanto la parte di soluzione espressa dalla (136). Questa parte di soluzione, inoltre, decade esponenzialmente nel tempo e il sistema quindi tende verso una condizione di regime che è descritta dal solo integrale particolare. Questo comportamento è evidente in Figura 2.18 in cui è riportato anche l’andamento nel tempo del contributo zp.


La rapidità con cui il contributo fornito dall’omogenea associata va riducendosi è descritta dal termine M=1/χ che rappresenta, pertanto, una misura della memoria del sistema. In prima approssimazione, e in analogia con quanto visto per il tempo di avviamento di una condotta, potremmo assumere che la memoria effettiva del sistema sia circa 3M, ovvero dopo un tempo pari a 3M il sistema non ricorda più le condizioni iniziali di partenza.

Nell’esempio appena illustrato si ha M=20 s (quindi 3M=60 s) e in effetti, come è immediato osservare in Figura 2.18, dopo circa 60 s, la soluzione completa z(t) e quella descritta dall’integrale particolare zp(t) risultano praticamente coincidenti.

Questo comportamento, caratteristico delle “oscillazioni forzate” viene frequente sfruttato nella soluzione di problemi per i quali è difficoltoso fissare le condizioni iniziali del sistema. In questi casi, stimata la memoria dl sistema, si estende opportunamente all’indietro il periodo di tempo simulato in modo che la soluzione nella finestra temporale che interessa non sia affetta dalle condizioni iniziali, poste arbitrariamente.

In molti problemi, inoltre, il transitorio descritto dalla soluzione dell’omogenea associata risulta breve e poco significativo mentre l’analisi del comportamento del sistema a regime riveste spesso un certo interesse. Queste condizioni si determinano soprattutto quando il sistema è fortemente dissipativo, cioè le dissipazioni di energia prevalgono rispetto agli effetti legati all’inerzia temporale.

In condizioni di regime si ha

[ ])cos(2)()()()2(

222222

2

ttsenazz p ψχψψψωψωχψ

ω−−

−+==

che può essere riscritta come segue

)()()2( 2222

2

ϕψψωχψ

ω+

−+= tsenaz e ⎟⎟

⎠

⎞⎜⎜⎝

⎛−

= 22

2arctanωψ

χψϕ (141)

Se confrontiamo l’andamento del livello nel pozzetto, espresso dalla (141) con l’andamento nel tempo del livello nel bacino si osserva una variazione di ampiezza e un ritardo di fase che dipendono dalle caratteristiche del sistema e dalla pulsazione ψ della forzante. E’ spesso interessante analizzare come varia il rapporto η tra le ampiezze nel pozzetto e nel bacino, e il ritardo di fase ϕ, al variare della pulsazione ω. Nei problemi dell’idraulica, inoltre, risulta di maggiore immediatezza considerare il ritardo temporale Δt anziché quello di fase ϕ, le due grandezze essendo legate dalla relazione Δτ=-ϕ/ψ.

Nel problema in esame si ha dunque

2222

2

)()2( ψωχψ

ωη−+

= e ⎟⎟⎠

⎞⎜⎜⎝

⎛−

−=Δ 22

2arctan1ωψ

χψψ

t (142)

Gli andamenti di questi due parametri al variare della pulsazione ψ (e del rapporto tra la frequenza della forzante e quella propria del sistema, ψ/ω) sono illustrati in Figura 2.19. Con riferimento al rapporto η tra l’ampiezza della forzante e quella delle oscillazioni nel pozzetto, si osserva un’attenuazione (η<1) solo per le frequanze maggiori (ψ>0.6 s-1) mentre per quelle minori l’ampiezza della risposta z(t) si amplifica con un valore massimo di η superiore a 4 in condizioni prossime a quelle di risonanza (ψ/ω≈1). Per quanto riguarda il ritardo di fase si osserva che anche per


valori molto piccoli della pulsazione ψ, ovvero per oscillazioni del livello nel bacino molto lente, il ritardo che caratterizza le oscillazioni nel pozzetto, pur piccolo, non si annulla, a causa dell’inerzia temporale.

Figura 2.19

Si noti, che qualsiasi sia l’andamento della grandezza assegnata, si può pensare di descrivere questo andamento mediante uno sviluppo in serie di Fourier. Essendo il problema lineare posso studiare la soluzione relativa ad una singola e generica armonica e poi sovrapporre gli effetti.

Soluzione numerica del problema di oscillazione forzata. Nel seguito si illustra una semplice soluzione numerica del problema utilizzando un metodo di integrazione di tipo esplicito. E’ da osservare, come per i casi precedenti, che la soluzione con un metodo esplicito, in assenza di dissipazioni di energia, risulta essere instabile mentre in presenza di dissipazioni di energia la stabilità può essere assicurata scegliendo un passo di integrazione temporale sufficientemente piccolo.

Procediamo come già visto negli altri esempi presentati in questo capitolo. Dividiamo l’asse dei tempi in intervalli Δt tali per cui

nktkt ,,3,2,1 K=Δ=

ed indichiamo con l’apice k la generica variabile all’istante t. Discretizzando le due derivate temporali che compaiono nelle equazioni (124) e

(125) alle differenze finite, e con riferimento ad un metodo di integrazione esplicito,si trova:


kkk

tsenazgvv

dLf

Lg

tvv

⎥⎥⎦

⎤

⎢⎢⎣

⎡−+−=

Δ−+

)(2

1

ψ

kkk

vAt

zzΩ

=Δ−+1


kk tsenazgvv

dLft

Lgvv

⎥⎥⎦

⎤

⎢⎢⎣

⎡−+Δ−=+ )(

21 ψ

kkk vAtzzΩ

Δ+=+1

(143)

Il sistema è completato dalle condizioni iniziali e al contorno per il particolare problema che si vuole studiare.

Figura 2.20


Con riferimento ai dati dell’esempio illustrato in questo paragrafo, la Figura 2.20 mostra il confronto tra la soluzione numerica, ottenuta con un passo di integrazione temporale Δt=0.05 s e quella analitica delle equazioni linearizzate assumendo una velocità caratteristica vM=2 m/s. Si osserva un ottimo accordo tra le due soluzioni ed è facile verificare che tale accordo è poco influenzato dal valore assunto dalla velocità vM a patto che questa si mantenga in un intervallo di valori verosimili.


2.4 Colpo d’ariete Il fenomeno del colpo d’ariete sarà descritto con maggior dettaglio nel capitolo 1

relativo ai problemi che richiedono una trattazione in ipotesi elastiche. Si tratta di un fenomeno determinato da brusche variazioni delle condizioni di moto in grado di produrre importanti stati di sovrapressione o forti depressioni. Il fenomeno, come sarà possibile vedere analizzandolo in ipotesi elastiche, è aggravato dal fatto che queste variazioni positive e negative di pressione si alternano con frequenze molto elevate.

In questo capitolo viene analizzato il fenomeno del colpo d’ariete nella forma semplificata consentita dalle ipotesi anelastiche ricavando una soluzione che può essere ritenuta accettabile solo quando le variazioni del campo di moto sono relativamente lente e, conseguentemente, modeste le variazioni di pressione che così si determinano.

Consideriamo il sistema schematicamente illustrato in Figura 2.21 costituito da un serbatoio e da una lunga condotta di scarico munita, al termine, di un ugello regolabile (otturatore). Assumiamo inoltre, per comodità, che il moto nel sistema sia stazionario. Note le caratteristiche della condotta e il dislivello h0. è facile stabilire la velocità v0 di regime nella condotta. Immaginiamo di chiudere istantaneamente l’ugello terminale e di determinare come evolve il moto nel tempo. Si tratta, qualitativamente, di un processo opposto a quello dell’avviamento di una condotta.

Figura 2.21

Integrando l’equazione dinamica tra i punti 1 e 2, e trascurando i termini di accelerazione temporale e di dissipazione di energia nel tratto 1-1’, come fatto nel paragrafo 2.1, si trova:

Ljdtdv

gLEE −−=− 12 (144)

in cui L la lunghezza della condotta. Per comodità assumiamo che durante il processo il livello nel serbatoio rimanga costante e, per semplicità, trascuriamo le dissipazioni di energia ammettendo, come spesso accade, che la velocità in condotta e la lunghezza L siano piccole.

In queste ipotesi, l’equazione (144) si riduce alla seguente


dtdv

gLh

gvh −=−+ 0

2*

2 (145)

in cui h* è la quota piezometrica in corrispondenza alla sezione 2 (il pedice 2 è stato omesso per comodità). Una ulteriore semplificazione, adottata come si vedrà anche nella trattazione elastica del problema, consiste nel trascurare il termine cinetico a fronte di quello piezometrico. L’equazione (145), nella sua forma semplificata finale, è dunque

dtdv

gLhh −=− 0

* (146)

Nella precedente equazione le grandezze variabili nel tempo sono due: la quota piezometrica h* e la velocità in condotta v. E’ necessario quindi associare un’ulteriore relazione per rendere risolubile il problema. 2.4.1 Riduzione lineare della velocità

Se immaginiamo una manovra di chiusura dell’otturatore, potremmo assumere, in prima approssimazione, che questa chiusura determini una riduzione della velocità nella sezione 2 (e quindi in tutta la condotta) dal valore v0 a zero nel tempo di manovra Tc. Potremmo altresì assumere che questa variazione di velocità sia lineare nel tempo e segua quindi la legge:

)/1(0 cTtvv −= (147)

In tal caso, sostituita la (147) nella (146) si trova

cTgvL

hh 100

* =− (148)

ovvero

⎟⎟⎠

⎞⎜⎜⎝

⎛+=

c

a

TT

hh 10* (149)

La soluzione appena illustrata prevede che durante la chiusura dell’otturatore, in corrispondenza alla sezione terminale si determini una sovrapressione costante nel tempo il cui valore è fornito dalla (149). A chiusura avvenuta, il termine dv/dt si annulla e si ristabilisce immediatamente la condizione h*=h0. E’ da osservare che al tendere a zero della durata della manovra la sovrapressione che si viene a creare tende ad infinito. Prima che tale condizione sia raggiunta, evidentemente intervengono nel fenomeno gli effetti legati alla deformabilità della condotta e alla comprimibilità del fluido. Per tempi di manovra brevi, quindi, le ipotesi anelastiche sono inaccettabili. Nella pratica, la soluzione fornita dalla (149) può ritenersi accettabile fino a valori Ta/Tc=0.10÷0.15.

Nell’ipotesi di riduzione lineare delle velocità, il problema presenta un semplice soluzione anche quando siano mantenuti il termine dissipativo e l’energia cinetica nella sezione 2. In questo caso, infatti, l’equazione (144), associata alla (147) si scrive


⎟⎠⎞

⎜⎝⎛ +

−−=−

dLf

gTtv

Tv

gLhh c

c

12

)/1( 2200

0* per t<Tc

VELOCITA’ SENZA SEMPLIFICARE LE DISSIPAZIONI

2.4.2 Chiusura lineare dell’otturatore In alternativa all’ipotesi precedente, si può assumere una riduzione lineare

dell’area della sezione 3 dal valore iniziale AuMAX (otturatore completamente aperto) al valore zero. In tal caso la trattazione risulta un po’ più articolata.

.Trascurando gli effetti inerziali e le dissipazioni di energia nel tratto 2-3, l’equazione dinamica integrata nello spazio (vedi eq. (144)) si riduce all’uguaglianza E2=E3. Quest’ultima, trascurando, come già visto, il carico cinetico nella sezione 2 si scrive:

gv

h2

23* = (150)

L’equazione di continuità, applicata tra le sezioni 2 e 3 fornisce

uAvAv ⋅=⋅ 3 (151)

in cui A è l’area della sezione della condotta e Au è l’area della sezione terminale dell’ugello (sezione 3) in un generico istante. Combinando le equazioni (150) e (151) per eliminare la variabile v3, si trova:

*2 hgAA

v u= (152)

Inoltre, nelle condizioni di regime iniziali, quando l’otturatore è completamente aperto (Au=AuMAX), la velocità in condotta è v0, e la quota piezometrica in corrispondenza alla sezione 2 è h*=h0, dalla (152) si ha

00 2 hgA

Av uMAX= (153)

Dal rapporto tra le equazioni (152) e (153), si trova

0

*

0 hh

AA

vv

uMAX

u= (154)

Indicati con η=Au/AuMAX il grado di chiusura dell’ otturatore e con 0* / hh=ζ , la

(154) si scrive

)()(0 ttvv ζη ⋅⋅= (155)

L’espressione (155) per la velocità può essere sostituita nella (146) che diventa


⎟⎠⎞

⎜⎝⎛ ⋅+⋅−=−

dtd

dtd

gvL

hh ζηζη00

* (156)

Ovvero, dividendo entrambi i membri per h0 e ricordando l’espressione del tempo di avviamento:

⎟⎠⎞

⎜⎝⎛ ⋅+⋅−=−

dtd

dtdTa

ζηζηζ 12 (157)

Una volta fissata la legge di chiusura η(t), l’equazione (157) risulta contenere la sola variabile dipendente ζ. Se assumiamo, ad esempio, una legge di chiusura lineare del tipo

cTt /1−=η (158)

L’equazione (157) diventa

⎟⎟⎠

⎞⎜⎜⎝

⎛−−=−−

ca

c

a

Tt

dtdT

TT

112 ζζζ (159)

Si tratta di un’equazione differenziale non lineare e a variabili separabili. La (159), può infatti essere scritta come segue

( ) 1)/(/1)/(

2 −−−=

− ζζζ

cac

a

c

c

TTd

TT

TtTtd

(160)

Posto λ=Ta/Tc, la soluzione analitica della precedente, con la condizione ζ=1 per t=0 è

( ) ( )

( )( )

λλ

λλ

λλ

λλλλζ

/4

2

2

22

2

1)(4242

)(4)(4

21

+

⎟⎟⎠

⎞⎜⎜⎝

⎛−=

++−

+−−=

−⋅++−⋅+−

=

cTttBA

AtBAtB

(161)

Un esempio dell’andamento nel tempo della soluzione è illustrato in Figura 2.22.


Figura 2.22

D’altra parte noi siamo interessati soprattutto alla stima dei valori massimi della sovrapressione, che si raggiungono, praticamente, in tempi brevi e sensibilmente inferiori al tempo di manovra (a rigore, il valore massimo lo si ha per t=Tc). Si raggiungono, infatti, in un tempo confrontabile con quello necessario all’avviamento di una condotta già visto nel paragrafo 2.1.

Come si osserva dalla soluzione illustrata in Figura 2.22, in corrispondenza ai massimi valori della sovrapressione la variazione nel tempo della quota h* e quindi della variabile ζ è trascurabilmente piccola. Posto dζ/dt=0, possiamo ridurre l’equazione (159) ad una semplice equazione algebrica per i valori estremi di ζ.

012 =−− MAXc

aMAX T

Tζζ (162)

L’equazione (162) presenta, come unica soluzione fisicamente accettabile, la seguente

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡+⎟⎟

⎠

⎞⎜⎜⎝

⎛+= 4

21

2

c

a

c

aMAX T

TTT

ζ (163)

La precedente relazione può essere semplificata ricordando che le ipotesi anelastiche, come visto nel paragrafo 2.4.1, possono essere ritenute accettabili solo quando Ta/Tc«1. In tal caso (Ta/Tc)2 è trascurabile rispetto a 4 e la (163) si riduce a

c

aMAX T

T211+=ζ (164)

Ricordando la definizione di ζ, si può scrivere 22

2

0

*

411

211 ⎟⎟

⎠

⎞⎜⎜⎝

⎛++=⎟⎟

⎠

⎞⎜⎜⎝

⎛+==

c

a

c

a

c

aMAX

MAX

TT

TT

TT

hh

ζ (165)


Ricorrendo alla stessa semplificazione utilizzata per la (163), si trova infine

⎟⎟⎠

⎞⎜⎜⎝

⎛+=

c

aMAX T

Thh 10

* (166)

Si ottiene, in pratica, lo stesso risultato conseguito con la trattazione vista nel paragrafo 2.4.1.

Moto vario in ipotesi elastiche

- 55 -

3 MOTO VARIO IN IPOTESI ELASTICHE Nel paragrafo precedente si è visto che la trattazione anelastica presenta seri limiti

quando le brusche variazioni delle caratteristiche del moto determinano forti e rapide variazioni di pressione. In queste situazioni l’effetto della pressione consiste in non trascurabili deformazioni della condotta e variazioni di massa legate alla variazione della densità del fluido.

Per queste situazioni si rende necessario ricorrere alla trattazione elastica. Almeno nella sua forma semplificata che sarà illustrata e descritta in questo capitolo.

3.1 Descrizione qualitativa del fenomeno del colpo d’ariete Prima di introdurre le equazioni del moto, si ritiene opportuno analizzare

qualitativamente il fenomeno del colpo d’ariete come si sviluppa nella realtà. A tale scopo consideriamo il sistema schematico illustrato in Figura 3.1a. La situazione è del tutto analoga a quella schematizzata in Figura 2.21 e analizzata sotto le ipotesi anelastiche. Immaginiamo di chiudere istantaneamente l’otturatore all’istante t=0. In corrispondenza della sezione terminale, quindi, la velocità deve annullarsi e la pressione cresce essendo la forza di pressione l’unica forza che si oppone all’inerzia del fluido che sta fluendo verso l’otturatore. L’incremento di pressione determina un aumento della densità del fluido (legge di Hook) e la dilatazione della condotta. Dopo pochi istanti, tutta la parte terminale della condotta è caratterizzata da velocità nulla e da pressioni superiori a quella iniziale (Figura 3.1b). Questo stato di moto si propaga verso monte con celerità a. Mentre nel tratto di monte della condotta, non ancora raggiunto dalla perturbazione di pressione, permangono le condizioni iniziali. Una volta raggiunto il serbatoio, l’onda viene riflessa e si propaga verso l’otturatore con la stessa celerità. Il passaggio dell’onda riflessa determina il ripristino delle pressioni ma produce un flusso caratterizzato dalla velocità v0 diretta verso il serbatoio (Figura 3.1c). All’istante t=2L/a l’onda di pressione riflessa ha raggiunto l’otturatore. Le pressioni sono tornate ovunque quelle iniziali ma adesso la condotta è animata da un flusso opposto a quello iniziale (Figura 3.1d). Trovando chiuso l’otturatore, in questa sezione si produce una nuova onda di pressione, questa volta negativa, uguale ed opposta a quella iniziale a causa del verso della velocità. Qualche istante più tardi (Figura 3.1e) si osserva l’onda che viaggia verso monte con celerità a e lascia alle sue spalle uno stato di depressione e velocità nulle. Quest’onda raggiunge il serbatoio e viene da questo riflessa. . Il passaggio dell’onda riflessa determina il ripristino delle pressioni e modifica ancora una volta il verso del flusso che ora è caratterizzato dalla velocità v0 diretta verso l’otturatore (Figura 3.1f). All’istante t=4L/a l’onda di pressione riflessa ha raggiunto l’otturatore e sia le pressioni che la velocità tornano ad essere quelle dell’istante iniziale (Figura 3.1a). Il fenomeno prosegue periodicamente esaurendosi lentamente per effetto delle dissipazioni di energia. Gli andamenti nel tempo delle sovrapressioni all’otturatore e delle velocità in condotta in corrispondenza dell’imbocco del serbatoio sono illustrati in Figura 3.2.


Figura 3.1


Figura 3.2 – Andamento nel tempo delle sovrapressioni all’otturatore (sinistra) e delle velocità

all’imbocco del serbatoio (destra)

Il tempo t=2L/a, necessario alla perturbazione di pressione prodotta all’otturatore per ritornare all’otturatore stesso dopo essere stata riflessa dal serbatoio è un tempo caratteristico che riveste, come vedremo, notevole importanza. Tale tempo caratteristico è detto “ritmo” della condotta ed è normalmente indicato con τ.

Nella realtà le manovre non sono istantanee ma caratterizzate da una durata finita Tc. Per valutare, ancora qualitativamente il moto che si sviluppa a seguito di una manovra non istantanea conviene rappresentare graficamente in modo diverso la propagazione delle onde di pressione. Essendo il fenomeno lineare, vale il principio di sovrapponibilità degli effetti. Possiamo quindi studiare il campo di moto complessivo sommando le onde che a passi temporali pari a τ/2 partono dall’otturatore o dal serbatoio. Per familiarizzare con questa rappresentazione, consideriamo ancora il caso di chiusura istantanea.

In Figura 3.3, nella colonna di sinistra, sono rappresentate le diverse onde mentre nella colonna di destra la loro somma, quindi lo stato di pressione effettivo. All’istante t=0, dall’otturatore parte un’onda positiva, cioè di sovrapressione (0<t<τ/2) che arrivata al serbatoio cambia segno (τ/2<t<τ). L’onda riflessa al serbatoio arriva all’otturatore all’istante t=τ e qui viene nuovamente riflessa mantenendo però lo stesso segno negativo (τ<t<1.5τ). Giunta al serbatoio, l’onda viene riflessa con segno opposto e torna quindi ad essere un’onda di sovrapressione (1.5τ<t<2τ). Quest’onda raggiunge l’otturatore all’istante t=2τ. Qui viene nuovamente riflessa come onda positiva, la situazione complessiva è identica a quella che caratterizza l’istante t=0 e il processo si ripete periodicamente. In pratica una variazione di pressione Δp (positiva o negativa) che dall’otturatore arriva al serbatoio viene da questo riflessa con segno opposto. Viceversa, una variazione di pressione Δp che dal serbatoio arriva all’otturatore viene da questo riflessa con lo stesso segno.


Figura 3.3

Nel caso di una manovra di durata finita, possiamo immaginare di scomporre la riduzione graduale di area all’otturatore in una successione (al limite infinita) di piccole riduzioni istantanee della sezione terminale di efflusso (al limite infinitesime) effettuate a passi temporali δt piccoli ma finiti. Ogni riduzione di sezione determina un’onda di sovrapressione, di intensità relativamente modesta (δp) essendo modesta la riduzione di velocità indotta dalla parziale chiusura dell’otturatore. Le onde di sovrapressione così prodotte, sovrapponendosi, si propagano verso monte come illustrato schematicamente in Figura 3.4.


Figura 3.4

Al termine della manovra di chiusura, se questa è sufficientemente veloce, all’otturatore si avrà una sovrapressione Δp=nδp essendo n=Tc/δt. Tale sovrapressione si manterrà inalterata nel tempo finchè l’onda negativa riflessa al serbatoio non avrà raggiunto l’otturatore. Al tendere a zero dell’intervallo δt, quindi per una manovra graduale, il fronte d’onda assumerà un profilo continuo anziché scalinato.

E’ possibile analizzare in questo modo, qualitativamente, lo stato di pressione che interessa la condotta a seguito di una manovra non istantanea.

In questa analisi è necessario distinguere il caso in cui sia Tc<τ da quello in cui Tc risulti superiore al ritmo della condotta. Nel primo caso, illustrato in Figura 3.5, la manovra è definita rapida. Si osserva infatti che, essendo Tc<τ, almeno all’otturatore si raggiunge la sovrapressione massima, al termine della manovra di chiusura, coincidente con quella che si determina nel caso di chiusura istantanea. Ciò è determinato dal fatto che l’onda negativa riflessa dal serbatoio, e quindi il suo effetto consistente in una riduzione di pressione, non raggiunge l’otturatore prima della completa chiusura.


Figura 3.5 – Chiusura rapida

Nel secondo caso, illustrato in Figura 3.6, la manovra è definita lenta (Tc>τ). Si osserva infatti che, essendo Tc>τ, l’onda negativa riflessa dal serbatoio, e quindi il suo effetto consistente in una riduzione di pressione, raggiunge l’otturatore all’istate t=τ, prima della completa chiusura, cioè prima che all’otturatore si raggiunga la sovrapressione massima. In questo caso, quindi, le massime sovrapressioni sono inferiori a quelle corrispondenti ad una manovra rapida.


Figura 3.6 – Chiusura lenta

Nella successiva Figura 3.7 sono illustrati gli andamenti nel tempo della pressione

all’otturatore determinati da una chiusura caratterizzata da tempi di manovra Tc differenti. Si osserva, in particolare, che per Tc<τ la massima (e la minima) pressione assumono lo stesso valore anche se il periodo durante il quale la condotta è sollecitata da questa pressione massima si riduce al crescere del tempo di manovra. Quando invece sia Tc>τ, la massima sovrapressione all’otturatore, per le ragioni su esposte, va riducendosi al cresce della durata della manovra di chiusura.

Nel paragrafo che segue sono introdotte le equazioni normalmente utilizzate per la

soluzione del moto in pressione in ipotesi elastiche.


Figura 3.7


3.2 Le equazioni semplificate del moto Consideriamo le equazioni del moto (1)1 e (2) sviluppate come indicato di seguito

jtv

ggvh

s−

∂∂

−=+∂∂ 1)

2(

2* (167)

0=∂∂

+∂∂

+∂∂

+∂∂

tA

tA

sQ

sQ ρρρρ (168)

Sviluppando la derivata spaziale nella (167) e sostituendo nella (168) Q=v.A, si trova

jtv

gsv

gv

sh

−∂∂

−=∂∂

+∂∂ 1*

(169)

0=∂∂

+∂∂

+∂∂

+∂∂

+∂∂

tA

tA

svA

svA

sAv ρρρρρ (170)

Nelle equazioni così scritte introduciamo una serie di semplificazioni che si basano sulle seguenti considerazioni. Immaginiamo che una qualsiasi delle caratteristiche dell’onda di pressione sia descritta dalla generica funzione F(s,t). Tale funzione, ad esempio, può essere la pressione, o la densità. Se ci spostiamo solidali all’onda di pressione e assumiamo di trascurare le dissipazioni di energia, vediamo la caratteristica F immutata sia nel tempo che nello spazio. In altre parole, per un osservatore in moto con la celerità a dell’onda di pressione deve essere dF/dt=0. Sviluppando, si ha

0=∂∂

+∂∂

=dtds

sF

tF

dtdF (171)

in cui la velocità ds/dt è proprio la celerità a dell’onda. Esplicitando la precedente relazione rispetto ad a=ds/dt e dividendo entrambi i

membri per la velocità v in condotta si ha

sFvtF

v ∂∂⋅∂∂

−=/

/a (172)

1 In realtà l’equazione (1) contiene già un’ipotesi semplificativa determinata dalla trascurabilità della variazione spaziale della densità. L’equazione di partenza dovrebbe infatti essere:

jtv

ggv

ssh

sp

−∂∂

−=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

+∂∂

+∂∂ 1

21 2

γ

Solo se è possibile assumere (1/γ)∂p/∂s≈∂(p/γ)/∂s le due equazioni coincidono. D’Atra parte si ha

⎟⎟⎠

⎞⎜⎜⎝

⎛−

∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

−∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

−∂∂

=∂∂

−∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

dpdp

sp

sp

dpdp

sp

sp

sp

sp

spp

sρ

ργρ

ργρ

ργγ

γγγ11111

2

Utilizzando la legge di Hook richiamata più avanti (equazione (178)),si ha (p/ρ)dρ/dp=p/E, essendo E il modulo elastico del fluido. Le pressioni in gioco sono dell’ordine di qualche MPa mentre per l’acqua si ha E=2GPa. Il termine p/E risulta pertanto trascurabile rispetto all’unità e l’ipotesi semplificativa può considerarsi senz’altro accettabile.


Come avremo modo di vedere, la celerità a è normalmente molto più grande della velocità v. Ad esempio, nel caso di una condotta in acciaio, la celerità è dell’ordine dei 1000 m/s mentre la velocità è spesso dell’ordine 1÷3 m/s. Dalla (172) segue che i termini del tipo v.∂F/∂s sono sensibilmente più piccoli dei termini del tipo ∂F/∂t, per qualsiasi variabile F, e possono essere trascurati. Accade così, nella (169), che il termine v.∂v/∂s possa essere trascurato rispetto al termine ∂v/∂t , essendo entrambi moltiplicati per 1/g. Raccolta a fattor comune la densità ρ nel primo e nel quarto termine della (170), il primo risulta trascurabile. Analogamente, raccolta a fattor comune l’area A nel terzo e nel quinto termine della (170), il terzo risulta trascurabile.

Le equazioni (169) e (170) si riducono pertanto alle seguenti

jtv

gsh

−∂∂

−=∂∂ 1*

(173)

0=∂∂

+∂∂

+∂∂

tA

tA

svA ρρρ (174)

E’ da osservare che la semplificazione introdotta nell’equazione dinamica (173) corrisponde a trascurare il carico cinetico in condotta.

L’equazione (174) necessita di laboriose manipolazioni. Assumiamo innanzitutto che il moto sia barotropico e che quindi non solo l’area A della sezione ma anche la densità ρ del fluido possano modificarsi solo per effetto della pressione: A=A(p) e ρ=ρ(p). In queste ipotesi, dividendo entrambi i membri per il prodotto ρA, l’equazione (174) si riscrive come segue

011=

∂∂

+∂∂

+∂∂

tp

dpd

tp

dpdA

Asv ρ

ρ (175)

Essendo h*=p/γ+h, con h quota geodetica della condotta, costante nel tempo, si ha

tp

tp

th

∂∂

−∂∂

=∂

∂ γγγ 2

* 1 (176)

nella quale l’ultimo termine a destra risulta essere trascurabile1. L’equazione (175) può quindi essere riscritta nella forma

011 *

=∂

∂⎟⎟⎠

⎞⎜⎜⎝

⎛++

∂∂

th

dpd

dpdA

Asv ρ

ργ (177)

Si tratta a questo punto di formulare le espressioni per i due termini contenuti tra parentesi.

1 L’equazione (176), utilizzando la legge di Hook nella forma data dalla (178), può essere riscritta come

⎟⎠⎞

⎜⎝⎛ −

∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

−∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

−∂∂

=⎟⎟⎠

⎞⎜⎜⎝

⎛∂∂

−∂∂

=∂

∂Ep

tp

tp

dpdp

tp

tp

tp

tp

tp

th 11111*

γρ

ργρ

ργγ

γγ

in cui il modulo elastico del fluido è molto più grande dei valori usualmente ammissibili per la pressione. Ovvero, risulta p/E«1 e può essere trascurato.


Dalla legge di Hook: dρ/ρ=dp/E, con E modulo elastico del fluido, si ha direttamente

Edpd 11

=ρ

ρ (178)

Assumendo che la condotta sia di sezione circolare e diametro D, si ha

dpdD

DdpDd

DdpdA

A2)4/(

4/11 2

2 ==π

π (179)

Dalla formula di Mariotte: pD=2sσ, con s spessore della condotta e σ sforzo di trazione, si ha1

σdDsdp 2

= (180)

Utilizzando ancora la legge di Hook: dD/D=dσ/Em, con Em modulo di elasticità del materiale, l’equazione (180) può essere scritta come

dDDEsdp m2

2= (181)

Sostituita questa espressione nella (179) si ha infine

mEsD

dpdA

A=

1 (182)

L’equazione (177) può quindi essere scritta come segue

01 *

=∂

∂⎟⎟⎠

⎞⎜⎜⎝

⎛++

∂∂

th

EEsD

sv

m

γ (183)

O meglio, esplicitata rispetto alla variazione temporale della quota piezometrica

sv

EE

sD

Egt

h

m

∂∂

⎟⎟⎠

⎞⎜⎜⎝

⎛+

−=∂

∂

1

/1* ρ (184)

Dall’equazione (172), in cui si assume F=h*, si ha

shth

∂∂∂∂

−=//

*

*

a (185)

Sostituite nella precedente, le espressioni per ∂h*/∂s e ∂h*/∂t date dalle (173), nell’ipotesi di trascurare le dissipazioni di energia, e (184), rispettivamente, si trova 1 Differenziando la formula di Mariotte e utilizzando la legge di Hook nella forma dD/D=dσ/Em, con Em modulo di elasticità del materiale, si ha

⎟⎟⎠

⎞⎜⎜⎝

⎛−=−=−=

mm Ed

Ds

Ed

Dsd

DsdD

Dsd

Dsdp σσσσσσσ 122222

2

in cui il rapporto σ/Em è trascurabilmente piccolo rispetto all’unità.


tvsv

EE

sD

E

m

∂∂∂∂

⎟⎟⎠

⎞⎜⎜⎝

⎛+

−=//

1

/ ρa (186)

Utilizzando ancora l’equazione (172) nella quale si ponga F=v, l’equazione (186) diventa

⎟⎟⎠

⎞⎜⎜⎝

⎛+

=

mEE

sD

E

1

/2 ρa ovvero

⎟⎟⎠

⎞⎜⎜⎝

⎛+

±=

mEE

sD

E

1

/ ρa (187)

Quindi l’equazione (184) può infine essere scritta nella forma

sv

gth

∂∂

−=∂

∂ 2* a (188)

Sommando membro a membro le equazioni (173) e (188),la prima moltiplicata per a=ds/dt, si ottiene

aa jtv

dtds

sv

gdtds

sh

th

−⎟⎠⎞

⎜⎝⎛

∂∂

+∂∂

−=∂∂

+∂

∂ **

(189)

Ovvero

aa jdtdv

gdtdh

−−=*

(190)

L’equazione (190) è l’equazione che governa i fenomeni di moto vario elastico nelle condotte. Nell’equazione (190), come risulta dalla soluzione della (187), la celerità a può assumere valori positivi o negativi in relazione alla direzione di propagazione dell’onda di pressione. Pertanto, se si fa riferimento al modulo della celerità, la (190) rappresenta in realtà un sistema di due equazioni.

3.2.1 Forma integrata delle equazioni del moto

Con riferimento all’impiego della (190), conviene, per comodità, fare riferimento al modulo della celerità specificandone, con il segno, la direzione di propagazione. Nel seguito, con il simbolo a, come nel precedente paragrafo 3.1, si indicherà il modulo della celerità.

Per un’onda che si propaga nella direzione positiva con riferimento all’ascissa curvilinea s si avrà a=ds/dt. Tale relazione descrive la cosiddetta “linea caratteristica positiva” indicata normalmente con C+. Lungo tale linea vale la (190) che riscriviamo come segue, utilizzando per le dissipazioni continue la formula di Darcy-Weisbach

)(2

*+−−= Ca

gvv

df

dtdv

ga

dtdh (191)

Per un’onda che si propaga nella direzione negativa sarà a=-ds/dt. Tale relazione descrive la cosiddetta “linea caratteristica negativa” indicata normalmente con C-. Lungo tale linea vale ancora la (190) che riscriviamo come segue


)(2

*−+= Ca

gvv

df

dtdv

ga

dtdh (192)

Nell’ipotesi di poter considerare costante la celerità a, le due linee caratteristiche, in forma integrata, si scrivono

tsa

ΔΔ

= (193)

in cui Δs è lo spazio, in modulo, percorso dal fronte d’onda in un intervallo di tempo Δt.

Le equazioni (191) e (192), integrate lungo le caratteristiche, si scrivono

[ ] )(2

),(),(),(),(),(),( ** +Δ−−Δ+Δ+−=−Δ+Δ+ Cs

gtsvtsv

dftsvttssv

gatshttssh

(194)

[ ] )(2

),(),(),(),(),(),( ** −Δ+−Δ+Δ−=−Δ+Δ− Cs

gtsvtsv

dftsvttssv

gatshttssh (195)

Tali relazioni forniscono, con le opportune condizioni al contorno, i valori delle variabili h* e v a passi discreti spaziali Δs e temporali Δt.

3.2.2 Esempio di applicazione

Nel seguito si illustra un semplice esempio dell’applicazione della forma integrata delle equazioni per il moto vario elastico. Consideriamo il problema illustrato qualitativamente all’inizio di questo capitolo e risolto in quello precedente mediante l’impiego delle equazioni anelastiche.

Consideriamo un passo spaziale Δs coincidente con la lunghezza totale L della condotta e un passo di integrazione temporale che, conseguentemente vale Δt=L/a. Assumiamo di trascurare le dissipazioni di energia.

Per comodità, possiamo riscrivere le equazioni (194) e (195) nella seguente forma.

[ ] )(),0(),(),0(),( 121*

2* +−−=− CtvtLv

gathtLh (196)

[ ] )(),(),0(),(),0( 121*

2* −−=− CtLvtv

gatLhth (197)

in cui gli istanti t2 e t1 sono ovviamente legati dalla relazione t2=t1+Δt. Si ricorda che in corrispondenza della sezione s=0 (serbatoio) si ha costantemente

h*(0,t)=h0. Studiamo in particolare l’andamento delle pressioni all’otturatore (s=L) che

immaginiamo di chiudere all’istante t=0. Fino al momento della chiusura (t=t2<0, t1<0) la velocità lungo la condotta è quella di regime, cioè v(L,t2)=v(0,t1)=v0. Dalla (196) si ha


[ ] 0),( 0002* =−−=− vv

gahtLh

Dunque, fino all’istante di chiusura la quota piezometrica all’otturatore vale h0. Appena effettuata la chiusura (t=0+) si ha: v(L,0+)=0 e v(0, 0+-Δt)=v0. Dalla (196) si

trova

[ ]gva

vgahLh 0

00* 0)0,( =−−=−+

Cioè all’otturatore si registra una sovrapressione pari a (av0/g). A partire da questo momento v(L,t) resta sempre nullo. E’ da osservare inoltre che la velocità nella sezione s=0 rimane v0 fintantochè la perturbazione partita dall’otturatore non raggiunge tale sezione. Quindi per un intervallo di tempo pari a Δt. Durante questo intervallo di tempo, pertanto, all’otturatore si mantiene la medesima sovrapressione. All’istante t=Δt+ la perturbazione, dopo aver raggiunto il serbatoio è appena riflessa. Dalla (197) con t2=Δt+ e t1=0+, si ha

[ ]0),0()( 000 −Δ=+− +tv

ga

gva

hh

Ovvero v(0,Δt+)=-v0. Dalla (196), invece, si ha

[ ]gva

vgahtLh 0

00* 0),( =−−=−Δ +

Pertanto all’otturatore si mantiene ancora la precedente condizione di sovrapressione. Inoltre le condizioni appena viste si mantengono finchè l’onda riflessa non ha raggiunto l’otturatore, cioè fino all’istante t=2Δt. All’istante immediatamente successivo t=2Δt+, l’onda è già stata riflessa dall’otturatore. Dalla (196) con t2=2Δt+ e t1= Δt+, essendo v(0,Δt+)=-v0, si ha

[ ]gva

vgahtLh 0

00* )(0)2,( −=−−−=−Δ +

Cioè all’otturatore si registra una depressione pari a (-av0/g). La velocità nella sezione s=0 rimane -v0 fintantochè la perturbazione partita dall’otturatore non raggiunge tale sezione. Quindi per un intervallo di tempo pari a Δt. Durante questo intervallo di tempo, pertanto, all’otturatore si mantiene la medesima depressione. All’istante t=3Δt+ la perturbazione, dopo aver raggiunto il serbatoio è appena riflessa. Dalla (197) con t2=3Δt+ e t1=2Δt+, si ha

[ ]0)3,0()( 000 −Δ=−− +tv

ga

gva

hh

Ovvero v(0,3Δt+)=v0. Dalla (196), invece, si ha


[ ]gva

vgahtLh 0

00* )(0)3,( −=−−−=−Δ +

Pertanto all’otturatore si mantiene ancora la precedente condizione di depressione.


3.3 Modalità di impiego delle equazioni del moto Consideriamo una condotta di lunghezza L, suddivisa in un numero n di tratti di

lunghezza Δs. Il sistema composto dalle equazioni (194) e (195) può essere utilizzato per determinare le pressioni e le velocità nelle sezioni, estremi dei tronchi con cui è stata suddivisa la condotta, a passi temporali Δt finiti. Per quanto visto, i passi di integrazione spaziale Δs e temporale Δt restano legati tra di loro dalla relazione Δs=a·Δt e non possono pertanto essere scelti arbitrariamente. In genere viene fissata la discretizzazione spaziale e stabilito, di conseguenza il passo Δt.

Figura 3.8 – Illustrazione schematica del procedimento di integrazione nel tempo delle equazioni

Con riferimento alle indicazioni contenute in Figura 3.8, possiamo riscrivere le equazioni per le caratteristiche positiva e negativa, rispettivamente, come segue

[ ] sg

tsvtsvdftsvttsv

gatshttsh kk

kkkk Δ−−Δ+−=−Δ+ −−−− 2

),(),(),(),(),(),( 0101

01001*

0* (198)

[ ] sg

tsvtsvdftsvttsv

gatshttsh kk

kkkk Δ+−Δ+=−Δ+ ++++ 2

),(),(),(),(),(),( 0101

01001*

0* (199)

Sommando le due precedenti relazioni si ottiene

[ ] [ ]

[ ]),(),(),(),(4

),(),(2

),(),(21),(

01010101

010101*

01*

0*

tsvtsvtsvtsvdgsf

tsvtsvg

atshtshttsh

kkkk

kkkkk

++−−

+−+−

−Δ

−

−−++=Δ+ (200)

la quale consente di determinare la quota piezometrica nella posizione sk all’istante t0+Δt note che siano le caratteristiche del moto all’istante precedente, t0.

Analogamente, sottraendo tra loro le relazioni (198) e (199) ed esplicitando il risultato rispetto al termine v(sk, t0+Δt) si trova:

[ ] [ ]

[ ]),(),(),(),(4

),(),(2

),(),(21),(

01010101

01*

01*

01010

tsvtsvtsvtsvdasf

tshtsha

gtsvtsvttsv

kkkk

kkkkk

++−−

+−+−

+Δ

−

−−++=Δ+ (201)


Le equazioni (200) e (201) consentono, a partire da una condizione iniziale nota, di ricostruire i valori della quota piezometrica e della velocità in tutti i punti interni della condotta, cioè nelle posizioni sk con k=1,2,..n-1.

Per la soluzione completa del problema è necessario fissare le condizioni iniziali (condizioni al contorno nel tempo) e le condizioni al contorno nello spazio. 3.3.1 Condizioni al contorno

Le condizioni iniziali consistono nei valori di quota piezometrica e velocità per tutti i punti della discretizzazione all’istante iniziale t0. Più articolato è invece il problema dell’assegnazione delle condizioni al contorno nello spazio. E’ evidente, innanzitutto, per come sono fatte le linee caratteristiche, che è necessario fissare una condizione a monte ed una a valle (cioè non è possibile fissare entrambe le condizioni nella stessa sezione).

In entrambe le sezioni di estremità (s0 e sn) la condizione può essere una quota piezometrica assegnata (es.: h*(s0,t)=f(t), con f(t) nota), oppure una velocità assegnata (es.: v(s0,t)=f(t), con f(t) nota), oppure un legame (possibilmente lineare) tra quota piezometrica e velocità (es.: f[h*(s0,t), v(s0,t)]=0). Tale condizione al contorno, messa a sistema con l’unica equazione caratteristica disponibile per i punti di estremità, consente di ricavare la soluzione del problema in detti punti. Nel seguito si riportano alcuni semplici esempi di condizioni al contorno, in relazione al fenomeno che si vuole studiare. 3.3.2 Intercettazione del flusso in una condotta (colpo d’ariete)

Consideriamo l’esempio discusso nel paragrafo 2.4 (Figura 2.21), nel quale sono trascurate le dissipazioni di energia. Assumiamo, ad esempio, le seguenti caratteristiche del condotto: L=1000 m, d=0.5 m e una celerità di propagazione a=1000 m/s. Sia inoltre h0=10 m. In tal caso, trascurando le dissipazioni di energia, la velocità allo sbocco (sezione 3 di Figura 2.21) in condizioni di regime stazionario vale v3=(2gh0)1/2=14.0 m/s. Consideriamo un ugello caratterizzato (quando è completamente aperto) da una sezione Au0=0.028 m2. Per continuità, la velocità lungo la condotta in regime stazionario vale v0=2.0 m/s.

Se discretizziamo la condotta in 10 tratti della lunghezza di 100m ciascuno (Δs=100m) il passo temporale di integrazione resta fissato in Δt=Δs/a=0.1s.

La condizione iniziale è rappresentata da

0* )0,( hsh k = e 0)0,( vsv k = per k=1, 2, …., 11 (202)

La condizione al contorno di monte è

00* ),( htsh =

mentre per la condizione al contorno di valle, simulando una chiusura praticamente istantanea, si può assumere:

0),( =tsv n per t>0


Mettendo a sistema la condizione al contorno di monte con l’equazione per la caratteristica negativa, si ha:

[ ] sg

tsvtsvdftsvttsv

gatshh Δ+−Δ++=

2),(),(

),(),(),( 0101010001

*0

ovvero

[ ] [ ]),(),(2

),(),(),( 010101*

00100 tsvtsvdasftshh

agtsvttsv Δ

−−+=Δ+ (203)

Mettendo invece a sistema la condizione al contorno di valle con l’equazione per la caratteristica positiva, si ha:

[ ] [ ]),(),(2

),(0.0),(),( 01010101*

0* tsvtsv

dgsftsv

gatshttsh nnnnn −−−−

Δ−−−=Δ+ (204)

L’applicazione delle relazioni (200) e (201), con le condizioni iniziali indicate dalla (202) e le condizioni al contorno imposte, contenute nelle equazioni (203) e (204) consentono, posto per ipotesi f=0, di ricostruire l’andamento nel tempo (a passi finiti Δt) di velocità e quote piezometriche per ogni punto sk della condotta. Alcuni di questi andamenti sono riportati in Figura 3.9.

Volendo considerare le dissipazioni di energia, la condizione iniziale si modifica in quanto la quota piezometrica non è più costante ma si riduce dal serbatoio all’otturatore. Assumiamo, per semplicità, un valore costante per la funzione di resistenza nella formula di Darcy-Weisbach e pari a f=0.02.

Il moto stazionario di partenza si ottiene da un bilancio di energia nel quale si trascura il carico cinetico allo sbocco (ipotesi contenuta nelle equazioni elastiche). Risulta:

gv

dskfhsh k 2

)0,(20

0* Δ

−= e 0)0,( vsv k = per k=1, 2, …., 11 (205)

Per rendere più agevole il confronto con il caso di assenza di dissipazioni, si assume che la velocità lungo la condotta in condizioni di regime sia ancora v0=2.0 m/s. In queste condizioni, dalla (202), la quota piezometrica h*(L,0), immediatamente a monte dell’ugello, si riduce ad appena 1.8 m, la velocità in corrispondenza della sezione di sbocco vale v3=[2g h*(L,0)]1/2≈6.0 m/s, conseguentemente dobbiamo ammettere che l’area dell’ugello (quando è completamente aperto) sia caratterizzato da una sezione Au0=0.065 m2.


Figura 3.9 – Andamento nel tempo delle sovrapressioni all’otturatore e delle velocità all’imbocco dal

serbatoio in assenza di dissipazioni di energia.

I risultati che si ottengono applicando le relazioni viste in precedenza sono illustrati in Figura 3.10.

Aggiungere andamenti della densità, del diametro e delle tensioni Nei due esempi precedenti è stata considerata una manovra di chiusura

istantanea (in realtà il processo di discretizzazione numerica ha consentito lo studio di una manovra effettuata in un tempo molto piccolo ma finito: Tc=Δt=0.1s).

Consideriamo ora l’effetto della durata della manovra assumendo una riduzione lineare della velocità v0 in corrispondenza dell’otturatore. Va osservato che per il sistema in esame il ritmo della condotta vale τ=2L/a=2s. L’andamento nel tempo della quota piezometrica in corrispondenza dell’otturatore per tempi di chiusura variabili tra 0 e 8 s (4τ) è illustrato in Figura 3.11.


Figura 3.10 – Andamento nel tempo delle sovrapressioni all’otturatore e delle velocità all’imbocco del

serbatoio in presenza di dissipazioni di energia.


Figura 3.11 – Andamento nel tempo delle sovrapressioni all’otturatore per diversi tempi di manovra.

Con riferimento ai risultati illustrati in Figura 3.11, si possono fare alcune

considerazioni: • La massima sovrapressione all’otturatore (corrispondente alla massima quota

piezometrica) si mantiene sostanzialmente invariata per tempi di manovra Tc inferiori al ritmo τ della condotta.

• Per tempi di chiusura superiori al ritmo, la massima sovrapressione si riduce all’aumentare di Tc.


• In questo ultimo caso la massima sovrapressione non si presenta al termine della manovra ma ad istanti che sono multipli del ritmo. Si osservi, in particolare, che il primo picco si sovrapressione non corrisponde necessariamente al valore massimo.

• La massima depressione all’otturatore (corrispondente alla minima quota piezometrica) si mantiene sostanzialmente invariata per tempi di manovra Tc inferiori al ritmo τ della condotta.

• Per tempi di chiusura superiori al ritmo, la massima depressione risulta essere inferiore al caso di chiusura rapida ma tale depressione non si riduce monotonicamente all’aumentare di Tc.

L’andamento dei massimi valori di sovrapressione e di depressione in

corrispondenza all’otturatore al variare del tempo di chiusura Tc è illustrato in Figura 3.12 (sinistra) nella quale, per confronto, sono riportati i valori dedotti utilizzando l’equazione (149) ricavata in ipotesi anelastiche.

In Figura 3.12 (destra) sono confrontati invece i risultati che si ottengono considerando o trascurando le dissipazioni di energia. Si osservano differenze assolutamente modeste che consentono di concludere che, almeno per il caso analizzato, l’influenza delle dissipazioni nella stima dei valori estremi della sovrapressione, è assolutamente trascurabile (e in parte legata al diverso valore iniziale della quota piezometrica a monte dell’otturatore) a patto che la velocità di regime, nei due casi, si mantenga la stessa.

Figura 3.12 – Sovrapressioni e depressioni massime in corrispondenza dell’otturatore al variare del

tempo di chiusura e confronto con la soluzione anelastica (a sinistra) e con il caso di dissipazioni nulle (a destra).

Nel caso si proceda come visto nel paragrafo 3.2.2 in cui è stato considerato un

passo di integrazione spaziale coincidente con la lunghezza della condotta e, di conseguenza, un passo di integrazione temporale corrispondente a mezzo ritmo, è possibile combinare opportunamente le equazione per le caratteristiche positiva e negativa ed ottenere una successione di relazioni, dette equazioni concatenate di Allievi, che forniscono l’andamento nel tempo della soprapressione in corrispondenza


dell’otturatore nota che sia la velocità in questa sezione. In particolare, assumendo una variazione lineare della velocità da v0 a zero, si trova (formula di Allievi-Michaud)

c

a

TT

hkLh 21),(

0

*

+=τ per k=1, 3, 5,…. (206)

la quale, ovviamente, fornisce gli stessi valori di massima sovrapressione, al variare del tempo di chiusura, illustrati in Figura 3.12 (a destra).

Consideriamo ora l’effetto della durata della manovra assumendo una riduzione lineare dell’area dell’ugello. In tal caso, rispetto alle situazioni appena esaminate, va modificata la condizione al contorno di valle (quella di monte resta espressa dall’equazione (203)). Si procede in modo analogo a quanto visto per il colpo d’ariete in condizioni anelastiche. In particolare risulta ancora valida l’equazione (152) che viene riscritta per evidenziare che velocità e quota piezometrica si riferiscono alla sezione prossima all’otturatore

),(2)(

),( 0*0

0 ttshgA

ttAttsv n

un Δ+

Δ+=Δ+ (207)

La precedente, a sistema con l’equazione per la caratteristica positiva

[ ] [ ]),(),(2

),(),(),(),( 010101001*

0* tsvtsv

dgsftsvttsv

gatshttsh nnnnnn −−−−

Δ−−Δ+−=Δ+ (208)

consente di scrivere

CA

ttAa

AttA

attsv uun +⎥⎦

⎤⎢⎣

⎡ Δ++⎥⎦

⎤⎢⎣

⎡ Δ+−=Δ+

402

20

0)()(),( con

[ ]⎭⎬⎫

⎩⎨⎧ Δ

−+⎥⎦

⎤⎢⎣

⎡ Δ+= −−−− ),(),(

2),(),(

)(2 01010101

*2

0 tsvtsvdgsftsv

gatsh

AttA

gC nnnnu

(209)

Le condizioni iniziali sono descritte dall’equazione (202) in assenza di dissipazioni e dall’equazione (205) quando invece si vogliano considerare le dissipazioni di energia.

Anche in questo caso si assumono due diversi valori per l’area dell’ugello alla massima apertura (Au0=0.028 m2 in assenza di dissipazioni e Au0=0.065 m2 in presenza di dissipazioni) in modo che la velocità di regime in condotta sia la stessa: v0=2 m/s sia in presenza che in assenza di dissipazioni di energia. A causa dei diversi valori che assume quota piezometrica immediatamente a monte dell’otturatore per effetto delle perdite, però, la legge con cui la velocità nella sezione terminale della condotta si riduce nel tempo sarà diversa e tale differenza ha un impatto notevole sulla soluzione. Alcuni esempi di come si riduce la velocità nel tempo per diversi valori della durata Tc della manovra sono illustrati in Figura 3.13.

Si osserva che, soprattutto in presenza di dissipazioni di energia, la velocità all’otturatore si mantiene pressoché invariata per buona parte del tempo di manovra per poi ridursi rapidamente a zero. E’ evidente, quindi, che questa legge di chiusura risulta essere sensibilmente più severa di quella esaminata in precedenza. Si


osserva inoltre che la chiusura risulta essere più brusca in presenza di dissipazioni. C’è da aspettarsi, pertanto, che in queste condizioni le massime sovrapressioni siano maggiori rispetto al caso di assenza di dissipazioni.

Figura 3.13 – Leggi di riduzione delle velocità a monte dell’otturatore in assenza e in presenza di

dissipazioni di energia per diversi tempi di chiusura. I tempi di chiusura possono essere letti sull’asse temporale in corrispondenza dell’annullarsi della velocità.

L’andamento nel tempo della sovrapressione all’otturatore per alcuni tempi di manovra sono riportati in Figura 3.14. Come prevedibile, le sovrapressioni massime in questo caso risultano, a parità di tempo di manovra, superiori a quelle calcolate assumendo una riduzione lineare nel tempo della velocità all’otturatore. E’ confermato inoltre che le sovrapressioni massime in presenza di dissipazioni sono sensibilmente superiori a quelle che si ottengono trascurando le perdite.

Da questi pochi esempi emerge chiaramente che l’entità delle massime sovrapressioni dipende dal valore della velocità di regime v0 è dalle modalità con cui la velocità all’otturatore si riduce a zero. E’ necessario pertanto, nella soluzione di problemi reali, porre molta cura nel fissare queste condizioni iniziali e al contorno in quanto la loro influenza sul risultato è decisamente importante.


Figura 3.14 – Andamento nel tempo delle sovrapressioni all’otturatore per diversi tempi di manovra:

confronto tra la soluzione in presenza (linea continua) e assenza (linea tratteggiata) di dissipazioni di energia.

3.3.3 Avviamento del flusso in una condotta (tempo di avviamento)

Consideriamo l’esempio discusso nel paragrafo 2.1 (Figura 2.1) relativo al calcolo del tempo di avviamento di una condotta. Assumiamo, come per il precedente esempio, le seguenti caratteristiche del condotto: L=1000 m, d=0.5 m e una celerità di propagazione a=1000 m/s. Sia inoltre h0=10 m e f=0.025, costante. Il bilancio di energia, in condizioni di moto stazionario e nell’ipotesi di trascurare il carico cinetico allo sbocco, fornisce v0=1.98 m/s. Il tempo di avviamento vale pertanto: Ta=Lv0/gh0=20.2 s.

Se discretizziamo la condotta in 10 tratti della lunghezza di 100 m ciascuno (Δs=100m) il passo temporale di integrazione resta fissato in Δt=Δs/a=0.1 s.

La condizione iniziale è rappresentata da

0* )0,( hsh k = e 0)0,( =ksv per k=1, 2, …., 11 (210)


La condizione al contorno di monte è

00* ),( htsh =

mentre la condizione al contorno di valle potrebbe essere, simulando un’apertura praticamente istantanea:

0),(* =tsh n per t>0


[ ] sg

tsvtsvdftsvttsv

gatshh Δ+−Δ++=

2),(),(

),(),(),( 0101010001

*0

ovvero

[ ] [ ]),(),(2

),(),(),( 010101*

00100 tsvtsvdasftshh

agtsvttsv Δ

−−+=Δ+ (211)


[ ] [ ]),(),(2

),(),(),(0 010101001* tsvtsv

dgsftsvttsv

gatsh nnnnn −−−−

Δ−−Δ+−=

ovvero

[ ] [ ]),(),(2

),(),(),( 010101*

010 tsvtsvdasftsh

agtsvttsv nnnn −−−

Δ−+=Δ+ (212)

L’applicazione delle relazioni (200) e (201), con le condizioni iniziali indicate dalla (210) e le condizioni al contorno imposte, contenute nelle equazioni (211) e (212) consentono di ricostruire l’andamento nel tempo (a passi finiti Δt) di velocità e quote piezometriche per ogni punto sk della condotta. In particolare, in Figura 3.15 è riportato l’andamento delle velocità nelle sezioni iniziale e terminale della condotta e in una sezione intermedia. Per confronto, nella stessa figura, è illustrata la soluzione analitica fornita dall’equazione (17).


Figura 3.15 – Andamento nel tempo della velocità in tre sezioni lungo la condotta e confronto con la

soluzione anelastica.

3.3.4 Arresto di una pompa

Consideriamo il sistema illustrato in Figura 2.12 che schematizza un impianto di sollevamento e immaginiamo che non sia presente la cassa d’aria posta immediatamente a valle della pompa. Assumiamo, come nell’esempio illustrato nel paragrafo 2.3.3 che la condotta di mandata sia lunga L=1000 m e sia caratterizzata da un diametro d=0.2 m. Assumiamo inoltre, per semplicità, che le dissipazioni di energia possano essere calcolate mediante la formula di Darcy-Weisbach nella quale la funzione di resistenza sia costante e pari a f=0.02.

Come per gli esempi precedenti discretizziamo la condotta in 10 tratti della lunghezza di 100 m ciascuno (Δs=100m) e consideriamo una celerità dell’onda di pressione pari ad a=1000 m/s. Il passo temporale di integrazione resta pertanto fissato in Δt=Δs/a=0.1s.

Se consideriamo come condizione iniziale quella rappresentata dal moto stazionario che si instaura quando la pompa solleva la portata Qp=0.05m3/s, i valori per la quota piezometrica h* e la velocità v0 sono:

gv

dskLfhsh k 2

)()0,(20

0* Δ−

−= e 4/

)0,( 20 dQ

vsv pk π

== per k=1, 2, …., 11 (213)

La condizione al contorno di valle è

mhtsh n 40),( 0* ==


La condizione al contorno di monte, da porsi per la prima sezione della condotta, è, in un problema di questo tipo, molto delicata. Si è visto, con riferimento agli esempi riportati nel paragrafo 3.3.2, che la descrizione della manovra di chiusura ha un impatto considerevole sulla soluzione. Lo stesso, ovviamente, si verifica in un problema di attacco o di arresto di una pompa. In quest’ultimo caso l’inquadramento del fenomeno risulta ancora più complesso in quanto è necessario considerare la dinamica della pompa e in particolare il funzionamento della girante. Nel momento in cui si arresta una pompa (cioè si interrompe l’alimentazione elettrica), la girante, per inerzia, continua a ruotare e quindi a fornire energia al fluido. Rapidamente l’inerzia della girante si esaurisce ed entra in gioco l’inerzia del fluido che, pur rallentato, continua a fluire dalla vasca di carico alla condotta di mandata attraversando la pompa. In queste condizioni la girante, che ostacola questo moto, è assimilabile ad una dissipazione localizzata di energia. Quasi sempre, poi, nel momento in cui il flusso inverte il verso di percorrenza della condotta di mandata, entra in azione una valvola di non ritorno (clapet) posta immediatamente a valle della pompa. Da questa descrizione sintetica ed approssimata del funzionamento degli organi posti nella stazione di pompaggio emerge chiaramente la complessità della condizione al contorno di monte. Una descrizione, anche approssimata, del funzionamento meccanico reale di una pompa non rientra nel carattere di questi appunti

Nel seguito assumiamo per semplicità di porre a monte una condizione al contorno sulla portata pompata che dal valore Qp=0.05 m3/s passa linearmente a zero in un tempo di arresto Tc:

⎭⎬⎫

⎩⎨⎧

−= 0;1max),( 0c

p

Tt

AQ

tsv (214)


[ ] sg

tsvtsvdftsvttsv

gatshttsh Δ+−Δ++=Δ+

2),(),(

),(),(),(),( 0101010001

*00

* (215)


[ ] [ ]),(),(2

),(),(),( 010101001*

0 tsvtsvdgsftsvttsv

gatshh nnnnn −−−−

Δ−−Δ+−=

ovvero

[ ] [ ]),(),(2

),(),(),( 0101001*

010 tsvtsvdasfhtsh

agtsvttsv nnnn −−−

Δ−−+=Δ+ (216)

L’applicazione delle relazioni (200) e (201), con le condizioni iniziali indicate dalla (213) e le condizioni al contorno imposte, contenute nelle equazioni (214) e (215) consentono di ricostruire l’andamento nel tempo (a passi finiti Δt) di velocità e quote piezometriche per ogni punto sk della condotta. In particolare, in Figura 3.16 sono riportati gli andamenti delle quote piezometriche immediatamente a valle della pompa e in una sezione intermedia (s=500 m) e delle velocità nelle sezioni iniziale,


intermedia e terminale della condotta. Tali andamenti sono stati ottenuti assumendo un tempo di manovra Tc=Δt (cioè il minimo tempo di manovra che la soluzione numerica consente).

Figura 3.16 – Andamento nel tempo della quota piezometrica (sinistra) e delle velocità (destra) in

alcune sezioni lungo la condotta, per il caso di riduzione lineare della portata nel tempo Tc=0.1 s.

Nella successiva Figura 3.17 sono riportati gli stessi andamenti nel tempo calcolati assumendo, in luogo della (214), una condizione al contorno di monte espressa dalla seguente relazione:

cTtp eA

Qtsv /

0 ),( −= (217)

Una riduzione esponenziale della portata pompata rappresenta, in prima approssimazione, la soluzione che si ottiene considerando l’inerzia della girante. E’ da osservare che entrambe le condizioni (214) e (217) forniscono la stessa velocità di riduzione della portata in un intorno di t=0.


alcune sezioni lungo la condotta, per il caso di riduzione esponenziale della portata con tempo Tc=0.1s.

Le differenze tra i risultati che si ottengono con le due ipotesi di variazione nel tempo della portata pompata sono, in questo caso, modeste. E’ invece da osservare che la soluzione ottenuta è una soluzione “ideale”. Infatti, le quote piezometriche scendono a valori minimi sensibilmente inferiori allo zero. Essendo presumibilmente


la quota geodetica della pompa solo qualche metro sotto lo zero (che in Figura 2.12 coincide con il livello nella cassa d’aria), si avrebbero pressioni minime inferiori alla tensione di vapore. In queste condizioni quindi il problema reale si complica sensibilmente per la formazione di sacche di vapore saturo che determinano la rottura della colonna liquida e seri problemi quando, dopo pochi istanti, la pressione torna ad assumere valori positivi. Anche i valori positivi di pressione sono tali da destare almeno qualche preoccupazione con riferimento alle sollecitazione cui risulta soggetta la condotta.

La situazione migliora, ma non sufficientemente se si considera un tempo Tc=1.5τ=3 s. Nelle successive Figura 3.18 e Figura 3.19 sono riportati i risultati del calcolo per le due diverse condizioni al contorno di monte. Anche in questo caso la pressione a valle della pompa presenta valori minimi estremamente bassi anche se non drammaticamente bassi come nel caso precedente.


alcune sezioni lungo la condotta, per il caso di riduzione lineare della portata nel tempo Tc=3 s.


alcune sezioni lungo la condotta, per il caso di riduzione esponenziale della portata con tempo Tc=3 s.

3.3.5 Arresto di una pompa con cassa d’aria

Consideriamo lo stesso esempio appena analizzato e immaginiamo che sia presente la cassa d’aria posta immediatamente a valle della pompa come illustrato in


Figura 2.12. Rispetto alla soluzione appena vista si modifica la condizione al contorno sulle velocità relative alla prima sezione. In particolare l’equazione (214) (ovvero la (217)) vanno sostituite con un’equazione di continuità di nodo:

AQ

AtQ

tsv cassap −=)(

),( 0 (218)

in cui Qcassa rappresenta la portata entrante nella cassa d’aria. Quest’ultima dipende dalle variazioni di livello nella cassa e quindi dalla quota piezometrica nel nodo.

L’equazione (218) può essere manipolata e riscritta segue nella seguente forma:

)(),()(),( 20*

10 ttctshttctsv Δ−+Δ−= (219)

in cui c1(t) e c2(t) sono opportune funzioni del tempo le cui espressioni sono ricavate nel seguito. L’equazione (219), messa a sistema con l’equazione per la caratteristica negativa, fornisce la seguente espressione per la velocità

)(1

),()()(),(),(

2),(

),(

1

0*

1

2000

0

ttcga

ttsshttcttcttssvttssv

dgsfttssv

ga

tsv

Δ−−

Δ−Δ+−Δ−Δ−

−Δ−Δ+Δ−Δ+Δ

−Δ−Δ+=

(220)

Per quanto riguarda le funzioni c1(t) e c2(t), si consideri la quota piezometrica h*(s0,t) nella quale z è la quota della superficie libera e la pressione p è espressa mediante l’equazione (108).

γγγatm

kM

kM

atm ptzz

zzpp

tztptztsh −−

−+

+=+=))((

1)()()()()(),( 00

0* (221)

Dalla precedente si ha immediatamente:

dtzdt

dttshd

)(),( 0

*

ξ= (222)

in cui si è posto:

100

))(()(

)(1)(+−

−++= k

M

kMatm

tzzzzpp

tγ

κξ (223)

Utilizzando la (222) e ricordando che la portata entrante nella cassa d’aria è Qcassa=Ac(dz/dt) essendo Ac l’area della sezione orizzontale della cassa, l’equazione (218) può essere riscritta nella seguente forma linearizzata nel tempo

dttshd

tAA

AtQ

tsv cpp ),()(

1)(),( 0

*

0 ξ−= (224)

La precedente relazione è linearizzata nel tempo sostituendo ξ(t-Δt) al posto di ξ(t) e discretizzata alle differenze finite:

tttshtsh

ttAA

AtQ

tsv cp

ΔΔ−−

Δ−−=

),(),()(

1)(),( 0

*0

*

0 ξ (225)


Dal confronto tra la (219) e la (233) si ha immediatamente:

tAttA

ttc c

ΔΔ−−=Δ−

)()(1 ξ

AtQ

tAttttshA

ttc pc )()(

),()( 0

*

1 +ΔΔ−Δ−

=Δ−ξ

(226)

Ad ogni passo di calcolo va aggiornato il livello z, utilizzando la forma discreta e linearizzata dell’equazione (222):

( )),(),()(

1)()( 0*

0* ttshtsh

ttttztz Δ−−

Δ−+Δ−=

ξ (227)

Figura 3.20 – Andamento nel tempo del livello z nella cassa d’aria e della portata Q lungo la condotta

di mandata nel caso di stacco della pompa. Confronto tra le soluzioni ottenute in ipotesi elastiche ed anelastiche.


Per il caso di arresto istantaneo della pompa, la soluzione in ipotesi anelastiche è quella illustrata in Figura 2.15 e discussa nel paragrafo 2.3.3. Lo stesso problema, risolto in ipotesi elastiche mediante il procedimento appena descritto, fornisce i risultati illustrati in Figura 3.2 e qui posti a confronto con quelli ottenuti in ipotesi anelastiche.

Si osservano differenze assolutamente modeste tra i risultati in ipotesi elastiche ed anelastiche, gran parte dovute al procedimento numerico e al diverso passo di integrazione temporale (necessariamente molto piccolo se si opera in ipotesi elastiche). Ciò conferma la validità, per alcune tipologie di problemi, dell’approccio semplificato anelastico.

Reti di condotte

- 89 -

4 RETI DI CONDOTTE In questo capitolo saranno sinteticamente illustrati alcuni aspetti del problema

della verifica del funzionamento di sistemi di condotte in pressione. La verifica può richiedere sia la soluzione a moto permanente che quella in

condizioni di moto vario. Apparentemente, la soluzione a moto permanente è più rapida ed agevole; in realtà, a causa delle forti non linearità presenti nelle equazioni dinamiche, la soluzione richiede necessariamente un processo iterativo che, viceversa, non è necessario nella soluzione a moto vario. Molto spesso, quindi, si preferisce determinare la soluzione di moto permanente come soluzione asintotica di un problema di moto vario nel quale siano poste condizioni al contorno costanti nel tempo. Per questo motivo, in questo capitolo ci occuperemo della più generale soluzione di una rete in condizioni di moto non stazionario.

Sarà considerato inizialmente il funzionamento di una rete di condotte quando si possano assumere valide le ipotesi di comportamento anelastico. Si esaminerà successivamente lo stesso problema quando sia necessario (o conveniente) ricorrere alle ipotesi elastiche.

4.1 Reti di condotte trattate in ipotesi anelastiche Consideriamo una generica rete di condotte, come illustrato in Figura 4.1

(sinistra), nella quale non siano presenti dispositivi di diversa natura (serbatoi, pompe, ecc.) e organi di regolazione (saracinesche, ecc.): consideriamo cioè una rete costituita da sole condotte in pressione.

Figura 4.1 – Esempio di una rete di condotte (sinistra) e sua rappresentazione schematica (destra).

Tale sistema può essere schematizzato come una serie di tronchi di condotta connessi tra loro mediante nodi (Figura 4.1 a destra). I nodi della rete sono posti necessariamente in corrispondenza dell’intersezione tra tre o più condotte e in corrispondenza dei punti dove andranno poste le condizioni al contorno. E’ però possibile introdurre nel sistema nodi addizionali (ad esempio i nodi A, B, C di Figura

Reti di condotte

4.1) laddove vi siano variazioni di sezione o per qualsiasi altra esigenza, anche semplicemente grafica.

Le variabili del problema sono costituite dai valori della quota piezometrica in corrispondenza delle due estremità di ciascun tronco di condotta e dalla portata fluente. Nel complesso quindi il numero di incognite è 3M, essendo M il numero di tronchi di condotta che costituiscono la rete. Per la determinazione di queste variabili possono essere scritte M equazioni dinamiche (una per ciascun tronco della rete) e 2M equazioni di “congruenza” in corrispondenza ai nodi della rete.

Tali condizioni di congruenza sono rappresentate da un’equazione di continuità (la somma delle portate entranti nel nodo deve essere uguale alla somma delle portate uscenti) e da m-1 equazioni che esprimano un legame tra le quote piezometriche che si instaurano in corrispondenza delle sezioni terminali delle m condotte che convergono nel nodo.

In generale, e per semplicità, si assume che in un intorno del generico nodo la quota piezometrica si mantenga la stessa cosicché le m-1 condizioni di congruenza esprimono l’uguaglianza tra le quote piezometriche di estremità prese a due per due. Con riferimento all’esempio illustrato in Figura 4.2 (a sinistra), per il quale è m=4, le condizioni nodali sono espresse dalle seguenti relazioni:

Figura 4.2 – Illustrazione schematica delle condizioni di congruenza in un nodo della rete.

0=−++ dcba QQQQ (228)

******dccbba hhhhhh === (229)

Nella rappresentazione schematica della rete (Figura 4.2, a destra), le estremità delle condotte che convergono in un nodo vengono fatte coincidere con il nodo stesso: la quota piezometrica nel nodo n corrisponde quindi alle quote piezometriche alle estremità a, b, c, d delle quattro condotte che convergono nel nodo n. Così facendo le condizioni di congruenza sulle quote piezometriche (229) sono implicitamente imposte e l’unica condizione di congruenza necessaria è rappresentata dall’equazione di continuità (228).

In definitiva, con questa schematizzazione del funzionamento della rete, le variabili incognite sono rappresentate dalle M portate fluenti lungo i tronchi di condotta e dalle N quote piezometriche in corrispondenza degli N nodi della rete.

Reti di condotte

Per queste M+N incognite è possibile costruire un sistema composto da N equazioni di continuità (228) e da M equazioni dinamiche, una per ciascun tronco di condotta, come descritto nel seguito.

Con riferimento ad un generico tronco di condotta, compreso tra i nodi i e j, l’equazione dinamica (3), integrata nello spazio, si scrive

Ljdtdv

gLEE ij −−=− (230)

in cui la lunghezza L, la velocità v e la dissipazione di energia per unità di lunghezza j si riferiscono, ovviamente, al tratto di condotta in esame.

Essendo costante la sezione del tratto di condotta tra i nodi i e j, per l’equazione di continuità sarà costante anche il carico cinetico. Essendo inoltre v=Q/A, l’equazione (230) può essere riscritta come segue

LjdtdQ

gALhh ij −−=− ** (231)

Il termine jL che compare nella (231) rappresenta la perdita complessiva (perdita continua ed eventuali dissipazioni localizzate) che si realizza lungo tronco di condotta e che, come è noto, dipende dalla portata. In forma esplicita, l’equazione dinamica si può quindi scrivere:

QQdtdQ

gALhh ij α−−=− ** (232)

in cui α è un opportuno coefficiente che, per il caso di moto turbolento di parete scabra, dipende dalle caratteristiche geometriche della condotta e dalla scabrezza, mentre, più in generale, dipende anche dalla viscosità e dalla velocità e quindi dalla stessa portata. Con riferimento alle formule più frequentemente utilizzate nella pratica, le espressioni per il parametro α sono:

221

4 AgRLf

H

=α (Darcy-Weisbach) 23/42 ARkL

HS

=α (Gauckler-Strickler)

Per poter essere risolta, l’equazione (232) va linearizzata e discretizzata, cioè riscritta in forma algebrica. Tra le diverse possibilità, un modo semplice ed efficace consiste nel discretizzare la derivata temporale alle differenze finite e linearizzare nel tempo il termine che descrive le perdite:

QQtQQ

gALhh ij '''** α−

Δ−

−=− (233)

in cui Δt è il passo di integrazione temporale e con l’apice sono denotate le grandezze riferite al passo temporale precedente.

L’equazione (233) può essere riscritta come segue:

tAgQLQQ

tgALhh ij Δ

−⎥⎦

⎤⎢⎣

⎡+

Δ−=−

'''** α (234)

in modo da evidenziare il legame lineare tra le variabili *ih , *

jh e Q.

Reti di condotte

La risoluzione, ad ogni passo temporale Δt, del sistema composto dalle M equazioni (234) e dalle N equazioni (228), unitamente alle condizioni al contorno ed iniziali, fornisce la soluzione del problema.

E’ importante osservare che, per il fatto che le incognite del problema non sono omogenee (sono incognite sia le portate sia le quote piezometriche), la matrice del sistema presenta valori dei coefficienti che possono essere anche sensibilmente diversi tra loro. Con riferimento all’equazione (234), ad esempio, mentre i coefficienti delle quote piezometriche sono unitari, il coefficiente della portata, cioè il termine tra parentesi quadre, può assumere valori anche dell’ordine di 105-106. Tale disomogeneità ha un’influenza negativa sulla precisione del calcolo che, soprattutto in sistemi complessi, può portare a soluzioni sensibilmente errate.

Per questo motivo è conveniente impostare il problema in modo leggermente diverso.

La portata Q, nell’equazione (234), è la portata fluente lungo la condotta che va dal nodo i al nodo j ed è quindi, al tempo stesso, la portata che entra nel nodo j e che esce dal nodo i. In altre parole, dalla (234), la portata Qj entrante nel nodo j può essere espressa in funzione delle quote piezometriche *

ih e *jh come segue:

ijjijiijj ChBhBQ +−= ** con

tAgQLBCQ

tgALB ijijij Δ

=⎥⎦

⎤⎢⎣

⎡+

Δ=

−';''

1

α (235)

mentre la portata Qi, entrante nel nodo i, è semplicemente Qi=-Qj:

ijjijiiji ChBhBQ −+−= ** (236)

Consideriamo ora l’equazione di continuità di nodo (228). Se nello scrivere questa condizione sostituiamo alla generica portata entrante/uscente la corrispondente espressione del tipo (235) o (236), otteniamo, per il generico nodo n, un’equazione algebrica lineare contenente, come incognite, la quota piezometrica del nodo n e le quote piezometriche dei nodi che costituiscono l’altro estremo dei tronchi di condotta che convergono in n.

Possiamo in tal modo ridurci ad un sistema di N equazioni (di nodo) contenenti le sole quote piezometriche come incognite. Così facendo la matrice dei coefficienti risulta più omogenea.

Determinate le quote piezometriche nodali, attraverso la soluzione del sistema, è quindi possibile calcolare le portate in ciascun tronco della rete, utilizzando la (235), semplicemente per sostituzione a ritroso.

E’ evidente l’ulteriore vantaggio di questo modo di procedere. Infatti, invece di dover risolvere (ad ogni passo temporale) un sistema di M+N equazioni, è sufficiente risolvere un sistema di sole N equazioni ed applicare poi, M volte, l’equazione (235).

Dal punto di vista operativo, si tratta di calcolare, per ogni tronco di condotta, le seguenti matrici locali (dei coefficienti e dei termini noti)

⎥⎦

⎤⎢⎣

⎡−

−=

ijij

ijijij BB

BBB ⎥

⎦

⎤⎢⎣

⎡−=

ij

ijij C

CC (237)

Reti di condotte

che discendono direttamente dalle (236), (235) e rappresentano le portate entranti nei nodi i e j espresse in funzione delle corrispondenti quote piezometriche. E assemblare poi le M matrici locali in quella globale e nel vettore dei termini noti ottenendo il sistema risolutivo finale

0* =+ CB h (238)

4.1.1 Condizioni iniziali e al contorno

Una volta scritto il sistema di equazioni, ottenuto assemblando le matrici locali, è necessario assegnare le condizioni iniziali e al contorno.

Condizioni iniziali. Le condizioni iniziali sono costituite dai valori di quota piezometrica in ciascun nodo e dai valori di portata in ciascun tronco della rete all’istante iniziale t=0.

E’ opportuno che la distribuzione di portate e quote piezometriche sia dinamicamente in equilibrio in modo da evitare la formazione e propagazione di onde di pressione determinate da una condizione iniziale squilibrata. Tra le possibili scelte, generalmente si adotta o una condizione di quiete per il sistema o una condizione di moto stazionario.

Nel primo caso a ciascun nodo della rete viene assegnata la stessa quota piezometrica (il cielo piezometrico risulta pertanto una superficie orizzontale) e a ciascun tronco della rete viene assegnata una portata iniziale nulla.

Nel secondo caso quote piezometriche e portate assegnate corrispondono ad una soluzione di moto stazionario. Per altro, tale soluzione viene generalmente determinata mediante una simulazione preliminare per la quale, a partire da una condizione iniziale di quiete, si assegnano condizioni al contorno costanti nel tempo e la simulazione viene protratta nel tempo finchè il sistema non si porta a regime.

Condizioni al contorno. E’ possibile assegnare, in corrispondenza di alcuni nodi della rete, una, ed una sola, delle seguenti condizioni al contorno. Con riferimento al generico nodo n, è possibile assegnare:

• l’andamento nel tempo della quota piezometrica: in tal caso, dal punto di vista operativo, ad ogni passo temporale di calcolo viene azzerata la riga n-esima della matrice dei coefficienti e si pone un coefficiente unitario in corrispondenza della diagonale principale. Al tempo stesso, al termine n-esimo del vettore dei termini noti viene assegnato il valore della quota piezometrica imposta (Cioè, al termine n-esimo del vettore C viene assegnato il valore della quota piezometrica imposta cambiato di segno);

• l’andamento nel tempo della portata entrante nel nodo n (per assegnare una portata uscente è sufficiente fornire valori negativi): è sufficiente, in questo caso, sommare al termine n-esimo del vettore dei termini noti il valore della portata assegnata;

• un legame lineare (o linearizzato) tra portata entrante e quota piezometrica: esprimiamo tale legame mediante la relazione: nnnn ChBQ += * , in cui Qn è la

Reti di condotte

portata entrante nel nodo, *nh la corrispondente quota piezometrica e Bn, Cn,

opportuni coefficienti. In tal caso la condizione viene imposta sommando Bn alla diagonale principale della matrice dei coefficienti in corrispondenza della riga n-esima, e Cn al termine n-esimo del vettore dei termini noti.

A conclusione di questo paragrafo, è solo il caso di sottolineare che, almeno in un

nodo della rete, è necessario assegnare una condizione al contorno del tipo “andamento nel tempo della quota piezometrica” o del tipo “legame lineare tra portata entrante e quota piezometrica”. In caso contrario il sistema ammette infinite soluzioni (un’unica distribuzione di portate e infinite distribuzioni di quote piezometriche che differiscono per un valore costante: la quota di riferimento)

4.1.2 Dispositivi ed organi di regolazione

In generale, in una rete di condotte sono presenti diversi dispositivi e organi di regolazione. Questi possono essere inclusi nel sistema come “tronchi di condotta a funzionamento speciale” o “nodi a funzionamento speciale”. Nel seguito vengono riportati solo alcuni esempi per illustrare la tecnica normalmente impiegata per descrivere questi dispositivi.

Dissipazione localizzata. Con riferimento allo schema di Figura 3.4 (sinistra), supponiamo che, ad esempio lungo un tronco della rete, sia presente un dispositivo che introduce una dissipazione localizzata di energia. In tal caso, nel tratto di rete in esame viene introdotto un tronco speciale (e due ulteriori nodi: i, j) come illustrato in Figura 3.4 (destra).

Figura 4.3

Per il nuovo tronco i-j di condotta possiamo scrivere:

QQAgQQ

hAgQQ

hEEEj

ji

ijiij κ=−−+=−=Δ 2*

2*

22

in cui il coefficiente κ dipende sostanzialmente dalle caratteristiche geometriche del dispositivo che determina la dissipazione localizzata.

Linearizzando nel tempo, si ha:

QQAAg

Qhh

ijji

⎥⎥⎦

⎤

⎢⎢⎣

⎡+⎟

⎟⎠

⎞⎜⎜⎝

⎛−=− '11

2'

22** κ

La precedente relazione può essere riscritta formalmente come la (235) e quindi, come la (235), assemblata nel sistema risolutivo:

Reti di condotte

ijjijiijj ChBhBQ +−= ** con 0;'112

'1

22 =⎥⎥⎦

⎤

⎢⎢⎣

⎡+⎟

⎟⎠

⎞⎜⎜⎝

⎛−=

−

ijij

ij CQAAg

QB κ (239)

In questo caso, dal punto di vista numerico, bisogna evitare che si annulli l’argomento tra le parentesi quadrate, ponendo artificialmente un limite inferiore al modulo della portata fluente al passo di calcolo precedente.

Nodo serbatoio. Per simulare la presenza di un serbatoio è possibile introdurre un “nodo speciale”. Di fatto, come si vedrà tra breve, un serbatoio rappresenta una immissione esterna di portata dipendente dalla quota piezometrica ed è quindi riconducibile ad una condizione al contorno del tipo “legame lineare tra portata entrante e quota piezometrica”.

Figura 4.4 – Schematizzazione di un serbatoio.

Consideriamo la situazione illustrata in Figura 4.4 nella quale, in corrispondenza del nodo n della rete è posto un serbatoio caratterizzato da una sezione orizzontale di area Ω (costante per semplicità). Con riferimento all’equazione dei serbatoi (19) si osserva che la portata netta che esce dal serbatoio (Qusc-Qentr) corrisponde alla portata Qn che dal serbatoio entra nel nodo n. Possiamo allora scrivere:

dtdh

dtdhQ n

n

*

Ω−=Ω−=

Nella precedente, trascurando eventuali dissipazioni di energia, si è assunto che il livello h nel serbatoio coincida con la quota piezometrica *

nh al piede del serbatoio. Discretizzando la derivata temporale alle differenze finite, possiamo scrivere:

( )thhQ nn

n Δ−

Ω−≅'**

in cui con l’apice è denotata la quota piezometrica del nodo n relativa al passo di calcolo precedente. La relazione appena scritta può essere posta nella forma che abbiamo visto per la condizione al contorno di tipo “legame lineare tra portata entrante e quota piezometrica”; si ha infatti:

nnnn ChBQ +≅ * con ( )'*; nnn ht

Ct

BΔΩ

=ΔΩ

−= (240)

Volendo introdurre eventuali dissipazioni localizzate di energia, in prima approssimazione è possibile includere il relativo tronco speciale come schematicamente illustrato in Figura 4.5.

Reti di condotte

Figura 4.5

Esempio. Consideriamo il problema dell’oscillazione di massa in un pozzo piezometrico visto nel paragrafo 2.3.2.

Il sistema può essere schematizzato come illustrato in Figura 4.6.

Figura 4.6 – Rappresentazione schematica del sistema illustrato in Figura 2.7.

Il sistema è composto da tre nodi (1, 2 e 3) di cui il terzo è un nodo speciale “serbatoio”, da un solo tronco di condotta, tra i nodi 1 e 2 (galleria di derivazione), e da un tronco speciale “dissipazione localizzata” tra i nodi 2 e 3. Nel nodo 1 è posta una condizione al contorno del tipo “quota piezometrica assegnata” che corrisponde al livello (costante) del serbatoio di monte. Nel nodo 2 è posta una condizione al contorno del tipo “portata assegnata” che corrisponde alla portata che scende attraverso la condotta forzata.

Volendo simulare, ad esempio, una manovra di apertura per la quale la portata Q2(t) derivata verso la condotta forzata è inizialmente nulla (Q2(t=0)=0), la condizione iniziale è data da:

0*3

*2

*1 === hhh 02312 == QQ

in cui il riferimento assunto coincide, evidentemente, con il livello del serbatoio di monte.

Ad ogni passo di calcolo vanno innanzitutto assemblati i contributi dei due tronchi di condotta (normale e speciale). Con riferimento al tronco 1-2 si ha:

matrici locali ⎥⎦⎤

⎢⎣⎡

−−

1212

1212

BBBB ⎥⎦

⎤⎢⎣⎡−

12

12

CC

matrici globali ⎥⎥

⎦

⎤

⎢⎢

⎣

⎡−

−

00000

1212

1212

BBBB

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡−

012

12

CC

Reti di condotte

in cui i coefficienti B12 e C12 sono dati dalla (235). Assemblando quindi il tronco 2-3 si ha:


⎢⎣⎡

−−

2323

2323

BBBB ⎥⎦

⎤⎢⎣⎡−

23

23

CC


⎦

⎤

⎢⎢

⎣

⎡

−−−

−

2323

23231212

1212

0

0

BBBBBB

BB

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡−

−

23

2312

12

CCC

C

in cui i coefficienti B23 e C23 sono dati dalla (239). Si assembla ora il contributo del nodo serbatoio rappresentato dalla portata Q3

entrante nel nodo 3 ed espressa dalle equazioni (240). La matrice globale diventa quindi:


⎦

⎤

⎢⎢

⎣

⎡

+−−−

−

32323

23231212

1212

0

0

BBBBBBB

BB

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

+−

−

323

2312

12

CCCC

C

Imponendo la condizione al contorno di livello assegnato sul nodo 1 ( 0*1 =h ) si ha:


⎦

⎤

⎢⎢

⎣

⎡

+−−−

32323

23231212

0

001

BBBBBBB

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

+−

323

2312

0

CCCC

Imponendo infine, in corrispondenza del nodo 2 la condizione al contorno sulla portata Q2(t) (che avrà evidentemente un valore negativo in quanto esce dal nodo) si avrà:


⎦

⎤

⎢⎢

⎣

⎡

+−−−

32323

23231212

0

001

BBBBBBB

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

++−

323

22312 )(0

CCtQCC i

in cui ti è l’istante di calcolo attuale. La soluzione del precedente sistema fornisce, all’istante attuale, le quote

piezometriche nei tre nodi della rete (livello nel serbatoio di monte, quota piezometrica al piede del pozzo e livello nel pozzo piezometrico). Sostituendo (a ritroso) i valori di quota piezometrica trovati, nelle due equazioni di tronco precedentemente assemblate (e corrispondenti alle equazioni (235) e (239)), si determinano i valori di portata:

12*212

*11212 ChBhBQ +−=

23*323

*22323 ChBhBQ j +−=

(241)

Si procede quindi a risolvere il passo temporale successivo

Reti di condotte

Nel caso particolare di questo esempio il sistema risolutivo è estremamente semplice è può essere risolto agevolmente per via analitica. In particolare, essendo

*1h noto, il sistema si riduce a due sole equazioni

0)(323

*11222312

*3

*2

32323

232312 =⎥⎦

⎤⎢⎣

⎡+

++−+⎥⎦

⎤⎢⎣

⎡⋅⎥⎦

⎤⎢⎣⎡

+−−−

CChBtQCC

hh

BBBBBB i

Si ricava, in particolare

( ) ( )( )

2312

*323

*11222312*

2

231232312

233123*11221223*

3

)(

)(

BBhBhBtQCC

h

BBBBBCCBChBtQCB

h

i

i

++++−

=

+−+++++

=

Figura 4.7 – Andamento nel tempo del livello nel pozzo piezometrico per manovre di apertura (in alto)

e chiusura (in basso) e confronto con la soluzione numerica ottenuta mediante un metodo esplicito ed illustrata nel paragrafo 2.3.2.

Reti di condotte

Con riferimento ai dati dell’esempio riportato nel paragrafo 2.3.2, l’applicazione del procedimento appena illustrato fornisce, con un passo di integrazione temporale Δt=1.0s, i risultati illustrati in Figura 4.7 (il livello z nel pozzo corrisponde alla quota piezometrica *

3h ). Nella stessa figura, per confronto, sono riportati i risultati ottenuti con il metodo

esplicito e un passo di integrazione temporale Δt=0.2s. E’ evidente la maggiore accuratezza che si ottiene col procedimento illustrato in

questo capitolo. In particolare, nel caso di una manovra di chiusura, l’ampiezza delle oscillazioni del livello nel pozzo piezometrico non soffrono dell’instabilità evidenziata per il metodo esplicito e si smorzano fino ad annullarsi completamente.

Pompa. L’inserimento di una pompa nel sistema può essere fatto mediante l’impiego di un tronco a funzionamento speciale “pompa” (Figura 4.8). Analiticamente, il funzionamento di tale tronco sarà rappresentato ancora da un’equazione lineare esprimente un legame tra la portata fluente nel tronco e le quote piezometriche agli estremi, formalmente identica alla (235) e quindi, come la (235), assemblata nel sistema risolutivo.

Figura 4.8

In relazione all’accuratezza e al dettaglio richiesti, l’equazione per il “tronco pompa” può essere formulata in diversi modi.

Nel caso più semplice si può ipotizzare che la portata pompata assuma un valore costante, QP (corrispondente, ad esempio, alla portata nominale della pompa), indipendente cioè dalle quote piezometriche nei nodi estremi, ovvero dalla prevalenza.

In alternativa si può far riferimento alla curva caratteristica della pompa che lega la portata sollevata QP alla prevalenza H. Quest’ultima, con buona approssimazione, è data dalla differenza **

ij hhH −= tra le quote piezometriche agli estremi del tronco. Facendo riferimento ai valori assunti dalle quote piezometriche all’istante di calcolo precedente è possibile stimare, ad ogni passo di calcolo e in modo esplicito, l’effettiva portata sollevata QP mediante la curva caratteristica anzidetta.

In entrambi i casi l’equazione per il “tronco pompa” si scrive:

ijjijiijj ChBhBQ +−= ** con Pijij QCB == ;0 (242)

Nel secondo caso, oltre al problema legato alla necessità di controllare costantemente che i valori della quota piezometrica si mantengano all’interno dell’intervallo per cui è disponibile il legame QP-H, c’è il pericolo che la soluzione numerica presenti fastidiose oscillazioni causate dalla determinazione esplicita della portata sollevata.

Reti di condotte

Per evitare queste oscillazioni si può ricorrere ad una formulazione implicita solo formalmente un po’ più complessa.

Esprimiamo, in forma generale, il legame tra la portata sollevata e la prevalenza come: QP=f(H) (generalmente tale legame è ben approssimato da un polinomio di secondo grado), e sviluppiamo tale legame in serie di Taylor, arrestata al primo ordine, in un intorno del valore H’ assunto dalla quota piezometrica all’istante di calcolo precedente:

)'()'('

HHHfHfQ

HHP −

∂∂

+≅=

(243)

Ricordando che **ij hhH −= , l’equazione per il “tronco pompa” si scrive:

ijjijiijj ChBhBQ +−= ** con ';)'('

HBCHfHfB ijij

HHij =

∂∂

−−==

(244)

Valvola clapét. La valvola clapét è un dispositivo che consente il flusso in un verso (dal nodo i verso il nodo j in Figura 4.9) ed impedisce il flusso nel verso opposto. Tale dispositivo può essere schematizzato mediante un tronco a funzionamento speciale, “tronco clapét”, per il quale l’equazione associata viene a dipendere dallo stato (aperto/chiuso) della valvola.

Figura 4.9

Quando il clapét è aperto, il dispositivo può essere assimilato ad una piccola dissipazione di energia localizzata e la sua equazione coincide pertanto con la (239), la quale, per essere Ai=Aj, si semplifica nella seguente:

ijjijiijj ChBhBQ +−= ** con 0;'

1=−= ijij C

QB

κ (245)

Al coefficiente vengono generalmente assegnati valori molto bassi essendo piccola la dissipazione prodotta dalla valvola aperta. Si osservi che non è possibile assumere una dissipazione nulla (κ=0), in questo caso infatti la portata fluente sarebbe indipendente dalla variazione di quota piezometrica.

Quando il clapét è chiuso si determina una sconnessione tra i nodi i e j. In questo caso il tronco non va assemblato o, che è lo stesso, la sua equazione diventa:

ijjijiijj ChBhBQ +−= ** con 0== ijij CB (246)

L’aspetto più interessante, e tutt’altro che banale, riguarda il criterio da adottare per stabilire se la valvola è aperta oppure chiusa. La regola generalmente adottata è la seguente.

Reti di condotte

Se al passo di calcolo precedente la valvola era aperta si controlla il verso della portata: se questo è positivo (cioè il flusso avviene dal nodo i verso il nodo j), la valvola si mantiene aperta, altrimenti si considera chiusa.

Se, invece, al passo di calcolo precedente la valvola era chiusa si controllano le quote piezometriche ai nodi di estremità: se **

ij hh > , la valvola rimane chiusa altrimenti, si considera aperta.

E’ necessario pertanto associare al “tronco clapét” un flag che ne descriva lo stato.

Cassa d’aria. Per simulare la presenza di una cassa d’aria è possibile introdurre un “nodo speciale”. Di fatto, come si vedrà tra breve, una cassa d’aria può essere schematizzata in modo analogo a quanto fatto per un serbatoio. Essa, infatti, rappresenta una immissione esterna di portata dipendente dalla quota piezometrica ed è quindi riconducibile ad una condizione al contorno del tipo “legame lineare tra portata entrante e quota piezometrica”.

Figura 4.10 – Schematizzazione di una cassa d’aria.

Consideriamo la situazione illustrata in Figura 4.10 nella quale, in corrispondenza del nodo n della rete è posta una cassa d’aria caratterizzata da una sezione orizzontale di area Ac (costante per semplicità). Con riferimento all’equazione dei serbatoi (19) si osserva che la portata netta che esce dalla cassa (Qusc-Qentr) corrisponde alla portata Qn che dalla cassa entra nel nodo n. Possiamo allora scrivere:

dtdzAQ cn −= (247)

in cui z è la quota della superficie libera nella cassa. Il livello z e la quota piezometrica *

nh sono tra loro legate come espresso dall’equazione (222). Linearizzando nel tempo tale legame e discretizzando la derivata temporale alle differenze finite, si può scrivere

( )thhA

Q nncn Δ

−−≅

'**

'ξ

in cui con l’apice sono state denotate quantità relative al passo di calcolo precedente. La relazione appena scritta può essere posta nella forma che abbiamo visto per la condizione al contorno di tipo “legame lineare tra portata entrante e quota piezometrica”; si ha infatti:

nnnn ChBQ +≅ * con ( )'*

';

' nc

nc

n ht

AC

tA

BΔ

=Δ

−=ξξ

(248)

Reti di condotte

Ovviamente durante il calcolo dev’essere aggiornata anche la quota z dell’acqua nella cassa. A tale scopo può essere utilizzata l’equazione di continuità (247) discretizzata ed esplicitata rispetto a z come segue

c

n

AQtzz Δ−= ' (249)

Restringimento/ugello. La presenza di un restringimento/allargamento ben raccordato o di un ugello terminale può essere simulata mediante l’impiego di un tronco a comportamento speciale: “tronco restringimento” (Figura 4.11) la cui equazione, essendo trascurabile l’inerzia, si basa su di un semplice bilancio di energia nel quale si trascura la dissipazione di energia localizzata.

Figura 4.11

Essendo:

2*

2*

22 jj

ii Ag

QQh

AgQQ

h +=+ (250)

si ha:

ijjijiijj ChBhBQ +−= ** con 0;11'

21

22 =⎟⎟⎠

⎞⎜⎜⎝

⎛−=

−

ijij

ij CAAQ

gB (251)

Anche in questo caso, come per l’equazione (239), bisogna porre artificialmente un limite inferiore al modulo della portata fluente al passo di calcolo precedente.

Reti di condotte

4.2 Reti di condotte trattate in ipotesi elastiche Nel caso sia necessario ricorrere alle ipotesi elastiche il problema della soluzione

di una rete in pressione diventa solo formalmente più complesso. La soluzione, da un punto di vista generale, si determina ancora strutturando un sistema di equazioni formalmente analogo a quello visto per il caso di moto in ipotesi anelastiche. In particolare le equazioni per i tronchi e per i nodi a funzionamento speciale restano le stesse viste nel paragrafo precedente. La differenza dal caso anelastico risiede sostanzialmente nel modo di scrivere le equazioni per le diverse condotte che compongono la rete.

Consideriamo un generico tronco di condotta, lungo L, collegante i nodi i e j della rete. Per risolvere il campo di moto lungo questa condotta è necessario (si veda il capitolo 3.3) suddividere la condotta in n tratti di lunghezza Δs in modo che sia soddisfatta la relazione Δs=a·Δt, essendo a la celerità di propagazione di una perturbazione di pressione (vedi equazione (187))

Figura 4.12 – Illustrazione schematica del procedimento di integrazione nel tempo delle equazioni

A partire dalla condizione nota relativa all’istante di calcolo precedente, mediante l’impiego delle equazioni per le caratteristiche C+ e C- è possibile determinare, come schematicamente illustrato in Figura 4.12, la soluzione all’istante di calcolo attuale per tutti i punti interni della condotta. Per questa valutazione si utilizzano, in particolare, le equazioni (200) e (201) che, riscritte in termini di portata, sono:

[ ] [ ]

[ ]),(),(),(),(4

),(),(2

),(),(21),(

010101012

010101*

01*

0*

tsQtsQtsQtsQdAgsf

tsQtsQgAatshtshttsh

kkkk

kkkkk

++−−

+−+−

−Δ

−

−−++=Δ+ (252)

Reti di condotte

[ ] [ ]

[ ]),(),(),(),(4

),(),(2

),(),(21),(

01010101

01*

01*

01010

tsQtsQtsQtsQAda

sf

tshtshaAgtsQtsQttsQ

kkkk

kkkkk

++−−

+−+−

+Δ

−

−−++=Δ+ (253)

Effettuato preliminarmente questo calcolo per tutte le condotte della rete, si passa alla scrittura (e quindi all’assemblaggio) delle equazioni di nodo la cui forma è ricavata qui di seguito.

Consideriamo dapprima l’estremità corrispondente al nodo j (s=L) e l’equazione per la caratteristica positiva (198), qui riscritta in termini di portata Qj entrante nel nodo j:

[ ] sAg

tsQtsQdftsQQ

Agatshh nn

njnj Δ−−−=− −−−− 2

01010101

**

2),(),(

),(),( (254)

in cui, per comodità e analogia con la scrittura utilizzata nel precedente paragrafo, si è indicato con *

jh il termine ),( 0* ttsh n Δ+ e con Qj il termine ),( 0 ttsQ n Δ+ .

Posto:

),(2

),(),(),(;;0 012

010101

* tsQsAg

tsQtsQdftsh

aAgC

aAgBB n

nnnjjjjji −

−−− +

⎥⎥⎦

⎤

⎢⎢⎣

⎡Δ−=−== (255)

l’equazione (254) può essere riscritta come segue:

jjjjjijij ChBhBQ ++= ** (256)

Consideriamo ora l’estremità corrispondente al nodo i (s=0) e l’equazione per la caratteristica negativa (199), qui riscritta in termini di portata Qi entrante nel nodo i:

[ ] sAg

tsQtsQdftsQQ

Agatshh ii Δ+−−=− 2

01010101

**

2),(),(

),(),( (257)

in cui, come prima, si è indicato con *ih il termine ),( 00

* ttsh Δ+ e con Qi il termine ),( 00 ttsQ Δ+− .

Posto:

),(2

),(),(),(;0; 012

010101

* tsQsAg

tsQtsQdftsh

aAgCB

aAgB niiijii −+

⎥⎥⎦

⎤

⎢⎢⎣

⎡Δ+==−= (258)

l’equazione (257) può essere riscritta come segue:

iijijiiii ChBhBQ ++= ** (259)

Nel complesso, quindi, le matrici locali per la generica condotta (corrispondenti alle (237) per il caso di moto anelastico) sono:

Reti di condotte

⎥⎦

⎤⎢⎣

⎡

jjji

ijii

BBBB

⎥⎦

⎤⎢⎣

⎡jj

ii

CC (260)

Va osservato che, essendo Bij=Bji=0, la matrice globale dei coefficienti si mantiene simmetrica come per il caso anelastico.

Infine, il modo di porre le condizioni al contorno è lo stesso visto per il caso anelastico mentre un po’ più complessa è la procedura da seguire per porre le condizioni iniziali in quanto non è sufficiente assegnare le quote piezometriche in tutti i nodi della rete e le portate nei tronchi ma è necessario fornire quote piezometriche e portate anche in tutti punti interni (in relazione alla discretizzazione spaziale in tratti Δs) di ciascuna condotta. Esempio. Consideriamo il problema dell’oscillazione di massa in una cassa d’aria visto nel paragrafo 2.3.3 (vedi anche 3.3.5).

Il sistema può essere schematizzato come illustrato in Figura 4.13.

Figura 4.13 – Rappresentazione schematica del sistema illustrato in Figura 2.12.

Il sistema è composto da due nodi (1 e 2) di cui il primo è un nodo speciale “cassa d’aria” e da un solo tronco di condotta, tra i nodi 1 e 2 (condotta di mandata). Nel nodo 2 è posta una condizione al contorno del tipo “quota piezometrica assegnata” che corrisponde al livello (costante) del serbatoio di recapito. Nel nodo 1 è posta una condizione al contorno del tipo “portata assegnata” che corrisponde alla portata pompata nel sistema (in questo caso non è necessario utilizzare un tronco speciale “pompa” in quanto la pompa è posta in corrispondenza di un nodo di estremità della rete.

Volendo simulare, ad esempio, una manovra di arresto della pompa per la quale la portata Q1(t) pompata nel sistema vale inizialmente 0.05 m3/s, la condizione iniziale è data da (si veda Figura 2.15):

mhmh 0.4091.52 *2

*1 == smQQ /05.0 3

2,121,12 ==

in cui il riferimento assunto coincide con il livello nella cassa d’aria in condizioni di sistema fermo (vedi paragrafo 2.3.3 e 3.3.5). Dovendo necessariamente seguire passo passo anche l’evoluzione temporale del livello z nella cassa d’aria è necessario fissare una condizione iniziale anche per tale parametro. Con riferimento all’esempio riportato nel paragrafo 2.3.3 si ha:

mz 149.0=

Reti di condotte

Si ricorda infine che, volendo risolvere il moto nel sistema seguendo un approccio elastico, è necessario fissare quote piezometriche e portate in ogni sezione in cui si suddivide ciascuna condotta. Nel caso in esame, con riferimento all’esempio illustrato nel paragrafo 3.3.5, immaginando di suddividere la condotta di mandata in 10 tronchi di lunghezza Δs=100m, le quote piezometriche iniziali nelle diverse sezioni interne sono fornite dalla relazione:

212

212*

2*12 2 Ag

Qd

sjfhh jΔ

+= j=1, 2, 3, …, 9

mentre i valori di portata nelle sezioni interne sono tutti uguali e corrispondono a Q=0.05 m3/s.

Ad ogni passo di calcolo vanno innanzitutto assemblati i contributi della condotta di mandata. Con riferimento al tronco 1-2 si ha:


⎢⎣⎡

22

11

00

BB ⎥⎦

⎤⎢⎣⎡

22

11

CC

matrici globali ⎥⎦⎤

⎢⎣⎡

22

11

00

BB ⎥⎦

⎤⎢⎣⎡

22

11

CC

in cui i coefficienti B11, B22, C1 e C2 sono dati dalla (258). Assemblando quindi il contributo del nodo “cassa d’aria” rappresentato dalla

portata Q1 entrante nel nodo 1 ed espressa dall’equazione (248), la matrice globale diventa:


⎢⎣⎡ +

22

111

00B

BB ⎥⎦⎤

⎢⎣⎡ +

22

111

CCC

Imponendo la condizione al contorno di livello assegnato sul nodo 1 ( mh 0.40*2 = )

si ha:


⎢⎣⎡ +

100111 BB ⎥⎦

⎤⎢⎣⎡ +

40111 CC

Imponendo infine, in corrispondenza del nodo 1 la condizione al contorno sulla portata Q1(t), corrispondente alla portata pompata (che avrà evidentemente un valore positivo in quanto entra nel nodo) si avrà:


⎢⎣⎡ +

100111 BB ⎥⎦

⎤⎢⎣⎡ ++

40)(1111 tQCC

La soluzione del precedente sistema fornisce, all’istante attuale, le quote piezometriche nei due nodi della rete (livello nel serbatoio di valle e quota piezometrica in corrispondenza della cassa d’aria). Sostituendo (a ritroso) i valori di quota piezometrica trovati, nelle due equazioni di tronco precedentemente assemblate, si determinano i valori di portata alle due estremità della condotta di mandata:

Reti di condotte

11*1111,12 ChBQ +=

22*2222,12 ChBQ +=

(261)

Analogamente, utilizzando l’equazione (248) è possibile valutare la portata Q1 che esce dalla cassa d’aria ed entra nel nodo 1

1*111 ChBQ += (262)

Ovviamente durante il calcolo dev’essere aggiornata anche la quota z dell’acqua nella cassa. A tale scopo può essere utilizzata l’equazione di continuità (247) discretizzata ed esplicitata rispetto a z come segue

cAQtzz 1' Δ−= (263)

Si procede quindi a risolvere il passo temporale successivo Nel caso particolare di questo esempio il sistema risolutivo è estremamente

semplice essendo, di fatto, costituito da una sola equazione. Si ha, in particolare

0.40

)(

*2

111

1111*1

=+

++−=

hBB

tQCCh

L’applicazione del procedimento appena illustrato fornisce, con un passo di integrazione

Terminare esempio

- 109 -

ComplementiDiIdraulica new 02 WEB - image.unipd.it€¦ · Moto vario in ipotesi anelastiche...

Documents

Transcript of ComplementiDiIdraulica new 02 WEB - image.unipd.it€¦ · Moto vario in ipotesi anelastiche...