格子玻尔兹曼模拟中幻数的解释？

Question

我最近遇到了 this fluid dynamics simulation 并且一直在尝试阅读格子玻尔兹曼方法以便更好地理解它。那里的大部分代码都非常透明，所以在阅读代码和阅读维基百科之间，我几乎已经弄明白了一切的作用……除了核心“碰撞”函数中的运动学计算。

我已经能够弄清楚1/9、4/9和1/36的因数与沿不同晶格方向连接单元格中心的向量的长度有关，我什至可以找到解释不同晶格类型使用哪些不同因素的资源（此代码中使用了 D2Q9 晶格）。但我还没有找到任何内容来解释如何从定义 Lattice Boltzmann Algorithm 的碰撞步骤的通用向量方程式到下面所示的具体九行实现算法。特别是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]);
    }
  }
}

Answer 1

代码（似乎是用 C# 编写的）开始于从粘度 计算弛豫频率 omega，这是通过强制执行某个雷诺数 Re:= U L /努。（通常这是一个常数，因此也可以提前计算，而不是在每个循环中计算。）

var omega = 1 / (3*viscosity + 0.5);

然后代码遍历整个域，省略每个方向的第一个和最后一个节点，这将是将由不同函数处理的边界条件。

for (var y=1; y<ydim-1; y++) {
   for (var x=1; x<xdim-1; x++) {`

为了存储种群，它依赖于每个晶格方向的全局一维数组 \alpha，因此转向线性索引以select 正确的节点。在这种情况下，其余节点的方向被命名为 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} [=86 来确定速度 =] 用于当前格子。在这种情况下，x 指向水平方向（东西），y 指向垂直方向（南北），其中东和北是正方向。这意味着 \alpha = {北，东北，西北，南，东南，西南}。

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是weights和one36thrho的组合horizontal/vertical 和对角线种群的密度分别为。这里令人困惑的是 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 中，因此足以 increment (+=) f_\alpha（数组 nNWSE）由 \omega*(f_\alpha^{eq} - f_\alpha)。括号中的第一项是平衡函数，而最后一项对应当前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) 对应最后一项。您必须考虑方向 i 上晶格速度 \vec e_\alpha 的相应速度投影 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 重新计算尽可能少的项，并进一步重写二次项。参考code that accompanies "The Lattice Boltzmann Method: Principles and Practice" for better written code. This should already increase performance by at least a factor of two. In order to simulate 3D on your machine you will have to apply several more difficult optimisations. Additionally there are better forums for this topic, for instance Palabos (University of Geneva) and OpenLB（卡尔斯鲁厄理工学院）

格子玻尔兹曼模拟中幻数的解释？

Explanation for magic numbers in lattice Boltzmann simulation?

simulation

physics

mathematical-lattices