0

我最近遇到了这种流体动力学模拟，并一直在尝试阅读格子玻尔兹曼方法，以便更好地理解它。那里的大部分代码都非常透明，所以在阅读代码和阅读 Wikipedia 之间，我已经大致了解了所有内容……除了核心“碰撞”函数中的运动学计算。

我已经能够弄清楚 1/9、4/9 和 1/36 的因子与沿不同晶格方向连接单元中心的向量的长度有关，我什至可以找到解释有什么不同的资源用于不同晶格类型的因素（此代码中使用了 D2Q9 晶格）。但是我还没有找到任何东西来解释你如何从定义格子玻尔兹曼算法的碰撞步骤的通用向量方程到下面显示的具体的九行实现算法。特别是 3、1.5 和 4.5 的那些因数是从哪里来的？

链接网页中使用的代码（稍作修改以删除自由变量并提高可读性）如下：

function collide(viscosity) { // kinematic viscosity coefficient in natural units
  var omega = 1 / (3*viscosity + 0.5); // reciprocal of relaxation time
  for (var y=1; y<ydim-1; y++) {
    for (var x=1; x<xdim-1; x++) {
      var i = x + y*xdim; // array index for this lattice site
      var thisrho = n0[i]+nN[i]+nS[i]+nE[i]+nW[i]+nNW[i]+nNE[i]+nSW[i]+nSE[i];

      // macroscopic horizontal velocity
      var thisux = (nE[i]+nNE[i]+nSE[i]-nW[i]-nNW[i]-nSW[i])/thisrho;

      // macroscopic vertical velocity
      var thisuy = (nN[i]+nNE[i]+nNW[i]-nS[i]-nSE[i]-nSW[i])/thisrho;

      var one9thrho = thisrho / 9;
      var one36thrho = thisrho / 36;
      var ux3 = 3 * thisux;
      var uy3 = 3 * thisuy;
      var ux2 = thisux * thisux;
      var uy2 = thisuy * thisuy;
      var uxuy2 = 2 * thisux * thisuy;
      var u2 = ux2 + uy2;
      var u215 = 1.5 * u2;

      n0[i]  += omega * (4/9*thisrho * (1                              - u215) - n0[i]);
      nE[i]  += omega * (  one9thrho * (1 + ux3       + 4.5*ux2        - u215) - nE[i]);
      nW[i]  += omega * (  one9thrho * (1 - ux3       + 4.5*ux2        - u215) - nW[i]);
      nN[i]  += omega * (  one9thrho * (1 + uy3       + 4.5*uy2        - u215) - nN[i]);
      nS[i]  += omega * (  one9thrho * (1 - uy3       + 4.5*uy2        - u215) - nS[i]);
      nNE[i] += omega * ( one36thrho * (1 + ux3 + uy3 + 4.5*(u2+uxuy2) - u215) - nNE[i]);
      nSE[i] += omega * ( one36thrho * (1 + ux3 - uy3 + 4.5*(u2-uxuy2) - u215) - nSE[i]);
      nNW[i] += omega * ( one36thrho * (1 - ux3 + uy3 + 4.5*(u2-uxuy2) - u215) - nNW[i]);
      nSW[i] += omega * ( one36thrho * (1 - ux3 - uy3 + 4.5*(u2+uxuy2) - u215) - nSW[i]);
    }
  }
}
4

1 回答 1

2
  • 代码（似乎是用 C# 编写的）首先从通过强制某个雷诺数 Re:= UL / nu 设置的粘度计算弛豫频率 omega 。（通常这是一个常数，因此也可以提前计算，而不是在每个循环中计算。）
var omega = 1 / (3*viscosity + 0.5);
  • 然后代码在整个域上循环，忽略每个方向的第一个和最后一个节点，这将是由不同函数处理的边界条件。
for (var y=1; y<ydim-1; y++) {
   for (var x=1; x<xdim-1; x++) {`
  • 对于种群的存储，它依赖于每个晶格方向α的全局一维数组，因此转向线性索引来选择正确的节点。在这种情况下，其余节点的方向命名为 0，其他种群根据相应的方向（北、东等）命名。通常代码使用编号约定来代替。
var i = x + y*xdim;
  • 它使用线性索引读取种群并将它们相加以获得密度（零阶动量）。
var thisrho = n0[i]+nN[i]+nS[i]+nE[i]+nW[i]+nNW[i]+nNE[i]+nSW[i]+nSE[i];
  • 然后它通过对当前晶格强制执行一阶动量 \rho*u_i = \sum_\alpha e_{i \alpha} f_\alpha 来确定速度。在这种情况下，x 指向水平方向（东西方向），y 指向垂直方向（南北方向），其中东和北是正方向。这意味着 α = {北，东北，西北，南，东南，西南}。
var thisux = (nE[i]+nNE[i]+nSE[i]-nW[i]-nNW[i]-nSW[i])/thisrho;

var thisuy = (nN[i]+nNE[i]+nNW[i]-nS[i]-nSE[i]-nSW[i])/thisrho;
  • 然后声明了几个变量，将在下一节中重复使用：one9thrho 和 one36thrho 分别是水平/垂直和对角种群的权重和密度的组合。这里令人困惑的是，ux3 表示 3*ux，而 ux2 旨在表示 ux^2。
var one9thrho = thisrho / 9;
var one36thrho = thisrho / 36;
var ux3 = 3 * thisux;
var uy3 = 3 * thisuy;
var ux2 = thisux * thisux;
var uy2 = thisuy * thisuy;
var uxuy2 = 2 * thisux * thisuy;
var u2 = ux2 + uy2;
var u215 = 1.5 * u2;

  • 最后一步处理碰撞f_\alpha^t（其中 t 代表临时性，因为它将被流式跟进）= f_\alpha + \omega*(f_\alpha^{eq} - f_\alpha)。每条线负责一个离散方向 α，因此对于向量方程的单个条目。笔记：

    • 在这种情况下，f_\alpha^t 再次存储在 f_\alpha 中，因此通过 \omega*(f_\alpha^{eq} - f_\alpha递增(+=) f_\alpha（数组 nNWSE）就足够了）。括号中的第一项是平衡函数，而最后一项对应于当前的 f_\alpha。

    • 不可压缩平衡函数为 f_\alpha^{eq} = w_\alpha \rho (1 + e_{i \alpha} u_i/c_s^2 + (e_i \alpha u_i)^2/(2 c_s^4) - u_i^2/(2 c_s^2)) 我们必须对包含索引 i 的所有项求和两次. 声速 c_s 由 1/\sqrt(3)*dx/dx 给出，因此 f_\alpha^{eq} = w_\alpha \rho (1 + 3 e_{i \alpha} u_i + 4.5 (e_i \alpha u_i)^2 - 1.5 u_i^2)。术语 one9thrho 和 one36thrho 对应于括号前的术语。ux3 和 uy3 之和对应于括号内的第二项 3 e_{i \alpha}。倒数第二项 4.5 (e_i \alpha u_i)^2 对应于水平或垂直方向的 4.5 u_x^2 或 4.5 u_y^2 和对角线方向的 4.5 (+-u_x +- uy)^2 作为两个分量存在，因此导致 4.5 (u_x^2 + u_y^2 +- u_x u_y)。减去 u215 = 1.5*u2 = 1.5*(ux+uy) 对应最后一项。您必须考虑晶格速度 \vec e_\alpha 在 i 方向上的对应速度投影 e_{i \alpha} 是 0 还是 +-1。

n0[i]  += omega * (4/9*thisrho * (1                              - u215) - n0[i]);
nE[i]  += omega * (  one9thrho * (1 + ux3       + 4.5*ux2        - u215) - nE[i]);
nW[i]  += omega * (  one9thrho * (1 - ux3       + 4.5*ux2        - u215) - nW[i]);
nN[i]  += omega * (  one9thrho * (1 + uy3       + 4.5*uy2        - u215) - nN[i]);
nS[i]  += omega * (  one9thrho * (1 - uy3       + 4.5*uy2        - u215) - nS[i]);
nNE[i] += omega * ( one36thrho * (1 + ux3 + uy3 + 4.5*(u2+uxuy2) - u215) - nNE[i]);
nSE[i] += omega * ( one36thrho * (1 + ux3 - uy3 + 4.5*(u2-uxuy2) - u215) - nSE[i]);
nNW[i] += omega * ( one36thrho * (1 - ux3 + uy3 + 4.5*(u2-uxuy2) - u215) - nNW[i]);
nSW[i] += omega * ( one36thrho * (1 - ux3 - uy3 + 4.5*(u2+uxuy2) - u215) - nSW[i]);

像这样的代码不太可能非常有效。为了加快代码速度，您应该为 f_\alpha 和 f_\alpha^t 使用两个数组，并一步而不是两步执行碰撞和流式传输。此外，可以通过结合 omega 和 oneXXthrho 进一步重写平衡函数以重新计算尽可能少的项，并进一步重写二次项。请参阅《格子玻尔兹曼方法：原理与实践》随附的代码，以获得更好的编写代码。这应该已经将性能提高了至少两倍。为了在您的机器上模拟 3D，您将不得不应用几个更困难的优化。此外，这个主题还有更好的论坛，例如Palabos（日内瓦大学）和OpenLB（卡尔斯鲁厄理工学院）。

于 2019-10-25T08:01:12.873 回答