simulation - 格子ボルツマンシミュレーションにおける魔法数の説明?

Question

私は最近、この流体力学シミュレーションに出くわし、それをよりよく理解するために格子ボルツマン法について調べようとしています。そこにあるコードのほとんどはかなり透明なので、コードを読んでウィキペディアを読む間に、私はすべてが何をしているのかほとんど理解できました...コアの「衝突」関数の運動学的計算を除いて。

1/9、4/9、および 1/36 の因数は、異なる格子方向に沿ってセルの中心を結ぶベクトルの長さに関連していることを突き止めることができました。異なるラティス タイプに使用する係数 (このコードでは D2Q9 ラティスが使用されています)。しかし、格子ボルツマン アルゴリズムの衝突ステップを定義する一般的なベクトル方程式から、以下に示す特定の 9 行の実装演算に至るまでの過程を説明するものを見つけることができませんでした。特に、3、1.5、4.5 という係数はどこから来るのでしょうか?

リンクされた Web ページで使用されているコードは次のとおりです (自由変数を削除して読みやすくするために少し編集しています)。

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]);
    }
  }
}