-
Notifications
You must be signed in to change notification settings - Fork 4
Expand file tree
/
Copy pathfluide_slow.py
More file actions
138 lines (108 loc) · 5 KB
/
fluide_slow.py
File metadata and controls
138 lines (108 loc) · 5 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
# Ecoulement d'un fluide
# Méthode de Boltzmann
# D'après Philip Mocz
# https://medium.com/swlh/create-your-own-lattice-boltzmann-simulation-with-python-8759e8b53b1c
# Version non optimisée mais plus lisible
import numpy as np
import matplotlib.pyplot as plt
Nx = 400 # Résolution en x
Ny = 100 # Résolution en y
Nv = 9 # Nombre de vitesses
rho0 = 100 # Densité moyenne
tau = 0.6 # Temps de collision
Nt = 50 # Nombre de pas de temps
X, Y = np.meshgrid(range(Nx), range(Ny), indexing='ij') # Grille de coordonnées (x,y)
poids = np.array([4/9,1/9,1/36,1/9,1/36,1/9,1/36,1/9,1/36]) # Somme des poids vaut 1
# Vecteurs vitesse v = (vx, vy) indexé de 0 à 8
# v = 0 : (0,0)
# v=1 : (0,1) ; v=2 : (1,1) ; v=3 : (1,0) ; v=4 : (1,-1)
# v=5 : (0,-1) ; v=6 : (-1,-1) ; v=7 : (-1,0) ; v=8 : (-1,1)
vx = np.array([0, 0, 1, 1, 1, 0,-1,-1,-1])
vy = np.array([0, 1, 1, 0,-1,-1,-1, 0, 1])
F = np.ones((Nx,Ny,Nv)) # Initialisation de la fonction de distribution à F(x,y,v) = 1 partout
np.random.seed(31415) # Graine du générateur aléatoire
F += 0.01*np.random.randn(Nx,Ny,Nv) # Ajout d'un petit bruit aléatoire
F[:,:,3] += 2 * (1+0.2*np.cos(2*np.pi*X/Nx*4)) # Ajout d'un flux vers la droite sur la composante v=3 cad v = (1,0)
rho = np.sum(F,2) # Calcul de la densité rho(x,y) = somme sur v de F(x,y,v)
for i in range(Nv):
F[:,:,i] *= rho0 / rho # On multiplie F(x,y,v) par rho0/rho(x,y) pour tout x,y
# Ainsi la densité moyenne en (x,y) est égale à rho0
# Obstacle : vaut 1 (True) si on est dans l'obstacle, 0 (False) sinon
obstacle = (X - Nx/4)**2 + (Y - Ny/2)**2 < (Ny/4)**2 # Cylindre de rayon Ny/4 centré en (Nx/4,Ny/2)
# Test de l'obstacle
# print(obstacle.shape)
# plt.scatter([0],[0], c='r', s=100)
# plt.scatter([100],[50], c='b', s=100)
# plt.imshow(obstacle.T, cmap='gray', alpha=0.3)
# ax = plt.gca()
# ax.invert_yaxis()
# plt.show()
fig = plt.figure(figsize=(10,5), dpi=80)
FF = np.zeros((Nx,Ny,Nv)) # Copie de F pour éviter les effets de bords
# Boucle principale
for it in range(Nt):
# Mouvement naturel : une mini-particule suit sont vecteur vitesse
for x in range(Nx):
for y in range(Ny):
for i in range(Nv):
FF[(x+vx[i])%Nx,(y+vy[i])%Ny,i] = F[x,y,i] # Mini-particules suivent leur vecteur vitesse (vx[i],vy[i])
# pass
# Exemple
# F[(x+0)%Nx,(y+1)%Ny,1] = F[x,y,1] # Mini-particules suivent leur vecteur vitesse (0,1)
F[:,:,:] = FF # Copie de FF dans F
# Densité rho et quantité de mvt rho u = rho (ux, uy) et vitesse moyenne u = (ux, uy)
rux = F[:,:,2] + F[:,:,3] + F[:,:,4] - F[:,:,6] - F[:,:,7] - F[:,:,8] # Quantité de mouvement en x : somme des vx*F(x,y,v)
ruy = F[:,:,1] + F[:,:,2] - F[:,:,4] - F[:,:,5] - F[:,:,6] + F[:,:,8] # Quantité de mouvement en y : somme des vy*F(x,y,v)
rho = np.sum(F,2) # Densité
ux = rux / rho # Vitesse moyenne en x
uy = ruy / rho # Vitesse moyenne en y
# Fonction d'équilibre
Feq = np.zeros((Nx,Ny,Nv))
for i in range(Nv):
Feq[:,:,i] = poids[i] * rho * (1 + 3*(ux*vx[i] + uy*vy[i]) + 9/2*(ux*vx[i] + uy*vy[i])**2 - 3/2*(ux**2 + uy**2))
# Formule principale d'itération
FF[:,:,:] = F -(1.0/tau) * (F - Feq)
# FF = F
# Bords réfléchissants
for x in range(Nx):
for y in range(Ny):
if obstacle[x,y]:
FF[x,y,0] = F[x,y,0]
FF[x,y,1] = F[x,y,5] # mini-particule allant vers la haut rebondit vers le bas
FF[x,y,2] = F[x,y,6]
FF[x,y,3] = F[x,y,7] # mini-particule allant vers la droite rebondit vers la gauche
FF[x,y,4] = F[x,y,8]
FF[x,y,5] = F[x,y,1]
FF[x,y,6] = F[x,y,2]
FF[x,y,7] = F[x,y,3]
FF[x,y,8] = F[x,y,4]
F[:,:,:] = FF
# Affichage
if (it%10 == 0) or (it == Nt-1):
print('t =',it)
plt.cla()
ux[obstacle] = 0
uy[obstacle] = 0
vortex = np.zeros((Nx,Ny))
# vortex = (np.roll(ux, -1, axis=0) - np.roll(ux, 1, axis=0)) - (np.roll(uy, -1, axis=1) - np.roll(uy, 1, axis=1))
for x in range(Nx):
for y in range(Ny):
vortex[x,y] = (uy[(x+1)%Nx,y] - uy[(x-1)%Nx,y]) - (ux[x,(y+1)%Ny] - ux[x,(y-1)%Ny])
if obstacle[x,y]:
vortex[x,y] = np.nan
# vortex = (np.roll(uy, -1, axis=0) - np.roll(uy, 1, axis=0)) - (np.roll(ux, -1, axis=1) - np.roll(ux, 1, axis=1))
# vortex[obstacle] = np.nan
# vortex[obstacle] = np.nan
vortex = np.ma.array(vortex, mask=obstacle)
# print(vortex.shape)
# print(vortex)
plt.imshow(vortex.T, cmap='bwr')
plt.imshow(~obstacle.T, cmap='gray', alpha=0.3)
plt.clim(-.1, .1)
ax = plt.gca()
ax.invert_yaxis()
# ax.get_xaxis().set_visible(False)
# ax.get_yaxis().set_visible(False)
ax.set_aspect('equal')
plt.pause(0.001)
plt.show()