钝体绕流在实际工程中有着广泛的背景和工程价值，具有一定粘性的流体经过钝体结构后在钝体的尾流中会产生各种尺度的涡，在脉动载荷作用下，结构升阻力、表面压力等都发生变化，振动和破坏由此产生。而当钝体被强迫以指定规律运动时，尾流场及受力都会产生变化，旋转圆柱正是如此。圆柱的旋转运动能在一定程度上降低激流振荡，修正柱体尾流。Prandtl^[1]进行了旋转圆柱的受力试验研究，结果显示均匀流中圆柱纵横比增加会引起整体升力系数的增加。Nobari等^[2]基于有限元的CBS算法发现超过临界无量纲旋转速度时，涡流的脱落被高度抑制。Lam^[3]基于PIV技术进行试验分析了剪切层产生的大尺度涡的特性。当前对于旋转圆柱的研究主要集中在有限元、有限差分以及实验等方法。

近几年，一种介观方法—格子玻尔兹曼(lattice Boltzmann method, LBM)应运而生，已应用于湍流、多组分流等多个研究领域^[4]。不同于有限体积等方法，LBM基于微观统计力学的思想从介观的角度，引入BGK模型，用碰撞和迁移求解宏观流场，可以解决连续介质假说不成立或难以用方程数值求解的复杂系统问题，在边界处理上采用Bounce-Back格式简化了传统数值求解时的映射和坐标转换过程。但LBM需要使用边界内外两侧的网格逐步迭代表示边界，对于更加复杂的边界的模拟比较困难。沉浸边界法(immersed boundary method, IBM)使用Lagrangian点和Eulerian点来描述边界和处理流场信息，对于外形复杂、多尺度以及大变形边界问题的处理更加简单。IBM等都使用笛卡尔网格，因此将二者结合既能很好地满足物理量的守恒性，又能在处理边界上节省资源^[5]。IB-LBM最初是由Feng等^[6]提出和发展的，他们通过引入IBM取代了以往多排格子的计算方式，但其结果受限于胡克定律中弹性系数的取值。Cheng等^[7]引入非定常不均匀力以及源项，改进了BGK方程的额外项，提出了新的体力计算方法，使得收敛性增加到二阶。Wu等^[8]基于带外力的格子Boltzmann方程，将宏观速度划分为中间速度和修正速度，使得无滑移条件得到了很好的满足。之后Wu等^[9]又将恢复力作为未知量，按照强制无滑移条件的方法对IB-LBM进行了改进，解决了TLLBM在非均匀格子计算中的权重系数计算量太大的问题。Suzuki等^[10]为解决速度梯度不连续问题提出了一种通过逐步迭代使计算精度提高到二阶的高阶格式。IB-LBM因其适合处理动边界及外形复杂的问题得到了迅速的发展和应用。Takada等^[11]采用格子Boltzmann方法研究了雷诺数在200和500时以给定转速突然旋转的二维圆柱的流场和受力，并与已有的试验值和FVM结果进行了对比，证明LBM方法可以适用于移动的曲面边界的问题求解。Fallah等^[12]基于多松弛时间的格子Boltzmann方法研究了雷诺数为100时非牛顿流体中不同旋转速度和幂律指数情况下的二维旋转圆柱的流场和受力，模拟结果说明LBM在模拟曲面边界的非牛顿流体时是可行的。郑海成^[13]基于LBM方法模拟了并列双圆柱旋转绕流。

以往的研究多集中在二维的数值模拟上，但是实际工程中往往涉及三维计算，因此，本文采用IB-LBM进行均匀来流条件下三维旋转圆柱的研究，分析流场结构及升阻力变化。

1 数值计算方法 1.1 格子玻尔兹曼法

LBGK模型是目前为止应用较为广泛的LBE模型，这是一种单松弛时间单相模型，其演化方程为：

$ \begin{array}{l} {f_{\rm{i}}}\left( {\mathit{\boldsymbol{x}} + {\mathit{\boldsymbol{e}}_i}{\delta _t}, t + {\delta _t}} \right) - {f_{\rm{i}}}(x, t) = \\ - \omega \left[ {{f_i}(x, t) - f_i^{(eq)}(x, t)} \right] \end{array} $

(1)

式中：f_i(x, t)为t时刻 x 处的粒子分布函数；δ_t为时间步长；ω=1/τ为松弛参数；e_i为离散速度。对于三维模型的计算，这里选择D3Q19模型，相比较于二维计算，三维问题涉及的离散方向更为复杂，由原来的平面上9方向变为空间的19方向，如图 1，相应的宏观速度组 e_i可以定义为：

$ {\mathit{\boldsymbol{e}}_i} = c \\ \left[ {\begin{array}{*{20}{c}} 0&1&{ - 1}&0&0&0&0&1&{ - 1}&1&{ - 1}&1&{ - 1}&{ - 1}&1&0&0&0&0\\ 0&0&0&1&{ - 1}&0&0&1&{ - 1}&{ - 1}&1&0&0&0&0&1&{ - 1}&1&{ - 1}\\ 0&0&0&0&0&1&{ - 1}&0&0&0&0&1&{ - 1}&1&{ - 1}&1&{ - 1}&{ - 1}&1 \end{array}} \right] $

式中：f_i^(eq)(x, t)为平衡态分布函数，即：

$ f_i^{(eq)}(x, t) = {w_i}\rho \left[ {1 + \frac{{{\mathit{\boldsymbol{e}}_i} \cdot \mathit{\boldsymbol{u}}}}{{c_s^2}} + \frac{{{{\left( {{\mathit{\boldsymbol{e}}_i} \cdot \mathit{\boldsymbol{u}}} \right)}^2}}}{{2c_s^4}} - \frac{{{\mathit{\boldsymbol{u}}^2}}}{{2c_s^2}}} \right] $

(2)

式中：w_i为权重系数；c_s为与声速有关的参数；w₀=1/3, w_1-6=1/18, w_7-18=1/36，$c_{s}=c / \sqrt{3}$；密度ρ是固定的常数，由此得到宏观的压力和速度为$\rho = \sum\limits_\alpha {{f_\alpha }} , \rho \mathit{\boldsymbol{u}} = \sum\limits_\alpha {{e_\alpha }} {f_\alpha }$。这里流体的粘性与松弛时间有关υ=(2τ-1)/6。

1.2 沉浸边界-格子Boltzmann方法

本文基于Suzuki和Inamuro^[10]提出的高阶格式进行IB-LBM程序编写，即对分布函数进行分解，由LBM基本方程得到未修正的分布函数，由力源项引起分布函数的变化得到分布函数的修正值。

如图 2，已知边界的点的位置 x_k^b，假定强迫运动的规律已知，即边界处的速度为 U_k^b，Suzuki和Inamuro方法是采用ghost点计算力源项，即边界上的力源项由边界内部ghost点 x_k^g_m和边界上 x_k^b处的值插值得到。而ghost点的位置和速度是由同一法向方向上的流体点和边界点值插值得到的。流体点 x_f^b可以由边界点位置和法向方向n_k(图 2箭头方向)得到，即 x_f^b= x_k^b+dxn_k，流体点速度可以用迁移碰撞计算得到的未修正速度插值得到，即：

$ \mathit{\boldsymbol{u}}_i^*\left( {\mathit{\boldsymbol{x}}_k^f} \right) = \sum\limits_x {\mathit{\boldsymbol{u}}_i^*} (\mathit{\boldsymbol{x}})W\left( {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{x}}_k^f} \right){({\rm{d}}x)^d} $

(3)

式中：W(x - x_k^f)为权重函数，由于计算的是三维问题，因此计算涉及3个方向的插值，即：

$ W(x, y, z) = \frac{1}{{{\rm{d}}x}}\phi (x/{\rm{d}}x) \cdot \frac{1}{{{\rm{d}}x}}\phi (y/{\rm{d}}x) \cdot \frac{1}{{{\rm{d}}x}}\phi (z/{\rm{d}}x) $

(4)

这样可以得到ghost点的位置和速度，即：

$ \mathit{\boldsymbol{x}}_k^{{g_m}} = x_k^b - m{\rm{d}}x{\mathit{\boldsymbol{n}}_k} $

(5)

$ \mathit{\boldsymbol{u}}_k^{{g_m}} = \mathit{\boldsymbol{u}}_k^b + m\left[ {\mathit{\boldsymbol{u}}_k^b - {\mathit{\boldsymbol{u}}^*}\left( {\mathit{\boldsymbol{X}}_k^f} \right)} \right] $

(6)

进而分别求出ghost点和边界点上的体力，插值得到力源项即：

$ g\left( {\mathit{\boldsymbol{x}}_k^b, t + {\delta _t}} \right) = \frac{1}{{\delta x}}\left[ {u_k^b - {u^*}\left( {x_k^b, t + {\delta _t}} \right)} \right] $

(7)

$ g\left( {\mathit{\boldsymbol{x}}_k^{{g_m}}, t + {\delta _t}} \right) = \frac{1}{{\delta x}}\left[ {u_k^{{g_m}} - {u^*}\left( {\mathit{\boldsymbol{x}}_k^{{g_m}}, t + {\delta _t}} \right)} \right] $

(8)

$ \begin{array}{l} g\left( {x, t + {\delta _t}} \right) = \sum g \left( {\mathit{\boldsymbol{x}}_k^b, t + {\delta _t}} \right)W\left( {x - \mathit{\boldsymbol{x}}_k^b} \right){\delta _{{v_b}}} + \\ \sum g \left( {\mathit{\boldsymbol{x}}_k^{{s_m}}, t + {\delta _t}} \right)W\left( {x - \mathit{\boldsymbol{x}}_k^{{g_m}}} \right){\delta _{{V_{{g_m}}}}} \end{array} $

(9)

式中：$u_k^bu_k^{{g_m}}, {u^*}\left( {\mathit{\boldsymbol{x}}_k^b, t + {\delta _t}} \right), {u^*}\left( {\mathit{\boldsymbol{x}}_k^{{g_m}}, t + {\delta _t}} \right)$分别为实际边界、ghost点，未修正时边界点和未修正ghost点速度。δ_Vb和δ_Vgm为体积元。

整个IB-LBM计算流程为：

1) 初始化，指定u₀、ρ₀以及f_i(x, 0)，计算f_i^(eq)(x, t)，并定义边界的运动；

2) 根据演化方程进行碰撞迁移，得到未修正部分的速度和密度；

3) 采用多次迭代方式计算和修正体力；

4) 利用步骤3)修正的体力进行分布函数修正：

$ \begin{array}{l} {f_{\rm{i}}}\left( {x + {e_i}{\delta _t}, t + {\delta _t}} \right) = f_{\rm{i}}^*\left( {x + {e_i}{\delta _t}, t + {\delta _t}} \right) + \\ 3{\rm{d}}x{w_i}{e_i}g(x, t); \end{array} $

5) 重复步骤2)~4)直到满足收敛条件。

2 计算工况与边界条件

本文对不同相对转速(α=ωD/2U₀, 这里ω为圆柱的旋转角速度，U₀为来流速度)及雷诺数(Re=U₀D/v)条件下旋转圆柱进行了数值模拟。柱体直径D=0.1 m，柱体高H=4D，整个计算域设置为30D×20D×4D，这样流动的三维性得以发展，消除了周围边界的影响^[14]。柱体结构中截面的圆心取在坐标原点O处，圆心距离入口10D，距离出口20D。Inlet为速度边界，outlet为流出边界，其余各面均采用周期性边界，如图 3。流体沿x轴正向流动，速度为U₀，密度为ρ。计算空间步长dx取为0.006，时间步长δ_t为0.002，模拟步数100 000步。

通过改变Re和α，探讨不同参数条件对旋转圆柱流场及受力的影响。具体工况如表 1。

3 计算结果分析 3.1 数值验证

为了验证本文IB-LBM代码的准确性，首先对孤立的三维圆柱进行数值模拟。计算域及模型等设置与第2小节设置一致。表 2给出了本文和文献计算的Strouhal数与力系数，图 4、5分别给出了流线和旋度结果图^[15-16]。

这里定义阻力系数C_d=2F_x/(ρU₀²A)，升力系数C_l=2F_y/(ρU₀²A)，Strouhal数St=Df_v/U₀，其中F_x与F_y分别为x、y方向上的力，A为柱体在yoz面的投影面积，f_v为涡脱落频率。从表中可以看到本文的计算结果与文献中的结果较为吻合，因此用本文的IB-LBM算法计算圆柱旋转问题是可靠的。

3.2 结果分析 3.2.1 流场结构

首先讨论不同相对转速对漩涡脱落形态的影响。如图 6给出了雷诺数Re=300时，不同转速条件下的瞬时涡量等值线图和流场图。

从图 6中可以看到在转速比较低时(如图 6(a)、(b))，旋涡排列较规则，正负涡的形状大小基本一致，逆时针的旋转导致脱落的旋涡上移，涡稍微被拉伸。这是因为对于静止圆柱来说，在分离点(流速为0的点)后，边界层发生脱离，进而产生周期性的旋涡，而这个旋转作用导致的速度增加使原有的环向流速增大，在圆柱上表面，这个速度与均匀来流方向相反，在圆柱下表面与来流方向相同，这样圆柱上表面分离点提前，而下表面分离点滞后，这就导致了边界层在脱离后形成的剪切层与均匀来流形成一定夹角，即旋涡上移。当转速增加后(如图(c))，旋涡被大大拉伸，圆柱下表面脱落的正涡比上表面脱落的负涡拉长效果更加明显。这是由于转速较大，旋转导致的剪切层与流动方向间夹角增加，在接近自由流区域外流速比内流速大的情况下，内外流速形成一定夹角，因此旋涡形状被拉伸。转速再增大(如图(d))，涡的长度逐渐减小，正涡大部分出现在靠近圆柱表面附近，而负涡覆盖在正涡外侧，并向外延伸。当α等于5时，圆柱产生的正负涡层层缠绕于圆柱表面。从流线图看，随着α的增加，流线逐渐上移贴近圆柱。在α较小时，没有出现环绕圆柱的流线，在α增大后，出现环绕在圆柱表面的流线，这也是转速增加涡交替附着圆柱表面造成的。在α从2逐步增加到5时，停滞点也发生了改变，旋转加快，停滞点逐渐上移，这是因为在圆柱周围的环绕密集流线在逐渐远离圆柱表面。

如图 7给出了不同相对转速条件下，在指定的流向位置(X/D=0.5，1.0，1.5，2.0)上的归一化流向速度。可以看到在接近圆柱表面处的归一化的速度U/U₀在(-2, 2)变化，当X/D向下游更远处延伸后，U/U₀的变化范围减小，在(0, 1)变化，这是因为距离圆柱越远，速度的损失逐渐减小，速度受柱体的影响减弱，速度逐渐恢复而接近来流速度。在不同位置处α对U/U₀的影响不同，X/D较小时，随着转速的增加U/U₀的波动增加，这是因为α增加，涡脱落变得不规则，正负涡逐步环绕在圆柱表面，脱落现象消失，这与图 6的结果相吻合。当远离圆柱后，随着α增加，U/U₀的波动先有一个小幅度增加，然后逐渐减小，这说明在α逐渐增加的过程中，其对于流速的影响有一个临界点，在临界点前后变化规律相反，此α在1~2。

图 8是雷诺数为300时在圆柱垂向方向上不同位置的涡量等值线。

圆柱垂向中点为z₀，图中(a)、(b)是相对转速为0.5时在z₀和z₀+0.2D位置时的涡量等值线，(c)、(d)是相对转速为1.0时在z₀和z₀+0.2D位置处的涡量等值线。可见圆柱尾流呈现三维的流动状态，在不同的垂向位置，脱落的漩涡形状不同，在圆柱垂向中点处漩涡的脱落速度较远离圆柱垂向中点脱落的速度要慢。

图 9是3种雷诺数下相对转速为1时y截面上的wy涡量云图。在雷诺数较低时，尾部涡脱落较为整齐，在z方向上正负涡分布均匀。当雷诺数增加时，在垂向方向上，圆柱尾涡发生扭曲和变形，流场紊乱缺乏规律性，正负涡环绕在一起，并且随着雷诺数的进一步增大，这种现象更加明显。在圆柱表面附近存在4个涡，沿着z方向一正一负交替排列。在远离圆柱表面过程中，在雷诺数100时，正负涡变成2个，而当雷诺数增大到200时，涡有2个和4个2种情况，部分正涡被负涡包裹，当雷诺数增加到300时，涡的分布展现出更为复杂的三维湍流特性。

3.2.2 受力特点

为了进一步研究旋转对圆柱绕流特性的影响，如图 10分别给出了雷诺数在300情况下，不同相对转速时的阻力系数均值以及升力系数的时间历程曲线。

在α逐步增加时，阻力系数均值减小，但是当α>1时，阻力系数均值缓慢的增加，且旋转越快增加的越快。当α较小时，升力系数曲线稳定后呈现周期性振荡，其幅值随着α的增加而减小。当α增加一定程度后，升力系数的周期性不复存在，曲线整体呈现一种紊乱的波动，并且α越大波动越剧烈，这是因为旋转速度的增加导致脱落的频率受到St数和α的共同影响从而失去周期性。此外α增加后升力系数均值也有大幅度增加，这体现了马格努斯效应，旋转速度增加圆柱会产生额外的升力，且随着转速的增加这个额外的升力也变大。

图 11是不同Re条件下相对转速为1时圆柱的升阻力曲线。雷诺数对于升阻力的振荡幅值和振荡平衡点大小都有影响。在Re=100时，周期性振荡的阻力系数C_d的幅值较小，但是振荡平衡点对应的值较大，当Re=200时，阻力系数振幅增加，而其振荡平衡点对应的数值减小，当Re=300后，振荡幅值减小，但其振荡平衡点对应的值却进一步减小。升力系数均方根的变化规律与阻力系数均值相反，随着Re的增加而逐渐增加，但是升力系数的振荡幅值先增加而后减小。

图 12给出了不同转速条件下升力系数的功率谱密度图。图中在0处的PSD值是由于升力系数均值不为0造成的。可以看到，随着转速增加，漩涡脱落的频率增加，这与图 10结果相吻合，在相对转速小于2时，漩涡仍呈现周期性脱落，转速越大，脱落频率越大，当相对转速大于2后，漩涡不再呈现规则的周期性脱落，因此PSD值趋近于0。

图 13给出了圆柱稳定状态下表面压力的分布。这里压力系数是根据圆柱受力达到稳定状态时即150~200 s区域压力均值计算得到的。圆柱在转速较低时，压力分布较为均匀，波动较小，在圆柱表面压力在-1.5左右，当转速增加时，压力系数沿着周向的分布逐渐变化，在圆柱上表面压力系数与下表面不同，因此产生了一个压力差，并且转速越大，这种差值越大，这与图 10结果相吻合。此外，还可以观察到随着转速的增加，压力系数的分布逐渐出现2个极值，分别发生在90°和270°附近。

图 14绘制了升阻力系数相位图，其中升阻力系数是在流动完全发展后的一个周期内计算得到的。这里相位图的变化是由于漩涡在一个周期内的脱落状态造成的。在相对转速为0时，在脱落周期内，力系数相位显示的是2个对称的半循环状态，这是因为柱体不旋转时，正涡和负涡分别从柱体上下表面脱落，流动是规则的。当相对转速不为0时，对称性消失，力系数相位图变成一个近似于圆形的图形，并且随着转速的增大其对应的|C_L|值增加而C_D值减小。此外在相对转速一定时，也可观察到雷诺数越大，C_D-C_L的相位对应的C_D值越小。

4 结论

柱体旋转运动是一种流体控制的方法，旋转作用能改变柱体的尾流场，产生的马格努斯效应具有很重要的工程意义，本文基于IB-LBM方法对3种雷诺数下，相对转速为0.5，1.0，2.0，3.5，5.0时的圆柱进行了数值模拟，分析了圆柱尾部流场涡量场变化以及升阻力的变化，得到如下结论：

1) IB-LBM方法在进行旋转边界处理、流场特性捕捉上具有一定的优势，能够清晰的模拟出移动边界的流场，当前的计算结果与文献中数据符合较好，计算结果具有良好的可信度。

2) 旋转作用会导致柱体尾部流场与流动方向间产生夹角，影响尾涡的周期性规则脱落，在大的相对转速时会导致正负涡交替环绕而不脱落现象。

3) 雷诺数影响尾流场的三维性，增加的雷诺数会导致尾涡场复杂紊乱，涡的扭曲现象更加明显

4) 相对转速增大会改变柱体尾部流速的变化趋势，且对于不同截面，转速对流速的变化范围影响不同，在远离柱体表面的截面上，相对转速对于流速的影响存在一个在1~2的临界点

5) 旋转作用会影响升阻力系数曲线的变化趋势，大的相对转速使升阻力曲线的周期性减弱，曲线图紊乱，升阻力相位关系由2个半循环转变成类圆形。升阻力的变化趋势受雷诺数影响，雷诺数的增大有助于提高升力，降低阻力，使得升阻力系数相位图对应的C_D值减小。