UNIVERSITÀ DEL SALENTO

FACOLTÀ DI SCIENZE MATEMATICHE, FISICHE E NATURALI

Corso di Laurea in Fisica

L E Z I O N I S U L L A D I N A M I C A D I AT M O S F E R A E O C E A N I

Piero LionelloDipartimento di Scienze dei Materiali

Anno Accademico 2010/2011

I N D I C E

I E L E M E N T I D I F L U I D O D I N A M I C A 1

Introduzione 3

1 C I N E M AT I C A 51.1 Coordinate ed operatori vettoriali 51.2 Le coordinate geografiche 71.3 Coordinate Lagrangiane ed Euleriane 8

1.3.1 Esempi 81.4 Variazioni di volumi e lunghezze 91.5 Funzione di flusso e linee di flusso 101.6 Divergenza, deformazioni e vorticità 11

2 D I N A M I C A 152.1 Le forze nei fluidi 15

2.1.1 Il tensore degli sforzi 152.1.2 La pressione 17

2.2 L’equazione del moto 182.3 Sistemi di riferimento non inerziali 192.4 L’equazione di Navier-Stokes 202.5 La vorticità 212.6 L’integrale di Bernoulli 222.7 Il teorema di Kelvin 232.8 La vorticità potenziale 24

3 M O T I A G R A N D E S C A L A 273.1 Le equazioni in coordinate geografiche 273.2 Sistemi di riferimento cartesiani 283.3 Moti quasi-piani 293.4 Il bilancio geostrofico 313.5 Le coordinate isobariche 313.6 Il vento termico 323.7 Modello a due strati 33

3.7.1 Il modo baroclino 35

4 O N D E L I N E A R I 374.1 Onde di gravità 384.2 Onde di Poincarè–Sverdrup 404.3 Onde di Kelvin 414.4 Onde di Rossby 42

4.4.1 Dinamica dell’onda di Rossby 43

II M O D E L L I S T I C A N U M E R I C A 45

5 I N T E G R A Z I O N E N U M E R I C A D I P D E 475.1 Premessa: metodo alle differenze finite 47

i

Indice

5.2 Onde in una dimensione 485.2.1 Schemi di integrazione espliciti 495.2.2 Dispersione numerica 51

5.3 Metodi impliciti 515.4 Equazione di diffusione 52

5.4.1 Schema FTCS 525.4.2 Schemi impliciti di integrazione 525.4.3 Onde in presenza di diffusione 53

5.5 Staggering e griglie multidimensionali 535.5.1 Equazioni shallow-water monodimensionali 545.5.2 Onde in griglie multidimensionali 545.5.3 Diffusione in griglie multidimensionali 555.5.4 Equazioni di shallow-water, forza di Coriolis e griglie di Araka-

wa 555.6 Boundary value problems 565.7 Metodo spettrale 585.8 Metodo agli elementi finiti 58

ii

Parte I

E L E M E N T I D I F L U I D O D I N A M I C A

I N T R O D U Z I O N E

Vi sono alcuni comportamenti che caratterizzano i fluidi in contrapposizione ai solidi.Mentre i solidi hanno una forma propria, i fluidi si adattano senza opporre resistenzaa variazioni della forma del recipiente in cui sono contenuti, ma è necessario esercitareuna forza per variare il volume da essi occupato. È sempre possibile mescolare un fluido,ossia modificare la posizione di alcuni suoi elementi rispetto agli altri inducendovi unmoto continuo di scorrimenti e vortici, che diminuisce con la viscosità del fluido, manon si annulla per quanto piccola sia la forza applicata. Quando viene esercitata unaforza che tende a deformare un solido, esso oppone una resistenza proporzionale alladeformazione che pertanto aumenta fino al valore per cui la resistenza eguaglia la forzaapplicata, ed una condizione di equilibrio viene conseguita. Quando viene esercitatauna forza che induce uno scorrimento di uno strato di fluido rispetto ad un altro, unfluido oppone una resistenza che è proporzionale alla velocità dello scorrimento stessoche, quindi, continua indefinitamente col valore della velocità a cui la resistenza egua-glia la forza applicata. Pertanto non esiste un meccanismo in grado di interrompere ilmoto prodotto dalla forza, per quanto piccola essa sia, e non viene raggiunto l’equilibrio.Al contrario, se sottoposto ad una forza che tende a variarne il volume, ad esempio unacompressione isotropa, un fluido reagisce con una forza proporzionale alla variazionestessa e viene raggiunta una condizione di equilibrio dove la differenza rispetto al volu-me iniziale è proporzionale alla forza applicata. In conclusione sono fluidi le sostanzeche non possono sostenere senza deformarsi una sollecitazione che ne lasci inalterato ilvolume. Questo criterio non è privo di ambiguità, ed esistono sostanze, quali il vetroed il catrame, che hanno comportamenti duplici, e, a seconda dell’intensità degli sforzia cui sono sottoposte e delle scale dei tempi considerati, si comportano come solidi ofluidi.

Durante il moto, mentre in un solido, in assenza di rotture, le variazioni delle distanzegeometriche tra i vari elementi sono piccole, in un fluido esse sono enormi. Ne segueche la cinematica dei fluidi è molto più complessa di quella dei solidi. I gradi di libertà diun solido, nel caso limite che si comporti come un corpo rigido, si possono ridurre a 6, esono comunque in numero finito, mentre quelli di un fluido sono formalmente infinitied il suo moto viene descritto in modo completo solo specificando punto per punto levelocità locali di espansione–contrazione e rotazione.

Le diverse caratteristiche di fluidi e solidi conseguono da differenze strutturali: i soli-di hanno una struttura cristallina, per cui le loro molecole rispondono in modo unitarioalle sollecitazioni, mentre i legami fra le molecole dei fluidi sono deboli e non perma-nenti (nei liquidi) o assenti (nei gas). Tipicamente il potenziale fra due molecole biato-miche ha un minimo ad una distanza d0 di circa 10−10 m = 1Å. A distanze inferiori laforza è repulsiva, se le molecole non si combinano chimicamente, per effetto della for-za elettrica–quantistica esercitata fra i nuclei, mentre a distanze superiori è attrattiva pereffetto delle fluttuazioni statistiche delle distribuzioni di carica fra gli atomi. Nei solidi ladistanza di separazione tra le molecole è d0, mentre nei liquidi è di poco superiore e nei

3

Indice

gas è nell’ordine di 10d0 nelle abituali condizioni di pressione e temperatura. Le forzeintermolecolari sono quindi molto intense nei solidi, dove l’agitazione termica delle mo-lecole non consente loro di spostarsi in modo apprezzabile dalla posizione di equilibrioattorno alla quale oscillano. Nei liquidi le forze sono più deboli, non vi sono permanen-ti gruppi di molecole parzialmente legate fra loro, e le distanze intermolecolari varianonell’ordine di d0 per effetto dell’agitazione termica. Nei gas le forze intermolecolari agi-scono solo durante le collisioni fra molecole, che altrimenti si muovono liberamente edin modo disordinato. La libertà di movimento delle molecole, che consegue dall’assenzadi una struttura cristallina, rigida e ordinata, è la caratteristica microscopica (comune agas e liquidi) che contraddistingue i fluidi.

L’obiettivo della fluidodinamica è la descrizione del moto relativo di porzioni od ele-menti di fluido sotto l’azione delle forze presenti, attuata ignorando completamente lastruttura corpuscolare della materia, ed, invece, introducendo delle variabili macrosco-piche quali velocità e densità come funzioni continue dello spazio e del tempo. Non esisteuna “prova” che questa, introdotta come ipotesi del continuo, sia una descrizione “cor-retta” della realtà. L’ipotesi del continuo è in accordo con la nostra percezione sensorialedella realtà ed è operativamente giustificata dal suo successo nel descrivere il moto deifluidi. Essa è tuttavia falsa poiché la struttura di tutta la materia, e quindi anche deifluidi, è corpuscolare. La massa è concentrata nei nuclei degli atomi e le molecole simuovono disordinatamente per effetto dell’agitazione termica. Dal punto di vista mi-croscopico, sia densità che velocità non sono affatto funzioni regolari e continue. L’ipo-tesi del continuo è valida per il moto di un fluido se esso coinvolge un intervallo di scalespaziali che verificano simultaneamente due proprietà:

• sono sufficientemente grandi da contenere un numero di molecole sufficiente arendere piccole le varianze dei valori medi di energia, quantità di moto, e massatotale;

• sono sufficientemente piccole da poter essere considerate locali dal punto di vistamacroscopico.

Questi requisiti sono soddisfatti da una vastissima classe di moti. Un esempio è fornitodalla densità. La densità di un fluido dipende dal numero di molecole contenute nel vo-lume considerato. Se il volume è estremamente piccolo, le molecole in esso contenutesono poche e le fluttuazioni del loro numero dovute al moto sono significative nel calco-lo della densità. Se il volume aumenta le fluttuazioni smettono di essere significative e lavalutazione della densità non dipende dal volume. Se il volume aumenta ulteriormentevengono inglobate parti del fluido che hanno proprietà diverse e la densità diventa dinuovo dipendente dalla scelta del volume. Il volume richiesto per contenere un enormenumero di molecole è sorprendentemente piccolo. Si consideri un gas: in un cubo di

lato D sono contenute N =(

D10d0

)3molecole. Quindi, in un cubo di 10−5 m di lato sono

contenute 1012 molecole, numero che spiega la validità dell’ipotesi del continuo.

4

1C I N E M AT I C A

1.1 C O O R D I N AT E E D O P E R AT O R I V E T T O R I A L I

La posizione di ogni elemento di un fluido viene specificata in relazione ad un prefissatosistema di coordinate. Questo, in generale, viene definito introducendo tre famiglie disuperfici nello spazio tridimensionale. A ciascuna famiglia è associata una coordinata ilcui valore denota una superficie, in modo che il vettore xc = (xc

1, xc2, xc

3), selezionando tresuperfici, una per ciascuna famiglia, identifichi il punto dello spazio tridimensionale chene costituisce l’intersezione. Quando le superfici sono piane e mutuamente ortogonalile coordinate x = (x1, x2, x3) si dicono cartesiane ortogonali, ma in generale è sufficientedefinire una tripletta di funzioni xc

1(x1, x2, x3), xc2(x1, x2, x3), xc

3(x1, x2, x3), ciascuna dellequali rappresenta una superficie nello spazio tridimensionale, il cui jacobiano

J(x,xc)= ∣∣∣∣∣ ∂xi

∂xcj

∣∣∣∣∣ (1.1)

sia non nullo, J (x,xc ) 6= 0, perché si possa definire un sistema di coordinate curvilinee eduna trasformazione biunivoca dalle coordinate cartesiane a quelle curvilinee.

Se i versori, la cui direzione varia nello spazio, che identificano localmente le direzionidel sistema di coordinate curvilinee sono mutuamente ortogonali, le coordinate vengo-no dette curvilinee ortogonali. Le coordinate non hanno necessariamente le dimensionidi lunghezze fisiche, le quali sono ad esse legate da un fattore metrico hi , ossia

dsi = hi dxci , h2

i =3∑

j=1

(∂x j

∂xci

)2

, (1.2)

dove dsi è la lunghezza infinitesima di un arco ottenuto variando di dxci la coordinata

xci e lasciando invariate le altre due. Il volume dV = dx1dx2dx3, in coordinate curvilinee

ortogonali diventa

dV =3∏

i=1hi dxc

i . (1.3)

Nel formulare le leggi della fisica vengono utilizzati tre operatori: gradiente, diver-genza e rotore. L’operatore gradiente applicato ad uno scalare produce un vettore cheverifica la proprietà:

lim∆x→0

grad(θ) ·∆x

∆θ= 1, (1.4)

dove ∆x rappresenta uno spostamento infinitesimo. L’operatore divergenza applicatoad un vettore produce uno scalare tale che

div(u) = lim∆V →0

dΦ(u)

∆V, (1.5)

5

C I N E M AT I C A

dove dΦ(u), rappresenta il flusso di u attraverso la superficie dell’elemento infinitesimodi volume ∆V . L’operatore rotore applicato ad un vettore produce un vettore tale che

rot(u) · S = limdS→0

dC (u)

dΦ(u), (1.6)

dove dC (u) rappresenta la circuitazione di u lungo il contorno della superficie infinite-sima dS con versore normale S e dΦ(u) il flusso di u attraverso di essa. Il div(grad(θ))definisce un altro operatore di uso comune chiamato laplaciano.

È consuetudine rappresentare tali operatori mediante il simbolo ∇:

gradθ ≡∇θ, divu ≡∇·u, rotu ≡∇×u, divgradθ ≡∇2θ;

questa notazione non è ambigua ed è perfettamente accettabile, ma non deve nasconde-re il fatto che si tratta di tre operatori ben distinti la cui espressione dipende dal sistemadi coordinate utilizzato. In coordinate ortogonali curvilinee si ha:

∇θ =3∑

i=1

1

hi

∂θ

∂xci

xci ,

∇·u = 1

h1h2h3

[∂u1h2h3

∂xc1

+ ∂h1u2h3

∂xc2

+ ∂h1h2u3

∂xc3

],

∇×u =

∣∣∣∣∣∣∣h1xc

1 h2xc2 h3xc

3∂xc

1∂xc

2∂xc

3

h1u1 h2u2 h3u3

∣∣∣∣∣∣∣ ,

∇2θ = 1

h1h2h3

[∂

∂xc1

(h2h3

h1

∂θ

∂xc1

)+ ∂

∂xc2

(h1h3

h2

∂θ

∂xc2

)+ ∂

∂xc3

(h1h2

h3

∂θ

∂xc3

)].

(1.7)

Solamente per sistemi cartesiani ortogonali, per cui ovviamente h1 = h2 = h3 = 1, si puòformalmente identificare ∇ con il vettore

∇=(∂

∂x1,∂

∂x2,∂

∂x3

)ed eseguire con esso le operazioni vettoriali mantenendo sempre l’analogia formale.

Alcune utili identità vettoriali si possono derivare dalle proprietà degli operatori:

∇×∇θ = 0, (1.8a)

∇·∇×u = 0, (1.8b)

∇·θu = ∇θ ·u+θ∇·u, (1.8c)

∇·u×v = v ·∇×u−u ·∇×v, (1.8d)

∇×θu = ∇θ×u+θ∇×u, (1.8e)

∇× (u×v) = u(∇·v)−v(∇·u)+ (v ·∇)u− (u ·∇)v. (1.8f)

Esse saranno richiamate, assieme alle proprietà degli operatori gradiente, divergenza erotore, nella derivazione delle equazioni in fluidodinamica.

In forma integrale, le proprietà dell’eq. (1.5) e dell’eq. (1.6) si traducono nei teoremidi Gauss e di Stokes, rispettivamente:∫

V∇·udV =

∫S

u ·dS, (1.9)

6

1.2 L E C O O R D I N AT E G E O G R A F I C H E

(dove il primo membro rappresenta un integrale di volume ed il secondo un integralesulla superficie S che racchiude V , ossia il flusso di u attraverso S), e∮

γu ·dl =

∫S∇×u ·dS, (1.10)

(dove il primo membro rappresenta la circuitazione di u lungo la curva chiusa γ ed ilsecondo membro rappresenta un integrale su una generica superficie S limitata dallacurva γ).

1.2 L E C O O R D I N AT E G E O G R A F I C H E

Figura 1: Coordinate geo-grafiche.

Esistono molte scelte possibili di coordinate curvilinee, edesse possono, sia pure con un procedimento non banale,essere opportunamente definite in modo da adattarsi allamorfologia del dominio di integrazione delle equazioni edavere, nella loro versione discreta, risoluzioni spaziali diver-se in diverse parti del dominio. Un sistema di coordinateprivilegiato per la descrizione dei moti su scala planetariasono le coordinate geografiche. Esse sono quasi equivalentialle coordinate sferiche di comune uso in molti argomentidi fisica, ma ne differiscono perché la latitudine è definitaa partire dal piano equatoriale e non dall’asse polare. La lo-ro importanza in geofisica necessita di una loro trattazionedettagliata.

Le coordinate geografiche sono (λ,ϕ,r ), longitudine,latitudine e raggio rispettivamente, e sono legate allecoordinate cartesiane dalle relazioni:

λ= arctan

(x2

x1

), ϕ= arctan

x3√x2

1 +x22

, r =√

x21 +x2

2 +x23 , (1.11)

e

x1 = r cosϕcosλ, x2 = r cosϕsinλ, x3 = r sinϕ. (1.12)

I fattori metrici sono h1 = r cosϕ, h2 = r e h3 = 1, in modo che dV = r 2 cosϕdr dλdϕ.In coordinate geografiche gli operatori gradiente, divergenza, rotore e laplaciano as-

sumono la forma seguente:

∇θ = 1

r cosϕ

∂θ

∂λλ+ 1

r

∂θ

∂ϕϕ+ ∂θ

∂rr

∇·u = 1

r cosϕr

[∂

∂λ(r u1)+ ∂

∂ϕ

(r cosϕu2

)+ ∂

∂r

(r cosϕr u3

)]

∇×u =

∣∣∣∣∣∣∣r cosϕλ r ϕ r

∂λ ∂ϕ ∂r

r cosϕu1 r u2 u3

∣∣∣∣∣∣∣∇2θ = 1

r 2 cos2ϕ

∂2θ

∂λ2 + 1

r 2 cosϕ

∂

∂ϕ

(cosϕ

∂θ

∂ϕ

)+ 1

r 2

∂

∂r

(r 2 ∂θ

∂r

).

7

C I N E M AT I C A

1.3 C O O R D I N AT E L A G R A N G I A N E E D E U L E R I A N E

Vi sono due descrizioni del moto di un fluido, una lagrangiana ed una euleriana. Il pun-to di vista euleriano descrive il moto del fluido con variabili θ(x, t ) che restituiscono ladistribuzione spaziale ad ogni istante del moto, ossia il valore di θ nel punto x all’istantet . Il punto di vista lagrangiano descrive il moto del fluido seguendone ogni elemento.Le variabili θ(x0, t0) descrivono il valore di θ dell’elemento di fluido che ad un istan-te di riferimento t0 si trovava nella posizione x0. Le leggi della fisica, che ad esempiodescrivono la variazione di quantità di moto, temperatura, densità di un elemento difluido, si esprimono naturalmente nella descrizione lagrangiana. Il punto di vista eu-leriano fornisce invece una più immediata rappresentazione del moto del fluido, ed èquello utilizzato più spesso per presentare osservazioni e costruire modelli numerici deifluidi ottenendo, ad esempio, mappe meteorologiche, immagini satellitari, mappe chemostrano l’elevazione superficiale della superficie marina o correnti oceaniche.

La derivata lagrangiana DθDt , detta anche derivata comovente o derivata totale, espri-

me la variazione di θ per un dato elemento di fluido. L’elemento di fluido che occupavaal tempo t la posizione x con velocità u(x, t ) si trova all’istante t +∆t nella posizionex+∆x, dove ∆x = u(x, t )∆t . Quindi

∆θ

∆t= θ(x+∆x, t +∆t )−θ(x, t )

∆t

= θ(x+∆x, t +∆t )−θ(x+∆x, t )+θ(x+∆x, t )−θ(x, t )

∆t,

da cui, nel limite ∆t → 0,

Dθ

Dt= u ·∇θ+ ∂θ

∂t, (1.13)

dove è stata usata la proprietà del gradiente espressa dall’eq. (1.4). La relazione (1.13)raccorda il punto di vista lagrangiano, Dθ

Dt , con il punto di vista euleriano, che osserval’evoluzione locale ∂θ

∂t e le variazioni spaziali ∇θ della stessa grandezza θ.

1.3.1 Esempi

Esempi di descrizioni lagrangiane ed euleriane possono essere identificate nelle tecni-che osservative. Le tecniche che producono naturalmente dati in forma euleriana sono,forse, più comuni, perché si basano sulla raccolta di dati forniti da un insieme di stru-menti di rilevamento la cui collocazione è fissa nello spazio (come sensori collocati supiattaforme oceanografiche od in stazioni meteorologiche) o dalla realizzazione di undispositivo che determini la distribuzione spaziale di un osservabile. Gli esempi piùdiffusi sono gli strumenti fotografici o altri strumenti per telerilevamento, come radare radiometri. La fig. 2 mostra il campo di vento sopra la superficie del mare derivatodai dati rilevati dallo scatterometro montato sul satellite ERS-2. Lo strumento forniscela visione complessiva ed istantanea del campo di vento consentendo la visualizzazio-ne di un fronte che attraversa il Mare Tirreno e delle deformazioni del campo di ventodeterminate dalla presenza delle isole.

8

1.4 VA R I A Z I O N I D I V O L U M I E L U N G H E Z Z E

Figura 2: Campo di vento alle ore 9:24del 18 novembre 1996 fornito dalloscatterometro installato su ERS-2. Lalunghezza delle frecce è proporzionaleal modulo della velocità del vento.

Misure di tipo lagrangiano richiedono l’utilizzodi un dispositivo capace di muoversi assieme alfluido, che, in applicazioni oceanografiche può es-sere costituito da apparati di misura la cui densitàmedia è stata portata a coincidere con quella delfluido da analizzare. La figura 3 mostra il motooscillatorio dovuto ad un’onda a bassa frequenza,osservata nell’oceano Atlantico. Ogni figura del-la sequenza mostra la traiettoria seguita dagli ele-menti del fluido durante una settimana (i dati sonostati mediati per eliminare oscillazioni di frequen-za superiore che sono sovrapposte all’onda a bassafrequenza).

1.4 VA R I A Z I O N I D I V O L U M I E L U N G H E Z Z E I N

U N F L U I D O I N M O T O

Le relazioni che descrivono le variazioni nel tempodel volume e della lunghezza di un elemento fluido

consentono di esprimere in forma differenziale le leggi fisiche. Queste relazioni vengo-no chiamate integrali materiali per sottolineare che si intendono applicate a porzionifissate di fluido, che si deformano e spostano durante il moto.

Figura 3: Traiettorie seguite in un periodo dicirca 3 mesi da un insieme di boe rilasciate ori-ginariamente nella Corrente del Golfo, che con-sente di visualizzare il moto della corrente e lesue fluttuazioni.

Ogni elemento δS della superficie chedelimita un elemento di volume δV sisposta durante il moto con la velocità lo-cale u del fluido, determinando una varia-zione di volume u ·δS. Sommando su tut-ta la superficie dell’elemento di volume siottiene la variazione di un elemento di vo-lume δV in termini di δΦ(u), flusso di u:

DδV

Dt= δΦ(u) = δV (∇·u), (1.14)

dove è stata usata la proprietà espressadall’eq. (1.5) dell’operatore divergenza ∇·u = limδV →0

δΦ(u)δV .

I due estremi di un elemento lineare δlsi muovono con diverse velocità che diffe-riscono per la quantità δu = (δl ·∇)u, che

rappresenta, quindi, la velocità di variazione di δl:

Dδl

Dt= (δl ·∇)u. (1.15)

9

C I N E M AT I C A

Usando l’eq. (1.14), il tasso di variazione della massa di un elemento di fluido, DδMDt ,

risulta dal tasso di variazione della densità % e dalla divergenza del campo di velocità:

DδM

Dt= D(%δV )

Dt=

(D%

Dt+%∇·u

)δV . (1.16)

In assenza di aumento o perdita di massa, la conservazione della massa, DδMDt = 0, si

esprime in forma differenziale,

D%

Dt+%∇·u = 0, (1.17)

equazione denominata equazione di continuità: essa esprime il bilancio di massa perun elemento di volume che si muove con il fluido, secondo cui la variazione percentualedi densità di un elemento di fluido è bilanciata dalla variazione percentuale di volume(punto di vista lagrangiano). Se il fluido è incompressibile, allora D%

Dt = 0 e quindi

∇·u = 0, (1.18)

ossia il campo di velocità è indivergente. In presenza di una sorgente di massa SM perunità di volume si avrebbe DδM

Dt = SMδV e quindi D%Dt +%∇·u = SM .

In base all’identità vettoriale (1.8c), l’eq. (1.17) si può riscrivere come

∂%

∂t+∇· (%u) = 0, (1.19)

che esprime il bilancio per un elemento di volume fisso nello spazio fra la variazionelocale di densità e la divergenza del flusso di densità %u (punto di vista euleriano).

Le relazioni (1.13), (1.14), (1.15), (1.16) non esprimono leggi fisiche, ma conseguonodalla sola descrizione cinematica dei fluidi, e consentono di formulare dei tassi di va-riazione lagrangiani in forma euleriana. Occorrono la fisica e la dinamica del fluido perpoter specificare i flussi non avvettivi e le funzioni sorgenti e, quindi, il valore dei tassidi variazione lagrangiani, in modo da ottenere la forma differenziale delle leggi fisicheadatta allo studio del moto dei fluidi.

1.5 F U N Z I O N E D I F L U S S O E L I N E E D I F L U S S O

La linea di flusso è definita come una linea che è in ogni suo punto tangente alla velocitàu = (u, v , w) del fluido:

δl ·u = δl |u|⇔ δxu = δy

v = δzw , (1.20)

dove δl = (δx,δy ,δz) è un vettore infinitesimo localmente orientato come la linea diflusso. Un tubo di flusso è costituito dall’insieme delle linee di flusso che ad un datoistante passano per il contorno di una superficie fissata. Se il moto è stazionario, sia tubiche linee di flusso sono costanti ed il percorso seguito da un elemento di fluido, dettotraiettoria, coincide con una linea di flusso.

La condizione espressa dall’eq. (1.18) è automaticamente soddisfatta se si definisceun vettoreΨ, tale che

u =∇×Ψ

10

1.6 D I V E R G E N Z A , D E F O R M A Z I O N I E V O R T I C I TÀ

perché si ha che ∇ · (∇×Ψ) ≡ 0. Dalla definizione segue immediatamente che il flussodi u attraverso una superficie coincide con la circuitazione di Ψ lungo la curva che ladelimita.

Se il moto è bidimensionale, ossia u = (u, v ,0), allora Ψ = (0,0,ψ) e la funzione ψ

è detta funzione di flusso. Le isolinee di ψ sono linee di flusso, ossia sono tangenti allavelocità e la loro densità indica la velocità del moto, essa consente di visualizzare il motodel fluido. Valgono, infatti, le due proprietà seguenti.

1. Le isolinee di ψ sono tangenti alla velocità del fluido. Una curva lungo cui ψ ècostante soddisfa:

0 = δψ= ∂ψ

∂xδx + ∂ψ

∂yδy =−vδx +uδy

e quindi la condizione (1.20) è verificata.

2. Il flusso di della velocità fra due isolinee A e B diψ è costante e valeψ(B)−ψ(A). Ladimostrazione è immediata1:∫ B

A

(udy − vdx

)= ∫ B

A

(∂ψ

∂xdx + ∂ψ

∂ydy

)=ψ(B)−ψ(A).

Come le linee di forza di un campo visualizzano la direzione e l’intensità della forza,così le linee di flusso visualizzano l’intensità e la direzione della velocità del fluido. Ilgrafico che risulta è costante solo se il moto è stazionario, nel qual caso ψ non dipendedal tempo. La fig. 4 rappresenta la quota a cui la pressione atmosferica vale 500hPa,corrispondente approssimativamente al livello che divide verticalmente l’atmosfera indue parti di massa eguale. Nel regime dinamico caratteristico dell’atmosfera tale quotaè proporzionale alla funzione di flusso e le sue isolinee rappresentano la direzione el’intensità della circolazione atmosferica.

1.6 D I V E R G E N Z A , D E F O R M A Z I O N I E V O R T I C I TÀ

Localmente il moto di un fluido è costituito dalla somma di tre componenti con carat-teristiche diverse: una rotazione, un’espansione/contrazione isotropa, una deformazioneisocora senza deformazioni di taglio. Il campo di velocità del fluido può sempre esseredecomposto in queste tre componenti, che possono essere presenti sia contemporanea-mente che separatamente. Un moto viene detto irrotazionale se la rotazione è assente,incompressibile se l’espansione/contrazione isotropa è assente.

L’analisi del campo di velocità di un fluido si ottiene da un’espansione locale dellavelocità in funzione dello spazio. Al primo ordine, la velocità relativa∆u di due elementidi fluido separati da una distanza ∆x risulta2

∆ui = ui (x+∆x)−ui (x) ≈ ∂ui

∂x j∆x j , (1.21)

1 Occorre tenere presente che, nel caso bidimensionale, il flusso è attraverso curve; in particolare, dato il seg-mento orientato

(x′−x, y − y ′

), il flusso si ottiene moltiplicando scalarmente u per il vettore normale a tale

segmento e con lunghezza data dalla lunghezza del segmento, fornito dalla relazione n = (y ′− y ,−x′+x

).

2 È sottintesa la somma sugli indici ripetuti.

11

C I N E M AT I C A

Figura 4: Linee di flusso della circolazione atmosferica. La figura mostra il la quota acui la pressione vale 500hPa, detto geopotenziale a 500hPa, che rappresenta una adeguatarappresentazione della circolazione atmosferica media nella troposfera

che si decompone in una componente simmetrica ∆usi ed una antisimmetrica ∆ua

i :

∆ui =∆usi +∆ua

i = e si j∆x j +ea

i j∆x j , (1.22)

dove

e si j =

1

2

(∂ui

∂x j+ ∂u j

∂xi

), ea

i j =1

2

(∂ui

∂x j− ∂u j

∂xi

). (1.23)

La componente simmetrica è

∆usi =∆x j e s

i j =1

2

∂

∂∆xi

(∆xk∆x j e s

j k

)= ∂φ

∂∆xi, (1.24)

dove φ rappresenta una forma quadratica nello spazio 3-dimensionale, che, con unaopportuna rotazione del sistema di riferimento, si può sempre riscrivere come

φ(∆x) = 1

2

(∆xi e s

i i

)2 . (1.25)

Il contributo ∆usi = e s

i i∆xi rappresenta un moto in cui un elemento di fluido viene al-lungato, se

∑i e s

i i > 0, o schiacciato, se∑

i e si i < 0, nella direzione degli assi principali

dell’ellissoide (1.25). Nel caso di un moto incompressibile la divergenza si annulla, ossia∑i e s

i i = 0, ed allungamenti e schiacciamenti si compensano in modo da lasciare inalte-rato il volume dell’elemento. Pertanto, posto∆=∑

i e si i , la velocità∆us

i si divide ulterior-mente in una contrazione (espansione) isotropica ∆ue , caratterizzata dalla divergenza

∆, associata ad un tensoreδi j∆

3 , e in uno stiramento (schiacciamento) senza variazioni

di volume ∆ub associato al tensore e si j −

δi j∆

3 .

12

1.6 D I V E R G E N Z A , D E F O R M A Z I O N I E V O R T I C I TÀ

L’espressione del contributo antisimmetrico ∆uai coinvolge il vettore

ωk = εki j∂u j

∂xi= (∇×u)k .

Usando l’identità εi j kεlmk = δi mδ j l −δi lδ j m , tale che eai j = 1

2εi j kωk , l’espressione com-pleta di ∆ua

i diventa

∆uai = 1

2εi j kω j∆xk = 1

2(ω×∆x)i .

Il contributo antisimmetrico rappresenta pertanto una rotazione con velocità angolareω2 . Il vettoreω è chiamato vorticità del fluido.

Al primo ordine in una piccola regione attorno ad x la velocità di un fluido è costituitada tre distinti contributi: una rotazione, una compressione isotropa, una deformazio-ne (allungamento+ compressione) senza variazioni di volume. Indicando le rispettivevelocità come ∆us , ∆ue , ∆ub si ha:

u(x+∆x) = u(x)+∆us +∆ue +∆ub .

Ne segue che il campo di velocità è costituito da tre componenti fondamentali: ue , ub ,us , dove i contributi ue , ub sono entrambi irrotazionali, mentre us è solenoidale. Questadecomposizione non vale solo per il campo delle velocità. Un generico campo vettorialesi può sempre decomporre in questi tre campi costituenti (Teorema di Clebsh).

L’utilità di questa decomposizione non è soltanto l’analisi delle componenti fonda-mentali del moto, ma essa consente la ricostruzione del moto del fluido a partire daquantità dinamiche che possono essere determinate in altro modo. La componenteue viene determinata dalla distribuzione di sorgenti di volume interne al fluido, chepossono essere dovute ad effettive immissioni di fluido o a variazioni del suo stato ter-modinamico (espansioni e/o contrazioni). La componente us viene invece determina-ta dalla distribuzione della vorticità, la cui evoluzione nel tempo è determinata dallacorrispondente equazione (vedasi sez. 2.5).

In alcuni casi sono possibili semplificazioni. In geofisica si assume quasi sempre unmoto incompressibile, ossia ∆ = 0, e ue risulta nullo. Nel limite di moti a grande scalala vorticità del fluido è determinata dalla rotazione terrestre ed è quindi nota a priori.La vorticità caratterizza il moto dei fluidi e la sua presenza consente il moto anche sela superficie libera del fluido è immobile. Un esempio è costituito dalle onde corte diRossby (sez. 4.4). Il moto associato a ub ha invece la caratteristica opposta poiché èunivocamente determinato dal moto della sua superficie libera: esso è importante permolti moti a piccola scala con vorticità trascurabile. Un esempio è il moto delle onde digravità di alta frequenza (4.1).

13

2D I N A M I C A

2.1 L E F O R Z E N E I F L U I D I

Le forze che agiscono sugli elementi di un fluido appartengono a due categorie: forze divolume e forze di superficie.

Le forze di volume sono le forze a lungo raggio che penetrano all’interno degli elementidi fluido e sono proporzionali alla massa dell’elemento di fluido. Se f è la forza per unitàdi massa, allora sull’elemento di volume δV agisce una forza δF data da

δF(v) = %fδV .

Gli esempi più significativi sono la gravità e le forze fittizie in sistemi di riferimento noninerziali.

Le forze di superficie sono forze a corto raggio che, essendo dovute alle collisioni mo-lecolari, coinvolgono solo le molecole superficiali dell’elemento e sono quindi propor-zionali alla sua superficie. Se Σ(n) è la densità superficiale di forza – che dipende dall’o-rientazione della superficie rappresentata dal versore n = (n1,n2,n3) – e δS l’elementodi superficie, allora

δF(s) =Σ(n)δS

rappresenta la forza di superficie esercitata attraverso δS dal fluido verso cui n punta. Lapressione e la viscosità dei fluidi sono forze di superficie.

2.1.1 Il tensore degli sforzi

La densità di forza di superficie Σ soddisfa tre proprietà:

1. Σ è una funzione dispari di n;

2. Σ si può sempre scrivere come

Σi =σi j n j (2.1)

dove σi j è un tensore, detto tensore degli sforzi;

3. il tensore σi j è simmetrico, εi j kσ j k = 0, ossia σ j k =σk j .

Per provare la prima proprietà, si consideri la forza che agisce tra due elementi di fluido,1 e 2, separati da una superficie ∆S, con n versore ad essa normale. Alla forza F12 =Σ(n)∆S esercitata dall’elemento 1, verso cui n punta, su 2 corrisponde allora una forzaeguale e contraria F21 = −F12, esercitata da 2 su 1. Per definizione si ha F21 = Σ(−n)∆S.Quindi Σ(n) =−Σ(−n) .

15

D I N A M I C A

Figura 5: Elemento infinitesimo a forma di tetraedro a cui sono applicate le forze di superficie.

Per dimostrare la seconda proprietà, si consideri un elemento di fluido a forma ditetraedro come in fig. 5. Il bilancio delle forze risulta

%Du

DtδV = f%δV +Σ(n)δS +Σ(−n1)δS1 +Σ(−n2)δS2 +Σ(−n3)δS3.

Poiché

δS j = n j · nδS,

si ottiene(Du

Dt− f

)%δV =

[Σ(n)−

3∑j=1Σ(n j )n j · n

]δS.

Poiché i membri sono infinitesimi di ordine diverso, l’eguaglianza è valida indipenden-temente dalle dimensioni del tetraedro solo se il membro di destra è identicamentenullo:

Σ(n) =3∑

j=1Σ(n j )n j · n.

Il membro di sinistra è un vettore ed il membro di destra è il prodotto scalare fra un vet-tore ed un’altra grandezza, che deve essere un tensore di rango 2 affinché l’eguaglianzasia valida indipendentemente dal sistema di riferimento, ossia

Σi (n j )n j =σi j ,

da cui

Σi =σi j n j .

La terza proprietà si dimostra in modo analogo, considerando il bilancio del momen-to angolare per un elemento di fluido di forma cubica e dimensione L rispetto al suo

16

2.1 L E F O R Z E N E I F L U I D I

baricentro O. Le forze di volume risultano, a meno di infinitesimi di ordine superiore inL, in una forza applicata nel punto O stesso il cui momento angolare è nullo. Le forzedi superficie generano un momento angolare la cui componente verticale z è di ordine(σx y −σy x )L3. Il momento di inerzia del cubo è proporzionale alla massa del cubo ed alquadrato della sua dimensione, cioè ad L5. Pertanto l’accelerazione angolare è finita so-lo se σx y =σy x . Applicando lo stesso ragionamento alle altre componenti del momentoangolare si ottiene che σi j =σ j i per ogni i e j .

Per ogni direzione i esiste quindi un vettoreσi = (σi 1,σi 2,σi 3) la cui componente per-pendicolare alla superficie δS fornisce la densità superficiale della componente i dellaforza: F (s)

i =σi j n jδS.

2.1.2 La pressione

Poiché la traccia di σi j è indipendente dal sistema di riferimento, è possibile definireuno scalare p

p =− tr [σ]

3. (2.2)

Si consideri un sistema di riferimento opportunamente ruotato in modo che il tensoreσi j sia diagonale e l’azione delle forze di superficie si riduca a compressioni od espan-sioni σi j = 0 se i 6= j e quindi lo si decomponga in due parti

σ=

−p 0 00 −p 00 0 −p

+

σ11 +p 0 00 σ22 +p 00 0 σ33 +p

.

Il primo termine rappresenta una forza uniforme che agisce sulla superficie di un ele-mento di fluido determinando una compressione (od espansione) isotropa. La sua azio-ne determinerebbe necessariamente una variazione di volume a cui la struttura del flui-do è in grado di opporsi. Il secondo termine è associato ad almeno una coppia compres-sione + espansione poiché la sua traccia è 0. Esso quindi può deformare l’elemento difluido senza modificarne il volume, e deve quindi annullarsi se il fluido è in equilibrio.Ne consegue la seguente necessaria condizione per l’equilibrio: in un fluido in quiete laforza esercitata su ogni elemento è normale alla sua superficie ed indipendente dal suoorientamento, ossia Σ(n) = −pn. Lo scalare p è chiamato pressione. Esso rappresentail valore medio della forza normale ad una superficie nel fluido per unità di superficie,dove la media è intesa su tutte le possibili orientazioni della superficie.

Condizione necessaria per l’equilibrio del fluido è che la risultante delle forze applicatesia nulla,∫

δVfi% dV −

∫δS

pni ·δS = 0,

che, per il teorema di Gauss, diventa:∫δV

[fi%−∇· (pni )

]dV = 0.

17

D I N A M I C A

Includendo in un’unica equazione vettoriale le tre componenti si ottiene la condizionenecessaria per l’equilibrio:

f%=∇p . (2.3)

Se l’unica forza agente è la forza di gravità g =−g z, allora dall’eq. (2.3) si ottiene

%g =−∂p

∂z. (2.4)

La pressione così definita viene denominata pressione idrostatica. Se la densità è costan-te, ossia %= %0, l’eq. (2.4) si integra esplicitamente :

p = p0 −%0g z.

Tale condizione è generalmente soddisfatta nell’oceano, dove le variazioni di densitàsono nell’ordine dei pochi per mille.

La densità non è costante nell’atmosfera, dove, in condizione di temperatura costante,per la legge di Boyle vale

p

%= R?T = costante = g H , R? = R

m, con m peso molecolare medio.

e, sostituendo nell’eq. (2.3), si ottiene

p = p0 exp(− z

H

)dove H rappresenta la scala verticale dell’atmosfera.

2.2 L’ E Q U A Z I O N E D E L M O T O

La seconda legge di Newton, applicata ad un elemento di fluido, assume la forma

DδP

Dt= δF (2.5)

dove δP è la quantità di moto e δF è la risultante delle forze applicate. Poiché la massa èconservata, il tasso di variazione della quantità di moto è:

DδP

Dt= D

Dtδ(Mu) = δV %

Du

Dt.

Il contributo delle forze di superficie nella direzione ni si ottiene dall’integrazione dellacomponente normale del vettoreδF(s) sulla superficieδSi che delimita l’elemento di flui-do, ossia dal flusso del vettoreσi attraverso δSi . La somma delle forze di volume e delleforze di superficie dà la forza complessiva, le cui componenti sono δFi = %δV fi +Φ(σi ),per cui la legge di Newton (2.5) diventa %Dui

Dt = % fi +∇·σi . Utilizzando la notazione car-tesiana e separando la pressione dal restante componente di j dello stress, l’equazionedel moto risulta:

%Dui

Dt= % fi − ∂p

∂xi+ ∂di j

∂x j. (2.6)

18

2.3 S I S T E M I D I R I F E R I M E N T O N O N I N E R Z I A L I

Il termine a sinistra rappresenta la variazione di momento per unità di volume. Poichéla massa è conservata, la variazione di densità è compensata da una variazione di volu-me e non contribuisce alla variazione di momento. I termini di destra rappresentano larisultante delle forze di volume e delle forze di superficie, il secondo essendo determi-nato dalla differenza fra forze di superficie agenti su superfici opposte dell’elemento difluido.

2.3 S I S T E M I D I R I F E R I M E N T O N O N I N E R Z I A L I

Nella fluidodinamica geofisica il sistema di riferimento naturalmente utilizzato per os-servare e descrivere i moti del fluido è solidale con la Terra e pertanto non inerziale. Inesso appaiono forze fittizie di fondamentale importanza dinamica per i moti su scalaplanetaria.

Si consideri un sistema di riferimento che ruota con velocità angolareΩ, in modo chei suoi assi di riferimento ci varino nel tempo dci

dt = Ω× ci . La variazione nel sistemainerziale di un vettore v =∑

i vi ci , dove vi sono le coordinate nel sistema di riferimentonon inerziale, è data da

dv

dt

∣∣∣∣inerziale

=∑i

(dvi

dtci + vi

dc

dt

)= dv

dt

∣∣∣∣non inerziale

+Ω×v, (2.7)

dove il primo contributo rappresenta la variazione relativa al sistema di riferimento noninerziale e il secondo risulta dalla rotazione degli assi. Applicando questa formula alvettore x che rappresenta la posizione di un elemento di fluido si ottiene

d2x

dt 2

∣∣∣∣in.

= d

dt

(dx

dt

∣∣∣∣non in.

+Ω×x)= d2x

dt 2

∣∣∣∣non in.

+2Ω×Vr + dΩ

dt

∣∣∣∣in.

×x+Ω× (Ω×x),

che fornisce la relazione fra l’accelerazione nel sistema di riferimento rotante e in quelloinerziale. Se il sistema non inerziale subisce inoltre un’accelerazione, il membro di de-stra contiene un ulteriore contributo f0. Il membro di sinistra è l’accelerazione del fluidoche è proporzionale alla risultante delle forze applicate secondo l’equazione del moto.Quindi si conclude che

%dV

dt

∣∣∣∣non inerziale

+2%Ω×Vr +%dΩ

dt×x+%Ω× (Ω×x)+ f0 = %f+∇·σ.

La rotazione determina la presenza di tre forze: la forza di Coriolis 2%Ω×Vr , una forzadeterminata dalle variazioni della velocità angolare, %dΩ

dt ×x e la forza centrifuga %Ω×(Ω×x).

La forza di Coriolis dipende dalla componente della velocità perpendicolare all’assedi rotazione. Se il moto ha luogo nel piano x y , in un sistema di riferimento con velocitàdi rotazione orientata nella direzione z la forza di Coriolis determina una deviazioneverso destra rispetto alla direzione del moto.

La variazione della velocità angolare del sistema di riferimento è normalmente deltutto trascurabile per applicazioni meteorologiche ed ambientali.

19

D I N A M I C A

La forza centrifuga tende ad allontanare il fluido dall’asse di rotazione e si può espri-mere come il gradiente di un opportuno potenziale: essa viene perciò inglobata nell’e-spressione delle forze conservative, e ha come risultato il loro aumento con l’aumentaredella distanza dall’asse di rotazione (vedasi sez. 2.5).

2.4 L’ E Q U A Z I O N E D I N AV I E R - S T O K E S

L’equazione del moto eq. (2.6) non è risolvibile senza specificare l’espressione del tenso-re σi j . L’equazione di Navier-Stokes si ottiene dall’eq. (2.6) esplicitando la forma di σi j

e conseguentemente di di j , in funzione della velocità del fluido.Come visto nella sez. 2.1.2, il tensore σi j si scompone in una componente dovuta alla

pressione ed in un tensore di deviazione di j = σi j + pδi j che è simmetrico, ha traccianulla di i = 0, e si annulla se il fluido è in equilibrio. Si assume che di j sia lineare nellederivate della velocità:

di j = Ai j kl∂uk

∂xl.

dove Ai j kl è un tensore isotropo, ovvero della forma generale Ai j kl =µ1δi jδkl+µ2δi kδ j l+µ3δi lδ j k ,con µ1,µ2,µ3 costanti arbitrarie. Queste assunzioni corrispondono all’eviden-za che gli attriti sono nulli in assenza di gradienti della velocità e che in fluidi come l’ariae l’acqua non esistono direzioni preferenziali. Poiché di j è simmetrico, l’espressione sisemplifica in Ai j kl =µ1δi jδkl +µ2(δi kδ j l +δi lδ j k ). In precedenza è stato dimostrato che

∂uk

∂xl= e s

kl +1

2εklmωm

e quindi si ha

di j = [µ1δi jδkl +µ2(δi kδ j l +δi lδ j k )]

(e s

kl +1

2εklmωm

)= [µ1δi jδkl +2µ2δi kδ j l ]e s

kl ,

cioè

di j =µ1δi j e skk +2µ2e s

i j . (2.8)

Poiché di j ha traccia nulla, 3µ1 +2µ2 = 0 e l’espressione finale per di j è

di j =µ(∂ui

∂x j+ ∂u j

∂xi− 2δi j

3

∂uk

∂xk

). (2.9)

L’equazione del moto risulta

%Dui

Dt= % fi − ∂p

∂xi+ ∂

∂x j

[µ

(∂ui

∂x j+ ∂u j

∂xi− 2δi j

3

∂uk

∂xk

)](2.10)

che è l’equazione di Navier-Stokes. Generalmente µ è una costante che dipende dallanatura del fluido e l’equazione si semplifica in

%Dui

Dt= % fi − ∂p

∂xi+µ

(∂2ui

∂x j∂x j+ 1

3

∂2u j

∂xi∂x j

)(2.11)

20

2.5 L A V O R T I C I TÀ

che, per moti incompressibili si semplifica ulteriormente in

Dui

Dt= fi − 1

%

∂p

∂xi+ν ∂2ui

∂x j∂x j(2.12)

dove ν= µ% viene chiamata viscosità cinematica.

La variazione locale ∂ui∂t = Dui

Dt − u j∂ui∂x j

della velocità è dovuta alle forze esterne, algradiente di pressione, alla viscosità ed alle nonlinearità. Il numero di Reynolds R è ilrapporto fra queste due ultime quantità:

R = nonlinearità

viscosità=

∣∣∣u j∂ui∂x j

∣∣∣∣∣∣ν∂2ui

∂x j2

∣∣∣ ∼ velocità · lunghezza

ν.

Quindi R À 1 caratterizza i regimi non lineari, con presenza di vortici intensi e getti,viceversa R ¿ 1 caratterizza i regimi dove la viscosità prevale, ed il moto del fluido nonsubisce brusche variazioni di velocità.

2.5 L A V O R T I C I TÀ

La vorticità ω = ∇×u del fluido è un’importante quantità dinamica che caratterizza ilmoto di un fluido, e rappresenta il doppio della sua velocità locale di rotazione. In unsistema di riferimento rotante essa risulta dalla somma di una vorticità relativa ω eduna vorticità ambientale 2Ω: ωa = 2Ω+ω. La componente 2Ω è dovuta al moto delsistema di riferimento non inerziale, mentre la componenteω è dovuta al moto relativodel fluido, osservato nel sistema di riferimento inerziale.

Per derivare le equazioni che ne descrivono le variazioni si utilizzano tre identità vet-toriali:

u ·∇u = ∇|u|22

+ω×u (2.13a)

A× (B×C) = (A ·C)B− (A ·B)C (2.13b)

∇× (A×B) = A(∇·B)−B(∇·A)+ (B ·∇)A− (A ·∇)B (2.13c)

L’equazione del moto (2.11), riscritta in forma più compatta ed in un sistema di riferi-mento rotante trascurando il termine (piccolo per le applicazioni geofisiche) di variazio-ne della velocità angolareΩ, risulta

Du

Dt=−∇p

%+ f−2Ω×u−Ω× (Ω×x)+ν

(∇2u+ 1

3∇(∇·u)

). (2.14)

La (2.14) si semplifica ulteriormente, in quanto la forza centrifuga, usando l’eq. (2.13b),si può esprimere mediante un potenzialeΨc :

Ω× (Ω×x) = (Ω ·x)Ω−Ω2x =−Ω2x⊥ =−∇(Ω2x2

⊥2

)≡−∇Ψc ; (2.15)

21

D I N A M I C A

pertanto, se le forze di volume sono conservative, e.g. f =−g z =∇φ, i due potenziali Ψc

eφ si possono inglobare in un unico potenziale totaleΨ=φ+Ψc e, utilizzando anche la(2.13a), l’equazione del momento diventa

∂u

∂t+ω×u+∇|u|2

2=−∇p

%+∇Ψ−2Ω×u+ν

(∇2u+ 1

3∇(∇·u)

). (2.16)

L’equazione per la vorticità si ottiene prendendo il rotore di (2.16):

∂ω

∂t= ∇%×∇p

%2 +ν∇2ω− (2Ω+ω)(∇·u)+ [(2Ω+ω) ·∇]u, (2.17)

dove si è assunto cheΩ sia costante.La variazione della vorticità di un elemento di fluido è dovuta a quattro termini. Il

primo è presente se isolinee di pressione e di densità non sono parallele. In tal caso ilfluido si dice baroclino, altrimenti, se il termine è nullo, il fluido si dice barotropico. Inun fluido barotropico, poiché i gradienti di densità e pressione sono paralleli, l’una puòessere espressa in funzione dell’altra. Questo termine rappresenta le rotazioni indottein elementi di fluido dall’azione dello stesso gradiente di pressione su porzioni di fluidoaventi diverse densità o, viceversa, da un diverso gradiente di pressione su porzioni difluido aventi uguale densità. Il secondo termine rappresenta l’azione della viscosità, chetende ad omogeneizzare la vorticità all’interno del fluido. Per comprendere il significatodegli altri due termini, si consideri un sistema di riferimento il cui asse z sia orientatolungo la direzione locale della vorticità totaleωa = 2Ω+ω:

(ωa ·∇)u−ωa(∇·u) =ωa∂z(ux x+uy y+uz z

)−ωa z(∂x ux +∂y uy +∂z uz

)=ωa∂z

(ux x+uy y

)−ωa z(∂x ux +∂y uy

).

Il secondo termine rappresenta l’intensificazione della vorticità dovuta alla conver-genza del fluido in un piano normale alla direzione locale della vorticità, che conseguedalla conservazione del momento angolare. Se il moto è incompressibile, questo motocorrisponde all’allungamento degli elementi di fluido lungoωa ed è chiamato stiramen-to dei tubi di vorticità. Il primo termine rappresenta l’effetto di variazioni lungoωa dellavelocità perperdicolare aωa .

In assenza di baroclinicità e viscosità, considerando per semplicità un fluido incom-pressibile nel caso bidimensionale,ωa =ωa z, si ha

Dωa

Dt=ωa

∂uz

∂z→ dωa

ωa= duz

uz.

Questo termine esprime la tendenza dei filamenti di vorticità a inclinarsi con i filamentidi fluido in assenza di baroclinicità e di viscosità ed è chiamata inclinazione dei tubi divorticità.

2.6 L’ I N T E G R A L E D I B E R N O U L L I

L’eq. (2.16) consente di derivare un’importante legge di conservazione. Se % è costante,le forze sono conservative, il moto è stazionario (ossia ∂u

∂t = 0), e la viscosità trascurabile(ν= 0), allora consegue che (2Ω+ω)×u =∇H dove

H =Ψ− p

%− |u|2

2,

22

2.7 I L T E O R E M A D I K E LV I N

che implica H costante lungo linee di flusso e lungo isolinee della vorticità totale 2Ω+ω. Se la vorticità totale è nulla, allora H è costante su tutto il fluido. La funzione H èchiamata funzione di Bernoulli.

2.7 I L T E O R E M A D I K E LV I N

La divergenza sia della vorticità totaleωa che di quella relativaω è identicamente nulla.Infatti

ωa =∇×ua =∇× (u+Ω× r) =ω+2Ω

per cui

∇·ωa = 0, ∇·ω= 0.

Per tubo di vorticità si intende la superficie costituita dalle linee di vorticità che passanoper una linea chiusa. Tali linee sono le linee tangenti in ogni loro punto al vettore vor-ticità. Poiché il flusso della vorticità attraverso una superficie chiusa è nullo ed il flussoattraverso le pareti di un tubo di vorticità è nullo per definizione, ne segue che il flussodella vorticità attraverso una sezione del tubo non dipende dalla scelta della sezione. Sidefinisce forza di un tubo

Γ=∫

Sω ·dS =

∮γ

u ·dl

dove γ è una linea chiusa che giace sul tubo di vorticità. Utilizzando l’eq. (1.15) si ottiene

DΓa

Dt=

∮γ

dl · Dua

Dt+

∮γ

ua (dl ·∇)ua

dove il secondo integrale si annulla perché è un differenziale esatto (ua · (dl ·∇)ua = dl ·∇u2

a = du2a). Per il calcolo del primo integrale si utilizza l’eq. (2.14). Le forze conservative

danno un contributo nullo. Il contributo del gradiente di pressione è:

−∫γ

∇p

%·dl =

∫S

∇%×∇p

%2 ·dS

e quindi si annulla se il fluido è barotropico. Si può enunciare quindi il seguente

TEOREMA DI KELVIN 2.1 - Se la viscosità è nulla, se il fluido è barotropico e le forze sonoconservative, l’intensità di un tubo di vorticità totale è costante.

Dal teorema di Kelvin consegue che i tubi di flusso sono permanenti, ossia una su-perficie materiale che individua un tubo di flusso continua a individuarlo durante ilsuccessivo moto del fluido. Una linea di vorticità è il limite cui tende un tubo di flus-so al diminuire della sua sezione. Il teorema di Kelvin mostra che le linee di vorticità simuovono con il fluido, se il fluido è barotropico e non viscoso.

Così come la vorticità, anche la forza di un tubo può essere divisa in una componen-te ambientale ed una relativa: Γa = Γ+2Ω∆A⊥, dove ∆A⊥ rappresenta una sezione del

23

D I N A M I C A

tubo ortogonale adΩ, vorticità planetaria. Poiché Γa è conservata, vi è un meccanismoche incrementa/diminuisce Γ in presenza di diminuzioni/aumenti di∆A⊥, ossia di moticonvergenti/divergenti del fluido in un piano perpendicolare ad Ω. Il meccanismo se-condo cui tale conversione ha luogo è dovuto alla forza di Coriolis, che, deviando versodestra le traiettorie degli elementi di fluido, induce una rotazione antioraria/oraria permoti convergenti/divergenti in presenza di vorticità ambientaleΩ positiva.

2.8 L A V O R T I C I TÀ P O T E N Z I A L E

Il teorema di Kelvin non vale quando il vettore baroclino B = ∇p×∇%%2 è non nullo e quindi

modifica la forza di un tubo di vortice. Tuttavia è possibile, grazie al teorema di Ertel,identificare una diversa quantità, detta vorticità potenziale, che è conservata.

Se esiste uno scalare λ costante durante il moto del fluido, ossia corrispondente aduna quantità conservata, allora le superfici dove λ è costante corrispondono a superficimateriali permanenti, ossia che si muovono con il fluido, o, in altre parole, la porzionedi tubo di vorticità, inizialmente compresa fra due superfici corrispondenti al valore λ eλ+∆λ, resterà compresa fra di esse durante il moto. Se questo scalare è funzione solodella densità e della pressione del fluido1, allora il vettore baroclino è perpendicolare algradiente di λ,

B ·∇λ= ∇p ×∇%%2 ·

(∂λ

∂p∇p + ∂λ

∂%∇%

)= 0

ed il suo flusso attraverso una superficie dove λ è costante è sempre nullo. In questeipotesi la forza di un tubo di vorticità non è costante, ma è costante il flusso di vorticitàattraverso la superficie che delimita la porzione di tubo considerata. Se il tubo ha unasezione infinitesima δS allora

∂

∂t(ωa ·δS) = 0. (2.18)

Poiché la massa δM contenuta nella porzione di tubo è costante allora δS = δM%δl , dove

δl distanza fra le due superfici vale δl = ∆λ|∇λ| , con ∆λ costante allora l’eq. (2.18) risulta

un’equazione di conservazione:

DΠ

Dt= D

Dt

(ωa · ∇λ

%

)= 0.

La funzione Π è chiamata vorticità potenziale. Pertanto è stato dimostrato il seguenteteorema:

TEOREMA DI ERTEL 2.2 - Se DλDt = 0 e λ è una funzione solo di p e %, allora la vorticità

potenzialeΠ=ωa · ∇λ% è conservata: DΠDt = 0.

Il termine vorticità potenziale esprime la capacità di questo meccanismo di variare lavorticità del fluido secondo le variazioni del gradiente diλ, precisamente di fornire vorti-

1 Ad esempio la densità stessa, se il moto è incompressibile, o la temperatura potenziale in un caso piùgenerale.

24

2.8 L A V O R T I C I TÀ P O T E N Z I A L E

cità quando il gradiente diminuisce e di assorbirla quando esso aumenta, ossia, se il mo-to è incompressibile, quando gli elementi di fluido vengono rispettivamente allungati ocompressi. Nel caso in cui λ≡ %,Π=ωa · ∇%% .

Nel caso di un moto che si sviluppa in modo bidimensionale, u = (u, v ,0), il ruolo diλ è assunto semplicemente da z, Dz

Dt = 0. Perciò, ricordando che DδVDt = 0 e ponendo

%= δMδV come sopra

D

Dt(ωa ·δS) = D

Dt

(ωa

δM

%δl

)= D

Dt

(ωa

δV

∆z

)= δV

D

Dt

(2Ω+ω

H

)= 0,

Π= 2Ω+ωH

,

dove abbiamo indicato H ≡∆z.

25

3M O T I A G R A N D E S C A L A N E L L’ AT M O S F E R A E N E L L’ O C E A N O

Questo capitolo introduce le equazioni necessarie per descrivere i moti a grande scalache caratterizzano le strutture principali della circolazione generale dell’oceano e del-l’atmosfera. In particolare è discussa la semplificazione dinamica che deriva dalla geo-metria quasi-piana dei moti a grande scala in cui la velocità è quasi esclusivamente oriz-zontale, e sono descritte le caratteristiche del regime geostrofico, in cui forza di Coriolise gradiente di pressione si bilanciano.

3.1 L E E Q U A Z I O N I I N C O O R D I N AT E G E O G R A F I C H E

Per studiare i moti dell’atmosfera e dell’oceano rispetto alla superficie terrestre, le equa-zioni della sezione precedente vanno riscritte in coordinate geografiche r (distanza ra-diale), φ (latitudine), λ (longitudine). La componente zonale della velocità u è nelladirezione di λ, la velocità meridionale v nella direzione di φ e la componente verticalew , nella direzione di r . Per semplicità si assume che la risultante della gravità e dellaforza centrifuga sia sempre diretta lungo r. Si ha:

u = r cosφλ, v = r φ, w = r .

In coordinate geografiche la derivata totale di uno scalare ξ è

Dξ

Dt= ∂ξ

∂t+ u

r cosφ

∂ξ

∂λ+ v

r

∂ξ

∂φ+w

∂ξ

∂r

e la divergenza di un vettore v:

div (v) = 1

cosφ

[∂

∂λ

( vλr

)+ ∂

∂φ

(vφ cosφ

r

)+ 1

r 2

∂

∂r(r 2vr cosφ)

].

Pertanto la conservazione della massa diventa:

D%

Dt+ %

r cosφ

∂u

∂λ+ %

r cosφ

∂

∂φ(v cosφ)+ %

r 2

∂r 2w

∂r= 0.

L’equazione del momento ha una struttura complessa. Il termine 2Ω×u dà un contri-buto

2Ω×u = 2Ω(−v sinφ+w cosφ,u sinφ,−u cosφ).

Il contributo non lineare determina la presenza di numerosi termini quadratici che, acausa della curvatura terrestre, non sono confinati esclusivamente nella derivata totale.

27

M O T I A G R A N D E S C A L A

La loro presenza nelle equazioni deriva dalle variazioni dell’orientamento del sistema diriferimento durante il moto di un elemento di fluido. L’equazione risultante è

Du

Dt−

(2Ω+ u

r cosφ

)(v sinφ−w cosφ) =− 1

%r cosφ

∂p

∂λ+Dλ,

Dv

Dt+ w v

r+

(2Ω+ u

r cosφ

)u sinφ=− 1

%r

∂p

∂φ+Dφ,

Dw

Dt− u2 + v2

r−2Ωu cosφ= 1

%

∂p

∂r− g +Dr , (3.1)

dove il vettore (Dλ,Dφ,Dr ) rappresenta il contributo dei termini di attrito viscoso.

3.2 S I S T E M I D I R I F E R I M E N T O C A R T E S I A N I

Le equazioni in coordinate sferiche sono geometricamente complesse e la comprensio-ne della dinamica risulta più immediata in coordinate cartesiane. È ovvio che se ci silimita ad una porzione limitata della superficie terrestre il cui spessore sia piccolo ri-spetto al raggio della Terra si può utilizzare un sistema cartesiano il cui piano x, y siatangente alla superficie.

Riscriviamo le equazioni introducendo le coordinate curvilinee

x = R cosφ0λ, y = R(φ−φ0) z = r −R

dove il sistema di coordinate ha origine al punto λ = 0, φ = φ0, r = R. Introduciamo lescale L, H ,U ,W in modo che le variabili adimensionali

x ′ = x

L, y ′ = y

L, z ′ = z

H, u′ = u

U, v ′ = v

U, w ′ = w

W

siano di ordine uno ed assumiamo che:

H

R¿ 1,

(L

R

)2

¿ 1,L

Rtanφ0 ¿ 1. (3.2)

Queste condizioni indicano uno spessore di fluido piccolo rispetto al raggio terrestre,un’estensione in latitudine e longitudine ristretta ad una piccola porzione della superfi-cie terrestre ed infine che si è distanti dai poli.

Localmente le coordinate curvilinee sono equivalenti ad un sistema di riferimentocartesiano, il cui piano x y è tangente alla superficie sferica. Si noti che utizzando lecoordinate curvilinee, fermandosi al primo ordine nell’espansione di Taylor, si ha:

1

r≈ 1

R

(1− z ′ H

R

), cosφ≈ cosφ0

(1− tanφ0

y ′LR

), sinφ≈ sinφ0

(1+ 1

tanφ0

y ′LR

).

Per cui, nella stessa approssimazione, la derivata totale della densità è

D%

Dt≈ ∂%

∂t+ 1

R cosφ0

(1+ y ′L

R− z ′H

R

)u∂%

∂x+

(1− z ′H

R

)v∂%

∂y+w

∂%

∂z

≈ ∂%

∂t+u

∂%

∂x+ v

∂%

∂y+w

∂%

∂z.

28

3.3 M O T I Q U A S I - P I A N I

e la divergenza della velocità diventa

∂u

∂x

(1+ tanφ0

y ′LR

− z ′HR

)+ ∂v

∂y

(−z ′H

R

)+ ∂w

∂z+ v ′ L2

RUtanφ0 + 2HW

RH

(1− z ′H

R

)w ′ ≈

∂u

∂x+ ∂v

∂y+ ∂w

∂z,

cioè l’equazione della conservazione della massa in coordinate cartesiane, per la cuiderivazione, sono state adottate, implicitamente, per le scale le relazioni:

1

T≡ U

L= W

H,

seguendo lo stesso procedimento per le equazioni del momento, si ottiene il seguentesistema di equazioni:

Du

Dt− f v + f w + 1

%

∂p

∂x= 0,

Dv

Dt+ f u + 1

%

∂p

∂y= 0,

Dw

Dt− f u − u2 + v2

R+ 1

%

∂p

∂z+ g = 0,

1

%

D%

Dt+ ∂u

∂x+ ∂v

∂y+ ∂w

∂z= 0,

dove

f = f0 +βy , f0 = 2Ωsinφ0, β= 2Ωcosφ0

R, f = 2Ωcosφ0.

Queste equazioni costituiscono la cosiddetta approssimazione del piano β. Se si scegliedi porre β = 0 si parla di approssimazione del piano f . Alle medie latitudini si ha β =1.6 ·10−11 ms−1 e f0 = 10−4 s−1.

3.3 M O T I Q U A S I - P I A N I

Una caratteristica dei moti dell’atmosfera e dell’oceano è che si svolgono su uno stratoestremamente sottile. Ciò è stato in parte già sfruttato nella sezione precedente, maha ulteriori implicazioni. Infatti non solo il rapporto H

R è piccolo, ma anche il rapportoHL è piccolo: la scala orizzontale dei moti è molto più grande della scala verticale (la

profondità media dell’oceano è circa 4km mentre la larghezza è di migliaia di kilometri).Quindi H

L ≈ 10−3. Un’importante conseguenza è che la pressione è idrostatica anchequando il fluido è in moto. Usando le stesse scale del capitolo precedente introduciamole due condizioni:

H

L¿ 1, cotφ0

H

L¿ 1.

La seconda condizione indica che siamo sufficientemente distanti dall’equatore.

29

M O T I A G R A N D E S C A L A

La scala dei vari termini nell’equazione della componente x del momento è:

Du

Dt︸︷︷︸UT

− f v︸︷︷︸2Ωsinφ0U

+ f w︸︷︷︸2Ωcosφ0

U HL

+ 1

%

∂p

∂x︸ ︷︷ ︸P%L

= 0.

Abbiamo ipotizzato che U T = L ossia che∣∣∣∂u∂t

∣∣∣∼ |u ·∇u|. In molti casi le non linearitàsono invece minori ed U T ¿ L, ma le conclusioni di questa sezione non dipendono daquesto rapporto. Risulta allora che:∣∣ f w

∣∣∣∣ f v∣∣ ∼ H

Lcotφ0.

per cui l’unica componente della rotazione dinamicamente rilevante è quella perpendi-colare al piano del moto. Inoltre risulta che la scala della pressione è:

P = %U L max

f ,

1

T

.

La scala dei vari termini nell’equazione per la componente z è:

Dw

Dt︸︷︷︸U 2 H

L2

− f u︸︷︷︸2Ωcosφ0U

− u2 + v2

R︸ ︷︷ ︸U 2

R

+ 1

%

∂p

∂z︸ ︷︷ ︸U LH max f , 1

T

+ g︸︷︷︸g

= 0.

Confrontando i vari termini con quello contenente la pressione si ha che∣∣DwDt

∣∣∣∣∣ 1%∂p∂z

∣∣∣ ∼(

H

L

)2

,

∣∣ f u∣∣∣∣∣ 1

%∂p∂z

∣∣∣ ∼ H

Lcotφ0,

∣∣∣u2+v2

R

∣∣∣∣∣∣ 1%∂p∂z

∣∣∣ ∼ H

R

pertanto la gravità è l’unica forza che può bilanciare il gradiente di pressione, e vale ilbilancio idrostatico.

Si consideri uno strato di fluido di densità costante: allora p(x, y , z) = g%(η− z) doveη è il livello della superficie libera del fluido, ed il gradiente orizzontale non dipendedalla coordinata verticale. È plausibile allora sostenere che anche il moto risultante nondipende da z e le derivate secondo z vanno eliminate dalle equazioni del moto. Si hannoallora le equazioni per moti quasi-piani (shallow–water in inglese) non viscosi:

∂u

∂t+u

∂u

∂x+ v

∂u

∂y− f v + g

∂η

∂x= 0, (3.3a)

∂v

∂t+u

∂v

∂x+ v

∂v

∂y+ f u + g

∂η

∂y= 0, (3.3b)

∂

∂t(η−hb)+ ∂

∂xu(η−hb)+ ∂

∂yv(η−hb) = 0; (3.3c)

L’ultima equazione discende dalla relazione di continuità, supponendo D%Dt = 0 ed inte-

grando tra hb , altezza del fondale, ed η, livello della superficie:∫ η

hb

∂u

∂xdz +

∫ η

hb

∂v

∂ydz +

∫ η

hb

∂w

∂zdz = 0;

30

3.4 I L B I L A N C I O G E O S T R O F I C O

dal fatto che ∂∂x

∫ h(x)g(x) f (x, y)dy = ∂

∂x

∫ h(x)g(x)

∂ f (x,y)∂x dy+ f (x,h(x))∂h(x)

∂x − f (x, g (x))∂g (x)∂x , pos-

siamo scrivere∫ η

hb

∂u∂x dz = ∂

∂x

∫ ηhb

u dz−u ∂η∂x +u ∂hb

∂x ed analogamente per v . Poiché infineDηDt = ∂η

∂z +u ·∇η= w |η ed analogamente per hb , si ha

Dη

Dt− Dhb

Dt+ ∂

∂x

∫ η

hb

u dz −u∂η

∂x+u

∂hb

∂x+ ∂

∂y

∫ η

hb

v dz − v∂η

∂y+ v

∂hb

∂y= 0

che può approssimativamente riscriversi come la (3.3c), supponendo una debole dipen-denza di u e v da z.

3.4 I L B I L A N C I O G E O S T R O F I C O

Molte delle approssimazioni del par. 3.3 derivano dall’analisi delle caratteristiche geo-metriche del moto, senza considerarne la dinamica. L’approssimazione geostrofica èinvece valida per i moti, detti a grande scala, in cui la forza di Coriolis ha un ruolodominante, cioè sono fondamentalmente condizionati dalla rotazione terrestre.

I moti a grande scala sono caratterizzati da scale temporali paragonabili al periododi rotazione della Terra. L’effetto della rotazione terrestre è importante per questi moti,mentre è trascurabile per i moti a piccola scala normalmente coinvolti nelle attività uma-ne. La distinzione fra moti a grande e piccola scala viene resa più precisa introducendoil numero di Rossby εR . Se L è l’ampiezza del moto e U la sua velocità caratteristica allorail tempo caratteristico è T = L

U . Il numero di Rossby è definito come

εR = U

2ΩL.

Il numero di Rossby rappresenta simultaneamente il rapporto fra frequenza del moto efrequenza della rotazione terrestre, fra vorticità relativa e vorticità planetaria, fra deriva-ta totale e termine di Coriolis nell’equazione del moto. Un valore piccolo del numero diRossby caratterizza i moti a grande scala.

L’analisi dell’equazione del moto mostra che se il numero di Rossby è piccolo alloravale approssimativamente l’eguaglianza

2Ω×u =−∇p

%−%g, (3.4)

ossia il gradiente di pressione è bilanciato dalla forza di Coriolis. Il fluido che si trova inquesta condizione dinamica viene definito in bilancio geostrofico.

3.5 L E C O O R D I N AT E I S O B A R I C H E

Quando la pressione è idrostatica è spesso vantaggioso sostituire alla coordinata vertica-le la pressione ed ottenere un sistema di equazioni in cui non appaia esplicitamente ladensità del fluido. In questo nuovo sistema di coordinate (x, v , p), chiamate isobariche,ogni elemento di fluido è caratterizzato dalla densità %(x, y , p) e dalla velocità (u, v ,ξ)dove ξ= Dp

Dt .

31

M O T I A G R A N D E S C A L A

In coordinate isobariche la conservazione della massa per un elemento di fluido è

D%δV

Dt= D%δxδyδz

Dt= Dδxδyδp

Dt= 0, (3.5)

per cui la conservazione della massa diventa la conservazione del “volume” nello spaziodefinito dalle coordinate isobariche:

∂u

∂x+ ∂v

∂y+ ∂ξ

∂p= 0 (3.6)

ossia la “velocità” è sempre non-divergente.Per esprimere le equazioni del moto in coordinate isobariche occorre esprimere me-

diante esse il gradiente orizzontale della pressione. Si noti che (ignorando tempora-neamente la variabile y) se due spostamenti successivi ∆x e ∆z sono scelti in modo dariportarsi sulla stessa superficie isobarica si ha(

∂p

∂x

)z∆x +

(∂p

∂z

)x∆z = 0 ⇒

(∂p

∂x

)z=−

(∂p

∂z

)x

(∂z

∂x

)p=−g%

(∂z

∂x

)p

,

che, introducendo la funzione geopotenzialeΦ definita dalla relazione

dΦ=−dp

%= g dz (3.7)

diventa(∂Φ

∂x

)z= g

(∂z

∂x

)p

.

Le due equazioni per le componenti orizzontali del momento diventano, ignorando icontributi viscosi

Du

Dt− f v =−∂Φ

∂x, (3.8a)

Dv

Dt+ f u =−∂Φ

∂y. (3.8b)

Si noti che, a z costante, la pressione e il geopotenziale hanno massimi e minimi incorrispondenza nonostante il segno ‘−’ presente nella definizione (3.7); ciò può esse-

re visto dal fatto che 0 =(∂p∂x

)p=

(∂p∂x

)z+ ∂p

∂z

(∂z∂x

)p=

(∂p∂x

)z− g%

(∂z∂x

)p

e dunque(∂p∂x

)z=

g%(∂z∂x

)p= %

(∂Φ∂x

)p

.

3.6 I L V E N T O T E R M I C O

Perché il numero di Rossby rimanga permanentemente piccolo durante il moto del flui-do occorre che la vorticità relativa resti sempre trascurabile rispetto a quella planetaria.Pertanto l’equazione della vorticità deve semplificarsi in:

(2Ω ·∇)u−2Ω(∇·u) =−∇%×∇p

%2 (3.9)

32

3.7 M O D E L L O A D U E S T R AT I

ossia la produzione di vorticità per baroclinicità deve essere bilanciata dallo stiramentoe dall’inclinazione degli elementi di fluido.

In un sistema di riferimento con l’asse z parallelo adΩ, le componenti dell’equazione(3.9) si scrivono come:

2Ω∂u

∂z= 1

%2

(∂p

∂z

∂%

∂y− ∂p

∂y

∂%

∂z

), (3.10a)

2Ω∂v

∂z= 1

%2

(∂p

∂z

∂%

∂x− ∂p

∂x

∂%

∂z

), (3.10b)

2Ω

(∂u

∂x+ ∂v

∂y

)= 1

%2

(∂p

∂x

∂%

∂y− ∂p

∂y

∂%

∂x

). (3.10c)

Le prime due collegano la derivata verticale della velocità alle derivate orizzontali delladensità. Se la velocità del fluido ad un livello di riferimento è nota, queste equazioni con-sentono di ricavare i profili della velocità orizzontale da misure di densità, utilizzandoper la pressione l’espressione idrostatica. Poiché le variazioni di densità sono comu-nemente dovute a variazioni di temperatura le correnti (o i venti) descritti da (3.10a) e(3.10b) sono detti venti termici.

Le equazioni del vento termico assumono una forma particolarmente semplice incoordinate isobariche. Infatti dalla definizione di geopotenziale, derivando rispetto ap le equazioni (3.8) in bilancio geostrofico, si ottiene:

f∂v

∂p= ∂2Φ

∂x∂p= 1

%2

∂%

∂x, f

∂u

∂p= ∂2Φ

∂y∂p=− 1

%2

∂%

∂y.

oppure

f∂v

∂z=−g

%

(∂%

∂y

)p

, f∂u

∂z= g

%

(∂%

∂y

)p

.

relazioni che esprimono il tasso di variazione verticale della velocità in funzione dellavariazione orizzontale di densità lungo superfici isobariche.

3.7 M O D E L L O A D U E S T R AT I

Supponiamo di disporre di due fluidi di differente densità, immiscibili tra loro (o tali dapoter essere considerati immiscibili). Se %1 è la densità (costante) del liquido superiore,con velocità orizzontali u1 e v1, e %2 è quella (costante) del liquido inferiore, con velocitàorizzontali u2 e v2, supponiamo che l’interfaccia si trovi ad una profondità z = −H1 +h(x, y , t ) (all’equilibrio z =−H1) e che la superficie superiore (all’equilibrio individuatada z = 0) sia tale che z = η(x, y , t ). Allora p1 = %1g (η− z) per −H1 +h < z < η. Poichésupponiamo di essere in assenza di vorticità, le equazioni shallow-water diventano

∂u1

∂t=−g

∂η

∂x,

∂v1

∂t=−g

∂η

∂y, (3.11a)

∂(η+H1 −h)

∂t+H1

(∂u1

∂x+ ∂v1

∂y

)= 0; (3.11b)

33

M O T I A G R A N D E S C A L A

derivando la (3.11b) rispetto al tempo ed utilizzando le (3.11a)

∂2

∂t 2 (η−h) = H1g

(∂2

∂x2 + ∂2

∂y2

)η≡ g H1∇2

Hη. (3.12)

Analogamente per il secondo strato, p2 = %1g (η+ H1 − h) + %2g (−H1 + h − z) per z <−H1 +h, per cui, introducendo la gravità ridotta

g ′ = g%2 −%1

%2,

si ottengono le equazioni

∂u2

∂t=−g

%1

%2

∂η

∂x− g ′ ∂h

∂x,

∂v2

∂t=−g

%1

%2

∂η

∂y− g ′ ∂h

∂y, (3.13a)

∂h

∂t+H2

(∂u2

∂x+ ∂v2

∂y

)= 0, (3.13b)

detta z =−H2 la quota del fondale. Procedendo come prima si ottiene

∂2h

∂t 2 = H2∇2H

(gη− g ′η+ g ′h

). (3.14)

Le equazioni ottenute governano dunque l’evoluzione del sistema. Il problema può es-sere tuttavia semplificato di molto ponendo h(x, y , t ) =µη(x, y , t ), doveµ è una costante;perciò le equazioni (3.12) e (3.14) si riducono a

∂2η

∂t 2 = c2e∇2

Hη, (3.15)

dove

c2e = g H1

1−µ = H2g − g ′(1−µ)

µ⇒ c4

e − g Hc2e + g g ′H1H2 = 0,

dove H = H1 + H2, profondità totale all’equilibrio. Definendo c2e = g He , He profondità

equivalente, l’equazione diventa

g H 2e − g H He + g ′H1H2 = 0;

poiché negli oceani le variazioni di densità sono molto piccole (dell’ordine del 3%%), g ′g

è infinitesimo e risolvendo l’equazione si ottengono due diversi valori per c2e : il maggio-

re è c20 = g H

(1− g ′

gH1 H2

H 2 + . . .)

per cui ηh ≈ H

H2e u2

u1≈ 1 nel limite g ′

g → 0 (in cui il liqui-

do diventa a densità uniforme, modo barotropico1); la radice minore vale invece c21 =

g ′H1 H2

H

(1+ g ′

gH1 H2

H 2 + . . .)

e corrisponde ai rapporti approssimati ηh ≈ − g ′H2

g H e u2u1

≈ −H1H2

(modo baroclino, ovvero la pressione non è costante sulle superfici di densità costante).Per gli oceani solitamente c1 = 2−3m · s−1, per cui la profondità equivalente è di 0.5−1m.

1 Solitamente per barotropico si intende che la pressione è costante sulle superfici a densità costante, e dun-que è costante all’interfaccia: ciò è solo approssimativamente vero ma si è soliti usare questa espressioneegualmente.

34

3.7 M O D E L L O A D U E S T R AT I

3.7.1 Il modo baroclino

Il modo barotropico si ottiene semplicemente trattando il fluido come dotato di densitàuniforme. Ci sono due particolari approssimazioni da fare per ottenere, invece, il modobaroclino. La prima consiste nel supporre η molto piccola rispetto alle fluttuazioni sul-l’interfaccia (approssimazione a coperchio rigido). Così la terza delle equazioni (3.11b)diventa

−∂h

∂t+H1

(∂u1

∂x+ ∂v1

∂y

)= 0. (3.16)

Il secondo tipo di approssimazione consiste nel rimpiazzare %1

%2con l’unità nelle (3.13a)

(approssimazione di Bousinnesq)

∂u2

∂t=−g

∂η

∂x− g ′ ∂h

∂x,

∂v2

∂t=−g

∂η

∂y− g ′ ∂h

∂y. (3.17)

Sottraendo le (3.17) dalle (3.11a) si ottiene

∂u

∂t= g ′ ∂h

∂x,

∂v

∂t= g ′ ∂h

∂y, (3.18)

con u = u1−u2 e v = v1−v2 differenze di velocità tra i due strati; sottraendo ora 1H2

volte

la (3.13b) ad 1H1

volte la (3.16) si ottiene

∂u

∂x+ ∂v

∂y−

(1

H1+ 1

H2

)∂h

∂t= 0. (3.19)

Combinando la (3.19) con la (3.18) otteniamo l’equazione delle onde

∂2h

∂t 2 = c21∇2

H h, c21 = g ′H1H2

H1 +H2, quadrato della velocità del modo baroclino.

Il valore di c21 ottenuto corrisponde a quello già trovato nel limite g ′

g → 0. Se He = c21

g ,allora la precedente corrisponde all’equazione

1

g He= 1

g ′H1+ 1

g ′H2.

Tipicamente g ′ = 0.03m · s−2, H1 = 400m, H2 = 4000m, per cui He ≈ 1m.

35

4O N D E L I N E A R I

Le onde lineari sono ottenute dalle equazioni del moto ignorando i termini avvettivi.Questo è giustificato se U 2

L ¿ UT , dove L è la lunghezza d’onda, T il suo periodo e U la

velocità del fluido, cioè se la velocità del fluido è molto più piccola della velocità di fasedell’onda. Questo corriponde a onde in cui lo spostamento degli elementi di fluido U Tè molto più piccolo della lunghezza d’onda L, ossia onde poco ripide, nel caso che ilsegnale propagato sia l’elevazione.

Il materiale contenuto in questo capitolo si suddivide in due categorie: onde irrota-zionali ed onde in presenza di rotazione (limitatamente al solo regime barotropico).

Le onde irrotazionali si ottengono dalla soluzione dell’equazione del moto in un siste-ma di riferimento non rotante, in assenza di viscosità

Du

Dt=−∇p

%,

con velocità non divergente, ∇ ·u = 0, e densità costante. Esse sono onde con vorticitànulla, dove la velocità verticale è dello stesso ordine di quella orizzontale e pertanto nonvale il bilancio idrostatico. Sono onde di piccola scala (il numero di Rossby è εR =∞) lacui lunghezza raggiunge al massimo il chilometro ed il cui moto è prodotto dall’instabi-lità dell’interfaccia aria–mare in presenza di vento. Le forze di richiamo sono la gravitàe la tensione superficiale.

L’effetto della rotazione è analizzato identificando le soluzioni ondulatorie delle equa-zioni lineari per moti quasi-piani:

∂u∂t − f v =−g ∂η

∂x , (4.1a)

∂v∂t + f u =−g ∂η

∂y , (4.1b)

∂Hu∂x + ∂H v

∂y + ∂η∂t = 0, (4.1c)

dove si suppone hb indipendente dal tempo e η−hb ≈ −hb ≡ H , profondità del fluido.La linearizzazione dell’ultima equazione implica η¿ H , che è un requisito non neces-sariamente soddisfatto da onde di piccola ampiezza, ma richiede H À L, ossia una lun-ghezza d’onda non grande rispetto alla profondità. In questo tipo di onde la pressione èidrostatica e la vorticità non nulla.

In un piano f , in assenza di variazioni di profondità le uniche onde descrivibili dalle(4.1) sono le onde di Poincarè–Sverdrup, la cui frequenza angolare verifica la condizione$ ≥ f . La presenza di una costa permette una ulteriore soluzione detta onda di Kelvin.Se sono incluse variazioni della profondità o di f sono possibili onde di bassa frequenza($< f ) dette onde di Rossby.

37

O N D E L I N E A R I

4.1 O N D E D I G R AV I TÀ

Dalla (2.17) si ottiene che in un sistema di riferimento non rotante se la vorticità è nullainizialmente, essa resta nulla anche nella continuazione del moto. Si può quindi assume-re per moti a piccola scala che il moto sia irrotazionale. L’irrotazionalità è garantita se lavelocità è espressa come gradiente di una funzione φ, potenziale della velocità, u =∇φ.In tal caso l’equazione ∇·u = 0 diventa l’equazione di Laplace:

∂2φ

∂x2 + ∂2φ

∂z2 = 0, (4.2)

dove, per semplicità, l’asse x è stato orientato nella direzione della componente oriz-zontale della velocità. Tale equazione è valida fra il fondo del fluido z = −H e la suasuperficie z = η. Quest’ultima non è costante ma è parte della soluzione che richiedela determinazione di η e di u. Cercando una soluzione ondulatoria dell’eq. (4.2) nellaforma φ(x, z, t ) = Z (z)ei (kx−$t ) si ottiene Z (z) = c1ekz +c2e−kz . Imponendo come primacondizione al contorno, sul fondo del fluido, al livello z =−H , l’annullarsi della velocitàverticale ∂φ

∂z = 0, si ottiene:

φ(x, z, t ) = w0 cosh[k(z +H)]ei (kx−$t ).

Per completare la soluzione occorre specificare la relazione di dispersione che collegak ad $, la cui derivazione richiede altre due condizioni. La prima, detta cinematica,impone che durante il moto gli elementi di fluido che si trovano alla superficie η non sidiscostino da essa:

D

Dt(η− z) = 0 ovvero

∂η

∂t+u

∂η

∂x= w . (4.3)

La seconda, detta dinamica, deriva dall’equazione del moto eq. (2.16) che si semplifica,per moti non viscosi ed irrotazionali, in

∂u

∂t+∇|u|2

2=−∇p

%+∇Ψ. (4.4)

Riscrivendo l’eq. (4.4) in termini del potenziale della velocità si ottiene

∇(∂φ

∂t+ 1

2|∇φ|2 + p

%+ g z

)= 0,

da cui segue che alla superficie del fluido, dove p = p0, pressione atmosferica,(∂φ

∂t+ 1

2|∇φ|2

)z=η

+ gη= 0. (4.5)

La (4.3) e la (4.5) vengono linearizzate espandendo il potenziale della velocità attorno lalivello z = 0, ed ignorando i termini di ordine kη, ossia assumendo onde poco ripide. Siottiene la coppia di equazioni:

∂η

∂t= ∂φ

∂z

∣∣∣∣z=0

= k sinh(kH)ei (kx−$t ), −gη= ∂φ

∂t

∣∣∣∣z=0

=−i$cosh(kH)ei (kx−$t ), (4.6)

38

4.1 O N D E D I G R AV I TÀ

Figura 6: Propagazione delle onde di gravità verso la costa.

per cui, facendo il rapporto membro a membro

∂η

∂t=−i

g k

$tanh(kH)η⇔−∂

2η

∂t 2 = ig k

$tanh(kH)

∂η

∂t= g 2k2

$2 tanh2(kH)η,

che permette direttamente di trovare la relazione di dispersione

⇒$2 = g k tanh(kH). (4.7)

Le onde descritte dalla relazione (4.7) sono dispersive, con le onde corte che viaggianopiù lentamente delle onde lunghe. Nel limite di profondità infinita (kH →∞) si ha $2 =g k. Nel limite opposto, per acqua poco profonda, (kH → 0) si ha $2 = k2g H e le ondenon sono più dispersive, ma viaggiano tutte alla stessa velocità. Nei casi intermedi lavelocità di propagazione diminuisce con la profondità, per cui, durante l’avvicinamentoalla costa, se il profilo del fondale è sufficientemente regolare, le linee di cresta tendonoad allinearsi alle isolinee di profondità (fig. 6). L’espressione completa della soluzione è:

η(x, t ) = η0 cos(kx −$t ),

u(x, z, t ) = u0coshk(z +H)

coshkHcos(kx −$t ), dove u0 = η0

g k

$,

w(x, z, t ) = w0sinhk(z +H)

coshkHsin(kx −$t ), dove w0 = η0

g k

$,

p(x, z, t ) = p0coshk(z +H)

coshkHcos(kx −$t ), dove p0 = η0g%.

Si noti che mentre velocità orizzontale e pressione sono in fase con l’elevazione, la ve-locità verticale è in anticipo di π

2 . Il moto delle particelle (fig. 7) di fluido segue quindidelle ellissi, la cui ampiezza diminuisce con la profondità, con il rapporto fra i moduliW0 e U0 delle due componenti della velocità dato da:

W0

U0= tanh[k(z +H)] .

39

O N D E L I N E A R I

Figura 7: Traiettorie degli elementi di fluido in un’onda di garvità in acqua profonda ed acquabassa.

La pressione non è idrostatica ma si attenua esponenzialmente con la profondità (inacqua profonda) con un tasso di diminuzione più accentuato per le onde corte.

4.2 O N D E D I P O I N C A R È – S V E R D R U P

Un modo per ottenere le soluzioni ondulatorie consiste nel manipolare le equazioni eriscriverle nella sola variabile η. A tal fine si sommano la derivata rispetto a x della (4.1a)e la derivata rispetto a y della (4.1b) ottenendo:

∂

∂t

(∂u

∂x+ ∂v

∂y

)+ f

(∂u

∂y− ∂v

∂x

)=−g

(∂2η

∂x2 + ∂2η

∂y2

)e si sottrae la derivata rispetto a x della (4.1b) dalla derivata rispetto ad y della (4.1a)ottenendo l’equazione della vorticità

∂

∂t

(∂u

∂y− ∂v

∂x

)= f

(∂u

∂x+ ∂v

∂y

). (4.8)

Le due equazioni combinate assieme consentono di eliminare la vorticità ottenendoun’equazione per η e la divergenza della velocità. Se il fondo è piatto, cioè H = costante,la divergenza è proporzionale alla derivata dell’elevazione:

∂u

∂x+ ∂v

∂y=− 1

H

∂η

∂t.

Risulta quindi un’equazione per la sola η:

∂

∂t

[∂2η

∂t 2 + f 2η− g H

(∂2η

∂x2 + ∂2η

∂y2

)]= 0. (4.9)

Sostituendo nell’eq. (4.9) la soluzione ondulatoria η= η0ei (kx+l y−$t ) si ottiene la rela-zione di dispersione per le onde di Poincarè–Sverdrup:

$2 = f 2 + g H(k2 + l 2). (4.10)

40

4.3 O N D E D I K E LV I N

Le onde sono dispersive, hanno sempre frequenza maggiore di f , e quindi non rappre-sentano moti geostrofici. Se la loro lunghezza d’onda è molto inferiore al raggio di Ros-

sby RR =p

g Hf allora il parametro di Coriolis è trascurabile, le onde non risentono della

rotazione e non sono dispersive.

Per individuare il tipo di moto, si consideri un sistema di riferimento orientato nelladirezione di propagazione dell’onda (l = 0). In tale sistema le equazioni si semplificanoe sostituendo nelle equazioni la soluzione ondulatoria (η0,u0, v0)ei (kx−$t ) si ottiene:∣∣∣∣u0

v0

∣∣∣∣= $

f, arg(v0) = arg(u0)+ π

2.

Pertanto gli elementi di fluido percorrono in senso orario (se f > 0) delle ellissi con l’assemaggiore orientato nella direzione di propagazione dell’onda. Il rapporto fra gli assidell’ellisse vale $

f ed è eguale ad 1 quando le onde hanno la frequenza minima f (ondeinerziali). In tale caso il moto segue orbite circolari ed è determinato dalla sola forza diCoriolis, senza alcun contributo della pressione.

4.3 O N D E D I K E LV I N

La presenza di un limite sul dominio di integrazione, ossia di una costa, permette unasoluzione dalle caratteristiche particolari, detta onda di Kelvin, in cui la componentedella velocità parallela alla costa è identicamente nulla. Per semplicità, il sistema di rife-rimento si considera orientato con l’asse y parallelo alla costa, che costituisce il bordodestro del bacino. Se u ≡ 0 allora le equazioni del moto quasi–piano diventano:

f v = g∂η

∂x, (4.11a)

∂v

∂t=−g

∂η

∂y, (4.11b)

∂η

∂t+H

∂v

∂y= 0. (4.11c)

La soluzione cercata è del tipo η0φ(x)ei (k y−$t ), cioè è periodica in y , ma ha una strutturain x, che deve essere determinata dalla dinamica del moto.

Le eq.(4.11) si riducono ad una equazione nella sola η:

∂2η

∂t 2 = g H∂2η

∂y2 , (4.12)

da cui segue la relazione di dispersione

$2 = g Hk2. (4.13)

La relazione di dispersione dell’onda di Kelvin è eguale a quella delle onde di Poincarè–Sverdrup nel limite per raggio di Rossby grande rispetto alla lunghezza d’onda. Sono pos-sibili onde di Kelvin di qualsiasi frequenza e pertanto le onde di Kelvin rappresentano lerisposte a forzanti di frequenza troppo bassa per generare le onde di Poincarè–Sverdrup

41

O N D E L I N E A R I

e troppo alta per generare le onde di Rossby. I principali esempi sono le forzanti mareali,che hanno periodi concentrati in due intervalli centrati attorno alle 12h e 24h.

La struttura trasversale dell’onda di Kelvin viene determinata dalla (4.11a). Sostituen-do in essa la soluzione ondulatoria trovata si ottiene:

φ′(x) = f sign(k)√g H

φ(x),

da cui φ(x) = exp

(f sign(k)p

g Hx

). Mentre in direzione parallela alla costa l’onda di Kelvin è

sinusoidale, nella direzione trasversale cresce esponenzialmente e tale crescita è inter-rotta dalla presenza della costa. La soluzione con k < 0 diverge e non è accettabile. Seinvece il bacino fosse limitato da una costa sul bordo sinistro non sarebbe accettabilela soluzione con k > 0. Pertanto l’onda di Kelvin si propaga sempre tenendo la costa adestra.

La rotazione, anche se non appare nella relazione di dispersione, ha un ruolo fon-damentale nell’onda di Kelvin. La componente non nulla della velocità è infatti in bi-lancio geostrofico, e la tendenza della forza di Coriolis a indurre nel fluido una velocitàtrasversale è bilanciata esattamente dal gradiente di pressione.

4.4 O N D E D I R O S S B Y

L’onda di Rossby è estremamente importante in geofisica. Infatti essa mostra come sia-no possibili moti ondulatori in bilancio geostrofico e spiega l’evoluzione dei vortici digrande scala del fluido. Essa consegue dalle equazioni della quasi-geostrofia e la suadinamica è spiegata in modo semplice dalla conservazione della vorticità potenziale edalla divergenza del trasporto geostrofico nel piano β (vedasi prossima sezione).

Le onde di Poincarè-Sverdrup costituiscono la sola soluzione possibile delle equa-zioni del moto quasi-piano in assenza di coste. Se però il fondale non è piatto o siadotta l’approssimazione piano-β, allora delle nuove soluzioni diventano possibili. Ri-petendo il procedimento seguito per le onde di Poincarè–Sverdrup, ma considerandoH = H0 +γy , f = f0 +βy , si ottiene

∂

∂t

(∂u

∂x+ ∂v

∂y

)+ f

(∂u

∂y− ∂v

∂x

)+βu =−g

(∂2η

∂x2 + ∂2η

∂y2

),

∂

∂t

(∂u

∂y− ∂v

∂x

)−βv = f

(∂u

∂x+ ∂v

∂y

), (4.14)

che utilizzando

∂u

∂x+ ∂v

∂y=− 1

H0

∂η

∂t−γv

si riduce a

∂

∂t

[∂2η

∂t 2 + f 2η− g H0

(∂2η

∂x2 + ∂2η

∂y2

)]−βH0

(∂u

∂t+ f v

)+γ

(∂2v

∂t 2 + f 2v

)= 0. (4.15)

42

4.4 O N D E D I R O S S B Y

Per moti lenti, ossia $¿ f , la (4.15) si semplifica in

∂

∂t

[f 2η− g H0

(∂2η

∂x2 + ∂2η

∂y2

)]+H0 f

(γ f

H0−β

)v = 0. (4.16)

Introduciamo ora βT = β− f γH0

; osserviamo che la (4.16) contiene ancora l’incognita vche non è stata espressa in funzione di η. Per moti a bassa frequenza vale però il bilanciogeostrofico, e quindi v = g

f∂η∂x e quindi si ha un’equazione per la sola η:

∂

∂t

[f 2η− g H0

(∂2η

∂x2 + ∂2η

∂y2

)]−βT g H0

∂η

∂x= 0. (4.17)

Sostituendo nell’eq. (4.17) la soluzione ondulatoria si ottiene la relazione di dispersio-ne dell’onda di Rossby:

f 2 + g H0

(k2 + l 2 +βT

k

$

)= 0 =⇒$=−R2

RkβT

1+R2R (k2 + l 2)

. (4.18)

Vi sono due limiti della relazione di dispersione: per lunghezze d’onda molto più grandidel raggio di Rossby, 1

k2+l 2 À R2R , e per numeri d’onda molto piccoli. Per onde di Rossby

lunghe si ottiene

$=−R2RβT k;

per onde di Rossby corte si ottiene

$=−βTk

k2 + l 2 .

In entrambi i casi la velocità di fase è negativa, mentre la velocità di gruppo è positivaper onde corte e negativa per onde lunghe.

L’onda di Rossby può essere identificata anche dalla conservazione della vorticità po-tenziale, prendendone il limite lineare e cercandone soluzioni ondulatorie che verifichi-no il bilancio geostrofico.

4.4.1 Dinamica dell’onda di Rossby

Il moto dell’onda di Rossby si presta a due diverse interpretazioni corrispondenti ai duelimiti per onde lunghe e corte.

Si consideri il limite ad onde lunghe in presenza di un gradiente di vorticità ambienta-le, ad esempio l’aumento di vorticità planetaria in direzione meridionale. Nel limite on-de lunghe il ruolo dell’elevazione superficiale, ossia dello stiramento-compressione deitubi di vorticità è trascurabile e la superficie superiore del fluido si può assumere piana.Se una colonna di fluido nella posizione x0 viene spostata in direzione del gradiente divorticità orientato come l’asse y, poiché aumenta la propria vorticità ambientale, dimi-nuisce la propria vorticità relativa, cioè acquisisce la tendenza a ruotare in verso orario.Il moto che ne segue sposta le adiacenti colonne di fluido, in posizione x0+∆x e x0−∆x,in direzione parallela al gradiente della vorticità ambientale. Conseguentemente, tali

43

O N D E L I N E A R I

Figura 8: Meccanismi di evoluzione dell’onda lunga (a sinistra) e corta (a destra) di Rossby.

colonne di fluido acquisiscono una vorticità relativa di segno tra di loro opposto. Il mo-to che ne segue riporta la colonna x0 verso la sua posizione originaria lungo y e muovelungo la stessa direzione altre colonne adiacenti determinando un moto ondulatorio(vedasi fig.8).

Se viceversa le onde sono corte la deviazione della superficie del fluido dalla posizio-ne di equilibrio ha il ruolo fondamentale. Se, in presenza di un moto convergente sicrea, localmente, un’elevazione, per il bilancio geostrofico ad essa è associata una circo-lazione antioraria. A causa dell’effetto β la circolazione risultante è meno intensa nelladirezione dove la vorticità ambientale è maggiore. Pertanto si creano ai lati est ed ovestdell’elevazione una zona di divergenza e di convergenza, che producono una depressio-ne ed una elevazione della superficie. Le circolazioni geostrofiche associate a tali duedeformazioni, producono, per via dell’effetto β, una divergenza nella regione fra essacompresa, per cui l’elevazione che vi era in precedenza collocata scompare ed inizianouna serie di oscillazioni (vedasi fig.8).

44

Parte II

M O D E L L I S T I C A N U M E R I C A

5I N T E G R A Z I O N E N U M E R I C A D I P D E

5.1 P R E M E S S A : M E T O D O A L L E D I F F E R E N Z E F I N I T E

Sia data una funzione f (x) continua; definiamo fi = f (xi ), con xi successione mono-tona strettamente crescente in R (con passo eventualmente non costante). Si ipotizzache il passo sia tale da poter considerare lineare la funzione tra xi e xi+1 (poiché non cisono informazioni specifiche, non ha senso andare oltre il primo ordine). Pertanto perx ∈ (xi−1, xi ) si ha f (x ′) ≈ fi + fi− fi−1

∆x (x ′ − xi−1) = fi

(xi−x ′∆x

)+ fi−1

(1− xi−x ′

∆x

)≡ α(x ′) fi +

(1−α(x ′)) fi−1, dove abbiamo adoperato l’approssimazione lineare per stimare f ′. Per

valutare l’entità dell’errore, dallo sviluppo di Taylor f (x±ε) = f (x)+∑∞j=1

(±ε) j

j !∂ j f (x)

∂x j , per

cui fi−1 ≈ f (x ′)−∆x1 f ′(x ′)+ ∆x21

2 f ′′(x ′) mentre fi ≈ f (x ′)+∆x2 f ′(x ′)+ ∆x22

2 f ′′(x ′), ∆x1 =x ′ − xi−1, ∆x1 +∆x2 = ∆x, da cui si ottiene fi∆x1+ fi−1∆x2

∆x ≈ f (x ′)+ ∆x1∆x22 f ′′(x ′), per cui

l’errore è dell’ordine di ∆x2, è nullo agli estremi ed è massimo al centro dell’intervallo.

Le serie di Taylor permettono di disporre di una approssimazione della funzione e

delle sue derivate: da fi−1 = fi −∆x ∂ f∂x (xi )+ ∆x2

2∂2 f∂x2 (xi )+O(∆x2) si ottiene che ∂ f

∂x (xi ) =fi− fi−1

∆x + ∆x2∂2 f∂x2 (xi ), ovvero ∂ f

∂x (xi ) = fi−1− fi

∆x a meno di ordini superiori in ∆x.

Se consideriamo x = xi+xi−12 , allora fi− 1

2∓ 12≈ f (x)∓ ∆x

2 f ′(x)+ ∆x2

8 f ′′(x)∓ ∆x3

48 f ′′′(x) ⇒f ′(x) = fi− fi−1

∆x + ∆x3

24 f ′′′(x), ovvero l’errore va a zero più rapidamente e dunque possiamopensare di introdurre una griglia sfasata che permetta il calcolo delle derivate con unabuona approssimazione.

Se ora consideriamo la derivata seconda, fi±1 ≈ fi ±∆x f ′i + ∆x2

2 f ′′i ± ∆x3

3! f ′′′i + ∆x4

4! f (4)i ,

da cui f ′′i ≈ fi+1−2 fi+ fi−1

∆x2 − ∆x2

12 f (4)i , ovvero l’errore scala come ∆x2. Solitamente questa

approssimazione viene ritenuta sufficiente.

Per le derivate superiori possiamo formalmente scrivere

......

......

1 −2∆x 2∆x2 ...

1 −∆x ∆x2

2

...

1 0 0...

1 ∆x ∆x2

2

...

1 2∆x 2∆x2 ......

......

...

︸ ︷︷ ︸

A

fi

f ′i

f ′′i

f ′′′i

f (4)i

f (5)i

f (6)i...

=

...fi−2

fi−1

fi

fi+1

fi+2...

,

per cui il problema si riduce all’inversione della matrice A. Per procedere si può utiliz-zare, ad esempio, il metodo di Gauss: scrivendo la matrice (A‖1), dove 1 è la matrice

47

I N T E G R A Z I O N E N U M E R I C A D I P D E

identità, si può passare tramite una serie di combinazioni lineari tra righe, alla matrice(1‖B): allora B = A−1.

5.2 O N D E I N U N A D I M E N S I O N E

L’equazione

∂ f

∂t+ c

∂ f

∂x= 0 (5.1)

descrive la propagazione di un’onda piana unidirezionale in un mezzo le cui proprietàvariano su una scala alquanto maggiore della lunghezza d’onda, come può essere un’on-da marina in acqua bassa in un canale di profondità variabile in modo lento e regolare(in modo che la velocità di propagazione c possa essere supposta costante). Tale equa-zione rappresenta il più semplice esempio di PDE (Partial Differential Equation) iperbo-lica. Una soluzione analitica della (5.1) con condizione iniziale f (x,0) = Aei kx è datada f (x, t ) = Aei k(x−ct ). La soluzione numerica con i metodi alle differenze finite si basasulla sostituzione delle derivate in forma continua con rapporti incrementali di tipo

∂ f

∂x≈ ∆ f

∆x,

basati sulla sostituzione della variabile continua con una sua rappresentazione discre-ta, f (x, t ) → f n

j = f (x0 + j∆x, t0 +n∆t ). La soluzione consiste nel calcolare i valori di

f nj in successione per n = 1, . . . , N , a partire dalla condizione iniziale f (x,0) → f 0

j =f (x0 + j∆x,0), fino a coprire l’intero intervallo di integrazione. La rappresentazione del-le derivate può essere effettuata in modi diversi. Indicando con i l’indice sulla nostravariabile e con j gli indici sulle altre si ha:

C E N T E R E D∂ f∂x = fi+1, j− fi−1, j

2∆x +O(∆x2);

F O R W A R D∂ f∂x = fi+1, j− fi , j

∆x +O(∆x);

B A C K W A R D∂ f∂x = fi , j− fi−1, j

∆x +O(∆x);

analogamente si procede per t . La scelta della rappresentazione determina le proprietàdello schema numerico.

Una proprietà fondamentale per lo schema numerico è la sua stabilità, ossia la pro-prietà di fornire soluzioni finite (non divergenti) a seguito della sua iterazione. Il teore-ma di Lax assicura che uno schema stabile e consistente (in cui le rappresentazioni dellederivate convergono all’espressione analitica se il passo di integrazione tende a zero)converge alla soluzione corretta. La stabilità viene usualmente studiata con il criteriodi von Neumann, consistente nel sostituire alla soluzione l’espressione f n

j = ξnei j k∆x everificare che, per ogni valore di n, |ξ| < 1, il che garantisce che la soluzione resta fini-ta. Il fattore di amplificazione ξ caratterizza il comportamento del metodo numerico,cioè un’alterazione dell’ampiezza o della fase della componente dovuta alla proceduradi integrazione.

48

5.2 O N D E I N U N A D I M E N S I O N E

5.2.1 Schemi di integrazione espliciti

S C H E M A D I I N T E G R A Z I O N E F T C S Lo schema FTCS (Forward in Time, Centered inSpace) è il caso più semplice ma ha applicazioni pratiche limitate, in quanto è intrinse-camente instabile. Adottando tale schema per la (5.1) si ottiene

f n+1j − f n

j

∆t+ c

f nj+1 − f n

j−1

2∆x= 0,

ossia la formula esplicita per la soluzione è data da

f n+1j = f n

j −Cf n

j+1 − f nj−1

2, (5.2)

dove C = c ∆t∆x = velocità fisica

velocità numerica . Con un po’ di algebra si vede che |ξ| =√

1+C 2 sin2 k∆x >1 ∀k, dunque lo schema è sempre instabile. Per rimediare a questa instabilità, si puòadoperare lo schema di Lax, che è una variazione dell’FTCS: esso consiste nel sostituiref n

j con la sua media tra due punti adiacenti, in modo che la (5.2) diventi

f n+1j =

f nj+1 + f n

j−1

2−C

f nj+1 − f n

j−1

2, (5.3)

da cui ξ= cosk∆x − iC sink∆x e dunque |ξ| ≤ 1 se C ≤ 1. La stabilità, dunque, dipendedal cosiddetto criterio di Courant o di CFL (Courant–Friedrichs–Lévy):

C = c∆t

∆x= velocità fisica

velocità numerica< 1, (5.4)

che impone delle restrizioni sul rapporto ∆x∆t una volta fissata c. Il criterio CFL è una

condizione ricorrente per la stabilità degli schemi di integrazione numerica di PDE ditipo iperbolico. Il risultato, apparentemente sorprendente, ha una logica spiegazioneche è evidente se riscriviamo la (5.3) nella forma

f n+1j = f n

j +f n

j+1 −2 f nj + f n

j−1

2−C

f nj+1 − f n

j−1

2,

che rappresenta in uno schema FTCS l’equazione

∂ f

∂t+ c

∂ f

∂x= ν∂

2 f

∂x2 , ν= ∆x2

2∆t.

Lo schema di Lax introduce dunque, in definitiva, un termine diffuso che limita la cre-scita dell’instabilità.

S C H E M A D I I N T E G R A Z I O N E C T C S O L E A P F R O G Uno schema molto usato è il CTCS

(Centered in Space, Centered in Time) che produce

f n+1j − f n−1

j

2∆t+ c

f nj+1 − f n

j−1

2∆x= 0,

49

I N T E G R A Z I O N E N U M E R I C A D I P D E

da cui la formula risolutiva

f n+1j = f n−1

j −C(

f nj+1 − f n

j−1

),

che fornisce l’equazione ξ2 +2iCξsink∆x −1 = 0 ⇒ ξ=−iC sink∆x ±√

1−C 2 sin2 k∆x.Se C ≤ 1 il radicando è positivo, la radice reale e |ξ|2 = 1 per ogni valore di k. Se inveceC > 1 esistono valori di k per cui, essendo C 2 sin2 k∆x > 1, il radicando è negativo e laradice ha solo parte immaginaria. In tal caso

|ξ|2 =(−C sink∆x ±

√C 2 sin2 k∆x −1

)2= 2C 2 sin2 k∆x −1∓2C sink∆x

√C 2 sin2 k∆x −1

|ξ|2 > 1∓2C sink∆x√

C 2 sin2 k∆x −1,

e dunque almeno in un caso è maggiore di 1. Similmente alla correzione di Lax all’FTCS,lo schema CTCS è stabile se la condizione CFL è soddisfatta.

Lo schema CTCS presenta la particolarità che i punti pari e dispari del grigliato com-putazionale (assimilabili alle caselle bianche e nere di una scacchiera) sono disaccop-piati, ossia la loro integrazione numerica procede in modo indipendente determinandola possibilità che le due soluzioni numeriche divergano tra loro. Per evitare tale pro-

blema, usualmente si aggiunge all’equazione un termine diffusivo ε(

f nj+1 −2 f n

j + f nj−1

),

dove ε¿ 1 rappresenta un coefficiente arbitrario e piccolo.

S C H E M A D I I N T E G R A Z I O N E U P W I N D Un ulteriore schema è lo schema upwind. Sef rappresenta la concentrazione di uno scalare generico, c f ne rappresenta il flusso. Loschema upwind calcola il tasso di variazione della concentrazione di f nel punto j comedifferenza tra flusso entrante e uscente tenendo conto del verso secondo cui avviene lapropagazione

f n+1j − f n

j

∆t + cf n

j − f nj−1

∆x = 0 se c > 0,f n+1

j − f nj

∆t + cf n

j+1− f nj

∆x = 0 se c ≤ 0,

che, nel caso c > 0, si traduce nella formula risolutiva

f n+1j = (1−C ) f n

j +C f nj−1 ⇒|ξ| = 1+4C (C −1)sin

k∆x

2,

ovvero |ξ| > 1 se C > 1. La sostituzione nella prima delle precedenti dello sviluppo inserie di Taylor di f n+1

j , f nj+1, f n

j−1 al secondo ordine fornisce

∂ f

∂t+ c

∂ f

∂x=−1

2

∂2 f

∂t 2 ∆t − c

2

∂2 f

∂x2∆x,

che, se f è soluzione di (5.1), diventa

∂ f

∂t+ c

∂ f

∂x= (C −1)

C∆x

2

∂2 f

∂x2 .

L’equazione realmente risolta dallo schema contiene un supplementare termine diffusi-vo, che determina la presenza di una diffusione numerica nella soluzione calcolata.

50

5.3 M E T O D I I M P L I C I T I

5.2.2 Dispersione numerica

La discretizzazione spaziale introduce una dispersione numerica che altera la velocità acui un’onda si propaga. Questo è già evidente osservando che il fattore di amplificazioneξ risultante dall’analisi della stabilità è caratterizzato non solo da un modulo ma ancheda una fase, che, a parte il caso che speciale in cui ξ = 1, implica una dipendenza dak nella propagazione di onde che in base alla soluzione analitica sono non dispersive.Questo errore degli schemi numerici deriva dalla discretizzazione spaziale. Ad esempioscrivendo la (5.1) come

∂ f j

∂t+ c

f j+1 − f j−1

2∆x= 0, f j = F (t )ei k j∆x ,

risulta

∂ f j

∂t+ i kc f j = 0, c = c

sink∆x

k∆x, (5.5)

che mostra come solo le onde lunghe si propaghino con la velocità corretta c.

5.3 M E T O D I I M P L I C I T I

Gli schemi risolutivi presentati fino ad ora sono detti espliciti, in quanto il calcolo dellivello temporale incognito n+1 avviene mediante una semplice combinazione algebri-ca dei valori della funzione incognita al livello temporale noto n. Se si utilizza il livellotemporale incognito n +1 anche nel calcolo della derivata spaziale mediante la mediapesata dei livelli temporali n ed n +1 la situazione si complica:

f n+1j − f n

j

∆t+ c

ω(

f n+1j+1 − f n+1

j−1

)+ (1−ω)

(f n

j+1 − f nj−1

)2∆x

= 0,

dove se ω = 0 si ha lo schema esplicito FTCS, mentre se ω = 1 lo schema si dice com-pletamente implicito mentre se ω = 1

2 è detto semi-implicito. Il vantaggio degli schemiimpliciti è la loro stabilità senza la necessità di soddisfare il criterio di Courant. L’analisidi von Neumann nel caso completamente implicito fornisce

ξ= 1− iC2 sink∆x

1+ iC2 sink∆x

,

che, posto θk = arctan(C

2 sink∆x), diventa ξ= e−2iθk , ovvero |ξ| = 1. Comunque, poiché

la stabilità non implica l’accuratezza della soluzione numerica, il passo di integrazionetemporale non può essere scelto arbitrariamente grande. Inoltre l’integrazione con me-todi impliciti richiede la soluzione di un sistema algebrico lineare, ovvero l’inversione diuna matrice tridiagonale. Si consideri il caso semi-implicito, che si può scrivere come

f n+1j − f n

j +Cf n

j+1 − f nj−1 + f n+1

j+1 − f n+1j−1

4= 0

51

I N T E G R A Z I O N E N U M E R I C A D I P D E

ed in forma matriciale AF (n+1) =−BF (n), dove

A =

1 C

4 0 0 · · ·−C

4 1 C4 0 · · ·

0 −C4 1 C

4 · · ·...

......

.... . .

, B =

−1 C

4 0 0 · · ·−C

4 −1 C4 0 · · ·

0 −C4 −1 C

4 · · ·...

......

.... . .

,

che si risolve ponendo F (n+1) = −A−1BF (n), cioè mediante il calcolo dell’inverso dellamatrice A.

5.4 E Q U A Z I O N E D I D I F F U S I O N E

Consideriamo ora l’equazione che descrive un processo diffusivo in una dimensione

∂ f

∂t= ν∂

2 f

∂x2 . (5.6)

Essa rappresenta il caso più semplice di una categoria di equazioni denominate parabo-liche ed è fondamentale nei processi ambientali che descrivono la diffusione di caloreo di una generica quantità scalare per effetti molecolari o di turbolenza (salinità, com-ponenti della quantità di moto, concentrazioni di componenti chimiche o popolazionibiologiche. . . ).

5.4.1 Schema FTCS

Lo schema FTCS per la sua integrazione è

f n+1j − f n

j

∆t= ν

f nj+1 −2 f n

j + f nj−1

∆x2 ,

dove si è usata l’approssimazione ∂2 f∂x2 = f n

j+1−2 f nj + f n

j−1

∆x2 +O(∆x2); dall’equazione si ottienedunque

f n+1j = (1−2D) f n

j +D( f nj+1 + f n

j−1), D = ν ∆t

∆x2 .

La consueta analisi di stabilità di von Neumann restituisce ξ= 1−4D sin2 k∆x2 , per cui lo

schema è stabile se D ≤ 12 , cioè se il passo di integrazione è inferiore alla metà del tempo

caratteristico in cui il segnale diffonde lungo la distanza ∆x.

5.4.2 Schemi impliciti di integrazione

Uno schema completamente implicito

f n+1j − f n

j = D(

f n+1j+1 −2 f n+1

j + f n+1j−1

)fornisce ξ = 1

1+4D sin2 k∆x2

che ha modulo minore di 1 ∀D e dunque è sempre stabile. Lo

schema intermedio tra quello completamente implicito e quello esplicito si ottiene con

52

5.5 S TA G G E R I N G E G R I G L I E M U LT I D I M E N S I O N A L I

il calcolo della derivata seconda media fra livelli temporali n ed n + 1. Tale schema èincondizionatamente stabile ed è detto di Crank–Nicholson:

f n+1j − f n

j

∆t= ν

2

(f n+1

j+1 −2 f n+1j + f n+1

j−1

∆x2 +f n

j+1 −2 f nj + f n

j−1

∆x2

),

da cui, infatti, |ξ| = 1−2D sin2 k∆x2

1+2D sin2 k∆x2

≤ 1 ∀D . Come discusso in precedenza, si pone tuttavia il

problema dell’inversione della matrice tridiagonale.Il caso completamente implicito in una dimensione è un esempio di problema alge-

bricamente gestibile in modo poco dispendioso. Consideriamo

f n+1j − f n

j = a j

(f n+1

j+1 −2 f n+1j + f n+1

j−1

),

dove ai = νi∆t∆x2 sono gli elementi della matrice tridiagonale che dovrebbe essere invertita,

mentre νi indica una generica dipendenza del coefficiente di diffusione dallo spazio. Laprecedente si riscrive

f n+1i =αi f n+1

i+1 +βi , (5.7a)

αi = ai

1+ai +ai−1(1−αi−1), βi =

ai−1βi−1 + f ni

1+ai +ai−1(1−αi−1), (5.7b)

che, fissati i valori α1 e β1, fornisce una procedura sequenziale per il calcolo in ordinedi i crescente. L’equazione (5.7a) permette di calcolare f n+1

i in ordine di i decrescente el’equazione viene risolta evitando l’inversione algebrica di una matrice tridiagonale.

5.4.3 Onde in presenza di diffusione

L’equazione delle onde in presenza di diffusione è

∂ f

∂t+ c

∂ f

∂t= ν∂

2 f

∂x2FTCS−−−→

f n+1j − f n

j

∆t= ν

f nj+1 −2 f n

j + f nj−1

∆x2 − cf n

j+1 − f nj−1

2∆x

⇒ f n+1j = f n

j −Cf n

j+1 − f nj−1

2+D

(f n

j+1 −2 f nj + f n

j−1

). (5.8)

Dall’analisi della stabilità ξ = 1 − iC sink∆x − 2D(1 − cosk∆x) ⇒ |ξ| ≤ 1 purché sianosoddisfatte le condizioni C 2 ≤ 2D ≤ 1. La presenza di diffusione rende stabile dunqueuno schema che, applicato all’onda semplice, era instabile.

Uno schema implicito per la diffusione

f n+1j = f n

j −Cf n

j+1 − f nj−1

2+D

(f n+1

j+1 −2 f n+1j + f n+1

j−1

),

fornisce ξ= 1−iC sink∆x1+2D(1−cosk∆x) , per cui |ξ| ≤ 1 ∀D ,C .

5.5 S TA G G E R I N G E G R I G L I E M U LT I D I M E N S I O N A L I

Le equazioni considerate fino ad ora rappresentano dei paradigmi di riferimento chein realtà hanno limitato interesse pratico. Questa sezione ha l’obiettivo di avvicinarsi

53

I N T E G R A Z I O N E N U M E R I C A D I P D E

maggiormente a casi di interesse geofisico in cui siano richiesti e correntemente appli-cati schemi numerici. Le equazioni shallow-water (4.1) costituiscono da un lato la basesu cui costruire schemi tridimensionali più complessi, dall’altro si prestano a modellinumerici barotropici (intesi in modo semplicistico come modelli monostrato) di circo-lazione con cui si possono descrivere tsunami, inondazioni costiere, maree, e che sonostati anche utilizzati per i primi modelli di circolazione generale.

5.5.1 Equazioni shallow-water monodimensionali

Le equazioni shallow-water in una dimensione accoppiano l’evoluzione del livello dellasuperficie η del fluido a quella della velocità u in base al bilancio di volume e pressioneidrostatica. In assenza di termine di Coriolis esse sono

∂η

∂t+H

∂u

∂x= 0,

∂u

∂t+ g

∂η

∂x= 0, (5.9)

dove H è la profondità del fluido, g l’accelerazione di gravità; le precedenti fornisconocome soluzione onde non dispersive con velocità c =√

g H in entrambi i versi.L’integrazione delle equazioni (5.9) può essere effettuata con gli schemi di integrazio-

ne FTCS-Lax o CTCS. In entrambi i casi lo schema è stabile se il criterio CFL è soddisfat-to. L’introduzione della discretizzazione spaziale introduce una dispersione numericasimile a (5.5).

Il trasferimento delle equazioni shallow-water su una griglia computazionale eviden-zia i vantaggi di introdurre griglie sfasate (staggered) dove ogni variabile dipendente ècalcolata in punti diversi. Se η ed u sono calcolate su due griglie diverse, ciascuna conpasso∆x ma fra loro sfasate di ∆x

2 , le differenze centrate nello spazio possono essere cal-

colate comeun

j+1−unj

∆x eηn

j −unj−1

∆x nei punti j dove sono richieste le derivate temporali di ηed u. In questo modo si raddoppia la precisione nel calcolo delle derivate spaziali senzaaumentare il costo computazionale.

5.5.2 Onde in griglie multidimensionali

L’equazione delle onde in due dimensioni ha la forma

∂ f

∂t+ cx

∂ f

∂x+ cy

∂ f

∂y= 0;

l’analisi di von Neumann impone una forma di soluzione del tipo f nj ,l = ξnei (kx j+ky l )∆,

dove ∆=∆x =∆y . Lo schema di integrazione FTCS-Lax restituisce

f n+1j ,l =

f nj+1,l + f n

j−1,l + f nj ,l+1 + f n

j ,l−1

4−Cx

f nj+1,l − f n

j−1,l

2−Cy

f nj ,l+1 − f n

j ,l−1

2,

Cx = cx∆t

∆, Cy = cy

∆t

∆, (5.10)

da cui ξ= coskx∆+cosky∆

2 − i(Cx sinkx∆+Cy sinky∆

). Con un po’ di algebra si ha

|ξ|2 = 1−(coskx∆+cosky∆

)2

4−(

Cx sinkx∆+Cy sinky∆)2−(

sin2 kx∆+ sin2 ky∆)[1

2−

(C 2

x +C 2y

)],

54

5.5 S TA G G E R I N G E G R I G L I E M U LT I D I M E N S I O N A L I

per cui lo schema è stabile se C 2x +C 2

y ≤ 12 , ossia

p2√

c2x + c2

y ≤ ∆∆t , che è un caso partico-

lare della relazione generale in N dimensioni

C ≤ 1pN

, C = ∆t

∆

√√√√ N∑i=1

c2i , (5.11)

dove ci sono le componenti della velocità di propagazione delle onde. Si noti che rispet-to al caso monodimensionale, in N dimensioni la condizione di Courant è più restrittivaper via del termine

pN .

5.5.3 Diffusione in griglie multidimensionali

L’equazione di diffusione in due dimensioni diventa

∂ f

∂t= ν

(∂2 f

∂x2 + ∂2 f

∂y2

); (5.12)

la sua integrazione sulla base di uno schema esplicito fornisce

f n+1j − f n

j = D(δ2

x f nj ,l +δ2

y f nj ,l

),

D = ν∆t

∆2 , δ2x f n

j ,l = f nj+1,l −2 f n

j ,l + f nj−1,l , δ2

y f nj ,l = f n

j ,l+1 −2 f nj ,l + f n

j ,l−1.

L’analisi di von Neumann restituisce come criterio di stabilità D ≤ 14 , che nel caso N -

dimensionale diventa D ≤ 12N .

Utilizzando il metodo di Crank–Nicholson per integrare la (5.12) si ottiene

f n+1j − f n

j = Dδ2

x f n+1j ,l +δ2

x f nj ,l +δ2

y f n+1j ,l +δ2

y f nj ,l

2,

che sarebbe associata ad una matrice non più tridiagonale ma con termini più diffusi,detta pentadiagonale, che dovrebbe essere invertita.

Un’alternativa, che permette di ridurre il problema all’inversione di una matrice tri-diagonale, è quello di separare il passo di integrazione in due semipassi (se siamo in Ndimensioni, il passo viene diviso in N parti), dove viene calcolato in modo semimplicitoun solo contributo alla volta, mentre i rimanenti restano espliciti. Ad esempio in duedimensioni

f n+1j − f

n+ 12

j = Dδ2

x fn+ 1

2

j ,l +δ2y f n+1

j ,l

2.

Come lo schema di Crank–Nicholson, anche questo si prova essere incondizionatamen-te stabile.

5.5.4 Equazioni di shallow-water, forza di Coriolis e griglie di Arakawa

Un fattore fondamentale nelle equazioni per i fluidi geofisici è la forza di Coriolis, cheintrodotta nelle equazioni fornisce come abbiamo visto

∂η

∂t+H

(∂u

∂x+ ∂v

∂y

)= 0,

∂u

∂t− f v + g

∂η

∂x= 0,

∂v

∂t+ f u + g

∂η

∂y= 0.

55

I N T E G R A Z I O N E N U M E R I C A D I P D E

n

n12

n1

j j + ½ j j j + ½ j j j + ½ j

n

n12

n1

j j + ½ j j j + ½ j

u

v

η

AA BB CC

DD EE

Figura 9: Griglie di Arakawa.

La discretizzazione di queste equazioni su un grigliato a passo costante può essere attua-ta con cinque diverse modalità, chiamate griglie di Arakawa, indicate con A, B, C, D, E

(fig. 9). Ad esempio nella griglia A le tre variabili dipendenti sono calcolate in punti coin-cidenti, nella griglia B u e v sono sfasate di mezzo passo di griglia rispetto ad η eccetera.Le prestazioni di uno schema numerico dipendono dall’efficienza con cui esso consentedi calcolare le variabili dipendenti e le loro derivate (evitando di ricorrere ad interpola-zioni spaziali) che sono richieste per l’avanzamento temporale della soluzione. Nessunoschema consente di evitare interpolazioni e tutte implicano delle deviazioni della disper-sione numerica rispetto a quella esatta delle onde. Ad esempio la griglia B è ottimale peril calcolo dei termini di Coriolis, ma non per quello del gradiente di pressione; la grigliaC verifica esattamente le esigenze contrarie.

5.6 B O U N D A R Y VA L U E P R O B L E M S

I problemi affrontati fino ad ora vengono classificato come iperbolici o parabolici e ri-chiedono il calcolo simultaneo dei valori della funzione incognita ad un nuovo livellotemporale in funzione di quelli al livello temporale precedente e noto; essi sono inoltreproblemi caratterizzati da una condizione iniziale. La classe di problemi complementa-re viene detta delle equazioni ellittiche e richiede il calcolo di una soluzione che verifichil’equazione differenziale all’interno di un dominio sul cui bordo sono definite delle con-dizioni al contorno: in questo caso la soluzione deve essere calcolata simultaneamentesu tutti i punti del dominio. Nel caso dei problemi alle condizioni iniziali è in generaleimportante la stabilità della soluzione e l’ottimizzazione del passo temporale adottato.Nel caso di problemi alle condizioni al contorno è in generale importante la velocità concui l’algoritmo (di natura iterativa) converge verso la soluzione finale.

56

5.6 B O U N D A R Y VA L U E P R O B L E M S

Il paradigma delle equazioni di natura ellittica è rappresentato dall’equazione di Pois-son

∂2 f

∂x2 + ∂2 f

∂y2 = %(x, y), (5.13)

che in termini di differenze finite diventa

fi+1, j −2 fi , j + fi−1, j

∆x2 + fi , j+1 −2 fi , j + fi , j−1

∆y2 = %i , j +O(∆x2,∆y2). (5.14)

I metodi risolutivi di natura iterativa si basano su una procedura del tipo f k+1i , j = f k

i , j +cki , j ,

secondo la solita notazione con l’accortezza che k questa volta indicizza le iterazioni enon il tempo. La relazione permette di trovare la soluzione all’iterazione k +1 in funzio-ne del valore noto al passo precedente mediante un termine correttivo: la convergen-za è verificata quando il termine correttivo diventa minore di una quantità prefissatasufficientemente piccola.

I metodi iterativi più usati per risolvere le equazioni ellittiche sono i metodi di Jacobi,di Gauss–Seidel e SOR (Successive Over Relaxation)

M E T O D O D I J A C O B I L’equazione (5.14) può essere riscritta esplicitando il valore nelpunto i , j in funzione dei punti adiacenti

fi , j = 1

2

11∆x2 + 1

∆y2

(fi+1, j + fi−1, j

∆x2 + fi , j+1 + fi , j−1

∆y2 +%i , j

),

che può essere interpretata come una relazione per il calcolo di f all’iterata k + 1 infunzione dei valori all’iterata k. Il metodo così ottenuto è detto metodo di Jacobi:

f k+1i , j = 1

2

11∆x2 + 1

∆y2

(f k

i+1, j + f ki−1, j

∆x2 +f k

i , j+1 + f ki , j−1

∆y2 +%i , j

). (5.15)

M E T O D O D I G AU S S – S E I D E L Il metodo di Gauss–Seidel consiste in una variante delmetodo di Jacobi, in cui nel membro di destra vengono inseriti, al posto di f k , i valoricorrispondenti f k+1 se disponibili. Ad esempio, se l’integrazione procede lungo righe,cioè aumentando i a j fissato, i valori nei punti (i −1, j ) e (i , j −1) sono già stati calcolatiquando si raggiunge il punto (i , j ) e pertanto si può scrivere

f k+1i , j = 1

2

11∆x2 + 1

∆y2

(f k

i+1, j + f k+1i−1, j

∆x2 +f k

i , j+1 + f k+1i , j−1

∆y2 +%i , j

).

Un ulteriore schema si ottiene considerando l’espressione (5.15) come una stima virtua-le intermedia da utilizzare assieme al valore precedente in una media pesata di f k+1

i , j :

f k+1i , j = (1−ω)φk+1

i , j + ω

2

11∆x2 + 1

∆y2

(f k

i+1, j + f k+1i−1, j

∆x2 +f k

i , j+1 + f k+1i , j−1

∆y2 + %i , j

ω

),

dove φk+1i , j è l’espressione di f k+1

i , j restitutita dalla (5.15) ed ω è il peso assegnato al finedi accelerare la convergenza verso la soluzione (in generale 1 <ω< 2).

57

I N T E G R A Z I O N E N U M E R I C A D I P D E

5.7 M E T O D O S P E T T R A L E

Supponiamo di avere una funzione f (x) e di poterla sviluppare come f (x) =∑∞j=1 a jφ j (x),

dove φ j j è una base nello spazio L2(R). Se l’equazione differenziale ha la forma L f (x) =%(x), dove L è un operatore differenziale lineare, si ha che

r (x) =L

(N∑

j=1a jφ j (x)

)−%(x)

costituisce una misura della bontà dell’approssimazione che prevede di troncare ad Nlo sviluppo. Si può scegliere a questo punto di minimizzare ‖r‖ = ∫

r 2(x)dx, oppureporre r (xk ) = 0 per certi valori di un reticolo, o ancora che

⟨r ,φ j

⟩= ∫r (x)φ(x)dx = 0 ∀ j

(metodo di Galërkin), o, come variante,∫

r (x)θ j (x)dx = 0 per ogni elemento θ j di unafamiglia prefissata.

Nel caso di equazioni alle derivate parziali con variabili indipendenti x e t , i coefficien-ti a j sono supposti dipendenti dal tempo. Il metodo spettrale applicato a L = ∂t + c∂x

consiste nel considerare una base ortonormalizzataφ j ,⟨φ j ,φk

⟩= δ j k , conφ j = 1p2π

ei j x

e di adoperare il metodo di Galërkin

1

2π

∫e−i kx

(∂

∂t+ c

∂

∂x

)∑j

a j (t )ei j x dx = 0 ⇒ ∂t ak + i kcak = 0.

M O D E L L O D I S T O M M E L Il modello di Stommel prevede l’introduzione di un para-metro ε= ν

β in un’equazione nella forma

εφxx −φx +1 = 0, x ∈ [0,1].

In genere si utilizzano, come base, seni e coseni. Imponendo come condizioni al contor-no φ(0) =φ(1) = 0, allora φ(x) =∑

j a j sin(π j x). Sostituendo nell’equazione di Stommel,

tenendo in considerazione che∫ 1

0 sin( jπx)sin(kπx)dx = δ j k

2 ed imponendo il criterio diGalërkin si ottiene

εakk2π2

2+∑

ja j j k

1− (−1) j+k

j 2 −k2 = 1− (−1)k

kπ⇔ (1+ A)a =τ.

Il costo computazionale è legato al calcolo dell’inversa della matrice 1+ A.

5.8 M E T O D O A G L I E L E M E N T I F I N I T I

Supponiamo di avere l’equazione delle onde in presenza di diffusione ∂t u + c∂x u =

ν∂2xx u. Sia u = ∑

i uiφi , con ad esempio φ j =

0 se∣∣x −x j

∣∣>∆x,

1− |x−x j |∆x se

∣∣x −x j∣∣<∆x,

dove ∆x è

una certa ampiezza prefissata. Le funzioni non sono ortogonali:⟨φ j ,φk

⟩ = 23∆xδ j k +

∆x3

6 δ j±1,k . Adoperando la stessa tecnica precedente si arriva all’equazione

∂t ui−1 +4∂t ui +∂t ui+1

6+ c

ui+1 −ui−1

2∆x−νui+1 +2ui +ui−1

∆x2 = 0.

58