Hello, friends.
Recently I am interested with sand simulation, and I read several papers and a few examples.“Drucker-Prager Elastoplasticity for Sand Animation” and “Multi-species simulation of porous sand and water mixtures”.Both of them give me more ideas. And I noticted that, the Drucker-Prager model is used for the sand simulation. And it is refered by Mast thesis. But there is some simplify points, so is it still as same as the constitutive model in the soil mechanics. I am confused with the Project function.
And I’m also working on the algorithm development of MPM, I tried to add the sand model and set a hourglass model. The results show that it looks like the fluid (shown in the following figures). I just set the friction angle as 35° and use harding parameters and cohesion as 0.
[I used the GIMP shape function]
BTW, I noticted that there are two algorithms named as sand in the taichi, and the Project funtion are different.
Anyone have any idea about the constitutive models, how should I compare it with the soil mechanics. and why the sand simulation results looks like fluid~
Thanks a lot.
1 个赞
Hi,
Would you mind also posting complete references or links to the papers and the code you mentioned? That helps us better understand your question.
I read the paper:
(1) Drucker-prager elastoplasticity for sand animation
https://www.researchgate.net/publication/305218069_Drucker-prager_elastoplasticity_for_sand_animation
(2) Multi-species simulation of porous sand and water mixtures
https://www.researchgate.net/publication/318612741_Multi-species_simulation_of_porous_sand_and_water_mixtures
1 个赞
And the sand-water coupling test:
import taichi as ti
import numpy as np
import time
ti.init(arch=ti.cpu)
quality = 1
n_particles = 20000 * quality ** 3
n_s_particles = ti.field(dtype = int, shape = ())
n_w_particles = ti.field(dtype = int, shape = ())
n_grid = 128 * quality
dx, inv_dx = 1 / n_grid, float(n_grid)
dt = 2e-4 / quality
gravity = ti.Vector([0, 0, -9.8])
d = 3
# sand particle properties
x_s = ti.Vector.field(3, dtype = float, shape = n_particles) # position
v_s = ti.Vector.field(3, dtype = float, shape = n_particles) # velocity
C_s = ti.Matrix.field(3, 3, dtype = float, shape = n_particles) # affine velocity matrix
F_s = ti.Matrix.field(3, 3, dtype = float, shape = n_particles) # deformation gradient
phi_s = ti.field(dtype = float, shape = n_particles) # cohesion and saturation
c_C0 = ti.field(dtype = float, shape = n_particles) # initial cohesion (as maximum)
vc_s = ti.field(dtype = float, shape = n_particles) # tracks changes in the log of the volume gained during extension
alpha_s = ti.field(dtype = float, shape = n_particles) # yield surface size
q_s = ti.field(dtype = float, shape = n_particles) # harding state
# sand grid properties
grid_sv = ti.Vector.field(3, dtype = float) # grid node momentum/velocity
grid_sm = ti.field(dtype = float) # grid node mass
grid_sf = ti.Vector.field(3, dtype = float) # forces in the sand
# water particle properties
x_w = ti.Vector.field(3, dtype = float, shape = n_particles) # position
v_w = ti.Vector.field(3, dtype = float, shape = n_particles) # velocity
C_w = ti.Matrix.field(3, 3, dtype = float, shape = n_particles) # affine velocity matrix
J_w = ti.field(dtype = float, shape = n_particles) # ratio of volume increase
# water grid properties
grid_wv = ti.Vector.field(3, dtype = float) # grid node momentum/velocity
grid_wm = ti.field(dtype = float) # grid node mass
grid_wf = ti.Vector.field(3, dtype = float) # forces in the water
block_size = 16
block0 = ti.root.pointer(ti.ijk, n_grid // block_size)
block1 = block0.dense(ti.ijk, block_size)
block1.place(grid_sv, grid_sm, grid_sf, grid_wv, grid_wm, grid_wf)
# constant values
p_vol, s_rho, w_rho = (dx * 0.5) ** 3, 400, 400
s_mass, w_mass = p_vol * s_rho, p_vol * w_rho
w_k, w_gamma = 50, 3 # bulk modulus of water and gamma is a term that more stiffy penalizes large deviations from incompressibility
n, k_hat = 0.4, 0.2 # sand porosity and permeability
E_s, nu_s = 3.537e5, 0.3 # sand's Young's modulus and Poisson's ratio
mu_s, lambda_s = E_s / (2 * (1 + nu_s)), E_s * nu_s / ((1 + nu_s) * (1 - 2 * nu_s)) # sand's Lame parameters
mu_b = 0.75 # coefficient of friction
a, b, c0, sC = -3.0, 0, 1e-2, 0.15
# The scalar function h_s is chosen so that the multiplier function is twice continuously differentiable
@ti.func
def h_s(z):
ret = 0.0
if z < 0: ret = 1
if z > 1: ret = 0
ret = 1 - 10 * (z ** 3) + 15 * (z ** 4) - 6 * (z ** 5)
return ret
# multiplier
sqrt2 = ti.sqrt(2)
@ti.func
def h(e):
u = e.trace() / sqrt2
v = ti.abs(ti.Vector([e[0, 0] - u / sqrt2, e[1, 1] - u / sqrt2]).norm())
fe = c0 * (v ** 4) / (1 + v ** 3)
ret = 0.0
if u + fe < a + sC: ret = 1
if u + fe > b + sC: ret = 0
ret = h_s((u + fe - a - sC) / (b - a))
return ret
state = ti.field(dtype = int, shape = n_particles)
pi = 3.14159265358979
@ti.func
def project(e0, p):
e = e0 + vc_s[p] / d * ti.Matrix.identity(float, 3) # volume correction treatment
e += (c_C0[p] * (1.0 - phi_s[p])) / (d * alpha_s[p]) * ti.Matrix.identity(float, 3) # effects of cohesion
ehat = e - e.trace() / d * ti.Matrix.identity(float, 3)
Fnorm = ti.sqrt(ehat[0, 0] ** 2 + ehat[1, 1] ** 2 + ehat[2, 2] ** 2) # Frobenius norm
yp = Fnorm + (d * lambda_s + 2 * mu_s) / (2 * mu_s) * e.trace() * alpha_s[p] # delta gamma
new_e = ti.Matrix.zero(float, 3, 3)
delta_q = 0.0
if Fnorm <= 0 or e.trace() > 0: # Case II:
new_e = ti.Matrix.zero(float, 3, 3)
delta_q = ti.sqrt(e[0, 0] ** 2 + e[1, 1] ** 2 + e[2, 2] ** 2)
state[p] = 0
elif yp <= 0: # Case I:
new_e = e0 # return initial matrix without volume correction and cohesive effect
delta_q = 0
state[p] = 1
else: # Case III:
new_e = e - yp / Fnorm * ehat
delta_q = yp
state[p] = 2
return new_e, delta_q
h0, h1, h2, h3 = 35, 9, 0.2, 10
@ti.func
def hardening(dq, p): # The amount of hardening depends on the amount of correction that occurred due to plasticity
q_s[p] += dq
phi = h0 + (h1 * q_s[p] - h3) * ti.exp(-h2 * q_s[p])
phi = phi / 180 * pi # details in Table. 3: Friction angle phi_F and hardening parameters h0, h1, and h3 are listed in degrees for convenience
sin_phi = ti.sin(phi)
alpha_s[p] = ti.sqrt(2 / 3) * (2 * sin_phi) / (3 - sin_phi)
@ti.kernel
def substep():
# set zero initial state for both water/sand grid
for i, j, k in grid_sm:
grid_sv[i, j, k], grid_wv[i, j, k] = [0, 0, 0], [0, 0, 0]
grid_sm[i, j, k], grid_wm[i, j, k] = 0, 0
grid_sf[i, j, k], grid_wf[i, j, k] = [0, 0, 0], [0, 0, 0]
# P2G (sand's part)
for p in range(n_s_particles[None]):
base = (x_s[p] * inv_dx - 0.5).cast(int)
if base[0] < 0 or base[1] < 0 or base[2] < 0 or base[0] >= n_grid - 2 or base[1] >= n_grid - 2 or base[2] >= n_grid - 2:
continue
fx = x_s[p] * inv_dx - base.cast(float)
# Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2]
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2]
U, sig, V = ti.svd(F_s[p])
inv_sig = sig.inverse()
e = ti.Matrix([[ti.log(sig[0, 0]), 0, 0], [0, ti.log(sig[1, 1]), 0], [0, 0, ti.log(sig[2, 2])]])
stress = U @ (2 * mu_s * inv_sig @ e + lambda_s * e.trace() * inv_sig) @ V.transpose() # formula (25)
stress = (-p_vol * 4 * inv_dx * inv_dx) * stress @ F_s[p].transpose()
# stress *= h(e)
# print(h(e))
affine = s_mass * C_s[p]
for i, j, k in ti.static(ti.ndrange(3, 3, 3)):
offset = ti.Vector([i, j, k])
dpos = (offset.cast(float) - fx) * dx
weight = w[i][0] * w[j][1] * w[k][2]
grid_sv[base + offset] += weight * (s_mass * v_s[p] + affine @ dpos)
grid_sm[base + offset] += weight * s_mass
grid_sf[base + offset] += weight * stress @ dpos
# P2G (water's part):
for p in range(n_w_particles[None]):
base = (x_w[p] * inv_dx - 0.5).cast(int)
fx = x_w[p] * inv_dx - base.cast(float)
# Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2]
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2]
stress = w_k * (1 - 1 / (J_w[p] ** w_gamma))
stress = (-p_vol * 4 * inv_dx * inv_dx) * stress * J_w[p]
# stress = -4 * 400 * p_vol * (J_w[p] - 1) / dx ** 2 (special case when gamma equals to 1)
affine = w_mass * C_w[p]
# affine = ti.Matrix([[stress, 0], [0, stress]]) + w_mass * C_w[p]
for i, j, k in ti.static(ti.ndrange(3, 3, 3)):
offset = ti.Vector([i, j, k])
dpos = (offset.cast(float) - fx) * dx
weight = w[i][0] * w[j][1] * w[k][2]
grid_wv[base + offset] += weight * (w_mass * v_w[p] + affine @ dpos)
grid_wm[base + offset] += weight * w_mass
grid_wf[base + offset] += weight * stress * dpos
# Update Grids Momentum
for i, j, k in grid_sm:
if grid_sm[i, j, k] > 0:
grid_sv[i, j, k] = (1 / grid_sm[i, j, k]) * grid_sv[i, j, k] # Momentum to velocity
if grid_wm[i, j, k] > 0:
grid_wv[i, j, k] = (1 / grid_wm[i, j, k]) * grid_wv[i, j, k]
# Momentum exchange
cE = (n ** 3 * w_rho * gravity[2]) / k_hat # drag coefficient
if grid_sm[i, j, k] > 0 and grid_wm[i, j, k] > 0:
sm, wm = grid_sm[i, j, k], grid_wm[i, j, k]
sv, wv = grid_sv[i, j, k], grid_wv[i, j, k]
d = cE * sm * wm
M = ti.Matrix([[sm, 0], [0, wm]])
D = ti.Matrix([[-d, d], [d, -d]])
V = ti.Matrix.rows([sv, wv])
G = ti.Matrix.rows([gravity, gravity])
F = ti.Matrix.rows([grid_sf[i, j, k], grid_wf[i, j, k]])
A = M + dt * D
B = M @ V + dt * (M @ G + F)
X = A.inverse() @ B
grid_sv[i, j, k], grid_wv[i, j, k] = ti.Vector([X[0, 0], X[0, 1], X[0, 2]]), ti.Vector([X[1, 0], X[1, 1], X[1, 2]])
elif grid_sm[i, j, k] > 0:
grid_sv[i, j, k] += dt * (gravity + grid_sf[i, j, k] / grid_sm[i, j, k]) # Update explicit force
elif grid_wm[i, j, k] > 0:
grid_wv[i, j, k] += dt * (gravity + grid_wf[i, j, k] / grid_wm[i, j, k])
if grid_sm[i, j, k] > 0:
if i < 3 and grid_sv[i, j, k][0] < 0: grid_sv[i, j, k][0] = 0 # Boundary conditions
if i > n_grid - 3 and grid_sv[i, j, k][0] > 0: grid_sv[i, j, k][0] = 0
if j < 3 and grid_sv[i, j, k][1] < 0: grid_sv[i, j, k][1] = 0
if j > n_grid - 3 and grid_sv[i, j, k][1] > 0: grid_sv[i, j, k][1] = 0
if k < 3 and grid_sv[i, j, k][2] < 0: grid_sv[i, j, k][2] = 0
if k > n_grid - 3 and grid_sv[i, j, k][2] > 0: grid_sv[i, j, k][2] = 0
if grid_wm[i, j, k] > 0:
if i < 3 and grid_wv[i, j, k][0] < 0: grid_wv[i, j, k][0] = 0 # Boundary conditionw
if i > n_grid - 3 and grid_wv[i, j, k][0] > 0: grid_wv[i, j, k][0] = 0
if j < 3 and grid_wv[i, j, k][1] < 0: grid_wv[i, j, k][1] = 0
if j > n_grid - 3 and grid_wv[i, j, k][1] > 0: grid_wv[i, j, k][1] = 0
if k < 3 and grid_wv[i, j, k][2] < 0: grid_wv[i, j, k][2] = 0
if k > n_grid - 3 and grid_wv[i, j, k][2] > 0: grid_wv[i, j, k][2] = 0
# G2P (water's part)
for p in range(n_w_particles[None]):
base = (x_w[p] * inv_dx - 0.5).cast(int)
fx = x_w[p] * inv_dx - base.cast(float)
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
new_v = ti.Vector.zero(float, 3)
new_C = ti.Matrix.zero(float, 3, 3)
for i, j, k in ti.static(ti.ndrange(3, 3, 3)):
dpos = ti.Vector([i, j, k]).cast(float) - fx
g_v = grid_wv[base + ti.Vector([i, j, k])]
weight = w[i][0] * w[j][1] * w[k][2]
new_v += weight * g_v
new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
J_w[p] = (1 + dt * new_C.trace()) * J_w[p]
v_w[p], C_w[p] = new_v, new_C
x_w[p] += dt * v_w[p]
# G2P (sand's part)
for p in range(n_s_particles[None]):
base = (x_s[p] * inv_dx - 0.5).cast(int)
if base[0] < 0 or base[1] < 0 or base[2] < 0 or base[0] >= n_grid - 2 or base[1] >= n_grid - 2 or base[2] >= n_grid - 2:
continue
fx = x_s[p] * inv_dx - base.cast(float)
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2]
new_v = ti.Vector.zero(float, 3)
new_C = ti.Matrix.zero(float, 3, 3)
phi_s[p] = 0.0 # Saturation
for i, j, k in ti.static(ti.ndrange(3, 3, 3)): # loop over 3x3 grid node neighborhood
dpos = ti.Vector([i, j, k]).cast(float) - fx
g_v = grid_sv[base + ti.Vector([i, j, k])]
weight = w[i][0] * w[j][1] * w[k][2]
new_v += weight * g_v
new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
if grid_sm[base + ti.Vector([i, j, k])] > 0 and grid_wm[base + ti.Vector([i, j, k])] > 0:
phi_s[p] += weight # formula (24)
F_s[p] = (ti.Matrix.identity(float, 3) + dt * new_C) @ F_s[p]
v_s[p], C_s[p] = new_v, new_C
x_s[p] += dt * v_s[p]
U, sig, V = ti.svd(F_s[p])
e = ti.Matrix([[ti.log(sig[0, 0]), 0, 0], [0, ti.log(sig[1, 1]), 0], [0, 0, ti.log(sig[2, 2])]])
new_e, dq = project(e, p)
hardening(dq, p)
new_F = U @ ti.Matrix([[ti.exp(new_e[0, 0]), 0, 0], [0, ti.exp(new_e[1, 1]), 0], [0, 0, ti.exp(new_e[2, 2])]]) @ V.transpose()
vc_s[p] += -ti.log(new_F.determinant()) + ti.log(F_s[p].determinant()) # formula (26)
F_s[p] = new_F
@ti.kernel
def initialize():
n_s_particles[None] = 10000 * quality ** 3
for i in range(n_s_particles[None]):
x_s[i] = [ti.random() * 0.25 + 0.4, ti.random() * 0.25 + 0.4, ti.random() * 0.4 + 0.01]
v_s[i] = ti.Matrix([0, 0, 0])
F_s[i] = ti.Matrix([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
c_C0[i] = -0.01
alpha_s[i] = 0.267765
n_w_particles[None] = 0
@ti.kernel
def update_jet():
if n_w_particles[None] < 20000 - 50:
for i in range(n_w_particles[None], n_w_particles[None] + 50):
x_w[i] = [ti.random() * 0.03 + 0.92, 0.485 + ti.random() * 0.03, ti.random() * 0.03 + 0.5]
v_w[i] = ti.Matrix([-1.5, 0, 0])
J_w[i] = 1
n_w_particles[None] += 50
@ti.func
def color_lerp(r1, g1, b1, r2, g2, b2, t):
return int((r1 * (1 - t) + r2 * t) * 0x100) * 0x10000 + int((g1 * (1 - t) + g2 * t) * 0x100) * 0x100 + int((b1 * (1 - t) + b2 * t) * 0x100)
color_s = ti.field(dtype = int, shape = n_particles)
color_w = ti.field(dtype = int, shape = n_particles)
@ti.kernel
def update_color():
for i in range(n_s_particles[None]):
color_s[i] = color_lerp(0.521, 0.368, 0.259, 0.318, 0.223, 0.157, phi_s[i])
for i in range(n_w_particles[None]):
color_w[i] = color_lerp(0.2, 0.231, 0.792, 0.867, 0.886, 0.886, v_w[i].norm() / 3.0)
pos_s = ti.Vector.field(2, dtype = float, shape = n_particles)
pos_w = ti.Vector.field(2, dtype = float, shape = n_particles)
@ti.kernel
def update_pos():
for i in range(n_s_particles[None]):
pos_s[i] = ti.Vector([x_s[i][0], x_s[i][2]])
# pos_s[i] = ti.Vector([x_s[i][0], x_s[i][1]])
for i in range(n_w_particles[None]):
pos_w[i] = ti.Vector([x_w[i][0], x_w[i][2]])
# pos_w[i] = ti.Vector([x_w[i][0], x_w[i][1]])
initialize()
project_view = False
frame = 0
gui = ti.GUI("2D Dam", res = 512, background_color = 0xFFFFFF)
while True:
for e in gui.get_events(ti.GUI.PRESS):
if e.key == gui.SPACE: project_view = not project_view
elif e.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]:
exit()
update_jet()
for s in range(5):
substep()
update_pos()
if project_view:
gui.circles(pos_w.to_numpy(), radius = 1.5, color = 0x068587)
colors = np.array([0xFF0000, 0x00FF00, 0x0000FF], dtype = np.uint32)
gui.circles(pos_s.to_numpy(), radius = 1.5, color = colors[state.to_numpy()])
else:
update_color()
gui.circles(pos_w.to_numpy(), radius = 1.5, color = color_w.to_numpy())
gui.circles(pos_s.to_numpy(), radius = 1.5, color = color_s.to_numpy())
# gui.show(f'{frame:06d}.png')
gui.show()
frame += 1
目前 taichi 提供的 svd 函数实现有问题，返回的矩阵都是 det 为 +1 的，而且输入必须是可逆矩阵，所以这可能是一个原因。我们会在 1.0 的 release 以后重构这部分代码。
先试试使用你原来的实现方式？
感谢您的回复与建议。我暂时不确定问题的具体原因，但是我认为不是svd函数导致的。目前我是自己根据开源库中的svd改成的device函数。目前我怀疑的原因是关于Project函数中的操作或者参数设置等。我再检查一下~另外想问下，是否有关于下图沙漏模型的算例：
请问一下，您找到这个沙漏模型的的算例了吗？