Page 55 - 《应用声学》2025年第2期

第 44 卷 第 2 期 杨雨等： 微通道内声传播格子 Boltzmann 建模及其特性研究 315

(1) 首先根据具体问题确定物理模型和模拟区 对角化矩阵为
 
域大小，然后均匀地划分网格。因为本文只考虑二   1 0 1
2
维的情况，所以只需要对流场进行二维划分，然后将 c −c s ρ 0  2c 2 s 2c 2 
s
s 
  
 1 1 
流场变成离散的空间。 P x = 0 0 1  , P x −1 =  0  ,


 −

   2ρc s 2ρc s 
(2) 设置模拟参数，当 t = 1 时，将流场的宏观 c 2 c s ρ 0  
量和分布函数初始化并相应地赋初始值，根据粒子 s 0 1 0
演化过程得到粒子的分布函数。  1 1 
  0
2
(3) 在t时刻执行碰撞迁移步骤。 c 0 −c s ρ  2c 2 s 2c 2 s 
s
   
 
(4) 根据具体的实际问题，选择适合的边界条 P y = 0 1 0  , P y −1 =  0 1 0  .


件，然后处理边界节点。  2   


c 0 c s ρ  1 1 
s
(5) 根据具体问题求解节点上宏观密度和速度 − 0
2ρc s 2ρc s
等宏观量，不同模型的宏观量计算有所不同。 (14)
(6) 判断计算收敛与否，若收敛直接输出宏 式(10)与P x 相乘得到
观量及所需结果，若不收敛需要返回第 (3) 步继续 ∂n ∂n
计算。 P x ∂t + Λ x P x ∂x = 0. (15)
上述 LBM 的模拟步骤中，前两步主要为初始 式(10)对x求导可以表示为
化数据，后四步主要为模拟，其中碰撞和迁移、边界  
I x,1
处理最重要，这些步骤将施加在模拟区域的所有的 ∂n ∂n  


X = P x −1 Λ x P x = P x −1 I x , I x =  I x,2  ,
粒子上。 ∂x ∂x  
I x,3
1.2 声学无反射边界
(16)
由于声学仿真常用的边界会导致声波反射，从
其中，I x 表示特征向量，I x 组成部分是
而使压力场受到干扰，所以学者提出了无反射边界 [ ∂ρ ]
克服了这一问题。 I x,1 = (u x − c s ) c 2 s − c s ρ ∂u x ,
∂x ∂ x
首先假设边界上的守恒方程是近似欧拉的，即 ∂u y
不考虑黏性的影响。对于没有外力二维情况，欧拉 I x,2 = u x ∂x ,
[ ]
守恒方程可以用向量和矩阵形式表示为 2 ∂ρ ∂u x (17)
I x,3 = (u x + c s ) c
s ∂x + c s ρ .
∂n ∂n ∂n ∂ x
+ X + Y = 0, (10) 在x边界上，用I 表示特征向量I x ：
′
∂t ∂x ∂y x

其中，流体变量向量是 n = (ρ, u x , u y )，矩阵 X、Y  I x , 输出特性,
I = (18)
′
分别为 x  0, 输入特性.
   
u x ρ 0 u y 0 ρ 在 x 边界处，有几种不同的方法来演化流体变
   
 2   
X = c /ρ u x 0  , Y = 0 u y 0  . (11) 量，可以将式(12)修改为
s
   
2
0 0 u x c /ρ 0 u y ∂n = −P −1 ′ ∂n , (19)
I − χY
s
∂t x x ∂y
矩阵经过对角化后为 X = P x −1 Λ x P x 和 Y = 其中，χ = 3/4 [25] 。
P y −1 Λ y P y ，其中 Λ x 和 Λ y 表示含有特征值 X 和 Y 为了确定时间导数 ∂n/∂t，必须根据式 (18) 和
的对角矩阵，对应地表示为
式(19)估计特征向量，这需要估计空间导数。对于x
Λ x =diag (λ 1 , λ 2 , λ 3 )=diag (u x − c s , u x , u x + c s ) , 边界上的 x 导数，可以通过单边二阶有限差分方法
(12) 来实现：
∓3n i (x) ± 4n i (x ± ∆x) ∓ n i (x ± ∆x)
∂n i
Λ y =diag (µ 1 , µ 2 , µ 3 )=diag (u y − c s , u y , u y + c s ) . (x) ≈ ,
∂x 2∆x
(13)
(20)