Algoritmo per il calcolo dei profili di corrente in presenza di manufatti in serie e studio delle soluzioni multiple

1 Introduzione
1.1 Visione d’insieme
Nella prima parte dell’elaborato viene illustrato il funzionamento del codice di calcolo svilup-
patonelcorsodeidueannidistudiopressolamagistralediIngegneriaperl’AmbienteeilTerritorio.
Si mostrerà come l’approccio utilizzato sia fedelmente basato sulla teoria dell’idrodinamica. Esso
distingue i manufatti (oggetti puntuali attraverso i quali si applica un bilancio di energia) dal resto
della canaletta, in cui viene integrata l’equazione dei profili. Si illustrerà l’algoritmo dinamico che
consente di determinare le influenze dei manufatti e che rende quindi il codice in grado di risol-
vere un problema del tutto generale di manufatti in serie. Infine mostrerà come è stato possibile
costruire un algoritmo in grado di risolvere un qualunque problema di Deflusso di portata e un
problema di Imbocco.
Nella seconda parte dell’elaborato viene approfondito il fenomeno poco conosciuto delle soluzioni
multiple allo stato stazionario che si possono realizzare nei problemi a pelo libero in un canale
rettangolare in presenza di manufatti come Luci di fondo, Restringimenti, Salti di fondo ecc...
Verrà illustrato un secondo algoritmo che consente di ricercare soluzioni alternative per il medes-
imo problema stazionario. Infine si parlerà anche del tema della stabilità delle soluzioni multiple,
e del fenomeno che regola l’instaurarsi di una soluzione piuttosto che di un’altra.
1.2 Oggetto della tesi
Il presente elaborato finale del corso di laurea magistrale in Ingegneria per l’Ambiente e il
Territorio tratta del moto a pelo libero nei canali di forma rettangolare allo stato stazionario.
In particolare lo studio si focalizza sul moto a pelo libero in presenza di manufatti "in serie" ossia
che interagiscono idro-dinamicamente l’uno con l’altro. É stato realizzato un codice di calcolo in
linguaggio Python nel corso dei due anni di magistrale che rappresenta il cuore della tesi.
Il codice sviluppato contiene un primo algoritmo che consente di trovare una soluzione del moto
in presenza contemporaneamente di: manufatti in serie, afflussi di portata, variazioni graduali di
larghezza e un deflusso di portata o un problema di Imbocco.
Per fare questo è stata utilizzata una programmazione orientata agli oggetti, ogni elemento sopra
elencato ha una classe associata che contiene i metodi necessari per la creazione della geometria e
delle condizioni al contorno e per la risoluzione del problema idrodinamico dell’oggetto stesso (nel
caso dei manufatti puntuali dato il tirante di moto indisturbato consiste nel trovare il tirante di
monte e di valle). In particolare è risultato particolarmente complesso creare un algoritmo che fosse
in grado di risolvere un qualunque problema di imbocco e ancor di più un problema di deflusso di
portata.
La risoluzione del problema di deflusso in caso di corrente lenta è molto complessa. Tale comp-
lessità nasce dal fatto che la condizione al contorno è dettata da valle e a valle non si conosce a
priori la portata, che è funzione dell’altezza del profilo di corrente che si instaura all’interno della
finestra di deflusso. In altre parole per conoscere la condizione al contorno di valle ho bisogno di
aver risolto il deflusso, ma per risolvere il deflusso ho bisogno di conoscere la condizione al contorno
di valle.
Per tale ragione la soluzione non può che trovarsi in maniera iterativa, ipotizzando una portata a
valle del deflusso, da questa ricavare il valore di tirante (cioè la C.C.) ( Per fare ciò è necessario
costruire la scala di deflusso della sezione di valle ) risolvere la finestra di deflusso e verificare a
posteriori l’ipotesi sulla portata di valle in base a quanta acqua viene sfiorata.
Un problema ancora più complesso si ha quando le condizioni al contorno sono due. Vi è cioè una
condizione di corrente veloce dettata da monte e una condizione di corrente lenta dettata da valle.
Anche in questo caso la soluzione non può che essere di natura iterativa, ma la maggior complessità
7
nasce nella costruzione della scala di deflusso della sezione di valle (può essere necessario ipotizzare
di avere transizione nei manufatti di valle) e nel determinare la condizione al contorno che monte
genera nella sezione iniziale della finestra di deflusso (anche in questo caso può essere necessario
ipotizzare di avere transizione nei manufatti di monte).
Anche il problema di Imbocco può essere particolarmente complesso nella costruzione della scala
di deflusso (può essere necessario fare ipotesi di transizione).
La capacità di risolvere un problema di deflusso e un problema di imbocco sono le principali abilità
in più che distinguono questo primo algoritmo da un software commerciale (Hec-ras).
Un’altra qualità del codice sviluppato nella presente attività di tesi è la praticità d’uso, cioè è uno
strumento molto più speditivo di Hec-ras per l’analisi di problemi stazionari in canali rettangolari.
In generale i canali rettangolari non esistono in natura, quindi si parla di canali costruiti dall’uomo.
É tipico trovare canali rettangolari con presenza di manufatti (anche in serie) in contesti urbaniz-
zati. Può essere usato anche come strumento di controllo di simulazioni su fiumi e torrenti naturali
che sono stati analizzati con software più complessi.
Inoltre è uno strumento molto teoricamente basato, che integra l’equazione delle correnti nelle zone
che non presentano manufatti, e utilizza invece un approccio energetico per risolvere quelle sezioni
che presentano manufatti. Uno speciale algoritmo dinamico consente di risolvere il moto nel canale
rettangolare anche in presenza di manufatti in serie.
Il codice contiene inoltre un secondo algoritmo (concettualmente diverso dal primo) che consente di
trovare tutte (per lo meno la maggior parte) le soluzioni ammesse allo stato stazionario in presenza
di manufatti in serie, variazioni graduali di larghezza, afflussi di portata.
Viene inoltre riportato nella tesi un articolo scientifico "Extended Theory of Hydraulic Hysteresis
in Open-Channel Flow" che parla della molteplicità di soluzioni e del fenomeno di isteresi nel caso
di un Salto di fondo. Questo articolo è stato tradotto dall’inglese, riassunto e riportato nella tesi
utilizzando una notazione per le formule che fosse conforme al resto della tesi.
Esso è stato molto utile perché ha contribuito a dare vita a una serie di riflessioni sulla stabilità
delle soluzioni multiple e sul fenomeno dell’isteresi.
Le soluzioni multiple sono un argomento di cui si parla e si conosce veramente poco in letteratura.
Eppure sono una realtà che può manifestarsi. Esse sono strettamente correlate alla presenza di
manufatti, che interrompono la continuità dell’equazione dei profili. Basta la presenza di un solo
manufatto per doversi domandare se più di una soluzione è ammissibile allo stato stazionario, ed è
sempre più probabile avere soluzioni multiple al crescere del numero dei manufatti, in particolare
se interagiscono tra loro idro dinamicamente.
1.3 struttura dell’elaborato
Nella prima parte dell’elaborato si fanno dei richiami teorici sulle equazioni che regolano il
moto a pelo libero dell’acqua, partendo dalle equazioni in Shallow water per una corrente monodi-
mensionale (o equazione di De Saint Venant) per arrivare fino all’equazione dei profili in canale
rettangolare. Seguono dei richiami teorici sul metodo di Eulero, utilizzato per l’integrazione nu-
merica di tale equazione differenziale.
In seguito si fa una panoramica globale del modello di calcolo (elementi dell’interfaccia grafica,
introduzione agli algoritmi, potenzialità del programma). In particolare si spiega il funzionamento
dell’algoritmo associato al pulsante START SIMULATION.
Successivamente si ha una trattazione teorica (qualora sia necessaria) e pratica (cioè come vengono
gestiti dal codice) dei singoli oggetti inseribili: Luci di fondo, Restringimenti, Salti di fondo, Gra-
dini, Brusche variazioni di larghezza, cambio di pendenza, cambio di scabrezza, afflusso di portata,
variazione graduale di larghezza, problema di imbocco e problema di deflusso. La trattazione teor-
ica precede la presentazione dell’algoritmo perché è propedeutica per la comprensione dello stesso.
8
In particolare è interessante e particolarmente complesso lo studio effettuato sulla soluzione
del problema di imbocco e del problema di deflusso (che sono poco trattati in letteratura). Per
esempio per entrambi è necessario costruire la scala di deflusso prima di poterli risolvere. Proprio
per questo motivo il codice non è in grado di gestire un problema con più di un deflusso o con un
deflusso e un problema di imbocco.
Segue un’analisi della molteplicità delle soluzioni ammesse allo stato stazionario. Nella prima
parte viene spiegata la teoria sugli ordini di influenza che è stato necessario sviluppare per ottenere
un algoritmo maneggevole (cioè con tempi di calcolo rapidi). Poi si da una spiegazione sul fun-
zionamento dell’algoritmo per la ricerca di soluzioni alternative, e si riporta un esempio esplicativo
di problema con Luci di fondo in successione che ammette più di una soluzione.
Infine si discute il tema della stabilità di tali soluzioni multiple, definendo in primo luogo il
concetto stesso di stabilità, e si dedica l’ultimo capitolo al fenomeno dell’isteresi, cioè al fenomeno
che governa l’instaurarsi di una piuttosto che di un’altra soluzione.
9
2 Modello di calcolo
2.1 Equazioni
Il problema è di trattazione di una corrente monodimensionale in acqua bassa, quindi dove i
terminidivariazionedellegrandezzedelmotosonotrascurabiliindirezionex(rispettoalledirezioni
y e z):
δ
δx
<<
δ
δy
,
δ
δz
(1)
Sotto tali ipotesi le equazioni monodimensionali a fondo fisso (dette anche equazioni di De
Saint Venant) assumono la forma:







δQ
δx
+
δΩ
δt
= 0
δQ
δt
+
δ
δx
(β
Q
2
Ω
2
) +gΩ(
δh
δx
+j) = 0
(2)
Dove:
• β =
R
y
ˆ u(y)ydy
U
2
Ω
e rappresenta quindi il rapporto tra la quantità di moto vera e la quantità di
moto media della sezione. Poiché il profilo di velocità lungo y presenta significative disomo-
geneità (cosa che invece non accade per il profilo verticale logaritmico), β rappresenta una
quantità spesso non assimilabile al valore unitario.
• j =
Q
2
B
C
2
Ω
3
g
rappresenta il termine di dissipazione energetica
• C è la costante di Chezy
• B è la larghezza massima della sezione, in corrispondenza del pelo libero
• h è la quota della superficie libera
Si prende ora il secondo termine dell’equazione del moto. Sotto l’ipotesi che
δβ
δx
’ 0 (fermo
restando che β è diverso da 1)
δ
δx
(β
Q
2
Ω
) =
2Q
Ω
β
δQ
δx
−
Q
2
Ω
2
β(
δΩ
δx
|
h
+b
∂h
∂x
)
L’equazione del moto può quindi essere riscritta:
δQ
δt
+
2Q
Ω
β
δQ
δx
−
Q
2
Ω
2
β(
δΩ
δx
|
h
+b
∂h
∂x
) +gΩ(
δh
δx
+j) = 0
Da che se si raccogli il termine
δh
δx
l’equazione prende il nome di equazione delle correnti:
δQ
δt
+
δh
δx
(gΩ−
Q
2
Ω
2
bβ) +
2Qβ
Ω
δQ
δx
−
Q
2
Ω
2
β
δΩ
δx
|
h
+gΩj = 0 (3)
Si introduce quindi l’ipotesi di moto stazionario: (
δ
δt
= 0). Dall’equazione di continuità si
ricava:
δQ
δx
= 0. Semplificando quindi questo termine più tutti i termini di derivata parziale nel
tempo nell’equazione delle correnti si trova:
δh
δx
(1−
bβQ
2
gΩ
3
)−
βQ
2
gΩ
3
δΩ
δx
|
h
+j = 0
Da cui si può ricavare l’equazione dei profili:
10
δh
δx
=
β
F
2
r
b
δΩ
δx
|
h
−j
1−βF
2
r
(4)
L’equazione dei profili può essere riscritta in termini di tirante Y:
δY
δx
=
δ
δx
(h−Z
f
) =
βF
2
r
b
δΩ
δx
|
y
+i
f
−j
1−βFr
2
(5)
Avendo sostituito nell’ equazione (5):
δΩ
δx
|
h
=
δΩ
δx
|
Y
+bi
f
Alveo rettangolare Si introduce ora l’ipotesi di alveo rettangolare. L’equazione dei profili in
termini di tirante (equazione (5)) si semplifica nell’equazione (6), avendo approssimato β’ 1.
δY
δx
=
F
2
r
b
Y
δb
δx
+i
f
−j
1−Fr
2
(6)
2.1.1 Correnti con variazioni di portata
Sono correnti caratterizzate da una variazione graduale di portata lungo x. Si noti come
la condizione di variazione graduale (ovvero lenta) sia coerente con l’ipotesi di variazioni delle
grandezze del moto lente lungo x (condizione (1)). In generale si tratta di correnti con afflussi e
deflussi. Gli afflussi sono generati da immissioni di affluenti oppure da acqua che precedentemente
derivatadalfiumestesso, vienereintrodottaavalleoppuredacascateverticali. Ideflussisonolegati
a prelievi idrici a scopo idroelettrico oppure ad esondazioni a seguito di alluvioni. La trattazione
fino ad ora svolta non comprende le variazioni graduali di portata. Occorre dunque riscrivere
entrambe le equazioni del sistema (2). Si mostra uno schema del problema in Figura 1
Figura 1: Schema correnti con variazione di portata a seguito di afflussi e deflussi
L’equazione di continuità può essere riscritta come:
∂Ω
∂t
=−
Q
out
−Q
in
dx
=
Q(x +dx) +q
u
(x)dx−q
e
(x)dx−Q(x)
dx
=
∂Q
∂x
−q
u
+q
e
∂Ω
∂x
+
∂Q
∂x
=qe−qu (7)
11
L’equazione della conservazione della quantità di moto vede riscritto un termine:
∂
∂x
(U
2
Ω) =
∂
∂x
(U
2
Ω) +q
u
u
ux
−q
e
u
ex
Si introducono le seguenti ipotesi:
• u
ex
= 0 L’acqua che entra nel fiume ha velocità nulla lungo x. Questa è una buona schema-
tizzazione di afflussi laterali circa perpendicolari al fiume oppure dei "doccioni" ovvero delle
cascate verticali che immettono acqua nel fiume.
• u
ux
=U L’acqua che esce dal fiume ha la stessa velocità in x della corrente. Si può pensare
a un canale defluente circa parallelo al all’asse longitudinale del fiume o a una finestra di
deflusso laterale.
Partendo dall’equazione del moto convenientemente riscritta si ricava l’equazione corrispon-
dente delle correnti, semplificata per moto stazionario:
2
Q
Ω
(qe−qu)
δh
δx
(gΩ−
Q
2
Ω
2
bβ)−
Q
2
Ω
2
β
δΩ
δx
|
h
+gΩj +q
u
U = 0
Da cui si ricava l’equazione dei profili, espressa in termini di tirante:
δY
δx
=
βF
2
r
b
δΩ
δx
|
y
+i
f
−j−
U
gΩ
(2q
e
−q
u
)
1−βFr
2
(8)
Alveo rettangolare L’equazione dei profili, espressa i termini di tirante, sotto l’ipotesi di
alveo rettangolare, in correnti che ammettono variazioni graduali di portata lungo x (e variazioni
di larghezza lungo x) assume dunque la forma:
δY
δx
=
F
2
r
b
Y
δb
δx
+i
f
−j−
Q
gΩ
2
(2qe−qu)
1−Fr
2
(9)
L’equazione (9) è l’equazione implementata dal codice di calcolo per la risoluzione dei moti a
pelo libero in alveo rettangolare.
Risoluzione dell’equazione dei profili











Q =Q(X) (nota a priori)
Y (X
0
) =Y
0
δY
δx
=
F
2
r
b
Y
δb
δx
+i
f
−j−
Q
gΩ
2
(2qe−qu)
1−Fr
2
(10)
Nella condizione in cui la geometria della canaletta rettangolare sia fissata (noto b(x) e
δb
δx
(x)) e
sia nota la pendenzai
f
, la risoluzione dell’equazione dei profili rappresenta un problema di Cauchy.
Questo richiede due condizioni al contorno: la portata deve essere nota a priori su tutto il dominio
e occorre conoscere un valore di tirante da cui "partire" con l’integrazione. Dal punto di vista
matematico un valore di tirante Y (X
o
) ammette soluzioni sia a monte sia a valle di esso (cioè
matematicamente esiste soluzione) ma dal punto di vista fisico le condizioni al contorno vengono
integrate verso monte se sono condizioni di corrente lenta e verso valle se sono di corrente veloce.
12
2.2 Schemi adottati dal modello
Come già detto si tratta di risolvere il problema di Cauchy illustrato nel sistema (10). Comune-
mente per l’integrazione di tale equazione si adotta il metodo di Eulero. Il metodo di Eulero è
una procedura numerica del primo ordine per risolvere equazioni differenziali ordinarie (ODE) una
volta fornito un valore iniziale. Si tratta del più basilare dei metodi espliciti per l’integrazione nu-
merica di equazioni differenziali ordinarie. In equazione (11) si riporta la struttura di un’equazione
(ODE).





d
dx
y(x) =f(y(x))
y(x) =y
o
(11)
2.2.1 Metodo di Eulero
Si supponga che la funzione y sia di classe C
2
. Allora si può scrivere l’espansione in serie di
Taylor:
y(¯ x + Δx) =y(¯ x) +y
0
(¯ x)Δx +y
00
(¯ x)
Δx
2
2
+o(Δx
3
)
avendo indicato con:
•
¯
X un certo punto per cui è nota la soluzione
• y(¯ x) la soluzione in ¯ x
• y
0
la derivata prima della funzione y
• y
00
la derivata seconda della funzione y
• o(Δx
3
) indica una funzione con ordine di grandezza inferiore a Δx
3
Il metodo di Eulero considera l’espansione di Taylor troncata al primo ordine, ovvero:
y(¯ x + Δx) =y(¯ x) +y
0
(¯ x)Δx +o(Δx
2
)
Ciò significa che ad ogni aggiornamento spaziale della grandezza y si commette un errore
proporzionale a Δx
2
. Invece il numero delle iterazioni è inversamente proporzionale a δx: N =
λ/ΔX. Dunque l’errore che si genera nell’integrazione dell’equazione ODE sul dominio, dovuto
alla propagazione dell’errore commesso ad ogni iterazione spaziale, può scriversi:
e = Δx
2
λ
ΔX
=o(ΔX)
Quindi il metodo è lineare, cioè l’errore si propaga linearmente in base alla dimensione dello step
spaziale di discretizzazione del dominio. Questo può rendere il metodo di Eulero dispendioso a
livello computazionale, dovendo in alcuni casi scegliere un passo estremamente piccolo. Il metodo
di Eulero può dare inoltre problemi di stabilità, cioè la soluzione può divergere soprattutto in caso
di equazioni in cui sono coinvolte grandezze con ordini di grandezza molto diversi, per cui può
rendersi necessario scegliere un passo particolarmente piccolo.
Si vuole ora riscrivere il metodo di Eulero applicato ai profili di corrente nella stessa notazione
di equazione (9):
Y (x +dx) =Y (x) +
δY
δx
dx
Dove:
• Y (x) è stata precedentemente calcolata (nota)
•
δY
δx
è calcolato puntualmente con equazione (9)
13
2.2.2 Metodo di Eulero in alveo rettangolare cilindrico a Q=cost
Si analizza a titolo di esempio la situazione in cui si abbia un canale rettangolare cilindrico
(
δb
δx
= 0) con portata costante Q = cost in x. In tale situazione i profili che si realizzano sono i
classici profili M1, M2, M3 in alveo fluviale e S1, S2, S3 in alveo torrentizio. Si analizza ora la
qualità della soluzione ottenuta con il metodo di Eulero per ciascuno di questi profili.
Figura 2: Profili che si realizzano in alveo fluviale cilindrico a Q=cost
Profili M
Profilo M3 Il profilo M3 è un profilo di corrente veloce che si realizza in alveo fluviale
al di sotto dell’altezza critica. Essendo un profilo di corrente veloce (F
2
r
> 1) viene risolto per
integrazione verso valle. In Figura 2 si riporta la forma che assume il profilo M3. Il profilo presenta
un primo tratto di salita con concavità verso il basso seguito da un tratto in salita a pendenza circa
costante per terminare con un ultimo tratto in salita con concavità verso l’alto. Nella prima parte
il valore del coefficiente angolare decresce al crescere del valore del tirante. Dunque, la soluzione
di Eulero, dopo un primo discostamento, tenderà a convergere alla soluzione reale. Infatti, dopo
il primo punto si produce una sovrastima di tirante a cui è associato un coefficiente angolare
maggiore dell’analogo punto del profilo M3 nella stessa posizione x. Dunque, il discostamento nel
successivostepvieneridotto(sesiassumeunintervallodidiscretizzazionesufficientementepiccolo).
Nella parte più o meno rettilinea il profilo M3 e la soluzione di Eulero praticamente coincidono e
infine nella parte terminale c’è un discostamento, poiché il valore del coefficiente angolare cresce
al crescere del valore del tirante. Si riporta una schematizzazione di questo concetto in Figura 3
Figura 3: differenza tra profilo M3 e soluzione di Eulero
14
Profilo M2 Il profilo M2 è un profilo di corrente lenta che si realizza al di sopra dell’altezza
critica e al di sotto del tirante di moto uniforme in alveo fluviale. É un profilo dunque di corrente
lenta e viene quindi integrato verso monte. Si può formalizzare il concetto espresso per il profilo
M3 dicendo che in generale se il valore assoluto del coefficiente angolare cresce con il procedere
dell’integrazione la soluzione di Eulero non è una buona approssimazione del profilo. Integrando
il profilo M2 verso monte, tuttavia, tale situazione non si realizza e quindi il metodo di Eulero è
un buon metodo per questo profilo.
Profilo M1 Il profilo M1 è un profilo di corrente lenta che si realizza al di sopra del tirante
di moto uniforme. Essendo un profilo di corrente lenta viene integrato verso monte. Non presenta
particolari concavità, in particolare il coefficiente angolare si mantiene praticamente costante e
pari a:
δy
δx
’ i
f
salvo in prossimità del tirante di moto uniforme, dove comunque non si hanno
variazioni rapide. Per questo motivo il metodo di Eulero è un buon metodo per il calcolo di questo
profilo
Figura 4: Profili che si realizzano in alveo torrentizio cilindrico a Q=cost
Profili S
Profilo S3 Il profilo S3 è un profilo di corrente veloce che si realizza al di sotto del tirante di
moto uniforme in alveo torrentizio. Per questo motivo viene integrato verso valle. Poiché il modulo
del coefficiente angolare decresce nella direzione di integrazione, facendo riferimento allo studio
fatto per i profili M, si può dire che Eulero è un buon metodo, realizzando una discretizzazione
sufficientemente piccola.
Profilo S2 Il profilo S3 è un profilo di corrente veloce che si realizza al di sotto dell’altezza
critica in alveo torrentizio. Per questo motivo viene integrato verso valle. Poiché il modulo del
coefficiente angolare decresce nella direzione di integrazione, facendo riferimento allo studio fatto
per i profili M, si può dire che Eulero è un buon metodo, realizzando una discretizzazione suffi-
cientemente piccola.
Profilo S1 Il profilo S3 è un profilo di corrente lenta che si realizza al di sopra dell’altezza
critica in alveo torrentizio. Per questo motivo viene integrato verso monte. Poiché il modulo del
coefficiente angolare cresce nella direzione di integrazione, facendo riferimento allo studio fatto
per i profili M, si può dire che Eulero non è il metodo migliore per determinare la forma del
profilo, soprattutto in prossimità della critica dove si assiste a un rapido incremento del modulo
del coefficiente angolare.
15