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))
}