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

NoiseFloor · 发表于 2022-4-23 17:02

本次要复现的是一篇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 · 发表于 2022-4-23 17:03

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

johnsoncodehk · 发表于 2022-4-23 17:06

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

		自动登录	找回密码
密码			立即注册

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

本帖子中包含更多资源