移動項 境界 ベクトル場

KABIRA2011-01-11

  • 拡散方程式のシミュレーション
  • 境界からの流出入をなくす
    • ソースのなかで (i != 1)や(i != Nx)の論理値 をかけることで境界では流出入が0になる
    • 総量の保存について
      • 初期条件によることであるが t=0 で36、 t=300 で 34.2 になった
      • 端で流出している
  • 今回は真ん中まで集まってくる
  • やはりまだスケーリングを気にしないといけない
  • 移動量の計算を3回繰り返しているの部分があるので、計算を工夫したほうがいいかもしれない
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(Nx,Ny,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<-array(c(runif(Nx*Ny)/10,rep(0,Nx*Ny*(Nt-1))),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.4
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*(max(0,V[max(i-1,1),j,1]*(i != 1))*U[max(i-1,1),j,t]-min(0,V[min(Nx,i+1),j,1]*(i != Nx))*U[min(i+1,Nx),j,t]-abs(V[i,j,1])*U[i,j,t])*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*(max(0,V[i,max(j-1,1),2]*(j != 1))*U[i,max(j-1,1),t]-min(0,V[i,min(j+1,Ny),2]*(j != Ny))*U[i,min(j+1,Ny),t]-abs(V[i,j,2])*U[i,j,t])*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,0.3)))
#image(U[,,t],col=topo.colors(50))
print(sum(U[,,t]))
}