Corso di Laurea in Fisica - infn.it coinvolte pongono severi problemi computazionali, ... (che si...

43
UNIVERSIT ` A DEGLI STUDI DI PAVIA FACOLT ` A DI SCIENZE MM. FF. NN. Corso di Laurea in Fisica Equazioni di Navier-Stokes e automi cellulari Tesi di Laurea di Giovanni Balossini Relatore Oreste Nicrosini Anno Accademico 2002/2003

Transcript of Corso di Laurea in Fisica - infn.it coinvolte pongono severi problemi computazionali, ... (che si...

UNIVERSITA DEGLI STUDI DI PAVIAFACOLTA DI SCIENZE MM. FF. NN.

Corso di Laurea in Fisica

Equazioni di Navier-Stokese automi cellulari

Tesi di Laurea di Giovanni Balossini

Relatore Oreste Nicrosini

Anno Accademico 2002/2003

Indice

Indice i

Introduzione 1

Capitolo 1. Le equazioni fondamentali della fluidodinamica 31.1. Gradiente di velocita e deformazioni 31.2. Tensore degli sforzi 41.3. Descrizione di Eulero e derivata sostanziale 61.4. Equazioni di bilancio 71.5. Equazioni di Navier-Stokes 91.6. Risoluzione analitica: moti paralleli incomprimibili 10

Capitolo 2. Generalita sugli automi cellulari 132.1. Definizione di automa cellulare 132.2. Life 142.3. Classi di Wolfram 162.4. Formulazione matematica 172.5. Esempi e applicazioni 182.6. Automi cellulari e equazioni differenziali 19

Capitolo 3. Gas reticolare 213.1. Caratteristiche generali 213.2. Microdinamica 223.3. Modello HPP 233.4. Modello FHP-1 243.5. Altri modelli 26

Capitolo 4. Limite macrodinamico 294.1. Distribuzione di equilibrio 294.2. Equazioni della fluidodinamica 31

Capitolo 5. Simulazioni di gas reticolari 355.1. Fluido bidimensionale attorno ad un ostacolo 355.2. Instabilita di Kelvin-Helmholtz 365.3. Aggregazione di particelle 37

Bibliografia 39

Ringraziamenti 41

i

Introduzione

In anni recenti si e assistito ad un crescente interesse nella cosiddetta“Scienza della Complessita”nelle sue piu varie accezioni. Questo dato difatto trae origine da almeno due considerazioni: da un lato, l’interessescientifico che i “sistemi complessi” suscitano di per se, dall’altro ladisponibilita di risorse di calcolo sempre piu potenti, capillari e a costisempre piu accessibili.

In questo contesto, un tema di particolare interesse e la dinamica deifluidi. Questa puo, infatti, essere considerata un paradigma di comples-sita, in particolare con riferimento ai regimi turbolenti. Le equazionicoinvolte pongono severi problemi computazionali, che possono essereaffrontati sia con tecniche tradizionali sia con procedure computazio-nali innovative. Scopo di questa tesi e quello di illustrare il legametra descrizione dinamica tradizionale e moderno approccio ad automicellulari.

Questo lavoro si articola in cinque capitoli. In particolare, i primidue saranno propedeutici agli altri, trattando separatamente le due par-ti fondamentali costituenti l’oggetto della tesi: le equazioni di Navier-Stokes e gli automi cellulari. Nelle rimanenti sezioni si studieranno piuin dettaglio particolari automi cellulari, i modelli a gas reticolare, e leloro applicazioni.

Il primo capitolo e dedicato, come gia detto, alle equazioni fonda-mentali della fluidodinamica. La derivazione classica che viene propostaparte da considerazioni sul tensore degli sforzi, per poi introdurre la de-rivata sostanziale nell’ambito della descrizione di Eulero del moto deifluidi. Lo studio delle equazioni di bilancio permette, poi, di arrivaread ottenere le equazioni di Navier-Stokes. Queste, come verra ricorda-to piu volte in seguito, non hanno una soluzione analitica generale, masolo in casi particolari, come il regime alla Poiseuille, la cui illustrazionechiude il capitolo.

Ad una descrizione generale sugli automi cellulari e sulle loro ca-ratteristiche salienti e dedicato il secondo capitolo. Dopo una defini-zione formale, si procedera ad un semplice esempio: il “gioco” Life diJohn Horton Conway. Ad un accenno alle classi di Wolfram ed allaformulazione matematica seguiranno altri esempi ed una conclusioneriguardante il rapporto tra automi cellulari ed equazioni differenziali,che e, in un certo senso, introduttivo alle parti successive.

1

2 INTRODUZIONE

Gli ultimi tre capitoli sono strettamente legati a particolari automicellulari, i modelli a gas reticolare. Il terzo ricopre un ruolo introdutti-vo, con l’illustrazione delle caratteristiche generali, della microdinamicae di alcuni tra i modelli piu importanti. La relazione con le equazio-ni della fluidodinamica sara sviluppata nel quarto capitolo, piuttostotecnico, in cui sara utilizzato un approccio basato in primo luogo sullameccanica statistica.

La tesi si conclude con una breve rassegna di esempi applicativi deimodelli a gas reticolare.

CAPITOLO 1

Le equazioni fondamentali della fluidodinamica

In questo primo capitolo si intende presentare la derivazione classicadelle equazioni di Navier-Stokes e in generale delle equazioni su cui sibasa la meccanica dei fluidi.

Per far questo si rende, prima di tutto, necessario porsi nel limitedi mezzo continuo, trattando, cioe, la materia non come aggregato diparticelle (che si possono considerare elementari), bensı supponendo lamassa distibuita con continuita nello spazio. A questo punto si puoprendere in considerazione una porzione di fluido le cui dimensioni so-no decisamente inferiori a quelle caratteristiche a livello macroscopico,rimanendo, tuttavia, piu grandi di quelle della dinamica microscopica,cosı che si possano considerare, nello sviluppo in serie di Taylor di ognigrandezza, solo i termini lineari. A questa porzione si da il nome diparticella fluida.

1.1. Gradiente di velocita e deformazioni

Per studiare la cinematica di una particella fluida si sceglie unaterna di assi cartesiani che all’istante t0 ha origine nel baricentro Gdella particella stessa. La velocita del baricentro e uG ed e tale che, senell’intervallo di tempo ∆t questo si e spostato dal punto G al puntoG′, si ha GG′ = uG∆t.

Un generico punto P non si sposta, in generale, con velocita uG,tuttavia si puo operare lo sviluppo1:

u = uG +∂u∂xi

∣∣∣∣G

xi +∂2u

∂xi∂x j

∣∣∣∣G

xix j + . . . .

Considerando solo i termini lineari si nota, allora, che lo spostamentodi un generico punto di una particella fluida e descritto dal tensoregradiente di velocita nel baricentro. Omettendo il pedice G si scrive:

∇u =

∂u1∂x1

∂u1∂x2

∂u1∂x3

∂u2∂x1

∂u2∂x2

∂u2∂x3

∂u3∂x1

∂u3∂x2

∂u3∂x3

.

1Qui e nel seguito si adottera la convenzione di Einstein: dove un indice appareripetuto e sottointesa la somma su quell’indice.

3

4 1. LE EQUAZIONI FONDAMENTALI DELLA FLUIDODINAMICA

Il tensore gradiente di velocita puo essere decomposto nella somma di

un tensore antisimmetrico ~~r e di uno simmetrico ~~e:

r i j =12

(∂ui

∂x j− ∂u j

∂xi

)

e

ei j =12

(∂ui

∂x j+

∂u j

∂xi

).

Se il moto e descritto dal solo tensore simmetrico la velocita, in unsistema di riferimento solidale con il baricentro, e

u =~~e·xe il moto risulta di deformazione e non di rotazione, in cui il tensoreantisimmetrico interviene in modo analogo.

1.2. Tensore degli sforzi

In un fluido, sia ∆S una superficie e P un punto su ∆S. Sia, poi,n il versore normale in P (si veda la figura 1.1). In base al verso di nsi individuano la faccia positiva e quella negativa della superficie. Sidefinisca, infine, ∆F la forza di superficie esercitata dal fluido che sitrova sul lato positivo su quello che si trova sul lato negativo.

P

n

∆F ∆S

Figura 1.1. Elemento di superficie.

Per come e stato definito il mezzo continuo, la funzione

τn.= lim

∆S→0

∆F∆S

esiste continua e di classe Cn. τn si chiama sforzo in P relativo allanormale n. E immediato notare che τ−n =−τn.

1.2. TENSORE DEGLI SFORZI 5

Si consideri, ora, il tetraedro di Cauchy (figura 1.2). Indicata con Mla massa del fluido contenuto nel tetraedro, la condizione di equilibrioe data da:

τiSi− τnSn = Mg.

Inoltre

Si = Sn cosαni = Snni M = ρhSn

3,

dove h e l’altezza relativa alla faccia Sn, quindi:

τini− τn = ρh3

g.

Facendo tendere a zero il volume del tetraedro mantendolo sempresimile a se stesso, quindi facendo tendere h a zero, si ha:

τn = n ·~~τ .

c3

c1

c2

nSn

S2

S3

S1

Figura 1.2. Tetraedro di Cauchy.

Si puo dimostrare ([2]) che:

(1.1) τmn =(−p+λ

∂uk

∂xk

)δmn+2µemn,

dove p e la pressione, mentre λ e µ sono i coefficienti di viscosita. Sinoti che la 1.1 lega il tensore degli sforzi al tensore gradiente di velocitanei fluidi detti newtoniani.

6 1. LE EQUAZIONI FONDAMENTALI DELLA FLUIDODINAMICA

1.3. Descrizione di Eulero e derivata sostanziale

La descrizione del moto del baricentro delle particelle fluide puoessere effettuata in due differenti modi:

lagrangiano: seguendo il moto di ciascuna particella;euleriano: valutando cio che accade in un fissato punto dello

spazio al passare del tempo relativamente a particelle differen-ti.

Riferendosi a quest’ultimo modo, si definisca la derivata sostanziale diuna grandezza caratteristica A (temperatura, pressione, velocita, . . . )come:

DA

D t.= lim

∆t→0

A(P′)−A(P)∆t

.

Con P si e indicato il punto di coordinate x occupato all’istante t dauna data particella, che, dopo un intervallo di tempo ∆t si e spostatain P′, di coordinate x+∆x (vedi figura 1.3).

PQ

c3

c2

c1

Figura 1.3. Rappresentazione del moto delle particellefluide P e Q.

Sia, ora, Q il punto occupato all’istante t da una particella fluida cheal tempo t +∆t si trova, invece, nella posizione P. Le derivate parzialirispetto al tempo e alle coordinate spaziali risultano essere:

∂A∂t

∣∣∣∣x=cost

= lim∆t→0

A(x, t +∆t)−A(x, t)∆t

∂A∂xi

∣∣∣∣t=cost

= lim∆x→0

A(x+∆x, t)−A(x, t)∆x

.

1.4. EQUAZIONI DI BILANCIO 7

La variazione locale (cioe in P≡Q′) di A e dunque:

A(Q′

)−A(P) =∂A∂t

∆t .

Poiche

A(P′

)−A(P) ={

A(P′

)−A(Q′

)}+

{A

(Q′

)−A(P)}

e

A(P′

)−A(Q′

)=

∂A∂x j

∆x j

si ottiene, quindi:

(1.2)DA

D t=

∂A∂t

+∂A∂x j

dx j

dt=

∂A∂t

+∂A∂x j

u j .

1.4. Equazioni di bilancio

1.4.1. Teorema del trasporto di Reynolds.Si consideri ora un volume di fluido V, racchiuso da una superficie S,

che contiene sempre le stesse particelle. Si calcoli la derivata sostanzialedell’integrale di A(x, t) su V:

DD t

V(t)AdV = lim

∆t→0

∫V(t ′) A(x, t ′)dV− ∫

V(t) A(x, t)dV

∆t=

= lim∆t→0

∫V(t) A(x, t ′)dV +

∫∆V A(x, t ′)dV− ∫

V(t) A(x, t)dV

∆t=

= lim∆t→0

V(t)

[A(x, t ′)−A(x, t)]∆t

dV + lim∆t→0

∆V

A(x, t ′)∆V

dV ,

posto t ′ .= t + ∆t e ∆V.= V (t ′)−V (t). In particolare quest’ultimo e (si

veda la figura 1.4):

∆V =∫

Su ·n∆tdS,

quindi

dV = u ·n∆tdS.

Poiche, inoltre:

lim∆t→0

A(x, t ′

)= A(x, t) ,

si puo scrivere:

DD t

V(t)AdV =

V

∂A∂t

dV +∫

SAu ·ndS.

8 1. LE EQUAZIONI FONDAMENTALI DELLA FLUIDODINAMICA

L’ultimo termine e chiaramente il flusso di A attraverso S. Si notiche il volume V ha perso la dipendenza dal tempo, divenendo fisso ecoincidente con V (t). Dal teorema di Green si ricava ancora:

(1.3)DD t

V(t)AdV =

V

(DA

D t+A∇ ·u

)dV .

n

u∆t

V t( )

V t+ t( )∆

S t( )

S t+ t( )∆

Figura 1.4. Moto di una particella fluida ed elementirappresentativi.

1.4.2. Equazione di bilancio di massa.La massa M del volume V rimane costante nel tempo:

M =∫

V(t)ρdV = cost ,

cioe:DM

D t= D

D t

V(t)ρdV = 0.

Dalla (1.3) si ottiene, allora:

(1.4a)∫

V

[∂ρ∂t

+∇ · (ρu)]

dV = 0

e

(1.4b)∫

V

[DρD t

+ρ∇ ·u]

dV = 0.

1.5. EQUAZIONI DI NAVIER-STOKES 9

Affinche le (1.4) siano verificate per ogni V occorre che:

(1.5a)∂ρ∂t

+∇ · (ρu) = 0

cioe, applicando la (1.2),

(1.5b)DρD t

+ρ∇ ·u = 0.

1.4.3. Equazione di bilancio della quantita di moto.Se Q e la variazione della quantita di moto del fluido contenuto nel

volume V composto sempre dalle stesse particelle, si ha:

DQ

D t= F

eDD t

V(t)ρudV =

VρfdV +

S

~~τ ·ndS,

con f forza di massa per unita di volume. Applicando la (1.3) si ricaval’equazione:

V

[D(ρu)

D t+ρu∇ ·u

]dV =

V

[ρf +∇ ·~~τ

]dV

cioe ∫

V

[u

(DρD t

+ρ∇ ·u)

+ρDu

D t

]dV =

V

[ρf +∇ ·~~τ

]dV .

Per la (1.5b), il termine in parentesi tonde al primo membro deve es-sere nullo. Affinche, quindi, la relazione trovata valga per ogni volumearbitrario, e necessario che sia:

(1.6a) ρDu

D t= ρf +∇ ·~~τ

che, per componenti, si scrive

(1.6b) ρDui

D t= ρ fi +

∂τi j

∂x j.

1.5. Equazioni di Navier-Stokes

Inserendo la 1.1 nella (1.6b) e supponendo costanti i coefficienti diviscosita si ottiene:

ρDui

D t= ρ fi +

∂∂xk

[(−p+λ

∂um

∂xm

)δik +µ

(∂ui

∂xk+

∂uk

∂xi

)]=

= ρ fi− ∂p∂xi

+(λ+µ)∂

∂xi

∂uk

∂xk+µ

∂2ui

∂xk2 .

10 1. LE EQUAZIONI FONDAMENTALI DELLA FLUIDODINAMICA

Per i liquidi incomprimibili ∇ ·u = 0, dunque si ricavano le equazionidi Navier-Stokes:

(1.7a)Dui

D t= fi− 1

ρ∂p∂xi

+ν∂2ui

∂xk2

(1.7b)Du

D t= f− 1

ρ∇p+ν∇2u ,

con ν .= µρ .

1.6. Risoluzione analitica: moti paralleli incomprimibili

Le (1.7) sono risolvibili analiticamente solo in pochi casi particolar-mente semplici. Nella maggior parte dei casi realistici, di effettivo inte-resse pratico, bisogna procedere, invece, con una trattazione numerica.Questo avviene non solo per i moti turbolenti, in cui non si possonotrascurare i termini non lineari delle equazioni di Navier-Stokes, maanche la dove tale operazione e possibile. La risoluzione numerica, conqualunque metodo venga fatta, comporta necessariamente un processodi discretizzazione del tempo e dello spazio.

Di seguito viene illustrato, a titolo esemplificativo, il piu semplicetra i casi di risoluzione analitica: il moto alla Poiseuille. Le ipotesi sonoquelle di fluido incomprimibile e di forze di massa dotate di potenzialef =−∇φ, con φ = gz. Inoltre si assume:

u1 6= 0, u2 = u3 = 0.

Dalla (1.5a) risulta, allora:

∂u1

∂x1= 0,

quindi u1 = u1(x2, x3). Proiettando la (1.7a) sui tre assi si ha:

∂u1

∂t=−g

∂∂x1

(p

ρg+z

)+ν

(∂2u1

∂x22 +

∂2u1

∂x32

)

0 =−g∂

∂x2

(p

ρg+z

)

0 =−g∂

∂x3

(p

ρg+z

).

Indicando con γ .= ρg il peso per unita di volume si definisce quotapiezometrica la quantita h

.= pγ +z, funzione unicamente di x1. La pen-

denza piezometrica e, invece, i.= − ∂h

∂x1. Il fluido si muove sempre con

pendenza i > 0.

1.6. RISOLUZIONE ANALITICA: MOTI PARALLELI INCOMPRIMIBILI 11

Se il moto avviene a velocita costante nel tempo in un piano si hache, poiche u = u(x2),

(1.8) −ig = νd2u1

dx12

e, di conseguenza, entrambi i membri devono essere costanti. Proce-dendo all’integrazione della (1.8) si ottiene:

νu1 = (−gi)x2

2

2+k1x2 +k2 ,

ove si puo notare che la distribuzione della velocita e di tipo parabolico.Se il moto avviene tra due piastre parallele indefinite distanti b l’una

dall’altra si possono imporre, ad esempio, le condizioni al contorno:

u1(x2 = 0) = u1(x2 = b) = 0

e determinare le costanti k1 e k2. La distribuzione della velocita e,dunque:

u1 =− 12µ

dpdx1

(b−x2)x2 .

CAPITOLO 2

Generalita sugli automi cellulari

Gli automi cellulari sono stati introdotti in forma pressoche com-piuta a meta del Ventesimo secolo in alcuni lavori di Johann Ludwigvon Neumann, nel tentativo di sviluppare un modello astratto di au-toriproduzione in biologia. Nel 1947, il matematico ungherese, manaturalizzato americano, propose un modello in tre dimensioni, sem-plificato, poi, nel 1951, su suggerimento di Stanislaw Ulam, in unoa sole due dimensioni. Negli ultimi cinquant’anni gli automi cellula-ri sono stati utilizzati per simulare un gran numero di processi fisicicome, ad esempio, diffusione, propagazione delle onde, meccanismi dicrescita, transizioni di fase su reticoli di spin, evoluzione di sistemitermodinamici e dinamica dei fluidi.

2.1. Definizione di automa cellulare

Un automa cellulare e una semplice idealizzazione matematica diun sistema fisico [13]. Nella sua forma piu elementare consiste in unasequenza di punti su una linea, punti a cui si possono assegnare i valori0 o 1. L’evoluzione nel tempo di questi valori e, in linea di principio,assolutamente deterministica1, secondo alcune regole concernenti uni-camente i valori dei primi vicini. Piu in generale, i punti di un automacellulare possono essere disposti su un reticolo regolare k-dimensionalee possono assumere qualunque insieme discreto di valori [12].

Da un punto di vista formale, un automa cellulare e definito [7]come un array di celle soddisfacente le seguenti proprieta:

spazio discreto: array N-dimensionale di celle;stati discreti : numero finito di bit per ogni cella;tempo discreto: ogni cella si aggiorna ad ogni “ticchettio” del-

l’orologio;aggiornamento sincrono: le celle si aggiornano tutte insie-

me;regole omogenee (isotrope): le regole di aggiornamento sono

le stesse ovunque;

1Il requisito di determinicita puo anche cadere, almeno in parte, conl’introduzione di parametri aleatori

13

14 2. GENERALITA SUGLI AUTOMI CELLULARI

regole spazialmente locali : si riferiscono solo alle celle conti-gue;

regole temporalmente locali : solo lo stato presente (al piuquello immediatamente precedente) condiziona quello succes-sivo.

Nonostante queste regole possano sembrare piuttosto restrittive, inrealta consentono la modellizzazione di una grande varieta di fenomeni.Si consideri, infatti, un reticolo di celle ciascuna con k primi vicini,compresa essa stessa (quindi k = 3, 9, 27 per reticoli rispettivamente1, 2, 3-dimensionali). Si assuma che ogni cella possa assumere di volta

in volta uno tra n possibili stati diversi. E immediato rilevare che ilvicinato puo assumere nk combinazioni diverse di stati e per ognuna diqueste la cella considerata puo assumere uno degli n stati. Sono, cosı,

possibili nnkregole. Nel caso di Life (vedi paragrafo 2.2) questo numero

e 229 ≈ 10154.La definizione di vicinato non e univoca. Si possono, infatti, ad

esempio, avere il vicinato di Neumann e quello di Moore (si veda lafigura 2.1).

NW

W

SW

NE

E

SE

N

S

(a)

W E

N

S

(b)

Figura 2.1. Tipi di vicinato: di Moore (a) e di Neumann(b).

2.2. Life

Uno dei piu noti esempi di automa cellulare e stato sviluppato dalmatematico John Horton Conway e prende il nome di Life.

Si consideri una scacchiera abbastanza grande e su di essa si dispon-ga una configurazione di pedine nere chiamate “organismi”(vedi figura2.2). Le regole proposte da Conway per l’evoluzione del sistema, oltrea soddisfare i requisiti brevemente illustrati nel paragrafo 2.1, devonoessere anche tali che [5]:

• non esista una configurazione iniziale per cui si possa banal-mente dimostrare un aumento illimitato della popolazione;

2.2. LIFE 15

• esistano, tuttavia, configurazioni iniziali che evidentemente au-mentino oltre ogni limite;

• esistano semplici configurazioni iniziali che crescono e si mo-dificano per un periodo di tempo finito, prima di estinguersicompletamente, oppure raggiungere una configurazione stabile(immutabile), oppure entrare in una fase oscillatoria che ripeteun ciclo infinito di due o piu periodi.

In altre parole, il comportamento della popolazione iniziale deve essereimprevedibile. Le regole scelte da Conway sono, infine, solo quattro:

(1) ogni organismo con due o tre vicini (considerando il vicinatodi Moore) sopravvive fino alla generazione seguente;

(2) ogni organismo con quattro o piu vicini muore (la sua pedinaviene rimossa) per sovrappopolazione;

(3) ogni organismo con un solo o nessun vicino muore per isola-mento;

(4) in ogni cella vuota con esattamente tre vicini nasce un organi-smo (vi si mette una nuova pedina).

Se si assegna un valore numerico allo stato di un organismo (ad esempio0 =morto e 1 =vivo), le regole sono di tipo totalistico. Il valore z diuna data cella dopo un passo evolutivo dipende, cioe, dalla somma Sdei valori delle celle vicine:

(1) se S= 2, 3, allora z= 1;(2) se S≥ 4, allora z= 0;(3) se S= 0, 1, allora z= 0.

Figura 2.2. Life: evoluzione di una semplice configu-razione di cinque organismi in un reticolo 10×10. Taleconfigurazione e detta “aliante” [5].

La popolazione varia continuamente, arrivando, talvolta all’estin-zione, talvolta a configurazioni stabili (“nature morte”), talvolta a con-figurazioni oscillanti all’infinito, concordemente alle richieste fatte sulleregole. In altri casi, si suppone possa crescere all’infinito, sebbene Co-nway avesse avanzato la congettura che una configurazione compostada un numero finito di organismi non possa accrescersi oltre un limitemassimo (finito). Questo assunto e stato confutato, almeno, per quantoriguarda scacchiere infinite.

Una delle configurazioni piu interessanti e quella detta ad “aliante”(vedi figura 2.2), che raggiunge nuovamente la disposizione iniziale in

16 2. GENERALITA SUGLI AUTOMI CELLULARI

quattro mosse, dopo, tuttavia, essersi spostata sulla scacchiera. Lapopolazione rimane invariata in numero, senza, pero, entrare in unciclo vero e proprio. In generale, i casi piu rilevanti sono quelli conuna popolazione iniziale di almeno quattro organismi (“tetromini”, vedifigura 2.3).

Figura 2.3. Esempio di tetromino: dopo nove mosse siraggiunge la configurazione di “semaforo”, costituita daquattro “lampeggiatori” che oscillano con periodo 2.

2.3. Classi di Wolfram

Stephen Wolfram, in un lavoro del 1983 ([13]), ha proposto di di-stinguere gli automi cellulari in quattro classi. Una prima definizione,empirica e qualitativa, di questa suddivisione e la seguente:

classe 1: lo stato iniziale evolve verso uno stato omogeneo, ca-ratterizzato dal fatto che tutte le celle hanno lo stesso valore;

classe 2: lo stato iniziale evolve verso uno stato con strutturesemplici, stabili o periodiche, separate tra loro nel reticolo;

classe 3: lo stato iniziale evolve verso uno schema completa-mente caotico;

classe 4: lo stato iniziale evolve verso uno stato con strutturecomplesse (frattali).

Un’analisi piu approfondita del grado di previdibilita dell’evolutodi un automa cellulare permette una definizione rigorosa delle classi

2.4. FORMULAZIONE MATEMATICA 17

di Wolfram, ma non verra qua portata avanti. Si noti, comunque,che per la classe 1 la predizione e assolutamente banale. Per la classe2, invece, gli effetti dovuti al particolare valore di una data cella sipropagano ad una distanza finita da essa; si ha, cioe, un’evoluzionelocale, ben delimitabile. La conoscenza di uno stato evoluto dipende,quindi, soltanto da quella dei valori iniziali di un insieme finito di celle.Al contrario, per la classe 3 il valore di una particolare cella dipende,a tempi crescenti, da quello di un sempre maggior numero di altre.In ultimo, gli automi cellulari della classe 4 sono caratterizzati da unelevato grado di imprevedibilita. Da queste premesse si puo intuire chegli automi cellulari della classe 2 agiscono come una sorta di filtro perselezionare particolari configurazioni e farle sopravvivere, a discapitodi altre che vengono annullate.

2.4. Formulazione matematica

Gli automi cellulari possono essere efficacemente rappresentati construmenti matematici. Negli ultimi anni e stato sviluppato un com-plesso formalismo matematico, essenzialmente basato sull’algebra esulla teoria dei gruppi (si veda [11]), tuttavia in questa sezione sivuole solamente presentare la semplice formulazione matematica dellecaratteristiche essenziali degli automi cellulari.

Si definisca, innanzi tutto, lo spazio delle celle:

L.= {(i, j) | i, j ∈ N, 0≤ i < n, 0≤ j < m} ,

dove i e j sono indici rispettivamente di riga e di colonna, mentre me n danno la massima estensione del reticolo. Il vicinato di Moore edefinito come:

Ni j.= {(k, l) ∈ L | |k− i| ≤ 1, |l − j| ≤ 1} ,

e in modo analogo si puo definire quello di Neumann.E evidente il modo per rappresentare le regole totalistiche: ponen-

dosi, per semplicita, nel caso di un automa cellulare unidimensionale,in cui ogni cella puo assumere valori binari (0 o 1), siano ζα un insiemedi costanti intere e numerate secondo il pedice α, tali che, considerandoun vicinato di Moore, 0≤ ζα ≤ 2. Allora il valore della cella i-esimaall’istante di tempo t +1 e:

zi (t +1) =

1, se [zi−1(t)+zi+1(t)] = ζ1

. . .

0, altrimenti .

Un’interessante restrizione delle regole totalistiche e data dalle regolelegali. Queste sono state introdotte da Wolfram per impedire che da

18 2. GENERALITA SUGLI AUTOMI CELLULARI

un reticolo completamente nello stato 0 risulti una qualunque confi-gurazione non nulla. Nel caso di un automa cellulare unidimensionalebinario con due vicini per ogni lato della cella (vicinato di Moore este-so) si hanno 32 regole legali su 1024totali. Queste sono identificate daun numero di codice, che si ricava come illustrato di seguito. Sia

zi (t +1) = f

(2r

∑j=−2r

zi+ j (t)

)

con

f : {0, 1, . . . , 2r +1}→ {0, 1} ,

allora si puo definire un codice per tutte le regole legali:

Cf.=

2r+1

∑j=0

f(zj

)2 j .

Nel caso preso in considerazione, le 32 regole legali sono individuate datutti i numeri dispari compresi tra 0 e 64.

2.5. Esempi e applicazioni

Life e poco piu di un gioco, ma gli automi cellulari hanno ancheapplicazioni di una notevole rilevanza pratica e teorica. Come verra il-lustrato nel dettaglio nei prossimi capitoli, vengono utilizzati per simu-lare il comportamento microscopico dei fluidi, in un’idealizzazione cheporta a vedere ogni singolo automa come una particella che interagiscemeccanicamente con le altre.

Un’altra nota applicazione degli automi cellulari, anche se non sem-pre viene vista sotto questo profilo, e costituita dal modello di Ising peril ferromagnetismo. In questo caso ogni automa rappresenta lo spin diun piccolo magnete, con orientazione variabile e dipendente da quelladei vicini.

In ultimo, si vuole spendere qualche parola su un problema piu ge-nerale. Si accennato, all’inizio del capitolo, al fatto che lo studio degliautomi cellulari ha avuto inizio con le ricerche sull’autoriproduzionecondotte da Johann Ludwig von Neumann. Il problema di von Neu-mann era quello di scoprire quale tipo di organizzazione logica di unautoma fosse sufficiente a produrre auto-organizzazione. La risposta futrovata in una macchina universale di Turing inserita in uno schemadi celle aventi 29 stati e cinque vicini. Un sistema di questo tipo echiamato costruttore universale, in quanto e capace di costruire mac-chine descritte dal segnale d’ingresso, ma anche a riprodurre il segnale

2.6. AUTOMI CELLULARI E EQUAZIONI DIFFERENZIALI 19

d’ingresso stesso ed il meccanismo di costruzione. E. F. Codd semplifi-co l’automa di von Neumann, introducendo una macchina richiedentesoltanto otto stati.

2.6. Automi cellulari e equazioni differenziali

Per concludere questa breve rassegna sulle caratteristiche salientidegli automi cellulari, si accennera solamente un aspetto che avra unimportante sviluppo nel seguito.

La maggior parte della fisica e basata sull’utilizzo di equazioni dif-ferenziali. Ad esempio, come si e visto nel capitolo 1, per illustrareil comportamento di un fluido in movimento e necessario ricorrere adequazioni alle derivate parziali, che hanno una soluzione analitica espli-cita solo in pochissimi casi, il piu delle volte di scarsa rilevanza pratica.E, quindi, prassi comune ricorrere a metodi di approssimazione numeri-ca. Questi constano, innanzi tutto, di una discretizzazione dello spazioe del tempo, perdendo cosı, implicitamente, il requisito di continuitache si impone usualmente ai sistemi fisici.

Per tale motivo risulta evidente l’utilita, in certe circostanze, del-l’utilizzo di automi cellulari, discreti per loro stessa natura. E questol’argomento dei prossimi capitoli, che mostreranno come e possibile il-lustrare le stesse caratteristiche viste nella derivazione classica delleNavier-Stokes mediante particolari tipi di automi cellulari. Il vantag-gio di un approccio computazionale di questo tipo, differente da quelli“ordinari” sta nella sua efficienza. Si puo, infatti, procedere alla paral-lelizzazione dei processi, a dispetto della serialita di metodi come quellibasati sulla risoluzione di sistemi lineari. In altre parole, il lavoro vienediviso tra piu elaboratori, al limite uno per ciascun automa, che, invirtu delle proprieta degli automi cellulari, ha in ogni istante sempretutte le informazioni necessarie per governare un singolo passo dell’e-voluzione, per poi scambiare dati con gli altri computer e ricominciareil processo. In tal modo vengono vantaggiosamente ridotti i tempi dicalcolo.

CAPITOLO 3

Gas reticolare

I flussi idrodinamici presentano un’ampia varieta di fenomeni com-plessi, tanto da poter essere considerati un paradigma della complessitastessa dei fenomeni naturali. Classicamente il moto di un fluido e de-scritto dalle equazioni di Navier-Stokes. Queste equazioni non linearialle derivate parziali sono del tutto indipendenti, nella loro struttura,dalla dinamica microscopica del sistema, cosı che e possibile costruire“micromondi” stilizzati che, nel limite macroscopico, riescono a descri-vere bene il comportamento di un fluido. Modelli di questo tipo vengo-no detti a gas reticolare e sono in tutto e per tutto automi cellulari, conalcune regole aggiuntive che li rendono adatti alla rappresentazione difenomeni fisici

3.1. Caratteristiche generali

Un primo approccio ai gas reticolari permette di considerarli comesemplici modelli di un fluido, caratterizzati dall’essere completamentediscreti. Viene, cioe, negato il requisito di continuita introdotto per laderivazione “classica” delle equazioni dinamiche. La discretizzazione edata dal fatto che si assumono particelle fittizie, localizzate in regionifinite di un reticolo regolare di Bravais. Queste particelle si muovono aintervalli di tempo regolari da un nodo all’altro del reticolo, con regoledi collisione locali e isotrope, cosı come assunto nel definire gli automicellulari.

L’equilibrio statistico a livello macroscopico e garantito da due tipidi condizioni. La prima riguarda le proprieta di conservazione (massa,momento, energia, . . . ) ed e necessaria, ma non sufficiente. Occorreanche che il reticolo appartenga ad un appropriato gruppo spazialecristallografico.

E altresı importante il principio di esclusione: due particelle asso-lutamente identiche non possono trovarsi simultaneamente nello stessonodo. Di conseguenza, se esistono b diversi tipi di particelle, lo stato diun dato nodo e completamente caratterizzato da b variabili booleane.

In maniera piu formale, un gas reticolare e una particolare classedi automi cellulari dove i singoli automi sono localizzati sui nodi diun reticolo di Bravais in uno spazio euclideo di dimensione D. Tale

21

22 3. GAS RETICOLARE

corrispondenza e biunivoca, cosı che d’ora in poi si chiameranno nodigli automi stessi. Il numero di possibili stati interni di ciascun nodo e2b con b intero. b e chiaramente il numero di nodi vicini, il cui statoinfluisce sull’evoluzione di quello del nodo considerato.

Il processo evolutivo e distinto in due fasi: collisione e propagazio-ne. Nella collisione lo stato di ogni nodo si modificano secondo regolelocali, che possono essere anche irreversibili e (parzialmente) non deter-ministiche. La propagazione e la fase in cui l’informazione sullo statodi un nodo si propaga a quelli vicini. Si suppone che segua sempreregole isotrope e deterministiche e vi si puo associare un vettore retico-lare ci (1≤ i ≤ b oppure 0≤ i ≤ b−1). Questo si puo pensare come lavelocita con cui una particella si muove lungo i “canali” che colleganotra loro i nodi. E, tuttavia, necessario tenere a mente che nei modellia gas di reticolo non c’e alcun moto fisico, bensı solamente un traspor-to di informazione da un nodo agli altri vicini nell’intervallo di tempo∆t, che puo sempre essere assunto unitario. Affinche sia garantito cheogni nodo sia raggiungibile a partire da un altro con una catena finitadi vettori e che, al contempo, non possano formarsi sotto-reticoli conevoluzione distinta dal resto del sistema, e necessario che tra i ci sianocompresi D generatori del reticolo.

3.2. Microdinamica

La dinamica completa del gas reticolare e completamente contenutanell’equazione microdinamica, cioe nell’espressione algebrica della rela-zione logica esatta che lega il campo booleano all’istante t +1 al campobooleano all’istante iniziale t. Formalmente:

n(t +1) = E(n(t)) ,

dove E e un operatore di evoluzione. Per il caso di un modello nondeterministico occorre considerare anche un insieme ξ di parametrialeatori:

n(t +1) = E(n(t) , ξ(t)) .

Espressioni come queste sono, tuttavia, piuttosto generiche. Perarrivare ad una forma piu esplicita si puo cominciare a notare che E

puo essere inteso come prodotto di un operatore di propagazione nonlocale P e di un operatore di collisione locale C:

E = P◦C .

In accordo con le regole di propagazione, deve essere:

(3.1) P(ni−c) = ni (r) ,

con i numero di canale.

3.3. MODELLO HPP 23

Si introduca, ora, una matrice di collisione A(s→ s′), quadrata 2b×2b, che contiene l’informazione completa sulle regole di collisione. Ognielemento di matrice rappresenta la probabilita di transizione dallo statobooleano s allo stato s′. Secondo questa interpretazione,la condizionedi normalizzazione deve essere:

∑s′

A(s→ s′

)= 1.

A questo punto l’operatore di collisione e definito da:

(3.2) C(ni)(r) = ∑(s,s′)

s′iξ(s→ s′

) b

∏j=1

n jsj n

sjj ,

con n j.= 1−n j e sj

.= 1−sj . ξ(s→ s′) e, invece, una matrice quadrata

2b×2b di valori casuali (0 o 1) che sono in media gli A(s→ s′) e taliche:

∑s′

ξ(s→ s′

)= 1.

Il prodotto ∏bj=1n j

sj nsjj e una sorta di delta di Kronecker δkr (n(r)−s).

Dalla (3.1) e dalla (3.2) si ricava l’equazione microdinamica com-pleta:

ni (r +ci , t +1) = ∑(s,s′)

s′iξ(s→ s′

) b

∏j=1

n jsj n

sjj .

Una forma alternativa e:

ni (r +ci , t +1) = ni + ∑(s,s′)

(s′i−si

)ξ(s→ s′

) b

∏j=1

n jsj n

sjj .

3.3. Modello HPP

Storicamente il primo modello e quello chiamato HPP (Hardy, de

Pazzis, Pomeau), introdotto nel 1973. E basato su un reticolo bidimen-sionale quadrato (si veda la figura 3.1). Ogni nodo ha, dunque, b = 4canali e sono possibili sedici stati interni. I vettori reticolari ci hannomodulo unitario e componenti:

c1 =(

10

), c2 =

(01

), c3 =

( −10

), c4 =

(0−1

).

Tornando a considerare le particelle fittizie, si assume che questeabbiano tutte massa unitaria, cosı che per il momento si ha pi = ci ,mentre l’energia cinetica e uguale a 1

2. Durante il processo evolutivo,poi, le particelle che collidono vengono sempre diffuse lungo la direzioneperpendicolare a quella di incidenza. In altre parole, se la collisionee nord-sud, il risultato e una propagazione est-ovest e viceversa. Si

24 3. GAS RETICOLARE

y

x

c1

c2

c3

c4

Figura 3.1. Reticolo HPP con i corrispondenti ci .

noti che solo due stati corrispondono ad una collisione effettiva, conlo stato “pre-collisionale” diverso da quello “post-collisionale”. Dunque,l’efficienza collisionale e:

∑s6=s′ A(s→ s′)∑s,s′ A(s→ s′)

= 12.5%.

La matrice A(s→ s′) e simmetrica e di dimensione 16× 16, conelementi 0 e 1, mentre l’equazione microdinamica risulta:

ni (r +ci , t +1)= ni−ni ni+1ni+2ni+3+nini+1ni+2ni+3 , per i = 1, . . . , 4.

3.4. Modello FHP-1

Un altro modello elementare, ma sufficientemente completo di gas direticolo e chiamato FHP-1 dalle iniziali di Frisch, Hasslacher e Pomeau,che l’hanno introdotto nel 1986. Il reticolo ha celle elementari di formatriangolare, cosı che ogni nodo e direttamente collegato a sei vicini (siveda la figura 3.2). Come nell’HPP i vettori reticolari ci hanno modulounitario e le loro componenti sono:

c1 =(

10

), c2 =

(12√3

2

), c3 =

(−1

2√3

2

),

c4 =( −1

0

)c5 =

(−1

2

−√

32

)c6 =

(12

−√

32

).

Per quanto riguarda l’evoluzione del sistema, mentre la propaga-zione e assolutamente analoga a quella del modello HPP, la collisioneintroduce un fattore stocastico. Come, infatti, si puo desumere dal-la figura 3.3, le leggi di conservazione permettono che due particelleprovenienti da direzioni opposte possano essere diffuse con uno stato

3.4. MODELLO FHP-1 25

1-i 4+i

3+i

2+i1+i

i

Figura 3.2. Celle elementari in un reticolo FHP-1.Sono indicati i sei canali.

finale ruotato di +60◦ o −60◦, con probabilita p e 1− p rispettivamen-te. Un’altra caratteristica che differenzia il modello FHP-1 dall’HPPe la possibilita di triple collisioni deterministiche: l’angolo tra le di-rezioni d’incidenza, cosı come quello tra le direzioni di diffusione, e120◦. Le collisioni efficaci sono, dunque, cinque su sessantaquattro,con un’efficienza del 7.81%.

c1

c4

c2

c6

c3

c5

t

c1

c4

c2

c6

c3

c5

t+1

c1

c4

c2

c6

c3

c5

Figura 3.3. [1] Collisione binaria in un reticolo FHP-1.

26 3. GAS RETICOLARE

L’equazione microdinamica mantiene la struttura formale gia vistain precedenza, solamente ampliata ai sei canali possibili:

ni (r +ci , t +1) = ni

+ni ni+1ni+2ni+3ni+4ni+5+

+ξ(+)nini+1ni+2ni+3ni+4ni+5+

+ξ(−)ni ni+1ni+2ni+3ni+4ni+5++ni ni+1ni+2ni+3ni+4ni+5++ nini+1ni+2ni+3ni+4ni+5 ,

per i = 1, . . . , 6.

ξ(±) = ξ(±) (r , t) sono insiemi di variabili aleatorie che possono assumerei valori 0 e 1 e sono tali che ξ(+) +ξ(−) = 1 ad ogni tempo t e ad ogninodo r .

3.5. Altri modelli

Il modello HPP e quello FHP non sono gli unici gas reticolari pro-posti. In questo paragrafo ci si limitera ad accennare agli altri, senzaentrare nei dettagli.

Per prime si presenteranno due varianti di FHP-1: FHP-2 e FHP-3. Il primo e sostanzialmente analogo all’originale, salvo permettere lapresenza di una particella a riposo nei nodi. Il numero dei canali e,dunque, b = 7, ma i vettori ci non sono tutti unitari. Infatti, |c0|= 0,corrispondendo c0 alla particella riposo. L’efficienza di collisione e del17.2%. Proprio allo scopo di incrementare questa grandezza e statointrodotto FHP-3: le regole di collisione sono modificate per permettereil maggior numero possibile di urti, cosı che l’efficienza e portata al59.4%.

Tutti i modelli presentati fino a questo momento sono caratterizzatidal coinvolgere particelle tutte con la stessa massa, che si differenzianotra loro solo per la velocita. Fisicamente, questo significa che si con-sidera sempre una singola specie chimica. Per studiare i fenomeni didiffusione, FHP-3 puo essere modificato per comprendere due diversespecie chimiche, costruendo CFHP (colored FHP). Nella sua versionepiu semplice, le particelle continuano ad essere meccanicamente ugualie chimicamente passive, cioe non sono possibili reazioni. Ai due diversitipi di particella sono assegnati due colori fittizi, come, ad esempio,rosso e blu. Le regole di collisione sono, poi, le stesse di FHP-3, conla differenza nel fatto che il colore e assegnato casualmente alle parti-celle uscenti, con l’unica restrizione di mantenere costante il numero diparticelle rosse e blu.

3.5. ALTRI MODELLI 27

Modelli piu complessi sono costruiti su reticoli tridimensionali. Tut-tavia si puo dimostrare che la dimensione minima per un modello vera-mente efficiente per rappresentare un sistema tridimensionale realisticoe D = 4. Un classe di modelli di questo tipo e FCHC (face centeredhyper-cubic), basato su un reticolo di Bravais definibile come quell’in-sieme di punti di R4 aventi coordinate intere con somma pari. La regoledi collisione sono numerosissime, a causa dell’elevata dimensionalita ingioco, e nella loro scelta risulta particolarmente importante considerarele condizioni di invarianza.

CAPITOLO 4

Limite macrodinamico

4.1. Distribuzione di equilibrio

Sia L un reticolo finito. Lo spazio delle fasi Γ e definito da tutti ipossibili stati s(·) = {si (r)}, dove i = 1, . . . , b e r individua un dato nodonel reticolo. Sia P(0, s(·)) una probabilita iniziale assegnata su Γ. L’e-voluzione di P e governata dall’equivalente probabilistico dell’equazionedi Liouville:

(4.1) P(t +1, Ps(·)) = ∑s′(·)

P(t, s′ (·))∏

r∈LA

(s′ (r)→ s(r)

).

Come e gia stato visto, P e un operatore di propagazione che muovela particella dal nodo r lungo il canale i-esimo. Tramite la (4.1) epossibile calcolare il valor medio di un generico osservabile, ricordandodi sommare su tutti gli stati s(·):(4.2) 〈 f (t)〉= ∑

s(·)f (s(·))P(t, s(·)) .

I valori medi calcolati con la (4.2) di alcune quantita notevoli sono:

(4.3a) Ni (t, r) = 〈ni (t, r)〉 ,

(4.3b) ρ(r) = ∑i

Ni ,

(4.3c) J(t, r) = ∑i

Nici ,

ovvero, rispettivamente, popolazione media, densita e corrente di mas-sa.

Continuando a procedere con i metodi della meccanica statistica, sidimostra [4] che puo essere introdotta un’entropia locale:

(4.4) S(r)≡−∑s(r)

P(t, s(r)) logP(t, s(r))

che non diminuisce mai in seguito a collisioni. L’entropia S e costante see solo se P e completamente fattorizzata sopra i link del sito r, ovvero:

(4.5) P =b

∏i=1

Nsii (1−Ni)

1−si

con Ni dato dalla (4.3a).

29

30 4. LIMITE MACRODINAMICO

Una soluzione t-indipendente della (4.1) e detta distribuzione diequilibrio, e indicata con Peq e si puo scrivere in forma completamentefattorizzata:

Peq(s(·)) = ∏r

P(s(r)) ,

con P(s(r)) data dalla (4.5). La corrispondente popolazione media Neqi

ha l’espressione di Fermi-Dirac:

(4.6) Neqi = [1+exp(h+q ·ci)]

−1 .

Il principio di esclusione sull’occupazione dei link fa sı che le particellepossano essere assimilate a fermioni, nonostante non sia stato definitoper esse un momento di spin, da cui il segno + davanti all’esponenziale.h e q sono i moltiplicatori di Lagrange associati alle quantita conservate.Per velocita molto inferiori a quella del suono, si puo assumere:

(4.7a) h = h0 +h2(u ·u)+ . . .

e

(4.7b) q = q1u+ . . . ,

ovvero h e q dipendono dalle potenze rispettivamente pari e dispari diu.

Per i modelli della classe generale a cui appartengono FHP e FCHCvalgono le relazioni:

(4.8) ∑i

ci = 0, ∑i

ciαciβ =(

bc2

D

)δαβ ,

dove b e sempre il numero di canali per sito, c= |ci |2 e D e la dimensionedel reticolo (D = 4 per FCHC). Sostituendo le (4.7) nella (4.6) e tenendopresente le (4.8) si ottiene:

(4.9) Neqi =

ρb

+(

ρDbc2

)uαciα +ρG(ρ)Qiαβuαuβ +O

(v3) ,

avendo posto

G(ρ) =

(D2

2bci4

)(b−2ρ)

b−ρe

(4.10) Qiαβ = ciαciβ−(

ci2

D

)δαβ .

4.2. EQUAZIONI DELLA FLUIDODINAMICA 31

4.2. Equazioni della fluidodinamica

Partendo dai risultati del paragrafo 4.1 e in particolare dall’espres-sione della popolazione media (4.9) si giungera ora a ricavare le equa-zioni di Navier-Stokes, in modo da poter completamente giustificarel’assunzione di un automa cellulare come modello microscopico di fluidoreale.

Si decomponga la popolazione media nella forma1

Ni = Neqi (ρ, u)+Nneq

i (∂ρ, ∂u) ,

dove Nneqi e il termine di non equilibrio che si assume essere propor-

zionale a O(ε). Si assume, in altre parole, che un generico sistemapossa essere costruito “perturbando” un sistema in equilibrio, per ilquale e definibile una distribuzione come la (4.5). Si e interessati allostudio dell’espansione lineare al prim’ordine in ε di Nneq

i (∂ρ, ∂u). Conconsiderazioni di simmetria si giunge dimostrare che il gradiente delladensita di massa, ∂ρ, non contribuisce e che l’unico contributo di ∂αJβa Nneq

i deriva dal prodotto con il tensore Qiαβ definito dalla (4.10). Indefinitiva, si puo scrivere:

Nneqi = ψ(ρ)Qiαβ∂αJβ +O

(ε2)

con ψ(ρ) da determinare.Si esprimano le leggi di conservazione di massa e quantita di moto

nella forma seguente:

(4.11a) ∑i

Ni (r +ci , t +1)−∑i

Ni (r , t) = 0

e

(4.11b) ∑i

Ni (r +ci , t +1)ci−∑i

Ni (r , t)ci = 0,

in accordo con l’assunzione di porre la massa totale uguale al numerodi celle occupate. Nel limite di piccoli numeri di Knudsen, ovvero con-siderando il fluido come un mezzo continuo, le differenze (4.11) possono

1Qui e nel seguito si fara uso del seguente formalismo:

∂x.= ε

∂∂x1

, ∂t.= ε

∂∂t1

+ ε2 ∂∂t2

+ . . . ,

dove ε e un parametro di piccolezza, interpretabile come un numero di Knudsenlocale:

ε v λ∇FF

,

con λ libero cammino medio di una particella e F generico campo macroscopico.Il numero di Knudsen fornisce una valutazione della rarefazione del fluido consi-derato. Il significato fisico di questa assunzione e che vari fenomeni avvengono adifferenti scale temporali e spaziali, cioe la piu piccola variazione della grandezzamacroscopica descritta da F si svolge su dimensioni comunque piu grandi del passoreticolare e su tempi piu lunghi di quelli di evoluzione del modello.

32 4. LIMITE MACRODINAMICO

essere sviluppate in serie di derivate spaziali e temporali. Eguagliandoi termini in ε e quelli in ε2 si ottiene

(4.12a) ∂tρ+∂ ·J = 0

e

(4.12b) ∂tJα +∂β

(ρG(ρ)Tαβγδuγuδ +

c2

Dδαβ

)+

+∂β

[(ψ(ρ)+

D2bc2

)Tαβγδ∂γ (ρuδ)

]= 0,

avendo posto

(4.12c) Tαβγδ = ∑i

ciαciβQiγδ

e utilizzato la (4.9). La (4.12a) e un’equazione di continuita. Nel limitedi densita costante si ha, in particolare, ∇ ·u = 0.

Per ricavare le equazioni di Navier-Stokes e essenziale l’isotropia deltensore Tαβγδ. Da proprieta di simmetria discende, dunque, che:

(4.13) Tαβγδ =bc4

D(D+2)

(δαγδβδ +δαδδβγ−

2D

δαβδγδ

),

espressione che vale sia per i modelli FHP sia per quelli FCHC, men-tre risulta inadeguata per HPP. Sostituendo la (4.13) nella (4.12b) siottengono le equazioni di Navier-Stokes nella forma:

∂tJα +∂βPαβ = ∂βSαβ ,

dove

Pαβ = ρc2

D

(1−g(ρ)

v2

c2

)δαβ +ρg(ρ)uαuβ

e

Sαβ = ν[

∂αJβ +∂βJα− 2D

(∂γJγ

)δαβ

].

Il fattore galileiano

g(ρ) =D

D+2b−2ρb−ρ

e la differenza piu notevole dalle equazioni classiche. Tuttavia, la sceltadi opportune unita di misura permette di eliminarlo, giungendo a delleequazioni di Navier-Stokes “adimensionate”. Per far questo occorreintrodurre una lunghezza e una velocita tipiche `0 e u0, per poi poterdefinire il tempo-scala tipico:

t0 =`0

g(ρ)u0.

Definendo, a questo punto, delle opportune variabili adimensionali:

T.= t

t0, R .= r

`0, U .= u

u0, Π .= P 1

g(ρ)ρ0u02

4.2. EQUAZIONI DELLA FLUIDODINAMICA 33

si ottengono, dopo opportuni passaggi algebrici, infine, le equazioni:

∇ ·U = 0

∂tU+U ·∇U =−∇Π+1Re

∇2U .

Gli operatori ∇ e ∇2 sono, ovviamente, presi rispetto a variabili spazialidi scala. In ultimo,

Re.=

u0`0g(ρ)ν(ρ)

svolge la funzione di numero di Reynolds.Si e, dunque, mostrato che le pseudo-particelle del gas reticolare

obbediscono ad equazioni del tipo di quelle di Navier-Stokes. Un fluidoreale puo, allora, essere simulato da un’assemblea di queste particelle,dopo aver fissato delle opportune regole di collisione. In altre parole sipuo utilizzare un automa cellulare.

CAPITOLO 5

Simulazioni di gas reticolari

In questo ultimo capitolo verranno presentati alcuni esempi di simu-lazioni effettuate con modelli a gas reticolari. Saranno, in particolare,studiati quei fenomeni turbolenti per cui la ricerca di soluzioni analiti-che per le equazioni di Navier-Stokes e estremamente ardua, in quantoentrano in gioco anche quei termini non lineari trascurati, ad esempio,nello studio del regime di Poiseuille.

5.1. Fluido bidimensionale attorno ad un ostacolo

Per il primo esempio di simulazione di gas reticolare si considera unfluido che si muove in una geometria a galleria del vento, attorno ad unabarra che attraversa tutto l’apparato. Lo studio viene portato al casobidimensionale, cosı che la barra si riduce ad un ostacolo rettangolareposto tra due piastre parallele.

Puo essere logico pensare che, se per rappresentare un fenomenotridimensionale e necessario un modello a quattro dimensioni, per unfluido bidimensionale si puo usare un gas a tre dimensioni. In realta,questo non e possibile perche modelli di questo tipo non soddisfano irequisiti di isotropia delle equazioni idrodinamiche. Gas bidimensio-nali sono, invece, corretti da questo punto di vista, ma pongono serielimitazioni sui coefficienti di viscosita rappresentabili. L’unica opzionepossibile rimane, quindi, quella della classe FCHC.

Figura 5.1. [8] Fluido incompressibile bidimensionaleattorno ad un ostacolo rettangolare in una geometria agalleria del vento. Simulazione effettuata con un modelloFCHC-8.

Le figure 5.1 e 5.2 illustrano i risultati ottenuti con FCHC-8. Ildominio della simulazione e un reticolo di 4096×1024nodi e ai confini

35

36 5. SIMULAZIONI DI GAS RETICOLARI

Figura 5.2. [8] Campo di vortici corrispondente alcampo di velocita in figura 5.1.

superiore e inferiore sono state imposte regole di riflessione. Il flusso haun momento di ingresso che viene annullato all’uscita, mantenendo unavelocita in linea di massima costante e con modulo ∼ 0.3. A sinistra, acento nodi dall’ingresso, c’e l’ostacolo, di dimensioni 32×100 nodi. Ilreticolo e, poi suddiviso, in celle quadrate 32×32, in cui la velocita e ladensita del fluido sono spazialmente (ma non temporalmente) mediate.L’effetto di questa media e ben visibile nei“pixel”quadrati che si notanoin figura 5.2. Le immagini rappresentano l’evoluzione del sistema dopo40000passi.

5.2. Instabilita di Kelvin-Helmholtz

Si considerino due strati fluidi che si muovono uno rispetto all’altrocon velocita parallele e separati da un’interfaccia inizialmente piana.Quando la velocita relativa diventa sufficientemente elevata, la superfi-cie di separazione diventa instabile e viene a crearsi una perturbazionesinusoidale. Di conseguenza si genera una serie di vortici “rotolan-ti” tra i due fluidi. Questo fenomeno e conosciuto come instabilita diKelvin-Helmholtz.

Flussi di questo tipo si possono bene simulare con un modello FHP-3, eventualmente permettendo anche due diverse specie di fluido, condiffusione o reazione-diffusione. In questi ultimi casi diventa ancora piuevidente la visualizzazione dei vortici. Si veda, ad esempio, la figura5.3, che mostra lo stato di una simulazione dopo 1200 passi: il fluidosuperiore si muove verso sinistra ed e piu concentrato in una dellespecie chimiche, viceversa quello inferiore, che si muove verso destra.L’interfaccia piana diventa chiaramente instabile dopo 1000 passi. Atitolo di confronto, si osservi un’istabilita di Kelvin-Helmoltz “vera” infigura 5.4.

Questo modello puo essere variato imponendo una regola di reazio-ne: se, ad un dato nodo, sono presenti in maggior numero particelledella specie A, allora le particelle della specie minoritaria B si trasfor-mano in A. In tal modo la simulazione mostra instabilita dopo unperiodo di tempo piu lungo (circa 4000 passi) e i contorni dei vorticirisultano piu nitidi (si veda la figura 5.5).

5.3. AGGREGAZIONE DI PARTICELLE 37

Figura 5.3. [8] Simulazione dell’istabilita di Kelvin-Helmholtz effettuata con un modello FHP-3 a duespecie.

Figura 5.4. [8] Fotografia di nuvole reali tra due strati atmosferici.

5.3. Aggregazione di particelle

Una semplice modifica alle regole di FHP-3 permette di studiareanche l’aggregazione di particelle solide. Per far questo si attribuisce

38 5. SIMULAZIONI DI GAS RETICOLARI

Figura 5.5. [8] Simulazione dell’instabilita di kelvin-Helmholtz effettuata con un modello FHP-3 a due speciee reazione-diffusione.

ad una particella a riposo scelta casualmente il ruolo di seme: e “ap-piccicosa”, nel senso che ogni altra particella che arriva in un sito traquelli vicini si ferma e diventa a sua volta appiccicosa. Il risultato euna struttura frattale come quella mostrata in figura 5.6, ottenuta inpochi secondi.

Figura 5.6. [8] Aggregazione frattale con un semelocalizzato al centro della figura.

Bibliografia

[1] R. Benzi, S. Succi, and M. Vergassola. The lattice Boltzmann equation: theoryand applications. Physics Reports, (222):145–197, 1992.

[2] Antonio Cenedese. Meccanica dei fluidi. McGraw-Hill, 2003.[3] D. d’Humieres, P. Lallemand, and U. Frisch. Lattice gas models for 3d

hydrodynamics. Europhysics Letters, (2:291), 1986.[4] U. Frisch, D. d’Humieres, B. Hasslacher, P. Lallemand, Y. Pomeau, and J-

P. Rivet. Lattice gas hydrodynamics in two and three dimensions. In G. D.Doolen, U. Frisch, B. Hasslacher, and S. Wolfram, editors, Lattice gas methodsfor partial differential equations. Addison-Wesley, 1990.

[5] Martin Gardner. Giochi matematici. Le Scienze, (33):102–105, Maggio 1971.[6] Martin Gardner. Giochi matematici. Le Scienze, (35):101–105, Luglio 1971.[7] David E. Hiebeler and Robert Tatar. Cellular automata and discrete physics. In

Lui Liam, editor, Introduction to nonlinear physics, pages 275–295. Springer,1997.

[8] Jean-Pierre Rivet and Jean Pierre Boon. Lattice gas hydrodynamics. CambridgeUniversity Press, 2001.

[9] S. Succi, R. Amati, and R. Benzi. Challenges in lattice Boltzmann computing.Journal of Statistical Physics, 81:5–16, 1995.

[10] John F. Wendt, editor. Computational fluid dynamics. Springer, 1996.[11] S. Wolfram, O. Martin, and A. M. Odlyzko. Algebraic properties of cellu-

lar automata. Communications in Mathematical Physics, (93):219–258, Marzo1984.

[12] Stephen Wolfram. Cellular automata as simple self-organizing systems.http://www.stephenwolfram.com/publications/articles/ca/82/cellular/, 1982.

[13] Stephen Wolfram. Cellular automata. Los Alamos Science, Fall 1983.[14] G. Zambotti. Introduzione alla meccanica statistica. Lezioni del corso di Mec-

canica Statistica I, Universita di Pavia - Dipartimento di Fisica Nucleare eTeorica, 2003.

39

Ringraziamenti

Questo lavoro e dedicato a molte persone: ai miei amici.A Massimiliano, anche se lo si vede poco.A Davide e Luca, i miei vicini, che eroicamente continuano a soppor-tarmi.Ai fagioli del Paradiso Giardino: Patrick, Mattia, Davide e Paolo.Ad Alessio.A Guglielmo, primus inter pares, per essermi sempre stato vicino, peravermi incoraggiato, spronato e sostenuto, nonche per avermi datopreziosi consigli.

Infine, voglio ringraziare anche Mamma e Papa, per non avermi maifatto mancare fiducia e affetto.

A tutti, grazie.

41