CN104951601A - 一种海冰-海洋结构相互作用的离散元高性能仿真系统 - Google Patents

Abstract

本发明属于计算机应用技术领域，一种海冰-海洋结构相互作用的离散元高性能仿真系统。海冰-海洋结构相互作用的离散元高性能仿真系统采用三个步骤实现对海冰与海洋平台结构作用的仿真，分别是建立模型，数值计算以及结果显示。用户可以直接在计算机终端设置海冰与海洋结构的参数并其相互作用，准确高效地获得结构的冰荷载和振动响应等力学性能和可视化的三维结果动画。从而使得基于GPU离散元数值分析对冰区海洋平台的结构设计和疲劳分析具有实际的指导性意义。系统适用于环境、气象、水文、工业、科研等多个涉及海冰与海洋结构相互作用的研究领域，界面友好，设计规范，具有良好的用户体验。

Description

一种海冰-海洋结构相互作用的离散元高性能仿真系统

技术领域

本发明属于计算机应用技术领域，涉及一种针对海冰-海洋结构相互作用的数值仿真方法，特别涉及通过数值仿真系统基于GPU的并行高性能数值算法实现海冰与海洋结构耦合作用的大规模离散元计算方法。通过构造平整冰、碎冰、冰脊和堆积冰等不同的冰类型，对导管架平台、自升式平台和浮式海洋平台结构的冰荷载及结构振动响应进行高效、精确的数值分析，为冰区海洋平台的结构设计和疲劳分析提供参考依据。

背景技术

在冰区的油气开发过程中，海冰作为主要的环境荷载对工程结构的安全有着重要的影响。因此，冰区海洋工程结构的设计，需要对冰荷载问题有着明确的了解，例如冰的类型，冰与结构作用时的破坏模式及冰载荷情况等。

在海冰与海洋结构相互作用过程中，冰载荷是影响平台结构震动响应和疲劳寿命的重要因素。海洋结构的冰载荷可以通过现场测量、模型试验、理论分析和数值方法得到。其中，数值方法具备研究成本低、周期短等优势而被广泛应用于海冰与海洋结构作用的冰载荷分析。在数值分析中常用有限元方法计算海冰与海洋结构作用，该方法可以处理相对复杂的力学模型，但不能合理地模拟海冰的动力破坏特性。而受到计算机硬件条件的制约，大型离散元模型很难用CPU并行计算实现。因此，发展高效的大规模离散元并行计算是解决海冰与海洋结构作用的一种有效方法。

发明内容

本发明基于离散元法模拟海冰，计算海冰与海洋结构相互作用，得到冰荷载和结构响应参数，具有计算结果的动态显示和绘制冰力时程曲线等功能，并通过采用GPU并行计算的方式大幅提高计算效率。

本发明的技术方案：

一种海冰-海洋结构相互作用的离散元高性能仿真系统，采用三个步骤实现对海冰与海洋平台结构作用的仿真，分别是建立模型，数值计算以及结果显示，具体如下：

(一)建立模型：建立的模型包括海冰的离散元模型和海洋结构模型

海冰的离散元模型：海冰分为平整冰、浮冰和冰脊，海冰的主要参数包括海冰区域的长、宽、厚度、位置、温度、盐度、压缩强度、弯曲强度、冰速、水面高度、海冰密度和海水密度。浮冰的形状有多边形、圆形和矩形，其他参数包括密集度和单块浮冰的平均尺寸。

海冰的离散元模型由球型颗粒通过六边形的排列方式构成，球型颗粒的参数包括粒径大小、杨氏模量、颗粒间的摩擦系数、颗粒与结构摩擦系数、颗粒与结构的回弹系数、颗粒间切向和法向的刚度比和阻尼比、宏观冰的弯曲强度与微观颗粒法向粘结强度比以及宏观冰的压缩强度与微观颗粒切向粘结强度的比。

海洋结构模型：包括直立腿、椎体和复杂结构，其中复杂结构需要调入。海洋结构模型的参数包括桩腿直径、高度、个数、中心位置、质量、刚度、阻尼系数、转动惯量、转动刚度和转动阻尼系数。

仿真模型中海冰离散元模型和海洋结构模型的参数根据海冰和海洋结构的尺寸和力学性能设定，以此建立海冰与海洋结构的数值模型。

(二)数值计算：在采用离散单元模型计算海冰的动力作用及破碎过程中，将海冰离散为具有一定质量和大小的颗粒单元，颗粒间具有相应的粘接作用。这里对海冰离散元模型中的接触力模型、粘接模型及粘接强度的设定进行介绍。

(a)颗粒间的接触力模型

在两个颗粒相互碰撞过程中，考虑颗粒间因相对速度和弹性变形而引起的黏弹性作用力，采用Mohr-Coulomb摩擦定律计算接触颗粒之间的接触力，如图1所示。其中，M_A和M_B为颗粒A和B的质量，K_n和K_s分别是法向和切向刚度系数，C_n和C_s是法向和切向阻尼系数，μ是摩擦系数。

在线性接触模型中，两个颗粒间的作用力包括弹性力和粘滞力两部分，又可分为法向力和切向力，按下式计算：

$F_n=K_n{x}_n-C_n{\stackrel{\·}{x}}_n---\left(1\right)$

式中，x_n和分别为颗粒的法向变形和应变率。

两个颗粒间的切向力也由弹性和粘滞两部分组成，且满足Mindline理论和Mohr-Coulomb摩擦定律，切向接触力为

$F_s^{*}=K_s{x}_s-C_s{\stackrel{\·}{x}}_s---\left(2\right)$

$F_s=\min (F_s^{*},\mathrm{sign}\left(F_s^{*}\right)\mathrm{\μ}{F}_n)---\left(3\right)$

式中，x_s和分别为颗粒的切向变形和应变率。

两个颗粒碰撞的法向有效刚度系数为

$K_n=\frac{2k_n^A{k}_n^B}{k_n^A+k_n^B}---\left(4\right)$

式中，和分别为球单元A和B的刚度系数。

法向阻尼系数按下式计算，即

$C_n={\mathrm{\ζ}}_n\sqrt{2{\mathrm{MK}}_n}---\left(5\right)$

这里无量纲法向阻尼系数为

${\mathrm{\ζ}}_n=\frac{-\ln e}{\sqrt{{\mathrm{\π}}^2+{\ln }^{2}e}}---\left(6\right)$

式中，M为两颗粒单元的有效质量，e为回弹系数。切向和法向刚度、阻尼系数有如下关系：K_s＝αK_n，C_s＝βC_n，这里取α＝0.5，β＝0.0。

在线性接触模型中，计算步长一般取时间步长为二元接触时间的1/50。该二元接触时间定义为

$T_{\mathrm{bc}}=\frac{\mathrm{\π}}{\sqrt{\frac{{2K}_n}{M}(1-{{\mathrm{\ζ}}_n}^2)}}---\left(7\right)$

式中，T_bc为二元接触时间,即两个球单元从碰撞到分离的接触时间。在线粘弹性模型中，它是一个与颗粒大小和材料性质相关的常数。

(b)平行粘接模型

平行粘接就是在两个粘接颗粒间设定一个弹性粘接圆盘，如图2所示。圆盘可以传递两个颗粒间的粘接力和力矩，即拉力、剪力、弯矩和扭矩。并且以上粘接力F和力矩M都可以分解为法向分量和切向分量

F＝F_n+F_s  (8)

M＝M_n+M_s  (9)

式中和分别是力和力矩的法向分量和切向分量。

在平行粘接模型中，粘接圆盘上的最大拉应力σ_max和剪应力τ_max依据梁的拉伸、扭转和弯曲理论有

${\mathrm{\σ}}_{\max }=\frac{\left|F_n\right|}{A}+\frac{\left|M_s\right|}{I}R---\left(10\right)$

${\mathrm{\τ}}_{\max }=\frac{\left|F_s\right|}{A}+\frac{\left|M_n\right|}{J}R---\left(11\right)$

式中，A、J和I分别为平行接触圆盘的面积、极惯性矩和惯性矩，有A＝πR²，J＝1/2πR⁴，J＝1/4πR⁴，其中R为粘接圆盘的半径。当最大拉应力和剪应力超过其拉伸强度和剪切强度时粘接颗粒将断开。

(c)粘结强度的确定

在海冰材料的离散单元模型中，颗粒的粘接强度是影响计算结果的关键参数。考虑计算中颗粒的大小对粘接力的影响，这里设定颗粒间的粘接强度为

${\mathrm{\σ}}_b=\frac{F_n^c}{A}---\left(12\right)$

式中，σ_b为粘接强度，为粘接单元破碎时的最大拉力。这里取粘接面积为A＝πR²。

试验结果表明，海冰单轴压缩强度是海冰卤水体积(温度、盐度)、加载速率等因素的函数^[6,18]。在主要考虑卤水体积影响的情况下，可将颗粒间的粘接强度用颗粒间最大粘接强度表示，即

${\mathrm{\σ}}_b=\mathrm{\β}\left(v_b\right){\mathrm{\σ}}_b^{\max }---\left(13\right)$

式中，β(v_b)为卤水体积v_b影响下的海冰强度折减系数。颗粒间的最大粘接强度可通过海冰单轴压缩强度的敏度分析进行确定。考虑海冰的压缩和弯曲强度与卤水体积有相似的对应关系^[19]，则有

$\mathrm{\β}=e^{-4.29\sqrt{v_b}}---\left(14\right)$

式中，v_b可设为海冰温度和盐度的函数，即

式中，T为海冰温度(℃)，S为海冰盐度(‰)。

数值计算过程主要包括颗粒搜索、接触判断、内力叠加和位置更新四个步骤。首先根据颗粒的位置，判断颗粒与颗粒之间的接触对，计算每个颗粒接触对之间的接触力及颗粒与结构的接触力，再将上述的接触力叠加到每个颗粒和结构上，通过接触力计算颗粒和结构的速度和位移，并更新颗粒和海冰结构的位置。在每个时间间隔内循环上述过程，直到计算时间达到预设的计算总时长，循环停止。根据预设的结果输出频率，记录颗粒和结构的速度、加速度、位移和内力。

其中，颗粒间的接触判断采用网格单元法(Uniform Grid)，首先对颗粒所占空间划分网格单元，当颗粒被分配到某一个单元内时，只有在同一单元或者直接相邻的单元内的颗粒之间才进行接触判断。该算法基于GPU并行计算环境设计，与传统的网格法相比，该网格法允许一个网格内有多个颗粒。每个网格的搜索范围是以该网格为中心的3×3×3个网格，所以网格大小不能小于颗粒直径。这种方法有效地提高了离散元法中颗粒接触判断的效率，更加适合大运算量的并行环境。

(三)结果显示：海冰与海洋结构相互作用的计算输出结果包括海冰颗粒的速度、加速度、位移和颗粒间的接触力以及结构的速度、加速度、转速、位移和颗粒对结构的作用力。按时间顺序显示颗粒和结构的计算结果，可以得到海冰与海洋结构相互作用过程中的破坏模式和运动规律。此外还可得到冰力时程曲线，即海冰对结构的作用力随时间变化的曲线。

本发明的效果和益处是，用户可以直接在计算机终端设置海冰与海洋结构的参数并其相互作用，准确高效地获得结构的冰荷载和振动响应等力学性能和可视化的三维结果动画。从而使得基于GPU并行算法的离散元数值分析对冰区海洋平台的结构设计和疲劳分析具有实际的指导性意义。系统适用于环境、气象、水文、工业、科研等多个涉及海冰与海洋结构相互作用的研究领域，界面友好，设计规范，具有良好的用户体验。

附图说明

图1是两个颗粒单元间的接触力模型。

图2是两个颗粒单元间的平行粘结模型。

图3是本发明的海冰与海洋平台相互作用的程序流程框图。

具体实施方式

以平整冰海冰与单个直立腿结构相互作用为例，首先打开前处理模块(Preprocess)设置参数，点击海冰(Sea Ice)选择平整冰(Level)，海冰的长宽为10m×10m，冰厚0.2m，冰层数为1，在边界条件(Boundary Condition)中设置边界速度0.2m/s，边界与颗粒的刚度比为1，在强度(Strength)设置中压缩强度为2.4MPa，弯曲强度为1.5MPa，海冰和海冰的密度分别为920Kg/m³和1035Kg/m³，水深3m，冰面右边缘的位置(-1,0)，水流速度0.2m/s。然后点击结构(Structure)选择圆柱直立腿结构(Cylinder)，桩腿的直径为0.75m，高度为6m，点击桩腿数目及位置选项(Multi_leg Distribution)选择1个桩腿，设置桩腿下底面中心位置为(0,0,0)。为了显示桩腿上压力分布可以对桩腿设置局部压力网格(Local Pressure)将桩腿沿轴向和周向方向划分若干等份，局部网格的高度为2m，轴线的数量是10，周向上每个扇形等份的角度为30°。然后设置桩腿的质量中心(0,0,3)，质量1.6×10⁷，刚度7.5×10⁶，阻尼系数0.07，转动惯量4.0×10⁹，转动刚度1.0×10⁸，转动阻尼系数0.07。

点击“Element”设置离散单元颗粒性质，弹性模量1.0×10⁹，颗粒间摩擦系数0.1，颗粒间回弹系数为0.3，颗粒与结构间摩擦系数为0.2，颗粒与结构回弹系数0.3，颗粒的法向与切向刚度比0.1，法向和切向的阻尼比0.5，微观拉伸强度与宏观弯曲强度的比0.2，微观剪切强度与宏观压缩强度的比0.32。

完成模型参数的输入后，点击预览(Preview)，检查建立的海冰与结构模型，确定无误后点击完成(Finish)，继续点击计算(Calculate)开始数值计算。计算开始前要设置计算时间长度(Total Time)25s，输出显示动画的次数(PrintStep)30，输出冰载荷的频率(Force Frequency)100Hz。继续点击计算(Calculate)，弹出程序运行窗口。根据计算机配置选择显卡，单一显卡时输入“0”点回车键，程序开始运行，可在窗口中查看程序运行进度。程序运行中切勿关闭窗口，待计算结束时窗口自动关闭。

点击“后处理”按钮，弹出输入动画文件的对话框，选择计算得到的三维动画数据ANIMATION_3D.DAT，具体文件格式要求，可以在主界面的“帮助”选项中查看。调入动画文件后，点击开始按钮查看动画。通过调节动画播放控制按钮，实现动画播放，动画复位，动画快退，动画暂停，动画快进，保存当前画面图片(jpg，bmp格式)，保存动画视频(avi格式)，动画自动旋转，冰载荷结果时程曲线，结构位移时程曲线及结构振动时程曲线。在“Colors”中可以选择颗粒颜色代表的不同物理意义包括速度，颗粒直径，内力，切向力和法向力等，点击颜色条，更换颜色。颗粒颜色代表数值的大小，颜色越偏向红色，所代表的数值越大；反之，越偏向蓝色表示数值越小。“Scaling”调节视图的位置和画面的放大缩小，“Timestep”调节动画的播放速度，“BeginningAngle”是动画播放时的初始角度。同时还可通过单击鼠标左键可以在调节图像位置区域或直接拖拉图像，选择合适的观察角度。由此实现对海冰与海洋结构相互作用破碎过程和冰载荷的数值模拟。

整个软件系统可视为由模型建立，数值计算，结果显示等三大应用模块构成，采用上述的离散元算法，为用户提供了综合型集成化的海冰与海洋结构相互作用的分析平台。

Claims (1)

1.一种海冰-海洋结构相互作用的离散元高性能仿真系统，其特征在于，包括建立模型、数值计算和结果显示，步骤如下：

(一)建立模型：建立的模型包括海冰的离散元模型和海洋结构模型

海冰的离散元模型：海冰分为平整冰、浮冰和冰脊，海冰的参数包括海冰区域的长、宽、厚度、位置、温度、盐度、压缩强度、弯曲强度、冰速、水面高度、海冰密度和海水密度；浮冰的形状有多边形、圆形和矩形；还包括密集度和单块浮冰的平均尺寸；

海冰的离散元模型由球型颗粒通过六边形的排列方式构成，球型颗粒的参数包括粒径大小、杨氏模量、颗粒间的摩擦系数、颗粒与结构摩擦系数、颗粒与结构的回弹系数、颗粒间切向和法向的刚度比和阻尼比、宏观冰的弯曲强度与微观颗粒法向粘结强度比以及宏观冰的压缩强度与微观颗粒切向粘结强度的比；

海洋结构模型：包括直立腿、椎体和复杂结构，其中复杂结构需要调入；海洋结构模型的参数包括桩腿直径、高度、个数、中心位置、质量、刚度、阻尼系数、转动惯量、转动刚度和转动阻尼系数；

仿真模型中海冰的离散元模型和海洋结构模型的参数根据海冰和海洋结构的尺寸和力学性能设定，以此建立海冰与海洋结构的数值模型；

(二)确定海冰的离散元模型的接触力模型、平行粘接模型及粘接强度

(a)接触力模型

在两个颗粒相互碰撞过程中，考虑颗粒间因相对速度和弹性变形而引起的黏弹性作用力，采用Mohr-Coulomb摩擦定律计算接触颗粒之间的接触力；

在接触力模型中，两个颗粒间的作用力F_n包括弹性力和粘滞力两部分，分为法向力和切向力，按下式计算：

$F_n=K_n{x}_n-C_n{\stackrel{\·}{x}}_n---\left(1\right)$

式中，M_A和M_B为颗粒A和B的质量，K_n和K_s分别是法向和切向刚度系数，C_n和C_s是法向和切向阻尼系数，μ是摩擦系数；x_n和分别为颗粒的法向变形和应变率；

两个颗粒间的切向力由弹性和粘滞两部分组成，且满足Mindline理论和Mohr-Coulomb摩擦定律，切向接触力为

$F_s^{*}=K_s{x}_s-C_s{\stackrel{\·}{x}}_s---\left(2\right)$

$F_s=\min (F_s^{*},\mathrm{sign}\left(F_s^{*}\right)\mathrm{\μ}{F}_n)---\left(3\right)$

式中，x_s和分别为颗粒的切向变形和应变率；

两个颗粒碰撞的法向有效刚度系数K_n为

$k_n=\frac{2k_n^A{k}_n^B}{k_n^A+k_n^B}---\left(4\right)$

式中，和分别为颗粒A和B的刚度系数；

法向阻尼系数C_n按下式计算，即

$C_n={\mathrm{\ζ}}_n\sqrt{{2\mathrm{MK}}_n}---\left(5\right)$

这里无量纲法向阻尼系数ζ_n为

${\mathrm{\ζ}}_n=\frac{-\ln e}{\sqrt{{\mathrm{\π}}^2+{\ln }^{2}e}}---\left(6\right)$

式中，M为两颗粒单元的有效质量的和，e为回弹系数；切向和法向刚度、阻尼系数有如下关系：K_s＝αK_n，C_s＝βC_n，这里取α＝0.5，β＝0.0；

在线性接触模型中，计算步长取时间步长为二元接触时间的1/50；该二元接触时间T_bc定义为

$T_{\mathrm{bc}}=\frac{\mathrm{\π}}{\sqrt{\frac{{2K}_n}{M}(1-{{\mathrm{\ζ}}_n}^2)}}---\left(7\right)$

式中，T_bc为二元接触时间，即两个颗粒从碰撞到分离的接触时间；在线粘弹性模型中，它是一个与颗粒大小和材料性质相关的常数；

(b)平行粘接模型

平行粘接模型：在两个粘接颗粒间设定一个弹性粘接圆盘，圆盘传递两个颗粒间的作用力和力矩，即拉力、剪力、弯矩和扭矩；力和力矩分解为法向分量和切向分量

${\stackrel{\→}{F}}_i={{\stackrel{\→}{F}}_i}^s+{{\stackrel{\→}{F}}_i}^n---\left(8\right)$

${\stackrel{\→}{M}}_i={{\stackrel{\→}{M}}_i}^s+{{\stackrel{\→}{M}}_i}^n---\left(9\right)$

式中，和分别是法向分量和切向分量；

在平行粘接模型中，粘接圆盘上的最大拉应力和剪应力依据梁的拉伸、扭转和弯曲理论：

${\mathrm{\σ}}_{\max }=\frac{-{\stackrel{\‾}{F}}^n}{A}+\frac{\left|{\stackrel{\‾}{M}}_i^s\right|}{I}\stackrel{\‾}{R}---\left(10\right)$

${\mathrm{\τ}}_{\max }=\frac{\left|{\stackrel{\‾}{F}}_i^s\right|}{A}+\frac{\left|{\stackrel{\‾}{M}}^n\right|}{J}\stackrel{\‾}{R}---\left(11\right)$

式中，A、J和I分别为平行接触圆盘的面积、极惯性矩和惯性矩，有 其中为粘接圆盘的半径；当最大拉应力和剪应力超过其拉伸强度和剪切强度时粘接单元将断开；

(c)粘结强度

在海冰材料的离散单元模型中，颗粒间的粘接强度是影响计算结果的关键参数；考虑颗粒的大小对粘接力的影响，设定颗粒间的粘接强度为

${\mathrm{\σ}}_b=\frac{F_n^c}{A}---\left(12\right)$

式中，σ_b为粘接强度，为粘接单元破碎时的最大拉力，粘接面积为A＝πR²；

将颗粒间的粘接强度用颗粒间最大粘接强度表示，即

${\mathrm{\σ}}_b=\mathrm{\β}\left(v_b\right){\mathrm{\σ}}_b^{\max }---\left(13\right)$

式中，β(v_b)为卤水体积v_b影响下的海冰强度折减系数，颗粒间的最大粘接强度通过海冰单轴压缩强度的敏度分析确定；考虑海冰的压缩和弯曲强度与卤水体积有相似的对应关系，则有

$\mathrm{\β}=e^{-4.29\sqrt{v_b}}---\left(14\right)$

式中，v_b可设为海冰温度和盐度的函数，即

$v_b=s(0.532+\frac{49.185}{\left|T\right|})$ -0.5℃≥T≥-22.9℃       (15)

式中，T为海冰温度℃，S为海冰盐度‰；

确定海冰离散元模型的接触力模型、粘接模型及粘接强度包括颗粒搜索、接触判断、内力叠加和位置更新四个步骤；首先根据颗粒的位置，判断颗粒与颗粒之间的接触对，计算每个颗粒接触对之间的接触力及颗粒与结构的接触力，再将上述的接触力叠加到每个颗粒和结构上，通过接触力计算颗粒和结构的速度和位移，并更新颗粒和海冰结构的位置；在每个时间间隔内循环上述过程，直到计算时间达到预设的计算总时长，循环停止；根据预设的结果输出频率，记录颗粒和结构的速度、加速度、位移和内力；

其中，颗粒间的接触判断采用网格单元法，首先对颗粒所占空间划分网格单元，当颗粒被分配到某一个单元内时，只有在同一单元或者直接相邻的单元内的颗粒之间才进行接触判断；该算法基于GPU-CUDA并行计算环境设计，该网格法允许一个网格内有多个颗粒；每个网格的搜索范围是以该网格为中心的3×3×3个网格，网格大小不能小于颗粒直径；

(三)结果显示：海冰与海洋结构相互作用的计算输出结果包括海冰颗粒的速度、加速度、位移和颗粒间的接触力以及结构的速度、加速度、转速、位移和颗粒对结构的作用力；按时间顺序显示颗粒和结构的计算结果，得到海冰与海洋结构相互作用过程中的破坏模式和运动规律；此外，还得到冰力时程曲线，即海冰对结构的作用力随时间变化的曲线。