1 Flusso di calore.
In questo tutorial scriveremo un codice per risolvere una equazione differenziale alle derivate parziali di tipo parabolico, come quella che descrive l'andamento nello spazio e nel tempo della temperatura di un certo sistema unidemensionale:
Come visto nelle dispense (pag.) il metodo FTCS si ottiene discretizzando le derivate della precedente equazione e riordinando i termini si arriva alla seguente relazione:
Come si nota, è quindi possibile trovare la temperatura in tutti i punti del nostro sistema unidimensionale ad un certo istante , conoscendo la temperatura al tempo precedente nello stesso punto e in quelli adiacenti . È perciò necessaria la conoscenza della temperatura in certo istante iniziale, per calcolare la sua evoluzione ad istanti successivi e dobbiamo inoltre fissare le condizioni al contorno, quindi decidere cosa accade alle estremità del nostro sistema. L'implementazione dello schema FTCS consiste nel tradurre in codice la relazione precedente e iterarla all'interno di un ciclo for un certo numero di volte. Vediamo il codice di seguito:
1 2 3 | for time in range(1,tstep):
for i in range(1,N-1):
T[i,time]=T[i,time-1]+0.5*tau/ts*(T[i-1,time-1]+T[i+1,time-1]-2*T[i,time-1])
|
- la temperatura del sistema ai vari istanti è memorizzata nella matrice T, dove l'indice di riga i rappresenta la posizione mentre l'indice di colonna time rappresenta l'istante temporale;
- abbiamo due cicli for: uno sul tempo che parte da 1 quindi dall'istante successivo all'istante iniziale (per il quale conosciamo già T); l'altro sullo spazio che va dalla posizione 1 alla N-1, quindi esclude le posizioni ai bordi che stiamo cosi fissando costanti al valore iniziale, implementando di fatto le condizioni di Dirichlet.
Al solito, per poter eseguire correttamente le righe precedenti dobbiamo definire le condizioni iniziali:
1 2 3 4 5 6 7 8 9 | N=100
tstep=1000
T=zeros((N,tstep))
#Temperatura iniziale
T[50,0]=10.
ts=1
tau=0.9*ts
|
Fissiamo quindi le dimensioni del sistema con N e la durata temporale tstep della nostra simulazione, che risultano anche essere le dimensioni della matrice T. Decidiamo anche la costante ts e l'intervallo temporale tau in maniera che esso sia leggermente piu piccolo, per avere la stabilità dello schema di integrazione.
La temperatura iniziale è stata posta a 10 in un singolo punto al centro del sistema, simulando cosi uno spike di temperatura. Tutto il resto del sistema è a temperatura nulla.
Siamo ora pronti per eseguire il nostro programma e vedere come evolve la temperatura del nostro sistema. Una volta eseguito dovremo fare un grafico della temperatura in funzione del tempo e dello spazio perciò necessitiamo di un grafico 3D. Per otternere questo ci servono le seguenti linee di codice:
1 2 3 4 5 | from mpl_toolkits.mplot3d import Axes3D
ax = gca(projection='3d')
gridx , gridy = meshgrid(range(tstep),range(N))
ax.plot_wireframe(gridx,gridy,T,cstride=10,rstride=2)
|
- importiamo il tool Axes3D per poter disegnare grafici in 3D;
- creiamo la figura con i 3 assi X,Y,Z;
- mediante meshgrid() creiamo due matrici gridx e gridy che ci danno le coordinate x,y di tutti i punti di cui vogliamo rapresentare la coordinata z; in input alla funzione diamo due sequenze di lunghezza pari al numero di righe e colonne, N e tstep nel nostro caso, che rappresentano i punti negli assi x e y.
- infine, disegnamo il grafico mediante plot_wireframe() dando in input le matrici gridx, gridy e T. Se la matrice da rappresentare è grande come nel nostro caso, è conveniente limitare il numero di punti realmente disegnati usando cstride e rstride che fissano ogni quante x e y disegnare.
Per regolare meglio il range di valori sia in x che in y, quindi posizione e tempo nel nostro caso, e avere un grafico più comprensibile, ho usato le seguenti righe:
1 2 3 4 5 | from mpl_toolkits.mplot3d import Axes3D
ax = gca(projection='3d')
gridx , gridy = meshgrid(range(2,200),range(0,N,2))
ax.plot_wireframe(gridx,gridy,Temp[::2,2:200])
|
Con range(2,200) ho creato selezionato un range temporale che mi esclude i due istanti iniziali e che considera solo i successivi 200 istanti. Con range(0,N,2) disegno un solo punto su due lungo l'asse x.
Si nota come il picco di temperatura iniziale si abbassa e si allarga lungo x con il passare del tempo, segno che il calore inialmente generato al centro del sistema, diffonda in entrambe i lati ed fuoriesca dal sistema ai bordi (condizioni di Dilichlet). Perciò, per tempi molto lunghi, la temperatura torna a zero su tutto il sistema.
Proviamo ora ad implementare le condizioni di Neumann, per le quali la derivata della temperatura rispetto a x negli estremi e' nulla:
Queste condizioni impongono un flusso di calore nullo agli estremi, per cui se supponiamo di avere un certa quantità di calore all'istante iniziale concentrata in una zona del sistema in esame, questa si distribuisce sino a che la temperatura raggiunge un valore costante per tutto il sistema. Per implementare queste condizioni al contorno dobbiamo imporre che, aggiungendo queste due righe al codice precedente:
1 2 3 4 5 | for time in range(1,tstep):
#ciclo for su i
...
T[0,time]=T[1,time]
T[99,time]=T[98,time]
|
Quindi una volta calcolata la temperatura su tutti i punti interni, aggiorniamo gli estremi imponendo la temperatura uguale a quella dei loro punti adiacenti.
Ora vediamo cosa succede se imponiamo una temperatura iniziale più alta in una zona centrale del sistema lunga 30 punti e lasciando andare la simulazione per 3000 passi:
1 2 3 4 | N=100
tstep=3000
Temp=zeros((N,tstep))
Temp[35:65,0]=500.
|
Come vediamo il calore iniziale si distribuisce aumentanto la temperatura delle zone inizialmente a zero gradi, sino a che questa non diventa uniforme per tutto il sistema.
1 | ax.plot_surface(gridx,gridy,Temp,cmap=cm.coolwarm,vmax=250,linewidth=0,rstride=2, cstride=100)
|
- cmap=cm.coolwarm, permette di aggiungere un colore alle singole caselle della superficie in base al valore della loro coordinata z; la particolare colorazione è data da cm.collwarm, una delle moltissime palette di colori a disposizione;
- vmax=250, il valore massimo sopra il quale il colore non cambia;
- linewidth=0, le linee di confine delle caselle della superficie hanno spessore 0 (appaiono in bianco);
- cstride e rstride, limitano i punti realmente disegnati e quindi la grandezza delle caselle.