移動項

KABIRA2011-01-10

  • 拡散方程式のシミュレーション
  • 拡散と移動項による変化
    • このときのCは移動項が入っている
    • その他の項についてはこちらで考えている
  • 移動項について
    • ベクトル場にそって移動するようにする
    • ベクトル場の情報はarray Vに入れてある
      • V[x,y,1]が(x,y)でのベクトルのx成分
      • V[x,y,2]が(x,y)でのベクトルのy成分
    • 今回は中心へ渦を描くようなベクトル場
      • O \hspace{8} (O_x,\hspace{4} O_y)を中心とする同心円の接線方向に対して中心へ向かって\phiの角度をなすような向き
  • 拡散しつつ、反時計回りに回る
  • 斜めの移動は入っていないが、全体としては斜め方向への移動をしているようにみえる
  • 境界条件が入っていないので、フィールドの外へ拡散してしまっている
  • まだ修正が必要
X<-30    #空間
Y<-30
T<-100    #時間


dx<-1
dy<-1
dt<-0.5  #時間間隔

Nx<-X/dx   #格子の数
Ny<-Y/dy
Nt<-T/dt   #計算する回数

phi<-pi/8
Ox<-X/2
Oy<-Y/2
V<-array(0,c(X,Y,2))

f<-function(x,y,Ox,Oy,phi){
#r<-sqrt((x-Ox)^2+(y-Oy)^2)
r<-1/3
th<-atan2(y-Oy,x-Ox)
v<-c(-r*sin(th+phi),r*cos(th+phi))
v
}

for(i in 1:Nx){
for(j in 1:Ny){
V[i,j,]<-f(i,j,Ox,Oy,phi)
}
}

U<-array(0,c(Nx,Ny,Nt))
U[5:10,5:10,1]<-1     #初期分布
c<-0.5
#c1<-0.1
#c2<-0.1
#nux<-(1-c1^2)*dx^2/dt
#nuy<-(1-c2^2)*dy^2/dt
 
nux<-0.1
nuy<-0.4

for(t in 1:(Nt-1)){
for(i in 1:Nx){
for(j in 1:Ny){
kx<-length(max(1,i-1):min(Nx,i+1))
ky<-length(max(1,j-1):min(Ny,j+1))

dux<-nux/2*(sum(U[max(1,i-1):min(Nx,i+1),j,t])-kx*U[i,j,t])/(dx^2)*dt-c*V[i,j,1]*(U[min(i+1,Nx),j,t]-U[max(i-1,1),j,t])/((kx-1)*dx)*dt
duy<-nuy/2*(sum(U[i,max(1,j-1):min(Ny,j+1),t])-ky*U[i,j,t])/(dy^2)*dt-c*V[i,j,2]*(U[i,min(j+1,Ny),t]-U[i,max(j-1,1),t])/((ky-1)*dy)*dt

du<-dux+duy

U[i,j,t+1]<-U[i,j,t]+du
}
}
}
min(U)

for(t in 1:Nt){
persp(U[,,t],col=6,theta=20,phi=60,zlim=(c(0,1)))
#image(U[,,t],col=topo.colors(10))
#print(sum(U[,,t]))
}