王文全，王金霖，骆佳玲

(1.昆明理工大学 建筑工程学院，云南 昆明 650500；2.四川大学 水力学与山区河流开发保护国家重点实验室，四川 成都 610065；3.四川大学 水利水电学院，四川 成都 610065)

自1972年PESKIN[1]提出浸入边界方法(immersion boundary method,IBM)以来，该方法就广泛用于模拟生物力学、流固耦合和多相流动等问题。IBM的基本思想是将浸入流体中的固体视为边界，并且通过边界上的恢复力来表示流体与固体之间的相互作用。根据恢复力的不同计算方法又形成了不同的浸入边界方法。如根据施力点的范围可以分为扩散界面法和界面一致法。在扩散界面法中，通过使用光滑函数将奇异力分布到周围的背景网格节点上，从而消除浸入边界。扩散界面法可以进一步分类为经典IBM[2-3]，直接力法[4-6]和惩罚力法[7-8]，具体取决于浸没固体影响流体的虚拟力的计算方式。传统扩散界面法的一个固有特征是，由于使用了离散函数或光滑函数，浸入边界不是清晰的，而是横跨多个网格节点的模糊边界。消除这种缺陷的界面一致方法被提出[9]，但是传统界面一致浸入边界法绝大多数都是显式求解方法，对于动边界问题，插值运算十分复杂和耗时，同时求解添加了动量源项的Navier-Stokes(N-S)方程时，数值解不易收敛甚至会发散，而且压力求解效率较低，而格子Boltzmann方法(lattice boltzmann method,LBM)是一种可快速求解流场和压力场的数值方法，提高了计算效率。浸入边界-格子Boltzmann方法(immersion boundary-lattice boltzmann method,IB-LBM)在求解流固耦合(fluid-structure interaction,FSI)问题中取得了可观的进展[10-11]，例如，侧重于拍打柔性箔的仿生模拟[12-13]和多个弹性结构的拍打[14]，颗粒-流体相互作用流[15]，固液相变问题[16]。为此，本文有机结合IBM与LBM的各自优点，实现了一种流体固体相互作用的计算方法，采用隐式扩散界面法来处理流固界面间的相互作用，界面力通过严格满足无滑移边界条件而获得，很好地消除界面流线穿透和边界模糊等缺陷。

1 数学模型与数值方法

1.1 数学模型

浸入边界法的黏性不可压缩流的控制方程可写为，

(1)

▽·u=0

(2)

式中：u是流体速度；p是流体压力；雷诺数定义为Re=U∞L/μ，μ是动力黏度，U∞是自由流速度，L是固体的特征长度，ρ是流体密度。浸没固体施加给流体的力可表示，

(3)

式中：x和X分别代表流场的网格点坐标(整个计算域)和IB点(固体边界的拉格朗日坐标)；F(X(s),t)是固体边界力密度；δ(x-X(s,t))是格子网格和拉格朗日网格之间的信息传递函数，其具体形式将在下一部分中进行讨论。为了求解方程式(1)和(2)，基于浸入边界的格子Boltzmann方程可写为，

(4)

(5)

(6)

ρ=∑αfα，ρu=∑αfαeα+f△t

(7)

压力p可以根据状态方程来计算，

(8)

有效格子松弛时间τLB与运动黏度νLB有关，如式(9)所示：

(9)

1.2 数值实现

等式(7)可以进一步写为，

u=u′+△u

(10)

式中：u′是求解中间速度，ρu′=∑αfαeα，△u=f△t。格子点处的力密度f的分布函数通过公式(3)从边界力F的分布函数而来，其离散形式可表示为，

(11)

(12)

(13)

(14)

将(11)式代入(10)式得,

(15)

将(15)式代入(14)式可得,

(16)

2 数值算例

2.1 圆柱绕流

(17)

δ(r)=

(18)

由图1(a)可见，高阶光滑分段函数计算效率略低，但对阻力系数的收敛速度没有影响,见图1(b)。使用高阶的光滑函数可以有效消除数值振荡[19]，因此本文采用高阶光滑分段函数进行流固间的信息交换。如表1所示，Re=40时阻力系数和涡旋长度与文献相吻合。图2为不同雷诺数下的流线图。当雷诺数小于40时，在圆柱体的尾部产生两对稳定的对称涡旋，并且随着雷诺数的增加，旋涡长度会增加。当Re=80和200时，涡交替脱落破坏了旋涡的对称性。流线没有穿透圆柱体表面，表明无滑移边界条件得到精确的满足。

表1 圆柱绕流的阻力系数和尾涡长度的比(Re=40)

图1 不同光滑函数下时间和阻力系数演化(Re=40)

图2 不同雷诺数下的流线分布图

2.2 NACA0012机翼绕流

翼型边界由160个节点离散，雷诺数Re=500。图3给出了NACA0012机翼的横截面，图4显示了x方向的速度轮廓以及翼型周围的流线，从图中可以看出无滑移边界条件得到很好的满足。图5(a)和图5(b)分别表示在x=-0.5, -0.25, 0.0, 0.25和0.5断面x方向速度和y方向速度分布。从图中可以发现，机翼表面附近的速度接近于零。

图3 NACA0012机翼的横截面

图4 NACA0012机翼绕流x方向速度分布和流线图

图5 NACA0012在不同位置处不同方向的速度

定义压力系数为，

(19)

式中：p∞是进口压力；pB是沿机翼表面的压力，可以从以下公式获得，

(20)

沿翼型表面的压力系数分布如图6所示，本文结果与IMAMURA等[23]和WU等[11]的数值结果吻合较好。

图6 沿NACA0012机翼的压力系数

2.3 转子的被动运动

图7 转子示意图

(21)

式中：Fi表示作用在坐标为(xi,yi)的第i个浸入边界点上的外力；fxi和fyi分别表示x方向和y方向上的外力；di表示第i个浸入边界点和旋转点之间的距离；m代表离散边界点的数量，△s是两个相邻离散边界点之间的弧长。

根据角动量定律，外部转矩可以表示为,

(22)

式中：LO和JO分别表示对旋转点O的角动量和旋转惯量；ω为角速度；α为角加速度；转子的材料密度取为7.85ρf，ρf为流体密度。阻力系数和升力系数定义为，

(23)

图8为不同Re下转矩、角速度和角位移随时间的变化。当Re=3和70时，旋转不稳定，当Re=1 000时，可保持恒定的角速度稳定地旋转。图9为不同Re下的阻力系数和升力系数的变化。不同Re下阻力系数保持恒定，但Re=3和70时，升力系数波动失稳，Re=1 000时，升力系数在零值附近呈周期性波动，转子可以保持稳定旋转，因此升力对诱导转子旋转具有重要作用。图10为不同Re下三个不同时刻转子附近的涡结构，Re=3和Re=70时涡脱落不稳定，导致升力大幅波动。Re=1 000时脱落涡强度更大且涡结构更复杂，但稳定的周期性涡脱落有利于转子稳定的被动运动。

图8 不同雷诺数下转矩、角速度和角位移的演变

图9 三个雷诺数的阻力系数和升力系数

图10 三个雷诺数在三个时刻的涡旋结构

3 结论

结合隐式浸入边界和格子Boltzmann方法，实现了流固耦合的模拟。采用格子Boltzmann方法快速获得预测的速度场和压力场，耦合界面力由精确的无滑移边界条件获得，通过高阶光滑函数实现流固间信息交换。该方法具有易于实施，计算效率高，并且减少了非物理振荡的优势。固定圆柱和NACA0012翼型绕流的数值结果与文献相吻合，验证了该方法的有效性。对被动旋转转子进行仿真，在低雷诺数Re=3和70下转子无法保持稳定旋转，当雷诺数Re=1 000时，转子进行稳定的旋转运动。