格子玻尔兹曼方法（LBM）中的MRT作用力模型
格子玻尔兹曼方法（LBM）MRT作用力模型
格子玻尔兹曼方法搞流动模拟的老司机都知道，MRT（多松弛时间）模型可比单松弛时间模型（BGK）香多了。这玩意儿最大的特点就是数值稳定性强，边界条件处理灵活，特别适合折腾复杂几何形状的流场。今天咱们就扒一扒MRT里那个让人又爱又恨的作用力模型实现。
先来看个经典场景：二维Poiseuille流动模拟。咱们用D2Q9格子结构，碰撞步的实现代码大概长这样：
def collide_mrt_force(f, force, rho, u, tau_vec): m = np.zeros(9) meq = np.zeros(9) M = np.array([ [1, 1, 1, 1, 1, 1, 1, 1, 1], [-4, -1, -1, -1, -1, 2, 2, 2, 2], [4, -2, -2, -2, -2, 1, 1, 1, 1], [0, 1, 0, -1, 0, 1, -1, -1, 1], [0, -2, 0, 2, 0, 1, -1, -1, 1], [0, 0, 1, 0, -1, 1, 1, -1, -1], [0, 0, -2, 0, 2, 1, 1, -1, -1], [0, 1, -1, 1, -1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 1, -1, 1, -1] ]) # 转换到矩空间 m = np.dot(M, f.T) # 计算平衡态矩 meq[0] = rho meq[1] = rho*(-2 + 3*(u[0]**2 + u[1]**2)) meq[2] = rho*(1 - 3*(u[0]**2 + u[1]**2)) # ...其他矩分量计算... # 添加力项（关键部分！） S = np.diag(1/tau_vec) force_term = np.dot(np.linalg.inv(M), (np.eye(9) - 0.5*S) @ np.dot(M, force)) # 碰撞过程 m_post = m - S @ (m - meq) + force_term # 转换回速度空间 return np.dot(np.linalg.inv(M), m_post)这段代码里藏着几个魔鬼细节。首先是转换矩阵M的设计，这玩意儿把分布函数从速度空间转到矩空间，不同矩对应不同的物理量。比如第一个矩就是质量守恒，第三个矩对应能量项。作用力模型的关键在force_term的计算，这里用到了Guo提出的修正项，比传统方法能更好处理非均质力场。
特别注意这个(np.eye(9) - 0.5*S)的操作，这可不是随便写的玄学系数。它确保了力项在离散化过程中的二阶精度，少了这个0.5因子的话，模拟结果会出现明显的数值扩散。曾经有个哥们儿把这个系数写成了1.0，结果模拟的涡街直接变成抽象画，debug了三天三夜才找到这个坑。
格子玻尔兹曼方法（LBM）MRT作用力模型
再说说松弛时间tau的选择。MRT的精髓就在于不同的矩有不同的松弛率，比如对于剪切相关的矩（对应黏性系数），tau_v通常取1.0附近；而高阶矩的tau可以调大些，这样数值稳定性更好。调参的时候记得遵守这个原则：奇数阶矩的松弛时间要大于0.5，否则模拟直接崩给你看。
实际跑代码时，处理边界条件的姿势也很重要。比如用非平衡外推法处理曲面边界时，作用力项需要跟着分布函数一起做外推。这里有个骚操作：可以把力项先转到矩空间再做外推，收敛速度比传统方法快两倍不止。具体实现就是在边界节点处：
# 边界节点处理示例 boundary_node = get_boundary_nodes() m_boundary = M @ f[boundary_node] force_m = M @ force[boundary_node] # 非平衡外推 m_neighbor = M @ f[neighbor_node] force_neighbor = M @ force[neighbor_node] m_boundary_post = 2*m_neighbor - m_boundary + (force_m + force_neighbor)*0.5*dt这套组合拳打下来，复杂几何边界的流场也能稳如老狗。不过要注意内存开销，MRT比BGK多用了近一倍的存储空间，特别是三维模拟的时候，显存分分钟爆炸。这时候就得祭出内存优化大法——把不用的矩分量及时清零，或者用稀疏矩阵存储转换矩阵。
最后说个实战技巧：用MRT模拟非牛顿流体时，可以通过动态调整松弛时间来反映黏度变化。比如卡森流体模拟中，这样改松弛参数：
tau_v = 0.5 + (mu_eff / c_s**2) / dt tau_vec[4] = tau_v # 调整剪切相关松弛时间 tau_vec[6] = 1.2 # 高阶矩保持固定这种动态调整策略能让MRT模型适应各种复杂的本构方程，比传统LBM方法灵活得多。不过要注意时间步长dt和松弛时间的耦合关系，调不好的话整个流场会像果冻一样抖个不停。
总之，MRT作用力模型就像瑞士军刀，功能强大但需要精细操作。摸透这些门道后，处理复杂流动问题就跟玩似的，管它什么微通道流动还是血液细胞运输，统统拿下不费劲。