5.5.4 代码分析
参考5.4节的代码，这里重点对有变动的代码进行注释说明。
第一部分，参数初始化：
这部分代码主要设置通道的尺寸，定义雷诺数Re，根据雷诺数和特征速度，以10为特征长度(小圆柱直径)，可以计算得到运动黏度nu，进一步获得松弛因子τ的倒数omega。第7～10行定义格子参数，包括声速cs、方向权重t，方向向量ex和ey。
第二部分，变量初始化：
第15行，meshgrid()是MATLAB中用于生成网格坐标序号的函数，函数格式为[X，Y]=meshgrid(xgv，ygv)。meshgrid()函数生成的是X、Y大小相等的矩阵，xgv、ygv是两个行向量。X的值通过将xgv复制length(ygv)行(严格意义上是length(ygv)-1行)得到；Y的值通过先对ygv进行转置得到ygv'，再将ygv'复制(length(xgv)-1)次得到。
第17行，对数组中所有元素进行并行操作。代码“bnds=bnds+(((x'-Cen_r(i，1)).^2+(y'-Cen_r(i，2)).^2)＜=Cen_r(i，3)^2)；”是 一 个“比较”操作，对于满足条件的元素返回结果1，若不满足条件则返回结果0。举一个例子，令“A=[1 2 4 5 7 8]”，那么“M=(A＞=5)”返回的结果即“M=[0 0 0 1 1 1]”。本行代码的内涵是计算所有点距圆心的距离，将距离小于等于半径的点标识为1，否则标识为0。由于本例中设置了两个圆柱(两个圆柱无公共点)，因此可以用累加方式“bnds=bnds+…”获得两个圆柱所有节点标识的结果。
第19行，直接将通道的上下壁面处节点均标识为1，这样就设置了固体壁面。
第20行，把二维向量转换为一维向量，并获得固体节点在一维向量中的位置序号。举个例子，令
若要返回4的整倍数的元素位置给变量b，那么首先需获得标识矩阵，其MATLAB代码为“b=mod(a，4)==0”，即
之后，采用“[bx，by]=find(b)”获得值为1的元素位置坐标，可得bx=2，3，4，by=1，2，3。其中bx为元素1的行号、by对相应列号。但如果返回值只有一个，即采用代码“obs=find(b)”，则返回的结果obs为2、8、14。这个结果如何来的呢?看看图5.19中箭头指示方向，沿着指示的方向就可以发现，在第2、8、14个位置的值为1。因此在MATLAB中，采用的是“列优先”的方式进行二维数组元素的排序。
同时应注意到，此处的obs和bx、by之间有换算关系“obs=(by-1)*5+bx”。有了以上的基础，就可以更好地理解第23～30行中关于m和mb的操作。
第21行，给定反弹序列，导致在反弹操作前后，节点的概率密度函数方向发生如图5.20所示的变化。
图5.19 二维数组一维化的列优先原则示意
图5.20 反弹操作前后概率密度函数的分布
(a)反弹前；(b)反弹后
因此，对于流场内部的固体节点，反弹是在其节点上给了一个阻止流动的趋势，从而达到模拟无滑移、无穿透效果的目的。需要注意的是，对于本书使用的LBM程序架构，用于反弹的概率密度函数无须经过碰撞操作，因此在进行碰撞操作之前，应先存储需要反弹的概率密度函数。在所有节点进行碰撞操作之后，替换相应固体节点上的概率密度函数值。
接下来，分析第23～30行代码。这段代码在整体上是一个双重循环，外面的一层循环变量为i，内部一层的循环变量为j。所以，位于内层循环的表达式“n=n+1”返回的是两层循环的累计循环次数结果。
首先来看内层的循环的“m(n)=obs(j)+(i-1)*lx*ly；”。注意：obs是从lx×ly的数组中获得的元素序号序列，因此obs中的最大值不超过lx×ly。那么，在后面为何要加“(i-1)*lx*ly”呢?其原因是概率密度函数F()是一个三维数组，从第24行代码“F(：，：，i)=rho*t(i)”可以看出，F()实际上包括9个维度为lx×ly的二维数组。如果将每个二维数组当成一个二维平面，那么F()的结构就可以表示成如图5.21所示的形式。
【思考】已知点p在F(：，：，5)中的obs=1088(px=180，py=8)，求其在全局中的值m?
由于这里lx=300，ly=60，因此答案为：m=1 088+(5-1)×300×60=19 088。得到这个m值后，在MATLAB中通过在代码中用“F(m)”就可直接定位到F(180，8，5)指代的位置。因此，采用一个数值，即可定位到F()的三维空间指定坐标。而在“mb(n)=obs(j)+(oppo(i)-1)*lx*ly”中，mb存储了反弹方向的位置，因此，令F(m)=F0(mb)，就可以精准地将F0中存储的反弹方向的概率密度函数赋予F中对应位置，以实现反弹操作(参见第三部分的第56行代码)。这种做法的优点是：计算过程不涉及无须参与反弹的流体节点，因此可以避免冗余计算。
图5.21 概率密度函数数组结构图
第三部分，LBM主循环体：
第三部分是程序的核心部分，这里大部分代码在介绍LBM方法时就已经介绍了，此处不再赘述。本部分程序与5.4.2节的主程序Poiseuille.m比较，有两个主要变化。变化一，边界条件设置。首先，入口被设为速度边界，出口被设为自由出口边界；其次，处理上下壁面的微观边界条件时，直接用反弹边界代替非平衡外推方法，因此所用的代码只有两行(第59、60行)，而在Poiseuille.m中所用的代码有4行(第52～55行)。变化二，增加了反弹条件设置，用于构造固体壁面或固体区域，主要包括三处：①第41、42行，对固体节点的速度进行矫正；②第54行，在碰撞之前备份反弹量；③第56行，执行反弹操作。
第四部分，计算结果可视化：
其中，第64～66行显示速率场，第67～70行显示涡量场。若需显示涡量图，则取消语句前的注释符并注释屏蔽第64～66行代码即可。
对比5.4节和5.5节的程序内容，可以发现大部分代码重复，为了获得所需的通道结构，只需对基础的通道流程序进行局部修改。至于是否改得正确，则通常需要进行定量验证，只有定量验证通过了，才可确认代码的修改是正确的。比如，在5.4节中，将通道流速度的数值解与解析解进行对比(图5.11)，就是一种定量化的验证方法。