低雷诺数变截面细长柔性结构流固耦合能量交换特征分析

邓秀兵, 于曰旻, 庞玺源

(1. 浙江省水利水电勘测设计院有限责任公司,杭州 310000; 2. 海南大学 土木建筑工程学院, 海口 570228; 3. 上海交通大学 船舶海洋与建筑工程学院,上海 200240)

在土木工程领域,风、浪、流等环境动力荷载作用下的细长结构物受损甚至破坏屡见不鲜.比如,悬索桥中的主缆索作为典型的细长柔性结构,在来风作用下将产生涡激振动,甚至发生疲劳损坏;又如,在主缆索施工过程中,缆索截面形式发生非对称性变化,当风速超过临界值后,将产生空气动力负阻尼,使得缆索振动逐渐增强,甚至超负荷振幅而破坏.因此,科学评估和有效抑制细长结构的流致动力响应是风、浪、流敏感细长柔性结构设计必须考虑的内容.然而,当前细长柔性结构的流固耦合效应方面研究不足,尤其精准模拟和分析细长结构流致振动响应的力学模型和算法方面仍较缺乏,且亟待深入理解其流固耦合机理并建立系统的分析方法与手段[1].

波浪形变截面圆柱体结构具有良好的流动减阻效果.近年来,国内外学者对该类型结构的绕流特性开展了较多研究,发现波浪形变截面能够延后剪切层的相互作用,从而有效减小阻力.Ahmed等[2]采用实验方法研究了波浪形圆柱的边界层分离线和尾流结构,发现波浪形圆柱流动分离线呈现明显的三维特性,在分离节点附近形成流向涡,且发生边界层“上卷”现象,从而延迟或抑制剪切层中湍流的生成和发展.Lam等[3]利用多种流动显示技术对波浪形圆柱的近尾迹进行了实验研究,给出了平均速度和波动速度分量沿流向、展向和横向的分布特征.实验结果表明,波浪形柱涡旋平均形成长度比大于光滑圆柱.湍流统计分析也表明,光滑圆柱尾迹中的涡街更为规则,而由于波浪形圆柱后面的涡具有较强的三维效应,使得由湍流掺混增强而尾迹表现出更强的非相干流动结构.

目前变截面结构绕流分析主要集中在静止或弹性支承的刚性波浪形结构的减阻和涡激振动特性方面.然而在实际土木工程中,缆索等细长结构的刚度系数通常较小,因此这种结构的柔性作用对流动减阻作用以及涡激振动响应的影响需要进一步研究,包括对初始动力激励的敏感性方面也需要深入分析.在文献[4-7]的基础上,针对这一问题,本研究采用基于高精度谱单元方法的流固耦合分析方法,建立了细长变截面柔性圆柱体结构的涡激振动力学模型,对其在低速均匀流作用和驻波初始动力激励下的流致振动机理进行了量化分析,获得了包括波浪形结构的尾流特性、结构动力响应特性、能量传递规律、涡脱频率展向变化特征,并对其在驻波扰动下的减阻、减振机理进行了深入探讨,对波浪形变截面细长结构的工程设计和应用提供参考.

1 数值算法

细长结构的流致振动是一种典型的流固耦合问题.求解流固耦合的基本方法通常有两大类,分别为界面匹配法和非界面匹配法,包括经典的任意拉格朗日-欧拉(ALE)方法和浸入边界法(IBM).本文采用另一种由Dimas等[8]提出的随体坐标系(Body-fitted Coordinates)方法,该方法在惯性坐标系中求解N-S方程,再通过坐标转换至非惯性坐标系中,从而不需要动网格和浸泡边界近似.

N-S方程和连续性方程在惯性坐标系(x′,y′,z′)内可描述为

(1)

(2)

由惯性坐标系至非惯性坐标系的坐标转换采用如下关系式:

x=x′-ζx(z,t),y=y′-ζy(z,t)

(3)

式中:ζx(z,t)和ζy(z,t)分别表示结构在顺流向和横流向的位移.相应地,速度项和压力项采用以下变换:

(4)

因此,将式(3)和(4)代入N-S方程式(1)和(2)可得:

A(u,p,ζ)

(5)

(6)

(7)

在此定义:

(8)

此外,本文采用线性张力梁模型描述细长柔性结构的动力学行为,该模型采用小变形假设,可写成:

(9)

式中:ρc,k和T分别表示结构的单位长度的质量、阻尼比和张力.需要指出的是,张力T的大小将影响相速度c,c=(T/ρc)1/2;ζ(z,t)=(ζx,ζy)表示结构在顺流向和横流向的位移;F(z,t)是作用于结构上的流体力,通过对压力和黏性力项沿结构固壁表面积分获得,

(10)

n为指向结构外向的法线单位向量;s为结构表面微分.此外,假设结构的动力响应满足沿展向的周期性条件,则对流场和结构变量可采用傅里叶级数(Fourier Expansion)表示,即

(11)

(12)

式中:β=2π/Lz表示展向的波数,Lz为结构长度;M为展开式中傅里叶模态个数;m为各阶模态.将式(11)和(12)代入式(5)和(6)可获得解耦后的二维模态方程:

(13)

(14)

(15)

对结构运动的解耦方程式(15),采用二阶Newmark-β方法进行求解.对不定常流场模态方程式(13)和(14),采用Karniadakis等[9]提出的高阶分步格式进行时间离散.具体地,对每一时间迭代步对速度和压力进行解耦计算.

(16)

式中:αq和βq均为与强稳定化积分(Stiffly Stable Integration)有关的参数.

第2步,将考虑压力梯度的作用修正速度场,并施加连续性约束条件和纽曼边界条件:

(17)

第3步,考虑黏性项更新下一步的速度场:

(18)

式中:γ0为强稳定计算过程中向后差分系数(Backwards Differentiation Coefficient).

对上述时间离散后的模态方程,根据Karniadakis等[9-10]将(x,y)平面的二维空间域离散为四边形有限元单元,并采用高斯-洛巴托-勒让德(GLL)高阶正交拉格朗日多项式作为形函数进行空间离散.

2 数值模型

2.1 物理模型

考虑在均匀来流作用下的波浪形变截面圆柱(见图1),其直径沿展向变化由下式确定:

图1 均匀流下截面直径沿展向余弦变化的波浪柱几何示意图

(19)

式中:Dz为展向相应位置的圆截面直径.由于波浪形柱直径在展向呈余弦变化,平均直径为D=(Dmax+Dmin)/2.基于平均直径和均匀来流速度的雷诺数取值为100,因此本文考虑的流动为层流流动.定义波浪形柱直径最大截面所在的展向位置为几何节点,波浪柱直径最小截面所在的对应位置为几何鞍点.A代表波状表面的波高,取值范围为0.1D～0.3D;展向波长设置为λ=4πD.A=0对应于光滑圆柱,也作为基准工况用于与波浪形柱工况进行对比分析.参照 Newman等[11]的研究,可以通过规定圆柱的初始振幅和速度来确定其初始条件,本研究中初始扰动为驻波:ζy(z,t)=acos(ωt)cos(2πz/Lz),其中a为振幅,ω为振动频率,ω=2π/(cLz).

二维平面网格划分如图2所示.计算域几何形状为C形,由圆柱上游半径为20D的半圆弧和下游矩形组成,其流向长度为30D,在横流向长度设置为40D,圆柱放置在计算域横流向的中心位置.在(x,y)平面内采用了877个基础网格;由于对每一基础单元再用6阶正交多项式(P=6)对速度和压力场进行近似,所以每个基础网格再细分为5×5个高阶网格.此外,根据文献[11-12],沿展向采用64个傅里叶模态,满足本文雷诺数下的精度要求.

图2 二维网格拓扑结构

2.2 模型验证

为了验证力学模型和数值模拟方法的适用性,首先获得了均匀来流作用下细长光滑圆柱涡激振动模拟结果,并与已有文献结果进行对比分析.本模拟中质量比设置为m*=6.0,而其他模拟参数均与Newman等[11]相同.图3所示为光滑圆柱在驻波和行波初始激励下的横向振幅(ξy)、阻力系数(Cd)和升力系数(Cl)模拟结果,与Newman等[11]结果吻合良好.尤其在驻波初始激励下本模拟结果与Newman等[11]结果在展向节点(响应振幅最大的点)和反节点(响应振幅最小的点)位置坐标完全相同.在行波初始激励下的行波响应幅值和行波速度等均与Newman等[11]结果一致,故说明本文力学模型和模拟方法的可靠性.

图3 柔性光滑圆柱涡激振动拟结果

接着正交多项式阶数对计算结果的影响进行验证分析.通过变化多项式阶数P改变高阶网格的疏密;在展向上通过改变傅里叶模态个数,改变展向网格的疏密.由计算结果表明:P=6,M=64时模拟结果足够收敛至精确解.

3 结果与讨论

3.1 运动响应与水动力系数

图4 ξy、Cd和Cl沿展向的均方根值分布

图5为在A=0.10D,0.20D时的结构横向振动响应时程模拟结果.由图可知,两种工况下节点和反节点的展向位置与光滑截面结构相同,即节点位于z/D=0, 2π, 4π,反节点位于z/D=π, 3π处,表明该工况下结构振动模式均与光滑表面结构相同.阻力系数在A=0.10D时其随时间演化特征与光滑表面结构基本相同,均在节点处取得最大值,在反节点处取最小值.然而A=0.20D工况下的阻力系数时空分布则与光滑表面结构完全不同,在阻力系数分布中看不到明显的驻波效应.在z/D=0, 4π时阻力系数并未取得最大值,且最大阻力系数仅有光滑表面结构的50%,表明A=0.20D扰动波高具有良好的减阻效果.进一步观察升力系数时程可知,A=0.10D工况下波形柱在展向π

图5 结构响应与水动力系数计算结果

3.2 尾流特性

为了进一步了解扰动表面波高对结构振动响应的影响机理,需要对柔性波形柱和柔性光滑柱近尾迹三维涡结构进行对比分析.图6给出了不同扰动波高时展向涡量等值面(ωz=±1.0)计算结果.其中,图6(a)～6(c)为柔性光滑柱展向涡量等值面透视图,图6(d)～6(f)为俯视图,与Newman等[11]计算结果吻合良好.驻波初始激励下柔性光滑柱展向涡量中存在明显的交织结构.这种交织涡结构的展向特征是,在z/D=0, 2π截面处形成交错脱落的卡门涡街结构;而在z/D=π, 3π截面处上、下表面同时脱落的完全对称型涡结构.图6中分别给出A=0.10D时模拟结果,显然在z/D=π, 3π截面处的涡结构也呈现对称分布.值得注意的是,在此工况下几何节点处的展向涡量幅值较大从而升力系数幅值也随之增大,而几何鞍点处涡量幅值较小从而升力系数幅值也随之减小,这与图5所示结果一致.图6中分别给出A=0.20D时模拟结果,在该工况下交织涡结构完全消失,弯曲涡管以交错方式从波浪形表面上、下侧分别脱离,并且在尾流中迅速消失.更重要的是,与柔性光滑柱相比,此时波形柱表面剪切层的卷曲和相互作用明显较弱,使得涡形成区长度进一步增大.Lin 等[13]用数值方法研究了亚临界雷诺数条件下具有相对较大展向波长的刚性波形柱周围的流动特征,也得到了类似的结果.

图6 展向涡量瞬时等值面(ωz=±1.0)的透视图和横流向俯视图

图7所示为不同扰动波高条件下的涡量场顺流向(ωx=±1.0)与横流向(ωy=±1.0)分量瞬时等值面云图.从该图可清晰分辨出顺流向涡结构和横流向涡结构在尾迹附近的演化特征.显然,与柔性光滑柱和A=0.10D扰动波高工况相比,A=0.20D时结构附近涡结构完全不同且涡量强度相对较低,并且沿下游耗散较快.Lin等[13]指出,附加流向涡在波形柱表面几何鞍点和节点处分别诱导出上涌流和下涌流,使得在几何鞍点处产生宽尾流,而在几何节点处产生窄尾流.这与本文模拟结果非常一致,表明波形结构的涡激振动抑制和减阻效应的力学机理有相同之处.

图7 3种工况下的顺流向(ωx=±1.0)与横流向(ωy=±1.0)瞬时等值面

为了进一步研究流向涡对结构响应的影响,图8中给出了3种工况下(x,y)平面内流向涡量切片图;其中图7中分别为柔性光滑柱和在A=0.10D,0.20D时柔性波形柱的模拟结果.该结果表明,柔性光滑柱两侧和A=0.10D波形柱两侧均形成一对同向旋转涡结构,而沿结构展向上、下侧分布的两对涡旋转方向刚好相反.然而,在A=0.20D情况下,在同一展向位置处柔性波形柱两侧生成一对较强的反向旋转涡,在这对反向旋转涡的外侧还分布着一对与各自强涡旋转方向相反的弱涡.Lin 等[13]指出,最优控制波形柱附加产生的反向流向涡对剪切层稳定性具有重要作用,能够防止剪切层与强涡结构的相互作用,使得涡旋形成长度增大.在A=0.20D的情况下观察到的反向旋转涡旋也稳定了剪切层,拉长了涡旋形成区长度,与本文之前观察到的现象一致.

图8 3种工况下(x, y)平面内圆柱周围的瞬时等值面和相应的流向涡量示意图

3.3 能量转换特征与频谱特性

结构的运动响应与流体与结构之间的能量转化特征密切相关,根据 Newman 等[11]的研究,一个旋涡脱落周期的无量纲时均能量E(z)可定义为

(20)

式中:Γ代表无量纲脱落周期;Wl代表升力产生能量的功率,即Wl=∂ζy/∂t.当E(z)为正值时,代表能量从流体转移至结构,E(z)为负值则转移方向相反.Newman等[11]和Zhu等[12]研究指出对于驻波响应,E(z)取值与展向位置有关,而与无量纲时间无关.同时,E(z)在一个脱落周期内沿整数波长的积分应该为0,表明在一个旋涡脱落周期内,结构从流体中得到的能量和在流体中做功耗散的能量相同,在没有外界能量输入的前提下,整个流体和结构系统满足能量守恒.因此,时均能量代表结构与流体进行能量交换的强度,其值越大则相应结构与流体之间的能量交换越频繁,反映到实际工程中,就会引起结构的疲劳,加速结构的老化.

图9所示为所有工况下时均能量随展向分布计算结果.由图9可知,A=0.10D柔性波形柱的时均能量比柔性光滑柱显著增大,当扰动波面高度增至A=0.15D其时均能量相比于柔性光滑柱维持在同一水准.而当A≥0.20D时,E(z)始终在0附近小幅振荡,说明在上述3种工况下波形柱与流体之间的能量交换接近于0,从能量角度验证了波形柱与流体之间的动力响应被显著抑制.同时,在z/D=0.37π～1.62π以及z/D=2.37π～3.62π之间柔性光滑柱E(z)为正值,在其余位置E(z)为负值;而A=0.10D波浪形结构在两端E(z)为正值,而在跨中附近E(z)为负值,这与图5中升力系数的时空演化模拟结果相互一致.

图9 细长柔性结构涡激振动时的无量纲时均能量分布

完全消除初始瞬态效应后,对3种典型工况下横向位移和升力系数时程进行了频域分析.

图10中分别为柔性光滑柱横向位移和升力系数功率谱密度(PSD)展向分布计算结果,可以据此确定旋涡脱落频率为fD/U∞=0.163,此外,在PSD的展向分布云图中可观察到,谱峰幅值与圆柱展向波长密切相关.Bourguet等[14]研究了长径比为200的细长柔性结构涡激振动问题,并也报道了类似现象.由此可以进一步证实,结构的振动响应与波长相关.此外,阻力频率(这里没有显示)是升力频率的两倍,说明尾迹是经典的2S模式.

图10 光滑圆柱无量纲时均能量分布

图11中分别为A=0.10D柔性波浪柱横向位移和升力系数PSD展向分布模拟结果,可以得知其主频为fD/U∞=0.161,略小于光滑柱时相应频率.从图11(b)可以看到,与光滑柱相比,A=0.10D柔性波形柱振动时程在跨中对应的PSD值明显减小,表明该位置附近升力的波动强度显著减小,与图5显示的现象一致.

图11 A=0.1D波浪柱无量纲时均能量分布

图12中分别为A=0.20D柔性波形柱横向位移和升力系数PSD的展向分布计算结果.对于A=0.20D柔性波形柱,除了主频fD/U∞=0.161外,在fD/U∞=0.126处也可识别出次峰频率,在图12用白色虚线表示.这种低频次生频率,在PSD-Cl中比在PSD-ζy中更强.次生频率产生的原因涉及流固耦合的非线性问题,其机理有待进一步探讨.

图12 A=0.2D波浪柱无量纲时均能量分布

4 结论

采用基于高精度谱单元方法的流固耦合分析方法,建立了细长柔性波形柱体结构的涡激振动力学模型,对其在低速均匀流作用和驻波初始动力激励下的流致振动机理进行了量化分析.结果显示,当扰动波高A≥0.20D时,柔性波形柱横向振幅相比于光滑柱明显减小,升阻力系数显著降低,水动力特性得到明显改善,表明在合适波高下,柔性波形柱具有良好的流致振动抑制作用.进一步比较扰动波高分别为A=0.10D,0.20D的结构横向振幅和水动力系数的时空演变特征,大致确定出波形柱抑制涡激振动的控制失效范围(A<0.10D)以及控制有效范围(0.20D≤A<0.30D).

为了探究扰动表面波高对结构振动响应的影响机理,对扰动表面结构和无扰动表面结构近尾迹三维涡结构进行了对比分析.结果表明,与柔性光滑柱相比,A=0.20D工况下交织涡结构完全消失,使得波形扰动截面剪切层的卷曲和相互作用明显较弱,涡形成区长度进一步增大.对柔性光滑柱在A=0.10D,0.20D时波形柱的流向涡结构瞬时等值面进行了模拟比较,结果表明在A=0.20D的情况下观察到的反向旋转涡旋也稳定了剪切层,拉长了涡旋形成区长度.比较了所有工况下时均能量的展向分布计算结果,发现当A≥0.20D时,E(z)始终在0附近小幅振荡,从能量转移角度进一步验证柔性波形柱对涡激振动的抑制作用.在完全消除初始瞬态效应后,对3种典型工况下横向位移和升力系数时程进行了频域分析,确定了3种工况下的漩涡脱落频率,对于A=0.20D波形柱工况下观察到的次生频率,其产生的原因涉及流固耦合的非线性问题,有待对其机理进行进一步探讨.