Il Lattice Boltzmann Method (LBM) è una tecnica potente e versatile per la simulazione numerica dei fluidi. Si distingue per la sua capacità di gestire complesse condizioni al contorno e per la relativa facilità di implementazione, rispetto ad altri metodi computazionali. Esploreremo, brevemente, la metodologia LBM attraverso una specifica implementazione: la simulazione del flusso bidimensionale attorno a un cilindro.
Il cuore del metodo LBM è l’equazione di Boltzmann discretizzata nello spazio, nel tempo e nelle velocità. In particolare, il modello D2Q9, che utilizza nove direzioni di velocità in uno spazio bidimensionale, è una scelta comune per via della sua semplicità e accuratezza. L’equazione di Boltzmann discretizzata è implementata per evolvere la funzione di distribuzione delle particelle nel reticolo.
import numpy as np
from numpy import fromfunction, roll
from numpy.linalg import norm
import matplotlib.pyplot as plt
from matplotlib import cm
import imageio
import os
# Parametri del problema
maxIter = 200000 # Numero totale di iterazioni temporali
Re = 220.0 # Numero di Reynolds
nx, ny = 520, 180 # Dimensioni della griglia
ly = ny - 1.0
q = 9 # Numero di popolazioni in ciascuna direzione
cx, cy, r = nx // 4, ny // 2, ny // 9 # Coordinate del cilindro
uLB = 0.04 # Velocità in unità di griglia
nulb = uLB * r / Re # Viscosità in unità di griglia
omega = 1.0 / (3.0 * nulb + 0.5) # Parametro di rilassamento
# Costanti del reticolo
c = np.array([(x, y) for x in [0, -1, 1] for y in [0, -1, 1]]) # Velocità del reticolo
t = 1.0 / 36.0 * np.ones(q) # Pesi del reticolo
t[np.asarray([norm(ci) < 1.1 for ci in c])] = 1.0 / 9.0
t[0] = 4.0 / 9.0
noslip = [c.tolist().index((-c[i]).tolist()) for i in range(q)] # Indice delle condizioni di non scorrimento
i1 = np.arange(q)[np.asarray([ci[0] < 0 for ci in c])] # Parete destra
i2 = np.arange(q)[np.asarray([ci[0] == 0 for ci in c])] # Parete centrale verticale
i3 = np.arange(q)[np.asarray([ci[0] > 0 for ci in c])] # Parete sinistra
# Funzioni di supporto
def sumpop(fin):
return np.sum(fin, axis=0)
def equilibrium(rho, u): # Funzione di equilibrio
cu = 3.0 * np.dot(c, u.transpose(1, 0, 2))
usqr = 3.0 / 2.0 * (u[0] ** 2 + u[1] ** 2)
feq = np.zeros((q, nx, ny))
for i in range(q):
feq[i, :, :] = rho * t[i] * (1.0 + cu[i] + 0.5 * cu[i] ** 2 - usqr)
return feq
# Inizializzazione
obstacle = fromfunction(lambda x, y: (x - cx) ** 2 + (y - cy) ** 2 < r ** 2, (nx, ny))
vel = fromfunction(lambda d, x, y: (1 - d) * uLB * (1.0 + 1e-4 * np.sin(y / ly * 2 * np.pi)), (2, nx, ny))
feq = equilibrium(1.0, vel)
fin = feq.copy()
# Percorsi delle immagini e della GIF
image_paths = []
# Loop principale
for time in range(maxIter):
# Condizione di uscita
fin[i1, -1, :] = fin[i1, -2, :]
rho = sumpop(fin)
u = np.dot(c.transpose(), fin.transpose((1, 0, 2))) / rho
# Parete sinistra: calcolo della densità
u[:, 0, :] = vel[:, 0, :]
rho[0, :] = 1.0 / (1.0 - u[0, 0, :]) * (sumpop(fin[i2, 0, :]) + 2.0 * sumpop(fin[i1, 0, :]))
# Condizione di Zou/He
feq = equilibrium(rho, u)
fin[i3, 0, :] = fin[i1, 0, :] + feq[i3, 0, :] - fin[i1, 0, :]
# Collisione
fout = fin - omega * (fin - feq)
for i in range(q):
fout[i, obstacle] = fin[noslip[i], obstacle]
# Streaming
for i in range(q):
fin[i, :, :] = roll(roll(fout[i, :, :], c[i, 0], axis=0), c[i, 1], axis=1)
# Visualizzazione e salvataggio immagini
if time % 100 == 0:
plt.clf()
img_path = "vel." + str(time // 100).zfill(4) + ".png"
plt.imshow(np.sqrt(u[0] ** 2 + u[1] ** 2).transpose(), cmap=cm.Reds)
plt.savefig(img_path)
image_paths.append(img_path)
# Creazione della GIF
with imageio.get_writer('flow_simulation.gif', mode='I') as writer:
for img_path in image_paths:
image = imageio.imread(img_path)
writer.append_data(image)
# Rimozione delle immagini
for img_path in image_paths:
os.remove(img_path)
Il codice Python fornito implementa la simulazione LBM per il flusso attorno a un cilindro in 2D. Esaminiamo il codice passo dopo passo per vedere come i concetti di LBM sono tradotti in programmazione.
Inizializzazione
Il codice inizia con l’importazione dei moduli necessari e la definizione dei parametri di simulazione:
from numpy import *; from numpy.linalg import * import matplotlib.pyplot as plt; from matplotlib import cm maxIter = 200000 # Numero totale di iterazioni temporali. Re = 220.0 # Numero di Reynolds.
La scelta di maxIter e Re determina rispettivamente la durata della simulazione e il regime di flusso. Il numero di Reynolds è particolarmente critico poiché definisce la transizione tra flusso laminare e turbolento.
Definizione del Reticolo e Parametri di Flusso
nx = 520; ny = 180; ... # Dimensioni del reticolo e popolazioni. uLB = 0.04 # Velocità in unità di reticolo. nulb = uLB*r/Re; omega = 1.0 / (3.*nulb+0.5); # Parametro di rilassamento.
La definizione delle dimensioni del reticolo (nx, ny) e del parametro uLB stabilisce il quadro della simulazione. Il parametro di rilassamento omega è derivato dalla viscosità e gioca un ruolo cruciale nell’equazione di evoluzione.
Configurazione del Reticolo e Condizioni al Contorno
Il cilindro viene modellato e le condizioni al contorno vengono stabilite:
obstacle = fromfunction(lambda x,y: (x-cx)**2+(y-cy)**2<r**2, (nx,ny)) vel = fromfunction(lambda d,x,y: (1-d)*uLB*(1.0+1e-4*sin(y/ly*2*pi)), (2,nx,ny)) feq = equilibrium(1.0, vel); fin = feq.copy()
obstacle rappresenta la posizione del cilindro nel reticolo, mentre vel stabilisce il profilo di velocità in ingresso, fondamentale per l’innesco del flusso.
Ciclo Principale di Simulazione
Il nucleo della simulazione è un ciclo che aggiorna le funzioni di distribuzione secondo l’equazione di Lattice Boltzmann:
for time in range(maxIter):
... # Condizioni al contorno
... # Passo di collisione
... # Passo di streaming
Qui, la dinamica del flusso è simulata con passi sequenziali che includono la gestione delle condizioni al contorno (come l’uscita del flusso e il movimento del muro), la collisione (aggiornamento delle funzioni di distribuzione verso l’equilibrio), e lo streaming (spostamento delle funzioni di distribuzione attraverso il reticolo).
L’implementazione fornita illustra chiaramente come i principi teorici di LBM siano incorporati in un contesto di simulazione pratica. Questa specifica simulazione LBM ci fornisce un esempio concreto di come metodi computazionali avanzati possano essere utilizzati per investigare fenomeni fisici complessi in modo relativamente semplice e intuitivo. Attraverso un’adeguata sintonizzazione dei parametri e delle condizioni iniziali, il metodo Lattice Boltzmann si rivela uno strumento potentissimo nell’analisi e nella visualizzazione dei flussi fluidi, sia in contesti accademici che industriali.
.
della funzione 
.
come esempio.
da 
a 
.
.
, without change of shape.
. Then the exact solution of the equation is 
. 
![x_i \in [0,2]](https://s0.wp.com/latex.php?latex=x_i+%5Cin+%5B0%2C2%5D&bg=ffffff&fg=000&s=0&c=20201002)
. We have define a variable 
is given as 
in the interval 
and 
everywhere else in 
(i.e., a hat function). Here, we use the function 
for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.
timestep, using once again the NumPy function 