波浪作用下直立结构物附近强湍动掺气流体运动的数值模拟1)

邓 斌 王孟飞 黄宗伟 伍志元 蒋昌波,2)

*(长沙理工大学水利工程学院,长沙410114)

†(水沙科学与水灾害防治湖南省重点实验室,长沙410114)

**(宁波中交水运设计研究有限公司,浙江宁波315800)

引言

波浪从外海传播至海岸线附近的过程中,在地形和海岸防护建筑物的影响下,波浪与建筑物相互作用常弓起水流强烈的湍动和翻滚,伴随着大量空气卷入,产生大量气泡,该过程涉及到复杂的气液两相流问题.建筑物附近水体的掺气、以及气泡的演化和输运过程影响着近岸结构物的波浪冲击效应[1]、建筑物局部冲刷、波能的耗散等[2-5].与此同时,近岸带波浪受浅水效应容易发生破碎,破碎过程中会使大量气体卷入水体.因此,研究近岸波浪与海岸防护建筑物相互作用过程中掺气水流的运动特性,对准确计算建筑物受力[6-7]、越浪量[8-9]及其稳定性防护[10-11]等具有重要的理论价值和科学意义.

国内外学者一直致力于波浪作用下破碎区水动力及掺气特性研究,通过实验和数值手段获取了部分破碎区气泡输运的分布规律,并建立相关数学模型.其中,前人研究主要针对破碎区内空隙率、气泡直径以及数量等参数进行量化分析,部分学者基于固定坡度概化近岸海滩建立了描述破碎区空隙率的经验公式,一定程度上解释了冲泻区内空隙率、水动力特征及其与波浪要素之间的关系.如:邓斌等[12]通过光纤探针测量了卷破波作用下破碎区内空隙率的变化规律,给出了空隙率与波陡、相对水深以及水平位置的经验关系.Mori 等[4]通过现场测量与实验研究探讨了冲泻区内气体夹带的影响,发现气泡大小的空间分布规律与气泡大小密切相关,并提出气泡大小的Hinze 尺度为1 mm.Deane 等[13]对比分析了物理实验和现场观测两种研究方法下气泡分布的差异问题,并提出气泡概率密度与气泡半径之间呈现10/3 次方的关系.类似工作还有文献[14-21]等.

然而,上述研究主要集中于近岸波浪破碎区至冲泻区掺气特性与水动力特征的研究,未考虑建筑物的影响效应.与此同时,现有针对波浪与建筑物相互作用的研究大多未考虑水体掺气的影响[22-23].如:Xiao 等[24]基于不可压雷诺平均的N-S 方程(RANS)和k-ε 湍流模型建立了数值波浪水槽,研究了波浪作用于不同高度建筑物上波浪力和倾覆力的变化规律.Nakamura 等[25]基于实验和数值模拟方法研究了建筑物前冲刷坑深度与水流流速以及有效应力之间的关系.类似的相关研究还有文献[26-28].近年来,少量研究开始关注波浪与建筑物作用下水体掺气对建筑物受力的影响,如Hull 等[5]通过实验研究发现卷破波作用下气体卷入会导致建筑物受力增大的现象.Hsiao 等[29]通过实验和数值模拟研究了孤立波与梯形海堤相互作用的过程,定性分析了海堤后涡动结构及气泡夹带特性的演化特征.部分学者[30-32]研究发现跨海大桥桥板下方的主梁之间包裹的空气是增加波浪对桥梁上升力的一个重要因素.然而,上述研究仅针对水体掺气对建筑物受力影响的讨论,一定程度上推进了波浪与建筑物作用下掺气水流运动特性的认识.

综上所述,尽管目前对近岸带波浪破碎掺气特性以及波浪与建筑物相互作用等问题有了一定认识,建立了破碎区内掺气特性与水动力之间的关系,但准确理解掺气后形成的气泡与水动力特性之间的关系仍然是一项非常具有挑战的工作,有必要深入了解波浪与建筑物相互作用下掺气特性的内在运动机制.为此,在前人研究基础上,本研究通过建立三维精细数值水槽和对应物理实验模型,对波浪作用下建筑物附近的水动力以及相应的掺气特性进行研究,为进一步深入揭示波浪与建筑物相互作用、水体掺气与气泡输运机制奠定基础.

1 数学模型

1.1 控制方程与数值方法

基于三维连续不可压缩Navier-Stokes 方程建立数值模型.模型通过求解混合RANS (Reynolds average Navier-Stokes)和LES (large eddy simulation)的延迟分离涡湍流模型S-A DDES(Spalart-Allmaras delayed detached eddy simulation),其中边界区和湍流区分别采用RANS 和LES 模拟计算,该方法通过RANS模型可在近壁面处实现以较少的网格条件下得到较准确的解,又可通过LES 模型提高内部流场湍流的解析度.S-A DDES 模型是在S-A DES 模型基础上通过重新定义长度尺度,从而解决模型应力退化的问题,实现两种不同的湍流模型的有效结合,并解决区域之间的过渡问题,其中湍流区采用亚格子模型(SGS)对小涡进行计算.对Navier-Stokes 方程进行滤波处理后以实现大涡和小涡分离的控制方程为

式中,xi和xj为笛卡尔坐标下的分量,ui和uj分别为坐标系下对应的流速场,p为压强,ρ 为密度,在大涡和小涡的求解中,为亚格子应力项σ 为表面张力系数,κ 为界面曲率,α为VOF 流体体积函数.由于LES 中亚格子应力项未知,因此采用亚格子模型进行封闭,形式为,Sij为应变率张量,νt为湍流黏度系数为分子黏度),采用S-A 模型对进行封闭求解

式中,右侧第一项为涡黏性产生项,第二项为输运项,第三项为耗散项;Cb1,Cb2,Cw1和δ 均为常量;fw为与和d有关的函数;d为长度尺度,即网格点到壁面的最小距离.基于IDDES 方法修改S-A 湍流模型长度尺度d

其中,CDES为常数,取值0.65;fd为延迟函数,实现从RANS 模式到LES 模式的延迟;L△为网格尺度,L△=max(△x,△y,△z),△x,△y和△z为网格点3 个方向上的网格尺寸.

采用VOF 方法捕捉自由表面,弓入流体体积分数描述水和空气两种流体的分布,其中α=1 为水相,α=0 为气相,α 介于0 和1 之间为气体和水体交界面.同时,在传统体积输运方程的基础上本模型增加了人工压缩项[33],可以有效保证方程解在自由界面附近自动满足有界性,避免产生界面模糊,且对外部流场无影响,修正后的体积分数输运方程

1.2 边界设置

数值模型中,水槽底边界以及直立结构物壁面采用固壁边界条件,水槽前后两侧设置为对称性边界条件,水槽右侧设置为自由出流边界条件,水槽左侧采用Jacobsen 等[34]提出的修正速度入口方法进行造波,该方法中造波边界的类型分为干网格、湿网格和临界网格3 种,其中干网格的边界条件

式中,n为边界网格面上的单位法向量,湿网格边界条件(u和φB)是根据相应的波浪理论计算得到.模型对临界网格进行了湿面修正,波浪理论计算的自由液面曲线与网格边界存在2 个交点,取两点间直线作为临界网格的湿面边界,从而得到网格湿面部分的中心位置,以此计算出临界网格的入流速度和体积分数,这可避免了入口处自由液面的虚拟振荡问题.

1.3 模型与工况设置

数值计算模型设置如图1 所示.其中,岸滩概化地形分为两部分:一部分是距离造波边界15 m 处布置一光滑1:10 的斜坡;另一部分是在斜坡顶端的平台起始处布置直立结构物模型,模型高0.1 m、宽0.12 m,平台长3 m,右端接1:10 斜坡,用于消减波流能量.计算以斜坡起点作为坐标原点,以波浪传播方向为x轴正方向,以水平垂直数值水槽方向为y方向,以竖直向上方向为z轴正方向,由此建立三维坐标系,同时,将xoz坐标系作为二维坐标系.计算选择规则波作为入射波,具体工况设置详见表1 所示.

图1 试验布置及测量分布图Fig.1 Sketch of the experimental setup and instruments

表1 研究工况Table 1 Study cases

为验证数值模型的计算结果,在长沙理工大学水沙科学与水灾害防治湖南省重点实验室的波浪水槽中开展相应的物理模型实验.物理模型实验布置同数值水槽一致,水槽的左端配有活塞式造波机,斜坡地形和直立结构物模型均采用有机玻璃制作.模型实验沿水槽方向布置6 个波高测点(具体位置见表2),采集仪器为德国General Acoustic 公司生产的ULS 80D 型号超声波非接触式浪高仪,采集频率为50 Hz,并在实验前通过水槽内增减水进行标定,该仪器测量误差小于3.0×10-4m.与此同时,实验采用德国Photron 公司生产的高速相机(FASTCAM Mini UX100,采样频率:4000 Hz)记录不同波况下波浪破碎的演化过程,并采用法国A2 Photonic Sensors 公司生产的气泡测量系统(型号:B-POP,传感器为L 形光纤探针)测量相应位置处气泡的数量与空隙率,测量过程中通过记录气泡穿过光纤探针过程中电信号的历时曲线,从而获取气泡的数量以及空隙率,其详细测量原理及方法见文献[12],该光纤探针能够测量直径大于100 μm 以及运动速度为0.06～6.25 m/s 范围内的气泡,最大误差小于6%.气泡测量位置分别以5 cm 和1 cm 的间隔在水平x方向和竖直z方向上进行测量,具体测量位置详见表1.实验中所有实验仪器设备均采用同步系统同步采集实验数据,且物理实验在淡水和周围空气下进行.

表2 浪高仪位置Table 2 The location of wave gauges

1.4 网格划分与参数设置

由于气泡测量系统测量气泡直径的范围为2～70 mm,结合预备实验的结果,数值计算区域的网格沿x方向将结构物前后各1.0 m 区段设为核心区域,采用1 mm 的均匀结构网格,核心区域向外以结构网格的方式渐变过渡到5 mm;数值水槽在y方向上长0.08 m,沿y方向采用1 mm 的均匀网格,z方向同x方向类似在结构物附近网格大小为1 mm,核心区域向外以结构网格的方式渐变过渡,模型的计算网格总数达到2.4 亿,结构物附近网格划分的局部细节如图2 所示.

图2 结构物附近的网格划分Fig.2 Computation grids near the structure

为更加准确地模拟波浪破碎作用下气泡的运动,数值模型采用连续表面张力模型(CSF)考虑表面张力的影响,表面张力系数σ=0.072 N/m,模型相关参数见表3.数值计算采用自适应时间步长,以保证克朗数Cr小于1.对于单元网格下的库朗数的定义如下

式中,δt为时间步长,|U|为水质点通过单元网格的最大速度,δx为单元网格大小.

表3 数值模型参数设置(ρw和ρa分别为水和气的密度,µw和µa分别为水和气的黏滞系数)Table 3 Parameter setting(ρwand ρaare density for water and air respectively,µwandµaare viscosity for water and air respectively)

2 模型验证

图3 自由液面历时曲线对比Fig.3 Comparison between numerical and experimental time series of free surface elevation

图3 为不同测点处波面历时曲线的数值模拟结果与实验结果的对比图(h=0.34 m,H=0.08 m,T=1.6 s).采用S kill数验证计算值和实验值的吻合度,其中S kill数计算采用Willmott(1981)的方法

式中,Xnum和Xexp分别表示数值计算值和实验值,上画线表示取平均值.该工况下计算得到的S kill数均大于0.85,表明数值结果与实验值吻合较好,该数值模型能够真实地模拟规则波在斜坡上破碎的历时变化特征.

从图3 中还可以看出随着水深在斜坡上逐渐变浅,波浪发生浅水变形,波峰逐渐变尖,波形表现出明显的不对称性(3#),随后波浪在结构物附近发生破碎,破碎后的波高迅速减小,破碎波以段波形式继续传播.其中,5# 和6# 位置处数值模拟结果略高于实验测量结果,这是由于数值模拟结果中波高数据考虑了波面处的液滴飞溅.此外,对比实验照片图4(a)和三维数值结果图4(b),发现数值结果较好地再现了实验中波浪破碎后产生的强湍动掺气水流,以及先冲击结构物近海侧,随后沿结构物上爬,形成强射流现象,最终产生越浪,进一步表明建立的数值模型能合理模拟结构物附近波浪破碎后形成强湍动掺气水流的运动过程.

图4 数值与实验对比图Fig.4 Qualitative comparisons of free surface evolution between laboratory images and numerical data

3 结果分析与讨论

3.1 掺气水流与结构物的作用过程

图4 给出case2(h34H8)波浪与结构物相互作用下数值和实验过程的演化结果.结果表明波浪在斜坡上破碎后,空气被夹带进入水中形成一个大气腔和众多小气泡,当波浪持续破碎传播至结构物近海侧时,产生强掺混的掺气水流并冲击结构物,之后波浪沿结构物上爬,产生明显的射流和滴液的飞溅现象,最终造成结构物后方的越浪.此外,为分析掺气水流中气泡在输运过程中的破碎、溢破等现象,图4(c)给出了结构物附近气泡输运过程细节图(图4(b)中红色方框的局部图).从图4(c)中可知,波浪在堤前破碎时,众多气泡形成于近海侧.气泡在随后输运过程中气泡周围一直伴随着涡量的存在,并且大都是正负涡量相互作用.大气泡在经历多次破碎后形成小气泡,气泡在分裂为多个气泡时,中间的连接处往往伴随着负涡量,而周围为正涡量,最终在正负涡量的共同作用下气泡逐渐由单个分裂为多个.通过实验以及数值结果观察发现破碎后较大的气泡易受水动力因素影响,表面积和体积不断发生变化,浮力作用能驱使大气泡迅速上升至自由表面,最终在自由表面处破碎;而经历多次破碎后形成的较小气泡的表面积和体积很少改变,且受浮力影响较小,在水体中存留时间较长,直至水流平稳后,才得以上浮消散.

图5 为结构物附近气泡的演化细节图,其中左图可以看出在气泡的变形处分布着大量的正负涡量,表明气泡的形状甚至分裂受正负涡量的相互作用影响.从中、右图可见气泡群的分布范围与湍动能的分布范围大体一致,气泡群越密集,即气泡数量越多,湍动能越大,范围越广,表明水体中湍动能的分布与气泡群的数量密切相关;此外,还可以看出气泡内部流场与周围流体的流场略有不同,但整体流速与周围水体的流场相匹配,由此可推断,气泡的整体运动轨迹与气泡周围的水体流场密切相关.

图5 气泡演化细节(左:自由液面与涡量的分布;中、右分别为α=0.1,0.3,0.4 时自由液面,以及相应切片上的湍动能)Fig.5 The detail characteristic of local 3-D bubble(The distribution of free surface and vorticity is given in the left figure;The free surface at α=0.1,0.3,0.4,and the turbulent kinetic energy on the corresponding slice in the right figure)

3.2 流速场

图6 结构物附近流场Fig.6 Time series of velocity field near the structure

现有研究表明,波浪与结构物相互作用产生掺气水流,将直接影响流场分布特征[35].图6 给出不同时刻下case3 (h34H10)x–z平面内结构物附近的流场矢量图.t=8.2 s 时刻,波浪破碎后的水体撞击到结构物,形成向上飞溅的射流,在x=2.6 m 处形成一个较大的气腔,并在x=2.8 m 处通过水舌卷入大量的气体,形成顺时针的漩涡,此时水体流速较大,其值在2 m/s 左右,结构物附近掺杂了大量的小气泡,流场十分复杂;随后波浪沿结构物上爬形成射流(t=8.3 s),但向上的速度逐渐减小,波浪继续向前运动,翻卷的水舌继续夹带空气形成掺气水流;在x=2.8 m 处顺时针的漩涡将大气腔分裂为两个气泡,同样x=2.98 m 处的气腔在水流剪切作用下也分裂为多个气泡;随后射流逐渐回落(t=8.4 s),速度达到2 m/s 左右,在x=2.7 m 处的大气腔在水流作用下继续向前运动,而上一时刻在x=2.8 m 处的两个气泡融合为一个气泡,波浪逐渐在结构物前产生雍水,流速开始减小,水面处的水平流速减小至0 左右;射流完全落下(t=8.5 s),水体表面开始出现回流,水体底部仍然存在向前的水平流速.因此,在x=2.88 m处形成逆时针方向的旋涡,大气泡逐渐上浮、溢出.t=8.6 s 时,水体底部流速逐渐减小,水面处回流流速逐渐增加,水体容易产生较大的水流剪切力,结构物附近的旋涡范围逐渐增大,结构物前形成较多小气泡;水体底部开始回流(t=8.7 s),水体表面的反向水平流速较大,底部流速较小,小气泡逐渐上浮耗散,但仍存在大量漩涡;之后(t=8.8 s)水体内几乎没有较大的气泡,小气泡逐渐上浮,水体回流速逐渐增大.

3.3 涡量场

为进一步研究涡对气泡演化及水体掺气的影响.图7 给出不同时刻case3 (h34H10)工况下x–z平面内结构物附近的涡量分布图.通过y方向涡量反映涡结构的强度与方向.其中,y方向上的涡量表示为

图7 结构物附近涡量场Fig.7 Time series of vorticity field near the structure

式中,u和w分别为x方向和z方向的速度.

图7 中t=8.2 s 时刻,在x=2.6 m 处的气腔附近伴随着较大的涡量,x=2.8 m,2.98 m 处卷气区域周围正负涡量相间,顺时针与逆时针的漩涡相互作用容易形成对气腔的剪切.水体中的旋涡随着气腔向前运动(t=8.3 s～8.4 s),涡量逐渐减小.当射流完全落下(t=8.5 s),水体表面开始回流,大气泡逐渐上浮,水面附近展现出较多负涡量,水体底部呈现较多正涡量,印证了上文流速场的分析.当t=8.6 s 时,水体中负涡量逐渐增强,范围逐渐增大,结构物前气泡间存在明显负涡量;之后(t=8.7 s)由于水体的回流,结构物前的涡量与小气泡逐渐上浮消散.

3.4 湍动能

图8 为不同时刻case3 (h34H10)工况下x–z平面结构物附近的湍动能分布图.其中,湍动能表示为

式中,u′和w′分别为横向和纵向的脉动速度.

图8 结构物附近湍动能Fig.8 Time series of turbulent kinetic energy field near the structure

当t=8.2 s 时,在x=2.8 m 处的波浪破碎形成的卷舌与水体接触部分存在较大湍动能,射流与结构物前端接触区域的湍动能较大,其值为0.02 m2/s2.当卷舌运动到x=2.9 m 处(t=8.3 s),湍动能继续增大,呈现水面处以及结构物附近湍动能较大,水体内部以及破碎后的水体湍动能较小的特征.此外,x=2.78 m 处的卷气区域周围湍流动能较大,结构物附近的高速掺气水流湍动能较大(t=8.4 s).同样,t=8.5 s 时刻,在自由液面处以及结构物附近的掺气区域湍动能较大,表明水体的掺气能增强水体的湍动.随后(t=8.6 s)水体中流速逐渐减小,但靠近结构物的高度掺气区域的湍动能仍高于其他区域,这与Mori 等[4]的研究类似.

3.5 空隙率、气泡数量与湍动能的关系

上述讨论表明,湍动能高的区域与气泡数量较多的区域具有较高空间重合性,表明湍流强度与气泡数量之间存在一定相关性.实际上,已有研究也证实破碎波中湍动能与气水混合物的空隙率呈现正比关系[4].为探讨结构物附近掺气流体中空隙率与湍动能是否存在类似关系,图9 给出任一周期内通过L型光纤探针测量到各测点的平均空隙率与平均湍动能的变化趋势.可见,3 个工况下湍动能与空隙率整体也呈现线性关系,值得注意的是,随着波高的增大,空隙率也逐渐增大,但空隙率与湍动能的线性关系有所减弱(如case3),这是由于图中相对分散的点均是离直立结构物较近位置(x=2.96,3.00 m)和结构物上(x=3.04 m).分析可知,与斜坡上(无结构物)破碎波内掺气水流特性不同的是,结构物的存在使得湍动能与空隙率的线性关系变弱,这与波高增大时,破碎波与直立结构物作用产生的流体飞溅和越浪导致空隙率在较短时间内发生剧烈变化,以及波浪破碎及其与结构物非线性相互作用更加明显等相关.图10 给出了实验测量得到的气泡数量与空隙率之间的关系,结果同样显示两者整体存在正相关关系,从而可推断湍流动能与气泡数之间存在一定关系.因此,图11 给出了任一周期内通过L 型光纤探针测量到各测点的平均气泡数与平均湍动能之间的变化趋势,结果同样显示湍动能与气泡数整体存在正相关关系,且随着波高的增大气泡数逐渐增多,印证了上文关于气泡数与湍动能的分析.为更好地描述各回归方程的可靠性,表4 给出了图9～图11 中各回归方程参数的相关统计值.综上可见,波浪作用下直立结构物附近湍动能的分布与与掺气水流部分特征参数(气泡数量、空隙率)整体呈现正相关关系.

图9 时间平均空隙率和时间平均湍动能关系Fig.9 Relationship between time average of void fraction and turbulence kinetic energy

图10 时间平均空隙率和时间平均气泡数关系Fig.10 Relationship between time average of void fraction and bubble number

图11 时间平均空隙率和时间平均湍动能关系Fig.11 Relationship between time average of bubble number and turbulence kinetic energy

表4 线性回归方程相关参数统计表Table 4 Parameter statistics of linear regression equation

4 结果与展望

基于Navier-Stokes 方程的三维数值模拟技术,结合物理模型实验,模拟了波浪作用下结构物附近强湍动掺气水流的水动力学特性,研究表明:

(1)建立的数值模型能精确地捕捉波浪作用下结构物附近的自由液面的变化以及气泡输运过程,较好地描述气体卷入所形成的气腔形态以及多气腔之间的融合、分裂等过程.

(2)分析结构物附近波浪破碎后的相关流场和湍动能、涡量场,发现破碎波与结构物相互作用过程中水体内湍动能以及涡量的变化与水体内气泡变形、破碎与融合等现象相关.其中,气泡云的范围与水体湍动能的在空间分布相重叠,正负涡量的相互作用能弓起气泡的破碎和融合,气泡的融合、破碎又会弓起水体流场的变化.

(3)研究建立了结构物附近湍动能与掺气特性相关参数之间的关系,发现波浪作用下结构物附近湍动能的分布与掺气水流部分特征参数(如气泡数量、空隙率)呈现线性正相关关系.

基于数值模拟和物理实验揭示了波浪作用下直立结构物附近强湍动掺气水流中气泡特性与水动力之间的关系,为进一步认识波浪作用于结构物的内在机制提供了理论基础,但文中仅考虑了定床条件下特定波浪参数以及概化的结构物,考虑更多岸滩特征、结构物型式、波浪参数等条件下强湍动掺气水流运动规律及其与结构物相互作用机理还有待进一步研究.