1 粒子法とは

粒子法とは，流体を粒子の集合とみなしてシミュレーションを行う手法である．

2 シミュレーション

3 初期設定

まず，初期設定を行う

L0 <- 0.1
nu <- 0.04
ratio <- 0.3
DT <- 0.005
tend <- 2.0
col_dist <- 0.6*L0
col_dist2 <- col_dist^2
ce <- 0.6
save_interval <- 0.05
movie_interval <- 0.1
err <- FALSE
g <- 9.8
re <- 2.25*L0
re2 <- re^2

height <- ceiling(1/L0)
width <- ceiling(ratio/L0)
nparticle <- height*width
nwall <- (3*height+8)*4
nall <- nparticle + nwall

stepnum <- ceiling(tend/save_interval)
u <- array(rep(0,nall*6*stepnum), dim=c(nall,6,stepnum))
dimnames(u) <- list(NULL, c("rx","ry","vx","vy", "ax", "ay"), NULL)
type <- rep(1,nall)  # type=1 if usual, 2 if realwall, 3 if dummywall

distance2_matrix <- array(rep(0,nall*nall), dim=c(nall,nall))
w_matrix <- array(rep(0,nall*nall), dim=c(nall,nall))
forceDT_matrix <-  array(rep(0,nall*nall), dim=c(nall,nall))
a_matrix <- array(rep(0,nall*nall), dim=c(nall,nall))
b_vector <- rep(0,nall)
pressure_vector <- rep(0,nall)
lap_matrix <- array(rep(0,nall*nall), dim=c(nall,nall))

3.1 初期配置

次に，粒子を初期配置に並べる．

#####
# initial state
## particles
for(i in 1:height){
  for(j in 1:width){
    u[(i-1)*width+j,"rx",1] <- j * L0
    u[(i-1)*width+j,"ry",1] <- i * L0
  }
}
## left and right wall
for(i in 1:(height+4)){
  for(j in 1:4){
    u[nparticle+(i-1)*4+j,"rx",1] <- (1-j)*L0
    u[nparticle+(i-1)*4+j,"ry",1] <- (i-4)*L0
    type[nparticle+(i-1)*4+j] <- ifelse(j<=2,2,3)
    
    u[nparticle+(height+4)*4+(i-1)*4+j,"rx",1] <- (height+j)*L0
    u[nparticle+(height+4)*4+(i-1)*4+j,"ry",1] <- (i-4)*L0
    type[nparticle+(height+4)*4+(i-1)*4+j] <- ifelse(j<=2,2,3)
  }
}
## below wall
for(i in 1:4){
  for(j in 1:height){
    u[nparticle+(height+4)*8+(i-1)*height+j,"rx",1] <- j*L0
    u[nparticle+(height+4)*8+(i-1)*height+j,"ry",1] <- (i-4)*L0
    type[nparticle+(height+4)*8+(i-1)*height+j] <- ifelse(i>=3,2,3)
  }
}

k <- 1
plot(u[,"rx",k],u[,"ry",k],  xlim=c(-0.1,1.1), ylim=c(-0.1,1.1), col=ifelse(type==1,"blue","black"))

3.2 初期情報設定

初期情報の設定と関数定義．

#####
# initial density
i0 <- floor(height/2)*width+floor(width/2)
distance2_vector <- (u[,"rx",1] - u[i0,"rx",1])^2 + (u[,"ry",1] - u[i0,"ry",1])^2
w_vector <- ifelse(distance2_vector<re2 & distance2_vector>0,re/sqrt(distance2_vector)-1,0)
n0 <- sum(w_vector)
lambda0 <- sum(distance2_vector*w_vector)
rho0 <- 1000

distance2 <- function(x,y){(x-y)^2}

3.3 計算

t <- 0
k <- 1
while(t<=tend){
  if(t>=k*save_interval){
    u[,,k+1] <- u[,,k]
    k <- k + 1
  }
  t <- t + DT

  distance2_matrix <- outer(u[,"rx",k],u[,"rx",k],distance2)+outer(u[,"ry",k],u[,"ry",k],distance2)
  w_matrix <- ifelse(distance2_matrix<re2 & distance2_matrix>0,re/sqrt(distance2_matrix)-1,0)
  u[,"ax",k] <- nu*4/lambda0*colSums(w_matrix*outer(u[,"vx",k],u[,"vx",k],"-"))
  u[,"ay",k] <- nu*4/lambda0*colSums(w_matrix*outer(u[,"vy",k],u[,"vy",k],"-")) - g
  u[type>1,"ax",k] <- 0
  u[type>1,"ay",k] <- 0
  u[,"vx",k] <- u[,"vx",k] + u[,"ax",k]*DT
  u[,"vy",k] <- u[,"vy",k] + u[,"ay",k]*DT
  u[,"rx",k] <- u[,"rx",k] + u[,"vx",k]*DT
  u[,"ry",k] <- u[,"ry",k] + u[,"vy",k]*DT
  
  diff_rx <- outer(u[,"rx",k],u[,"rx",k],"-")
  diff_ry <- outer(u[,"ry",k],u[,"ry",k],"-")
  distance2_matrix <- diff_rx^2+diff_ry^2
  diff_vx <- outer(u[,"vx",k],u[,"vx",k],"-")
  diff_vy <- outer(u[,"vy",k],u[,"vy",k],"-")
  forceDT_matrix <- ifelse(distance2_matrix<col_dist2 & distance2_matrix>0,
                           ce*(diff_vx*diff_rx+ diff_vy*diff_ry)/distance2_matrix,
                           0)
  u[,"ax",k] <- colSums(forceDT_matrix*diff_rx)
  u[,"ay",k] <- colSums(forceDT_matrix*diff_ry)
  u[type>1,"ax",k] <- 0
  u[type>1,"ay",k] <- 0
  u[,"vx",k] <- u[,"vx",k] + u[,"ax",k]
  u[,"vy",k] <- u[,"vy",k] + u[,"ay",k]
  u[,"rx",k] <- u[,"rx",k] + u[,"ax",k]*DT
  u[,"ry",k] <- u[,"ry",k] + u[,"ay",k]*DT
  
  diff_rx <- outer(u[,"rx",k],u[,"rx",k],"-")
  diff_ry <- outer(u[,"ry",k],u[,"ry",k],"-")
  distance2_matrix <- diff_rx^2+diff_ry^2
  w_matrix <- ifelse(distance2_matrix<re2 & distance2_matrix>0,re/sqrt(distance2_matrix)-1,0)
  ns_vector <- colSums(w_matrix)
  inner_particle <- (ns_vector>0.97*n0 & type<=2)
  cof <- 1/rho0*4/lambda0
  a_matrix <- ifelse(outer(inner_particle,inner_particle,"&"),
                     ifelse(distance2_matrix>0,-cof*w_matrix,cof*ns_vector),
                     ifelse(distance2_matrix>0,0,1))
  b_vector <- ifelse(inner_particle, 1 / (DT^2) * (ns_vector/n0 -1), 0)
  pressure_vector <- try(solve(a_matrix,b_vector))
  if(class(pressure_vector) == "try-error"){err <- T; break }
  
  lap_matrix <- w_matrix*outer(pressure_vector,pressure_vector,"-")
  lap_matrix <- ifelse(distance2_matrix>0,lap_matrix/distance2_matrix,0)
  u[,"ax",k] <- -1/rho0*2/n0*colSums(lap_matrix*diff_rx)
  u[,"ay",k] <- -1/rho0*2/n0*colSums(lap_matrix*diff_ry)
  u[type>1,"ax",k] <- 0
  u[type>1,"ay",k] <- 0
  u[,"vx",k] <- u[,"vx",k] + u[,"ax",k]*DT
  u[,"vy",k] <- u[,"vy",k] + u[,"ay",k]*DT
  u[,"rx",k] <- u[,"rx",k] + u[,"ax",k]*DT^2/2
  u[,"ry",k] <- u[,"ry",k] + u[,"ay",k]*DT^2/2
}

3.4 動画出力

データを動画で出力する．

  for(k in 1:stepnum){
    plot(u[,"rx",k],u[,"ry",k],  xlim=c(-0.1,1.1), ylim=c(-0.1,1.1), col=ifelse(type==1,"blue","black"))
  }