前言：寻找写作灵感？中文期刊网用心挑选的波浪柱流致振动高保真数值模拟，希望能为您的阅读和创作带来灵感，欢迎大家阅读并分享。
摘要波浪柱具有良好的减阻减振效果,在海洋工程和桥梁工程中应用广泛。为了研究低雷诺数均匀流下细长柔性波浪柱的流致振动特性,选取了3种不同波高的截面直径呈正弦变化的波浪柱,分别给予其驻波和行波的初始激励,对其动力特性和尾流特征进行分析,并与同等条件下的光滑直圆柱进行了对比。研究结果表明,柔性波浪柱在驻波和行波初始激励下的流致振动特性有很大差别:在驻波激励下,当波高≥0.2D时,波浪柱振动得到了很好的控制,这主要归功于波浪柱的特殊结构防止了剪切层与强度较高的涡结构发生掺混,增大了漩涡生成长度;而在行波初始激励下,波浪柱的控制效果有限,局部雷诺数主导了结构响应、尾流模式和展向变化特征。
关键词低雷诺数;柔性波浪柱;流致振动;驻波;行波;动力特性;尾流特征
波浪柱由于具有良好的减阻和抑制涡激振动的效果,在实际工程中得到越来越广泛的应用。近年来,针对波浪柱的流致振动特性和减振机理,学者们展开了一系列的研究。Lam等[1-3]采用大涡模拟,对亚临界雷诺数下波浪柱波长、波幅、入流角对减阻效果的影响进行了分析,发现在合适的波长和波幅条件下,波浪柱在亚临界雷诺数下减阻效果依然优秀,同时波浪柱减阻效果与偏转角成反比;Garba-ruk等[4]采用准三维方法将波浪柱的全三维流动结构建模转化为一系列代表截面的二维流动特性问题,并对预测模型进行了稳定性分析,结果表明,稳定性分析结果与数值模拟结果吻合较好,为以后对波浪柱流场进行准三维模拟提供了可靠的模型框架;Zhang等[5]数值模拟研究了初始激励下Re=5000时波浪柱的流致振动特性;给圆柱施加了一个横向正弦形式的激励,并且固定了激励的振幅,结果表明,随着初始激励频率的增加,升力系数的分布发生显著变化,但初始激励与波浪柱及其流场之间的相互作用机理还待进一步探明;安苗等[6]利用风洞试验,对波浪柱在高雷诺数下的气动性能进行了分析,结果表明在高雷诺数下波浪柱的减阻效果有限;赵桂欣等[7]使用了大涡模拟对不同波长波幅组合的有限长波浪柱进行了研究,发现在亚临界雷诺数下大部分波长波幅组合的波浪柱具有明显减阻效果。目前研究大都以固定或者弹性支承的刚性波浪柱为对象,然而在实际工程中,柔性波浪柱具有更广阔的应用前景。本次研究在Re=100的均匀来流下,以截面直径呈正弦规律变化的柔性波浪柱为研究对象,对其在行波和驻波初始激励下的尾流特性、动力响应以及背后的耦合效应机制进行了探究。
1计算模型的建立及网格划分
1.1数值计算模型
均匀流中的波浪柱模型的示意图如图1所示。该圆柱直径沿展向变化由下式确定。其中:Dz为展向相应位置的圆截面直径;D代表波浪柱的平均直径,由于波浪柱直径在展向呈正弦变化,平均直径D=(Dmax+Dmin)/2,定义波浪柱直径最大截面所在的展向位置为几何节点,波浪柱直径最小截面所在的展向位置为几何鞍点;A代表波状表面的波高,共设置了3种波高,分别为0.1D、0.2D、0.3D,展向波长设置为λ=4πD。为了方便对比,在相同来流条件下对直径为D的光滑圆柱也进行了数值模拟。此处需要说明的是:以下讨论中所说的圆柱,只要没有具体指出是波浪柱还是光滑圆柱,则两种类型圆柱都包含在内。
1.2数值计算方法
采用Serson等[8]详细介绍的坐标映射方法并结合式(1)来构建波浪柱。设定均匀流的来流速度为U∞,与之相应的雷诺数定义为Re=U∞D/ν(ν为流体的运动粘度系数),将雷诺数设定为100,使流动保持层流状态。展向的计算域设置为4πD,与波浪柱一个波长的大小相同。参考Prasanth等[9]的设置,将圆柱的质量比设置为10,并给出圆柱质量比的计算公式m*=4ρc/ρfπD2(ρc和ρf分别是单位长度的圆柱质量和流体密度)。限制圆柱在顺流向的运动,只考虑圆柱在横流向的运动,因为顺流向运动的振幅通常仅为横流向振幅的10%~15%[10]。采用张力梁模型来表征圆柱的动力学响应,张力梁模型的控制方程为ρc∂2ξy∂t2+b∂ξy∂t=T∂2ξy∂z2+F(z,t),(2)其中:b为圆柱的结构阻尼系数,设置为0,允许产生最大的横向振幅振荡;T、ξy分别表示作用于两端的拉力和横向振幅。由于我们只对锁定状态感兴趣,结构固定频率被设置为2πSt,以保证其与漩涡脱落频率相近,其中St是在相同雷诺数下均匀流流过光滑圆柱的斯特罗哈尔数(Strouhalnumber)。然后,根据Newman等[11]计算出相应的圆柱张力T=m*L2S2t,F(z,t)是沿展向作用在圆柱上的力,由压力和粘滞应力项的积分计算。三维流动由开源软件Nektar++中的不可压缩Navier-Stokes求解器控制,该软件基于傅立叶谱元法[12]。其中:u是速度;p是压强。Nektar++中采用的流固耦合方法是(x,y)平面上的Jacobi-Galerkin公式和z方向上的傅立叶展开的混合格式。将(x,y)平面的二维空间域离散为四边形有限元单元,在每个四边形有限元单元内,再采用高斯-洛巴托-勒让德(GLL)插值点关联的高阶拉格朗日多项式作为形状和权重函数进行离散。关于空间离散化算法的细节可参考文献[12-13]。时间离散化基于高阶刚性稳定分裂格式,允许原始变量在每一个时间步长内独立处理。类似应用于流体流动控制方程的傅立叶变换,对结构的运动施加了一个在全长范围内周期性的假设。然后将运动变量和力表示为复傅立叶级数,将拉伸梁模型的偏微分方程解耦为一组常微分方程,用二阶Newmark-β方法求解。
1.3网格设置
网格设置如图2所示。图2(a)为二维平面基础网格的示意图,圆柱被放置在一个C形计算域中,C形计算域由圆柱上游半径为20D的半圆形计算域和圆柱下游的矩形计算域组成,矩形计算域在流向长度为30D,在横流向长度为40D。在(x,y)平面上采用了877个基础网格,同时使用P=6对网格进行差分,如图2(b)所示,每个基础网格被分成了(6-1)×(6-1)个高阶网格。此外,根据Newman等[11]和Zhu等[14]的研究,将计算域沿展向划分了64个傅立叶面。计算域的边界条件设置如下:入口边界为均匀来流条件,来流速度为(U∞,0,0);在出口边界上采用纽曼边界条件,令出口边界上以及p=0;在顶部和底部边界上指定(u,v,w)=(U∞,0,0)的远场边界条件;在圆柱固壁采用无滑移和不可渗透边界条件;通过假设无限长波浪柱的展向周期性,在波浪柱的两端实施了只允许横向运动的自由边界条件[11]。此外,Lam等[1]证明了只有单个展向波长的周期边界条件对于计算波浪柱绕流足够精确。
1.4网格验证
为了确保模拟使用的网格具有足够的精度,并且保证在现有的网格设置下继续增大网格密度对模拟结果影响较小,需要对测试算例进行网格独立性验证。采用在(x,y)平面改变多项式的阶数P来改变平面网格的疏密,在展向上通过改变傅立叶面数目来改变z方向网格的疏密。不同网格数目下的测试算例结果对比见表1。测试算例均是在雷诺数Re=100条件均匀来流下的光滑圆柱,所有测试算例光滑圆柱的质量比m*=10。由表1可得,在不同的网格数目下,模拟结果的差异非常小。为了兼顾模拟的精度与计算速度,在主要模拟中使用的网格精度设置为P6Nz64。
2驻波初始激励下模拟结果分析
2.1动力响应与振幅
横向振幅(ξy)、阻力系数(Cd)和升力系数(Cl)沿展向的均方根值分布如图3所示。阻力系数Cd=2Fd/ρfU2∞D,升力系数Cl=2Fl/ρfU2∞D,其中Fd和Fl分别代表圆柱受到的流向和横向的截面力的分量。为了能够更加精确地捕捉到升阻力和振幅随波浪柱波高变化的规律,在此处添加了A=0.15D和A=0.25D两种工况。由图3(a)可知,随着波浪柱波高A的增加,波浪柱的横向振幅逐渐减小,尤其是在A≥0.2D的工况下,波浪柱的振动几乎完全被抑制。图3(b)中的阻力系数均方根分布也显示出同样的规律,并且在A≥0.2D的3种工况下,阻力系数相比于光滑圆柱有了明显降低。而图3(c)中升力系数的均方根分布显示,在A<0.2D的两种工况下,升力系数相比于光滑圆柱反而有所升高,继续增加波浪柱波高至A≥0.2D后,升力系数才会显著下降。为了进一步阐明波浪柱波高与振动抑制之间的关系,对横向振幅及升阻力系数的时空演变规律进行了分析,此处取光滑圆柱、A=0.1D和A=0.2D3种代表工况进行分析。光滑圆柱和波浪柱横向振幅(ξy)的时空演变云图见图4。由图4可以看到,A=0.1D工况的节点和反节点位置与光滑圆柱相同,节点位于z/D=0、2π、4π,反节点位于z/D=π、3π。同时在任何无量纲时间(tU∞/D)下,A=0.2D工况下波浪柱的横向振幅均较光滑圆柱大幅降低,而A=0.1D工况下波浪柱抑制振动的效果较差。阻力系数(Cd)的时空演变云图见图5。由图5可以看到,波高A=0.1D工况下阻力系数的时空分布和大小与光滑圆柱的基本相同;而A=0.2D工况下的阻力系数时空分布则与光滑圆柱完全不同,在阻力系数分布中看不到明显的驻波效应,在z/D=0、4π时阻力系数并未取得最大值,并且在不同无量纲时间处阻力系数最小的点是沿展向变化的,同时A=0.2D波浪柱的最大阻力系数仅有光滑圆柱的一半,说明A=0.2D波浪柱具有良好的减阻效果。波浪柱升力系数(Cl)的时空演变云图见图6。由图6可以看到,A=0.1D工况下波浪柱在展向π<z/D<3π范围内的升力系数较光滑圆柱有所下降,而在其他展向位置升力系数较光滑圆柱有所升高,这是由于局部雷诺数的增加引起的。而A=0.2D的波浪柱仅在展向的两侧有较高的升力系数,其他位置升力系数非常小。根据波浪柱在驻波初始激励下对横向振幅和升阻力的抑制效果,可以将波浪柱分成控制有效(A=0.2D、A=0.3D)和控制失效(A=0.1D)波浪柱两类。
2.2尾流特性
为了进一步了解不同波高导致波浪柱不同动态响应的内在机制,需对波浪柱和光滑圆柱近尾迹的三维涡模式展开研究。3种代表工况下展向涡量瞬时等值面的透视图和横流向俯视图见图7。从图7可以看出,驻波初始激励下光滑圆柱的展向涡量中存在明显的交织结构,这与Newman等[11]所观察到的相同。在A=0.1D工况下也可以大概看到交织的涡结构,并且在几何节点(波浪柱截面最大展向位置,处在波浪柱的两端)的展向涡强度较大,而几何鞍点(波浪柱截面最小展向位置,处在波浪柱的跨中)处的展向涡强度较小,强度较强的涡会在两端产生较高的截面升力系数,也解释了图6(b)的现象。在A=0.2D工况下交织涡结构完全消失,弯曲涡管以交错的方式在波浪柱上下表面分别脱离,并且在尾流中迅速消失。更重要的是,与光滑圆柱相比,A=0.2D波浪柱剪切层的稳定性较差,并且向上卷动较弱,形成涡旋。涡旋在下游进一步分离,使得涡形成区长度进一步增大。Lin等[15]用数值方法研究了亚临界流动区域中具有相对较大展向波长的刚性波浪圆柱周围的流动,也发现了类似的结果。
3行波初始激励下模拟结果分析
Newman等[11]的研究表明,对于无量纲长度超过4π的光滑直圆柱而言,结构的振动模式最终发展为行波响应。然而,对于长度为4π的光滑直圆柱,经过行波-驻波响应的竞争,最终结构的振动模式发展为驻波形式,这表明了结构的长度是影响结构动力响应模式的关键因素之一。对波浪柱在行波初始激励下的动力响应特性和尾流特征展开研究,赋予结构的行波初始激励为其中:a是振幅;ω是振动频率,ω=2πcLz,c=T/ρ表示相速度。
3.1动力响应与振幅
行波初始扰动下光滑圆柱的横向振幅和升阻力系数云图见图8。由图8(a)所示的光滑圆柱横向振幅可以看出:随着时间的演化,结构响应呈现出明显的行波响应特征,即响应由一端传递至另外一端,任意处的结构响应最大幅值一致且同一时刻结构展向不同位置处存在连续相位差,这与Newman等[11]的模拟结果相同,进一步证明了本模拟的准确性。行波初始扰动下3类典型工况结构的横向振幅时空演化云图见图9。总体而言,不同工况下结构的最大振幅响应基本一致,相较于光滑圆柱的横向振幅均有所增加。当A=0.1D时,同样出现了明显的行波响应,但结构两端(z/D=0,4π)和中间(z/D=2π)出现了相对较大振幅响应。随着波高A的继续增加,较大振幅响应位置逐渐集中至z/D=π,3π两位置处,如图9(b)和图9(c)所示。该现象呈现出了驻波-行波两类响应叠加的特征。行波初始扰动下波浪柱的阻/升力系数时空分布见图10。由图10可以看到,A=0.1D工况下作用于结构上的流体动力同样不再是单一的行波响应,而是呈现了明显的驻波-行波叠加响应。New-man等[11]发现柔性管道在无量纲时间为500~570时结构响应由驻波演化为行波,Facchinetti等[16]的研究同样表明,行波响应是该类结构的稳定状态。对于本次研究而言,初始激励为行波激励,但行波响应并不是稳定状态,这与结构的波浪形有关,即结构为一种固定形式(可当作驻波),与行波叠加,从而呈现这一特殊现象。此外,A=0.1D工况下阻力和升力幅值与光滑圆柱相比均有所提升。A=0.2D工况除了呈现与A=0.1D工况一致的驻波-行波叠加响应之外,结构两端(z/D=0,4π)阻升力明显大于展向其他位置。这与结构的形式有关,在该工况下结构两端的直径为1.4D,跨中(z/D=2π)为0.6D,根据局部雷诺数(Reloc=ρU∞dloc/μ)计算,升阻力与光滑圆柱相比产生了明显变化(最大阻升力提升40%)。A=0.3D工况除了呈现驻波-行波叠加响应外,流体动力响应的最大特征是失去了对称性,而该非对称特征随时间演化并未消失,而是成为稳定状态。该现象的出现是由于研究中结构本身为超柔性体,结构自激动力响应呈现强非线性,与展向变化的结构形式以及流体的耦合相互作用,使得整个结构系统呈现出强烈的非线性特征,从而动力响应呈现了分叉特征,当结构响应偏向一侧时,响应会朝着该方向继续发展。此外,由于波浪柱的各工况均呈现出驻波-行波叠加响应,以工况A=0.2D为例做出解释。A=0.2D工况无量纲时间tU∞/D=0~300的结构横向振幅时空演化图见图11。由图11(a)可见,在初始行波激励下,结构呈现行波特征,在tU∞/D=40时行波特征消失,随后在tU∞/D=90时再次出现行波响应;由图11(b)可见,随着时间的演化振幅不断增大(超过0.6D),且行波响应特征增强;由图11(c)可见,至tU∞/D=240~260,变为驻波-行波响应叠加稳定状态。结合图9(b),驻波-行波响应竞争结束后,结构横向振幅稳定至0.54D。因此,认为结构横流向呈现的驻波-行波响应是两种振动响应模式相互竞争、相互结合的结果。该现象通过涡结构开展进一步讨论。
3.2尾流特性
行波初始扰动下4类典型工况展向涡量等值面图(ω=±0.5)见图12,图中黄、蓝色分别表示展向正负值涡量。由图12(a)可以看出,结构展向呈现了规则的行波结构,这与前述以及Newman等[11]的研究结果一致,漩涡从结构两侧交替脱落,与管道结构呈现固定角度。此外,随着波浪柱波高A值的增加,驻波响应特征出现。由图12(b)可见,工况A=0.1D时,漩涡不再呈固定角度从管道结构上脱落,位于管道两端的涡尺度相对更大,至A=0.2D时(见图12(c)),两端涡强度显著大于中间部位。对于工况A=0.3D而言,行波或驻波特征通过涡量等值图难以识别,结合动力响应分析结构认为是驻波-行波叠加响应的结果。为更加清楚地表明圆柱不同展向位置后侧尾流形式的变化规律,提取了各工况的两个展向位置(z/D=π,2π)的展向涡量云图见图13。对于光滑圆柱而言(见图13(a)和图13(b)),两展向位置仅有相位差,尾流形式一致;当A=0.1D时(见图13(c)和图13(d)),相对于z/D=π位置处的尾流,z/D=2π处的尾流具有更宽的横流向涡核间距,该现象意味着位于z/D=2π的结构振幅更大,这与图8(b)的结果一致;当A=0.2D时(见图13(e)和图13(f)),z/D=π位置处的横流向涡核间距相对更大,认为这与显著增强的驻波响应有关;对于工况A=0.3D而言,z/D=π位置处的尾流呈现了转捩特征,而z/D=2π处的尾流出现了流动不分离的特征,认为这与局部雷诺数有关,在该工况下结构节点处的局部雷诺数高达160(对于圆柱绕流而言,层流转捩发生在雷诺数180左右),因此出现了早期转捩现象。而该工况时最小局部雷诺数为60(对于圆柱绕流而言,流动失稳发生在雷诺数50左右),因此出现了位于z/D=2π处下游区的稳定流动特征。
4结论
在驻波初始激励下合理地改变波高,可有效减小波浪柱的振动幅值和作用在其上的阻力。根据波浪柱在驻波初始激励下对横向振幅和升阻力的抑制效果,将波浪柱分成控制有效(A=0.2D、A=0.3D)和控制失效(A=0.1D)波浪柱两类。对于控制失效波浪柱,其动力特性与光滑圆柱基本一致,而对于控制有效波浪柱,其动力响应和受到的流体动力均显著降低。对于驻波响应具有良好控制效果的波浪柔性柱并未对行波响应起到良好控制效果,即截面尺寸沿展向正弦变化的波浪形圆柱对行波响应的控制效果有限,在行波响应下局部雷诺数主导了结构响应、尾流模式的展向变化特征;波浪形细长柔性柱在行波初始扰动下的动力响应与波高显著相关,呈现两类不同振动形式:当A≲0.1D时,为典型的行波振动模式;当A>0.1D时,经历驻波、行波振动竞争后,呈现出了驻波-行波叠加振动模式。
参考文献:
[1]LamK,LinYF.LareEdgdySimulationofFowlAroundWavyCylindersataSbcriticalRueynoldsNumber[J].InternationalJur-onalofHatandFuidFowell,2008,29(4):1071-1088.
[2]LamK,LinYF.EffectsofWavelengthandAmplitudeofaWavyCylinderinCoss-flowatLrowReynoldsNumbers[J].JournalofFuidMlechanics,2009,620:195-220.
[3]LamK,LinYF,ZouL,etal.InvestigationofTurbulentFowlPastaYawedWavyCylinder[J].JournalofFuidsandSruc-lttures,2010,26(7-8):1078-1097.
[4]GarbarukA,CrouchJD.Quasi-threeDmensionalAinalysisofGlobalIstabilitiesn:OnsetofVortexSeddinhgBehindaWavyCylinder[J].JournalofFuidMlechanics,2011,677:572-588.
[5]ZhangK,KatsuchiH,ZhouD,etal.NumericalSudtyofFowlPastaTansverselryOscillatingWavyCylinderatRe=5000[J].OceanEngineering,2018,169:539-550.
[6]安苗,刘庆宽,孙一飞,等.波浪形圆柱的气动力特性试验研究[J].振动与冲击,2021,40(8):209-215.
[7]赵桂欣,桂洪斌,王晓聪.有限长波浪形圆柱绕流数值模拟[J].哈尔滨工业大学学报,2021,53(6):163-170.
[8]SersonD,MeneghiniJR,SherwinSJ.Velocity-correctionS-chemesfortheIcomnpressibleNavier-StokesEquationsinGeneralCordinateoSystems[J].JournalofComputationalPhysics,2016,316:243-254.
[9]PrasanthTK,MittalSVortex-inducedVbrationsofaCrcu-.iilarClinderatLyowReynoldsNumbers[J].JournalofFuidlMechanics,2008,594:463-491.
[10]EvangelinosC,KarniadakisGE.DynamicsandFowlStruc-turesintheTurbulentWakeofRigidandFexibleClinderslySubjecttoVortex-InducedVbrationsi[J].JournalofFuidMle-chanics,1999,400:91-124.
[11]NewmanDJ,KarniadakisGE.ADirectNumericalSmulationiStudyofFowlPastaFeelryVibratingCable[J].JournalofFluidMechanics,1997,344:95-136.
[12]KarniadakisG,SherwinSS.pectral/hpElementMethodsforComputationalFuidDlynamics[M].OxfordUiversitnyPress,2013.
[13]KarniadakisGE,IsraeliM,OrszagSA.High-orderSlittinpgMethodsfortheIcomnpressibleNavier-StokesEquations[J].JournalofComputationalPhysics,1991,97(2):414-443.
[14]ZhuHongbo,WangR,BaoY,etal.FlowOveraSymmetri-callyCurvedCrcularClinderwiththeFeeSreamiyrtParalleltoThePaneofCrvatureatLluowReynoldsNumber[J].Jour-nalofFuidsandSructureslt,2019,87:23-38.
[15]LinYF,BaiHL,AlamMM,etal.EffectsofLrageSanwisepWavelengthontheWakeofaSnusoidalWiavyCylinder[J].JournalofFuidsandSructureslt,2016,61:392-409.
[16]FacchinettiML,deLnagreE,BiolleyF.Vortex-inducedTav-rellingWavesAonlgaCblea[J].EuropeanJurnalofMoechan-ics-B/Fluids,2003,23(1):199-208.
作者:庞玺源 朱宏博 张春晓 包艳 单位:上海交通大学船舶海洋与建筑工程学院 内蒙古路桥集团有限责任公司