图形学论文复现-Vortex Fluid for Gaseous …

NoiseFloor

本次要复现的是一篇2005年的论文，Vortex Fluid for Gaseous Phenomena。复现了2D版本，编程语言：C++
这篇文章使用涡粒子法来模拟烟雾，引入了一个vortex sheet边界层模型来处理边界条件。实现之前可以先看看大佬们的一些文章:
以及games201第六讲。robert brison那本流体书也有一章讲涡方法的。
常见的流体模拟方法比如sph这种都是以速度场作为追踪变量，涡方法解的则是N-S方程的旋度版，以涡量作为追踪变量。
N-S方程：



对上式两边作旋度运算，然后省略亿点点推导，得到一个旋度版的方程：



这篇文章就是在围绕求解这个旋度版的方程展开的。本文使用拉格朗日粒子来离散流体元素，因此对流项会被隐式的求解。
算法步骤：使用split算法求解这个方程。

• 在网格上用有限差分计算外力的旋度，更新粒子的涡量为\omega_1
  
• 计算粘性项
  
• 计算Stretching项
  
• 通过Biot-Savart法则从涡量场中恢复出速度场
  
• 处理边界问题（无穿透边界和无滑移边界）
  

1.外力的计算
烟雾模拟的主导力是空气浮力，受空气密度和温度的影响。这里没有实现浮力场的建模，会在下次去复现浮力场建模的论文。

2.粘性项
暂时没有复现（感觉这个粘性不是很重要，先偷个懒。。。）。这里大概是用一种随机分布函数去对涡粒子进行扰动。

3.Stretching计算
2d情况下Stretching直接为0，不用算，真是太高兴了。

4.速度场的恢复
每一个涡粒子都有自己的一个涡量，涡量的计算公式：



里面那个东西的表达式：



其中D是维数，\sigma是一个自己指定的光滑系数。然后有一些细节推导就不写了，通过下面这个公式可以把速度场恢复出来：



这几个公式都很好离散。对于2D版本的涡粒子法，如果不考虑外力，其涡量的来源只有两个：1，手动指定；2，边界的扰动。到这里其实就足够复现出这样的效果了（论文截图）：



浮力场的建模会在下一篇文章里复现。这里可以写个vortex ring leapfrogging效果。（详见最上面的文章链接）



可以把涡量为0的粒子叫做tracer粒子，涡量不为0的叫做涡粒子。初始化4个涡粒子（大圆），手动指定其涡量，然后若干tracer粒子（灰色区域），然后让它进行演化。

这里先提供一个伪代码，完整代码我正在整理，有一些命名不规范的地方：
//公式10的离散
Vector2D computeUSingle(const Vector2D& pos, int i)const {
    auto position = particles->positions();
    auto gamma = particles->gamma();
    auto r2 = (pos - position).getLengthSquared();
    auto uv = Vector2D(position.y - pos.y, pos.x - position.x);
    return gamma * uv / (PI * r2) * 0.5 * (1.0 - pow(kE, -r2 / (eps * eps)));
}

//初始化涡粒子
void FoamVortexSolver::emitParticles() {
    auto& n = data->numberOfParticles();
    auto& pos = data->positions();
    auto& gamma = data->gamma();

    n = 4;
    Vector2D A(0, 1.2);
    Vector2D B(0, 1.1);
    Vector2D C(0, 1.0);
    Vector2D D(0, 0.8);

    pos.push(A);
    pos.push(B);
    pos.push(C);
    pos.push(D);

    gamma.push(1.0);
    gamma.push(1.0);
    gamma.push(-1.0);
    gamma.push(-1.0);

}

//时间积分
int n = particle->numberOfParticles()；
auto position = particle->positions();
for(int i = 0;i<n;++i){
   Vector2D vortexVel;
   for(int j = 0;j<n;++j){
        if(i!=j){
          vortexVel+=computeUSingle(positions, j);
    }
      
    newPositons = positions+dt*vortexVel;
}
positions=newPositons;

看一下效果：


再增加一些tracer粒子：



<b>5.处理边界问题（无穿透边界和无滑移边界）
涡方法的边界处理比较特殊，因为它是追踪涡量的，所以不能简单的对速度场作处理。这篇文章引入了一个边界层模型，对于无粘流来说，流体会光滑的滑过边界，这时只需要保证在边界附近的速度场沿着法线的分量为0，这样流体就不会穿透边界。本文使用一个vortex sheet boundary模型，假设边界上分布着无数个涡量，它在空间中会引发一个场，这个场是由vortex sheet的\Gamma来决定的，因此就可以通过联立方程，求出一个\Gamma来使得它能够消去法向分量。

首先是边界的建模，这里以一个绕柱流为例，2D下的圆是用多边形逼近的：
RegularPolygon::RegularPolygon(int edgeNum_, Vector2D center_, double r_) :
    _edgeNum(edgeNum_),
    _center(center_),
    _r(r_) {

     //_data里是多边形的每条边
    _data.resize(edgeNum_);

    double ca = 0;
    auto aiv = 360.0 / edgeNum_;
    auto ata = kPiD / 180;
    vector<Vector2D> list(edgeNum_);
    for (int k = 0; k < edgeNum_; k++) {
        auto x = cos(ca * ata) * r_ + center_.x;
        auto y = sin(ca * ata) * r_ + center_.y;
        Vector2D temp(x, y);
        list[k] = temp;
        ca += aiv;
    }
    int i = 0;
    SurfaceElement2 temp1;
    for (; i < edgeNum_ - 1; ++i) {
        temp1.start = list;
        temp1.end = list[i + 1];
        auto mid = (temp1.start + temp1.end) * 0.5;
        temp1.normal = (mid - center_).getNormalize();
        _data = temp1;
    }

    temp1.start = list[0];
    temp1.end = list;
    auto mid = (temp1.start + temp1.end) * 0.5;
    temp1.normal = (mid - center_).getNormalize();
    _data = temp1;
   
}

这里用panel method来对vortex sheet进行数值离散，即每条线段是一个panel，然后其中点代表panel上的所有点。
边界会在空间中引发一个场，先考虑单个panel对粒子的影响：



它做了一个坐标变换，首先把世界坐标系下的点变到这个局部坐标系中，变换代码：
Vector2D vel_to_world(const Vector2D vel_local,  const Vector2D normal_i) {
    double cos_theta = normal_i.y;
    double sin_theta = normal_i.x;
    Vector2D vel = Vector2D(cos_theta * vel_local.x + sin_theta * vel_local.y,
        -sin_theta * vel_local.x + cos_theta * vel_local.y);
    return vel;
}
然后在局部坐标系下计算由panel引入的速度，注意计算出来后要把结果变回到世界坐标系：



Vector2D FoamVortexSolver::computeSingleVelocityFromPanels(int index) {
    auto& pos = _foamVortexData->positions();
    auto panel = _foamVortexData->panelSet;
    auto panelSize = panel->size();
    auto& gama1 = _foamVortexData->strength;

    Vector2D result;

    for (int i = 0; i < panelSize; ++i) {

        auto start = panel->lookAt(i).start;
        auto end = panel->lookAt(i).end;
        auto normal = panel->lookAt(i).normal;

        Vector2D temp2;

        double beta = 0.0;
        auto vec_r1 = start - pos[index];
        auto vec_r2 = end - pos[index];
        auto r1 = vec_r1.getLength();
        auto r2 = vec_r2.getLength();
        auto temp1 = vec_r1.dot(vec_r2) / (r1 * r2);
        beta = acos(std::min(std::max(temp1, -1.0), 1.0));
        if (isnan(beta)) {
            beta = kPiD;
        }

        temp2.x = beta / (2.0 * kPiD);
        temp2.y = log((r2 + fv_eps) / (r1 + fv_eps)) / (2.0 * kPiD);

        double tempgamma = 0.0;
        if (gama1.size() > 0)
            tempgamma = gama1;
        temp2 = temp2 * tempgamma;
        temp2 = vel_to_world(temp2, start, normal);
        result += temp2;
    }

    return result;
}

然后通过这个公式，把所有panel对粒子的贡献加起来，就是整个边界引入的速度场：



但现在\Gamma是未知的，因此要联立一个方程组，求解出满足no-through边界条件的\Gamma来:



这个就是no-through条件的约束，这是一个N*N的稠密矩阵求解：



然后还需要满足一个条件：



这样的话就成了一个（N+1）*N的方程组求解，我这里使用的是Eigen的最小二乘求解。
void FoamVortexSolver::computeBoundaryMatrix() {
    auto panels = _foamVortexData->panelSet;
    auto panelSize = panels->size();
    Eigen::MatrixXd& A = _foamVortexData->A;
    Eigen::MatrixXd& B = _foamVortexData->B;

    auto sizex = panelSize + 1;
    auto sizey = panelSize;
    A.resize(sizex, sizey);
    B.resize(sizex, sizey);

    for (int j = 0; j < panelSize; ++j) {
        auto normal = panels->lookAt(j).normal;
        for (int i = 0; i < panelSize; ++i) {
            auto mid_j = panels->midPoint(j);
            auto u_ji = computeUnitVelocityFromPanels(i, mid_j).dot(normal);
            A(j, i) = -u_ji;
        }
    }

    for (int i = 0; i < panelSize; ++i) {
        A(panelSize, i) = 1.0;
    }

   
}

对于静止的刚体边界，这个矩阵是不变，只需要每次组装新的b即可：
void FoamVortexSolver::vortexSheetSolve(double timeIntervalInSeconds) {
    auto data = _foamVortexData;
    auto n = data->numberOfParticles();
    auto panels = data->panelSet;
    auto panelSize = panels->size();
    auto& pos = data->positions();
    Eigen::MatrixXd& A = _foamVortexData->A;
    Eigen::VectorXd& x = _foamVortexData->strength;

    //组装b
    Eigen::VectorXd b(panelSize + 1);
    for (int i = 0; i < panelSize; ++i) {
        auto normal = panels->lookAt(i).normal;
        auto pos = panels->midPoint(i);
        auto vec = vs_vec;

        Vector2D temp;
        //for (int j = 0; j < n; ++j) {
            //temp += computeUSingle(pos, j);
        //}
        //vec是背景流场，temp是涡粒子在边界上引发的速度场，no-through边界条件无需计算这一项
        //no-slip是需要加上这一项的
        b = (vec + temp).dot(normal);
    }

    b[panelSize] = 0;

    x = A.colPivHouseholderQr().solve(b);

}
这样就可以把\Gamma求出来了，然后用这个求出边界附近正确的速度分布。看下效果：



流体光滑的滑过了边界。
对于无滑移边界条件，和上面是同样的道理，只是建立的是对速度场切向分量的约束。解出一个消去切向分量的\Gama后，这里使用的是random walk的方法。即每隔几个时步，就从每个panel的中点向空间中发射一个涡粒子，这个涡粒子的涡量就是解出来的\Gamma计算得到的。然后这些涡粒子就会从边界出来，搓出一个卡门漩涡来：



这里它说这个\alpha是panel的长度，我傻傻的信了他的鬼话，用panel两个端点相减得到一个panelLength，结果怎么也做不出卡门漩涡来，后来发现他这里只是做了一个简化的处理，只是粗略的给涡量乘上一个系数，谁知道他的单位是啥，我直接给涡量乘了30，果然卡门漩涡就出来了:



这种做法其实不是很精确，边界附近的涡量好像没了。下一篇论文里会使用一种更为精确的边界处理。
代码我最近会整理一下然后上传，里面有些我自己的实验需要分开一下。
总结：这里面其实有一个很不理解的地方，就是边界的坐标系到底该怎么定，按照图上给的，每一个坐标系都是以panel的起始点为原点，end-start方向为x轴，法线方向为y轴，应该是需要保证y.cross(x)>0才行。但是实际编程的时候我发现这样做算出来是错的，得搞成<0才行，但也不能全部<0,最后一个闭合的panel那个的y.cross(x)>0才可以。有空再思考一下这个问题。

请大家多多批评指正！

IT圈老男孩1

粘性项我没理解错的话，应该就是求一下vorticity的二阶导数乘以粘性系数就行吧。你说的random扰动是想实现烟雾的随机扩散吗？加油更新，共同学习[赞同][赞同]

johnsoncodehk

你说粘性项里的扰动吗？论文里说是用了一种高斯分布函数，避免向空间中发射粒子，我没有细看[捂脸]