如何使用基于 Jos Stam 论文的模拟使液体停留在容器中？

问：

我正在根据 Jos Stam 的论文在 JS/WebGL 中进行 3D 流体模拟。我相信一切都在工作，除了边界条件、平流或超出我专业知识的东西。

我的模拟会扩散，对速度的变化做出响应，但在到达边界后会继续超出其容器。

产卵

运动

不保持液体

在这些图片中，透明部分是空的单元格和容器边界。

我担心我正在做的事情有根本性的错误。我希望看到流体生成、移动并与其容器交互。也许模拟的性质或我的网格是有缺陷的。如果我能提供更多细节，请告诉我。

GitHub 存储库

以下是我的边界、扩散和平流例程。我的网格是矩形的，尺寸为 N + 额外的 * N * N，而不是像 Stam 的论文那样的 （N + 1） ^3。因此，边界位于：

0, j, k // horizontal parallel to page
N + extra - 1, j, k

i, 0, k // vertical parallel to page
i, N - 1, k

i, j, 0 // perpendicular to page
i, j, N - 1


boundary(value, flag) {

        for (let i = 0; i < N; i++) {
            for (let j = 0; j < N; j++) {
                // yz parallel horizontal plane
                if (flag == 1) {
                    value[idx(0, i, j)] = -value[idx(1, i, j)];
                    value[idx(N + extra - 1, i, j)] = -value[idx(N + extra - 2, i, j)];
                } else {
                    value[idx(0, i, j)] = value[idx(1, i, j)];
                    value[idx(N + extra - 1, i, j)] = value[idx(N + extra - 2, i, j)];
                }
            }
        }

        for (let i = 0; i < N + extra; i++) {
            for (let j = 0; j < N; j++) {
                // xz parallel vertical plane
                if (flag == 2) {
                    value[idx(i, 0, j)] = -value[idx(i, 1, j)];
                    value[idx(i, N - 1, j)] = -value[idx(i, N - 2, j)];
                } else {
                    value[idx(i, 0, j)] = value[idx(i, 1, j)];
                    value[idx(i, N - 1, j)] = value[idx(i, N - 2, j)];
                }

                // xy perpendicular plane
                if (flag == 3) {
                    value[idx(i, j, 0)] = -value[idx(i, j, 1)];
                    value[idx(i, j, N - 1)] = -value[idx(i, j, N - 2)];
                } else {
                    value[idx(i, j, 0)] = value[idx(i, j, 1)];
                    value[idx(i, j, N - 1)] = value[idx(i, j, N - 2)];
                }
            }
        }

        value[idx(0, 0, 0)] = .33 * (value[idx(1, 0, 0)] + value[idx(0, 1, 0)] + value[idx(0, 0, 1)]);
        value[idx(0, N - 1, 0)] = .33 * (value[idx(1, N - 1, 0)] + value[idx(0, N - 2, 0)] + value[idx(0, N - 1, 1)]);

        value[idx(N + extra - 1, 0, 0)] = .33 * (value[idx(N + extra - 2, 0, 0)] + value[idx(N + extra - 1, 1, 0)] + value[idx(N + extra - 1, 0, 1)]);
        value[idx(N + extra - 1, N - 1, 0)] = .33 * (value[idx(N + extra - 2, N - 1, 0)] + value[idx(N + extra - 1, N - 2, 0)] + value[idx(N + extra - 1, N - 1, 1)]);
        
        value[idx(0, 0, N - 1)] = .33 * (value[idx(1, 0, N - 1)] + value[idx(0, 1, N - 1)] + value[idx(0, 0, N - 2)]);
        value[idx(0, N - 1, N - 1)] = .33 * (value[idx(1, N - 1, N - 1)] + value[idx(0, N - 2, N - 1)] + value[idx(0, N - 1, N - 2)]);
        
        value[idx(N + extra - 1, 0, N - 1)] = .33 * (value[idx(N + extra - 2, 0, N - 1)] + value[idx(N + extra - 1, 1, N - 1)] + value[idx(N + extra - 1, 0, N - 2)]);
        value[idx(N + extra - 1, N - 1, N - 1)] = .33 * (value[idx(N + extra - 2, N - 1, N - 1)] + value[idx(N + extra - 1, N - 2, N - 1)] + value[idx(N + extra - 1, N - 1, N - 2)]);

        return value;
    }


    linear_solve(flag, value, value_old, rate, c) {
        for (let k = 0; k < 10; k++) {
            for (let x = 1; x < N + extra - 1; x++) {
                for (let y = 1; y < N - 1; y++) {
                    for (let z = 1; z < N - 1; z++) {
                        let sum =   value[idx(x - 1, y, z)] + value[idx(x + 1, y, z)] +
                                    value[idx(x, y - 1, z)] + value[idx(x, y + 1, z)] +
                                    value[idx(x, y, z - 1)] + value[idx(x, y, z + 1)];

                        value[idx(x, y, z)] = (value_old[idx(x, y, z)] + rate * sum) / c;
                    }
                }
            }
        }

        value = this.boundary(value, flag);
        return value;

    }

    diffuse(flag, value, value_old, viscosity, dt) {
        let rate = dt * viscosity * (N + extra) * N**2;

        let c = 1 + 6 * rate;

        value = this.linear_solve(flag, value, value_old, rate, c)
        return value;

    }

    advect(flag, value, value_old, dt) {
        let dtOX = dt * (N + extra);
        let dtO = dt * N;

        for (let i = 1; i < N + extra - 1; i++) {
            for (let j = 1; j < N - 1; j++) {
                for (let k = 1; k < N - 1; k++) {

                    let x = i - dtOX * this.u[idx(i, j, k)];
                    let y = j - dtO * this.v[idx(i, j, k)];
                    let z = k - dtO * this.w[idx(i, j, k)];

                    if (x < 0.5) {
                        x = 0.5;
                    }
                    if (y < 0.5) {
                        y = 0.5;
                    }
                    if (z < 0.5) {
                        z = 0.5;
                    }


                    if (x > N + extra - 1.5) {
                        x = N + extra - 1.5;
                    }
                    if (y > N - 1.5) {
                        y = N - 1.5;
                    }
                    if (z > N - 1.5) {
                        z = N - 1.5;
                    }

                    let i0 = Math.floor(x);
                    let j0 = Math.floor(y);
                    let k0 = Math.floor(z);

                    let i1 = i0 + 1;
                    let j1 = j0 + 1;
                    let k1 = k0 + 1;

                    let s1 = x - i0;
                    let t1 = y - j0;
                    let u1 = z - k0;

                    let s0 = 1 - s1;
                    let t0 = 1 - t1;
                    let u0 = 1 - u1;

                    value[idx(i, j, k)] = 
                    
                    s0 * (
                            t0 * u0 * value_old[idx(i0, j0, k0)] +
                            t1 * u0 * value_old[idx(i0, j1, k0)] +
                            t0 * u1 * value_old[idx(i0, j0, k1)] +
                            t1 * u1 * value_old[idx(i0, j1, k1)]
                        ) +

                    s1 * (
                            t0 * u0 * value_old[idx(i1, j0, k0)] +
                            t1 * u0 * value_old[idx(i1, j1, k0)] +
                            t0 * u1 * value_old[idx(i1, j0, k1)] +
                            t1 * u1 * value_old[idx(i1, j1, k1)]
                        );
                }
            }
        }
        value = this.boundary(value, flag);
        return value;
    }