0 引　言

IB-LBM方法，即浸没边界法^[1]（Immersed Boundary Method）与格子Boltzmann方法^[2]（Lattice Boltzmann Method）耦合方法，最早由Feng等^{[3 − 4]}结合，是一种用于模拟流体流动和流固耦合现象的数值方法。它结合了IBM方法在处理复杂边界条件方面的优势和LBM方法在流体模拟上的高效性，能够有效地模拟具有复杂几何形状和边界条件的流体流动问题。另外，由于浸没边界法和格子Boltzmann方法均使用直角网格，在计算运动物体时无需进行网格重构，具有较高的运算效率。

IB-LBM方法的核心思想是通过在计算域中引入一个浸没边界来模拟固体边界，在传统的IB-LBM方法中，浸没边界通过在流体演化方程中增加一个边界作用力项来实现。Shu等^[5]提出了速度修正概念，通过边界速度对其附近流体网格点速度进行修正，由修正前后的速度差值计算得到边界作用力，从而确定浸没边界对流场的影响。相较于传统的IB-LBM，该方法中的速度修正过程使边界上的速度等于同一位置流场中的速度，严格保证了无滑移的边界条件。受到该思想的启发，Wu等^[6]提出了一种隐式修正的狄拉克插值修正方法，王星^[7]对其进行了改进，采用拉格朗日多项式插值进行速度修正，相较于狄拉克插值具有更高的精度。王富海^[8]将多松弛时间模型（MRT）^{[9 − 10]}推广到速度修正算法中，多松弛模型中流体的体黏性和剪切黏性由不同的松弛时间决定，相对于单松弛模型具有更高的灵活性和计算稳定性。Gsell等^[11]通过引入半解析误差修正机制对边界作用力进行校准，有效抑制了流固界面的滑移与渗透现象。在数值方法创新方面，Qin等^[12]将水平集函数与IB-LBM框架相融合，通过在浸没界面施加法向分量跳跃条件，数值实验证实该方法在衰减流场模拟中展现出二阶空间精度特征。Yang等^[13]学者从近场动力学角度出发，构建了结合衰减核函数与表面效应修正的新型计算模型，该模型显著提升了IB-LBM在复杂流动中的数值稳定性。此外，在非牛顿流体领域，基于Giesekus-Oldroyd本构关系的新型IB-LBM框架^[14]被成功开发，为固液两相流中悬浮颗粒的动态特性研究提供了有效的数值模拟工具。

在拉格朗日速度修正IB-MRT-LBM方法上进行改进，针对浸没边界拉格朗日点密集的情况提出了改进的速度修正IB-MRT-LBM，并对计算域进行了局部网格细化以提升计算效率。翼型作为螺旋桨叶片、船舵等船舶部件的本质，在船舶领域得到了广泛应用，具有较高的研究价值。以尾部拉格朗日点密集的NACA0012翼型为研究对象，分别研究了其在不同攻角下的静态绕流性能和不同振荡频率下的推进性能，从静态和动态两个角度验证了算法的有效性，为船舶领域中具有翼型特征部件的流场仿真提供了一种有效的方法。

1 数值方法和改进算法 1.1 多松弛时间格子Boltzmann模型

在LBM中常见的碰撞算子模型包括单松弛时间模型（Single Relaxation Time，SRT）和多松弛时间（Multiple Relaxation Time，MRT）模型，其中SRT模型使用单一的松弛时间来控制平衡态的松弛过程，具有实现简单、计算效率高的特点，但数值稳定性较差；而MRT模型采用多个松弛时间参数控制松弛过程，粒子分布函数$f$向平衡态分布函数${f^{eq}}$弛豫的过程不再由单一的松弛时间控制，而是由碰撞矩阵$\boldsymbol{\Lambda} $控制，计算更为稳定^{[9 − 10]}，能够精确模拟复杂流动。MRT模型的演化方程可以表示为：

$ {f_i}(x + {c_i}\delta t,t + \delta t) - {f_i}(x,t) = - {\boldsymbol{\Lambda }}[{f_i}(x,t) - f_i^{eq}(x,t)]。$

(1)

本文中采用多松弛时间下的D2Q9（二维九个离散速度方向）模型进行二维数值模拟计算，该模型采用正方形格子，如图1所示，每个格点上的粒子有9个离散方向的速度，各方向离散速度由下式给出：

$ {{c_i} = \left\{{\begin{aligned} &{ (0,0),i = 0}，\\ & { \left( {\cos \left[ {\displaystyle\frac{{(i - 1){\text{π}} }}{2}} \right],\sin \left[ {\displaystyle\dfrac{{(i - 1){\text{π}} }}{2}} \right]} \right),i = 1,2,3,4}，\\ &{\sqrt 2 \left( {\cos \left[ {\displaystyle\dfrac{{(i - 1){\text{π}} }}{2} + \frac{{\text{π}} }{4}} \right],\sin \left[ {\displaystyle\dfrac{{(i - 1){\text{π}} }}{2} + \frac{{\text{π}} }{4}} \right]} \right),i = 5,6,7,8} 。\end{aligned}} \right.} $

(2)

在D2Q9模型中，格子间声速${c_s} = 1/\sqrt 3 $，每个离散速度方向的权重系数$w$通过下式计算：

$ {w_i} = \left\{ {\begin{aligned} &{\displaystyle\frac{4}{9},i = 0}，\\ &{\displaystyle\frac{1}{9},i = 1,2,3,4} ，\\ &{\displaystyle\frac{1}{{36}},i = 5,6,7,8}。\end{aligned}} \right. $

(3)

MRT模型的碰撞过程在矩空间内进行，根据变换矩阵$M$，可以确定分布函数从速度空间到矩空间的相互转换关系：

$ m(x,t)=\boldsymbol{M}f(x,t)，$

(4)

$ { \boldsymbol M = \left( {\begin{array}{*{20}{c}} 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} \end{array}}\right)。} $

(5)

因此，式(1)可以表示为：

$ {{f_i}(x + {c_i}\vartriangle t,t + \vartriangle t) - {f_i}(x,t) = - {\boldsymbol M^{ - 1}}S[m(x,t) - {m^{eq}}(x,t)]。} $

(6)

在D2Q9-MRT中，$ \boldsymbol{S} $为松弛时间矩阵；$m$为矩空间；${m^{eq}}$为矩空间的平衡态：

$ S={\boldsymbol{M}}\Lambda{\boldsymbol{M}}^{-1}={\rm{diag}}\left(0,s_e,s_{\varepsilon},0,s_q,0,s_q,s_{\nu},s_{\nu}\right)，$

(7)

$ m=\left(\begin{array}{*{20}{c}}\rho,e,\varepsilon,j_x,q_x,j_y,q_y,p_{xx},p_{xy}\end{array}\right)\mathrm{^T}，$

(8)

$ \begin{split}{m^{(eq)}} = &\rho \Big( 1, - 2 + 3({u^2} + {v^2}),\alpha + \beta ({u^2} + {v^2}),\\ &u, - u,v, - v,{u^2} - {v^2},uv\Big)^{\rm{T}}。\end{split} $

(9)

式(7)中，${s_v}$与SRT模型中松弛频率的取值相同，${s_e}、{s_\varepsilon }、{s_q}$选取比1稍大的数即可。式(8)中，$\rho $为流体密度；$e$为内能；$\varepsilon $与内能$e$的平方有关；${j_x}、{j_y}$分别为动量通量在水平方向与垂直方向的分量；${q_y}$为内能通量在垂直方向分量；${p_{xx}}、{p_{xy}}$分别为应力张量的对角线和非对角线分量。式(9)中，$u、v$分别为流体粒子在水平方向和垂直方向的速度分量；$\alpha 、\beta $为自由参数。

剪切黏性$\nu $和体黏性系数$\varsigma $分别为：

$ \nu = c_s^2\left( {\frac{1}{{{s_\nu }}} - \frac{1}{2}} \right)\delta t，$

(10)

$ \varsigma = c_s^2\left( {\frac{1}{{{s_e}}} - \frac{1}{2}} \right)\delta t。$

(11)

1.2 速度修正IB-LBM的改进

速度修正IB-LBM方法的核心思想是通过建立插值公式得到欧拉网格点和拉格朗日边界点之间的速度关系，从而进行边界作用力的计算，本文主要研究拉格朗日插值速度修正^{[7 − 8]}的IB-LBM算法。

如图2(a)所示，P点表示拉格朗日边界点，A、A₁、B₁为欧拉网格点，假设每个拉格朗日点只影响附近一个格点内的速度，对于垂直网格上拉格朗日点的附近格点A，需要对其进行速度修正。本文中采用拉格朗日插值进行速度修正，修正后的速度用$ u_A^c $来表示：

$ \begin{split}u_A^c = &{u_P}\frac{{({x_A} - {x_{A_1}})({x_A} - {x_{B_1}})}}{{({x_P} - {x_{A_1}})({x_P} - {x_{B_1}})}} + {u_{A_1}}\frac{{({x_A} - {x_P})({x_A} - {x_{B_1}})}}{{({x_{A_1}} - {x_p})({x_{A_1}} - {x_{B_1}})}} +\\& {u_{B_1}}\frac{{({x_A} - {x_{A_1}})({x_A} - {x_P})}}{{({x_{B_1}} - {x_{A_1}})({x_{B_1}} - {x_P})}}。\\[-1pt] \end{split}$

(12)

式中：$u$为水平方向速度；$x$为格点纵坐标。对于修正了2次的修正点，取2次修正的平均值。

而对于拉格朗日点较为密集（内部网格节点数小于3）的边界，如图2(b)所示，其中P₁、P₂为拉格朗日点，A₁、B₁为欧拉网格点，采用上述速度修正算法在处理图2(b)中的密集点浸没边界问题时由于浸没体内部网格点过少会发生发散现象。对于修正点P，本文根据A₁、P₁、P₂、B₁建立拉格朗日插值公式，对其进行三阶精度的速度修正：

$ \begin{split} u_A^c =& {u_{P_1}}\displaystyle\frac{{({y_A} - {y_{A_1}})({y_A} - {y_{P_2}})({y_A} - {y_{B_2}})}}{{({y_{P_1}} - {y_{A_1}})({y_{P_1}} - {y_{P_2}})({y_{P_1}} - {y_{B_2}})}} +\\ &{u_{A_1}}\displaystyle\frac{{({y_A} - {y_{P_1}})({y_A} - {y_{P_2}})({y_A} - {y_{B_2}})}}{{({y_{A_1}} - {y_{P_1}})({y_{A_1}} - {y_{P_2}})({y_{A_1}} - {y_{B_2}})}} +\\ & {u_{P_2}}\displaystyle\frac{{({y_A} - {y_{A_1}})({y_A} - {y_{A_2}})({y_A} - {y_{B_2}})}}{{({y_{P_2}} - {y_{A_1}})({y_{P_2}} - {y_{P_1}})({y_{P_2}} - {y_{B_2}})}} +\\ &{u_{B_2}}\displaystyle\frac{{({y_A} - {x_{A_1}})({y_A} - {y_{P_2}})({y_A} - {y_{P_1}})}}{{({y_{B_2}} - {y_{A_1}})({y_{B_2}} - {y_{P_2}})({y_{B_2}} - {y_{P_1}})}} 。\end{split} $

(13)

图3给出了图2所示情况中2种修正方式的效果对比图，其中固体边界为静止边界$({u_{P_1}} = {u_{P_2}} = 0)$，流场速度为0.1，横轴$y$为各点垂直坐标位置，纵轴$u$为速度大小。可以看出，相较于二阶精度修正，三阶精度修正同时考虑了上下2个边界点对A点的作用，这样得出的结果无疑更合理，因此该修正方式可确保计算稳定性，而对外部流体格点统一使用二阶精度速度修正，从而保证外部格点修正方式的统一性，防止出现因修正方式差异而导致的流线穿透。

计算出边界点附近速度后，边界作用力计算式为：

$ {{f}} = 2\rho \frac{{({u^c} - {u^ * })}}{{\delta t}}。$

(14)

式中：${u^c}$为目标节点修正后速度；${u^*}$为修正前速度。对于多松弛时间的D2Q9离散速度模型，F为离散的边界作用力，可通过式(15)求得，式中$ {f}_x、{f}_y $分别为边界作用力在水平方向和垂直方向的分力。

$ \left\{\begin{split} & F_0=0, \\ & F_1=6\left(1-\frac{s_e}{2}\right)(u{f}_x+v{f}_y), \\ & F_2=-6\left(1-\frac{s_{\varepsilon}}{2}\right)(u{f}_x+v{f}_y)\cdot{f}, \\ & F_3={f}_x, \\ & F_4=-\left(1-\frac{s_q}{2}\right){f}_x, \\ & F_5={f}_y, \\ & F_6=-\left(1-\frac{s_q}{2}\right){f}_y, \\ & F_7=2\left(1-\frac{s_v}{2}\right)\left(u{f}_x-v{f}_y\right), \\ & F_8=\left(1-\frac{s_v}{2}\right)\left(u{f}_y-v{f}_x\right)。\end{split}\right. $

(15)

1.3 改进算法的优势

选取翼长为$D$的NACA0012翼型为研究对象，对均匀绕流中翼型的阻力系数进行计算，雷诺数$Re = 500$，计算区域大小为$10D \times 5D$，采用网格步长为1的均匀网格，边界条件均采用非平衡外推格式。根据不同的取点间隔$ （{\delta }_{y}^{l}=1、2） $设定2种取点方式。2种取点方式如图4所示，其中虚线表示翼型尾部边界，虚线上黑色点表示拉格朗日点。

由于$\delta _y^l = 1$的拉格朗日点在翼型尾端分布较为密集，为了测试2种算法的稳定性，分别使用Matlab软件对文献[8]算法与改进算法进行了编程，首先在$\delta _y^l = 1$时分别使用2种算法进行计算，并对迭代过程中边界内部的网格节点a、b、c的速度变化进行监测。监测结果如图5所示，可以看出，在$t = 32$时文献[8]算法的监测点速度值发生大幅度发散，而改进算法则未发生发散，计算可以稳定进行，说明改进算法加强了处理密集点边界时的稳定性。接下来对$\delta _y^l = 1$情况使用改进算法进行计算，由于$\delta _y^l = 1$时使用文献[8]算法会发散，因此将改进算法计算结果与$\delta _y^l = 2$情况下的文献[8]计算结果进行比较。

计算结果如表1所示，文献[8]算法的阻力系数计算误差达到了11.14%，同时发生多处流线穿透现象（见图6(a)圆圈处），这是由于拉格朗日点数量较少，难以严格保证无滑移的边界条件。而改进算法通过实现密集拉格朗日点的稳定计算，更好地还原出了浸没边界的几何特征，计算结果更为精确，速度边界的无滑移条件得到了保证（见图6(b)）。

2 MRT-LBM二重网格加密

对流体演化剧烈的区域采用二重网格加密技术^{[16 − 18]}以提高计算效率与计算精度。二重网格加密结构如图7所示，加密网格边界为ABC，粗糙网格边界为DEF，设置$n$为2种网格尺寸$\delta x$之比，则2种网格演化的时间步长$\delta t$之比也为$n$，本文中用上标$c$表示粗网格中的参数，上标$f$表示细网格中的参数，则有:

$ n = \frac{{\delta {x^c}}}{{\delta {x^f}}} = \frac{{\delta {t^c}}}{{\delta {t^f}}}。$

(16)

粗网格区和细网格区的交界处，须满足流体力学参数$\delta \rho 、u、v$连续，因此有$\delta _\rho ^c = \delta _\rho ^f，j_x^c = j_x^f，j_y^c = j_y^f$，并且两区域内各矩的平衡态须一致，则有${m^{eq,c}} = {m^{eq,f}}$，所以仅需要确定粗网格区，和细网格区之间不守恒矩的非平衡态部分的关系。

不同区域的剪切黏性和体黏性须保持一致，可推导出加密区域松弛参数和粗网格区域松弛参数之间的关系为：

$ \frac{1}{{s_e^f}} = \frac{1}{2} + n\left( {\frac{1}{{s_e^c}} - \frac{1}{2}} \right) ，$

(17)

$ \frac{1}{{s_\nu ^f}} = \frac{1}{2} + n\left( {\frac{1}{{s_\nu ^c}} - \frac{1}{2}} \right)。$

(18)

松弛参数${s_\varepsilon }$和${s_q}$可根据需要自行确定，本文为方便计算，将粗细网格中的松弛参数关系也设置为:

$ \frac{1}{{s_\varepsilon ^f}} = \frac{1}{2} + n\left( {\frac{1}{{s_\varepsilon ^c}} - \frac{1}{2}} \right),\frac{1}{{s_q^f}} = \frac{1}{2} + n\left( {\frac{1}{{s_q^c}} - \frac{1}{2}} \right)。$

(19)

粗细网格之间的数据传递方程如下：

$ {f^c} = {M^{ - 1}}{m^{eq,c}} + {M^{ - 1}}{\tau ^f}\left( {{m^f} - {m^{eq,f}}} \right)。$

(20)

式中：${\tau ^f} = {\rm{diag}}\left( {1,\displaystyle\frac{{ns_e^f}}{{s_e^c}},\displaystyle\frac{{ns_\varepsilon ^f}}{{s_\varepsilon ^c}},1,\displaystyle\frac{{ns_q^f}}{{s_q^c}},1,\displaystyle\frac{{ns_q^f}}{{s_q^c}},\displaystyle\frac{{ns_\nu ^f}}{{s_\nu ^c}},\displaystyle\frac{{ns_\nu ^f}}{{s_\nu ^c}}} \right)$。通过上式可实现粗细网格直接的数据交换，需要注意的是，粗网格向细网格进行数据传递时，边界上会缺少所对应的点，因此需要对粗网格进行二阶或高阶精度的空间插值来获得相应点的数据，并且由于2种网格时间尺度不同，在$\delta {t^c} + m\delta {t^f}(0 < m < n)$时刻进行计算时还需对其进行时间插值。

3 验证和模型计算

为了验证加密网格及改进速度修正算法的有效性，分别对均匀来流下的圆柱绕流、翼型绕流、翼型拍动进行计算验证，本文计算结果中阻力系数$C_d = \dfrac{{{F_x}}}{{\dfrac{1}{2}\rho U_\infty ^2D}}$，升力系数$C_l = \dfrac{{F_y}}{{\dfrac{1}{2}\rho U_\infty ^2D}}$，其中F_x、F_y分别为浸没体受到的水平方向和垂直方向作用力；$D$为特征长度；$\rho $为流体密度；${U_\infty }$为来流速度。值得一提的是，由于LBM算法的介观特性适用于低雷诺数下粘性效应显著的情况，在低雷诺数流场中表现出良好的数值稳定性，而高雷诺数流场运动粘性较低容易导致LBM程序计算发散，因此本文中的仿真均在低雷诺数流场中进行。

3.1 加密网格数值验证

对$Re = 40$流场中的圆柱绕流进行了仿真，流场左侧均匀来流，流速${U_\infty } = 0.1$，圆柱周围方形区域为加密区域，流场网格与加密区域网格尺寸比为2。分别设置圆面直径$ D=10、20、40 $，划分3种网格方案，其中方案1圆面直径与网格数量最小，网格较粗，方案3圆面直径与网格数量最大，网格较细。计算域尺寸如图8所示。在$Re = 40$时，流场为定常流动，当流场达到稳定状态时，在圆柱后方就会出现一对固定不动的对称涡（见图9）。定义从圆柱尾缘点到涡的尾端的距离与圆柱直径的比值为无量纲化循环区域长度$L$。

计算结果如表2所示，其中$C_d$为阻力系数；$L$为无量纲化长度。3种网格方案阻力系数的计算结果均与文献[7]结果吻合较好，中等网格（方案2）、细网格（方案3）2种网格方案的无量纲循环长度与文献结果吻合较好，证明了该网格加密方法的准确性，而粗网格方案下的无量纲循环长度误差较大，表明该网格加密方法需要足够的网格数量来保证计算准确性。

3.2 翼型绕流计算

在$Re = 1\;000$的情况下，对不同攻角的NACA0012翼型在流场中的受力情况进行了计算，流场采用二重网格加密，粗网格尺寸与细网格尺寸之比$n = 2$，粗网格区域的网格数为$1\;000 \times 500$，翼型周围为细网格区域，网格数为$500 \times 200$，翼型长度$D = 200$，计算攻角为$ {0}^{\circ }、{3}^{\circ }、{6}^{\circ }、{8}^{\circ }、{10}^{\circ }、{12}^{\circ }、{15}^{\circ }、{20}^{\circ } $。计算域尺寸如图10所示，其中流场左边界为速度入口边界，设定入口速度${U_\infty } = 0.1$；流场右边界为压力出口边界，设定出口密度$\rho = 1$；流场的上下边界均为无滑移的壁面。以上边界条件均采用非平衡外推模式进行计算。

图11给出了5个计算攻角下的流场速度分布图和流线图。可以看出，小攻角$ （AOA < 8^\circ） $时，绕翼型的流动是层流，流场中几乎不出现旋涡；而在大攻角$ （AOA\geqslant 8^\circ） $时，在翼型的背部会明显出现旋涡，旋涡的大小与攻角的大小有关，流体的流动也变为非定常流动，期间分离涡会周期性地从翼型背部脱落、传播。流线图中未发生流线穿过流体的现象，说明该算法可以保证速度边界的无滑移条件。传统的拉格朗日速度修正在翼型尾端计算时会发生发散，因此翼型计算采用改进算法，计算过程中未发生发散现象，证明了改进算法在进行静态浸没边界计算时的稳定性。

图12为阻力系数$C_d$和升力系数$C_l$随攻角变化的情况。可以看出，翼型的升力系数和阻力系数随着翼型攻角的增加而增大，变化趋势与文献[19]结果一致，在小攻角的翼型绕流中，计算结果较为稳定，与文献结果吻合也较好；而在大攻角的翼型绕流中，升阻力随着旋涡的脱落而发生周期性变化。计算结果与文献的误差主要由网格方案差异导致，总体而言相差较小，说明计算结果较为精确。

3.3 拍动翼型数值模拟

对不同拍动频率的NACA0012翼型拍动进行了数值模拟，拍动翼型翼长为$D$，计算区域为$12D \times 6D$，网格数为$600 \times 300$，加密区域为$2.0D \times 2.4D$，网格数为$200 \times 240$，粗细网格尺寸比为2。拍动翼型的运动包括升沉运动与俯仰运动，表达式为：

$ h(t) = {h_m}\cos(2{\text{π}} \omega t) ，$

(21)

$ \theta (t) = {\theta _m}\sin(2{\text{π}} \omega t)。$

(22)

式中：$h(t)$为$t$时刻翼型旋转点垂向位置；${\theta _m}$为翼型的升降幅值；$\theta (t)$为$t$时刻翼型旋转角；${\theta _m}$为翼型的俯仰角幅度；$\omega $为振动频率，与无量纲频率$S_t$的转化关系如式(23)所示。表3给出了本文拍动翼型仿真中各参数的符号与取值。

$ S_t = \frac{{\omega D}}{{{U_{_\infty }}}}。$

(23)

翼型拍动所产生的推力用平均推力系数$Ct$为：

$ C_t=\displaystyle\frac{-\displaystyle\frac{1}{T}\displaystyle\int_0^TF_x(t)\, \mathrm{d}t}{0.5\rho DU_{\infty}^2}。$

式中：$T$为翼型拍动周期；${F_x}$为$t$时刻翼型水平方向的瞬时力。

图13给出了$S_t = 0.3$时翼型在一个拍动周期内的流场速度分布图，图14(a)表示一个拍动周期内的推力系数变化。可以看出翼型拍动的前1/4个周期内行程向下并发生逆时针旋转，上表面出现涡脱离，上下表面的压力差增大导致推力逐渐增大；接下来的1/4个周期内翼型继续向下运动，发生顺时针旋转直至水平，上表面脱离涡逐渐消散，翼型上下表面压力差减小，推力随之降低。向上冲程的半个周期内，涡从翼型下表面产生、脱离直至消散，推力系数也随之发生变化。图14(b)给出了不同$S_t$下一个运动周期的平均推力系数$ Ct_{\mathrm{mean}} $，可以看出在$S_t = 0.2 \sim 0.3$的区间内，平均升力系数随$S_t$的增大而不断增大，这是由于更高的频率对应着更高的拍动速度。涡的发生更加剧烈，所产生的推力也就随之增加。同时，计算结果与文献[20]吻合良好，说明改进算法在进行运动物体的流场模拟时具有较高的精度。

4 结　语

本文针对拉格朗日点分布密集的浸没边界，对速度修正IB-MRT-LBM提出了改进，为了验证改进算法的计算稳定性，采用二重网格加密技术对流场网格进行局部细化以提升精确度和计算效率，通过对NACA0012翼型进行流场仿真，得到以下结论：

1）提出一种适用于拉格朗日点密集的浸没体模型的速度修正IB-MRT-LBM，通过在IB-MRT-LBM中引入三阶拉格朗日插值方法，避免了传统二阶插值处理密集边界点时的发散问题，增强了算法的稳定性，拓展了算法的使用范围。

2）使用改进算法进行了$ Re=1\;000 $流场中不同攻角翼型的升阻力系数，计算结果表明，小攻角$（AOA < 8^\circ） $时绕翼型流动为定常流动，大攻角$（AOA \geqslant 8^\circ） $时绕翼型流动为非定常流动，升阻力系数随着攻角的增加而逐渐变大，改进的速度修正的IB-MRT-LBM表现出良好的稳定性和计算精度，能够较好地处理静态浸没边界问题。

3）使用改进算法进行了$Re = 500$流场中6种振动频率下的拍动翼型仿真，拍动翼型的平均推力系数与随着频率的增大而增大，计算结果较为精确，证明改进的速度修正IB-MRT-LBM适用于运动边界的流场模拟。