限制水域船舶水动力数值计算研究*
2019-12-27陈明达程细得
陈明达 程细得
(武汉理工大学高性能船舶技术教育部重点实验室1) 武汉 430063) (武汉理工大学交通学院2) 武汉 430063)
0 引 言
当船舶驶入某些通航建筑物内的狭窄航道时，船舶的受力情况与在无限水域中有较大的差别.受阻塞效应影响船舶阻力值会增加，如果航道两侧相对船舶形状不对称，左右两侧流场不对称会导致船舶受横向力与首摇力矩的影响，因此，探究船舶在该类航道中航行时的水动力变化趋势，对通航建筑物的设计、船舶航行的操纵及控制问题具有理论意义和实用价值.
在过去几十年间，针对限制水域船舶水动力的研究有许多学者做了大量的模型实验，数值计算与分析.早年由于计算机能力限制，相关的研究还是主要以基于势流理论的细长体理论和三维面元法为主.Hess[1]针对近岸航行船舶的横向力计算的提出了一套理论模型.熊新民等[2]采用三维Rankine源面元法计算了考虑自由面影响下船舶平行于岸壁航行时的水动力.后来计算机技术愈发成熟，计算机性能提高导致基于黏性流的CFD方法在限制水域船舶水动力研究上得到了广泛应用.桑腾蛟等[3]以KVLCC2船型为对象,利用RANS方法计算其在不对称岸壁及浅水中斜航水动力,研究不对称岸壁、水深,以及其联合作用对船舶斜航水动力以及水动力导数的影响.Hoydonck等[4]对比不同的CFD计算条件分析了自由液面、螺旋桨、流体黏性(粘流与势流理论计算结果对照)和不同船岸距离与水深对船岸效应的影响.
尽管针对限制航道船舶水动力的研究开始较早，但是多数学者研究的还是理论上的规则均匀几何航道，而通航建筑物内的航道形状往往是不均匀的.对此，采用基于RANS(reynolds average navier-stokes)的CFD方法，针对船舶直航驶入船闸，应用流体分析商业软件FLUENT对船舶的黏性绕流场进行数值模拟，并计算船舶的水动力.通过对比计算结果与实验结果验证数值计算的可行性，在此基础上开展数值计算，分析船舶在航道的不同位置处水动力的变化趋势，分析航道形状对船舶进出船闸限制水域水动力的影响规律.
1 数学模型
1.1 控制方程
使用RANS的方法求解NS方程中的连续性方程和动量守恒方程：
(1)
(2)
用张量形式表达三个方向的动量方程，动量方程中速度的物理量都是时均值，带有顶线的物理量是雷诺应力项，船舶流场的马赫数远小于临界值，认为是流体不可压缩的，故方程中省去了密度的偏导项.
1.2 湍流模型
湍流模型选用的是Menter[5]提出的两方程SSTk-ω模型，其混合了k-e与k-ω模型的优势，是较为常用的湍流模型.选择SSTk-ω模型理由有三：①张志荣等[6]认为在常用的湍流模型中，SSTk-ω模型是最适合船舶黏性流场数值模拟的湍流模型；②Lu等[7]在做了多组湍流模型的数值计算对比后之后认为SSTk-ω模型在限制水域船舶水动力计算中更具有优势；③由于在计算的中后半段船舶会驶入以较窄的航道，船舶附近的流场较为混乱复杂，故在近壁处理上，需要使用近壁面模型，SSTk-ω模型对于近壁剪切流计算亲和性较好.为了保证y+=1 ，参考Frank[8]提供的公式得到第一层边界层厚度，使船舶面网格附有5层较密的边界层，边界层的增长率是1.5，以捕捉船舶附近的流场信息.
1.3 自由液面
由于孟庆杰等[9]认为对于船舶在极为窄浅水域中航行时，自由面对船舶水动力的影响较大，所以虽然计算模型的船舶弗劳德数Fr在0.01～0.02，计算时仍考虑自由液面的影响.
VOF的方法就是通过研究单元网格内各种填充介质占网格总体积的比值F(体积分数)来确定自由面[10].VOF方法的追踪的是流体自由液面形状的变化，而不是像势流理论中的自由液面运动边界条件一样追踪液面上质点的运动.
(3)
VOF方法可以追踪复杂的自由液面形状、易实施、易拓展到三维空间.但是VOF只适用于压力基求解器；使用VOF方法的算例的计算域必须充满介质，对于某些高速问题不适用；VOF方法对于在自由液面的捕捉需要高密度的网格支持，变相增加计算压力.
2 数值求解
2.1 计算对象
选取12 000 TEU超巴拿马集装箱船进行计算，船模几何模型侧视图见图1，x方向相对位置坐标的参考点在首垂线处,主要参数见表1.
图1 12 000 TEU几何模型
表1 船模主要尺度
航道岸壁垂直水底，闸室长1.62Lpp,宽度为0.687 5 m.从闸门处向外航道左侧伸出一道长1.32Lpp引墙，右侧则是逐渐外扩，过渡的岸壁的x方向投影长度为0.44Lpp，引航道的宽度为2.725 m，船舶初始位置首垂线距离闸门1.48Lpp.航道几何模型见图2.
图2 航道几何模型
与随船坐标系不同，固定坐标系计算域无来流速度，依靠动网格技术控制船舶做空间移动，故除TOP面需要设置为对称面边界条件(与气相接触，是否有切向物理量变化与计算结果影响不大)，其余壁面均是无滑移固壁(计算域尾段距离船尾初始位置有1倍LPP，理论上计算域尾段对船舶尾流无影响，只需要在该边界上无质量穿透，为了方便计算续性收敛，故选择Wall边界条件).
在航道中心线附近布置运动域，通过interfaces滑移网格与两侧静止与进行数据交换，运动域分成三个部分，船舶近域网格随船舶面网格同步运动，首尾方向的网格随船舶运动渐变，具体分块见图3.
图3 计算域划分与网格重构示意图
刚体运动域与网格渐变域之间由interior面交界,见图4，之间不需要产生网格的差值，与重叠网格相比，需要差值的面网格更少，相对的精度与计算稳定性更高.
2.2 数值方法
数值计算空间离散采用的是有限体积法，空间上采用二阶精度的中心差分格式.针对NS方程的求解压力基求解器，计算过程中对压力进行修正以满足连续性和动量守恒，配合耦合隐式算法Coupled Implicit Solver(适用于非定常计算，全速度范围求解，收敛性能好，但是占计算资源高)实现计算的时间迭代.
2.3 计算工况
Vantorre等[11]在比利时弗兰德水动力研究中心的实验公布的数据，针对其中的Test A与Test B船舶进闸工况的数值模拟，用以验证CFD计算方法的准确性.由于实验条件下未能使船模保持匀速，分析航道形状变化对船舶水动力的影响时需要控制速度变量，故在Test A的基础上设计一组匀速航行模拟Test C，见表2.
表2 三组不同的工况条件
3 计算结果分析
3.1 Test A结果分析
图5 Test A计算结果与实验数据对比图
对比发现阻力的数值计算结果与实验数据虽有一定差值，但是相差较小且总体趋势一致，而横向力与首摇力矩的计算结果差值较大.但是可以发现船舶横向受力的计算结果与实验数据初始分流位置明显(x>-0.3)，此时船模大部分区域驶入引墙范围内，有理由怀疑是船尾的某个变量与引墙相互作用干扰的数值计算的结果，由于数值计算未考虑螺旋桨，而实验条件中船尾配有一定转速的螺旋桨，可能是螺旋桨产生的尾流受引墙反射影响到船体进而影响船舶横向受力.
Kaidi[13]在研究螺旋桨对船模的岸壁效应影响做了多组对照仿真，发现顺时针旋转的螺旋桨对单右侧限制航道船舶水动力的影响较大.当量纲一的量的sbd(ship-bank dinstance)等于0.25时，螺旋桨的转动严重影响船舶去流段与岸壁之间的流场压力分布，使其出现一个明显的高压区，导致“岸推”现象加重.Test A引墙与船舶左舷的sbd为0.2，顺时针旋转的螺旋桨可能使得船舶去流段横向受到一个额外向左的吸力，导致实验横向力结果偏小.为了进一步验证计算方案的准确性，也为了证明当前的误差分析的合理性，故针对Vantorre的Test B做一组数值仿真.
3.2 Test B结果分析
Test B与Test A的主要区别在于Test B的水深只有0.209 m，且Test B进行过程中大部分时间螺旋桨是停止转动的，0转速的螺旋桨相对于船舶主体来说只相当于一个小构件附体，对水动力的影响较小，故Test B的数据结果能代表裸船体在进闸程中的水动力变化.
Test B数值计算结果与实验数据吻合较好，见图6.当x>0.6时船模速度剧烈震荡导致阻力计算结果与实验数据有一定出入，考虑到interface数据交换、网格拓扑能力以及几何模型质量等等不可避免的误差，可以认为该方案的计算结果是符合物理规律的.
Test B实验数据与计算结果的对比与在验证了CFD计算方案的准确性的同时，也在一定程度上能证明Test A的计算误差来源于实验条件中螺旋桨转动扰动流场带来的影响.
图6 Test B计算结果与实验数据对比图
3.3 Test C结果分析
Test A与Test B两组计算证明了CFD仿真方案的准确性，Test C是匀速工况，船模速度取Fr相似条件下的服务航速Vs=0.11 m/s.
3.3.1阻力分析
Test C船模阻力分布图见图7.由图7可知，当船舶位置x<0时阻力变化幅度较小，可以认为小阻塞比时单侧岸壁对裸船体的阻力影响不大(侧壁效应)，从船首部开始进入船闸时(x>0)船舶阻力才开始逐步增加.许立汀中间速度理论的船舶在限制航道内受到的阻力增加，是因为收缩的岸壁使得船舶航道的阻塞率上升，船舶的回流速度增加等于变相增加船舶航速，但是观察阻力计算数据发现，当航道收缩到一定程度，船舶阻力反而减小，说明阻塞率与船舶阻力并不是完全正相关，或者说还存在其他因素影响该航道内的船舶阻力.
图7 Test C船模阻力分布图
由于航道与船舶的相对形状是变化的，采用平均阻塞比的方法定义航道的阻塞效应，平均阻塞率为
(4)
式中：B为船宽；T为船舶吃水；ΔB为船舶首尾垂线平行范围内船岸平均距离；H为航道深度.
为了分析船舶在航行中阻力随阻塞率的增大反而减小的原因，将船舶阻力按成分分离，研究航道形状与阻塞率对船舶阻力的影响，见图8.
图8 航道不同位置阻塞率分布图与Test C船模不同成分阻力分布图
由于船舶的Fr(0.01
限制航道船舶的形状系数不仅与船舶自身形状有关，还与航道相对船舶形状有关.黏压阻力来源于船首尾的压力差，在Test C的剧烈变化区域(x>0)，船舶前半段附近的航道形状与后半段航道形状有明显差别.为了具体分析航道形状对船舶阻力影响，在CFD-POST 中故将船舶在船中处将船舶分割，分别计算船舶前、后半段段船舶阻力的变化趋势，见图9.
图9 Test C前、后半段船模x方向受力大小分布图
由图9可知船舶前半段船模x方向受力变化幅度大于后半段x方向受力变化幅度，即船舶总阻力主要受船舶前半段附近航道形状控制.当x>0，船首开始进入船闸室时，船舶前半段受力开始增加，当x>0.5时，此时船舶前半段完全进入船闸，在之后的航行中船舶前半段左右两侧的航道形状不在发生变化，而此时船模前半段x方向受力开始减小.需要注意的是，由于船闸室的长度(1.622Lpp)是有限的，当船舶逐渐驶入船闸室时，会不可避免的受下级闸门(不可穿透边界条件)的影响.
对于垂直船舶航行方向的壁面可能使船舶阻力的减小原因提一种解释.从能量的角度的看黏压阻力就是船舶单位时间将附进的流体向前推所需要做的功，当x越大，船舶离下级闸门越近，远前端受船舶作用的流体的质量就越小，则船舶对前方的流体单位时间内需要做的功就越少，船舶所受的黏压阻力就会减小.
在0.5
3.3.2横向受力分析
船模横向力初始稳定在-0.5处，当x>-1时，横向力开始增加;当x>-0.75时，横向力开始减小；当船舶完全进入引墙范围时(x>-0.25)，船舶横向力进入第一个波谷；之后横向力开始出现急速的上升；船舶前半段完全驶入闸室时(x>0.5)，横向力组件下降直至调转方向，最后当船舶完全驶入闸室时，横向力趋近于0，见图10.
图10 Test C船模横向力分布图
关于船舶横向力变化趋势分析，在后处理中分别对x=0.26，x=0.36与x=0.46处船模做iso surface分割处理，将船模表面按站位分割成20段，分别求解各个分段的横向力数值，将结果汇总见图11.
图11 Test C船模不同位置处横向力占比示意图
由图11可知，船模表面的横向力会集中在相邻的3～4个站上，且随着船模进入船闸(x增加)，横向力集中位置会相应的后移，且移动规律与船模运动规律反向且一致，说明船模横向力的集中与船模相对航道的某一个固定位置，分析横向力集中位置发现横向力的集中区域出现在闸门与船模交界处，该位置可能是引起船模横向力剧烈变化的关键.
通过分析x=0.46处z=0.185液面x方向速度云图分布(图12)，发现船舶右舷附近流场存在一处高回流速度区域.参考伯努利定理中对速度对压力的影响，切向速度越大法相压强越小，回流集中区域正是导致船舶右舷存在低压区的原因.
图12 x=0.46处z=0.185液面流场速度u云图
船舶驶入闸室的同时闸室内的水需要流出，当船模航行至处于，当船舶航行至x=0.46，闸门处断面系数最大.闸门处左右舷断面积相等，随着流体质点后移右舷断面积逐渐大于左舷，由于流体的连续性，闸室内的流体会倾向于从右舷流出，使得在闸门处右舷单位时间流过的流体体积大于左舷，即右舷的回流速度大于左舷，见图13.
图13 x=0.46处闸门位置剖面速度u云图 (面向x正方向观测)
船舶右舷的低压是由船舶自身形状与航道形状共同作用，当船舶开始驶入闸室时(x>0)，时船舶局部横剖面积变大，闸门处断面系数增大，左右舷回流速度差增大，右舷低压区更加明显，导致船舶向右的横向力增大，当船舶去流段开始于闸门位置接触时，船舶局部横剖面积减小，船舶横向力减小.图14为Test C船模首摇力矩分布图
图14 Test C船模首摇力矩分布图
低压区固定在于闸门处相近的船舶右舷，这意味着船舶的横向力合力位置会随着船舶航行慢慢从船舶前半段向后移动，这就可以解释当x>0.18时船舶首摇力矩开始减小甚至出现反向.
4 结 论
1) 基于黏性CFD理论，研究了某船舶进入巴拿马船闸的航行过程中的水动力进行了研究，对不同的物理条件进行了数值计算，并与实验结果进行比较，验证了本文计算方法合理性，为限制航道船舶操纵提供一定理论依据.
2) 船舶在船闸内航行时，会受左、右、前与底部岸壁的影响，使得船舶水动力发生变化，影响船舶航行的安全。对船舶而言，由于受航道阻塞率影响，船舶在驶入船闸过程中会出现明显的阻力增加现象，但是当船舶逐渐驶入船闸时，下级闸门的存在会使导致船舶阻力降低；对船舶横向受力而言，单侧的岸壁对船舶横向力的影响较小，但是航道形状的剧烈变化会导致船舶横向力的快速增加，且横向力的受力中心会随着船舶航行而后移，从而使得船舶在驶入船闸过程中所受首摇力矩也会出现明显的变化。
3) 仿真是基于文献[14]的基础上，针对其提出的展望，作出的一些计算上的调整，包括对于船模的速度的控制及考虑自由液面对水动力计算的影响，所以在计算精度上也有所提升，但是仿真仍未考虑船舶航行时的浮态变化.限制航道中船舶浮态变化大，对水动力结果有一定的影响，是下一步的研究方向.
