http://nbviewer.ipython.org/github/barbagroup/CFDPython/blob/master/lessons/15_Step_11.ipynb
J'ai découvert comment écrire le flux de cavité, mais ce n'est pas suffisant car il ne s'anime pas.
Vous trouverez ci-dessous le code qui affiche uniquement le résultat final. (Dans l'histoire originale, iPython Notebook est utilisé)
cavorg.py
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
nx = 41
ny = 41
nt = 5
nit=50
c = 1
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)
X,Y = np.meshgrid(x,y)
rho = 1
nu = .1
dt = .001
u = np.zeros((ny, nx))
v = np.zeros((ny, nx))
p = np.zeros((ny, nx))
b = np.zeros((ny, nx))
def buildUpB(b, rho, dt, u, v, dx, dy):
b[1:-1,1:-1]=rho*(1/dt*((u[1:-1,2:]-u[1:-1,0:-2])/(2*dx)+(v[2:,1:-1]-v[0:-2,1:-1])/(2*dy))-\
((u[1:-1,2:]-u[1:-1,0:-2])/(2*dx))**2-\
2*((u[2:,1:-1]-u[0:-2,1:-1])/(2*dy)*(v[1:-1,2:]-v[1:-1,0:-2])/(2*dx))-\
((v[2:,1:-1]-v[0:-2,1:-1])/(2*dy))**2)
return b
def presPoisson(p, dx, dy, b):
pn = np.empty_like(p)
pn = p.copy()
for q in range(nit):
pn = p.copy()
p[1:-1,1:-1] = ((pn[1:-1,2:]+pn[1:-1,0:-2])*dy**2+(pn[2:,1:-1]+pn[0:-2,1:-1])*dx**2)/\
(2*(dx**2+dy**2)) -\
dx**2*dy**2/(2*(dx**2+dy**2))*b[1:-1,1:-1]
p[-1,:] =p[-2,:] ##dp/dy = 0 at y = 2
p[0,:] = p[1,:] ##dp/dy = 0 at y = 0
p[:,0]=p[:,1] ##dp/dx = 0 at x = 0
p[:,-1]=0 ##p = 0 at x = 2
return p
def cavityFlow(nt, u, v, dt, dx, dy, p, rho, nu):
un = np.empty_like(u)
vn = np.empty_like(v)
b = np.zeros((ny, nx))
for n in range(nt):
un = u.copy()
vn = v.copy()
b = buildUpB(b, rho, dt, u, v, dx, dy)
p = presPoisson(p, dx, dy, b)
u[1:-1,1:-1] = un[1:-1,1:-1]-\
un[1:-1,1:-1]*dt/dx*(un[1:-1,1:-1]-un[1:-1,0:-2])-\
vn[1:-1,1:-1]*dt/dy*(un[1:-1,1:-1]-un[0:-2,1:-1])-\
dt/(2*rho*dx)*(p[1:-1,2:]-p[1:-1,0:-2])+\
nu*(dt/dx**2*(un[1:-1,2:]-2*un[1:-1,1:-1]+un[1:-1,0:-2])+\
dt/dy**2*(un[2:,1:-1]-2*un[1:-1,1:-1]+un[0:-2,1:-1]))
v[1:-1,1:-1] = vn[1:-1,1:-1]-\
un[1:-1,1:-1]*dt/dx*(vn[1:-1,1:-1]-vn[1:-1,0:-2])-\
vn[1:-1,1:-1]*dt/dy*(vn[1:-1,1:-1]-vn[0:-2,1:-1])-\
dt/(2*rho*dy)*(p[2:,1:-1]-p[0:-2,1:-1])+\
nu*(dt/dx**2*(vn[1:-1,2:]-2*vn[1:-1,1:-1]+vn[1:-1,0:-2])+\
(dt/dy**2*(vn[2:,1:-1]-2*vn[1:-1,1:-1]+vn[0:-2,1:-1])))
u[0,:] = 0
u[:,0] = 0
u[:,-1] = 0
u[-1,:] = 1 #set velocity on cavity lid equal to 1
v[0,:] = 0
v[-1,:]=0
v[:,0] = 0
v[:,-1] = 0
return u, v, p
fig = plt.figure(figsize=(11,7), dpi=100)
u, v, p = cavityFlow(nt, u, v, dt, dx, dy, p, rho, nu)
plt.contourf(X,Y,p,alpha=0.5)
plt.colorbar()
plt.contour(X,Y,p)
plt.quiver(X[::2,::2],Y[::2,::2],u[::2,::2],v[::2,::2])
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
Il semble que vous puissiez créer une animation gif en utilisant matplotlib.animation, alors je l'ai essayé, mais cela n'a pas fonctionné. La raison sera expliquée ci-dessous.
Vous pouvez créer une animation gif si imagemagick est inclus. Voir ci-dessous pour cela. http://kiito.hatenablog.com/entry/2013/11/27/233507
La méthode spécifique pour l'animation est la suivante.
Il semble y avoir deux façons d'animer.
animation.ArtistAnimation
... Dessinez les données préparées à l'avanceanimation.FuncAnimation
... Mettre à jour les données de temps en temps
from http://cflat-inc.hatenablog.com/entry/2014/03/17/214719
Autres exemples http://qiita.com/HirofumiYashima/items/5864fb9384053b77cf6f http://d.hatena.ne.jp/xef/20120629/p2 http://stackoverflow.com/questions/24166236/add-text-to-image-animated-with-matplotlib-artistanimation
Le premier code de simulation exécute une boucle for dans cavityFlow ()
, donc j'ai pensé que je cracherais des données pour chaque boucle et dessinerais avec ArtistAnimation, mais passer contourf au tableau de données de dessin est un bogue. sortir. Cela a été rapporté ailleurs, mais je ne suis pas sûr qu'il soit pris en charge.
http://matplotlib.1069221.n5.nabble.com/Matplotlib-1-1-0-animation-vs-contour-plots-td18703.html http://stackoverflow.com/questions/23070305/how-can-i-make-an-animation-with-contourf
Pour le moment, si vous le faites avec Func Animation au lieu de Artist Animation, il semble que contourf sera également une animation, mais il y a aussi des informations que vous pouvez faire avec Artist Animation. Sauf pour contourf, j'ai essayé de dessiner uniquement des lignes de courant et il a été créé avec succès, donc cela semble être un bogue causé par contourf.
Recommended Posts