冀 楠,黄浩东,罗 意,李浩然

(重庆交通大学 航运与船舶工程学院,重庆 400074)

0 引 言

丁坝是常见的防洪护岸水工建筑物,在长江等航道中分布广泛。丁坝迎水面和背水面出现壅水、旋涡、下潜等三维特性,其束窄作用也使得丁坝流场水流形态变得复杂。船舶驶经丁坝时,会受到丁坝产生的水动力干扰作用,对船舶的受力特性和航行的姿态产生较大的影响。不同的水流流速和船-丁坝相对位置是影响船舶水动力的重要因素,因此合理的预测船舶驶经丁坝时的水动力特性对于船舶的操纵乃至安全航行有非常大的意义。相比于试验研究,数值模拟的方法更具成本和时间效益,并且对空间没有固定的约束。

国内外对丁坝的水流特性以及船舶经过桥墩、船闸及限制水域等的水动力特性进行了较为丰富的研究,而鲜有研究船舶行驶过丁坝时的水动力特性。N.MUNETA等[1]采用定床实验得到丁坝流场不同水深处的横向和纵向流速,通过实验值和计算值的对比,来分析丁坝区的三维特性;R.MAYERLE等[2]采用三维数值模拟,并假设采用静水压力分布来解释丁坝附近水流模式的差异,结合丁坝坝头的静水压力的分布来分析对丁坝背水面涡旋分布的影响;程年生等[3]通过使用k-ε湍流模型模拟了丁坝区的二维流场特性,得出了丁坝流场湍动能的分布和湍动能的耗散率;许光祥等[4]通过实验数据对丁坝断面的水位进行分析,得到了计算丁坝临近断面水位变化值与丁坝束窄系数关系的经验计算公式,为预测丁坝整治效果提供了有效的参考;马峥等[5]通过对不同湍流模型在流场中的应用进行总结,讨论了CFD中船舶的湍流模型适用性以及选择方法。由于方形系数相对较小,对于集装箱类船,标准k-ε模型和SSTk-ω模型模拟效果相对较好;麻绍钧[6]通过分析近岸航行时船舶的水动力特性,对流场细节进行分析,提供了一系列水动力机理分析方法;张晨曦等[7-8]采用RNGk-ε湍流模型,分析了船舶经过桥墩以及船-船相遇时的受力特性,并改变船-桥墩以及船舶之间的横向距离对船舶的横向力和转艏力矩的影响,研究发现船舶经过桥墩时所受的横向力会出现明显的峰值,转艏力矩会出现由负到正的转换,并且随着船舶重心到桥墩的距离的减少,横向力及转艏力矩增加明显。船舶之间纵向距离为一倍船长左右时容易发生碰撞;李佳等[9]利用ADV(三维点式多普勒流速仪)检测丁坝流场,并通过测力天平和应变数据收集仪器,测量船模在丁坝流场不同水流流速和不同船-丁坝相对位置下船模的横向及纵向二维受力状况。试验表明船舶经过丁坝时受到的推力和吸力出现明显的峰值,船舶位于丁坝轴线位置处时受力较大。

笔者利用计算流体力学软件STAR-CCM+,用VOF法处理自由液面,采用不可压缩流动RANS方程以及Realizablek-ε湍流模型对丁坝的流场和船舶的受力状态进行数值模拟,并与试验对比以验证数值模拟的准确性。然后用该数值模型模拟研究了船舶在丁坝流场航行时,不同水流流速和船-丁坝相对位置时船舶的受力情况。

1 数值方法

基于RANS方法求解不可压缩黏性流体流场的控制方程为雷诺平均连续性方程和雷诺平均N-S方程为:

(1)

(2)

采用Realizablek-ε湍流模型封闭上述方程组,其中,湍流动能k和湍流耗散率ε的方程分别是:

(3)

(4)

式中:k为湍流动能;ε为湍流耗散率;μt为湍流粘性系数,其中μt表达式为:

(5)

式中:Cμ=0.09;C1=1.44;C2=1.9;σk=1.0,σε=1.2。

2 数值验证

2.1 丁坝流场模拟

2.1.1 计算区域

图1(a)为丁坝计算模型的坐标系及边界条件示意图,丁坝流场模拟的计算域设置参考J.JEON等[10]的实验,并与其进行对比验证。丁坝为直丁坝,坐标系OXYZ原点位于丁坝1/2厚度面、丁坝所在岸壁侧面和计算域底部平面的交点处,且X轴指向河道下游,Z轴垂直于底部平面朝上。为了更准确说明丁坝和船舶之间在X方向上的位置关系,特定义丁坝1/2厚度面和河道底面的交线为丁坝轴线。水槽模型长为6.5 m,宽为0.9 m,高为0.262 5 m。水槽的初始水深为0.21 m,丁坝高度H为0.262 5 m,丁坝厚度D为0.04 m,丁坎长度L为0.3 m,为非淹没式直丁坝。计算域水流入口位于丁坝轴线前1.5 m的位置,入口处水流流速为U0=0.144 m/s,进口水流量恒定为Q=0.027 m3/s。丁坝模型的入口和顶部为速度进口(velocity inlet),出口为压力出口(pressure outlet),其他边界都为壁面(wall)。

图1 网格及探针线布置示意Fig. 1 Schematic diagram of grid and probe line layout

2.1.2 计算模型验证

图1(b)为探针线的布置示意图和丁坝模型Z=0.105 m平面以及丁坝纵截面网格图。选取丁坝附近流速变化较大的位置设置4根探针线,分别位于X/L=-3.30、-0.90、1.67、3.33处,探测范围为Y=0～0.9 m。其中探针线布置从左至右依次是探针线1、2、3、4。

从图2中可以看出,s2和s3的模拟结果基本一致,整体的模拟结果要好于s1,在综合考虑计算精度与计算资源的限制后,最后采用s2的网格参数作为后续工况计算的网格设置。

图2 不同探针线处X向流速模拟值与实验值比较Fig. 2 Comparison between the simulated values of X-direction velocity at different probe lines and the experimental values

2.2 船舶阻力数值计算

2.2.1 计算模型

选取KCS船模为研究对象,该船模为2010年哥德堡研讨会[11]推荐的标准实验船型,船模的垂线间长Lpp为7.278 6 m,与实船的缩尺比为1∶31.6。KCS虽为海船船型,但其船型的长、宽、吃水与长江过闸船型集-18[12]的相应参数基本呈2∶1的尺度比,二者船型相似;另外,虽然KCS船型的方形系数0.65比集-18的方形系数0.53大,但根据船舶原理[13]中关于方形系数对船舶阻力影响的叙述可知,在下文低于1.2 m/s的水流流速范围内,方形系数差异造成的船型阻力变化仅为10%左右,因此,采用KCS船型进行内河船的水动力特性研究是可行的。

2.2.2 计算域与边界条件

计算域示意如图3。计算域的长度为4倍Lpp,左右边界和上下边界长度为3倍Lpp,水深为2倍Lpp,入口(inlet)、顶部(top)、底部(bottom)为速度入口,出口(outlet)为压力出口,左右边界(side、symmetry)为对称边界。在出口边界进行了消波。

图3 计算域Fig. 3 Computational domain

2.2.3 网格划分

计算域的网格划分如图4。为保证计算结果的精确性,在船周围进行局部加密。同时,为捕捉船舶的兴波波形,不仅对自由液面进行整体加密,还在尾流区进行了局部加密。网格划分Y+值控制在30～200内。

图4 计算网格Fig. 4 Computational grid

2.2.4 网格无关性与模型验证

表1 模拟总阻力系数与实验对比Table 1 Comparison of the simulated total resistance coefficients and experiment

(6)

式中:Cd为阻力系数;Fd为总阻力,N;v为水流流速,m/s;A为船舶的湿表面积,m2。

在综合考虑计算精度以及计算资源的限制后,最后采用S2的网格参数作为后续模拟的网格设置。

3 计算结果及分析

3.1 计算工况

分别设置了3个水流流速和9个船-丁坝相对位置,其中:3个水流流速分别为v=0.6、0.9、1.2 m/s(对应的傅汝德数Fr=0.135、0.203、0.271),船-丁坝相对位置定义为船舶重心距离丁坝轴线位置的距离,其值分别为X/Lpp=-2.00、-1.00、-0.50、-0.25、0、0.25、0.50、1.00、2.00,同时放开船舶升沉、纵倾、横倾3个自由度。图5为坐标系、边界条件和船只摆放位置示意。(以船舶位于X/Lpp=0位置处为例)。

3.2 丁坝后流场的时域分析

图6为水流流速0.9 m/s,计算时间t=350 s时,船舶位于不同船-丁坝相对位置处的瞬时流场图。如图6,船舶的存在会对丁坝后的流场产生扰动,无船时的纯丁坝绕流会在丁坝后形成稳定单一的漩涡,但有船时丁坝后的高流速带(上部区域)和低流速带(下部区域)不再稳定,发生横向交替摆动,流速带的摆动和漩涡变化分析如图7;在同一时刻点,船舶位于不同船-丁坝相对位置处的流场形态并不相同,说明船-丁坝相对位置会对丁坝后流场的摆动周期产生一定影响。

图6 v=0.9 m/s时船舶重心平面流速分布Fig. 6 The plane velocity distribution of the ship’s center of graving when v=0.9 m/s

图7 v=0.9 m/s,X/Lpp=0.25时坝后流场周期性摆动示意(t=373～423 s)Fig. 7 Schematic diagram of periodic oscillation of flow field behind the spur dike when =0.9m/s, v=0.9 m/s, X/Lpp=0.25,(t=373～423 s)

图7为水流流速为0.9 m/s,船舶位于X/Lpp=0.25工况点下,坝后高流速带和低流速带在一个周期T内的摆动情况示意。可以看出,在一个周期内,上部的高流速带中先形成高流速核心区,如图7(d),之后又逐渐消失;低流速带中的2个漩涡核心则先融合成一个与高流速核心区上下对应的单一漩涡,之后又分裂为2个漩涡。

3.3 船舶受力分析

3.3.1 船舶横向受力沿程分析

图8(a)为不同流速下,船舶在距离丁坝不同位置处的横向受力曲线图。船舶所受的横向力的正负规定为,沿Y轴正向的力为正,反之为负。从图8(a)可以看出,随着水流流速的增加,船舶在靠近丁坝的位置处(-1

图8 不同流速下船舶沿程受力及力矩变化情况Fig. 8 Changes in forces and moments of ships along the way under different flow speeds

从图9的压力及速度云图中可以看出,船舶从X/Lpp=0.5位置开始,逐渐从低压区向高压区移动,船舶位于X/Lpp=-0.25位置时,由于丁坝的绕流,此时船舶靠近丁坝一侧与远离丁坝一侧压强差最大,船舶的存在阻挡了水流的横向发展。随着水流流速的增加,横向力峰值出现的位置有向下移动的趋势,在0.6、0.9 m/s流速下,横向力的峰值出现在X/Lpp=-0.5附近,而流速增大为1.2 m/s时,横向力出现峰值的位置向丁坝轴线移动,出现在X/Lpp=-0.25附近,此时丁坝的绕流造成的横向水流流速较大。而X/Lpp=0位置相对X/Lpp=0.5位置,船舶左右两侧的压力数值均急剧减小,但船舶左侧高压区的面积减小更快,这使得虽然船舶仍受到数值为正的横向力,但船体两侧的压力差减小了,也就是船舶所受横向力减小了。

图9 v=1.2 m/s时船舶重心平面压力分布及流速分布(局部)Fig. 9 The plane pressure and velocity distribution of the ship’s center of gravity when v=1.2 m/s (Local)

3.3.2 船舶纵向受力沿程分析

图8(b)为不同流速下,船舶在距离丁坝不同位置处的纵向受力曲线图;纵向力的正负规定为,沿X轴正向为正,反之为负。图8(b)中可以看出,在3个水流流速下,船舶所受到的纵向力沿程变化趋势基本一致。在相同位置处,随着水流流速的增大,船舶受到的纵向力也随之增大。在丁坝前后一倍Lpp范围内(-1

3.3.3 船舶转艏力矩沿程分析

图8(c)为不同流速下,船舶在距离丁坝不同位置处的转艏力矩变化曲线图。转艏力矩的正负规定为沿Z轴正向逆时针旋转为正,反之为负。从图8(c)中可以看出,流速越大,转艏力矩的幅值越大。转艏力矩的变化曲线图有峰值和谷值,且转艏力矩在X/Lpp=-0.5附近由正变负。当水流流速为0.6、0.9 m/s时,转艏力矩峰值出现在X/Lpp=0.25附近处,而在1.2 m/s的水流流速下,转艏力矩的峰值移动到X/Lpp=0.5附近。

从图9(a)中可以看出,由于丁坝的绕流和束窄作用,导致丁坝轴线位置处横向水流流速较高,船舶位于X/Lpp=0.5位置处,船艏部位受到较高的水流冲击,水流流速在此处骤减而产生的低速高压区,使得船舶在此处受到的转艏力矩最大。而船舶位于X/Lpp=-0.5位置处时,船舶左侧水流只能通过船体尾部和丁坝之间的狭窄间隙向下游泄出,水流流速增大,相应的造成船舶左舷压力降低,船舶左、右舷的压强差导致船舶在此处受到的转艏力矩为负。而船舶继续向坝前移动,船舶左右两侧的压强差减少,故而船舶受到的转艏力矩逐渐减少。

3.4 船舶姿态沿程变化分析

3.4.1 船舶横摇分析

图10(a)为不同流速下,船舶的横摇变化曲线图。船舶横摇的正负规定为:船舶左倾为正,右倾为负。从图10(a)中可看出:当水流流速为0.6、0.9 m/s时,船舶的沿程横摇出现了2个明显的正负峰值,分别位于X/Lpp=0和X/Lpp=-0.5。当船舶位于X/Lpp=0处,丁坝前受丁坝阻挡而变成横向流动的水流,再次遇到船体阻挡,只能从丁坝和船体之间的空当向右侧流动,相应的水流流速更大,使得船体左侧的压力降低,船舶产生左倾,图11对船体3个界面的压力分析也可印证这个现象。图11中,船舶位于X/Lpp=0处,船舶的3个不同横截面流场压力图,船舶左侧下方出现明显的低压区,左右两侧的压强差使得船舶出现左倾,表现为横摇出现正向的峰值,这是由于丁坝的束窄作用,船舶左舷和丁坝之间的纵向水流高速发展而产生的低压区,船舶右舷流速相对较低,因此压强较高。同理对船舶位于X/Lpp=-0.5位置处,船体与丁坝之间的空当靠近船体尾部,高速水流产生的低压区也逐渐远离船体,船体左侧附近的水流因船体阻挡而流速降低,压力增大,而船舶右舷水流流速相对较高,因而产生的低压区使得船舶出现右倾,横摇值出现负向的峰值。而水流流速增至1.2 m/s,船舶位于X/Lpp=0.25左右时,丁坝的束窄作用导致坝后横向水流流速较大且高流速区范围较广,此时船舶左右舷水流流速差别不大,导致船舶在此处横摇值出现明显的凹值。

图10 不同流速下船舶横摇、纵摇、升沉沿程变化情况Fig. 10 Changes of the ship’s roll, pitch and heave along the way under different flow rates

3.4.2 船舶纵摇分析

图10(b)为不同流速下,船舶的纵摇变化曲线图。船舶的纵摇正负规定为:沿Y轴正向顺时针旋转为正,逆时针旋转为负。从10(b)中可看出:不同流速下,船舶的横摇沿程变化的规律基本一致,船舶的纵摇峰值出现在X/Lpp=-0.25附近处。随着水流流速的增加,船舶的纵摇峰值出现向坝前移动的趋势,且基本满足随水流流速增大,纵摇值增大的趋势。在靠近丁坝的范围内(-1

3.4.3 船舶升沉分析

图10(c)为不同流速下,船舶的升沉变化曲线图。船舶的升沉正负规定为沿Z轴正向移动为正,反之为负。从图10(c)可看出:在3个水流流速下,船舶的升沉值变化规律基本一致。由于丁坝的束窄作用导致河道的过水断面面积减少,流速越高,丁坝前的壅水量越多,丁坝前水位越高,船舶从丁坝轴线开始向坝前移动,由于水位不断上升导致船舶被抬升;船舶在坝前距离丁坝相同位置处,流速越高,上升值越大。由于丁坝前后的水面比降较大,在这一区域(-1

4 结 论

笔者通过CFD数值计算的方法,对不同水流流速下船舶沿设有丁坝的航道中线行驶时的水动力性能做了数值模拟,得出了船舶在驶过丁坝时的受力、力矩的变化情况,以及船舶的横摇、纵摇、升沉的变化趋势。通过研究得到以下结论:

1)船舶的存在使得丁坝后的流场发生不稳定的横向摆动。在一个摆动周期内,高流速带中先形成高流速核心区,之后又逐渐消失;低流速带中的2个漩涡核心则先融合成一个与高流速核心区上下对应的单一漩涡,之后又分裂为2个漩涡。

2)丁坝轴线附近(-0.5

3)船舶在丁坝轴线位置时会因为靠近丁坝侧的高速低压区而向丁坝发生倾斜,而在丁坝上游(X/Lpp=-0.5)则会因为横向推力的发展,加上丁坝的束窄作用导致远离丁坝侧流速高压强低而远离丁坝方向倾斜。由于丁坝的壅水作用,船舶从丁坝下游开始随水位的升高而被抬高,且流速越大,升沉的幅值越大。