流体描いてみたけどカルマン渦にならなかったの巻

vortex

#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()

gif
https://github.com/dekosuke/code-fragments/blob/master/test.gif