基于浸入边界-多松弛时间格子玻尔兹曼通量求解法的流固耦合算法研究∗

吴晓笛 刘华坪 陈浮

(哈尔滨工业大学能源科学与工程学院,哈尔滨 150001)

基于浸入边界-多松弛时间格子玻尔兹曼通量求解法的流固耦合算法研究∗

吴晓笛 刘华坪†陈浮

(哈尔滨工业大学能源科学与工程学院,哈尔滨 150001)

(2017年5月19日收到;2017年8月22日收到修改稿)

针对流固耦合问题,发展了基于浸入边界-多松弛时间格子玻尔兹曼通量求解法(immersed boundary method multi-relaxation-time lattice Boltzmann flux solver,IB-MRT-LBFS)的弱耦合算法.依据多尺度Chapman-Enskog展开,建立不可压宏观方程状态变量和通量与格子玻尔兹曼方程中粒子密度分布函数之间的关系;采用强制浸入边界法处理流固界面使固壁表面满足无滑移边界条件,根据修正的速度求解动量方程力源项;结构运动方程采用四阶龙格-库塔法求解.格子模型与浸入边界法的引入使流固耦合计算可以在笛卡尔网格下进行,无需生成贴体网格及运用动网格技术,简化了计算过程.数值模拟了单圆柱横向涡激振动、单圆柱及串列双圆柱双自由度涡激振动问题.结果表明,IB-MRT-LBFS能够准确预测圆柱涡激振动的锁定区间、振动响应、受力情况以及捕捉尾流场结构形态,验证了该算法在求解流固耦合问题的有效性和可行性.

10.7498/aps.66.224702

流固耦合,浸入边界法,格子玻尔兹曼通量求解法,涡激振动

1 引 言

流固耦合问题广泛存在于航空航天、海洋水利工程、石油运输及生物工程等领域中.依据耦合机理流固耦合问题可以分为两类[1]:流体、固体区域部分或完全重叠在一起,如渗流问题;流体和固体的耦合作用仅发生在两种介质交界面处,如涡激振动.依据控制方程解法流固耦合问题又可分为强耦合和弱耦合[2,3].强耦合将流体、固体域和耦合作用构造在同一控制方程中直接求解;弱耦合是在单一时间步内分别求解流体域和固体域,通过两相介质界面的数据交换实现流固耦合的计算,基于其概念简单,易于编程实现,因此得到广泛应用[4].

目前,对于流固耦合这种动边界问题常采用边界贴合方法处理,即根据每个时间步物体运动变形生成一套新的贴体网格,该方法能够很好地捕捉附面层特性,有利于模拟高雷诺数流动问题,但是网格生成的计算量将大大增加,降低整体程序的计算效率[5].近年来,非边界贴合方法越来越受到关注.该类方法往往通过对边界周围流场变量的恰当处理来实现边界条件,因此不严格要求网格与物体边界贴体,可以采用直角网格等规则的网格,使得网格的生成更容易.当物体运动时,可以直接在固定的网格上求解流动控制方程,从而无需在每个时间步上对网格进行变换、映射、重新生成等处理,大大提高了计算效率[6].因此,非边界贴合方法在模拟复杂流场和运动边界的流动时能够真正做到简便高效.

浸入边界法(immersed boundary method)作为非边界贴合方法之一,在流固耦合问题中得到了广泛的应用.该方法需要两套网格,流场使用固定直角坐标网格来求解,而浸入边界则用拉格朗日点来标识,物体与流场的作用通过两套网格之间的信息交互传递来完成.传统的浸入边界法主要分为直接力法[7]、惩罚力法[8]以及动量变化法[9].王文全等[10,11]应用浸入边界法实现了对刚体定轴转动的预测及水轮机活动导叶旋转摆动绕流的数值计算;明平剑等[12,13]采用浸入边界法及非结构化网格求解了低雷诺数下振荡圆柱绕流问题;李师尧等[14]采用浸没边界-格子玻尔兹曼格式模拟贯流式水轮机三维瞬变流.

若采用浸入边界法实现流固耦合问题的数值模拟,则需要一种有效的流场求解器.现阶段流场计算主要是通过求解Navier-Strokes(N-S)宏观方程或者格子玻尔兹曼方法(LBM).针对不可压控制方程,由于连续方程中压力变量的缺失,求解N-S方程的困难是对压力场的预测,需引入压力泊松方或者压力修正方法来求解,会影响整体的计算效率.而LBM中压力场可以由粒子分布函数直接给出,且具有算法简单、易于并行的优势,但还存在一些缺点,由于格子速度模型的均匀性且网格尺度与时间推进步长的关联性,造成该方法不能有效地计算工程应用问题,如计算湍流问题,边界层问题等;LBM中松弛时间项与黏性系数有关,因而该方法只能求解黏性流动问题,而N-S方程既可以求解黏性流动又可以求解无黏流动.基于以上分析,Shu等[15]开发了一种基于单松弛时间格式的格子玻尔兹曼通量求解器,其思想是基于有限体积法空间离散宏观方程,状态变量及通量由格子玻尔兹曼模型中粒子密度分布函数重新构造.该方法有效结合了两种求解器的优点.基于该思想,本文发展了一种基于多松弛时间格式的玻尔兹曼通量求解器(multi-relaxation-time lattice Boltzmann flux solver).将该求解流场的方法与Wang等[16,17]开发的强制浸入边界法相结合,开发了一种基于浸入边界-多松弛时间格子玻尔兹曼通量求解法(IB-MRT-LBFS)的弱耦合算法.该方法求解动边界问题无需使用动网格技术,节约了网格重新生成的计算量;强制浸入边界法的应用严格保证了无滑移边界条件;多松弛格子玻尔兹曼求解法改善了传统流体动力学计算方法必须花费大量时间求解泊松方程而得到压力场的限制,并且通量表达式直接由平衡分布函数给出,简化了流场的计算过程.

2 数值方法

2.1 2.1 Chapman-Enskog展开分析

在低马赫数条件下,借助Chapman-Enskog展开可以完成宏观Navier-Stokes方程与介观Boltzmann方程的相互转化,基本思想是将时空变量用多层次时空尺度表示,使其在各级尺度上物理量的量级保持一致.

宏观不可压Navier-Stokes控制方程表达式为:

其中,ρ,u,p和µ分别代表流体密度、速度、压力及动力黏性系数.

采用 Bhatnagar-Gross-Krook(BGK)碰撞模型的多松弛时间格子玻尔兹曼方程表示如下:

其中,α表示不同离散速度粒子;eα是粒子速度;fα是沿α方向的离散速度密度分布函数;r代表物理位置;M是正交转换矩阵,其作用是将速度空间密度分布函数fα及其平衡分布函数转换为多松弛时间尺度下的动量空间分布函数mα=Mfα及其平衡分布函数S=diag(s0,s1,s2,s3,s4,s5,s6,s7,s8)为对角碰撞松弛矩阵,由于在碰撞过程中为了保证质量和动量守恒,s0=s3=s5=0,s7=s8=1/τ,松弛时间τ是与动力黏性系数µ相关的物理量,取其中,cs为LBGK模型的声速.

采用D2Q9格子玻尔兹曼模型,离散速度定义为

其中,c通常取值为1.平衡分布函数定义为

其中,w0=4/9,w1=w2=w3=w4=1/9,

宏观方程中密度与动量可以用密度分布函数表示为

通过多尺度展开,分布函数、时间及空间导数可以表示为:

其中,ε为与克努森数成正比的参量.

将方程(3)进行二阶泰勒展开,得到以下方程:

将方程(7)—(10)代入方程(11)得到

其中I为单位矩阵.根据方程(12)和(13)推导出:

通过对方程(13)和(14)零阶矩和一阶矩的求和得到以下方程:

其中

通过以上公式的推导,建立了宏观方程通量与格子Boltzmann分布函数的关系,如(17)和(18)式所示.

2.2 流体控制方程

基于浸入边界法的宏观不可压Navier-Stokes控制方程表达式为

其中F是由浸入边界法求解的力源项.

根据2.1节中建立的宏观方程状态变量和通量与格子玻尔兹曼方程中密度分布函数之间的关系((6),(17),(18)式),采用IB-MRT-LBFS的不可压控制方程表达式为:

其中:

求解流体控制方程(23)和(24)分为两个步骤:

1)预测步,不考虑物面边界,采用格子Boltzmann通量求解法(LBFS)求解流场变量,即密度ρn+1和瞬时速度u∗,

2)修正步,采用强制浸入边界法(IBM)求解修正后的速度变量un+1,进一步求解边界作用力,

其中,网格界面通量重构与边界处理是数值计算中两个关键部分.

2.2.1 格子Boltzmann通量求解法

为了求解方程(27)和(28)中ρn+1和u∗,采用LBFS计算相邻单元界面的通量.基于有限体积法的控制方程离散形式:

其中W=[ρ,ρu]T;dVi为控制体的体积,dSk为控制体第k个控制面,n为控制面上的法向向量.

基于格子玻尔兹曼D2Q9模型的界面通量重构示意图见图1.已知相邻两个网格单元和界面的物理位置,分别定义为ri,ri+1以及r,则有

其中ψ代表流体变量ρ或u.

图1 界面通量重构示意图Fig.1.Local flux reconstruction at cell interface.

基于泰勒展开计算图1中离散点的速度,重构界面处密度和动量,表达式为:

通量表示为

其中,平衡项可以根据界面处宏观变量求解,而非平衡项可以用平衡项近似求解,表达式为:

2.2.2 强制浸入边界法

浸入边界法在数学建模过程中采用欧拉变量描述流场网格点,采用拉格朗日变量表示物面边界,如图2所示.应用强制浸入边界法进行流场速度的修正,使得物面满足无滑移边界条件.u∗为预测步求解的流场瞬时速度,∆u为修正速度,则修正后的欧拉点速度表达式为

图2 (网刊彩色)计算域二维示意图Fig.2.(color online)A solid boundary immersed in a two-dimensional computational domain.

欧拉点上修正速度可以通过Dirac delta函数插值求解,表达式如下:

其中rj为欧拉点坐标;分别为拉格朗日点坐标以及修正速度;D是连续Kernel函数,表达式为:

为了保证无滑移边界条件,拉格朗日点的速度等价于同一位置流场中的速度,通过插值来求解拉格朗日点的流动速度,即

其中,l=1,···,N;i=1,···,M;h为欧拉点网格尺寸,N和M分别为拉格朗日点与欧拉点的个数,则拉格朗日点上的速度为

将上述方程写成矩阵的形式:

其中,

求解矩阵方程AX=B获得拉格朗日点的修正速度通过(38)式求得欧拉点的修正速度∆u,进而求得修正后的流场瞬时速度un+1及边界作用力F.

2.3 IB-MRT-LBFS算法流程

IB-MRT-LBFS流固耦合算法步骤归纳如下:

1)给定初始状态值,选取迁移时间步长δt,计算碰撞松弛矩阵S;

3)根据(32)和(33)式求解界面处流场宏观变量,进而计算

4)根据(34)式计算界面通量,求解预测步中网格中心宏观变量;

5)依据不同情况下刚体运动方程,计算修正步中拉格朗日点的运动速度及坐标解矩阵方程(43)获得拉格朗日点修正速度;

6)根据(38)式求解流场修正速度,将其代入(37)式计算出修正后的流场宏观速度变量;根据(29)式求解边界作用力,完成修正步的计算;

7)重复步骤2)—6)至计算程序收敛.

3 数值计算与分析

利用IB-MRT-LBFS流固耦合算法开展单圆柱及串列双圆柱涡激振动问题的数值模拟研究.设圆柱为弹性支撑,双自由度无量纲运动方程表达式:

其中,X和Y分别为x,y方向上的无量纲位移,其一阶导数和二阶导数分别表示其相应方向上的无量纲速度和加速度;ζ为结构阻尼系数;质量比mr=4m/(πρD2);折合固有频率fr=fnD/U∞,fn为结构固有频率;折合速度Ur=U∞/(fnD);升阻力系数{CL,CD}={FL,FD}/(0.5ρu2D).升阻力由应力积分法求解得表达式为

然而,当采用浸入边界法处理运动刚体的物面条件时,由于刚体内部受到内部流体的作用力,因此应该考虑消除内部作用力的影响[18],(49)式修正为

3.1 单圆柱横向涡激振动

单圆柱钝体结构横向涡激振动问题作为经典验证算例已有较多的文献报道,其中Ahn和Kallinderis[19]采用ALE动网格技术数值模拟了流固耦合问题;Borazjani等[20]采用浸入边界法中的直接力法求解曲线边界的流固耦合问题;Jiang等[21]采用外插法处理运动边界问题;Wang等[22]采用ALE动网格与双重网格相结合的方法求解动边界的问题.应用IB-MRT-LBFS模拟雷诺数为150的单圆柱单自由度涡激振动,取质量比mr=8/π,为了使结构振动幅度最大,阻尼系数设置为ξ=0,计算折合速度为3≤Ur≤8.计算域及边界条件如图3所示,计算域入口距圆柱中心距离为40D且为Dirichlet型边界条件,密度为ρ=1.0,速度为U∞=0.1;上下边界距圆柱中心位置均为40D且为无滑移边界条件;出口距圆柱中心位置均为60D且为Neumann型边界条件;圆柱表面为无滑移边界条件,用均匀分布150个拉格朗日点代替,浸没在160×240的均匀欧拉网格区域中;计算域总直角网格数为547×681.

图3 单圆柱横向涡激振动计算域及边界条件示意图Fig.3.The computational domain and boundary conditions for a VIV cylinder in transverse direction.

图4给出了在不同折合速度下IB-MRT-LBFS流固耦合算法计算的无量纲横向最大振幅与文献数值结果的比较.由图可知,本文结果与现有文献[19—22]中采用不同数值方法的计算结果趋势相一致.当Ur=4时,单圆柱单自由度涡激振动进入“锁定区间”,尾流脱落涡频率接近圆柱固有频率,发生共振现象,振幅增加;当Ur=8时,单圆柱涡激振动超出“锁定区间”,圆柱振幅降低.

图4 (网刊彩色)不同折合速度单圆柱横向涡激振动最大振幅Fig.4.(color online)Maximum amplitude against reduced velocity for a VIV cylinder in transverse direction.

不同的数值计算方法对单圆柱涡激振动横向最大振幅结果影响较大,尤其是在非锁定区间,例如在Ur=8时文献[19,20]中计算结果相对误差达到40%以上.由图4可以看出,数值算法模拟结果介于不同文献数值结果之间,且可以准确地捕捉锁定区间范围及预测共振区域的振幅值.表1给出了单圆柱横向涡激振动最大振幅的误差分析,可以看出本文数值结果与文献结果相比较,多数误差在10%以内,且锁定区间内的误差较小,说明IB-MRT-LBFS流固耦合算法在计算涡激振动问题上的有效性.

图5给出了不同折合速度下单圆柱横向涡激振动升力系数与无量纲最大振幅随时间历程的变化情况.当Ur=3,4,5时,升力与振幅相位相同;而当Ur=6,7,8时,二者相位相反.图6给出了不同折合速度单圆柱横向涡激振动瞬时涡量图,可以清晰地观察到尾流场脱落涡结构特性.当Ur=3时,尾流模态为单行涡街;当Ur=4时,尾流模态呈现不规律性;随着折合速度进一步增加到Ur=5时,尾流模态出现双行涡街;当Ur≥6时,尾流又恢复到单行涡街形态.综上所述,单圆柱横向涡激振动响应规律以及尾流特性与文献[23]中相一致,说明IB-MRT-LBFS流固耦合新型算法能准确的求解单圆柱单自由度涡激振动问题.

图5 (网刊彩色)不同折合速度单圆柱涡激振动横向升力系数与振幅的时程变化Fig.5.(color online)Lift coefficient and amplitude against reduced velocity for a VIV cylinder in transverse direction.

表1 单圆柱横向涡激振动最大振幅误差分析Table 1.The error analysis on maximum amplitude for a VIV cylinder in transverse direction.

图6 (网刊彩色)不同折合速度单圆柱横向涡激振动瞬时涡量Fig.6.(color online)Instantaneous vorticity against reduced velocity for a VIV cylinder in transverse direction.

3.2 单圆柱双自由度涡激振动

采用IB-MRT-LBFS数值算法模拟60≤Re≤140条件下单圆柱双自由度涡激振动,与文献[24]中的计算参数相一致,取质量比mr=10.0,阻尼系数ξ=0,计算域、边界条件及网格与单圆柱横向涡激振动相同.

图7给出了在不同雷诺数条件下振幅及作用力的数值计算结果,且与文献[24]结果相对比,可以看出各项结果符合良好;对比图7(a)与图7(b),可以看出横向振幅值远远大于轴向振幅,当Re=90时,单圆柱进入锁定区间,产生共振,振幅达到最大值;随着雷诺数的增加,振幅值逐渐减小;当Re=130时,圆柱驶出锁定区间.从图7(c)与图7(d)可以看出升阻力均方根值的变化情况,升力系数均方根值与横向振幅的变化规律不同,首先随雷诺数的增大而增加;当进入“锁定区间”时(Re=90)达到最大值,随之逐渐减小;当驶出“锁定区间”时(Re=130)时,升力均方根值突然增大;阻力系数均方根值与轴向振幅的变化规律相一致.表2给出了横向最大振幅的误差结果,可以看出相对误差结果均在10%以下,进一步说明IB-MRT-LBFS数值算法能准确地预测低雷诺数条件下圆柱涡激振动问题.

图8给出了不同雷诺数下单圆柱涡激振动瞬时涡量图.当Re=60,70和80时,尾流模态为单行涡街;当Re=90和100时,尾流模态为双行涡街;当Re≥110时,尾流又恢复到单行涡街形态;当雷诺数增加到140时尾流模态由单行涡街逐渐演化为双行形态.以上尾流场脱落涡结构特性与文献[24]中观察到的尾流形态相一致,说明IB-MRT-LBFS能够有效预测流场结构特性.

图7 (网刊彩色)振幅及作用力系数随雷诺数的变化情况 (a)横向最大振幅;(b)轴向振幅均方根;(c)升力系数均方根;(d)阻力系数均方根Fig.7.(color online)Amplitude and force coefficients against Reynolds number:(a)Max transverse amplitude;(b)the rms value of in-line amplitude;(c)the rms value o flift coefficient;(d)the rms value of drag coefficient.

表2 单圆柱双自由度涡激振动横向最大振幅误差分析Table 2.The error analysis on maximum transverse amplitude for a VIV cylinder.

图8 (网刊彩色)不同雷诺数下单圆柱涡激振动瞬时涡量Fig.8.(color online)Instantaneous vorticity against Reynolds number for a VIV cylinder.

3.3 串列双圆柱双自由度涡激振动

串列双圆柱双自由度涡激振动计算域及边界条件如图9所示,两圆柱的质量和直径相同,圆心间距L/D=5;雷诺数为150;质量比mr=8/π;阻尼系数ξ=0.边界条件与单圆柱涡激振动相同,两圆柱表面均为无滑移边界条件,固壁边界同样用150个拉格朗日点代替,浸没在640×320的均匀欧拉网格区域中;计算域总直角网格数为1060×740.分别计算折合速度3≤Ur≤11情况下串列双圆柱双自由度涡激振动问题.

图9 串列双圆柱涡激振动计算域及边界条件示意图Fig.9.The computational domain and boundary conditions for two VIV cylinders in tandem.

图10给出了不同折合速度下串列双圆柱系统上下游圆柱轴向及横向最大振幅结果;图11给出了不同折合速度下串列双圆柱系统上下游圆柱受力情况.由图可知,本文数值方法模拟结果与现有文献[25,26]中结果符合良好,说明IB-MRT-LBFS流固耦合算法能够准确地预测多钝体结构双自由度涡激振动问题.

由图10(a)可知,在折合速度3≤Ur≤11情况下,上游圆柱轴向振动最大振幅小于0.05,且在Ur=4时轴向位移最大;当Ur≥6时随着折合速度的增加,轴向位移呈现先减小后增加的趋势.由图10(c)可知,下游圆柱轴向位移无显著变化规律,这是由于上游圆柱形成的尾迹脱落涡结构作用在下游圆柱表面,导致下游圆柱所受阻力有所变化.当Ur≥9时,下游圆柱轴向振幅突然增加,最大值达到0.34;说明在较高折合速度条件下,上游圆柱尾迹作用使得下游圆柱轴向振幅加大.

由图10(b)可知,上游圆柱横向振幅随折合速度的增加先变大后减小直至趋于平缓;且最大振幅出现在Ur=5时,约为0.6;对比单圆柱横向振动结果,上游圆柱的振动响应与单圆柱相似,但是对于串列双圆柱系统,在下游圆柱的作用下,上游圆柱的锁定区域有所增加.由图10(d)可知,下游圆柱横向振幅随折合速度变化趋势与上游圆柱基本相同;当Ur≥6时下游圆柱在上游尾迹的作用下振幅增加,均大于上游圆柱振幅,这是由于此折合速度下,下游圆柱升力均方根值大于上游圆柱(见图11(c)和图11(f));最大振幅出现在Ur=8时,约为1.06;对于串列双圆柱系统,上游圆柱的作用推迟下游圆柱进入锁定区间,且扩大了下游圆柱锁定区间的范围.

图10 (网刊彩色)最大振幅随折合速度的变化情况 (a)上游圆柱(轴向);(b)上游圆柱(横向);(c)下游圆柱(轴向);(d)下游圆柱(横向)Fig.10.(color online)Max amplitude against reduced velocity:(a)Upstream cylinder(in-line direction);(b)upstream cylinder(transverse direction);(c)downstream cylinder(in-line direction);(d)downstream cylinder(transverse direction).

图11 (网刊彩色)作用力系数随折合速度的变化情况 (a)上游圆柱平均阻力系数;(b)上游圆柱阻力系数均方根;(c)上游圆柱升力系数均方根;(d)下游圆柱平均阻力系数;(e)下游圆柱阻力系数均方根;(f)下游圆柱升力系数均方根Fig.11.(color online)Force coefficients against reduced velocity:(a)The mean drag coefficient,upstream cylinder;(b)the rms value of drag coefficient,upstream cylinder;(c)the rms value of lift coefficient,upstream cylinder;(d)the mean drag coefficient,downstream cylinder;(e)the rms value of drag coefficient,downstream cylinder;(f)the rms value of lift coefficient,downstream cylinder.

由图11(a)和图11(b)可知上游圆柱阻力系数的变化规律:当3≤Ur≤5时,阻力系数均值及均方根随折合速度的增加而增大;当5＜Ur≤9时,阻力系数随折合速度的增加而减小;随着折合速度的继续增加,阻力系数趋于平缓且趋近于单圆柱阻力.由图11(c)可知上游圆柱升力系数随折合速度的变化规律:升力系数均方根随折合速度的增加呈现先增大后减小直至Ur=8;当Ur≥9时,升力系数均方根趋于平缓且趋近于单圆柱升力.由图11(d)—(f)可知,下游圆柱所受作用力的变化规律不同于上游圆柱,这是由于下游圆柱受到来自上游圆柱不同结构尾迹形式的作用引起的.对比上下游圆柱所受阻力系数平均值,上游圆柱最大出现在Ur=5,而下游圆柱为Ur=7.

通过图12不同折合速度串列双圆柱涡激振动瞬时涡量可以观察串列双圆柱尾迹模态.在圆心间距L/D=5的情况下,上游圆柱在不同折合速度条件下均可以形成周期性的脱落涡结构,但是由于折合速度与脱落涡频率有关,因此形成的脱落涡的频率与涡长有所改变,作用在下游圆柱表面的涡结构不同进而影响下游圆柱表面作用力的分布及大小;观察下游圆柱尾迹结构,当Ur=3时,下游圆柱尾迹为双行涡街;当Ur=4时下游圆柱尾迹模态为单行涡街;当Ur=5,6时尾流模态呈现不规律性;随着折合速度进一步的增加,尾流模态出现双行涡街,当Ur=11时尾迹结构呈不规律性.

图13为串列双圆柱系统上下游圆柱位移轨迹对比结果;红色曲线表示上游圆柱运动轨迹,黑色曲线表示下游圆柱运动轨迹.由图可知,当Ur=3时上下游圆柱运动轨迹均呈现“8”字形,在此折合速度条件下,上游圆柱未进入“锁频区”,振动较弱,而下游圆柱在上游圆柱形成的周期性尾迹的作用下振动有所加强,双自由度的幅值均大于上游圆柱;当Ur=4时,上下游圆柱运动轨迹均呈现不规律性,在此折合速度条件下,上游圆柱进入“锁频区”,发生共振现象,振幅增加,此时上游圆柱脱落涡接近圆柱固有频率,因此下游圆柱在尾迹的作用下振动加强;随着折合速度的增加至Ur=8,上游圆柱在“锁频区”且均呈现“8”字形轨迹;当Ur=7,8时,下游圆柱轨迹为“8”字形;随着折合速度进一步增加,上游圆柱驶出“锁频区”,振动减弱,而下游圆柱在上游圆柱尾迹的作用下保持较强的振动状态,运动轨迹呈现非规律性.

图12 (网刊彩色)不同折合速度串列双圆柱涡激振动瞬时涡量Fig.12.(color online)Instantaneous vorticity against reduced velocity for two VIV cylinders in tandem.

图13 (网刊彩色)串列双圆柱上下游圆柱位移轨迹Fig.13.(color online)Displacement trajectories of the upstream and downstream cylinder centers with various reduced velocities.

4 结 论

发展了一种基于浸入边界-多松弛时间格式格子玻尔兹曼通量求解法的弱耦合流固算法(IBMRT-LBFS).对于流体控制方程,基于多尺度Chapman-Enskog展开,建立宏观方程状态变量及通量与格子玻尔兹曼方程中密度分布函数之间的关系,开发了多松弛时间格式的格子玻尔兹曼通量求解器;在流体方程中引入力源项,在保证无滑移边界条件下采用强制浸入边界法求解边界作用力大小;对于结构振动方程,采用四阶龙格-库塔法求解.该算法采用有限体积法进行空间离散,使得计算过程可以在非均匀网格下进行,对比传统的格子Boltzmann方法减少了网格量;控制方程中状态变量及压力项由粒子密度分布函数直接给出,网格界面通量重构使得宏观方程中对流项与黏性项可以同时由平衡分布函数求得,通量求解过程得以简化;强制边界法处理运动边界问题时,使得计算网格是固定的,无需再每一时间步更新网格的尺寸及位置,减少网格重新生成的计算时间.

数值计算和分析单圆柱横向涡激振动、单圆柱及串列双圆柱双自由度涡激振动问题.对于单圆柱涡激振动,横向振动响应大于轴向振动,当振动发生在锁定区间,尾流脱落涡频率接近圆柱固有频率,产生共振,振幅较大.对于串列双圆柱涡激振动,在较大折合速度条件下,下游圆柱振动响应大于上游圆柱,上游圆柱的作用力接近于单圆柱振动情况;上游圆柱推迟下游圆柱进入锁定区间且增加区间范围.通过对比现有文献数值结果,表明IB-MRT-LBFS能够准确预测圆柱涡激振动的锁定区间、振动响应、受力情况以及捕捉尾流场结构形态,验证了文中所提出的数值算法在计算流固耦合问题上的准确性与可行性.

[1]Xing J T,Zhou S,Cui E J 1997Adv.Mech.27 19(in Chinese)[邢景棠,周盛,崔尔杰 1997力学进展 27 19]

[2]Qian R J,Dong S L,Yuan X F 2008Spat.Struct.14 3(in Chinese)[钱若军,董石麟,袁行飞 2008空间结构 14 3]

[3]Guo P,Liu J,Wu W H 2013Chin.J.Theor.Appl.Mech.45 283(in Chinese)[郭攀,刘君,武文华 2013力学学报45 283]

[4]Zhou D,He T,Tu J H 2012Chin.J.Theor.Appl.Mech.44 494(in Chinese)[周岱,何涛,涂佳黄 2012力学学报44 494]

[5]Zhong G H,Liang A,Sun X F 2007J.Eng.Thermophys.28 399(in Chinese)[钟国华,梁岸,孙晓峰 2007工程热物理学报28 399]

[6]Liu Q Y 2012M.S.Dissertation(Nanjing:Nanjing University of Aeronautics and Astronautics)(in Chinese)[刘齐迎2012硕士学位论文(南京:南京航空航天大学]

[7]Luo H X,Dai H,Ferreira D S,Paulo J S A,Yin B 2012Comput.Fluids56 61

[8]Feng Z,Michaelides E 2004J.Comput.Phys.195 602

[9]Chen Y,Cai Q D,Xia Z H,Wang M,Chen S Y 2013Phys.Rev.E88 013303

[10]Wang W Q,Zhang G W,Yan Y 2017Trans.Beijing Inst.Tech.37 151(in Chinese)[王文全,张国威,闫妍2017北京理工大学学报37 151]

[11]Wang W Q,Su S Q,Yan Y 2015Chin.J.Comput.Mech.32 560(in Chinese)[王文全,苏仕琪,闫妍2015计算力学学报32 560]

[12]Ming P J,Zhang W P 2009Chin.J.Aeronaut.22 480

[13]Ming P J,Zhang W P,Lu X Q,Zhu M G 2010Chin.J.Hydrodyn.25 321(in Chinese)[明平剑,张文平,卢熙群,朱明刚2010水动力研究与进展25 321]

[14]Li S Y,Cheng Y G,Zhang C Z 2016J.Huazhong Univ.Sci.Tech.(Natural Science Edition)44 122(in Chinese)[李师尧,程永光,张春泽2016华中科技大学学报(自然科学版)44 122]

[15]Shu C,Wang Y,Teo C J,Wu J 2014Adv.Appl.Math.Mech.6 436

[16]Wang Y,Shu C,Teo C J,Wu J 2015J.Fluids Struct.54 440

[17]Wang Y,Shu C,Yang L M,Sun Y 2017Int.J.Numer.Meth.Fluids83 331

[18]Suzuki K,Inamuro T 2011Comput.Fluids49 173

[19]Ahn H T,Kallinderis Y 2006J.Comput.Phys.219 671

[20]Borazjani I,Ge L,Sotiropoulos F 2008J.Comput.Phys.227 7587

[21]Jiang R J,Lin J Z,Chen Z L 2013Phys.Rev.E88 023009

[22]Wang C L,Tang H,Duan F,Yu S C M 2016J.Fluids Struct.60 160

[23]Han Z L,Zhou D,Tu J H 2014J.Eng.Mech140 04014059

[24]Prasanth A K,Mittal S 2008J.Fluids Mech.594 463

[25]Bao Y,Huang C,Zhou D,Tu J H,Han Z L 2012J.Fluids Struct.35 50

[26]Yu K R,Etienne S Scolan,Y M,Hay A,Fontaine E,Pelletier D 2016J.Fluids Struct.60 37

PACS:47.63.mf,47.11.–j,47.11.Qr,43.40.–rDOI:10.7498/aps.66.224702

*Project supported by the National Natural Science Foundation of China(Grant No.51306042).

†Corresponding author.E-mail:hgdlhp@163.com

A method combined immersed boundary with multi-relaxation-time lattice Boltzmann flux solver for fluid-structure interaction∗

Wu Xiao-DiLiu Hua-Ping†Chen Fu

(School of Energy Science and Engineering,Harbin Institute of Technology,Harbin 150001,China)

19 May 2017;revised manuscript

22 August 2017)

This paper performs a newly developed method,which combines the immersed boundary method(IBM)with multi-relaxation-time lattice Boltzmann flux solver(MRT-LBFS),for solving fluid-structure interaction problems.Finite volume discretization is used to solve the macroscopic governing equations with the flow variables de fined at cell centers.Based on the multi-scale Chapman-Enskog expansion analysis,LBFS builds a relationship between the variables and fluxes in incompressible Navier-Stokes equations and density distribution functions in lattice Boltzmann equation.In order to ensure no-slip boundary condition,boundary condition-enforced immersed boundary method is used to treat the fluid-structure interface.The restoring force can be resolved by making a velocity correction in the flow field.The four-stage RungeKutta scheme is used to solve the motion equation of structure.Using the lattice model and immersed boundary method, fluid-structure coupling calculation can be implemented in a Cartesian grid,without generating the body- fitted mesh and using moving mesh technique.Therefore,the computational process is considerably simpli fied.In order to verify the validity and feasibility of IB-MRT-LBFS to solve fluid-structure interaction problems,both oneand two-degree of freedom vortex-induced vibrations(VIV)of a circular cylinder and two-degree of freedom VIV of two cylinders in a tandem arrangement are simulated by this proposal method.For a VIV cylinder system,the transverse vibration response is much stronger than the axial response.When the vibration occurs in the range of lock-in regime,the shedding vortex frequency of the wake is close to natural frequency of the cylinder so that resonance appears,consequently causing larger amplitude.For two VIV cylinders in a tandem arrangement,the dynamic behavior of each cylinder is signi ficantly di ff erent from that of a single cylinder.The gap spacing between the two cylinder centers is a signi ficant parameter which e ff ects vibration characteristics and the spacing is fixed in the simulations of two tandem cylinders.With the e ff ects of upstream cylinder wake,the axial and transverse amplitudes of downstream cylinder obviously increase with adding the reduced velocity.The downstream cylinder is delayed,coming into lock-in regime,and the range of lock-in regime is expanded under the e ff ects of the wake of the upstream cylinder.As the reduced velocity is relatively large,the vibration response of the upstream cylinder is close to a single cylinder and the vibration response of the downstream cylinder is more intense than the upstream cylinder.Compared with the existing literature results,our result illustrates that IB-MRT-LBFS owns the ability to correctly predict the lock-in regime,dynamic response and the forces of vortex-induced vibrations of cylinders.And this method can accurately capture the wake structures.

fluid-structure interaction,immersed boundary method,lattice Boltzmann flux solver,vortex-induced vibration

∗国家自然科学基金(批准号:51306042)资助的课题.

†通信作者.E-mail:hgdlhp@163.com