更新空いてしまいました。
Parabolic Bowl
先週末から平面二次元の計算を始めました。
手始めにparabolic bowlの図化をpythonで書いてみました。
気が向いたら回してみてください。DBに置いておきます。
さすがに遅すぎるので実計算には使えそうにないです。動画っぽくしたいなら、ImageMagickでgifにするほうが良いかもしれないです。
ソースコードはこんな感じです。ほぼグラフのところで、解析解の計算は数行です。やっぱりI/Oは手間ですね。
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
class MyGragh:
def __init__(self,X,Y,zb,wL,idX,idY):
self.X,self.Y,self.idX,self.idY = X,Y,idX,idY
self.fig = plt.figure(figsize=(12, 10))
self.ax = self.fig.add_subplot(2, 1, 1, projection='3d')
self.ax.set_zlim3d(0, 14)
self.ax.view_init(elev=33, azim=38)
Xm, Ym = np.meshgrid(self.X, self.Y)
h = list(map(list,zip(*zb)))
self.ax.plot_surface(Xm, Ym, h, color='k', zorder=1, rstride=1, cstride=1,linewidth=0.1, antialiased=False,alpha=0.3)
h = list(map(list,zip(*wL)))
self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1)
self.ax2 = self.fig.add_subplot(2, 2, 3)
plt.title('X = '+str(self.X[self.idX]))
h = zb[self.idX,:]
plt.plot(self.Y,h, color='k')
h = wL[self.idX,:]
self.lineX = self.ax2.plot(self.Y,h, color='b')
self.ax3 = self.fig.add_subplot(2, 2, 4)
plt.title('Y = '+str(self.Y[self.idY]))
h = zb[:,self.idY]
plt.plot(self.X,h, color='k')
h = wL[:,self.idX]
self.lineY = self.ax3.plot(self.X,h, color='b')
def plotWL(self,wL,t,u,v):
self.surf.remove()
Xm, Ym = np.meshgrid(self.X, self.Y)
h = list(map(list,zip(*wL)))
self.surf = self.ax.plot_surface(Xm, Ym, h, cmap=cm.jet, zorder=2, rstride=1, cstride=1,linewidth=0.1, antialiased=False, alpha=1)
for c in self.lineX : c.remove()
h = wL[self.idX,]
self.lineX = self.ax2.plot(self.Y,h, color='b')
for c in self.lineY : c.remove()
h = wL[:,self.idY]
self.lineY = self.ax3.plot(self.X,h, color='b')
self.fig.suptitle('t='+str(t)+'[sec]'+' U='+str(round(u,3))+'[m/s]'+' V='+str(round(v,3))+'[m/s]')
plt.pause(0.001)
def __enter__(self):
return self
def __exit__(self, exc_type, exc_value, traceback):
plt.savefig('output.png')
return False
nx, ny = 21, 21
X = np.arange(-10, 11, 1)
Y = np.arange(-10, 11, 1)
a = 7
h0 = 4
zb = np.array([[h0/a**2*(X[i]**2+Y[j]**2) for j in range(ny)] for i in range(nx)])
eta = a/3
tPriod=60
omega = 2*np.pi/tPriod
wL = lambda t, x, y: eta*h0/a**2*(2*x*np.cos(t*omega)+2*y*np.sin(t*omega)-eta)+h0
wLr = lambda t: np.array([[ wL(t,X[i],Y[j]) if wL(t,X[i],Y[j]) - zb[i,j] > 0 else zb[i,j] for j in range(ny)] for i in range(nx)])
u = lambda t: -eta*omega*np.sin(t*omega)
v = lambda t: -eta*omega*np.cos(t*omega)
idX=10
idY=10
tmax = 600
with MyGragh(X,Y,zb,wLr(0),idX,idY) as g:
for t in range(1,tmax):
if t%10==0 : g.plotWL(wLr(t),t,u(t),v(t))
とりあえず、矩形断面ではできたっぽくて、テンションが上がってます。パッと見せれる絵がないので勉強会のときにでも話します。
