船体波浪载荷的预报贯穿于整个船舶设计期，是结构校核的根本依据，也是船舶航行安全的重要保障；因此，波浪载荷预报中应考虑船体整个寿命期内可能遭受的极端海浪环境及响应．近年来，在节能减排、绿色航运的号召下，船舶大型化成为必然发展趋势．大型船舶的刚度较低、航速较高、船艏常具有大外飘结构，在高海况航行中，大幅船波相对运动会使船体受到波浪的强烈冲击，结构易产生“颤振”现象[1]，对船体的极限承载能力提出了更高要求．2007年MSC Napoli号集装箱船破坏事故后，国际船级社协会(IACS)提出了《集装箱船总纵强度准则》(UR S11A)[2]统一要求，明确要求考虑砰击颤振对船体载荷的影响．中国船级社(CCS)也发布了《波激振动和砰击颤振对船体结构疲劳强度影响计算指南》[3]，2018年进一步修订增加了极限强度评估，完成了《船体结构波激振动和砰击颤振直接计算评估指南》[1]的编写．然而UR S11A的设计载荷仍是基于经验公式，无法反映真实船型信息及结构特征．而《船体结构波激振动和砰击颤振直接计算评估指南》仅给出了非线性载荷的计算要求及流程，并没有规定算法及计算工具，要精准预报船体载荷仍须结合数值及试验手段．近年来，三维时域水弹性方法被广泛应用于船舶与海洋工程的流固耦合分析[4-8]．然而，现有数值计算多为弱非线性或部分非线性方法，不能反映船体的砰击力．文献[9]在弱非线性时域方法的基础上，采用二维Wagner模型计算砰击压力，预报了船舶的砰击颤振响应，数值结果与试验的运动及垂向弯矩符合较好，但部分工况的压力差异较明显．文献[10]应用间接时域法，采用动量定理计算船体的砰击力，预报了不规则波下船舶的运动及波浪载荷，其中垂向运动结果的超越概率分布与试验结果符合较好，但高海况下垂向弯矩存在一定差异．本研究提出了有航速船舶的全非线性时域水弹性数值方法，采用铁木辛柯梁理论分析船体的模态，基于全非线性边界条件，在时域中更新船体湿表面及自由液面，模拟船艏艉的出入水现象，能够更加真实模拟船体航行状态，精准预报船体的波浪载荷．基于超大型集装箱船经济效益及远洋恶劣环境的考虑，本研究提出了分别预报极限工况(对应5 kn低航速)与服务工况(对应经济航速18 kn)设计载荷的方法，预报两种航速下1×10-8概率水平的长期值．根据不同波高的全非线性水弹性数值结果及波浪载荷长期值，确定非线性设计波参数及载荷水平，结合UR S11A规范与分段模型试验结果，综合探究船体的设计载荷．1 规范公式根据UR S11A，对于满足其适用范围的集装箱船，沿船长分布的船体垂向波浪弯矩中拱值MW-Hog与中垂值MW-Sag可分别由下式计算MW-Hog=+1.5fRL3CCW(B/L)0.8fNL-Hog；(1)MW-Sog=-1.5fRL3CCW(B/L)0.8fNL-Sag，(2)式中：fR=0.85；L为船的长度；C为波浪参数；CW为船体结构吃水处的水线面系数；B为船的宽度；fNL-Hog与fNL-Sag分别为中拱与中垂状态的非线性修正系数．波浪参数C可由下式计算C=1-1.50(1-L/Lref)2.2    (L≤Lref);1-0.45(1-L/Lref)2.2    (LLref),(3)式中Lref为参考长度，Lref=315CW-1.3；(4)CW=AW/(BL)，(5)其中AW为吃水处的水线面面积．非线性修正系数fNL-Hog与fNL-Sag计算如下fNL-Hog=0.3(CB/CW)T；(6)fNL-Sag=4.5[(1+0.2fBow)/(CWCBL0.3)]，(7)式中：CB为结构吃水处的方形系数；T为结构吃水；fBow为艏部外飘相关的形状系数．又有fBow=(ADK-AWL)/(0.2Lzf)，(8)式中：ADK为最上层甲板(包括艏楼甲板)在水平面的投影面积，从0.8L开始向前到船艏的部分都计入；AWL为结构吃水处的水线面面积，范围也是从0.8L到船艏部；zf为在艏垂线处最上层甲板(包括艏楼甲板)到水线的垂向高度．2 数值方法2.1　线性波浪载荷的长期预报工程中常采用长期预报极值的方法确定船舶营运期间所受到波浪载荷的最大值．通常将长期预报作为一系列短期平稳随机过程的组合来处理．对于线性载荷，其波浪及载荷响应均满足瑞利分布，载荷幅值的概率密度函数及概率分布函数可表达为：f(x)=(2x/E)exp(-x2/E)；(9)F(x)=1-exp(-x2/E)，(10)式中E为方差的2倍，E=2m(Hs,Tz,β,V)．对于线性响应，船体对波浪的响应视为线性时不变系统，有    Sw(ω,Hs,Tz,V,β+θ)=H2(ω,V,β+θ)Sζ(ω,Hs,Tz,θ),(11)式中：Sw(ω,Hs,Tz,V,β+θ)为波浪载荷的谱密度函数；H(ω,V,β+θ)为幅频响应函数的模，指单位波幅下的载荷幅值；Sζ(ω,Hs,Tz,θ)为海浪谱密度函数；β为浪向角；θ为组合波与主浪向的夹角．根据谱密度函数的性质，方差可以表示为    m(Hs,Tz,β,V)=∫-π/2π/2∫0∞Sw(ω,Hs,Tz,V,β+θ)dωdθ=∫-π/2π/2∫0∞H2(ω,V,β+θ)Sζ(ω,Hs,Tz,θ)dωdθ.(12)一旦求得方差，式(10)的概率分布函数即可确定．那么载荷幅值大于某一定值的超越概率为    Q(x)=P(X≥x)=∑i∑j∑kpi(Hs,Tz)pj(β)pk(V)exp(-x2/E).(13)通常规定船舶航行25 a能够遭遇的总载荷循环次数n=1×108，长期预报超越概率Q=1/n=1×10-8．只要确定了船舶的航行区域对应的海况及概率水平，即可确定对应的载荷特征最大值Xmax．2.2　全非线性时域水弹性数值方法2.2.1　流场定义及初边值条件为了真实描述船舶在流域中的航行，引入三个笛卡尔右手坐标系，即固定坐标系、随船平动坐标系及固船坐标系，其二维示意图如图1所示，固定坐标系固定在水面上，随船平动坐标系以航速U0随船体行进，固船坐标系原点固结在船体重心处且随船体一起摇荡运动．10.13245/j.hust.241213.F001图1流场及坐标系二维示意图对于无黏、无旋、不可压缩的理想流体，速度势满足控制方程∇2ϕ=0．(14)将流场中的总速度势分离为入射波势ϕI及扰动势ϕD，可得：ϕ=ϕI+ϕD;ζ=ζI+ζD,(15)其中ϕD也满足控制方程．要确定其唯一解，须要给出扰动速度势在物面、自由水面、远方边界等满足的条件．由于船体表面是不可渗透的，在物面上扰动速度势满足    ∂ϕD/∂n=(U0+η˙1,2,3+η˙4,5,6×rb+u˙)∙n-∂ϕI/∂n,(16)式中：η为船体六自由度位移矢量；u为结构弹性变形．在随船平动坐标系中，自由面的运动学条件和动力学条件可表示为：    ∂ζ/∂t=∂ϕ/∂z-(∂ζ/∂x)∙(∂ϕ/∂x-U0)-(∂ζ/∂y)(∂ϕ/∂y);∂ζD/∂t=∂ζ/∂t-∂ζI/∂t;(17)∂ϕD/∂t=-(∇ϕD∙∇ϕD/2)-gζD+(∂ϕD/∂x)U0;δϕD/δt=∂ϕD/∂t+(∂ϕD/∂z)(∂ζ/∂t),(18)式中δ/δt为物质导数．在无穷远处，一般认为扰动速度势的影响为零，则∇ϕD→0．(19)Rankine源方法不满足辐射条件，本研究在自由面远端布置人工阻尼层吸收波浪，消除反射波的影响，则边界条件可表示为∂ζD/∂t=g(ϕD,ζD,t)-μ(r)ζD；(20)∂ϕD/∂t=f(ϕD,ζD,t)-μ(r)ϕD，(21)式中：μ(r)为阻尼强度；r为计算点到设置阻尼区起始处的距离[11]．假定初始自由面是静止的，则ϕD=0,ζD=0 ．(22)根据格林第二定理，速度势的控制方程可以通过以下边界积分方程来求解    A(p)ϕD(p)=∫{G(p,q)[∂ϕD(q)/∂nq]-ϕD(q)[∂G(p,q)/∂nq]}dsq,(23)式中：A(p)为源点p处的立体角；格林函数G=1/Rpq，Rpq为源点q和场点p之间的距离．在初边值条件已知后，通过数值方法求解该边界积分方程，即可得到船体上各点的扰动速度势[12]．2.2.2　广义运动方程及其解根据三维有限元法可将船体结构离散为有限个自由度的系统，其结构动力学方程为MU¨+CU˙+KU=P+F+G，(24)式中：M，C，K分别为系统的总质量矩阵、总阻尼矩阵和总刚度矩阵；U为系统总节点位移列阵；P，F，G分别为结构所承受的分布外力、集中力和体积力的等效节点力列阵．节点位移U=Dp=∑r=1mDrpr(t)，(25)式中：D为系统的振形，本研究将船体视为两端自由的非均匀铁木辛柯梁，通过迁移矩阵法计算得到船体固有频率及振形［13］；p={p1(t),p2(t),…，pm(t)}T为广义主坐标列阵，pr(t)(r=1,2,…，m)对应其第r阶模态的主坐标分量．将广义运动微分方程两端左乘DT，并将上式代入，可以得出结构离散系统的主坐标运动方程式ap¨+bp˙+cp=Z+Δ+Q，(26)式中：a,b,c分别为结构广义质量矩阵、广义阻尼矩阵、广义刚度矩阵；Z,Δ,Q分别为广义表面分布力、广义集中力和广义体积力列阵．根据伯努利方程，作用于船体瞬时湿表面上的流体压力可表达为p=-ρ(ϕt+∇ϕ∇ϕ/2+gz)．(27)而作用于物体的第r阶广义流体力可以写为Zr=-∬Sbp(n⋅ur)ds    (r=1,2,…,m)．(28)通常船舶结构系统所承受的体积力为重力，广义重力可以写为Qr=-∭Vρbgwrdv    (r=1,2,…,m)，(29)式中：ρb为结构质量密度；V为结构总体积．本研究中模型属无约束动力系统，无广义集中力作用．已知广义力，式(26)的主坐标即可通过数值手段求解，本研究采用Newmark-β法进行迭代求解[14]．根据模态叠加原理，可得到船体结构的位移w(x,y,z,t)、弯矩M(x,y,z,t)和剪力V(x,y,z,t)：w(x,y,z,t)=∑r=7mwr(x,y,z)pr(t);M(x,y,z,t)=∑r=7mMr(x,y,z)pr(t);V(x,y,z,t)=∑r=7mVr(x,y,z)pr(t).(30)当船体中部处于波峰上，其重力和浮力在纵向分布不对称，船舯浮力大于重力、艏艉重力大于浮力，这种甲板受拉、底部受压、船体中部拱起的纵向弯曲状态称为中拱，反之若船舶处于波谷上，船舯重力大于浮力、艏艉浮力大于重力引起的船体中部下垂的弯曲状态称为中垂．2.2.3　数值结果验证本研究以一艘10 000TEU集装箱船为分析船型[15]，结合分段模型试验开展对数值算法的验证分析．算法自身的收敛性分析和低工况对比验证结果已公开发表[16]．图2展示了航速5 kn下波高4，6，9，12，15，17 m，与航速18 kn下波高4，6，9，12 m工况，船舯垂向弯矩的数值及试验结果随波高的变化趋势，本研究数值结果与试验结果在趋势及数值上符合较好，且曲线反映出中拱弯矩与中垂弯矩的非对称分离现象．10.13245/j.hust.241213.F002图2全非线性水弹性结果与试验结果随波高的变化曲线图3为V=18 kn，λ/L=1.0，H=12 m工况下全非线性时域水弹性方法与试验方法预报的船舯垂向弯矩时历曲线，数值结果的中拱弯矩比试验结果高1.9%，中垂弯矩绝对值比试验值低11.5%，垂10.13245/j.hust.241213.F003图3V=18 kn，λ/L=1.0，H=12 m工况船舯垂向弯矩数值与试验时历结果对比向弯矩幅值误差为8.3%，整体幅值的预报较为近．根据时历曲线，数值与试验结果中高频发生的时刻并不一致，说明此时两组弯矩结果的总振动频率存在差异．已知弯矩曲线反映船体的湿模态，根据该工况的试验录像，船体发生了剧烈的砰击及出入水现象，甲板上浪显著，相当于增加了船体的附连水质量，不仅增大了船体的局部及整体载荷，同时也改变了船体的振动频率．而本文数值方法未考虑甲板上浪影响，推测这是图3中二者差异的主要来源之一．图4为船舯垂向弯矩峰谷值时刻的湿表面压力分布图，艏艉具有明显的出入水现象，船舶在砰击作用下产生如图3所示的高频非线性响应．本研究提出的全非线性时域数值方法预报的湿表面压力时空分布特征可以作为分析船舶水动力现象的重要依据．10.13245/j.hust.241213.F004图4船体出入水压力示意图2.3　非线性等效设计波在船体结构强度评估中，各载荷分量之间的组合较为复杂、预报困难，等效设计波法是指把实际的不规则波问题转换成一个等效的规则波，假定组合载荷最大值出现在某一规则波工况下，用幅频响应函数最大值对应的浪向与频率作为等效设计波的参数，设计波波幅ζa由下式计算[17]，ζa=Xmax/Ma，(31)式中：Xmax为1×10-8超越概率水平下的长期值，对应25 a可能出现一次的最大载荷；Ma为单位波幅下的船舯垂向弯矩最大值．本研究采用三维线性波浪载荷预报软件Compass-Walcs-Basic[18]对5 kn极限工况航速及18 kn经济航速不同浪向下船体的幅频响应函数进行预报，浪向角选为0°～180°，相邻浪向间隔为30°．频率为0.01～1.00 rad/s，间隔0.02 rad/s，预报的船舯垂向弯矩结果如图5所示，两航速下垂向弯矩最大值均发生在0°浪向(即迎浪)条件，频率为0.48 rad/s．10.13245/j.hust.241213.F005图5两种航速下的幅频响应函数基于北大西洋海况，采用ISSC双参数海浪谱，通过Compass-Walcs-Basic软件进行波浪载荷的长期分析，确定不同超越概率P下的长期值X(X1为5 kn航速的长期值，X2为18 kn航速的长期值)，如表1所示．10.13245/j.hust.241213.T001表1不同超越概率下的长期预报结果PX1/(GN∙m)X2/(GN∙m)1×10-10.810.791×10-21.571.571×10-32.392.391×10-43.233.251×10-54.104.141×10-64.985.031×10-75.865.931×10-86.726.811×10-97.567.68表2为三维线性波浪载荷预报软件Compass-Walcs-Basic计算的单位波幅下船舯垂向弯矩最大值Ma，1×10-8超越概率水平下的长期值Xmax，以及等效设计波浪向角β、频率ω及幅值ζa等参数．10.13245/j.hust.241213.T002表2长期预报及设计波参数V/knMa/(GN∙m)β/(°)ω/(rad∙s-1)Xmax/(GN∙m)ζa/m50.7100.486.729.5180.7500.486.819.1根据表2，航速18 kn与航速5 kn的线性计算结果最大值相差不大，因此计算得到的等效设计波波高也较接近．实际上对于超大型集装箱船，随着航速的增加，船舶在大幅波浪下航行的试验波浪载荷显示出较强的非线性，导致合成弯矩幅值远超波浪弯矩幅值，采用表2中的等效设计波波幅来预报非线性波浪载荷并不合理．因此，本研究用非线性设计波代替线性设计波[19]，避免因设计波波高过大造成预报的设计载荷值过大．图6为全非线性时域水弹性方法预报的不同波高下的船舯垂向弯矩幅值，用该航速下1×10-8超越概率水平对应的长期计算值(红色虚线)对该曲线进行截取，二者的交点对应的横坐标则代表非线性设计波波高．10.13245/j.hust.241213.F006图61×10-8超越概率水平得到的非线性设计波波高对于5 kn航速来说，非线性设计波波高为15.8 m，为线性设计波波高的83%．在18 kn航速下该值为12.0 m，仅为线性设计波的66%．在上述基础上，通过全非线性时域水弹性数值方法预报非线性等效设计波下的船体波浪载荷，并与规范结果及模型试验结果进行综合对比分析，最终可确定船体的设计载荷，本研究提出的船体设计载荷评估流程如图7所示．10.13245/j.hust.241213.F007图7大型船舶的设计载荷评估流程3 规范-数值-试验设计载荷综合分析表3列出了5 kn和18 kn航速下1×10-8超越概率水平对应的长期值，全非线性数值方法预报的非线性等效设计波结果，18 kn航速下非线性等效设计波的试验结果，UR S11A准则的结果对比情况，MH，MS，MA分别为中拱弯矩、中垂弯矩及幅值．根据长期分析结果，UR S11A准则的载荷值超过了1×10-8超越概率水平的长期值，5 kn航速下其载荷水平对应1×10-8.6超越概率水平，18 kn航速下对应1×10-8.5超越概率水平，而数值及试验的非线性设计波结果基于1×10-8超越概率水平的长期值．10.13245/j.hust.241213.T003表3不同方法得到的设计载荷预报方法MHMSMAMH/MSUR S11A5.299.227.260.36/0.645 kn长期值6.726.726.720.5/0.518 kn长期值6.816.816.810.5/0.55 kn非线性设计波数值结果5.318.406.850.39/0.6118 kn非线性设计波数值结果4.908.746.820.36/0.6418 kn非线性设计波试验结果5.159.497.320.35/0.65GN∙m本研究全非线性时域水弹性数值方法预报的18 kn经济航速下的非线性设计波的弯矩幅值略低于5 kn航速下的非线性设计波值，但中垂弯矩绝对值更大，中拱/中垂值与UR S11A准则一致，均为0.36/0.64，与该工况下分段模型试验的比值0.35/0.65较一致．尽管UR S11A准则对船体的非线性因素考虑较为全面，计入了艏部外飘的影响，预报的弯矩幅值较高，但比1×10-8超越概率水平下非线性等效设计波的试验值低0.8%，最终设计载荷采用非线性等效设计波的试验结果．4 结论本研究建立了有航速船舶的全非线性时域水弹性数值模型，能够较为准确地预报船体的波浪载荷．基于超大型集装箱船的经济效益与航行工况，提出采用全非线性水弹性时域数值方法预报服务工况(经济航速18 kn)及极限工况(航速5 kn)非线性设计波的方法，建立了一套适用于大型船舶的设计载荷评估流程，结合规范-数值-试验手段综合确定船舶的设计载荷．可得到以下结论：a．对于本研究的超大型集装箱船，服务工况下船体的非线性载荷更突出，其中垂弯矩极值高于极限工况结果，说明对于大型船舶高航速较大波高下的非线性载荷可能高于低速极端工况，更加危险；b．UR S11A准则与全非线性时域水弹性方法预报的设计载荷均能较好地反映船体梁垂向弯矩中拱与中垂的强非对称性；c．IACS的UR S11A准则对本研究超大型集装箱船极限工况结果具有一定安全裕度，服务工况下考虑不足，建议该准则设置安全系数．