Games101 大作业

发表于 | 分类于 ，

GAMES101 Final Project

选题：基于二维网格的流体动力学仿真（难度等级：2.0）

内容：编写一个二维的流体动力学模拟程序，以估算Navier-Stokes方程下不可压缩的均匀流体的行为。

NS方程是描述流体运动的方程。在研究该方程的历史中，首先被研究的就是不可压缩的流体运动。再后来提出了适用于可压缩流体的NS方程。现如今，三维空间中的NS方程组光滑解的存在性问题被设定为七个千禧年难题之一。

本次大作业要求是模拟不可压缩的均匀流体行为，所以本文就只探讨二维不可压缩NS方程的形式推导，还有应用。

NS方程内容

密度为$\rho$的不可压缩牛顿流体，受静水压力$p$和加速度$g$的作用，其运动可以描述为满足Navier-Stokes方程：

其中，$u$是速度矢量，可以表示为$u=u\vec{i}+v\vec{j}+w\vec{k}$，$t$是时间，$V$是体积。

NS方程推导

首先来看一下简单的方程$(2)$：（就当复习大物和微积分了）

$\nabla$是梯度算子，但是求梯度是用$\nabla{u}$，而求散度是用$\nabla\cdot{u}$，其实非常简单的区别，梯度是一个矢量，而散度是一个标量。散度等于零则表示这个点无源。

如何去理解呢？

我们假设有一个非常小的流体微元，由于不可压缩，所以通过一个闭合面的体积通量等于零，所以该点处的散度为零。那么该方程就显而易见了。

再来看一下复杂的方程$(1)$：

对于流体有欧拉方法描述和拉格朗日方法描述，我们这里采用欧拉方法。欧拉法是对欧氏空间中的每个点的速度和受力等情况的描述，但是该点对应的流体粒子可能会变更。

从某种角度来说，NS方程就是牛顿第二定律的一个特殊形式：

牛顿第二定律如下：

我们假设$u=u(x,t)$，代表了在时间$t$，位于位置$x$的粒子的速度。那么公式$(4)$可以写成：

把它带入$(4)$，可以得到：

其中，

所以就可以得到：

对于每一个非常小的微元，它们的质量是非常小的，在$(7)$右边除以质量$m$，就得到了

没错，就是NS方程中$(1)$的左边部分。

那么对于右半部分，很自然地想到，应该是表示该微元所受到的力的总和。我们还是先按照质量为$m$的液体考虑：

  1. 重力$mg$

  2. 液体受到周围流体的压力：

    对于压强求梯度，再乘以体积V，是$-V\nabla p$.

  3. 液体粘滞阻力：

    首先，如何定义流体的黏度？

    我们认为黏度是指流体对流动所表现的阻力。当流体流动时， 一部分在另一部分上面流动时， 就受到阻力， 这是流体的内摩擦力。

    在相邻两个液层之间，黏度导致的粘滞阻力会使液层具有速度差。

    那么在一个微元处的速度为$\nabla u$，我们再取该梯度的散度$\nabla \cdot(\nabla u)$，表示的是该速度场在此处的通量，如果通量为正，则说明在此处是有加速度，由于黏力导致形成速度差。反之亦然。

    所以，为了表示粘滞阻力我们定义粘度系数$\mu$，那么$\mu \nabla \cdot(\nabla u)$就是粘滞阻力。

综上，在一个微元处的合力就是：

我们再定义运动粘度$\nu=\mu/\rho$，除以$m$，就得到了：

这一项就是$(1)$的右半部分。

综上所述，我们由牛顿第二定律推出了NS方程的公式：

具体实现

具体实现

主函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def main():
    gui = ti.GUI('Stable-Fluid', (res, res))
    gen = MouseDataGen()
    paused = False
    while True:
        while gui.get_event(ti.GUI.PRESS):
            event = gui.event
            if event.key == ti.GUI.ESCAPE:
                exit(0)
            elif event.key == 'r':
                paused = False
                reset()
            elif event.key == 'p':
                paused = not paused

        if not paused:
            mouse_data = gen(gui)
            step(mouse_data)

        gui.set_image(dyes_pair.cur)
        gui.show()

主函数主要分为以下几个步骤：

  1. gui绘制
  2. 键鼠交互
  3. 更新数据

其中步骤1、2按照taichi的官方代码实现，步骤3根据网上资料和Games201课件得出。

接下来着重分析更新数据，也就是step()函数的实现。

step()函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def step(mouse_data):
    advect(velocities_pair.cur, velocities_pair.cur, velocities_pair.nxt)
    advect(velocities_pair.cur, dyes_pair.cur, dyes_pair.nxt)
    velocities_pair.swap()
    dyes_pair.swap()

    apply_impulse(velocities_pair.cur, dyes_pair.cur, mouse_data)

    divergence(velocities_pair.cur)
    for _ in range(p_jacobi_iters):
        pressure_jacobi(pressures_pair.cur, pressures_pair.nxt)
        pressures_pair.swap()

    subtract_gradient(velocities_pair.cur, pressures_pair.cur)

step()函数又由5个部分组成：

  1. advect()函数：根据上一帧的速度场（velocities_pair.cur）来得到当前帧的初始速度场（velocities_pair.nxt）和颜色场（dyes_pair.nxt）。
  2. apply_impluse()函数：根据鼠标的输入来调整速度场和颜色场。
  3. divergence()函数：根据速度场计算每个位置速度的散度。
  4. pressure_jacobi()函数：根据速度的散度来更新各个位置的压力，多次迭代找到平衡状态。
  5. subtract_gradient()函数：用加速度来更新速度。

其中swap()函数用来交换新旧缓冲。

advect()

1
2
3
4
5
6
@ti.kernel
def advect(vf: ti.template(), qf: ti.template(), new_qf: ti.template()):
    for i, j in vf:
        p = ti.Vector([i, j]) + 0.5
        p = backtrace(vf, p, dt)
        new_qf[i, j] = bilerp(qf, p) * dye_decay

使用半拉格朗日法对速度场进行更新，使用RK3的方法。具体见代码内backtrace()函数。在此不再赘述。

apply_impluse()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
@ti.kernel
def apply_impulse(vf: ti.template(), dyef: ti.template(),imp_data: ti.ext_arr()):
    g_dir = -ti.Vector([0, 1]) * 300
    for i, j in vf:
        omx, omy = imp_data[2], imp_data[3]
        mdir = ti.Vector([imp_data[0], imp_data[1]])
        dx, dy = (i + 0.5 - omx), (j + 0.5 - omy)
        d2 = dx * dx + dy * dy
        factor = ti.exp(-d2 / force_radius)

        dc = dyef[i, j]
        a = dc.norm()

        momentum = (mdir * f_strength * factor + g_dir * a / (1 + a)) * dt

        v = vf[i, j]
        vf[i, j] = v + momentum
        if mdir.norm() > 0.5:
            dc += ti.exp(-d2 * (4 / (res / 15)**2)) * ti.Vector([imp_data[4], imp_data[5], imp_data[6]])

        dyef[i, j] = dc

此函数主要规定了在鼠标互动时流体产生的方向，主要参考Taichi官方示例实现。

divergence()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
@ti.kernel
def divergence(vf: ti.template()):
    for i, j in vf:
        vl = sample(vf, i - 1, j)
        vr = sample(vf, i + 1, j)
        vb = sample(vf, i, j - 1)
        vt = sample(vf, i, j + 1)
        vc = sample(vf, i, j)
        if i == 0:
            vl.x = -vc.x
        if i == res - 1:
            vr.x = -vc.x
        if j == 0:
            vb.y = -vc.y
        if j == res - 1:
            vt.y = -vc.y
        velocity_divs[i, j] = (vr.x - vl.x + vt.y - vb.y) * 0.5

此函数求当前速度场的散度，其实对应了如下公式：

如果在边界，就把速度变反。

pressure_jacobi()

1
2
3
4
5
6
7
8
9
@ti.kernel
def pressure_jacobi(pf: ti.template(), new_pf: ti.template()):
    for i, j in pf:
        pl = sample(pf, i - 1, j)
        pr = sample(pf, i + 1, j)
        pb = sample(pf, i, j - 1)
        pt = sample(pf, i, j + 1)
        div = velocity_divs[i, j]
        new_pf[i, j] = (pl + pr + pb + pt - div) * 0.25

此函数求当前速度场的压强，对应了如下公式：

subtract_gradient()

1
2
3
4
5
6
7
8
@ti.kernel
def subtract_gradient(vf: ti.template(), pf: ti.template()):
    for i, j in vf:
        pl = sample(pf, i - 1, j)
        pr = sample(pf, i + 1, j)
        pb = sample(pf, i, j - 1)
        pt = sample(pf, i, j + 1)
        vf[i, j] -= 0.5 * ti.Vector([pr - pl, pt - pb])

此函数用来更新速度。对应公式：

运行结果

截取了一段运行视频，也在压缩包内。百度云网盘地址如下：

链接：https://pan.baidu.com/s/1_Ud9vUUOd8Ie6UlDT97TBA
提取码：1111

后记

这个大作业前前后后摸了好久，找了很多参考资料，发现确实还挺有难度的。特别是在代码实现这一块，以前用python写数值方法的时候都很麻烦，但是为了gui还是用python比较简单。Unity的话也有考虑过，但后来觉得放工程文件实在是太庞大了。其实在Games201上前几节课就讲了如何使用Taichi实现流体模拟，所以就想到用这个来实现大作业，也算是一个衔接吧。