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

Question

我最近遇到了这个流体动力学模拟，为了更好地理解它，我一直在尝试阅读晶格玻尔兹曼方法。那里的大多数代码都是非常透明的，所以在阅读代码和阅读维基百科之间，我几乎已经了解了一切。。。除了核心"碰撞"函数中的运动学计算之外。

我已经能够计算出1/9、4/9和1/36的因子与沿着不同晶格方向连接细胞中心的向量的长度有关，我甚至可以找到解释不同晶格类型使用哪些不同因子的资源(本代码中使用D2Q9晶格)。但我还没有找到任何东西来解释你是如何从定义Lattice Boltzmann算法碰撞步骤的通用矢量方程到下面显示的具体九行实现算法的。特别是，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#编写的)从根据粘度计算弛豫频率ω开始，该粘度是通过强制执行某个雷诺数Re:=U L/nu设置的。(通常，这是一个常数，因此也可以提前计算，而不是在每个循环中计算。)

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

• 然后代码在整个域上循环，忽略每个方向的第一个和最后一个节点，这将是由不同函数处理的边界条件

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

• 对于种群的存储，它依赖于每个晶格方向的全局1D阵列\alpha，因此转向线性索引来选择正确的节点。在这种情况下，其余节点的方向命名为0，其他种群则根据相应的方向(nord、east等)命名。通常，代码使用编号约定

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_\alphae_{i\alpha}f\alpha来确定速度。在这种情况下，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分别是水平/垂直和对角线总体的权重和密度的组合。这里令人困惑的是，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_\aalpha+\omega*(f_\albha^{eq}-f_\alpha)。每一行都负责一个离散方向\alpha，因此对于向量方程的单个条目也是如此。注：

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

  • 不可压缩平衡函数由f_\alpha^{eq}=w_\alph\rho(1+e_{i\alpha}u_i/c_s^2+(e_i\alpha u_i)^2/(2c_s^4)-u_i^2/(2 c_s^2))给出，其中我们必须对包含索引i的所有项进行两次求和。声的晶格速度c_s由1/\sqrt(3)*dx/dx给出，因此f_\alpha^｛eq｝=w_\alph\rho(1+3e_{i\alpha}u_i+4.5(e_i\alpha u_i)^2-1.5 u_i^2)。术语one9thrho和one36thrho对应于括号前的术语。ux3和uy3的和对应于括号3e_｛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使用两个数组，并在一步中执行冲突和流式处理，而不是两步。此外，平衡函数可以进一步重写，通过组合ω和oneXXthrho来重新计算尽可能少的项，并进一步重写二次项。请参阅"格子Boltzmann方法：原理与实践"中的代码，以获得更好的代码。这应该已经将性能提高了至少两倍。为了在您的机器上模拟3D，您必须应用几个更困难的优化。此外，还有更好的论坛来讨论这个话题，例如Palabos(日内瓦大学)和OpenLB(卡尔斯鲁厄理工学院)。