Skip to content

Instantly share code, notes, and snippets.

@dekosuke
Created November 6, 2011 13:28
Show Gist options
  • Save dekosuke/1342876 to your computer and use it in GitHub Desktop.
Save dekosuke/1342876 to your computer and use it in GitHub Desktop.
流体描いてみたけどカルマン渦にならなかったの巻
#http://ssrs.dpri.kyoto-u.ac.jp/~nakamichi/exercise/presentation.pdf
#を実装したけど渦にならない・・・
datadir = "YOUR SAVE DIR PATH"
import random,copy
def createMesh(xs,ys):
return [[0 for y in range(ys)] for x in range(xs)]
X = 400
Y = 100
Nu = 0.02
U=1
dt = 0.001
dx = 0.1
omega = createMesh(X,Y)
psi = createMesh(X,Y)
def isSpot(i,j):
if 45<i<55 and 45<j<55: return True
else: return False
def isSpot2(i,j):
if 45<i<55 and 45<j<55: return True
else: return False
def isNearSpot(i,j):
if 45<=i<=55 and 45<=j<=55: return True
else: return False
#initial condition
for i in range(0,X):
for j in range(0,Y):
#if isSpot(i,j):continue
psi[i][j]=j*U*dx+U
if 10<i<20 and 10<j<Y-10 and not isSpot(i,j):
psi[i][j]+=(random.random()-0.5)*0.0000001
pass
for i in range(1,X-1):
for j in range(1,Y-1):
omega[i][j]=-(psi[i+1][j]+psi[i][j+1]+psi[i-1][j]+psi[i][j-1]-psi[i][j]*4)/(dx*dx)
def getPsi(i,j):
global psi
if isSpot(i,j):return 0
return psi[i][j]
def u(i,j):
global U,omega,psi
if j==1 or j==Y-2 or i<=3:return U
if i==X-1:return 0
else:return (getPsi(i,j+1)-getPsi(i,j-1))/(2*dx)
def v(i,j):
global U, omega, psi
if j==1 or i==1 or i==X-2 or j==Y-2:
return 0
else: return -(getPsi(i+1,j)-getPsi(i-1,j))/(2*dx)
def getOmega(i,j):
global omega
if isSpot(i,j): return 0
return omega[i][j]
def dxw(i,j):
#if isSpot(i,j): return 0
if j==1 or j==Y-2: return 0
return (getOmega(i+1,j)-getOmega(i-1,j))/(2*dx)
def dyw(i,j):
#if isSpot(i,j): return 0
if j==1 or j==Y-2: return 0
return (getOmega(i,j+1)-getOmega(i,j-1))/(2*dx)
def lapomega(i,j):
#if isSpot(i,j): return 0
if i==1 or i==X-2 or j==Y-2: return 0
return (getOmega(i+1,j)+getOmega(i-1,j)+getOmega(i,j+1)+getOmega(i,j-1)-4*getOmega(i,j))/(dx*dx)
def lap(arr, i, j):
return (arr[i+1][j]+arr[i-1][j]+arr[i][j+1]+arr[i][j-1]-4*arr[i][j])/(dx*dx)
def update():
global X,Y,Nu,omega,psi
omega_n=copy.deepcopy(omega)
#get omega from psi
for i in range(1,X-1):
for j in range(1,Y-1):
#if isSpot(i,j):omega[i][j]=0
omega_n[i][j]+=(\
-u(i,j)*dxw(i,j)\
-v(i,j)*dyw(i,j)\
+Nu*lapomega(i,j)\
)*dt
omega=omega_n
psi_n=copy.deepcopy(psi)
#get psi from omega
for i in range(1,X-1):
for j in range(1,Y-1):
psi_n[i][j] = 0.25 * (\
(getPsi(i+1,j)+getPsi(i-1,j)+getPsi(i,j+1)+getPsi(i,j-1)) ) + 0.25*getOmega(i,j)*dx*dx
psi=psi_n
def mean2(anarray):
return sum([sum([x*x for x in line]) for line in anarray])/float(X*Y)
def drawEnergy():
e = createMesh(X,Y)
S = X*Y
es = []
for i in range(1,X-1):
for j in range(1,Y-1):
#print i,j,u(i,j),v(i,j)
es.append( (omega[i][j],i,j) )
es=sorted(es)
for k in range(len(es)):
(o,i,j) = es[k]
if k>0 and o==es[k-1][0]:
es[k]=(o,i,j,es[k-1][3])
else:
es[k]=(o,i,j,k)
for k in range(len(es)):
(o,i,j,c) = es[k]
if isSpot(i,j):continue
color = c*255/S
#print "kijc",k,i,j,color
e[i][j]=color
return e
#graphical representation
import PIL
import Image
img = Image.new('RGB', (X,Y))
for l in range(10000):
nonzero=0
nonzero_o=0
color=0
e=drawEnergy()
for i in range(X):
for j in range(Y):
c = e[i][j]
if c!=0:color+=1
img.putpixel( (i,j), (255,c,255-c) )
if psi[i][j]!=0.0:
nonzero+=1
if omega[i][j]!=0:
nonzero_o+=1
print "nonzero_p =",nonzero,"nonzero_o =",nonzero_o,"color =",color
print "mean2_p =",mean2(psi),"mean2_o =",mean2(omega)
#save img
img.save(datadir+"hoge"+str(l)+".png", "PNG")
#update
for t in range(100):
update()
@dekosuke
Copy link
Author

dekosuke commented Nov 6, 2011

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment