传统势流理论满足微幅波浪线性假设，并且仅考虑静水刚度的线性部分，忽略不同自由度之间耦合作用．随着海况的增加，一方面半潜式平台波浪载荷和运动响应的非线性特征不断增强，另一方面横撑或者浮筒出水会引起流体载荷的急剧变化．国内外学者对海上浮式结构物的非线性刚度开展了一系列研究工作．Al-Solihat等[1]对单柱式浮式结构物在大转角下的静水力载荷进行理论推导，并得出矩形柱体和圆形柱体非线性回复力矩的计算公式．Ullah等[2]针对三种构型的浮式风机，计算水线面变化引起的非线性回复刚度并进行时域模拟，结果表明考虑非线性刚度后纵摇和垂荡响应可大幅降低．Wei[3]等开展半潜式平台在规则波条件下的模型试验和数值模拟，研究发现当波浪激励接近垂荡固有周期时平台会发生参数激励纵摇．Kvaleid等[4]针对GG5000半潜式平台开展时域水动力分析来考虑瞬时湿表面变化对运动的影响，并指出垂荡和纵摇运动耦合所引起的马休不稳性及随之出现的纵摇在入射波浪下的异相运动是导致该平台气隙显著减小的主要原因．Wei等[5]研究横撑出水和耦合运动对共振运动的影响，发现横撑出水会导致共振频率增加，而耦合运动则会改变响应谱的峰值．目前海洋工程领域主流采用线性势流理，非线性回复力矩的研究相对较少，并且缺乏系统的理论分析和数值对比，本研究针对ExWave半潜式平台，通过理论研究、模型试验和数值模拟来综合分析非线性回复力矩对纵摇运动的影响．1 模型试验与频域分析结果对比ExWave平台浮体几何构型见图1，图中MP，MC和MB对应浮筒、立柱和横撑的Morison单元，M1~M9表示单元分段情况，浮筒长107.5 m，宽14.3 m，高9.5 m，立柱截面长宽均为12.5 m，横撑直径5.0 m．选定工况吃水23.0 m，横摇和纵摇的固有周期分别为52.1 s和56.5 s．该平台的水池试验在挪威Sintef海洋研究院水池完成，模型缩尺比为1∶50．不规则波采用Torsethaugen谱简化形式[6]，操作和自存工况下的有义波高分别为5.5 m和15.0 m，谱峰周期分别为8.5 s和 16.0 s．试验测得迎浪条件下纵摇统计结果见表1，滤波截止频率为0.1 rad/s．10.13245/j.hust.238950.F001图1平台及其Morison单元示意图10.13245/j.hust.238950.T001表1模型试验结果波浪条件低频纵摇波频纵摇纵摇总值方差/(°)²最大值/(°)方差/(°)²最大值/(°)方差/(°)²最大值/(°)操作2.074.880.161.952.265.66自存1.684.632.145.593.878.34低频纵摇运动预报的关键在于黏性阻尼的确定．OTG-13[7]推荐采用4%临界阻尼系数，自由衰减试验[8]测得该平台纵摇的临界阻尼系数为4.4%．Gao[9]和Pan等[10]研究表明，高海况下半潜式平台的阻尼可能被严重低估．本研究构建Morison单元并进行随机线性化来获得等效线性黏性阻尼，曳力系数参考Ma等[11]研究成果．表2为临界阻尼系数统计结果，可以看出黏性阻尼随着海况的增加急剧增大．基于线性化黏性阻尼，通过二阶势流计算得到的二阶差频纵摇传递函数R2-(μ,μ+ω)见图2，图中：ω为差频值；μ为波浪频率．可以看出：当差频接近纵摇固有频率时，后者结果明显小于前者．10.13245/j.hust.238950.T002表2临界阻尼系数波浪条件垂荡横摇纵摇操作海况2.15.85.9自存海况9.326.025.0%10.13245/j.hust.238950.F002图2二阶纵摇传递函数参考OTG-13低频响应频域计算方法得到的纵摇结果见表3．对比表1和3可知，自存海况下通过频域分析得到的方差和最大值均明显高于试验值．频域方法得到自存和操作海况下方差的比值为2.77，而试验测得的比值为0.81，这揭示了高海况条件下半潜式平台低频纵摇受到一定的抑制作用．频域分析方法无法考虑二阶以上成分和其他非线性因素，因此本研究通过理论推导和时域分析研究非线性回复力矩对低频纵摇运动的影响．10.13245/j.hust.238950.T003表3频域分析统计结果波浪条件方差/(°)²最大值/(°)操作海况1.135.67自存海况3.159.452 基于势流理论的纵摇回复力矩定义平台总体坐标系GCS：O-XYZ，以及随体坐标系BCS：o-xyz．平台六自由度运动X={η1,η2，η3,ζ1,ζ2,ζ3}T，转动欧拉角顺序为ζ1,ζ2,ζ3．BCS中任意一点位置向量为r，该点在GCS中的位置向量为rG=r0+η+Rr，其中：r0为初始BCS原点在GCS中的位置向量；η={η1,η2,η3}T为物体平移；ζ={ζ1,ζ2,ζ3}T为物体转动，旋转矩阵为R=R11R12R13R21R22R23R31R32R33=c2c3-c2s3s2c1s3+s1s2c3c1c3-s1s2s3-s1c2-c1s2c3+s1s3c1s2s3+s1c3c1c2,si=sin ζi,ci=cos ζi    (i=1,2,3)．浮体在静水中的纵摇回复力矩可表示为My=My1+My2；My1=-ρg[Vη1-Axfη3-L12ζ1+(Vzb+L11)ζ2-Vybζ3];My2=ρg[(η1ζ1-η3ζ3)Ayf-η1ζ2Axf-η1ζ2AyfZo+12ζ12+32ζ22AxfZo+η1η3A-η3ζ2AZo+ζ1ζ3(L11-L22)+2ζ2ζ3L12+12(ζ22+ζ32)Vxb],式中：ρ为水密度；g为重力加速度；(xb,yb,zb)为浮心；(xf,yf,zf)为漂心；A为水线面面积；V为排水体积；Lij为水线面面积矩，Lij=∬Axyds；Zo为BCS原点距离水线面高度；式中各项准确至二阶，上标(1)和(2)分别为一阶和二阶项．若BCS原点与平台浮心在同一垂线上，则有xb=yb=0．对于Exwave平台，xf=yf=0，L12=L21=0，静水纵摇回复力矩可整理为My=-ρg[-Vη1+(Vzb+L11)ζ2]+ρg[η1η3A-η3ζ2AZo+ζ1ζ3(L11-L22)].垂荡与横荡和纵荡之间耦合项的影响可以忽略．若初始BCS原点位于静水面，则有Zo=0．L11与L22非常接近，因此艏摇与横摇和纵摇的耦合项的数值很小．综上所述，静水纵摇回复力矩中的二阶项可以忽略．当浮体发生大幅运动时，还须考虑其他非线性因素．3 大转角纵摇回复力矩大转角条件下，回复力矩的成因可归为浮心和重心位置变化MBGG、立柱水线面变化MwpG、横撑出水MbrcG、浮筒出水MptnG四类，试验条件下浮筒未发生出水，因此本研究忽略该项．势流理论考虑MBGG和MwpG中的一阶和二阶部分，这里进一步分析MbrcG和MptnG及其占比随转角的变化趋势．3.1　立柱水线面改变图3为立柱水线面改变引起的回复力矩，B为浮心．立柱坐标系LCS：o'x'y'z'下，水线面与立柱交线上任意一点的高度可表示为h(x', y')=-[R(rt+ro,i+{x', y', z'}T)]z'N∙n，式中：rt={0,0,η3}T为浮体平移；ro,i={xo,i,yo,i，zo,i}T为在LCS原点在BCS中的位置向量；N和n分别为GCS和BCS的垂向单位向量，且N∙n=R33．从而上式可整理为：h(x',y')=-R31R33x'-R32R33y'-R31xo,i+R32yo,i+R33η3R33=ax'+by'+c;a=-R31R33;b=-R32R33;c=-R31L+R32B+R33η3R33.10.13245/j.hust.238950.F003图3立柱水线面变化引起的回复力矩水线面变化导致单根立柱排水体积的变化量为Vi'=∬Ac,i'g(x', y')dA=cAc,i' ,式中Ac,i'为立柱投影面积．rb,i'={xb,i', yb,i', zb,i'}T为Vi'浮心在LCS中的位置向量，可以表示为{xb,i',yb,i',zb,i'}=∬Ac,ih(x', y'){x', y', z'}dAV=aI'cAc', bI'cAc', (a2+b2)I'2cAc'+c2,式中I'为立柱投影面积的面积矩，由对称性可知I'=Ixx'=Iyy'，Ixy'=0．浮心在BCS和GCS下的位置向量分别为rb,i=ro,i+rb,i'和Rrb,i．浮力向量FB,i={0, 0, ρgcAc'}T，其引起的回复力矩Mwp,iG=Rrb,i×FB,i．总的回复力/力矩为FwpG=∑FB,i和MwpG=∑Mwp,iG．该平台∑x=∑y=∑xy=0，整理回复力和力矩公式为：Fwp,zGMwp,xGMwp,yG=ρga33000a44a450a54a55AI11I22;a33=-η3；a44=s13c3-2s1c3+c13s2s3-c1s2s32c12；a45=-s1s22c3-c1s23s3+2c1s2s32c12c22；a54=s13s3-2s1s3-c13s2c3+c1s2c32c12；a55=-s1s22s3+c1s23c3-2c1s2c32c12c22．图4为立柱水线面变化引起回复力矩与线性部分之比，r44和 r55分别定义为：r44=-Mwp,xGρgI11ζ1;r55=-Mwp,yGρgI22ζ2．10.13245/j.hust.238950.F004图4立柱回复力矩与线性值之比3.2　横撑出水横撑圆截面出水回复力矩示意图见图5，图中：di为圆心o'至水线面距离；Ri为圆截面半径．当-Ri≤di≤Ri时，出水面积Ai(di)和出水面积形心B至圆心的距离Li(di)可分别表示为：Ai(di)=πRi2-acos(di/Ri)Ri2+diRi2-di2；Li(di)=-2(Ri2-di2)3/23Ai(di)．由几何关系可知di=Zi/cos α，其中cos α=J∙j=R22，J和j分别为GCS和BCS下纵向单位向量，Zi为单元圆心距离水线面的瞬态高度．Zi的计算公式为Zi=-xiζ2+yiζ1+zini,i+η3-Si，式中：zini,i为静平衡时圆心距离水线面初始高度；(xi,yi,zi)为圆心o'在BCS下坐标；Si为各单元圆心处的瞬态波面升高．rb,i'可表示为rb,i'=cos ζ20sin ζ2010-sin ζ20cos ζ2{0, 0,  Li(di)}T．10.13245/j.hust.238950.F005图5横撑圆截面出水回复力矩示意图横撑出水浮力FB,i={0, 0, -ρgAi(di)li}T，其中li为圆截面单元长度，进而求得相应的回复力矩Mbrc,iG．通过积分可获得横撑出水引起总的回复力FbrcG和力矩MbrcG，横撑出水引起的纵摇回复力矩见图6，图7为纵摇回复力矩的不同成分随转角的变化．不难看出，横撑出水是纵摇回复力矩呈现非线性变化的主要原因．10.13245/j.hust.238950.F006图6横撑出水引起的纵摇回复力矩10.13245/j.hust.238950.F007图7回复力矩随转角的变换4 横撑出水对纵摇运动的影响4.1　横撑出水的载荷特征以自存海况下模型试验结果为研究对象，根据已建立的公式，计算得到前横撑b1与后横撑b2出水回复力矩Mb1G和Mb2G的时间历程如图8所示．此试验条件下的纵摇极值发生在2 614.4 s，选取2 000~3 000 s时间段内数据进行分析．图9为该时间段内平台的纵摇运动时历ζ2tot、滤波后的波频ζ2WF和低频ζ2LF成分对比，其中低频部分占主导地位．图10为横撑出水和二阶波浪激励纵摇力矩(My(2))的对比，前后横撑出水诱导回复力矩与波浪激励二阶纵摇力矩处于同一量级，其中横撑出水回复力矩持续时间短但变化速度大．10.13245/j.hust.238950.F008图8横撑出水回复力矩时程10.13245/j.hust.238950.F009图9纵摇不同成分对比10.13245/j.hust.238950.F010图10横撑出水和二阶波浪激励纵摇力矩的对比仅考虑可引起横撑最大回复力矩2%以上的出水事件，3 h内前后横撑分别发生136次和57次出水，平均出水间隔分别为79.4 s和189.5 s．为研究横撑出水的特征，继续对前后横撑在不同出水时长Td、不同出水间隔Ts和不同回复力矩比率M¯yG条件下出水次数进行分析．回复力矩比率定义为M¯yG=|MyG|/(ρgVLb)，式中：V为单个横撑排水体积；Lb为横撑中心线距离平台中横剖面的距离．图11为不同条件下前后横撑出水发生次数，可以看出：前后横撑出水时长均集中于1~5 s，但出水间隔分布较广；前横撑最频繁的出水间隔为10 s左右，但总数只有14次；前后横撑回复力矩比率主要分布于40%以下，出水次数最大值出现在5%以下．Kvaleid等[4]和Wei等[5]针对各自研究的半潜式平台均发现了横撑和浮筒入水会引发参数激励纵摇，其显著特征为纵摇周期接近两倍垂荡固有周期，该现象进一步导致了横撑垂向运动幅度增大和出水频率增加．从横撑出水时长和出水间隔判断，该平台的运动特征与上述两平台存在明显不同，它并未发生参数激励纵摇．该结论可从图12的纵摇响应谱得到验证，其最大值发生在纵摇固有频率ωn(ζ2)附近，且远离垂荡固有频率ωn(η3)及谱峰频率ωp的一半．运动特征差异主要是由于横撑与静水面距离不同造成的，本平台该值达到8.00 m，而上述两平台只有2.95 m和1.95 m．该平台的横撑离水线面较远，出水概率更低，一般认为偶尔出水并不会引起垂荡与纵摇之间持续的耦合效应；另一方面，纵摇固有周期相对较长而横撑出水时间极短，并且其引起的载荷变化相较于一阶波浪载荷是小量，因此其对波频纵摇的影响非常有限．综合上述分析，横撑出水呈现非连续长间隔的载荷特征，并且与低频运动密切相关．10.13245/j.hust.238950.F011图11不同条件下前后横撑出水发生次数10.13245/j.hust.238950.F012图12纵摇响应谱4.2　横撑出水对纵摇运动的影响通过时域仿真研究横撑出水对纵摇运动的影响，该模型中浮体参数和系泊设置等与模型试验[8]保持一致．时域运动方程参考文献[11]，本研究仅考虑阻尼项线性部分，忽略Morison单元非线性激励力和阻尼力，载荷项表示为F(t)=Fw(1)(t)+Fw(2)(t)+Fm+Fb，式中：Fw(1)(t)为一阶波浪激励力；Fw(2)(t)为二阶波浪激励力；Fm为系泊系统回复力；Fb为横撑出水诱导载荷．采用B0和B1两种数值模型，分别忽略和考虑Fb．图13和14对比两种模型在10个随机数下纵摇运动的极值ζ2max和方差ζ2var．对比ζ2max可以看出：横撑出水主要影响纵摇低频部分，对波频部分也具有一定的作用，后者可能由于前后横撑多次连续出水造成．考虑横撑出水作用后，纵摇的正向极值增加而负向极值减小，这说明横撑出水对纵摇运动既有抑制作用又有激励作用，同时反映了高海况条件下考虑横撑出水的必要性．对比ζ2var可知：横撑出水主要影响低频部分，但作用比较有限．10.13245/j.hust.238950.F013图13不同随机种子数下纵摇极值对比10.13245/j.hust.238950.F014图14不同随机种子数下纵摇方差对比假定不同种子数下的极值符合耿贝尔(Gumbel)分布，整理得到纵摇的时域分析结果见表4，数据显示横撑出水导致低频纵摇增大，与预期结果相反．该时域计算结果低于表3中的频域分析结果，但仍然高于表1中的模型试验结果，这表明仅仅补充回复力矩非线性效应并不能解决纵摇数值计算结果偏保守的问题．一方面，当前建立数值模型未考虑波浪的非线性及平台的辐射/绕射作用，也未考虑横撑出水和入水砰击载荷、黏性载荷和附加质量等变化；更重要的是，在极限海况下非线性因素对平台运动的影响逐渐凸显，例如细长杆构件载荷非线性、浮筒波浪力随垂向位置改变等，数值模型未计入其他非线性因素．与其他文献相比，Wei等[5]在恶劣海况下的数值模拟与模型试验结果同样存在较大差异，Kvaleid等[4]尝试同时调整线性阻尼和Morison单元参数来实现模型试验与数值模拟的响应谱保持一致，但未给出计算结果的对比．为了解决纵摇时域模拟结果偏保守的问题，在考虑横撑出水的同时，还须在数值模型中兼顾其他非线性效应．10.13245/j.hust.238950.T004表4时域分析统计结果数值模型低频纵摇波频纵摇纵摇总值方差/(°)²最大值/(°)方差/(°)²最大值/(°)方差/(°)²最大值/(°)B02.757.162.567.005.5611.27B13.078.252.477.005.8211.865 结语理论推导和数值计算结果表明，横撑出水是半潜式平台在大幅运动时非线性回复力矩的主要成因．在建立非线性回复力矩计算方法的基础上，通过开展时域分析研究半潜式平台在恶劣海况下回复力矩的载荷特征及其对纵摇运动的影响．计算结果显示，所研究平台的横撑出水呈现非连续长间隔的特征并与低频运动具有很高的相关性．对比考虑和忽略横撑出水两种数值模型的纵摇统计结果可知：横撑出水主要影响纵摇的低频部分，并且导致正向极值增加而负向极值减小．自存海况下，改进时域方法计算得到的低频纵摇比频域方法减小15%，结果更接近于试验值，但两者之间仍存在一定差异，因此时域模型中还应进一步考虑其他非线性因素．