3dでの拡散
- こちらのつづき
- 3dでの拡散のシミュレーション
- ナナメ方向も移動する
- 近傍のうち、中心間キョリが の18マスについて移動がある
- ナナメ方向も移動する
- 線維方向は考えていない
library(rgl) Nx<-35 Ny<-35 Nz<-35 Nt<-30 U<-tempU<-array(0,c(Nx,Ny,Nz)) #トーラスを作る A<-array(0,c(Nx,Ny,Nz)) r1<-10 r2<-5 for(x in 1:Nx){ for(y in 1:Nx){ for(z in 1:Nz){ l<-sqrt((sqrt((x-Nx/2)^2+(y-Ny/2)^2)-r1)^2+(z-Nz/2)^2) if(((r2-2)^2<=l^2)&&(l^2 <= r2^2)){A[x,y,z]<-1} }}} #ここの点の温度を固定して熱が伝わる U[20:30,20:30,0:30]<-1 U<-U*A #移動量が温度差に比例するときの係数 KX<-KY<-KZ<-array(0.1,c(Nx,Ny,Nz)) KXY<-KYZ<-KZX<-array(0.1,c(Nx,Ny,Nz)) KxY<-KyZ<-KzX<-array(0.1,c(Nx,Ny,Nz)) #マスの間の移動なので係数を移動する両側の平均値にする kX<-(KX[-Nx,,]+KX[-1,,])/2*A[-Nx,,]*A[-1,,] kY<-(KY[,-Ny,]+KY[,-1,])/2*A[,-Ny,]*A[,-1,] kZ<-(KZ[,,-Nz]+KZ[,,-1])/2*A[,,-Nz]*A[,,-1] kXY<-(KXY[-1,-1,]+KXY[-Nx,-Ny,])/2*A[-Nx,-Ny,]*A[-1,-1,] kYZ<-(KYZ[,-1,-1]+KYZ[,-Ny,-Nz])/2*A[,-Ny,-Nz]*A[,-1,-1] kZX<-(KZX[-1,,-1]+KZX[-Nx,,-Nz])/2*A[-Nx,,-Nz]*A[-1,,-1] kxY<-(KxY[-Nx,-1,]+KxY[-1,-Ny,])/2*A[-Nx,-1,]*A[-1,-Ny,] kyZ<-(KyZ[,-Ny,-1]+KyZ[,-1,-Nz])/2*A[,-Ny,-1]*A[,-1,-Nz] kzX<-(KzX[-1,,-Nz]+KzX[-Nx,,-1])/2*A[-1,,-Nz]*A[-Nx,,-1] #ナナメ方向を1/2に補正 kXY<-kXY/2 kYZ<-kYZ/2 kZX<-kZX/2 kxY<-kxY/2 kyZ<-kyZ/2 kzX<-kzX/2 for (t in 1:Nt){ #各方向の移動流の計算 kdUX<-(U[-1,,]-U[-Nx,,])*kX kdUY<-(U[,-1,]-U[,-Ny,])*kY kdUZ<-(U[,,-1]-U[,,-Nz])*kZ kdUXY<-(U[-1,-1,]-U[-Nx,-Ny,])*kXY kdUYZ<-(U[,-1,-1]-U[,-Ny,-Nz])*kYZ kdUZX<-(U[-1,,-1]-U[-Nx,,-Nz])*kZX kdUxY<-(U[-Nx,-1,]-U[-1,-Ny,])*kxY kdUyZ<-(U[,-Ny,-1]-U[,-1,-Nz])*kyZ kdUzX<-(U[-1,,-Nz]-U[-Nx,,-1])*kzX #移動した分の足し算と引き算 #18マス分の計算 U[-Nx,,]<-U[-Nx,,]+kdUX U[-1,,]<-U[-1,,]-kdUX U[,-Ny,]<-U[,-Ny,]+kdUY U[,-1,]<-U[,-1,]-kdUY U[,,-Nz]<-U[,,-Nz]+kdUZ U[,,-1]<-U[,,-1]-kdUZ U[-Nx,-Ny,]<-U[-Nx,-Ny,]+kdUXY U[-1,-1,]<-U[-1,-1,]-kdUXY U[,-Ny,-Nz]<-U[,-Ny,-Nz]+kdUYZ U[,-1,-1]<-U[,-1,-1]-kdUYZ U[-Nx,,-Nz]<-U[-Nx,,-Nz]+kdUZX U[-1,,-1]<-U[-1,,-1]-kdUZX U[-Nx,-1,]<-U[-Nx,-1,]-kdUxY U[-1,-Ny,]<-U[-1,-Ny,]+kdUxY U[,-Ny,-1]<-U[,-Ny,-1]-kdUyZ U[,-1,-Nz]<-U[,-1,-Nz]+kdUyZ U[-1,,-Nz]<-U[-1,,-Nz]-kdUzX U[-Nx,,-1]<-U[-Nx,,-1]+kdUzX #温度を固定している点に1を入れ直す U[20:30,20:30,0:30]<-1 U<-U*A plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),col=rainbow(1000)[ifelse(U < 0.001,1,U*1000)],alpha=ifelse(U < 0.001,0,1)) # plot3d(slice.index(U,1),slice.index(U,2),slice.index(U,3),col=rainbow(1000)[ifelse(U < 0.001,1,-log(U)*150)],alpha=ifelse(U < 0.001,0,1)) rgl.viewpoint(theta=30,phi=-10) # print(sum(U));print(max(U)) }