流体模拟入门指南

# mikeash.com: 流体模拟傻瓜教程 来源：https://www.mikeash.com/pyblog/fluid-simulation-for-dummies.html ## 流体模拟傻瓜教程 2005年春夏之际，我撰写了一篇关于高性能实时3D流体模拟与体渲染的硕士论文。我所采用的流体模拟基本原理虽然直截了当，但理解起来却异常艰难。现有的参考资料都非常出色，但对当时的我来说，物理和数学成分似乎偏重了些。由于找不到专门针对我这种思维模式的内容，我决定写下这篇一年前就希望看到的文章。本着这一目标，我将逐步展示如何实现简单的3D流体模拟，并尽可能将重点放在实际编程上。 本文介绍的模拟代码和思路基于Jos Stam的论文《Real-Time Fluid Dynamics for Games》（http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf）。如果你想更深入地了解原理，那是绝佳的参考。你还可以在我的硕士论文（http://www.mikeash.com/?page=thesis/）中阅读关于如何并行化模拟以及如何在3D中渲染输出的完整内容。 ### 基础知识 流体模拟建立在纳维-斯托克斯方程之上，这些方程非常重要，但如果没有良好的物理和微分方程基础，理解起来会很困难。因此，我打算在很大程度上忽略它们，仅作非常简要的解释。 将流体想象成由许多小格子或单元组成的集合。每个格子都有各种属性，比如速度和密度。这些格子与其相邻格子相互作用，影响彼此的速度和密度。现实世界中的流体可以看作是由数量庞大到极小的格子组成，它们每秒与邻居交互无数次。当然，真正的计算机无法每秒处理无数次交互，也无法处理无数个小格子，因此我们必须做出妥协。我们将流体分解成数量合理的格子，并尝试每秒进行若干次交互。 为了简化问题，我们仅研究不可压缩流体。空气是可压缩流体的一个例子：你挤压它，它体积会变小。水是不可压缩流体的一个例子：你挤压它，它会反弹，体积不会缩小。不可压缩流体模拟起来更简单，因为其密度和压力始终保持恒定。 当然，如果水里没有任何东西，流动的水就太无聊了，因为你根本看不到它在动！所以我们要添加一些染料，这样才能实际看到水的流动。水的密度处处相等，但有些部分染料浓度更高，这种差异让我们能看到东西在移动。请记住，每当我提到“密度”时，实际上指的是*染料*的密度，而不是*流体*的密度。（这一点我花了大约六个月才搞明白，所以我才如此强调。） ### 数据结构 这些格子将用多个3D数组来表示。当然，由于C语言不擅长多维数组，我们将使用一维数组并模拟出额外的两维。为此，我们有了第一段代码：一个 `#define` 宏，它接收三个坐标，并输出一维数组中对应的索引： ```c #define IX(x, y, z) ((x) + (y) * N + (z) * N * N) ``` 接下来，我们具体看看如何表示一个流体立方体。我们需要以下信息：其大小（此代码仅处理各边长度相等的完美立方体），时间步长，扩散系数（物质在流体中扩散的速度），以及粘度（流体的稠度）。我们还需要一个密度数组和三个速度数组，分别对应x、y、z方向。此外，每个数组都需要一个临时存储空间，以便在计算新值时保留旧值。综合以上所有元素，我们得到了流体立方体结构体。 ```c struct FluidCube { int size; float dt; float diff; float visc; float *s; float *density; float *Vx; float *Vy; float *Vz; float *Vx0; float *Vy0; float *Vz0; }; typedef struct FluidCube FluidCube; ``` 为了确保所有内容正确初始化，我们需要一个函数来填充整个结构体： ```c FluidCube *FluidCubeCreate(int size, int diffusion, int viscosity, float dt) { FluidCube *cube = malloc(sizeof(*cube)); int N = size; cube->size = size; cube->dt = dt; cube->diff = diffusion; cube->visc = viscosity; cube->s = calloc(N * N * N, sizeof(float)); cube->density = calloc(N * N * N, sizeof(float)); cube->Vx = calloc(N * N * N, sizeof(float)); cube->Vy = calloc(N * N * N, sizeof(float)); cube->Vz = calloc(N * N * N, sizeof(float)); cube->Vx0 = calloc(N * N * N, sizeof(float)); cube->Vy0 = calloc(N * N * N, sizeof(float)); cube->Vz0 = calloc(N * N * N, sizeof(float)); return cube; } ``` 最后，我们需要能够在完成后销毁这个结构体： ```c void FluidCubeFree(FluidCube *cube) { free(cube->s); free(cube->density); free(cube->Vx); free(cube->Vy); free(cube->Vz); free(cube->Vx0); free(cube->Vy0); free(cube->Vz0); free(cube); } ``` 由于这仅初始化了一个静止、无染料的立方体，我们需要方法来添加染料： ```c void FluidCubeAddDensity(FluidCube *cube, int x, int y, int z, float amount) { int N = cube->size; cube->density[IX(x, y, z)] += amount; } ``` 以及添加速度： ```c void FluidCubeAddVelocity(FluidCube *cube, int x, int y, int z, float amountX, float amountY, float amountZ) { int N = cube->size; int index = IX(x, y, z); cube->Vx[index] += amountX; cube->Vy[index] += amountY; cube->Vz[index] += amountZ; } ``` ### 模拟流程 我们使用三个主要操作来执行模拟步骤。 **diffuse（扩散）** – 将一滴酱油滴入水中，你会注意到它不会静止不动，而是会扩散开来。即使水和酱油都完全静止，这种现象也会发生。这被称为扩散。我们不仅在显而易见的情况（让染料扩散）中使用扩散，还在不太明显的情况（让流体的速度扩散）中使用它。 **project（投影）** – 还记得我说过我们只模拟不可压缩流体吗？这意味着每个格子中的流体量必须保持恒定。也就是说，流入的流体量必须恰好等于流出的流体量。其他操作往往会破坏这一点，导致某些格子净流出，某些净流入。此操作会遍历所有格子，修复它们，使一切达到平衡。 **advect（平流）** – 每个格子都有一组速度，这些速度使得物质移动。这被称为平流。与扩散一样，平流既适用于染料，也适用于速度本身。 此外，还有两个子程序被这些操作使用。 **set_bnd** – 这是“设置边界”的缩写，用于防止流体从立方体中泄漏。尽管它实际上不会泄漏，因为只是内存中的模拟，但没有墙壁会严重破坏模拟代码。通过将最外层单元视为墙壁来添加墙壁。基本上，与这个最外层相邻的层中的所有速度都会被镜像。因此，当次外层中有朝向墙壁的速度时，墙壁会获得一个完全抵消它的速度。 **lin_solve** – 老实说，我并不确切知道这个函数的作用和工作原理。我只知道它被同时用于 **diffuse** 和 **project**。它在求解某种线性微分方程，尽管具体如何求解以及求解什么对我来说并不完全清楚。 好了，现在我们可以看看主模拟函数了： ```c void FluidCubeStep(FluidCube *cube) { int N = cube->size; float visc = cube->visc; float diff = cube->diff; float dt = cube->dt; float *Vx = cube->Vx; float *Vy = cube->Vy; float *Vz = cube->Vz; float *Vx0 = cube->Vx0; float *Vy0 = cube->Vy0; float *Vz0 = cube->Vz0; float *s = cube->s; float *density = cube->density; diffuse(1, Vx0, Vx, visc, dt, 4, N); diffuse(2, Vy0, Vy, visc, dt, 4, N); diffuse(3, Vz0, Vz, visc, dt, 4, N); project(Vx0, Vy0, Vz0, Vx, Vy, 4, N); advect(1, Vx, Vx0, Vx0, Vy0, Vz0, dt, N); advect(2, Vy, Vy0, Vx0, Vy0, Vz0, dt, N); advect(3, Vz, Vz0, Vx0, Vy0, Vz0, dt, N); project(Vx, Vy, Vz, Vx0, Vy0, 4, N); diffuse(0, s, density, diff, dt, 4, N); advect(0, density, s, Vx, Vy, Vz, dt, N); } ``` 总结一下，该函数执行以下步骤： - 扩散所有三个速度分量。 - 修正速度，使其保持不可压缩性。 - 根据流体的速度移动速度本身（是不是有点晕了？）。 - 再次修正速度。 - 扩散染料。 - 根据速度移动染料。 以上就是基本原理。接下来我们将深入每个函数的实现。 ### set_bnd 如前所述，此函数设置立方体边缘的边界单元，使其完美抵消邻居的影响。这里有一点奇怪：这个 `b` 参数是什么？它可以取 0、1、2 或 3，每个值都有特殊含义，但这一点并不明显。答案在于哪些类型的数据可以传递给这个函数。我们有四种不同类型的数据（x、y、z 方向速度，以及密度）在流动，它们中的任何一个都可以传入 set_bnd。你可能已经注意到，这恰好是此参数可以取值的个数，这并非巧合。 让我们看看一个边界单元（糟糕的 ASCII 艺术警告）： ``` +---+---+ | |^ | | | / | | |/ | +---+---+ ``` 这里我们处于立方体的左边缘。左侧单元是需要抵消其邻居（右侧单元）的边界单元。右侧单元有一个向左上方向的速度。边界单元的 **x** 速度需要与邻居*相反*，但其 **y** 速度需要与邻居*相同*。这将产生如下结果： ``` +---+---+ | ^|^ | | / | / | |/ |/ | +---+---+ ``` 这会抵消导致流体穿过墙壁的运动，同时保留其余运动。因此，你可以看到采取什么行动取决于我们处理的是哪个数组；如果处理 x 速度，则必须将边界单元的值设置为与之邻居相反的值，而对于其他所有情况，则设置为相同值。这就是神秘参数 `b` 的作用。它告诉函数当前处理的是哪个数组，从而决定每个边界单元的值应该等于还是相反于其邻居的值。 此函数还处理角点。操作很简单：将每个角单元的值设置为其三个邻居的平均值。闲话少说，上代码： ```c static void set_bnd(int b, float *x, int N) { for(int j = 1; j < N - 1; j++) { for(int i = 1; i < N - 1; i++) { x[IX(i, j, 0 )] = b == 3 ? -x[IX(i, j, 1 )] : x[IX(i, j, 1 )]; x[IX(i, j, N-1)] = b == 3 ? -x[IX(i, j, N-2)] : x[IX(i, j, N-2)]; } } for(int k = 1; k < N - 1; k++) { for(int i = 1; i < N - 1; i++) { x[IX(i, 0 , k)] = b == 2 ? -x[IX(i, 1 , k)] : x[IX(i, 1 , k)]; x[IX(i, N-1, k)] = b == 2 ? -x[IX(i, N-2, k)] : x[IX(i, N-2, k)]; } } for(int k = 1; k < N - 1; k++) { for(int j = 1; j < N - 1; j++) { x[IX(0 , j, k)] = b == 1 ? -x[IX(1 , j, k)] : x[IX(1 , j, k)]; x[IX(N-1, j, k)] = b == 1 ? -x[IX(N-2, j, k)] : x[IX(N-2, j, k)]; } } x[IX(0, 0, 0)] = 0.33f * (x[IX(1, 0, 0)] + x[IX(0, 1, 0)] + x[IX(0, 0, 1)]); x[IX(0, N-1, 0)] = 0.33f * (x[IX(1, N-1, 0)] + x[IX(0, N-2, 0)] + x[IX(0, N-1, 1)]); x[IX(0, 0, N-1)] = 0.33f * (x[IX(1, 0, N-1)] + x[IX(0, 1, N-1)] + x[IX(0, 0, N)]); x[IX(0, N-1, N-1)] = 0.33f * (x[IX(1, N-1, N-1)] + x[IX(0, N-2, N-1)] + x[IX(0, N-1, N-2)]); x[IX(N-1, 0, 0)] = 0.33f * (x[IX(N-2, 0, 0)] + x[IX(N-1, 1, 0)] + x[IX(N-1, 0, 1)]); x[IX(N-1, N-1, 0)] = 0.33f * (x[IX(N-2, N-1, 0)] + x[IX(N-1, N-2, 0)] + x[IX(N-1, N-1, 1)]); x[IX(N-1, 0, N-1)] = 0.33f * (x[IX(N-2, 0, N-1)] + x[IX(N-1, 1, N-1)] + x[IX(N-1, 0, N-2)]); x[IX(N-1, N-1, N-1)] = 0.33f * (x[IX(N-2, N-1, N-1)] + x[IX(N-1, N-2, N-1)] + x[IX(N-1, N-1, N-2)]); } ``` ### lin_solve 如前所述，这个函数很神秘，但它执行某种求解。它通过遍历整个数组，将每个单元设置为邻居的某个组合值来实现。它会执行多次；迭代次数越多，结果越精确，但运行速度越慢。在上面的 step 函数中，使用了四次迭代。每次迭代后，它都会重置边界，以防止计算发散。代码如下： ```c static void lin_solve(int b, float *x, float *x0, float a, float c, int iter, int N) { float cRecip = 1.0 / c; for (int k = 0; k < iter; k++) { for (int m = 1; m < N - 1; m++) { for (int j = 1; j < N - 1; j++) { for (int i = 1; i < N - 1; i++) { x[IX(i, j, m)] = (x0[IX(i, j, m)] + a*( x[IX(i+1, j , m )] + x[IX(i-1, j , m )] + x[IX(i , j+1, m )] + x[IX(i , j-1, m )] + x[IX(i , j , m+1)] + x[IX(i , j , m-1)] )) * cRecip; } } } set_bnd(b, x, N); } } ``` ### diffuse Diffuse 非常简单；它只是预先计算一个值，然后将所有工作交给 lin_solve。也就是说，虽然我知道它的作用，但我并不真正了解其原理，因为所有工作都在那个神秘的函数中。代码如下： ```c static void diffuse (int b, float *x, float *x0, float diff, float dt, int iter, int N) { float a = dt * diff * (N - 2) * (N - 2); lin_solve(b, x, x0, a, 1 + 6 * a, iter, N); } ``` ### project 这个函数的具体工作原理也有些神秘，但它会进一步遍历数据并设置值，中间还穿插了对 lin_solve 的调用。代码如下： ```c static void project(float *velocX, float *velocY, float *velocZ, float *p, float *div, int iter, int N) { for (int k = 1; k < N - 1; k++) { for (int j = 1; j < N - 1; j++) { for (int i = 1; i < N - 1; i++) { div[IX(i, j, k)] = -0.5f*( velocX[IX(i+1, j , k )] - velocX[IX(i-1, j , k )] + velocY[IX(i , j+1, k )] - velocY[IX(i , j-1, k )] + velocZ[IX(i , j , k+1)] - velocZ[IX(i , j , k-1)] )/N; p[IX(i, j, k)] = 0; } } } set_bnd(0, div, N); set_bnd(0, p, N); lin_solve(0, p, div, 1, 6, iter, N); for (int k = 1; k < N - 1; k++) { for (int j = 1; j < N - 1; j++) { for (int i = 1; i < N - 1; i++) { velocX[IX(i, j, k)] -= 0.5f * ( p[IX(i+1, j , k )] - p[IX(i-1, j , k )]) * N; velocY[IX(i, j, k)] -= 0.5f * ( p[IX(i , j+1, k )] - p[IX(i , j-1, k )]) * N; velocZ[IX(i, j, k)] -= 0.5f * ( p[IX(i , j , k+1)] - p[IX(i , j , k-1)]) * N; } } } set_bnd(1, velocX, N); set_bnd(2, velocY, N); set_bnd(3, velocZ, N); } ``` ### advect 此函数负责实际移动物质。为此，它依次查看每个单元。在该单元中，它获取速度，沿着该速度回溯时间，查看最终落点。然后，它取落点周围单元值的加权平均，并将其应用于当前单元。代码如下： ```c static void advect(int b, float *d, float *d0, float *velocX, float *velocY, float *velocZ, float dt, int N) { float i0, i1, j0, j1, k0, k1; float dtx = dt * (N - 2); float dty = dt * (N - 2); float dtz = dt * (N - 2); float s0, s1, t0, t1, u0, u1; float tmp1, tmp2, tmp3, x, y, z; float Nfloat = N; float ifloat, jfloat, kfloat; int i, j, k; for(k = 1, kfloat = 1; k < N - 1; k++, kfloat++) { for(j = 1, jfloat = 1; j < N - 1; j++, jfloat++) { for(i = 1, ifloat = 1; i < N - 1; i++, ifloat++) { tmp1 = dtx * velocX[IX(i, j, k)]; tmp2 = dty * velocY[IX(i, j, k)]; tmp3 = dtz * velocZ[IX(i, j, k)]; x = ifloat - tmp1; y = jfloat - tmp2; z = kfloat - tmp3; if(x < 0.5f) x = 0.5f; if(x > Nfloat + 0.5f) x = Nfloat + 0.5f; i0 = floorf(x); i1 = i0 + 1.0f; if(y < 0.5f) y = 0.5f; if(y > Nfloat + 0.5f) y = Nfloat + 0.5f; j0 = floorf(y); j1 = j0 + 1.0f; if(z < 0.5f) z = 0.5f; if(z > Nfloat + 0.5f) z = Nfloat + 0.5f; k0 = floorf(z); k1 = k0 + 1.0f; s1 = x - i0; s0 = 1.0f - s1; t1 = y - j0; t0 = 1.0f - t1; u1 = z - k0; u0 = 1.0f - u1; int i0i = i0; int i1i = i1;