Introduzione
L'idea di base dello studio è l'analogia fra moto dei fluidi e flusso di traffico. Il traffico stradale, su larga scala (veicolo = particella), si può modellizzare come il flusso di un fluido monodimensionale comprimibile, densità e velocità che variano nello spazio e nel tempo. Con ispirazione alle note equazioni di Navier–Stokes, fluidodinamica (vedi anche: Python - Matplotlib - esempio idrodinamica, correnti fluviali 1D), arriviamo in questo caso ad un modello più pertinente, che si chiama LWR (Lighthill–Whitham 1955, Richards 1956), di fatto una legge di conservazione scalare (Navier-Stokes contiene anche un termine viscoso, descrive fluidi in cui le particelle reagiscono a ciò che sta davanti e dietro, mentre gli automobilisti reagiscono essenzialmente solo a ciò che hanno davanti quindi il modello LWR è più indicato in questo caso). Nello specifico, abbiamo:
∂ρ/∂t + ∂q/∂x = 0, con q(x,t) = ρu, dove:
- ρ(x,t) = densità [veicoli/km]
- u(x,t) = velocità media [km/h]
- q(x,t) = flusso [veicoli/h]
Occorre una relazione u = f(ρ), che definiamo tramite la cosiddetta relazione di Greenshields, lineare:
- U(ρ) = umax * (1 − ρ/ρmax)
- q(ρ) = umax * ρ * (1 − ρ/ρmax)
L'analogia è buona poiché rispetta l'intuizione strada libera = velocità libera (umax), strada satura (ρmax) = velocità nulla.
Nota: esistono anche modelli più complessi, che prevedono ad esempio l'introduzione di un termine di secondo ordine per la velocità (Payne–Whitham, Aw–Rascle–Zhang (ARZ)). Nel caso in esame continuiamo a trattare il modello LWR.
Codice Python
Vediamo il seguente codice Python, scritto con l'aiuto dell'IA (Claude Opus 5). Viene applicato il modello numerico Godunov (upwind conservativo, adatto alle leggi di conservazione non lineari perché gestisce correttamente gli urti), con dominio periodico, analogia anello autostradale.
import numpy as np
import matplotlib.pyplot as plt
# ---------------- Parametri fisici ----------------
L = 10.0 # lunghezza tratto [km] (anello)
rho_max = 180.0 # densità di ingorgo [veic/km]
u_max = 100.0 # velocità a flusso libero [km/h]
T_fin = 0.35 # tempo simulato [h] ~ 21 min
N = 400 # celle
dx = L / N
x = (np.arange(N) + 0.5) * dx
CFL = 0.9
dt = CFL * dx / u_max # condizione CFL: max|q'(rho)| = u_max
nt = int(T_fin / dt)
# ---------------- Modello (Greenshields) ----------------
def U(rho): # velocità media
return u_max * (1.0 - rho / rho_max)
def q(rho): # flusso, concavo, max in rho_c = rho_max/2
return rho * U(rho)
rho_c = rho_max / 2.0
def godunov_flux(rL, rR):
"""Flusso numerico di Godunov per flusso CONCAVO q(rho)."""
F = np.where(rL <= rR,
np.minimum(q(rL), q(rR)), # onda di rarefazione/urto
np.maximum(q(rL), q(rR))) # urto
# se il massimo di q cade dentro l'intervallo (caso rL > rR)
sonic = (rL > rR) & (rR <= rho_c) & (rho_c <= rL)
return np.where(sonic, q(rho_c), F)
# ---------------- Condizione iniziale: ingorgo localizzato ----------------
rho0 = 20.0 + 130.0 * np.exp(-((x - 3.0) ** 2) / (2 * 0.4 ** 2))
rho = rho0.copy()
RHO = np.zeros((nt + 1, N))
RHO[0] = rho
# ---------------- Integrazione ----------------
for k in range(1, nt + 1):
rL = rho # stato a sinistra di ogni interfaccia
rR = np.roll(rho, -1) # stato a destra (dominio periodico)
F = godunov_flux(rL, rR) # F[i] = flusso all'interfaccia i+1/2
rho = rho - dt / dx * (F - np.roll(F, 1))
rho = np.clip(rho, 0.0, rho_max) # solo pulizia numerica
RHO[k] = rho
VEL = U(RHO)
# ---------------- Controllo di conservazione ----------------
n0, n1 = RHO[0].sum() * dx, RHO[-1].sum() * dx
print(f"dt = {dt*3600:.2f} s, passi = {nt}")
print(f"Veicoli iniziali: {n0:.2f} | finali: {n1:.2f} | errore relativo: {abs(n1-n0)/n0:.2e}")
# ---------------- Diagrammi spazio-tempo ----------------
fig, ax = plt.subplots(1, 2, figsize=(13, 5), sharey=True)
ext = [0, L, T_fin * 60, 0] # tempo in minuti, t=0 in alto
im0 = ax[0].imshow(RHO, cmap='Reds', aspect='auto', extent=ext,
vmin=0, vmax=rho_max)
ax[0].set_title("Densità ρ(x,t)")
ax[0].set_xlabel("spazio [km]"); ax[0].set_ylabel("tempo [min]")
fig.colorbar(im0, ax=ax[0], label="veic/km")
im1 = ax[1].imshow(VEL, cmap='Greys_r', aspect='auto', extent=ext,
vmin=0, vmax=u_max)
ax[1].set_title("Velocità media u(x,t)")
ax[1].set_xlabel("spazio [km]")
fig.colorbar(im1, ax=ax[1], label="km/h")
plt.suptitle("Modello LWR – diagramma spazio-tempo (schema di Godunov)")
plt.tight_layout()
plt.show()
# ---------------- Diagramma fondamentale ----------------
r = np.linspace(0, rho_max, 200)
plt.figure(figsize=(6, 4))
plt.plot(r, q(r))
plt.axvline(rho_c, ls='--', lw=1)
plt.title("Diagramma fondamentale di Greenshields")
plt.xlabel("ρ [veic/km]"); plt.ylabel("q [veic/h]")
plt.tight_layout(); plt.show()
Risultati e considerazioni finali
I risultati della simulazione sono:
- dt = 0,81 s
- numero passi = 1555
- veicoli iniziali: = 330,34 (nota: numero decimale, dipende poi dalla scala usata dal modello, non è da intendersi 1 veicolo = 1 unità)
- veicoli finali = 330,34
- errore relativo = 3,44e-16
Il modello è valido per il caso in esame, poiché:
- l'ingorgo si propaga all'indietro (concettualmente ha senso, le auto vanno avanti, l'ingorgo va indietro, contrario al verso di marcia)
- urto a monte, rarefazione a valle: l'ingorgo si deforma in modo asimmetrico (un metodo numerico conservativo come Godunov è necessario in questo caso, il fronte di ingresso in coda resta netto mentre l'uscita dalla coda è sfumata)
- la velocità è l'immagine speculare della densità (questo è un modello semplice, per modellazioni più complesse si fa uso di altri modelli, come già accennato - es. formazione della coda e smaltimento, le velocità reali sarebbero diverse - come ipotesi semplificativa comunque va bene)
- verifica di conservazione (la verifica di conservazione del metodo numerico mostra un errore trascurabile (dell'ordine di 1e-16)
Vediamo infine i risultati nell'immagine che segue.
