Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
  • ISSN 1001-1455  CN 51-1148/O3
  • EI、Scopus、CA、JST收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊

一种模拟动边界绕流的锐利界面浸入边界法

郭涛 张晋铭 张纹惠 王文全

郭涛, 张晋铭, 张纹惠, 王文全. 一种模拟动边界绕流的锐利界面浸入边界法[J]. 爆炸与冲击, 2022, 42(8): 084201. doi: 10.11883/bzycj-2022-0342
引用本文: 郭涛, 张晋铭, 张纹惠, 王文全. 一种模拟动边界绕流的锐利界面浸入边界法[J]. 爆炸与冲击, 2022, 42(8): 084201. doi: 10.11883/bzycj-2022-0342
GUO Tao, ZHANG Jinming, ZHANG Wenhui, WANG Wenquan. A sharp-interface immersed boundary method for simulating flows around bluff body with moving boundary[J]. Explosion And Shock Waves, 2022, 42(8): 084201. doi: 10.11883/bzycj-2022-0342
Citation: GUO Tao, ZHANG Jinming, ZHANG Wenhui, WANG Wenquan. A sharp-interface immersed boundary method for simulating flows around bluff body with moving boundary[J]. Explosion And Shock Waves, 2022, 42(8): 084201. doi: 10.11883/bzycj-2022-0342

一种模拟动边界绕流的锐利界面浸入边界法

doi: 10.11883/bzycj-2022-0342
基金项目: 国家自然科学基金(51969009,52179087,52069010)
详细信息
    作者简介:

    郭 涛(1983- ),男,博士,副教授,guotaoj@126.com

    通讯作者:

    王文全(1977- ),男,博士,教授,wwqquan@126.com

  • 中图分类号: O351

A sharp-interface immersed boundary method for simulating flows around bluff body with moving boundary

  • 摘要: 为避免复杂贴体网格的更新和畸形对动边界流场计算效率、精度的影响,以充分掌握结构场的受力特性,采用一种改进的锐利界面(sharp-interface)浸入边界法模拟具有动边界绕流的流动问题。该方法将计算域中的固体视为流体,固体边界离散为若干个拉格朗日网格点,通过在界面单元处插值重构流动参数(速度),将其直接作为流动求解器的边界条件,由此来反映固体边界的影响。即通过构造“虚拟点—受力点—垂足点”的计算结构,借助双线性插值得到虚拟点的速度,再通过强制满足固体边界的无滑移条件计算出受力点的速度,以此为边界条件,最终求解基于浸入边界法的耦合系统方程,实现复杂动边界的流动数值模拟。采用C++编写该浸入边界法的数值程序,以单圆柱绕流为验证算例,通过与文献和实验结果的对比,验证了该方法的准确性和可靠性。在此基础上,对主动运动椭圆柱绕流问题进行了精细计算,探讨了不同轴长比(AR)、不同攻角(θ)下的椭圆柱对尾涡结构分布特征和水力不稳定现象的影响。捕捉到了反对称S型、“P+S” Ⅰ型、“P+S” Ⅱ型尾涡脱落模态,漩涡强度、涡脱频率和升阻比随ARθ的变化规律,以及确定了升阻比临界攻角(25°)。
  • 绕流是自然界以及实际工程运用中一个十分普遍的物理现象,也是一个经典的流体力学问题。例如:机翼在流动空气中的振动问题[1-3];导弹无侧滑大攻角飞行时,背风面分离涡的不对称运动导致的附加偏航力问题;水轮机导叶、叶片在水流作用下的摆动和振动问题[4-5];鱼在水中通过鱼鳍的摆尾提供推力[6];昆虫利用翅膀的拍动使其在空中飞行[7]等等。这些运动及工程都有涉及到翼型、椭圆在流场中的绕流问题。当流体绕过不具有流线型的钝体结构时,在钝体尾部会发生周期性的漩涡脱落,产生周期性的上下拖曳力,也就是升、阻力。随着来流情况的变化,当涡脱频率接近结构的自振频率时,漩涡脱落和结构振动相互锁定,形成共振,也就是涡激振动,其在实际生活中产生的危害十分常见并且影响也十分广泛。例如:1940年刚通车数月的塔科马大桥发生颤振风毁事件[8],震惊世界;以及近期的虎门大桥涡激振动事件[9]等。因此,流固耦合问题一直是计算流体力学领域较为关注的一类重要问题。

    一般地,研究流固耦合问题时,通常采用贴体网格,这类网格质量良好且生成速度快,但是在实际计算中,为了适应不断变化的边界条件,在每一个时间步中都要重新划分网格,这样就大大增加了计算成本。鉴于这些原因,Peskin[10]提出了浸入边界法,并应用到模拟人类心脏血液流动中的瓣膜运动中。浸入边界法(immersed boundary method,IBM)在整个物理计算区域中采用一套笛卡尔网格,避免了在每一个时间步更新网格的繁琐,使用力场来代替固体边界,使流固耦合问题中的动边界处理起来更加简单和精确。因此,力源项的获得是浸入边界法求解流体问题的关键。Peskin方法,早期是将固体边界离散为一组由弹簧连接的单元,当边界发生变形时,就会产生一个意图使其回到原位置的回复力,通过边界上的回复力结合delta函数求得力源项。但该方法涉及到固体的弹性变形,主要适用于柔性体,而且该变形不易通过流场物理量求得,不具有一般性。Goldstein等[11]在此基础上发展了虚拟边界法来求解刚性边界问题,主要使用反馈力法,配合谱方法来渐进地施加所需的边界条件,该法也称反馈力法。在浸入边界法的发展过程中,Jamaludin等[12]提出了直接力法,即通过直接构造受力区的速度场来推导力源项,从而确保满足边界条件,基本实现了对边界的准确描述,成为学者们的主要研究方向之一。如Fadlum等[13]将其应用于流固耦合研究,获得了与贴体网格下计算结果一致性较好的结果;Le等[14]提出了一种隐式直接力浸入边界法,Wu等[15]、王文全等[16]提出了一种对速度进行二次修正的隐式直接力浸入边界法;Uhlmann[17]等提出了显式直接力法。直接力法相比于Peskin的经典IBM法,无须通过繁琐的delta函数[18,19]完成拉格朗日点与欧拉点上信息交换,通过直接施加浸入边界处的速度边界,由此推导出力源项,不依赖于固体的本构关系,也不进行反馈调整,故名直接力法。

    按照是否使用delta函数或能否对界面进行清晰描述,可将IBM法分为扩散界面(diffused-interface)浸入边界法和锐利界面(sharp-interface)浸入边界法。理论上来讲,经典浸入边界法、反馈力法和直接力法均属于扩散界面浸入边界法。锐利界面法[20-23]不直接计算力源项来施加边界条件,主要是在界面处通过局部插值,重构边界网格点附近的流动参数(速度),以此作为边界条件将其直接施加到计算域进行求解。理论上可以具有高阶精度[24-25],在可压缩流体[26]以及激波冲击[27]问题中均得到了利用。本文将亦采用一种改进的锐利界面浸入边界法,将整个物理区域划分成纯流体区域以及包含固体的次流体区域。在流体次区域边界的法向,构造“虚拟点—受力点—垂足点”的计算结构。受力点速度由流体域内某一虚拟点和边界上垂足点之间的速度线性插值确定,实现复杂动边界的流动数值模拟。本文利用C++编程,通过圆柱绕流经典算例与其他方法的文献结果及试验结果进行对比,验证该方法的有效性和准确性。基于该方法,进一步探究不同轴长比AR、不同攻角θ下的椭圆柱绕流的涡结构分布特征、流线、压力等流场现象和升阻力系数、涡脱频率等水力不稳定现象。

    将整个物理区域(包括流体、固体)看作不可压缩牛顿流体的粘性流动。在整个计算域内,流体采用Euler描述,固体采用Lagrange描述。其连续性方程和动量方程可表示为

    u=0 (1)
    ut + (u)u=1ρp+μρ2u+f (2)

    式中:u代表速度矢量(u,v),p代表流体的压力,μ为动力黏度, ρ为密度,f表示力源项。

    采用二阶时间分布投影法对时间进行离散,来求解控制方程(1)、(2),具体步骤如下。

    (1)忽略力源项,计算考虑初始压力作用下的中间速度u*,即:

    u*unΔt=32Hn12Hn1pn (3)

    式中:

    H=uut+vut1Re(2ux2+2uy2) (4)

    (2)忽略压力项,考虑计算所得力源项作用下的中间速度u,即:

    uunΔt=32Hn12Hn1+fn+1 (5)

    式中:fn+1力源项采用如下格式进行离散:

    fn+1=ufu*Δt (6)

    式中:uf为考虑了力源项的Euler网格点速度。

    首先在次区域内通过算法判断和标识流体节点、固体节点,确定受力点。在任意的流体网格点周围找到相邻的4个网格点,若相邻的这4个网格点中只要有1个为固体点,那么认为该点为考虑固体边界影响的流体网格点,即受力点,如图1(a)所示;然后,构造如图1(b)所示的“虚拟点—受力点—垂足点”计算结构:过受力点F向固体边界做垂线交固体边界于垂足点B,反向延长垂线到虚拟点I,本文使点I到点F的距离等于受力点F到垂足B的距离h。首先,在局部插值得到虚拟点I的速度,再结合垂足点B的速度利用线性插值确定受力点F的速度uf,其中垂足点满足固体边界无滑移条件。

    图  1  固体边界处理
    Figure  1.  Treatment of the solid boundary

    虚拟点I的速度利用其周围4个Euler网格点速度进行双线性插值得到:

    uI=4n=1ϕiun (7)

    式中:un为距离虚拟点最近的4个Euler网格上的点(3点不可在同一直线上),ϕi为常系数:

    ϕiT = [xIyIxIyI1]A1 (8)

    式中:

    A = [x1y1x1y11x2y2x2y21x3y3x3y31x4y4x4y41] (9)

    式中:下标1~4分别指距离虚拟点I最近的四个Euler网格点。

    (3)考虑新的压力项,更新流场,求得下一步的速度u,即:

    un+1uΔt=pn+1 (10)
    2p = 1Δtu (11)

    式中:pn+1为通过解耦压力Poisson方程(11)得到的下一步压力值。

    对于对流项,设ϕ为流场内的任意变量,对于流体区域内部点,采用如下格式进行离散:

    ϕϕx=ϕi,j2Δx(ϕi+1,jϕi1,j)ϕ+λΔx(ϕi+1,j3ϕi,j+3ϕi1,jϕi2,j)ϕλΔx(ϕi+2,j3ϕi+1,j+3ϕi,jϕi1,j) (12)

    式中:ϕ+=12(ϕi,j+|ϕi,j|)ϕ=12(ϕi,j|ϕi,j|)λ=0.125。在边界位置,对流项采用一阶迎风格式离散。

    对于扩散项,在边界上采用一阶迎风格式,中间节点处采用中心差分格式,即:

    2φx2=φi+1,j2φi,j+φi1,jΔx2,2φy2=φi,j+12φi,j+φi,j1Δy2 (13)

    在对压力的离散过程中,采用五点差分离散压力泊松方程:

    pi+1,j2pi,j+pi1,jΔx2+pi,j+12pi,j+pi,j1Δy2=x(uΔt)+y(vΔt) (14)

    对于Neumann边界条件,在边界点处采用增设虚点的方法[28]。对式(14)等号右边项,内部节点采用中心差分格式:

    x(uΔt)+y(vΔt)=1Δt(ui+1,jui1,j2Δx+vi,j+1vi,j12Δy) (15)

    边界节点采用一阶迎风格式,如在左边界上i=0,有:

    x(uΔt)+y(vΔt)=1Δt(u1,ju0,jΔx+v0,j+1v1,j12Δy) (16)

    图2所示,计算域为一矩形区域,长×宽=15D×10D,流体次区域为边长1.5D的正方形,次区域左侧距流场入口的距离为4.25D,次区域中央浸没一个刚性圆柱,其约束方式为全固定边界,圆柱直径为D,圆心坐标为(0,0)。计算域左侧为均匀来流入口边界,采用Dirichlet边界条件,即u=U,v=0;上、下两侧均为无穿透边界;右侧为自由出流边界,采用Neumenn边界,即u/x=0v/x=0。整个流场采用一套间距为Δx=Δy=0.025的均匀四边形网格,固体边界采取与流体网格尺度相等的等间距离散,时间步长为Δt=0.002

    图  2  整体计算域及边界条件
    Figure  2.  Computational domain and boundary conditions

    图3t=160 s、雷诺数Re=300时圆柱绕流的流场速度分布(范围:1u/U1)。从图中可看出,速度等值线表现光滑,说明计算过程中对速度的更新、修正并没有造成流场的震荡。由此可知,采用设置中间速度并不断更新速度的方法是适用的。同时,固体边界周围的速度分布也符合真实的圆柱绕流流场分布规律。说明利用这种锐利界面法描述固体对流体的耦合作用和运用双线性插值方法计算虚拟点的速度是可行的。

    图  3  流向速度分布
    Figure  3.  Isolines of velocity

    图4为一个涡脱落周期T内圆柱尾迹涡随时间的演化过程(范围:8wzD/U8wz为涡量;实线为正值,虚线为负值)。从中可看出柱体表面的边界层分布、分离剪切层的卷起和圆柱后方的分离区域清晰可见。在一个周期内圆柱尾部出现明显的正、负漩涡交替脱落,向下游发展,呈现出典型的卡门涡街现象,与试验得到的卡门涡街现象一致。说明本文对边界附近的流动参数插值的计算方法及程序是可行的。

    图  4  一个周期内尾迹涡的演化
    Figure  4.  Isolines of vorticity against time

    图5Re=300时尾涡脱落稳定后的升力系数及阻力系数的时程曲线,可以看出,随着时间的推进,流场逐渐趋于稳定,在流场稳定后升阻力系数都随时间的发展而呈现出稳定的周期摆动。其阻力系数、Strouhal数(St)的全时域统计结果见表1所示。其中升力系数cL = 2Fy/ρU2D,阻力系数cD = 2Fx/ρU2D,Strouhal数St = fsD/U,fs为涡的脱落频率。对比表中结果可知,相比于文献[16,29]的数值解,本文的计算结果与文献[30]的实验值更为接近。验证了本方法及程序的准确性。

    图  5  升力、阻力系数时程曲线
    Figure  5.  Variations of the lift and drag coefficients with time
    表  1  本文结果与其他文献结果的对比(Re=300
    Table  1.  The results comparison of average drag coefficients and Strouhal number at Re=300
    算例平均阻力系数¯cDSt
    文献[16]1.240.215
    文献[29]1.270.21
    文献[30]1.400.20
    本文1.360.208
    下载: 导出CSV 
    | 显示表格

    目前,以圆柱为对象的钝体绕流的研究已较为充分,对椭圆柱的关注度相对要少一些,但椭圆柱作为介于圆柱和平板之间更具有代表性的钝体,其形状更具流线型,流动阻力低,具有较好的对抗流体的作用,而且剪切边界层的卷起、尾迹涡的脱落、转捩和耗散等也很复杂,其研究对工程实际具有重要的现实意义。取计算域尺寸、边界条件、网格尺寸、时间步长等均与验证算例一致,雷诺数Re=800。不同之处在于:流体次区域中央浸没的固体为一系列不同轴长比AR=b/a的椭圆柱(见图6),其绕枢轴点(原点O)旋转摆动时攻角θ(摆角)的变化规律由下式控制:

    图  6  椭圆柱绕流计算模型
    Figure  6.  Computational model of flow around an elliptical cylinder
    θ={θ0sin(12βωt)0tπβωθ0πβωtπ ω(21β)θ0sin[12βωt+π(1β)]πω(21β)t2πω (17)

    式中:θ0为最大攻角值,取80°;ω=2πf为振荡角频率,f为运动频率, β为曲线形状参数。该椭圆柱的摆动形式由水翼在正弦振荡的基础上演变而来(当 β=1时为正弦波的半个波长)。 β=2.5代表非正弦振荡,前后0.2 s为匀速开、关时间段,中间持续0.6 s,其1个开、关周期T内的方波曲线如图7所示。

    图  7  一个周期内椭圆柱的瞬时攻角变化曲线
    Figure  7.  Changes of the instantaneous angle of attack of the elliptical cylinder within a period

    图8图10为攻角θ=0时,平均阻力系数¯cD、最大升力系数|cL|max和涡脱频率 {f_{{\rm{s}}}} 随椭圆柱形状轴长比的变化规律,以及流场流线情况。

    图  8  平均阻力系数、最大升力系数随轴长比的变化
    Figure  8.  Variation of lift and drag coefficients with axis ratio
    图  9  涡脱频率随轴长比的变化
    Figure  9.  Variation of vortex shedding frequency with axis ratio
    图  10  不同轴长比下的流线图
    Figure  10.  The instantaneous streamlines with different axis ratio

    图8可以看出,随着AR的不断增大,平均阻力系数 \overline {{c_{{\rm{D}}} }} 以及最大升力系数 {\left| {{c_{{\rm{L}}} }} \right|_{\max }} 均逐渐增大。由于椭圆柱垂直于水流方向的投影面积不断增大,当流体流过椭圆柱时,对椭圆柱的冲击力增大,因此椭圆柱对流体的反作用力也会增大,椭圆柱所受到来自流体的阻力 \overline {{c_{{\rm{D}}} }} 必然会增大。但是,相比阻力而言,升力的增大幅度比较小。由图9可以看出,涡脱频率 {f_{{\rm{s}}}} 随着轴长比AR的变化趋势可以分成3个阶段。初期: {f_{{\rm{s}}}} 随着轴长比AR的增大而降低,当AR=0.2时,椭圆柱的涡脱频率 {f_{{\rm{s}}}} 最大,约为0.6 Hz;第二阶段:当AR进入某一区域(AR=0.5~0.7)后, {f_{{\rm{s}}}} 开始稳定下来,约为0.38 Hz;后期:当AR超出这一区间后,涡脱频率 {f_{{\rm{s}}}} 将再次随着AR的增大而递减,“稳定现象”消失,在AR=0.9时,fs达到最小值。

    图10可知,当AR=0.2时,椭圆柱更具流线型,流动阻力低,尾部并没有明显的分离涡,流场比较平稳;当AR=0.5时,尾部流线波动开始逐渐明显;随着轴长比的增大,流态波动更加剧烈;当AR=0.7时,在后方开始产生一个尺度较小的分离涡,其位置向上偏离椭圆柱的中心。

    图11Re=800时不同轴长比工况下,同一时刻下的瞬时涡量图( - 8 {\text{≤}} {{{w_{\textit z}}D}/ U} {\text{≤}} 8 {w_{\textit z}} 为涡量;实线为正值,虚线为负值)。可以看出:当AR=0.2时,尾迹涡以S型模态脱落,呈现出上侧负漩涡与下侧正漩涡交替脱落向下游演化发展的卡门涡街现象。脱落的正负漩涡之间的纵向间距y0、周期T较小,形成明显的2列,尾涡强度(尺寸)也较小,漩涡的分离点基本位于x轴上;随着轴长比AR的增大,卡门涡街规律依然明显,但漩涡强度(尺寸)和周期T逐渐增大,漩涡分离点偏离x轴,且脱落位置x0增大,正、负漩涡中心间距y0减小逐渐趋于一列。

    图  11  不同轴长比工况下的瞬时涡量
    Figure  11.  Variation of instantaneous vorticity with axis ratio

    图12给出了AR=0.2,Re=800时,不同攻角 \theta 下,升、阻力系数的变化规律。

    图  12  不同攻角对升、阻力系数的影响
    Figure  12.  Variations of lift and drag coefficients with angle of attack

    图12可知,当 \theta {0^ \circ } {90^ \circ } 发展时,平均阻力系数 \overline {{c_{{\rm{D}}} }} 整体呈现为上升的趋势,这是由于随着攻角的逐渐增加,垂直于来流方向的正投影面积也逐渐增加,必然导致阻力的随之增大。而最大升力系数 {\left| {{c_{{\rm{L}}} }} \right|_{\max }} 的变化规律不同于平均阻力系数 \overline {{c_{{\rm{D}}} }} {\left| {{c_{{\rm{L}}} }} \right|_{\max }} 的发展趋势分为两个阶段,以 \theta = {45^ \circ } 为临界攻角,当 \theta {\text{<}}{45^ \circ } 时,最大升力系数 {\left| {{c_{{\rm{L}}} }} \right|_{\max }} 随着 \theta 的增大而增大,而当椭圆柱摆动超出临界攻角,即当 \theta {\text{>}} {45^ \circ } 时, {\left| {{c_{{\rm{L}}} }} \right|_{\max }} 便随着 \theta 的增大而逐渐减小。

    最大升力系数与平均阻力系数的比值 {\left| {{c_{{\rm{L}}} }} \right|_{\max }}{\text{/}}\overline {{c_{\rm{D}}}} 也是呈先增加后减小趋势,但升阻比的临界攻角则为 \theta = {\text{2}}{{\text{5}}^ \circ } 。当攻角 \theta 小于临界攻角时,曲线呈现出随 \theta 的增大而逐渐增长的趋势,此时 {\left| {{c_{{\rm{L}}} }} \right|_{\max }}{\text{/}}\overline {{c_{\rm{D}}}} {\text{<}} 1 {\left| {{c_{{\rm{L}}} }} \right|_{\max }} 数值与 \overline {{c_{{\rm{D}}} }} 的数值相差不大;当 \theta = {\text{2}}{{\text{5}}^ \circ } 时,升阻比达到最大,为1.01;当 \theta 超过临界攻角后,升阻比则呈现随 \theta 的增加而逐渐减小的趋势,这时升力在 \theta = {45^ \circ } 之前虽然依旧呈现增长的趋势,但是增长速率却远小于阻力的增长速率,而超过最大升力系数的临界攻角后,由于升力开始逐渐减小而阻力却仍然呈现增长的趋势,所以升阻比必然会逐渐减小。为了讨论产生这种现象的原因,分析了位于临界攻角附近的4个工况( \theta = {0}^{\circ },\theta = {\text{25}}^{\circ },\theta = {\text{45}}^{\circ },\theta = {\text{75}}^{\circ } )的流场瞬时压力场,如图13所示。

    图  13  不同攻角下的瞬时压力场
    Figure  13.  Instantaneous pressure fields at different angles of attack

    图13中也可看出,当 \theta = {{\text{0}}^ \circ } 时椭圆柱上下表面的压力对称,压差接近零,升力系数很小;随着攻角增大,椭圆柱的上表面的压力,不管从数值和承压面积上来看都明显比下表面大,即上下表面的压差增大,故升力系数随攻角的增大而增大。当攻角发展到临界攻角 \theta = {\text{4}}{{\text{5}}^ \circ } 时,压差达到最大,之后升力系数的值随着攻角的增大而逐渐减小。

    图14漩涡脱频率曲线可以看到,涡脱频率{f_{\rm s}}随着攻角 \theta 的不断增大,整体呈现先增加再逐渐递减的趋势。其变化趋势可以分成3个阶段:初期, \theta =2°~5°阶段,{f_{\rm s}}随着攻角 \theta 的增大而增加, \theta = {2^ \circ } 时,涡脱频率达到最大,约为{f_{\rm s}}=0.61 Hz;第二阶段,当攻角进入 \theta =2°~10°区域后,{f_{\rm s}}有所稳定,递减速率较小,在0.6 Hz左右浮动;后期,当攻角 \theta >10°后,涡脱频率{f_{\rm s}}将随着攻角 \theta 增大而逐渐减小,后期{f_{\rm s}}曲线逐渐平缓,减小的速率逐渐减慢。

    图  14  不同攻角对涡脱频率的影响
    Figure  14.  Variation of vortex shedding frequency with angle of attack

    图15Re=800时不同攻角工况下椭圆柱尾迹涡的对比,由图可见,固体表面分布的边界层、分离剪切层的卷起和椭圆柱后侧的分离区,以及涡的脱落等现象都清晰可见。与圆柱绕流一样,椭圆柱绕流尾迹涡也是呈现出明显的一对正负漩涡从上下两侧交替脱落,向下游发展的卡门涡街现象。但是,椭圆柱摆动的攻角幅值对流态起决定性作用。随着攻角的增加,涡的脱落周期、尺寸、形状和演化形式等均有较大差异,涡脱落的反对称性越发明显。类似于水轮机活动导叶大攻角运行或导弹无侧滑大攻角飞行时,背流(风)面分离涡的不对称运动。攻角越大越容易产生不希望的附加偏航力,或加大水轮机活动导叶转动轴折断的风险。

    图  15  不同攻角下的瞬时涡量图
    Figure  15.  Variation of the instantaneous vorticity with the angle of attack

    随着椭圆柱摆动攻角的增大,漩涡脱落频率降低,周期T加长,漩涡尺寸加大。当攻角 \theta {\text{≤}} {30^ \circ } 时,在每个周期内,有两个符号相反的涡从椭圆柱上下侧脱落,在下游形成交替出现的涡对,是典型的反对称S型脱落模态,如图16(a)所示,类似于静止圆柱的涡脱落方式。当 \theta =45°~60°时,上下侧脱落的漩涡尾部均拖着一个方向相同的子涡(2个涡核);在向下游演化发展的过程中,上侧脱落的负漩涡的子涡与母涡融合,而下侧脱落的正漩涡的子涡则产生分裂,分裂后的次正漩涡位于最上侧,靠近负漩涡,最终尾迹涡表现为3排,在下游形成了一对涡和一个单涡交替脱落的“P+S”Ⅰ型模态,如图16(b)所示。当 \theta {\text{≥}} {75^ \circ } 时,则情况相反,上侧脱落的负漩涡的子涡发生分裂,分裂后的次负漩涡位于最下侧,靠近正漩涡,而椭圆柱下侧脱落的正漩涡的子涡则与母涡则融合,在下游形成了新的“P+S”Ⅱ型模态,如图16(c)所示。

    图  16  不同攻角下的涡脱落模态示意图
    Figure  16.  The diagram of vortex modes at different angles of attack

    采用一种锐利界面(sharp-interface)浸入边界法对二维椭圆柱绕流问题进行了数值计算,并讨论了不同轴长比AR、攻角 \theta 对涡结构分布特征及升阻力系数、涡脱频率等水力不稳定现象的影响,通过对比分析,得到以下结论:

    (1)在模拟钝体绕流时,通过与文献和试验结果的对比,本文的锐利界面浸入边界法结果与试验结果吻合较好,验证了本方法的准确性和有效性;

    (2)对于椭圆柱绕流问题,随着不同轴长比的增加,升阻力系数呈递增趋势,而流场流态波动也越来越剧烈,尾迹涡脱落位置、强度(尺寸)、间距、周期也越来越大。周期加长,涡脱频率则降低;

    (3)当攻角发生变化时,平均阻力系数随攻角的增大而呈递增趋势,而最大升力系数,以及最大升力系数和平均阻力比(升阻比),均表现为前期随着攻角的增大而增大,超过临界攻角后,再降低的趋势;升阻比的临界攻角为25°,而最大升力系数的临界角为45°;

    (4)当攻角发生变化时,涡脱频率总体呈随着攻角增大而逐渐减小的趋势。漩涡脱落频率降低,周期加长,漩涡尺寸加大;当攻角 \theta {\text{≤}} 3{{\text{0}}^ \circ } 时,尾涡呈反对称S型脱落模态,当 {45^ \circ } {\text{≤}} \theta {\text{≤}} 6{{\text{0}}^ \circ } 时,尾涡呈反对称“P+S”Ⅰ型脱落模态,当 \theta {\text{≥}} {75^ \circ } 时,尾涡呈反对称“P+S”Ⅱ型脱落模态。

  • 图  1  固体边界处理

    Figure  1.  Treatment of the solid boundary

    图  2  整体计算域及边界条件

    Figure  2.  Computational domain and boundary conditions

    图  3  流向速度分布

    Figure  3.  Isolines of velocity

    图  4  一个周期内尾迹涡的演化

    Figure  4.  Isolines of vorticity against time

    图  5  升力、阻力系数时程曲线

    Figure  5.  Variations of the lift and drag coefficients with time

    图  6  椭圆柱绕流计算模型

    Figure  6.  Computational model of flow around an elliptical cylinder

    图  7  一个周期内椭圆柱的瞬时攻角变化曲线

    Figure  7.  Changes of the instantaneous angle of attack of the elliptical cylinder within a period

    图  8  平均阻力系数、最大升力系数随轴长比的变化

    Figure  8.  Variation of lift and drag coefficients with axis ratio

    图  9  涡脱频率随轴长比的变化

    Figure  9.  Variation of vortex shedding frequency with axis ratio

    图  10  不同轴长比下的流线图

    Figure  10.  The instantaneous streamlines with different axis ratio

    图  11  不同轴长比工况下的瞬时涡量

    Figure  11.  Variation of instantaneous vorticity with axis ratio

    图  12  不同攻角对升、阻力系数的影响

    Figure  12.  Variations of lift and drag coefficients with angle of attack

    图  13  不同攻角下的瞬时压力场

    Figure  13.  Instantaneous pressure fields at different angles of attack

    图  14  不同攻角对涡脱频率的影响

    Figure  14.  Variation of vortex shedding frequency with angle of attack

    图  15  不同攻角下的瞬时涡量图

    Figure  15.  Variation of the instantaneous vorticity with the angle of attack

    图  16  不同攻角下的涡脱落模态示意图

    Figure  16.  The diagram of vortex modes at different angles of attack

    表  1  本文结果与其他文献结果的对比(Re =300

    Table  1.   The results comparison of average drag coefficients and Strouhal number at Re=300

    算例平均阻力系数 \overline {{c_{{\rm{D}}} }} {S_{\rm{t}}}
    文献[16]1.240.215
    文献[29]1.270.21
    文献[30]1.400.20
    本文1.360.208
    下载: 导出CSV
  • [1] 李仁年, 赵振希, 李德顺, 李银然, 陈霞, 于佳鑫. 风沙对风力机翼型绕流及其气动性能的影响 [J]. 农业工程学报, 2018, 34(14): 205–211; 303. DOI: 10.11975/j.issn.1002-6819.2018.14.026.

    LI R N, ZHAO Z X, LI D S, et al. Effect of wind sand on flow around airfoil wind turbine and its aerodynamic performance [J]. Transactions of Chinese Society of Agricultural Engineering, 2018, 34(14): 205–211; 303. DOI: 10.11975/j.issn.1002-6819.2018.14.026.
    [2] 贲安庆, 窦华书. 可压缩机翼绕流的数值模拟及其稳定性分析 [J]. 浙江理工大学学报, 2015, 33(9): 675–681.

    BEN A Q, DOU H S. Numerical simulation of compressible flow around the airfoil and its stability analysis [J]. Jounal of Zhejiang Sci-Tech University, 2015, 33(9): 675–681.
    [3] 刘雄, 梁湿. 风力机翼型在复合运动下的动态失速数值分析 [J]. 工程力学, 2016, 33(12): 248–256. DOI: CNKI:SUN:GCLX.0.2016-12-030.

    LIU X, LIANG S. Numerical investigation on dynamic stall of wind turbine airfoil undergoing complex motion [J]. Engineering Mechanics, 2016, 33(12): 248–256. DOI: CNKI:SUN:GCLX.0.2016-12-030.
    [4] KOIRALA R, NEOPANE H P, ZHU B S, et al. Effect of sediment erosion on flow around guide vanes of Francis turbine [J]. Renewable Energy, 2019, 136: 1022–1027. DOI: 10.1016/j.renene.2019.01.045.
    [5] LI W Z, WANG W Q, YAN Y, et al. Effects of pitching motion profile on energy harvesting performance of a semi-active flapping foil using immersed boundary method [J]. Ocean Engineering, 2018, 163(9): 94–106. DOI: 10.1016/j.oceaneng.2018.05.055.
    [6] 郝栋伟, 张立翔, 王文全. 流固耦合S-型自主游动柔性鱼运动特性分析 [J]. 工程力学, 2015, 32(5): 13–18.

    HAO D W, ZHANG L X, WANG W Q. Swimming patterns of an S-type self-propelled flexible fish in fluid-structure interaction [J]. Engineering Mechanics, 2015, 32(5): 13–18.
    [7] 向锦武, 孙毅, 申童, 李道春. 扑翼空气动力学研究进展与应用 [J]. 工程力学, 2019, 36(4): 8–23.

    XIANG J W, SUN Y, SHEN T, et al. Research progress and application of flapping wing aerodynamics [J]. Engineering Mechanics, 2019, 36(4): 8–23.
    [8] BOWERS A. Model tests showed aerodynamic instability of Tacoma narrows bridge [J]. Journal of Franklin Institute, 1941, 231(5): 470–470.
    [9] 颜大椿. 湍流、风工程和虎门大桥的风振 [J]. 力学与实践, 2020, 42(4): 523–525.

    YAN D C. Turbulence, wind engineering and wind vibration of Humen Bridge [J]. Mechanics in Engineering, 2020, 42(4): 523–525.
    [10] PESKIN C S. Flow patterns around heart valves: a numerical method [J]. Journal of Computational Physics, 1972, 10: 252–271. DOI: 10.1016/0021-9991(72)90065-4.
    [11] GOLDSTEIN D, HANDLER R, SIROVICH L. Modeling a no-slip flow boundary with an external force field [J]. Journal of Computational Physics, 1993, 105(2): 354–366. DOI: 10.1006/jcph.1993.1081.
    [12] JAMALUDIN M J. Combined immersed boundaries/ B-spline methods for simulations of flow in complex geometries [J]. Annual Research Briefs, Center for Turbulence Research, 1997, 161(1): 317–327.
    [13] FADLUN E A, VERZICCO R, ORLANDI P, et al. Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations. [J]. Journal of Computational Physics, 2000, 161(1): 35–60. DOI: 10.1006/jcph.2000.6484.
    [14] LE D V, KHOO B C, LIM K M. An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains [J]. Computer Methods in Applied Mechanics and Engineering, 2008, 197(25): 2119–2130. DOI: 10.1016/j.cma.2007.08.008.
    [15] WU J, SHU C. Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications [J]. Journal of Computational Physics, 2009, 228(6): 1963–1979. DOI: 10.1016/j.jcp.2008.11.019.
    [16] 王文全, 张国威, 闫妍. 模拟复杂流动的一种隐式直接力浸入边界方法 [J]. 工程力学, 2017, 34(2): 28–33; 93. DOI: 10.6052/j.issn.1000-4750.2015.07.0600.

    WANG W Q, ZHANG G W, YAN Y. An implicit direct force immersed boundary method for simulating complex flow [J]. Engineering Mechanics, 2017, 34(2): 28–33; 93. DOI: 10.6052/j.issn.1000-4750.2015.07.0600.
    [17] Uhlmann M. An immersed boundary method with direct forcing for the simulation of particulate flows [J]. Journal of Computational Physics, 2005, 209(2): 448–476. DOI: 10.1016/j.jcp.2005.03.017.
    [18] LAI M C, PESKIN C S. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity [J]. Journal of Computational Physics, 2000, 160(2): 705–719. DOI: 10.1006/jcph.2000.6483.
    [19] 郭涛, 郝栋伟, 李明华, 等. 基于浸入边界法研究超弹性红细胞在剪切流中的运动特性 [J]. 医用生物力学, 2015, 30(3): 243–248. DOI: 10.3871/j.1004-7220.2015.03.243.

    GUO T, HAO D W, LI M H, et al. Motion characteristics on hyper-elastic red cells in shear flow based on immersed boundary meth [J]. Journal of Medical Biomechanics., 2015, 30(3): 243–248. DOI: 10.3871/j.1004-7220.2015.03.243.
    [20] SOTIROPOULOS F, YANG X. Immersed boundary methods for simulating fluid-structure interaction [J]. Progress in Aerospace Sciences, 2014, 65(5): 1–21. DOI: 10.1016/.j.paerosci.2013.09.003.DOI:.
    [21] MOHAMMADI M H, SOTIROPOULOS F, BRINKERHOFF J. Moving least squares reconstruction for sharp interface immersed boundary methods [J]. International Journal for Numerical Methods in Fluids, 2019, 90(20): 57–80. DOI: 10.1002/fld.4711.
    [22] 郭涛,张纹惠,王文全,等. 基于IBM法的低雷诺数下涡激振动高质量比效应的研究 [J]. 工程力学, 2022, 39(3): 222–232. DOI: 10.6052/j.issn.1000-4750.2021.07.0566.

    GUO T, ZHANG W H, WANG W Q, et al. Effects of high mass and damping ration on VIV of a circular cylinder with low Reynolds number based on IBM [J]. Engineering Mechanics, 2022, 39(3): 222–232. DOI: 10.6052/j.issn.1000-4750.2021.07.0566.
    [23] XIE F T, QU Y G, ISLAM M A, et al. A sharp-interface Cartesian grid method for time-domain acoustic scattering from complex geometries [J]. Computers and Fluids, 2020, 202: 104498. DOI: 10.1016/j.compfluid.2020.104498.
    [24] YOUSEFZADEH M, BATTIATO I. High order ghost-cell immersed boundary method for generalized boundary conditions [J]. International Journal of Heat and Mass Transfer, 2019, 137(7): 585–598. DOI: 10.1016/j.ijheatmasstransfer.2019.03.061.
    [25] BRADY P T, LIVESCU D. Foundations for high-order, conservative cut-cell methods: stable discretizations on degenerate meshes [J]. Journal of Computational Physics, 2021, 426: 109794. DOI: 10.1016/j.jcp.2020.109794.
    [26] MONASSE L, DARU V, MARIOTTI C, et al. A conservative coupling algorithm between a compressible flow and a rigid body using an embedded boundary method [J]. Journal of Computational Physics, 2012, 231(7): 2977–2994. DOI: 10.1016/j.jcp.2012.01.002.
    [27] 张和涛, 宁建国, 许香照, 等. 一种强耦合预估-校正浸入边界法 [J]. 爆炸与冲击, 2021, 41(9): 094201. DOI: 10.11883/bzycj-2021-0129.

    ZHANG H T, NING J G, XU X Z, et al. A strong coupling prediction-correction immersed boundary method [J]. Explosion and Shock Waves, 2021, 41(9): 094201. DOI: 10.11883/bzycj-2021-0129.
    [28] 胡建伟, 汤怀民. 微分方程数值方法[M]. 2 版. 北京: 科学出版社, 2007: 79–105.
    [29] BALARAS E. Modeling complex boundaries using an external force field on fixed Cartesian grids in large-eddy simulations [J]. Computers & Fluids, 2004, 33(3): 375–404. DOI: 10.1016/S0045-7930(03)00058-6.
    [30] SCHLICHING H. Boundary-layer theory [M]. New York: Mcgraw-Hill Book Company, 1979: 19–21.
  • 期刊类型引用(1)

    1. 郭涛,张纹惠,王文全,罗竹梅. 基于IBM法的低雷诺数下涡激振动高质量比效应的研究. 工程力学. 2022(03): 222-232 . 百度学术

    其他类型引用(2)

  • 加载中
图(16) / 表(1)
计量
  • 文章访问数:  422
  • HTML全文浏览量:  392
  • PDF下载量:  49
  • 被引次数: 3
出版历程
  • 收稿日期:  2021-08-16
  • 修回日期:  2022-08-25
  • 网络出版日期:  2022-08-25
  • 刊出日期:  2022-09-09

目录

/

返回文章
返回