为增强寒区海洋平台结构的抗冰性能，锥体海洋平台得到了广泛的应用。锥体结构可使海冰由挤压破坏转变成弯曲破坏，从而有效地降低海洋平台结构冰载荷幅值[1－2]。然而现场监测表明锥体结构在海冰作用下仍然存在显著的动力问题[3－4]。由于锥体长期与海冰相互作用，需要对其强度、疲劳寿命以及冰激振动等因素进行研究。
目前，对锥体海洋平台结构与海冰相互作用过程，主要通过现场实测、室内试验以及数值模拟展开研究。在渤海，通过对锥体海洋平台的现场监测，对锥体平台的冰载荷、结构振动以及冰况进行统计分析[5]。在室内冰池内对冰-锥体相互作用过程也进行了一系列的研究[6]。然而，由于冰-锥体结构作用过程的复杂性，冰厚、冰速、海冰强度以及结构接触宽度等参数无法准确地得到。因此，很难通过现场实测和室内试验系统分析不同参数对平台结构冰载荷的影响。此外，数值模拟方法可以准确地改变计算参数[7－8]，对理解冰-海洋平台相互作用机制具有成本低、效率高的优势。
针对冰-锥体海洋结构相互作用过程的数值模拟，主要通过基于连续介质方法以及离散介质方法展开分析。武文华等[9]通过显式动力分析软件LS-DYNA模拟了海冰和锥体的相互作用。Sand等[10]将海冰假设为一种材料非线性模型，采用有限元方法计算海冰与平台结构相互作用的冰载荷。但由于海冰大多以离散碎块的形式漂浮于海面，具有很强的非连续性，且当海冰相互碰撞或与海洋平台结构相互作用时，会进一步发生破碎。基于非连续介质理论的离散元方法(DEM)在模拟连续介质材料的断裂方面体现了较好的力学性能[11]。目前，基于DEM的海冰模型发展迅速，可有效模拟海冰的破碎现象并能较合理地反映不同结构的冰载荷特性[12－13]。其中，通过构造粘结-破碎球体离散单元可模拟冰-锥相互作用以研究锥体结构的冰载荷[14]。为分析海洋平台结构在海冰作用下的振动特性，海洋平台结构可采用有限元(FEM)计算，海冰通过 DEM 构建，由此建立DEM-FEM耦合的方法。
目前，针对不同的工程问题存在不同的DEM-FEM 耦合方法。其中，Munjiza等[15]提出了模拟脆性材料断裂的DEM-FEM耦合方法。然而，该方法并不能模拟离散介质与连续介质接触的问题。因此，针对该问题的DEM-FEM耦合方法得到了发展[16]，并已应用到多个工程领域中，如轮胎与路面的相互作用[17－18]、列车载荷对有砟铁路地基的动力影响[19]、轨枕-道床-路基的耦合模型[20]、海底管线冲击问题[21]以及颗粒材料冲击问题[22]等。为满足实际工程需求，对DEM-FEM耦合方法的计算规模和效率提出了更高的要求。通过基于 GPU并行的 DEM-FEM 耦合方法可提高海冰-导管架平台相互作用的计算规模和效率。
本文在海冰DEM模型的基础上，建立了针对锥体海洋平台结构冰激振动的高性能DEM-FEM耦合模型，发展了GPU并行环境下DEM和FEM间参数传递算法以及球体离散元单元与壳单元在界面上的耦合算法。通过结构冰载荷、振动响应及冰压力和应力分布对耦合模型进行验证分析。
1 海冰离散元模型
本文通过平整冰离散元单元模型计算锥体结构的冰载荷。图1为平整冰与锥体结构相互作用的示意图，其中平整冰采用具有粘结-破坏的球体离散元单元构造。为模拟无限域的海冰，本文假设海冰模型在两侧的边界限制纵向(y方向)和垂向(z方向)的位移，并且在海冰的后侧施加流速边界，如图1所示。此外，在流体动力学方面，仅考虑水对离散单元的拖曳力以及浮力。
图1 平整冰-锥体结构相互作用示意图
Fig.1 Schematic of the interaction between the level ice and conical structure
1.1 海冰离散单元的运动方程
海冰离散元单元的运动可分为平动和转动，依据牛顿第二定律表示为：
式中：mi和Ii分别为i单元的质量和惯性矩；ui和θi 分别为 i单元的位移和转角； Fi和Mi分别为i单元受到的总外力和中心弯矩，可表示为：
式中：Fc为海冰颗粒单元的接触力；Fg为单元受到的重力；Fb和Fd分别为海冰单元受到的浮力和拖曳力，可表示为：
式中：ρw为水的密度；g为重力加速度；Vsub为单元浸入水中的体积；Cd为拖曳系数；vi为i单元的速度；vw为水流速度。
1.2 海冰离散单元的接触模型
基于粘弹性接触模型以及库仑摩擦准则计算颗粒与颗粒以及结构间的接触力[14]。图2为两球体离散单元接触力模型，主要分为法向作用力Fn和切向作用力Fs。其中，法向接触力可分为弹性力Fne和粘滞力Fnv，表示为：
式中：
为两颗粒单元的重叠量；
为两颗粒的法向相对速度；nij为颗粒的接触法向向量；kn和cn分别为接触的法向刚度和阻尼系数，可表示为：
式中：E为材料的弹性模量；Ri和Rj为单元的半径；Mp为单元质量；
为无量纲阻尼系数；ep为海冰颗粒回弹系数。
图2 球体离散单元接触力模型
Fig.2 Contact model in the spherical DEM
颗粒单元间的切向接触力可写作：
式中：Fs为切向力；
分别为颗粒切向的相对位移和速度；μp为颗粒间的摩擦系数；tij为切向向量；ks=kn/γ和cs=cn/γ分别为切向刚度系数和切向阻尼系数，这里γ=2[23]。
1.3 海冰离散单元粘结-破碎模型
每个粘结颗粒单元间具有一定的粘结强度以传递两个粘结单元间的作用力和力矩，如图3所示。粘结单元间的最大拉应力和剪应力可依据梁的拉伸、扭转和弯曲理论按下式计算：
式中：Fn和Fs分别为单元粘结面上的法向和切向作用力；Mn和Ms分别为粘结面的扭矩和力矩；A、
、J和I分别为颗粒接触的粘结面积、有效粘结半径、极惯性矩和惯性矩，即：
式中，Ri、Rj分别为颗粒i和j的半径。
图3 球体离散单元平行粘结模型
Fig.3 Parallel bonding model of the bonded spherical elements
通过多年现场实测数据发现渤海冰速主要集中在 0.2 m/s～0.8 m/s[24]，海冰主要发生脆性破坏[25]。因此本文海冰离散单元模型采用脆性破坏准则，即当单元间的最大拉伸σmax或剪切应力τmax超过单元间的法向粘结强度
或切向粘结强度
单元间的粘结会破坏失效。本文采用累积损伤的线性软化模型消除粘结颗粒破碎产生的能量[26]。本文离散元模型中，考虑锥体结构前海冰主要发生弯曲破坏[2]，取法向粘结强度和切向粘结强度相等
且与海冰的弯曲强度σf相关[23]，即：
海冰的弯曲强度与其卤水体积密切相关，可写作[27]：
式中：vb为海冰的卤水体积，可表示为温度和盐度的函数[28]。
式中：Tice/(℃)为海冰的温度；Sice/‰为海冰的盐度。
2 锥体导管架平台结构的有限元模型
本文采用渤海 JZ20-2 NW 单桩锥体海洋平台分析其在海冰作用下结构的动力响应。该平台主要由上部模块、导管架、桩基和破冰锥组成，如图4(a)所示。在有限元计算中，在保证主体结构几何形状以及振动频率和振型的真实性的前提下，对上部模块和桩基进行一定的简化，即结构的主体采用梁单元进行计算。破冰锥体是海洋平台抗冰的关键部分，安装于水线附近上的桩腿上。为分析海冰作用下锥体的强度，这里采用平板型壳单元进行有限元分析并用梁单元模拟其内部的加劲肋。综上所述，本文采用梁-壳组合单元构造海洋平台有限元模型,模型共有1026个节点和1530个单元，其中梁单元274个、壳单元1214个，计算模型如图4(b)所示。
图4 渤海锥体导管架海洋平台及其有限元模型
Fig.4 Conical jacket platform in the Bohai Sea (China)and FEM model
平板型壳单元可以看成是平面应力单元和平板弯曲单元的组合，如图5所示。对于各向同性的平板型壳体单元，这两部分变形是相互独立的，所以需要精细划分单元使其解收敛于原曲壳的解。为避免总体刚度矩阵由于单元共面或接近共面而出现的奇异性问题，这里采用具有附加刚体转角的三角形膜元 GT9[29]和三角形薄板元 TMT[30]组成的广义协调三角形平板型壳元。其单元刚度矩阵可由这两种单元的刚度阵组合而成，即：
式中：K为平板型壳单元的单元刚度阵；
为平面膜单元刚度系数；
为平板弯曲单元的刚度系数。
图5 三角形平板型壳元
Fig.5 Triangle plane shell element
这里采用Crank-Nicolson法求解结构的动力学方程，即：
式中：M、C和K分别为结构的质量、阻尼以及刚度矩阵；u、
和
分别为结构的位移、速度以及加速度矩阵；Δt为时间间隔；θ为计算参数，为保证计算的稳定性，取θ=0.5。其中，C采用瑞雷阻尼构造，质量矩阵与刚度矩阵的线性组合，即：
式中，参数α和β可分别通过结构任意两个模态阻尼比
和固有频率
确定，可表示为：
由于在冰载荷作用下结构的振动能量主要集中于低阶频率。因此
分别采用JZ20-2 NW平台1阶和10阶阻尼比。ωi和ωj分别为结构的1阶(1.07 Hz)和10阶频率(7.39 Hz)。
3 海冰与海洋平台相互作用的离散元-有限元耦合方法
这里采用DEM-FEM耦合方法计算锥体导管架海洋平台结构的冰激振动。图6为DEM-FEM耦合方法的主要计算过程，DEM部分是基于CUDA/C++在GPU并行环境下编译的[31]，FEM部分则是基于C++串行编译。将离散单元对结构的载荷转换为有限元单元节点载荷传递到有限元模型中，计算结构的动力响应，再进一步将更新后的有限单元的节点响应传递到离散元计算中作为位移边界条件。DEM-FEM在界面上传递参数是该耦合方法的关键问题。
图6 DEM-FEM耦合方法的参数传递
Fig.6 Coupled DEM-FEM scheme with a GPU-based parallel algorithm
3.1 GPU并行环境下的DEM-FEM参数传递
在 GPU并行环境下需要将离散单元与结构计算得到的接触力、作用位置以及接触数量等变量传递到有限元中进行结构的动力计算。然而，GPU是对线程进行并行计算，在离散元计算中每个线程代表一个离散单元，因此离散单元计算的结果在GPU内存中是随机分布的，这给参数的传递造成了困难。本文采用CUDA中的Thrust库对离散元的计算结果进行重新排序，它采用高度优化的基数排序算法，要比基于数据之间比较的排序算法速度快很多。基于 GPU并行环境下参数传递算法主要包括以下三个步骤。
首先，对结构有限单元进行并行搜索计算，每个有限单元代表GPU中的一个线程(threadx.x)。每个有限单元只搜索领域网格的离散单元。因为所有单元同时进行计算，所以这个过程是高效的。如果离散单元与结构单元发生接触，则需要将相关计算参数储存到数组中(ContactN，DEM_force，Sort_Num等)。“ContactN”储存每个有限单元与离散单元的接触数量，“DEM_force”储存每个有限单元与离散单元的接触力，“Sort_Num”储存线程号的负值。另外，当离散单元与结构没有接触时，令“Sort_Num”数组相应位置为1。
然后，采用 Thrust库中的 inclusive_scan和sort_by_key函数对耦合参数进行排序，如图7所示。图7(a)为通过inclusive_scan函数对颗粒与结构的接触数量进行排序，从图中可以看出排序后数组的最后一项为离散元与结构接触的总数。图7(b)是通过 sort_by_key函数对有效的接触力进行提取，需建立与线程号符号相反的数组。通过对该数据进行排序，可提取出连续有效的接触力、接触位置等数据，以便于传递到有限元计算中。
最后，将排序好的离散元计算参数传递到有限单元中。根据Thrust函数可以很容易地将接触数量、接触力和接触位置等参数传递到CPU的数组中。
图7 GPU并行环境下的参数传递方法
Fig.7 Parameter transfer scheme with the GPU
3.2 球体颗粒单元-三角形壳单元耦合的接触算法
通常球体单元作用于三角形单元的表面，如图8所示。有限元计算需要输入节点的载荷，需要将离散元的接触力通过插值等效到三角形单元的节点上。节点力的插值是基于虚功等效原理，可表示为：
式中：fi为三角形i节点的节点力向量；F为颗粒与三角形单元的接触力向量；Ni为三角形单元的形函数。
图8 三角形单元等效节点力计算模型
Fig.8 Sketch of the contact algorithm between DEM and FEM
这里采用面积坐标作为三角形单元的形函数。三角形的面积坐标可以直接通过接触点的坐标计算得到，而不必先求解广义坐标，然后再由广义坐标得到单元插值函数的复杂过程。对于线性三角形单元，面积坐标可表示为：
式中：Ai为离散元与三角面的作用点与另外两个单元节点围成的面积；A为三角形单元的整体面积。
3.3 DEM-FEM耦合方法的验证
为验证DEM-FEM耦合方法的精度，下面采用单个球体以一定的初始速度冲击平板，如图9所示。相应的几何参数，如表1所示。这里分别采用DEM-FEM耦合方法和LS-DYNA软件对其进行对比分析。在DEM-FEM耦合方法中，球体采用离散元单元构造，平板由平板型壳单元组成；LS-DYNA采用刚性球体冲击由壳单元构造的平板。
表1 颗粒冲击平板的主要计算参数
Table 1 Parameters of the particle impact on the plate
图9 刚性球与平板的相互作用
Fig.9 Schematic of interaction between the rigid ball and plate
这里分别对比平板受到的冲击力和板中心位移，如图10所示。从对比的结果可以发现本文提出的DEM-FEM耦合方法与通过LS-DYNA软件计算得到的结果相一致。因此，该DEM-FEM耦合方法可合理地模拟离散球体单元与连续体单元在界面的相互作用。
图10 DEM-FEM耦合方法与LS-DYNA的计算对比
Fig.10 Comparison of the coupled DEM-FEM method and Ls-dyna
4 锥体海洋平台结构冰激振动分析
渤海 JZ20-2 NW 锥体导管架平台安装了完备的冰激振动立体监测系统，如图11所示。该监测系统主要包括了海冰图像监测和振动加速度传感器等，实现了对该海域冰速、冰厚、冰向、海冰密集度等海冰参数以及结构振动加速度长期系统的现场监测，为我国研究导管架冰激振动提供大量的现场实测数据。
图11 渤海JZ20-2 NW平台冰激振动监测系统
Fig.11 Ice-induced vibration monitoring system of the platform
下面通过渤海JZ20-2NW单锥海洋平台结构的冰载荷和冰激振动对 DEM-FEM 耦合模型进行验证，并对锥体结构的应力分布进行分析。这里采用25种冰况，即不同冰速和冰厚，模拟海冰与单锥平台结构的相互作用，具体参数列于表2，海冰离散元计算参数在表3中列出。图12为冰速0.2 m/s和冰厚 0.2 m冰况下海冰与锥体平台相互作用的过程。从中可以看到平整冰在锥体作用下的破碎情况。海冰与锥体平台相互作用过程中，冰排呈现初次断裂、爬升、二次断裂和清除的过程，并由此引起交变动冰力。
表2 模拟中的冰速和冰厚
Table 2 Values of ice velocity and thickness in the simulations
4.1 锥体结构冰载荷
通过DEM-FEM耦合模型，得到锥体平台的冰载荷时程(v=0.2 m/s,h =0.2 m)，如图13所示。为减小冰力对比的随机性，这里将各冰况下冰载荷时程的幅值进行提取，并对其平均值进行对比分析。由冰荷载时程可以发现，锥体结构的冰荷载表现出很强的随机性并具有一定的周期脉冲特性。这与海冰的现场监测、室内模型试验以及数值模拟的结果相一致[7,32]。
表3 海冰离散元模拟中的计算参数
Table 3 Parameters in the DEM-FEM simulations
图12 海冰与JZ20-2 NW平台相互作用的数值过程
Fig.12 Numerical process of the sea ice-JZ20-2 NW platform interaction
图13 锥体冰力时程
Fig.13 Ice forces on the jacket platform
同时，对不同冰况下冰载荷的平均值(mean)进行了统计，如图14所示。从中可以发现，冰载荷的平均值随冰厚呈明显的增加，而对冰速变化并不敏感。由于离散单元间的粘结强度直接影响结构冰载荷的大小，而颗粒间的粘结强度又与颗粒的粒径相关。因此，冰厚是影响冰荷载的主要因素。另一方面，由于本文主要针对海冰的脆性破坏[33]，单元间的粘结强度并没有考虑加载速率的影响。因此计算结果受冰速的影响不大。这与室内模型试验和现场监测结果相一致[5,34]。
图14 不同冰速、冰厚下冰载荷均值
Fig.14 Mean ice forces with various ice velocities and thicknesses
为验证该耦合模型的合理性，这里将冰载荷峰值的均值与典型锥体静冰力公式进行对比，如图15所示，这里取冰速为0.2 m/s。这里典型锥体静冰力公式分别选用Ralstion、Kato、Yue和 Hirayama[35－38]，即：
式中：FR、FK、FY和FH分别为Ralstion、Kato、Yue和Hirayama锥体静冰力；A1、A2和A3为无量纲系数[36]；Ah和Bh通过模型实验的回归分析得到[37]；σf为海冰的弯曲强度；hi为海冰厚度；ρw为海水密度；ρi为海冰的密度；g为重力加速度；D为锥体直径；DT为锥顶直径；Lc为海冰的破碎长度。
图15 模拟的冰载荷峰值与经典冰载荷公式对比
Fig.15 Ice forces comparing the numerical model and function
从图中发现通过典型公式以及数值计算结果都与冰厚的平方呈线性增加的关系。但每个典型公式和模拟结果存在一定的差异，Ralston计算的冰载荷最大，通过DEM-FEM耦合方法计算的冰载荷最小。这主要是因为海冰和结构的参数选取存在很多限制和差异，从而影响海冰模型的基本假设和数值模拟的结果。同时，DEM 的计算参数，如颗粒层数、刚度以及粘结强度等，也需要进一步的确定。即便如此，当前通过模拟得到的冰载荷是合理的，可用于锥体平台的动力计算和强度分析。
4.2 锥体海洋平台结构的冰激振动
图16 v =0.2m/s,h=0.2m时平台的冰激振动加速度时程
Fig.16 IIVs acceleration of the platform
下面通过渤海 JZ20-2 NW 锥体平台结构的实测振动加速度数据对DEM-FEM耦合模型的合理性进行进一步的验证。图16(a)和图16(b)分别为相同冰况(v=0.2 m/s，h=0.2 m)下实测数据和本文模拟结果的冰激振动加速度时程。实测和数值结果都在-0.15 m/s2～0.15 m/s2范围内变化，具有较好的相似性。图17为模拟得到的冰速、冰厚与锥体平台结构冰激振动加速度均值的关系。结果表明平台结构的冰激振动加速度随冰速、冰厚的增加而增加，即冰速和冰厚是影响平台冰激振动的关键参数。此外，相对于冰速，冰厚对结构冰激振动呈现更明显的影响。
图17 冰速、冰厚与平台结构冰激振动加速度的关系
Fig.17 Relationship between platform vibration acceleration,ice velocity,and ice thickness
分别将不同冰厚下得到的振动加速度峰值与渤海实测数据进行对比分析，如图18所示。图中实心圆点表示实测振动加速度的最大值，空心圆点为通过 DEM-FEM 耦合方法得到的振动加速度峰值。从图中可以发现，模拟结果与实测结果的趋势基本保持一致，即冰厚与振动加速度呈二次非线性增加关系。但是，实测数据和模拟结果在数值方面存在一定的差异，造成这种差异的主要原因可能是海冰的模型参数与实际海冰参数不完全一致。由于现场实测环境复杂，渤海海洋平台只能监测冰速和冰厚两个参数，无法准确得到其他的海冰参数(海冰盐度、冰温和强度等)。而 DEM-FEM 耦合模型的海冰属性是根据相关试验结果得到，这与实际的海冰属性有一定的差异。即便如此，数值模拟的结果基本在实测数据的范围内，且具有相同的趋势。
对实测数据和模拟结果进行频谱分析，如图19所示。图中表明模拟结果与实测数据的功率谱密度(PSD)曲线在频率上基本一致，其主导频率集中在1.0 Hz附近。其中，实测数据PSD峰值对应的频率为1.14 Hz，模拟结果峰值对应的频率为1.01 Hz。考虑到JZ20-2 NW平台的基频为1.07 Hz，因此可以认为该平台的冰激振动主要由1阶模态提供。
4.3 锥体结构的强度分析
通过DEM-FEM耦合模型可得到海冰作用下锥体结构的Von-Mises应力。图20为冰速0.2 m/s，冰厚0.2 m冰况下锥体结构的Von-Mises应力时程。这里提取出最大应力时刻的应力和压力分布云图，如图21所示。图21(a)为结构的Von-Mises应力云图，图21(b)为相对应的法向压力云图。可以发现锥体结构应力分布与压力分布相对应，压力最大处也是应力最大的部分。
图18 JZ20-2NW平台冰激振动现场数据和数值模拟结果
Fig.18 IIV accelerations calculated by the proposed method and field measurements
图19 平台冰激振动加速度现场数据和数值模拟的频谱分析
Fig.19 PSD of the IIV acceleration on the platform via the field data and simulated results
图20 锥体结构Von-Mises应力时程
Fig.20 Von-Mises stress of the cone structure
本文对不同冰厚和冰速下锥体结构应力的最大值(Max)进行了统计分析，结果如图22所示。从中可以发现，冰速对应力的大小影响较小，冰厚则影响较大，这也与冰载荷的变化规律一致。此外，从图中可知本文冰况中最大应力为131.2 MPa并未超过钢的屈服强度(240 MPa)。因此，对于本文计算的所有冰况，线性分析是合理的。
图21 海冰作用下锥体结构的应力及压力云图
Fig.21 Stress and pressure contour of the cone structure
图22 不同冰厚下锥体结构的Von-Mises应力最大值
Fig.22 Maximum value of Von-Mises stress of the cone structure under different ice thicknesses
5 结论
本文提出了基于GPU并行离散元-有限元耦合方法的壳体结构冰激响应计算模型。采用具有粘结-破碎特性的离散单元模拟海冰的破碎过程，抗冰锥体结构采用平板型壳单元构造以分析结构的动力响应和应力分布，建立了离散元与平板型壳单元的耦合算法。同时，发展了 GPU并行环境下两个模型间计算参数的传递算法，通过与LS-DYNA对比验证了耦合方法的精度。基于上述耦合模型对渤海JZ20-2 NW单锥体平台的冰载荷、冰激振动以及锥体应力分布进行了详细讨论。通过与典型冰载荷计算公式对比表明了本文数值结果的合理性；结合渤海实测数据，对结构的冰激振动计算结果进行了对比，可较准确地描述结构冰激振动的基本规律；锥体结构的应力分布与冰载荷的作用位置相对应，且最大应力并未超过结构的屈服强度。
本文提出的离散元-有限元耦合模型不仅可在细观尺度上分析海冰与结构相互作用中的破坏状态，还可揭示结构冰激振动的内在机理，可为冰区结构的冰载荷、结构动力响应分析以及疲劳寿命估计提供一种有效的数值手段。但由于海冰的物理性质十分复杂，本文尚未考虑海冰强度、加载速率等因素。为更加合理地分析海冰与结构相互作用的问题，后续工作将从海冰的本构模型以及流体动力学方面发展相应的耦合模型。
[1]徐田甜,张晓,崔航.我国海上导管架平台抗冰锥体的设计实践[J].船舶工程,2010,32(2)：73－77.Xu Tiantian,Zhang Xiao,Cui Hang.Design practices of anti-ice structure on offshore jacket platform of our country [J].Ship Engineering,2010,32(2)：73－77.(in Chinese)
[2]岳前进,许宁,崔航,等.导管架平台安装锥体降低冰振效果研究[J].海洋工程,2011,29(2)：18－24.Yue Qianjin,Xu Ning,Cui Hang,et al.Effect of adding cone to mitigate ice-induced vibration [J].The Ocean Engineering,2011,29(2)：18－24.(in Chinese)
[3]Yue Q,Zhang L,Zhang W,et al.Mitigating ice-induced jacket platform vibrations utilizing a TMD system [J].Cold Regions Science ＆ Technology,2009,56(2/3)：84－89.
[4]岳前进,刘圆,屈衍,等.抗冰平台的冰振疲劳分析[J].工程力学,2007,24(6)：159－164.Yue Qianjin,Liu Yuan,Qu Yan,et al.Fatigue-life analysis of ice-resistant platforms [J].Engineering Mechanics,2007,24(6)：159－164.(in Chinese)
[5]岳前进,毕祥军,于晓,等.锥体结构的冰激振动与冰力函数[J].土木工程学报,2003,36(2)：17－19.Yue Qianjin,Bi Xiangjun,Yu Xiao,et al.Ice-induced vibration and ice force function of conical structure [J].China Civil Engineering Journal,2003,36(2)：17－19.(in Chinese)
[6]Huang Y,Ma J,Tian Y.Model tests of four-legged jacket platforms in ice：Part 1.Model tests and results [J].Cold Regions Science ＆ Technology,2013,95(11)：74－85.
[7]高良田,王键伟,王庆,等.破冰船在层冰中运动的数值模拟方法[J].工程力学,2019,36(1)：227－238.Gao Liangtian,Wang Jianwei,Wang Qing,et al.Numerical simulation method for motion of the icebreaker in level ice [J].Engineering Mechanics,2019,36(1)：227－238.(in Chinese)
[8]季顺迎,狄少丞,李正,等.海冰与直立结构相互作用的离散单元数值模拟[J].工程力学,2013,30(1)：463－469.Ji Shunying,Di Shaocheng,Li Zheng,et al.Discrete element modelling of interaction between sea ice and vertical offshore structures [J].Engineering Mechanics,2013,30(1)：463－469.(in Chinese)
[9]武文华,于佰杰,许宁,等.海冰与锥体抗冰结构动力作用的数值模拟[J].工程力学,2008,25(11)：192－196.Wu Wenhua,Yu Baijie,Xu Ning,et al.Numerical simulation of dynamic ice action on conical structure [J].Engineering Mechanics,2008,25(11)：192－196.(in Chinese)
[10]Sand B,Fransson L.Nonlinear finite element simulations of ice sheet forces on conical structures [C]// Proc.25th Int.Conf.on Offshore and Arctic Engineering (OMAE),Hamburg,Germany,2006.
[11]周先齐,徐卫亚,钮新强,等.离散单元法研究进展及应用综述 [J].岩土力学,2007(增刊1)：408－416.Zhou Xianqi,Xu Weiya,Niu Xinqiang,et al.A review of distinct element method researching progress and application.Rock and Soil Mechanics,2007(Suppl1)：408－416.(in Chinese)
[12]Hopkins M A.Onshore ice pile-up：a comparison between experiments and simulations [J].Cold Regions Science ＆ Technology,1997,26(3)：205－214.
[13]Polojärvi A,Tuhkuri J,Pustogvar A.DEM simulations of direct shear box experiments of ice rubble：force chains and peak loads [J].Cold Regions Science and Technology,2015,116：12－23.
[14]Ji S,Di S,Liu S.Analysis of ice load on conical structure with discrete element method [J].Engineering Computations,2015,32(4)：1121－1134.
[15]Munjiza A,Lei Z,Divic V,et al.Fracture and fragmentation of thin shells using the combined finite-discrete element method [J].International Journal for Numerical Methods in Engineering,2013,95：479－498.
[16]Chung Y C,Lin C K,Chou P H,et al.Mechanical behavior of a granular of a granular solid and its contacting deformable structure under uni-axial compression-Part1：Joint DEM-FEM modelling and experimental validation [J].Chemical Engineering Science,2016,144：404－420.
[17]Zheng Z M,Zang M Y,Chen H Y,et al.A GPU-based DEM-FEM computational framework for tire-sand interaction simulations [J].Computers and Structures,2018,209(15)：74－92.
[18]Michael M,Vogel F,Peters B.DEM-FEM coupling simulations of the interactions between a tire tread and granular terrain [J].Computer Methods in Applied Mechanics ＆ Engineering,2015,289(47)：227－248.
[19]Shao S,Yan Y,Ji S.Combined Discrete-Finite element modeling of ballasted railway track under cyclic loading[J].International Journal of Computational Methods,2016,14(5)：1750047.
[20]肖军华,张德,王延海,等.基于DEM-FDM耦合的普通铁路碎石道床-土质基床界面接触应力分析[J].工程力学,2018,35(9)：170－179.Xiao Junhua,Zhang De,Wang Yanhai,et al.Study on interface stress between ballast and subgrade for traditional railway based on coupled DEM-FDM [J].Engineering Mechanics,2018,35(9)：170－179.(in Chinese)
[21]邱长林,王菁,闫澍旺.冲击荷载作用下有碎石保护结构的海底管线 DEM-FEM 联合分析研究 [J].岩土工程学报,2015,37(11)：2089－2093.Qiu Changlin,Wang Jing,Yan Shuwang.Coupled DEM-FEM analysis of structure pipelines with rock armor berm under impact load [J].Chinese Journal of Geotechnical Engineering,2015,37(11)：2089－2093.(in Chinese)
[22]Dratt M,Katterfeld A.Coupling of FEM and DEM simulations to consider dynamic deformations under particle load [J].Granular Matter,2017,19(3)：49.
[23]Ji S,Di S,Long X.DEM simulation of uniaxial compressive and flexural strength of sea ice：parametric study [J].Journal of Engineering Mechanics,2016,143(1)：C4016010.(in Chinese)
[24]季顺迎,顾纵棋,王安良,等.基于海洋油气平台上雷达连续监测的海冰漂移特性分析[J].海洋通报,2016,35(3)：280－285.Ji S,Gu Z,Wang A,et al.Analysis of drifting characteristics of the sea ice based on continuous marine radar monitoring on the oil-gas offshore platform[J].Marine Science Bulletin,2016,35(3)：280－285.(in Chinese)
[25]Michel B,Toussaint N.Mechanisms and theory of indentation of ice plates [J].Journal of Glaciology,1977,19(81)：285－300.
[26]Oñate E,Rojek J.Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems [J].Computer Methods in Applied Mechanics ＆ Engineering,2004,193(27)：3087－3128.
[27]Ji S,Wang A,Su J,et al.Experimental studies on elastic modulus and flexural strength of sea ice in the Bohai Sea[J].Journal of Cold Regions Engineering,2011,25(4)：182－195.
[28]Frankenstein G,Garner R.Equations for determining the brine volume sea ice from -0.5℃ to -22.9℃ [J].Journal of Glaciology,1967,6(48)：943－944.
[29]Long Y,Xu Y.Generalized conforming quadrilateral membrane element with vertex rigid rotational freedom[J].Computers ＆ Structures,1994,52(4)：749－755.
[30]岑松,龙志飞.对转角场和剪应变场进行合理插值的厚板元[J].工程力学,1998,15(3)：1－14.Cen Song,Long Zhifei.A Mindlin triangular plate element with improved interpolation for the rotation and shear strain fields [J].Engineering Mechanics,1998,15(3)：1－14.(in Chinese)
[31]Nishiura D,Sakaguchi H.Parallel-vector algorithms for particle simulations on shared-memory multiprocessors[J].Journal of Computational Physics,2011,230：1923－1938.
[32]Qu Y,Yue Q,Bi X,et al.A random ice force model for narrow conical structures [J].Cold Regions Science and Technology,2006,45(3)：148－157.
[33]Timco G W,Weeks W F.A review of engineering properties of sea ice [J].Cold Regions Science and Technology,2010,60：107－129.
[34]Huang Y.Model test study of the nonsimultaneous failure of ice before wide conical structures [J].Cold Regions Science ＆ Technology,2010,63(3)：87－96.
[35]Ralston T D.Plastic limit analysis of sheet ice loads structure [C]// IUTAM Symposium on Physics and Mechanics of Ice,1980：289－308.
[36]Kato K.Simulation of ice load on a conical shaped structure：comparison with experimental results [C]//Proc of the 7th International Offshore and Polar Conference,1997：360－367.
[37]Yue Q J,Qu Y,Bi X J,et al.Ice force spectrum on narrow conical structures [J].Cold Regions Science and Technology,2007,49：161－169.
[38]Hirayama K,Obara I.Ice forces on inclined structure[C]// Proc of the 5th International Offshore Mechanics and Arctic Engineering,1986：515－520.