Julia MeshCat SPH 流体シミュレーション (PointCloud)
Julia MeshCat ライブラリを用いて、水柱崩壊のシミュレーションを行ってみました。
流体のシミュレーションには、SPH (Smoothed Particle Hydrodynamics) 法を
使用しています。
- MeshCat ( Julia bindings to the MeshCat WebGL viewer)
https://github.com/rdeits/MeshCat.jl
MeshCat は Three.js を用いた3D Wiewer です。簡単にmeshの描画ができるので、
MeshCat を利用しました。
- SPH
SPHによる流体シミュレーションは、様々なサイトで紹介されています。
ここで使用した数式と数値は、以下の2つのサイトを参考にしています。
1. SPH法による流体シミュレーション
https://soysoftware.sakura.ne.jp/archives/1559
2. 粒子法のプログラム
http://kamonama.blogspot.com/2009/02/blog-post_23.html
- 粒子の描画
粒子はPointCloudで描画しています。
PointCluodの描画
color = RGBA{Float32}(0.0, 1.0, 0.0, 1.0) material = PointsMaterial(color = color, size = 0.05, vertexColors = 0) setobject!(vis, PointCloud(positions), material)
PointCloudの色の設定、ポイントサイズの設定は、PointMaterialで行います。
- 画像のキャプチャー
MeshCat ウィンドウのControls から、save_image を使ってキャプチャーします。
Base.prompt()文(キー入力待ち)を挿入し、処理を停止させて、 1コマ1コマ
キャプチャーしています。動画は200コマを使用。
プログラム
# MeshCat-sph.jl using MeshCat vis = Visualizer() #open(vis) using Blink AtomShell.isinstalled() || AtomShell.install() open(vis, Window()) using ColorTypes using CoordinateTransformations using LinearAlgebra using GeometryTypes mutable struct Particle position::Point3 velocity::Point3 density pressure force::Point3 end # weight function w_rho(r, h) = r < h ? (h^2 - r^2)^3 : 0.0 w_pressure(r,h) = r < h ? (h - r)^2 : 0.0 w_viscosity(r,h) = r < h ? (h- r) : 0.0 # Gravity accel_g = [0.0, 0.0, -9.8] # constants mass = 0.00020543 pressure_stiffness = 1.0 rest_density = 600.0 h = 0.01 l0 = (mass / rest_density)^(1.0 / 3.0) # l0 = 0.0069958 simscale = 0.004 d = l0 / simscale * 0.86 # d = 1.5216 visc = 0.2 limit = 200.0 radius = 0.004 epsiron = 0.00001 extstiff = 10000.0 extdamp = 256.0 w_rho_coef = 315.0 / (64.0 * pi * h^9) w_pressure_coef = 45.0 / (pi * h^6) w_visc_coef = w_pressure_coef init_min = [ 0.0, 0.0, 0.0] pos_min = [0.0, 0.0, 0.0] pos_max = [10.0, 25.0, 20.0] nx = 5; ny = 5; nz = 10 particle_list = [] for iz = 0:nz+1, iy = 0:ny+1, ix = 0:nx+1 x = init_min[1] + d*ix; y = init_min[2] + d*iy; z = init_min[3] + d*iz position = Point3(x, y, z) velocity = Point3(0, 0, 0) density = 0.0 pressure = 0.0 force = Point3(0.0, 0.0, 0.0) push!(particle_list, Particle(position, velocity, density, pressure, force)) end particles = length(particle_list) println("number of particles: ", particles) function Simulation(DT) # neghiber_list neighbor_list = [] for i = 1:particles position_i = particle_list[i].position neighbors = [] for j = 1:particles j == i && continue position_j = particle_list[j].position distance = norm(position_j - position_i) * simscale if distance < h push!(neighbors, j) end end push!(neighbor_list, neighbors) #println("neighbors: ", neighbors) end # density(rho) for i = 1:particles position_i = particle_list[i].position sum = 0.0 map(neighbor_list[i]) do j position_j = particle_list[j].position distance = norm(position_j - position_i) * simscale sum += w_rho(distance, h) end # density particle_list[i].density = sum * mass * w_rho_coef #println("density: ", i, " ", particle_list[i].density) # pressure particle_list[i].pressure = pressure_stiffness * (particle_list[i].density - rest_density) #println("pressure: ", particle_list[i].pressure) end # force: pressure gradient, viscosity and adj(force at wall) for i = 1:particles pressure_gradient = [0.0, 0.0, 0.0] visc_term = [0.0, 0.0, 0.0] position_i = particle_list[i].position map(neighbor_list[i]) do j position_j = particle_list[j].position distance = norm(position_j - position_i) * simscale # diffji = (position_j - position_i) * simscale pweight = w_pressure_coef * w_pressure(distance, h) * diffji / distance pji = particle_list[j].pressure + particle_list[i].pressure rho_i = particle_list[i].density rho_j = particle_list[j].density pgradient = -0.5 * pji / (rho_i * rho_j) * pweight pressure_gradient .+= pgradient vweight = w_visc_coef * w_viscosity(distance, h) vji = particle_list[j].velocity - particle_list[i].velocity vterm = visc * vji / (rho_i * rho_j) * vweight visc_term .+= vterm end particle_list[i].force = pressure_gradient + visc_term #println("pressure gradient: ", particle_list[i].pressure_gradient) end # adj force accel_adj = Point3(0.0, 0.0, 0.0) diff_min = Point3(0.0, 0.0, 0.0) diff_max = Point3(0.0, 0.0, 0.0) map(particle_list) do particle accel = particle.force * mass speed = norm( accel )^2; if speed > limit^2 accel *= limit / sqrt(speed) end diff_min = 2.0 * radius .- ( particle.position - pos_min ) * simscale diff_max = 2.0 * radius .- ( pos_max - particle.position ) * simscale # Z if diff_min[3] > epsiron normal = Point3( 0.0, 0.0, 1.0 ) adj = extstiff * diff_min[3] - extdamp * dot( normal, particle.velocity ) accel += adj * normal end if diff_max[3] > epsiron normal = Point3( 0.0, 0.0, -1.0 ) adj = extstiff * diff_max[3] - extdamp * dot( normal, particle.velocity ) accel += adj * normal end # X if diff_min[1] > epsiron normal = Point3( 1.0, 0.0, 0.0 ) adj = extstiff * diff_min[1] - extdamp * dot( normal, particle.velocity ) accel += adj * normal end if diff_max[1] > epsiron normal = Point3( -1.0, 0.0, 0.0 ) adj = extstiff * diff_max[1] - extdamp * dot( normal, particle.velocity ) accel += adj * normal end # Y if diff_min[2] > epsiron normal = Point3( 0.0, 1.0, 0.0 ) adj = extstiff * diff_min[2] - extdamp * dot( normal, particle.velocity ) accel += adj * normal end if diff_max[2] > epsiron normal = Point3( 0.0, -1.0, 0.0 ) adj = extstiff * diff_max[2] - extdamp * dot( normal, particle.velocity ) accel += adj * normal end accel += accel_g # veloc = particle.velocity veloc = veloc .+ accel * DT # posit = particle.position posit = posit .+ veloc * DT / simscale particle.position = posit particle.velocity = veloc #println("position, velocity: ", particle.position, particle.velocity) end positions = map(particle_list) do particle particle.position end positions end # PointsCloud settings color = RGBA{Float32}(0.0, 1.0, 0.0, 1.0) material = PointsMaterial(color = color, size = 0.05, vertexColors = 0) N = 200 DT = 0.004 for i = 1:N println("loop: ", i) positions = Simulation(DT) positions .*= 0.1 setobject!(vis, PointCloud(positions), material) trans = Translation(0., -1, -0.5) rot = LinearMap(RotYZ(0.2, -0.2)) transform = compose(rot, trans) settransform!(vis, transform) sleep(0.01) #Base.prompt("continue?") # for capture end