移動項
- 拡散方程式のシミュレーション
- 昨日の記事に移動項を付け加える
- 拡散と移動項による変化
- 移動項について
- ベクトル場にそって移動するようにする
- ベクトル場の情報はarray Vに入れてある
- V[x,y,1]が(x,y)でのベクトルのx成分
- V[x,y,2]が(x,y)でのベクトルのy成分
- 昨日の記事のようにx方向とy方向で分けて考える
- 今回は中心へ渦を描くようなベクトル場
- 点を中心とする同心円の接線方向に対して中心へ向かっての角度をなすような向き
- 拡散しつつ、反時計回りに回る
- 斜めの移動は入っていないが、全体としては斜め方向への移動をしているようにみえる
- 境界条件が入っていないので、フィールドの外へ拡散してしまっている
- まだ修正が必要
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])) }