自船舶问世至今，高速船无疑是各类船型当中最引人瞩目的．高速船根据减阻原理的不同，又可分为表面效应船、水翼艇、滑行艇、多体船及复合船型等，其中滑行艇的发展历史是最长久的．这种船型结构简单，阻力性能优异，因此得到了非常广泛的应用．滑行艇高速航行时主要由水动升力支撑艇体，艇体尾倾且向上抬升，此时湿表面积相比于静浮状态显著减少，滑行艇摩擦阻力大幅降低．滑行艇在静水中的运动姿态直接影响了阻力性能．如果滑行艇纵倾角较小，会导致水动升力较小，艇体无法抬升，从而使湿表面积增大，阻力性能不佳；如果滑行艇纵倾角较大，可能导致滑行艇高速航行时产生海豚运动，对滑行艇的安全性十分不利．因此准确预报滑行艇在静水中的运动姿态具有重要意义．预报滑行艇在静水中的运动姿态主要包括三维黏流CFD(计算流体力学)方法和2D+t(二维半)势流方法．三维黏流CFD方法的核心是采用有限体积法求解RANS方程[1-5]．这种方法的优势在于能比较精确地捕捉滑行艇航行过程中的自由面变化，以及流体喷溅、方尾效应等流动分离现象，适用于各种复杂几何形状的艇型，对滑行艇受力预报比较准确．目前黏流CFD商业软件比较成熟，如Fluent，STAR-CCM+等，使用便捷，在工程领域被广泛应用．其缺点在于：为了捕捉精细流场，对VOF捕捉格式及网格质量要求非常高，相对于2D+t势流方法效率较低．2D+t势流方法的思路是将滑行艇直航问题转化为二维剖面入水砰击问题，建立滑行艇静水运动姿态的时域预报模型．在计算二维剖面入水砰击问题中，又可以采用动量定理[6-8]、简化砰击模型[9-11]、边界元[12-15]等不同的方法进行求解，这些方法的预报精度与效率也有所不同．2D+t势流方法计算效率高，缺点在于当将三维问题转化为二维问题时，忽略了首尾部的三维效应，对水动力的计算不够准确．本研究针对2D+t势流方法的不足之处提出了改进措施，并与黏流CFD结果进行了对比与分析．1 基本原理1.1　2D+t方法基本思想当滑行艇以速度U在静水条件下航行时，考虑在初始时刻，位于船前的一个固连于大地的横平面内，追踪该横平面内的流体质点，且此时平面内的流体质点不受艇体的扰动．随着时间推移，滑行艇的各个横剖面依次穿过该横剖面，从固连于大地的横平面来观察流动情况，一个二维物体以速度V=Uτ落入水中(τ为局部纵倾角)，从而将复杂的三维水动力问题简化为二维剖面入水问题，这便是2D+t理论的基本思想．这种方法的关键在于求解二维剖面的入水砰击力．求解二维剖面的入水砰击力是实现2D+t方法预报滑行艇水动力的关键步骤．基于势流理论的简化砰击模型相比于势流边界元方法及黏流CFD方法具有较高的计算效率，相对于动量定理方法具有更高的计算精度，因此将其应用到滑行艇水动力预报中是十分合适的．MLM(modified Logvinovch model)砰击模型[16]是简化砰击模型的一种，对局部斜升角较大的二维剖面砰击入水问题有较好的预报精度，适用于滑行艇的水动力预报．1.2　运动方程的建立基于2D+t理论与MLM砰击模型，可建立滑行艇的时域非线性运动方程．如图1所示：O-XYZ为大地坐标，坐标原点在静水面上；G-x1x2x3为船体坐标系，坐标系原点固连于艇体重心G；x1轴指向船艏方向，x3轴垂直于艇体龙骨线；o-xyz为随船平动坐标系；M为艇体质量；g为重力加速度．船体坐标系纵荡距离η1、垂荡距离η3和纵摇角度响应值η5是相对于随船平动坐标系给定的，其中η1与航速同向为正，η3向上为正，η5首倾为正．10.13245/j.hust.240640.F001图1滑行艇运动的坐标系忽略纵荡运动，只考虑重力及垂直于静水面的水动力N对运动的影响，可建立运动方程Mη¨3=N-Mg；I55η¨5=NLN，(1)式中：I55为艇体纵摇惯性矩；LN为水动力作用中心到重心的纵向距离．N可进一步分解为N=∫f3(t)+fs(t)dx，(2)式中：f3为流体和艇体剖面相互作用的动压力；fs为静压力；t为时间．f3可根据流体的动量变化率得到，有f3=-DDt(maVR)=-DVRDtma+DmaDtVR, (3)式中：ma为剖面附加质量，其与瞬时浸深相关，可通过势流方法或经验公式得到；D为维度；VR为剖面和流体的相对速度，与静水面垂直，剖面向上为正．VR又可以分解成两部分，包括定常滑行产生的垂向速度分量Uη5，以及摇荡运动产生的垂向速度分量．有VR=Uη5+η˙3-η˙5x．(4)由此可求解(DVR/Dt)ma．对于(Dma/Dt)VR这一项，其实际物理意义在于：从固连于大地的横平面看，船体剖面以速度VR入水而产生的二维入水砰击力可以采用MLM砰击模型来进行求解．静压力对剖面的作用也不可忽视．大地坐标系下无限水深条件下的静压力表达式为Ps=-ρgz，(5)式中ρ为水的密度．在瞬时湿表面上对静压力进行积分，即可得到剖面受到的总静压力．至此，艇体横剖面受力的各种成分都已得到，将其代入到运动方程式(1)中就可以求解得到滑行艇在静水中的垂荡与纵摇运动的时域解．1.3　针对2D+t方法固有缺陷的修正办法上述基于2D+t方法的滑行艇水动力预报模型忽略了首尾部流动分离产生的三维效应，对首尾处计算得到的压力与实际的三维情况并不完全相符；此外对静压力的求解也不准确，这些都是2D+t势流方法的固有缺陷，须要有针对性地进行修正．首先是针对艉部压力的改进修正．滑行艇航速很高，往往会产生方尾效应．水流在艇底尾端产生流动分离现象，尾板不与水接触，因此艇底尾端与大气连通，艉部剖面受到的压力也骤降到大气压力．而基于2D+t方法的滑行艇水动力模型由于未考虑方尾效应，艉部剖面处的压力要明显大于大气压．文献[7]依据单体滑行艇模型试验的尾端压力分布结果，提出了经验修正公式，即Ctr=tanh2.5αBFr1s，(6)式中：Ctr为方尾压力修正系数；α为无因次方尾修正系数，可由模型试验结果或者经验确定；B为艇宽；Fr1为宽度弗劳德数；s为横剖面到尾板的纵向距离．艇首也存在有三维效应，2D+t方法计算得到的首部压力偏大，也须要进行修正，这一现象由文献[17]发现．2D+t方法得到的首部峰值点发生在横剖面折角线处刚好产生流动分离的位置，动压力为f3，本研究采用取峰值压力修正系数aF，使动压力限制在fmax=aF f3．除了针对滑行艇首尾部三维效应的改进外，对静压力的修正也是须要考虑的．当滑行艇折角线处发生横向流动分离时，折角线与大气连通，压力也迅速降低至大气压，根据式(5)的求解方法得到的静压力相对于实际情况是偏大的，须要进行修正．设静压力修正系数为abf，则修正后的静压力为PS'，有PS'=abfPS=-abfρgz．(7)2 2D+t势流方法与黏流CFD方法计算结果对比2.1　约束模剖面力(矩)结果对比计算采用的滑行艇模型来自于文献[18]的标模试验．试验中滑行艇模型的特征参数有多种，选取其中一种模型进行数值计算．船模长度为1.143 m，船宽为0.229 m，横向斜升角为10°，质量为7.26 kg，重心纵向位置距离(相对于尾板)为0.457 m．对船模进行约束，约束状态下不同航速对应的船模姿态是采用黏流CFD商业软件STAR-CCM+计算三维自由模[19]确定的，如表1所示，表中Fr2为体积弗劳德数．针对约束模，采用四种方法对船体剖面受到的垂向力及力矩(相对重心)沿船长方向的10.13245/j.hust.240640.T001表1滑行艇计算工况Fr2纵倾角/(°)升沉/船宽2.8913.8290.0593.6143.5960.0914.3373.0300.110分布情况进行了计算．这四种方法包括：a．未考虑三维效应和静压力修正的原始2D+t势流方法；b．考虑三维效应和静压力修正的改进2D+t势流方法，且不同速度工况下的修正系数由表2给出；10.13245/j.hust.240640.T002表2改进2D+t势流方法修正系数Fr2αaFabf2.8910.20.60.83.6140.20.50.84.3370.20.50.8c．先采用黏流CFD方法求解二维滑行艇剖面匀速入水载荷，再将剖面载荷随时间的变化情况转化为三维滑行艇剖面载荷在船长方向的分布情况，即2D+t黏流CFD方法；d．三维黏流CFD方法．图2为当Fr2=2.891时，横剖面垂向力对比和横剖面纵倾力矩对比，图中：f1为横剖面垂向力；I1为横剖面纵倾力矩．这里以图2(a)为例，详细说明改进2D+t势流方法修正系数的确定方法．不同修正系数的作用范围有明显区别．从式(6)可知方尾修正系数Ctr随横剖面纵向位置(距船尾板)增大而迅速趋近于1.0，因此方尾修正只在靠近船尾区域(0～0.2 m)才会发挥主要作用．α决定船尾区域压力沿纵向变化的快慢．aF只影响首部区域(0.7～0.9 m)，对首部的压力峰值做截断处理．在船中区域(0.2～0.7 m)，横剖面压力基本只受到静压力修正系数的影响．由于不同修正系数的主要作用区域不同，因此可以调节三种修正系数，直到船尾、船中、船首区域的改进2D+t方法的压力分布结果都能够与三维黏流CFD结果符合．10.13245/j.hust.240640.F002图2横剖面垂向力对比和横剖面纵倾力矩对比(Fr2=2.891)船体剖面垂向力及力矩(相对重心)沿船长方向的分布情况如图3和4所示．原始2D+t势流方法和2D+t黏流CFD方法的结果十分接近．相比于三维黏流CFD的结果，这两种方法在滑行艇的首部会产生明显的压力峰值，在尾部的压力也较大，没有骤降到大气压，不符合方尾效应的实际情况．出现这种现象的原因是2D+t方法将三维滑行艇直航问题转化成二维剖面入水问题的过程中，忽略了首尾部的三维效应；而且由于首尾部压力到重心的力臂较长，压力的误差会造成首尾部相对重心的纵倾力矩存在较大的计算误差．此外，原始2D+t势流方法得到的艇体中后部(横剖面纵向位置小于0.5)的压力要略大于2D+t黏流CFD方法，这是因为2D+t黏流CFD方法考虑到了船体剖面横向的流动分离造成的折角线处静压力损失，结果更加合理．10.13245/j.hust.240640.F003图3横剖面垂向力对比和横剖面纵倾力矩对比(Fr2=3.614)10.13245/j.hust.240640.F004图4剖面垂向力对比和剖面纵倾力矩对比(Fr2=4.337)改进2D+t方法则考虑了三维效应对压力分布的影响，在首部对动压力进行了截断处理，在尾部采用了经验公式进行修正，使尾部压力能快速降低到大气压力，此外静压力也做了修正，因此剖面压力与纵倾力矩都能与三维黏流CFD结果符合得较好．2.2　自由模运动姿态结果对比约束模计算中通过压力曲线逼近得到的三种修正系数被继承到自由模计算中，作为改进2D+t方法的已知参数．根据以上的修正措施，开放自由度，可得到滑行艇运动姿态时域预报结果．将2D+t势流方法计算得到的滑行艇运动姿态结果与模型试验结果及三维黏流CFD结果进行对比，如图5和6所示．由于改进2D+t方法的修正系数是通过与三维黏流CFD的压力分布结果来确定的，因此改进2D+t方法与三维黏流CFD方法得到的自由模运动姿态结果基本符合．改进2D+t势流方法得到的纵倾角要明显大于原始2D+t势流方法，主要原因是原始2D+t势流方法对首尾部力矩的计算存在较大的误差．纵倾角直接影响到了求解二维剖面入水砰击问题中的入水速度，进而影响了滑行艇动压力的计算结果．动压力的预报精度关系到艇体能否向上抬升，因此改进2D+t势流方法得到的升沉结果要略大于原始2D+t势流方法，与模型试验结果比较接近．综上所述，改进2D+t势流方法对滑行艇在静水中运动姿态的预报精度要明显优于原始2D+t势流方法．10.13245/j.hust.240640.F005图5自由模纵倾角结果对比10.13245/j.hust.240640.F006图6自由模升沉结果对比3 结论a．对于约束模，改进2D+t势流方法可以明显改善滑行艇首尾部的压力和力矩，弥补2D+t方法的固有缺陷．b．对于自由模，改进2D+t势流方法对升沉运动和纵倾角的预报结果相对于原始2D+t方法有了明显提高，与模型试验结果十分接近，验证了改进2D+t势流方法的有效性．