PYTHON   25
ecoulement poiseuille
Guest on 14th March 2023 01:12:27 PM


  1. import matplotlib.pyplot as plt
  2. import matplotlib.animation as animation
  3. from mpl_toolkits.mplot3d import Axes3D
  4. import numpy as np
  5. import time
  6. import random as rn
  7.  
  8. class Application:
  9.     def __init__(self):
  10.         self.R = 1
  11.         self.p1 = 1.5
  12.         self.p2 = 1
  13.         self.eta = 0.1
  14.         self.L = 2
  15.         self.to = time.time()
  16.         self.t = time.time()
  17.  
  18.         self.scalex = 1
  19.         self.scaley = 1
  20.         self.scalez = 1
  21.         self.fig = plt.figure(figsize=(10,5))
  22.         self.ax = self.fig.add_subplot(111, projection='3d')
  23.         plt.subplots_adjust(left=0, right=1, top=1, bottom=0)
  24.  
  25.         self.plot = [self.ax.plot([0], [0], [0])]
  26.         self.create_boundaries()
  27.         self.scat = self.ax.scatter([],[],[], s=30)
  28.         self.man = ParticuleManager(self.R, self.L, self.p2, self.p1, self.eta)
  29.  
  30.     def create_boundaries(self):
  31.         i=0
  32.         while i <= 2:
  33.             self.plot[0] = self.ax.plot3D([0, self.L], [np.cos(np.pi*i), np.cos(np.pi*i)], [np.sin(np.pi*i), np.sin(np.pi*i)], c=[0,0,0,0.2])
  34.             i = i + 0.25
  35.         i = 0
  36.  
  37.         circle_angle = np.linspace(0, 2, 30)
  38.         circle_y = np.linspace(0, 0, 30)
  39.         circle_z = np.linspace(0, 0, 30)
  40.         while i < 30:
  41.             circle_y[i] = np.cos(circle_angle[i]*np.pi)
  42.             circle_z[i] = np.sin(circle_angle[i]*np.pi)
  43.             i = i + 1
  44.         i = 0
  45.  
  46.         while i <= self.L:
  47.             self.plot[0] = self.ax.plot3D(np.linspace(i, i, 30), circle_y, circle_z,  c=[0,0,0,0.2])
  48.             i = i + 0.5
  49.         self.ax.set_xlim(0, self.L)
  50.         self.ax.set_ylim(-self.R, self.R)
  51.         self.ax.set_zlim(-self.R, self.R)
  52.  
  53.  
  54.     def update(self, i):
  55.         t = time.time() - self.t
  56.         self.t = time.time()
  57.         self.man.move(t, self.p2, self.p1, self.eta, self.L, self.R)
  58.         self.man.spawn(self.L)
  59.         self.man.render(self.scat)
  60.  
  61. class Particule:
  62.     def __init__(self, r=1, an=0, x=0, col=0, v=0):
  63.         self.r = r
  64.         self.an = rn.uniform(0,2)
  65.         self.x = x
  66.         self.y = r*np.cos(self.an*np.pi)
  67.         self.z = r*np.sin(self.an*np.pi)
  68.         self.v = v
  69.         self.col = col
  70.  
  71. class ParticuleManager:
  72.     def __init__(self, R, L, p2, p1, eta):
  73.         self.time = time.time()
  74.         self.nb = 700
  75.         self.particules = []
  76.         for i in range(self.nb):
  77.             r=rn.uniform(0,R)
  78.             x=rn.uniform(0,L)
  79.             v=(p2-p1)*(r*r-R*R)/4/L/eta
  80.             self.particules.append(Particule(r=r, x=x, v=v))
  81.         #self.particules = [Particule(r=rn.uniform(0,R), x=rn.uniform(0,L)) for i in range(self.nb)]
  82.         self.positionsx = ([0] * self.nb)
  83.         self.positionsy = ([0] * self.nb)
  84.         self.positionsz = ([0] * self.nb)
  85.  
  86.         for i in range(self.nb):
  87.             self.positionsx[i] = self.particules[i].x
  88.             self.positionsy[i] = self.particules[i].y
  89.             self.positionsz[i] = self.particules[i].z
  90.  
  91.     def move(self, t, p2, p1, eta, L, R):
  92.         i=0
  93.         for p in self.particules:
  94.             p.x = p.x + t*p.v
  95.             self.positionsx[i]= p.x
  96.             i = i+1
  97.  
  98.     def spawn(self, L):
  99.         for i in range(self.nb):
  100.             if self.particules[i].x > L:
  101.                 self.particules[i].x = 0
  102.                 self.positionsx[i] = self.particules[i].x
  103.  
  104.  
  105.     def render(self, scat):
  106.         scat._offsets3d = (self.positionsx, self.positionsy, self.positionsz)
  107. a = Application()
  108.  
  109. ani = animation.FuncAnimation(a.fig, a.update, interval=40, blit=False)
  110. plt.show()
Raw Paste

Login or Register to edit or fork this paste. It's free.