0 引言
2016年全球风电新增的装机容量超过54.6 GW，累计容量达到486.7 GW [1]。随着海上风电的蓬勃发展，海上风力机逐渐增多，且多接近贸易繁忙的海运路线，故服役期内的海上风力机基础易受船舶撞击[2-4]。因此，探索有效的防护装置以减小船舶碰撞的危害，为风力机的安全运行提供保障，具有重要工程应用价值。
诸多学者对海上风力机支撑结构的船舶碰撞问题进行了研究。RAMBERG[5]模拟了不同质量船舶在不同角度下对海上风力机导管架基础连接点的碰撞，分析了不同状态下的导管架屈曲特性。BIEHL等[6]模拟海上风力机单立柱三桩基础的船舶碰撞过程，研究不同船舶速度下的结构抗撞特性及基础损伤情况，分析碰撞过程中的能量变化与结构损伤。DING等[7]通过模拟船舶以不同速度撞击重力式基础海上风力机，分析塔架结构的剪力弯矩特性以及重力式基础的应力应变。AMDAHL等[8]分析了单桩基础风力机在不同速度船舶撞击及风载荷作用下的塔架变形、塔顶位移响应。上述研究仅考虑海上风力机受撞损伤及动力响应，未对海上风力机结构的安全防护提出行之有效的方案。在海洋工程结构受碰撞过程中，防护装置主要运用在船-船碰撞、船-桥碰撞及船-石油平台碰撞的场合，用于海上风力机-船舶碰撞的较少。LEHMANN等[9]设计的钢结构防护装置能吸收碰撞过程中的巨大能量，通过试验与数值模拟验证了防护装置的有效性。文献[10]在研究船舶碰撞桥墩的过程中，提出了新型组合钢结构的防护装置；文献[11]指出，在船舶与海洋平台碰撞的过程中，橡胶护舷可以减小碰撞力和结构的损伤。
本文借鉴传统海洋工程结构安全防护装置的研究设计经验，针对单桩柱海上风力机，设计出4种型式的防护装置，模拟船舶在不同速度下与单桩柱基础的碰撞过程，分析其能量变化、塔架应力和撞深、接触力及塔顶风力机的响应，来验证防护装置的性能。
1 理论基础
1.1 控制方程
船舶与海上风力机碰撞运动控制方程为[12]
(1)
(2)
式中，M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；F(t)为作用在结构上的撞击力；H(t)为沙漏力；
分别为t时刻的加速度、速度和位移；m为船体质量；
为船体运动加速度；f(t)为t时刻的船舶撞击力。
1.2 材料本构模型
1.2.1 钢的本构方程
风力机塔架和船首材料为Q235B，单桩柱基础和防护装置外壳材料为Q345，参数如表1所示。由LS-DYNA提供的非线性弹塑性材料模型[13]是基于Cowper-Symonds关系式建立的[14], 可很好地模拟单桩柱基础撞击作用下的材料特性，如下式所示:
(3)
式中，σy为极限屈服应力；σ0为初始屈服应力；C、p为材料应变率常数，C=40.5，p=5；β为可调参数，塑性随动强化模型取0，等向强化模型取1；EP为塑性强化模量；
为塑性应变率。
1.2.2 橡胶的本构方程
碰撞过程中，各向同性橡胶材料的变形都是超弹性材料的均匀变形[15]。橡胶材料应变能函数有两种表达方式：由变形张量的3个不变量I1、I2、I3表示的应变能函数W(I1，I2，I3)；由主伸长比λ1、λ2和λ3表示的应变能函数H(λ1，λ2，λ3)。变形张量的不变量与主伸长比的关系式如下：
(4)
表1 部件的材料属性参数
Tab.1 Material properties of the parts
材料Q235BQ345密度ρ(kg/m3)78507870弹性模量E(GPa)210212泊松比υ0.270.31初始屈服应力σ0(MPa)235345切线模量Et(MPa)11801291可调参数β00应变率常数C(s-1)40.540.5应变率常数P55失效应变ε0.340.34
(5)
(6)
λi=1+εi
(7)
式中，εi为主轴方向的应变；i为张量序号，i=1，2，3。
橡胶属超弹性材料，其应变能函数由上述形式转化为多项式形式后，由应变偏量能和体积应变构成：
(8)
式中，W为应变势能；N为多项式阶数；J为弹性体积比；Di表示材料的压缩性；
为应变偏量；
为体积应变。
对于多项式中的应变偏量，橡胶的初始剪切模量G0、初始体积模量K0都取决于多形式的一阶系数，即G0=2(C10+C01)，K0=2/D1，其中，C10、C01、D1为材料常数。
对于式(8)，如果N=1，则仅保留线性部分的应变能，即Moony-Rivlin模型：
W=C10(I1-3)+C01(I2-3)+(J-1)2/D1
(9)
1.2.3 泡沫铝的本构方程
泡沫铝材料的屈服应力函数为
(10)
式中，Y为屈服应力；
为等效应力；σe为von Mises有效应力；σm为平均应力；α为屈服面形状参数；σp为平台应力；
为等效应变；εD为应变密度；γ为应力函数的材料系数，α2为硬化模型参数的工程应力，MPa；β为硬化系数。
泡沫铝材料隔垫的具体材料参数如表2所示。
1.3 接触定义
为更好地模拟船舶与海上风力机基础的碰撞过程，避免初始接触，将两者初始间隔距离设为0.6 m。在LS-DYNA中，将船首与海上风力机防护装置间设置为面面接触，塔架与桩柱之间设置为
表2 泡沫铝材料隔垫的材料参数
Tab.2 The material parameters for aluminum foam pad
材料参数数值材料参数数值密度ρ(kg/m3)178γ(MPa)1.87弹性模量E(GPa)370α2(MPa)93.5应变密度εD2.77σp(MPa)1.15可调参数β5.79
自动单面接触，避免船体穿透基础，橡胶、泡沫铝和外壳采用固连接触。
碰撞过程中，船首与防护装置接触时的摩擦力为
Fy=μ|fs|
(11)
μ=μd+(μs-μd)e-d|v|
(12)
式中，fs为节点法向接触力；μ为摩擦因数；μd为动摩擦因数；μs为静摩擦因数；d为衰减系数；v为接触面间相对速度。
船首和防护装置外壳之间的动摩擦因数、静摩擦因数取0.2，衰减系数取0，橡胶与塔架的静摩擦因数为0.9，泡沫铝与塔架的静摩擦因数为0.4。
2 研究对象与模型
2.1 单桩柱风力机
国内某海上风电场的风力机及塔架主要参数如表3所示。
表3 3 MW风力机部件主要参数
Tab.3 Main parameters of 3 MW offshore
wind turbine components
部件几何尺寸(m)质量(t)叶片2.4×44×4.211.2机舱13×5×6.24115轮毂5×528.5塔架上(3.09～3.9)×32.158中(3.9～4.5)×30.81596下4.5×14.8680单桩基础连接段(4.5～4.8)×7.7741管桩4.8×53283
单桩柱海上风力机上部结构由塔架、轮毂、机舱和叶片组成，其中桩基入土深度为41 m。单桩柱式海上风力机模型如图1所示。
图1 单桩柱式海上风力机模型
Fig.1 Monopile foundation offshore wind turbine model
海上风力机底部往往采取水泥浇筑等加固措施。因此，在ANSYS/LS-DYNA中将桩柱基础在泥土中的部分假定为刚性约束，忽略桩柱水平位移。塔顶部分较为复杂且网格划分数量巨大，占用较多计算资源，为提高效率，采用集中质量方法，用固定质量点代替风力机顶部(叶片、轮毂、机舱)的质量。
2.2 防护装置
为保证防护装置受船舶撞击后，不脱离单桩柱基础，模拟中的防护装置质量不能超过海上风力机总质量的15%[16]。防护装置安装在塔架下部，内部是内径4.5 m的橡胶，外部是一层10 mm厚的钢制防护壳，整体结构高6 m，质量70.7 t，占海上风力机总质量的9.6%，符合设计标准。泡沫铝密度更小，故设计表4所示的4种型式防护装置。
表4 防护装置材料和厚度
Tab.4 Material and thickness of protective device
类型ABCD内层材料橡胶泡沫铝橡胶泡沫铝厚度(m)0.30.30.150.15中间材料泡沫铝橡胶厚度(m)0.150.15外层材料钢板厚度(mm)10
对防护装置进行网格划分，橡胶和防护壳采用扫掠的网格划分方法，防护壳单元尺寸为0.1 m，橡胶单元尺寸为0.2 m，有限元模型如图2所示。
图2 防护装置有限元模型
Fig.2 Finite element model of protective device
2.3 船舶模型
固定式基础的海上风力机所处海域为近海，来往船只种类较多，本文选用较为常见的散货船，船首型式为飞剪型，主要尺寸参数如表5所示，实体模型如图3所示。模拟过程中，忽略碰撞船舶的变形，假定碰撞船舶为刚性体。
表5 船舶结构主要参数
Tab.5 Main parameters of ship structure m
总长型宽型深吃水102.410.625.84.12
图3 碰撞船舶模型
Fig.3 Collision ship model
碰撞船舶主要分为船首、船身和船尾，通过LS-DYNA中的关键字*ELEMENT_MASS控制船舶整体的重心和质量。为更好地分析防护装置对塔架主体结构的保护作用，本文模拟2 000 t、5 000 t的船舶以1 m/s、2 m/s、3 m/s的速度与单桩柱海上风力机的碰撞，船舶行驶速度沿x轴正方向。
碰撞过程中，船舶与海水的相互作用不可忽略，故通过流固耦合模型或附加质量模型表示船舶与海水的相互作用过程[17-18]。船舶与海水相互作用的耦合计算复杂且耗时[19-20]，因此本文中的船舶与水的相互作用选用附加质量法。
在船舶有限元模型网格划分中，单元厚度为20 mm，网格划分单元尺寸为0.2 m[21]。
3 结果与分析
3.1 无防护装置
3.1.1 全局能量分析
船舶撞击单桩柱海上风力机的过程中，船舶的初始动能转化为船舶和塔架的内能(变形能)、船舶的剩余动能，以及接触摩擦产生的滑移能、显式分析中采用缩减积分产生的沙漏能、塔架在阻尼作用下产生的阻尼能。数值模拟计算中，一般以沙漏能应小于总能量的5%作为判断数据准确性的依据 [14]。图4分别为5 000 t船舶以不同速度v撞击塔架下端全局的能量曲线图。
图4 5 000 t船舶撞击塔架的能量曲线图
Fig.4 Energy curve of 5 000 t ship impacts the tower
由图4可知，沙漏能低于总能量的1.5%，表明有限元模型合理，计算结果准确。由图4a可知，碰撞过程中，结构变形能最大值达到3.93 MJ，碰撞结束后，结构变形能在2.72～2.80 MJ之间波动；碰撞系统动能减小至接近于0后，反弹至1.29 MJ，然后在海水阻尼作用下缓慢减小。由图4b可知，碰撞过程中的结构变形能最大值达到3.63 MJ。由于防护装置中的橡胶的超弹性作用，动能先减小到0.52 MJ，然后迅速回升，此时橡胶受冲击压缩后反弹撞击船舶，减弱了对塔架主体的伤害。
3.1.2 塔架结构分析
图5是船舶冲击单桩柱基础碰撞区域的应力云图(船舶碰撞区域的正视图)及塔架撞深截取的云图。由图5可知，船舶速度增加时，结构损伤程度也随之增加。撞击速度v为1 m/s时，塔架最大应力在碰撞区域，且未发生变形；撞击速度为2 m/s、3 m/s时，船舶碰撞区应力迅速下降，碰撞区域中心出现应力最小值，开始产生明显变形；撞击速度为4 m/s时，塔架的变形δ最大，最小应力集中区域范围扩大，说明此时的塑性区域吸能特性使应力集中区域的能量释放。
图5 塔架下端结构应力及撞深
Fig.5 Stress and depth of tower structure
图6所示为不同速度下 5 000 t船舶碰撞单桩柱海上风力机的接触力曲线，为更直观分析接触力曲线，以速度4 m/s行驶的船舶为基准，并将接触初始点调至相同起点。
图6 不同速度船舶撞击接触力曲线
Fig.6 Contact force curve of ship impact at different speeds
由图6可知，碰撞初始阶段的船舶接触力曲线高度重合，且处于线性阶段，此时的碰撞属于弹性碰撞；随着船速的增加，接触力曲线快速分离，曲线的非线性波动特征显著，接触力曲线每次由波峰至波谷的过程存在力的卸载，表明船首与塔架变形损伤过程在不断变化；4 m/s速度时，接触力在1.5～1.9 s的卸载最显著，1.9～2.6 s的接触力并没有降低，结果说明，此时的塔架发生了较大的变形，初始动能主要转化为船舶的结构变形能，符合图5所示，随着船舶速度增加，塔架受撞击后的深度增加。
船舶以4 m/s速度碰撞时，0.50～0.54 s内，船舶开始与风力机塔架碰撞，接触力线性增长，此阶段的碰撞为弹性碰撞，塔架主要产生弹性变形；0.54～1.44 s时，接触力总体呈上升趋势，最大值为13.9 MN，此阶段碰撞为弹塑性碰撞；1.44 s后，曲线具有明显的非线性波动特征，直到碰撞结束(3.80 s)，此阶段的碰撞为塑性碰撞。对比分析每个阶段碰撞持续的时间可知，塑性碰撞的持续时间2.44 s明显长于弹性碰撞的持续时间0.04 s和弹塑性碰撞的持续时间0.9 s。
3.2 橡胶防护装置
3.2.1 塔架结构变形能分析
图7 塔架下端结构变形能
Fig.7 Structural deformation energy of the tower
图7所示为2 000 t船舶以不同速度撞击海上风力机时，塔架受碰区的结构变形能。由图7可知，碰撞过程中，有防护装置的最大结构变形能都小于无防护装置的，表明防护装置吸收了碰撞能量，减少了塔架的变形；随着船舶速度的增加，有防护装置的塔架在船舶速度分别为1 m/s、2 m/s、3 m/s的结构变形能相对于无防护装置分别减少了35.9%、26.6%和14.6%，说明该防护装置对于低速船舶的效果更好。
3.2.2 接触力分析
主从算法在求解的过程中，在垂直于主面的方向上施加一个作用力(接触力)以阻止从属节点的穿透。由图8可知，无防护装置时，在碰撞开始阶段，接触力曲线线性增长且斜率较大，此时的碰撞属于弹性碰撞；此后，曲线的非线性波动特征逐渐显著，这表示船首构件和塔架的变形或受损，此碰撞为弹塑性碰撞；增加防护装置后，塔架也存在由弹性变形到弹塑性、塑性变形的阶段，但碰撞持续的时间都在延长，在船舶速度为1 m/s、2 m/s、3 m/s时，最大接触力分别减少了45.6%、45.0%、36.1%。
图8 不同速度船舶撞击海上风力机接触力曲线
Fig.8 Contact force curve of ship impact on offshore wind turbine under different velocity
3.2.3 塔架应力云图
由图9可知，船舶以2 m/s速度撞击无防护装置塔架时，应力最大值达到278 MPa，超过屈服极限235 MPa，而船舶以3 m/s撞击有防护装置的塔架时，应力仅有193 MPa，由此可知，防护装置能很好地保护海上风力机基础结构并吸收船舶初始动能，减少碰撞载荷，避免应力集中。
图9 塔架应力云图
Fig.9 Von Mises stress of tower
3.3 三种防护装置的对比
为分析橡胶、泡沫铝及两种材料组合的防护装置性能，模拟5 000 t船舶以1 m/s速度正面撞击塔架下端。图10所示为3种型式防护装置下，海上风力机的动力响应。3种防护装置中，由塔架下端的结构变形能、接触力、塔顶位移及塔顶加速度响应参数对比分析可知，防护装置C即橡胶和泡沫铝组合，抑制塔顶动力响应的效果更好。
图10 3种防护装置的响应对比
Fig.10 Response characteristic curve of 3 protective devices
3.4 防护装置材料铺层顺序响应的对比
橡胶有很好的超弹性，泡沫铝有高阻尼减震性能及良好的冲击能量吸收率。为探究防护装置铺层顺序对撞击过程中风力机响应的抑制效果，设计型式C、D两种防护装置，模拟工况为5 000 t船舶以3 m/s速度正面撞击塔架下端。图11所示为C、D防护装置下，海上风力机的响应特性。D型防护装置的结构变形能、接触力、加速度响应稍小于C型防护装置，对塔顶位移响应抑制的效果明显。
图11 2种防护装置的响应对比
Fig.11 Response characteristic curve of 2 protective devices
4 结论
(1)无防护装置时，撞击过程中的塔架应力很快超过材料的屈服极限，撞深随船舶速度增加而增加。
(2)船舶速度在1 m/s、2 m/s、3 m/s时，A型防护装置的结构变形能分别是无防护装置的22.13%、23.80%和42.58%，最大接触力分别是无防护装置的54.38%、54.95%、63.92%，碰撞持续的时间都在增加。
(3)由响应参数的对比分析可知， C型(橡胶和泡沫铝组合)防护装置对响应的抑制效果较A型、B型更好。
(4)D型防护装置对塔顶位移响应的抑制效果明显，但在结构变形能、接触力、加速度响应与C型较为相近。
[1] Global Wind Energy Council. Global Wind Report[R]. Washington D C: Global Wind Energy Council, 2016.
[2] 丁勤卫，李春，叶柯华，等．风波流对多平台阵列浮式风机Spar平台运动特性的影响[J]. 农业工程学报，2016， 32(21)：223-229.
DING Qinwei, LI Chun, YE Kehua, et al. Effect of Wind, Wave and Current on Movement Characteristics of Array of Floating Wind Turbine Spar Platform[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(21): 223-229.
[3] 丁勤卫，李春，周国龙，等．海风力机动态响应对比[J]．动力工程学报，2016，36(1)：65-74.
DING Qinwei, LI Chun, ZHOU Guolong, et al. Comparison of Dynamic Response between Stationary and Floating Wind Turbine[J]. Journal of Chinese Society of Power Engineering, 2016, 36(1): 65-74.
[4] DING Qinwei, LI Chun. Research on the Influence of Helical Strakes on Dynamic Response of Floating Wind Turbine Platform[J]. China Ocean Engineering, 2017, 31(2): 131-140.
[5] RAMBERG H F. High Energy Ship Collisions with Bottom Supported Offshore Wind Turbines[D]. Trondheim: Norwegian University of Science and Technology, 2011.
[6] BIEHL F, LEHMANN E. Collisions of Ships with Offshore Wind Turbines: Calculation and Risk Evaluation[C]//25th International Conference on Offshore Mechanics and Arctic Engineering. Berlin: Springer, 2006.
[7] DING Hongyan, ZHU Qi, ZHANG Puyang. Dynamic Simulation on Collision between Ship and Offshore Wind Turbine[J]. Transactions of Tianjin University, 2014, 20(1): 1-6.
[8] AMDAHL B, HERVE L S, LOIC B, et al. Ship Collision Analysis on Offshore Wind Turbine Monopole Foundations[J]. Marine Structures, 2017, 51: 220-241.
[9] LEHMANN E, PESCHMANN J. Energy Absorption by the Steel Structure of Ships in the Event of Collisions[J]. Marine Structures, 2002, 15(4): 429-441.
[10] WANG Lili, YANG Liming, HUANG Dejin, et al. An Impact Dynamics Analysis on a New Crashworthy Device against Ship-bridge Collision[J]. International Journal of Impact Engineering, 2008,35(8):895-904.
[11] 鲍莹斌，李润培，顾永宁．船舶与海洋平台碰撞的动塑性分析[J]．上海交通大学学报，1997，31(7)：40-44.
BAO Yingbin, LI Runpei, GU Yongning. Dynamic Plastic Analysis of Ship Platform Collision[J]. Journal of Shanghai Jiaotong University, 1997, 31(7): 40-44.
[12] 白金泽. LS-DYNA 3D理论基础与实例分析[M]．北京：科学出版社，2005.
BAI Jinze. Basic Theory and Example Analysis of LS-DYNA 3D[M]. Beijing: Science Press, 2005.
[13] GLADMAN B. LS-DYNA Keyword User’s Manual[M]. Livermore: Livermore Software Technology Corporation, 2007.
[14] 石少卿，康建功，王敏，等．ANSYS/LS-DYNA在爆炸与冲击领域内的工程应用[M]．北京：中国建筑工业出版社，2011.
SHI Shaoqing, KANG Jiangong, WANG Min, et al. ANSYS/LS-DYNA in Explosion and Impact in the Field of Engineering Application[M]. Beijing: China Architecture & Building Press,2011.
[15] 李艳敏，马玉宏，罗佳润，等．考虑老化时间影响的隔震支座橡胶本构Moony-Rivlin模型常数研究[J]. 振动与冲击，2016，35(16)：164-169.
LI Yanmin, MA Yuhong, LUO Jiarun, et al. The Effect of Aging on the Material Constant of the Rubber Isolators Constitutive Model Moony-Rivlin[J]. Journal of Vibration and Shock, 2016, 35(16): 164-169.
[16] LIU Chunguang, HAO Ertong, ZHANG Shibo. Optimization and Application of a Crashworthy Device for the Monopile Offshore Wind Turbine against Ship Impact[J]. Applied Ocean Research, 2015, 51: 129-137.
[17] 王自力，蒋志勇，顾永宁．船舶碰撞数值仿真的附加质量模型[J]．爆炸与冲击，2002，22(4)：321-326.
WANG Zili, JIANG Zhiyong, GU Yongning. An Added Water Mass Model for Numerical Simulation of Ship Collision[J]. Explosion and Shock Waves, 2002, 22(4): 321-326.
[18] 姜华，王君杰，贺拴海．钢筋混凝土梁桥船舶撞击连续倒塌数值模拟[J]．振动与冲击，2012，31(10): 68-73.
JIANG Hua, WANG Junjie, HE Shuanhai. Numerical Simulation on Continuous Collapse of Reinforced Concrete Girder Bridge Subjected to Vessel Collison[J]. Journal of Vibration and Shock, 2012, 31(10): 68-73.
[19] PELERSEN M J. Dynamics of Ship Collisions[J]. Ocean Engineering, 2012, 9(4): 295-329.
[20] 李艳贞，胡志强，邹早建．海上风电站遭遇船舶侧向撞击时的结构动力响应分析[J]．振动与冲击，2010，29(10)：123-126.
LI Yanzhen, HU Zhiqiang, ZOU Zaojian. Structural Dynamic Response Analysis of an Offshore Wind Turbine Stricken by a Ship in Lateral Direction[J]. Journal of Vibration and Shock, 2010, 29(10): 123-126.
[21] Le SOURNE H, BARRERA A, MALIAKEL J B. Numerical Crashworthiness Analysis of an Offshore Wind Turbine Jacket Impacted by a Ship[J]. Journal of Marine Science and Technology, 2015，23 (5): 694-704.
