关于相对速度移动边界条件

类似刚体球体在液体中流动,胡老师说使用相对速度边界条件,想仔细了解一下其原理。请问大家有demo或者相关学习资料吗?

Coulomb friction

需要用库伦摩擦模型。

我写了个刚体球和弹塑性体碰撞的程序,也许可以参考下,用鼠标控制球的运动:

(主要看 friction_coeff = 0.4 这行上下文的部分。)

import taichi as ti
import numpy as np
import random

ti.init(arch=ti.gpu)  # Try to run on GPU

quality = 1  # Use a larger value for higher-res simulations
n_grid = (40 * quality, 80 * quality)
dx, inv_dx = 1 / n_grid[1], n_grid[1]
dt = 1e-4 / quality
p_vol, p_rho = (dx * 0.5)**2, 1
p_mass = p_vol * p_rho
E, nu = 0.15e4, 0.2  # Young's modulus and Poisson's ratio
mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / (
        (1 + nu) * (1 - 2 * nu))  # Lame parameters

max_num_particles = 1024 * 16
dim = 2
x = ti.Vector.field(dim, dtype=float, shape=max_num_particles)  # position
v = ti.Vector.field(dim, dtype=float, shape=max_num_particles)  # velocity
C = ti.Matrix.field(dim, dim, dtype=float, shape=max_num_particles)  # affine velocity field
F = ti.Matrix.field(dim, dim, dtype=float, shape=max_num_particles)  # deformation gradient
material = ti.field(dtype=int, shape=max_num_particles)  # material id
Jp = ti.field(dtype=float, shape=max_num_particles)  # plastic deformation

# ti.root.dense(ti.i, max_num_particles).place(x, v, C, F, material, Jp)
cur_num_particles = ti.field(ti.i32, shape=())

collider_radius = 0.1

grid_v = ti.Vector.field(dim, dtype=float,
                         shape=n_grid)  # grid node momentum/velocity
grid_m = ti.field(dtype=float, shape=n_grid)  # grid node mass


@ti.kernel
def substep(collider_x: float, collider_y: float, collider_vx: float, collider_vy: float):
    for i, j in grid_m:
        grid_v[i, j] = [0, 0]
        grid_m[i, j] = 0

    # p2g
    cur_num = cur_num_particles[None]
    for p in range(cur_num):
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[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]
        F[p] = (ti.Matrix.identity(float, 2) +
                dt * C[p]) @ F[p]  # deformation gradient update
        h = ti.exp(
            10 *
            (1.0 -
             Jp[p]))  # Hardening coefficient: snow gets harder when compressed
        h = min(h, 15.0)
        if material[p] >= 2:  # jelly, make it softer
            h = 0.3
        mu, la = mu_0 * h, lambda_0 * h
        if material[p] == 0:  # liquid
            mu = 0.0
        U, sig, V = ti.svd(F[p])
        J = 1.0
        for d in ti.static(range(2)):
            new_sig = sig[d, d]
            if material[p] == 1:  # Snow
                new_sig = min(max(sig[d, d], 1 - 2.5e-2),
                              1 + 4.5e-3)  # Plasticity
            Jp[p] *= sig[d, d] / new_sig
            sig[d, d] = new_sig
            J *= new_sig
        if material[
            p] == 0:  # Reset deformation gradient to avoid numerical instability
            F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J)
        elif material[p] == 1:
            F[p] = U @ sig @ V.transpose(
            )  # Reconstruct elastic deformation gradient after plasticity
        stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose(
        ) + ti.Matrix.identity(float, 2) * la * J * (J - 1)
        stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress
        affine = stress + p_mass * C[p]
        for i, j in ti.static(ti.ndrange(
                3, 3)):  # Loop over 3x3 grid node neighborhood
            offset = ti.Vector([i, j])
            dpos = (offset.cast(float) - fx) * dx
            weight = w[i][0] * w[j][1]
            grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
            grid_m[base + offset] += weight * p_mass

    for i, j in grid_m:
        if grid_m[i, j] > 0:  # No need for epsilon here
            grid_v[i,
                   j] = (1 / grid_m[i, j]) * grid_v[i,
                                                    j]  # Momentum to velocity

            grid_v[i, j][1] -= dt * 50  # gravity

            grid_pos = ti.Vector([i, j]) * dx
            collider_center = ti.Vector([collider_x, collider_y])
            collider_offset = grid_pos - collider_center
            if collider_offset.norm() < collider_radius:
                normal = collider_offset.normalized(1e-8)
                v_co = ti.Vector([collider_vx, collider_vy]) # collider
                v_rel = grid_v[i, j] - v_co

                v_n = v_rel.dot(normal)
                v_tan = v_rel - v_n * normal
                v_tan_len = v_tan.norm()
                friction_coeff = 0.4

                if v_n < 0:
                    v_tan_normalized = v_tan / (v_tan_len + 1e-5)
                    v_rel = max(v_tan_len + friction_coeff * v_n, 0) * v_tan_normalized
                grid_v[i, j] = v_rel + v_co

            if i < 3 and grid_v[i, j][0] < 0:
                grid_v[i, j][0] = 0  # Boundary conditions
            if i > n_grid[0] - 3 and grid_v[i, j][0] > 0: grid_v[i, j][0] = 0
            if j < 3 and grid_v[i, j][1] < 0: grid_v[i, j] = ti.Vector([0, 0])
            if j > n_grid[1] - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0

    # g2p
    for p in range(cur_num):
        base = (x[p] * inv_dx - 0.5).cast(int)
        fx = x[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, 2)
        new_C = ti.Matrix.zero(float, 2, 2)
        for i, j in ti.static(ti.ndrange(
                3, 3)):  # loop over 3x3 grid node neighborhood
            dpos = ti.Vector([i, j]).cast(float) - fx
            g_v = grid_v[base + ti.Vector([i, j])]
            weight = w[i][0] * w[j][1]
            new_v += weight * g_v
            new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
        v[p], C[p] = new_v, new_C
        x[p] += dt * v[p]  # advection


num_per_tetromino_square = 128
num_per_tetromino = num_per_tetromino_square * 4
staging_tetromino_x = ti.Vector.field(dim,
                                      dtype=float,
                                      shape=num_per_tetromino)


class StagingTetromino(object):
    def __init__(self, x):
        """
        Stores the info of the staging tetromino.

        Args
          x: a taichi field. Note this class just keeps |x| as a reference, and
             is not a taichi data_oriented class.
        """
        self._x_field = x
        self.material_idx = 0

        self._offsets = np.array([
            [[0, -1], [1, 0], [0, -2]],
            [[1, 1], [-1, 0], [1, 0]],
            [[0, -1], [-1, 0], [0, -2]],
            [[0, 1], [1, 0], [1, -1]],
            [[1, 0], [2, 0], [-1, 0]],
            [[0, 1], [1, 1], [1, 0]],
            [[-1, 0], [1, 0], [0, 1]],
        ])
        self._x_np_canonical = None

    def regenerate(self, mat, kind):
        self.material_idx = mat
        shape = (num_per_tetromino, dim)
        x = np.zeros(shape=shape, dtype=np.float32)
        for i in range(1, 4):
            # is there a more np-idiomatic way?
            begin = i * num_per_tetromino_square
            x[begin:(begin +
                     num_per_tetromino_square)] = self._offsets[kind, (i - 1)]
        x += np.random.rand(*shape)
        x *= 0.05
        self._x_np_canonical = x

    def update_center(self, center):
        self._x_field.from_numpy(np.clip(self._x_np_canonical + center, 0, 1))

    def rotate(self):
        theta = np.radians(90)
        c, s = np.cos(theta), np.sin(theta)
        m = np.array([[c, -s], [s, c]], dtype=np.float32)
        x = m @ self._x_np_canonical.T
        self._x_np_canonical = x.T


staging_tetromino = StagingTetromino(staging_tetromino_x)


@ti.kernel
def drop_staging_tetromino(mat: int):
    base = cur_num_particles[None]
    for i in staging_tetromino_x:
        bi = base + i
        x[bi] = staging_tetromino_x[i]
        material[bi] = mat
        v[bi] = ti.Matrix([0, -2])
        F[bi] = ti.Matrix([[1, 0], [0, 1]])
        Jp[bi] = 1
    cur_num_particles[None] += num_per_tetromino


def compute_tetromino_center(mouse):
    # TODO: better clipping...
    x = max(min(mouse[0], 0.45), 0.05)
    return np.array([x, 0.8], dtype=np.float32)


def main():
    gui = ti.GUI("Taichi Snow Effects",
                 res=(384, 768),
                 background_color=0x112F41)

    collider_center = np.array([0.5, 0.5])

    def gen_mat_and_kind():
        material_id = 1
        return material_id, random.randint(0, 6)

    staging_tetromino.regenerate(*gen_mat_and_kind())

    last_action_frame = -1e10
    last_mouse = gui.get_cursor_pos()
    for f in range(100000):
        if gui.get_event(ti.GUI.PRESS):
            ev_key = gui.event.key
            if ev_key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: break
            elif ev_key == ti.GUI.SPACE:
                if cur_num_particles[
                    None] + num_per_tetromino < max_num_particles:
                    drop_staging_tetromino(staging_tetromino.material_idx)
                print('# particles =', cur_num_particles[None])
                staging_tetromino.regenerate(*gen_mat_and_kind())
                last_action_frame = f
            elif ev_key == 'r':
                staging_tetromino.rotate()
        mouse = gui.get_cursor_pos()
        mouse = (mouse[0] * 0.5, mouse[1])
        staging_tetromino.update_center(compute_tetromino_center(mouse))

        collider_center = np.array(mouse)

        if f == 0:
            vx, vy = 0, 0
        else:
            vx = (mouse[0] - last_mouse[0]) / 2e-3
            vy = (mouse[1] - last_mouse[1]) / 2e-3
            last_mouse = mouse
        for s in range(int(2e-3 // dt)):
            mx, my = collider_center
            substep(mx, my, vx, vy)
        colors = np.array([
            0xA6B5F7, 0xEEEEF0, 0xED553B, 0x3255A7, 0x6D35CB, 0xFE2E44,
            0x26A5A7, 0xEDE53B
        ],
            dtype=np.uint32)
        particle_radius = 2.3
        gui.circles(x.to_numpy() * [[2, 1]],
                    radius=particle_radius,
                    color=colors[material.to_numpy()])

        if last_action_frame + 40 < f:
            gui.circles(staging_tetromino_x.to_numpy() * [[2, 1]],
                        radius=particle_radius,
                        color=int(colors[staging_tetromino.material_idx]))

        gui.circle(pos=collider_center * [2, 1], radius=collider_radius * 768)

        gui.text('Taichi Snow Effects', (0.07, 0.97), font_size=20)

        padding = 0.025
        segments = 20
        step = (1 - padding * 4) / (segments - 0.5) / 2
        gui.line(begin=(padding * 2, padding),
                 end=(1 - padding * 2, padding),
                 radius=2)
        gui.line(begin=(padding * 2, 1 - padding),
                 end=(1 - padding * 2, 1 - padding),
                 radius=2)
        gui.line(begin=(padding * 2, padding),
                 end=(padding * 2, 1 - padding),
                 radius=2)
        gui.line(begin=(1 - padding * 2, padding),
                 end=(1 - padding * 2, 1 - padding),
                 radius=2)

        gui.show()

if __name__ == '__main__':
    main()
2 个赞

刚体和流体的交互可以根据“流体是否会影响刚体的运动”分成两种:

  1. 流体不影响刚体。这种情况比较好处理,这时可以使用 @Aspen @yuanming 提到的摩擦力模型来控制碰撞行为。如果你使用的是粒子法,你可以认为是每个流体粒子与碰撞体发生刚体碰撞。网格法的处理方法可以看 Fluid Simulation for Computer Graphics 这本书里 boundary condition 相关的内容。
  2. 流体会影响刚体。则需要计算所有相关粒子对刚体产生的外力(冲量)总和。
1 个赞