基于重叠网格的桥墩防撞浮箱流场特性数值研究

2023-05-10舒麟棹钱志鹏万德成

重庆交通大学学报(自然科学版) 2023年3期

冀楠,杨光,舒麟棹,钱志鹏,万德成

(1. 重庆交通大学航运与船舶工程学院,重庆 400074; 2. 上海交通大学船海计算水动力学研究中心(CMHL),上海 200240)

0 引言

随着我国水路网络运输的发展,船撞桥事故频繁发生,若安装上防撞装置则能有效地保护桥墩安全[1]。桥墩防撞装置种类繁多,其中一种是自浮式防撞浮箱,它能随着水位升降而调整位置,特别适用于大变幅水位区域内的桥梁防撞保护,目前被广泛用于国内的桥梁防船撞设计中[2]。工程界针对防撞浮箱的设计一般是从撞击概率、结构变形、承载载荷和经济性等方面进行考虑[3],对防撞浮箱-桥墩系统则主要从碰撞角度研究浮箱防撞性能和桥墩结构损伤,极少会从水动力学角度考察防撞浮箱的流场特性及对桥墩影响[4-6]。因此,研究防撞浮箱在不同来流速度和水深条件下与桥墩之间的水动力性能及流场特性具有重要的理论意义和工程价值。

针对防撞浮箱与桥墩之间的研究可从两个方面展开:① 研究防撞浮箱的防碰撞性能,例如白凯等[7]采用显式有限元法模拟了船舶撞击防撞设施的过程,分析了防撞浮箱的耐撞性、变形及桩所受的力;单成林等[8]利用所建立的船舶、防撞浮箱、桥墩碰撞系统计算模型,对不同角度撞击该浮箱和桥墩的受力影响、能量转换、拨动船头程度进行了分析;单成林等[9]采用静力等效方法对SPS圆形防撞浮箱进行了受力分析。②研究裸桥墩的绕流特性,例如YU Peng等[10]基于N-S方程对桥墩周围的复杂流动进行了求解,模拟出了桥墩周围局部冲刷的动态演变过程;吴承伟等[11]对典型等直径45°错置双柱桥墩在不同柱间间距下的三维流场及水流力特征等进行了研究;F.M.CISTERNAS等[12]采用大涡模型对垂直圆形桥墩周围水流进行了三维数值模拟,研究了傅汝德数、雷诺数和几何参数(水深与桥墩直径之比)等3个参数对自由表面特性的影响。

笔者采用STAR-CCM+软件,基于重叠网格技术,以圆形自浮式浮箱-桥墩系统为研究对象,采用VOF方法捕捉自由液面,研究了水流流速和水深变化对防撞浮箱与桥墩间水动力相互作用和流场的影响,为防撞浮箱工程应用决策提供相应参考。

1 数值计算模型

1.1 控制方程

在空间固定坐标系下,数值模拟采用不可压缩流动RANS方程,连续性方程和雷诺平均方程分别如式(1)、式(2):

(1)

(2)

1.2 湍流模型

陈志乐等[13]研究表明:对于床面上直立圆柱绕流问题,RSM模型能更好地模拟出速度场、床面剪应力分布、柱侧分离角等流动现象。RSM模型对复杂流动有着更高的预测精度,摒弃了各向同性的涡黏性假设,通过附加雷诺应力运输方程和耗散率方程来封闭N-S方程,其控制方程如式(3):

(3)

式中:u′i、u′j分别为i、j方向速度;Dij为扩散项;Pij为产生项;Πij为压力应变项;εij为耗散项。

其表达式为:

(4)

(5)

(6)

(7)

2 模型验证与研究对象

2.1 模型验证

模型验证采用圆柱绕流实验[14],圆柱直径D=0.15 m,验证算例计算域长度取35D,圆柱上游入口距离圆柱中心为10D,下游出口距离圆柱中心为25D。计算域宽度为10D,左右两侧距圆柱中心都为5D,水深为0.2 m,水流流速为0.26 m/s。为了考察网格数量对计算结果的影响,通过改变网格基本尺寸,得到了3种质量不同的网格,其具体参数如表1。

表1 网格基本参数信息

实验测量了圆柱下游不同位置处的时均速度横向分布,如图1。图1中:时均流速可通过除以横向最大流速进行无量纲化,其中x/D表示圆柱下游距离圆柱柱面的相对位置,y/D表示圆柱左右两侧距离圆柱柱面的相对位置,z/H表示水深位置距离床面的相对高度。图1(a)～图1(c)表示圆柱下游不同位置近自由表面处的速度分布,图1(d)～图1(f)表示近床面处相应位置的速度分布。可以发现,网格数量对计算结果影响比较小,这3套网格数值计算结果与实验值非常吻合,仅图1(f)中的圆柱正后方位置(-0.6≤y/D≤0.6)速度分布误差较大,原因可能是床面的粗糙度设置与实验存在差别,导致近床面和圆柱远后方水流流态与试验结果存在一定误差。综合考虑计算精度与时间成本后,选择中等网格属性设置参数作为后续工况计算的网格划分参数。

图1 圆柱下游不同位置处的时均速度横向分布

2.2 研究对象

自浮式浮箱作为一种典型的桥墩防撞装置被广泛应用于实际工程中(图2)。自浮式防撞浮箱通过滑轮组件随河道水位涨落进行自动升降,从而达到在不同水位条件下保护桥墩的作用。

图2 自浮式防撞浮箱

笔者选择圆形自浮式浮箱-桥墩系统为研究对象,并进行数值建模计算。浮箱-桥墩几何模型如图3,其中图3(a)为防撞浮箱被静水面相切的圆截面,坐标原点O位于截面圆心处。实际内河航道桥梁直径大多在1～10 m,为了与实验验证工况保持一致,选择直径为7.5 m的大尺度桥梁为研究对象,采用1∶50的缩尺比,则桥墩模型直径D=0.15 m。浮箱实际尺寸参考文献[15],按缩尺比计算后得到圆形浮箱宽度L=0.07 m,高度h=0.06 m。为简化计算模型,浮箱与桥墩之间的滑轮组件和护舷均做简化处理,间距d=0.02 m,浮箱设计吃水为0.04 m。

图3 浮箱-桥墩系统

长江航道各段水流流速约为1～3 m/s,根据傅汝德相似定律,得到计算模型水流流速约为0.14～0.43 m/s,但实际河流中的水流流速是不断变化的,为尽可能覆盖较宽的流速范围,笔者对水流流速取0.13、0.26、0.52 m/s这3组工况。长江航道各段平均水深一般在2～10 m,以长江上游5000DWT过闸推荐船型为参考,其最大设计吃水为4.3 m,因此笔者将研究水深确定为5～10 m,按缩尺比计算后的3组水深工况分别设置为0.10、 0.15、 0.20 m。

2.3 计算域及边界条件

整个计算域分为背景域和包含防撞浮箱的重叠域,如图4。其中计算域的入口距离桥墩中心10D,计算域的出口距离桥墩中心30D,左右两侧距离桥墩中心各为10D,自由液面以上空气高度为0.1 m,采用VOF方法捕捉自由液面。计算域入口与顶部均设置为速度入口,均匀来流;出口设置为压力出口;计算域底部及左右两侧均设置为无滑移壁面。数值求解时间采用一阶离散;扩散项采用二阶中心差分格式离散;对流项采用二阶迎风格式,以提高计算精度;压力-速度耦合方程采用SIMPLEC算法进行求解。

图4 计算域

2.4 重叠网格划分

为捕捉防撞浮箱在不同外界条件下的升沉垂荡运动特性,笔者采用重叠网格法进行数值模拟。将防撞浮箱等效为一个浮体,在计算过程中只放开垂向自由度,计算域重叠网格划分如图5。

图5 计算域网格划分

网格划分采用切割体网格,边界层厚度取0.017D,对自由液面、桥墩周围进行体加密,同时对重叠域进行体加密,保证背景域与重叠域之间的网格平滑过渡。设置防撞浮箱运动域和桥墩流场背景域两套网格,两套网格相互独立,并通过交界面进行传值来实现数据交换,在计算过程中不会产生网格畸变和负体积等问题,能始终保持计算网格质量,从而保证计算结果准确性。

3 计算结果分析

3.1 防撞浮箱运动特性

防撞浮箱的升沉位移时程曲线如图6。由图6可知:水流流速主要影响防撞浮箱的升沉位移值,而水深主要影响防撞浮箱的垂荡历程。流速变化时,从图6(a)～图6(c)横向对比可看出:浮箱稳定后的升沉位移值变化很大。水深变化时,在这3个流速下的浮箱稳定后升沉位移值趋近相同,但其垂荡过程有所区别。当流速U=0.13 m/s时,防撞浮箱下沉量很小,且稳定后几乎不再垂荡;当流速U=0.26 m/s时,防撞浮箱下沉量增加,稳定后有小幅度垂荡;当流速U=0.52 m/s时,防撞浮箱下沉量进一步增加,且稳定后出现大幅度垂荡。当水深减小时,防撞浮箱在10～25 s区间内垂荡程度是逐渐加剧的;25 s后,浮箱垂荡逐渐平稳,但水深变化影响着垂荡历程。

图6 防撞浮箱升沉位移时程曲线

为更深入研究水深对防撞浮箱垂荡影响,笔者选择流速U=0.52 m/s时3个水深下的升沉位移,取计算稳定后的40～60 s垂荡数据,并通过快速傅里叶变换(FFT)求出频谱分析的结果,如图7。

图7 不同水深下防撞浮箱垂荡频谱分析

从图7中可看出:这3个水深下曲线均出现多个谱峰,说明有多个不同的垂荡频率。当水深H=0.20 m时,最高谱峰对应频率为0.2 Hz;当水深H=0.15、0.10 m时,最高谱峰对应频率均为0.1 Hz,说明水深减小时,防撞浮箱主要垂荡频率在减小。

3.2 防撞浮箱动力响应

为研究不同条件下防撞浮箱的动力响应特征,对防撞浮箱的横向力、纵向力和垂向力进行无量纲化处理,力系数Ci计算如式(8):

(8)

式中:ρ为流体密度;U∞为不受干扰的流体速度;S为防撞浮箱的湿表面积;Fi为防撞浮箱在不同方向上的力。

不同流速下防撞浮箱的时均力系数(纵向力系数Cx、横向力系数Cy、垂向力系数Cz)随水深变化的曲线如图8。

图8 不同流速下防撞浮箱时均力系数随水深的变化

由图8可知:随着水深增加,防撞浮箱的力系数基本呈下降趋势。当流速为0.13 m/s时,3个方向的力系数随着水深增加,其下降趋势较为平缓,Cx最大且随水深增加先略微增大后再减少,Cy和Cz最小;当流速为0.26 m/s时,3个方向的力系数大小居中,Cx随水深增加先略微增大后再减少;当流速为0.52 m/s时,Cx最小,Cy和Cz最大,力系数均随水深增加而下降明显。

当流速U=0.26 m/s时,防撞浮箱纵向力系数随水深变化的时程曲线和频谱分析结果如图9。从图9中可知:Cx随时间发展表现出强烈的非周期性脉动变化。从频谱分析结果可看出:3个水深下均出现了多个谱峰,且相邻谱峰间的能量峰值很接近,这也说明Cx无周期性的强烈脉动变化。

图9 防撞浮箱纵向力系数的时程曲线和频谱分析(U=0.26 m/s)

图10分别为3个水深下防撞浮箱Cz随流速变化的时程曲线和频谱分析结果,Cz的大小将直接影响着防撞浮箱的升沉垂荡。

在3个水深下,浮箱Cz的时间平均值随流速增加而变大,会导致流速增加时,防撞浮箱的升沉位移值增加,这一规律与上文关于防撞浮箱升沉位移和垂向力系数的分析一致。从频谱分析可以得出:水深H=0.10 m,流速在0.13、0.26 m/s时,频率谱有两个谱峰,一个主频率分别为2.20、2.15 Hz,一个副频率分别为1.70、1.35 Hz;当水深变化到0.15 m时,低流速(U=0.13 m/s)下的频率谱变成了单峰,而中流速(U=0.26 m/s)下的频率谱仍保持着两个谱峰;当水深增加到0.20 m时,低流速(U=0.13 m/s)下的频率谱也为单峰,但其频率减小到1.8Hz,说明水深增加后,低流速(U=0.13 m/s)下浮箱垂向力的脉动变化周期规律在增强,并且变化周期在增大。而中流速(U=0.26 m/s)下的频率谱出现了3个谱峰,并且最高峰值的能量在变大;当流速U=0.52 m/s时,3个水深下的Cz频域曲线都出现了多个谱峰,频率分布分散在0.50～4.55 Hz,这说明在高流速(U=0.52 m/s)下,防撞浮箱垂向力的周期性脉动规律完全消失,产生了强烈的垂荡。

图10 防撞浮箱Cz时程曲线和频谱分析

3.3 防撞浮箱-桥墩流场

用VOF方法捕捉的自由液面细节如图11。由图11可发现:浮箱与桥墩的间隙充满了水流,受浮箱干扰,桥墩尾部的流态紊乱。

图11 防撞浮箱-桥墩自由液面细节

图12为流速U=0.26 m/s、水深H=0.20 m时的防撞浮箱-桥墩三维马蹄涡。

由图12(a)可发现:防撞浮箱与桥墩共生成了3条马蹄涡,分别位于防撞浮箱自由液面处、防撞浮箱底平面处和桥墩床面处。由图12(b)可发现:桥墩床面处的马蹄涡涡量强度最大,呈粗管状形态,这也是引起桥墩床面冲刷的原因;防撞浮箱底平面处的涡量强度次之,形态扁平,这可能与防撞浮箱震荡有关;防撞浮箱自由液面处的涡量强度最小,呈纤细管状形态。

图13为流速U=0.26 m/s、水深H=0.20 m时防撞浮箱-桥墩在中纵剖面上的流线分布;图14为防撞浮箱-桥墩局部Q准则和流线。总体来看桥墩近后方流线呈螺旋运动形态,说明这一区域为回流区,随着流线逐渐往下游发展,螺旋形态漩涡消失,流线变为规则的左右摆动。

图12 防撞浮箱-桥墩三维马蹄涡

图13 防撞浮箱-桥墩中纵剖面流线

图14 防撞浮箱-桥墩局部Q准则和流线

由图13中的A、C区域和图14可看出:来流侧受到浮箱和桥墩的阻碍产生壅水,水面升高,水流速度降低,导致下游压力大于上游来流压力,产生了逆压梯度,一部分水流在逆压梯度作用下产生了回流,在防撞浮箱和桥墩来流侧形成了周向分布的流动分离漩涡,这就是防撞浮箱自由液面处和桥墩床面处产生马蹄涡的原因。其次,由图13中的B区域和图14可看出:遇到阻碍的水流还有一部分会分离下沉到浮箱底部,形成回流漩涡,这就是防撞浮箱底平面处产生涡量强度较高马蹄涡的原因。同时,防撞浮箱上下垂荡运动会产生加强和抑制漩涡的相反效果,使得浮箱底部马蹄涡形态不稳定。

图15分别展示了不同流速、不同水深共9个工况在自由液面(Z=0)处的二维瞬时涡量。

图15 自由液面处二维瞬时涡量

由图15可知:防撞浮箱垂荡主要产生两种类型漩涡:① 防撞浮箱后部边界层分离后形成的脱落涡,② 防撞浮箱升沉运动和水面相互作用形成的小尺度涡。当流速较小(U=0.13 m/s)时,小尺度涡分布在防撞浮箱周围,并在向下游运动过程中被拉长为条状涡,如图15(a)～(c);当流速增大,小尺度涡分布向桥墩后部的扇形区域转移,如图15(a)、(d)、(g);同时水深影响也很大,大流速加上小水深,小尺度涡的分布区域会大大扩展,如图15(i);由于受小尺度涡影响,在小流速(U=0.13 m/s)下,脱落涡还能呈现规律性的涡街形式,在中等流速(U=0.26 m/s)下,部分脱落涡被拉伸,呈现涡带,甚至涡丝形态,如图15(c)、(f);在大流速(U=0.52 m/s)下,脱落涡强度继续增大,但涡的脱落个数无明显规律。

图16分别展示了不同流速、不同水深的9个工况在半水深(Z=-H/2)处桥墩截面二维瞬时涡量。

图16 半水深处二维瞬时涡量

由图16可知:受水面处防撞浮箱垂向运动的影响,截面后部也存在大量小尺度漩涡,且漩涡的强度随流速增大和水深减小而逐渐增强;当水深较浅时(H=0.10 m),桥墩后部的小尺度漩涡与尾后的分离剪切层融合发展后,脱落为上下两对正旋涡和逆旋涡,且涡的强度随着流速的增大而增大,如图16(c)、(f)、(i),这也说明了图8(a)中纵向力系数随流速增加而增大的规律;当水深较大时(H=0.20 m),桥墩后的涡旋从强度较弱且被拉伸的涡带,逐渐过渡成强度较大的2S形漩涡,如图16(a)、(d)、(g)。值得注意的是,在水深H=0.10 m时,桥墩首部出现上下各一组条形漩涡,如图16(c)、(f)、(i),这是由于水深较浅时,半水深处桥墩截面距离防撞浮箱很近,受防撞浮箱底部的马蹄涡影响所导致。随着流速增大,条形漩涡被拉长,且在桥墩后部边界层分离点处又出现上下各一组较小的条形漩涡,如图16(f),随着流速继续增大,条形漩涡继续被拉长,直至与桥墩尾后分离剪切层发生融合,如图16(i)。

图17为防撞浮箱-桥墩三维瞬时涡结构等值面对比。

由图17可知:随着水深增大,涡量场的横向宽度在逐渐减小;当流速增大,马蹄涡逐渐与防撞浮箱后部漩涡发生融合;当流速较小(U=0.13 m/s),防撞浮箱后部漩涡融合程度较高,呈现尺度较大的涡带形态,在中等流速下(U=0.26 m/s),防撞浮箱后部漩涡形态多样,且尺度较小、数量增多,在大流速下(U=0.52 m/s),防撞浮箱后部漩涡整体呈现较规律的卡门涡街形态。