LBM格子玻尔兹曼方法模拟圆柱扰流及升阻力系数对比
格子玻尔兹曼方法（LBM）在流场模拟中有独特的魅力，特别是处理复杂边界时，圆柱绕流就是个典
型场景。今天咱们用Python实现一个2D圆柱扰流模拟，重点看看怎么抓取升力系数和阻力系数。
先整点核心代码。D2Q9模型是基础骨架，碰撞项用BGK近似搞定：
```python
def collide(f, omega):
rho = np.sum(f, axis=2)
u = np.dot(f, e).transpose(2,0,1) / rho
feq = equilibrium(rho, u)
f -= omega * (f - feq)
return rho, u
def equilibrium(rho, u):
cu = np.dot(e, u.transpose(1,0,2))
usq = np.sum(u**2, axis=0)
feq = w[:,None,None] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*usq)
return feq.transpose(1,2,0)
```
这里用了矢量化计算提升效率，注意速度投影cu的处理方式——把三个维度计算压缩成矩阵运算，比
逐点循环快至少20倍。
边界处理是绕流模拟的关键。圆柱用反弹格式处理，这里有个小技巧：用坐标掩码快速定位边界点
```python
def bounce_back(f, mask):
f[[0,1,2,3,4,5,6,7,8], mask] = f[[0,5,6,7,8,1,2,3,4], mask]
```
mask是提前生成的圆柱位置布尔矩阵，这种索引操作比逐点判断节省90%时间。
计算受力时，动量交换法比应力积分更稳定：
```python
Fx = np.sum( (f[[1,5,8], mask] - f[[3,6,7], mask]) * e[1,[1,5,8]] )
Fy = np.sum( (f[[2,5,6], mask] - f[[4,7,8], mask]) * e[2,[2,5,6]] )
