一种浮式结构的全耦合CFD-系泊非线性分析方法

技术领域

本发明涉及预测系泊浮式结构的大幅运动响应技术领域，具体涉及一种浮式结构的全耦合CFD-系泊非线性分析方法。

背景技术

海洋环境中浮式结构的动力响应一直是海洋工程领域关注的热点。波浪荷载的冲击影响其稳定性和使用寿命。此外，浮式结构的剧烈运动可能导致系统故障，而设备维修和部件更换会大大增加运营成本。因此，准确预测浮体结构在波浪中的运动响应是至关重要的。

目前，波浪荷载作用下浮体结构运动响应的耦合分析主要采用势流和计算流体力学CFD数值方法。其中，基于势流理论的数值方法由于计算效率高的优点，早已被用于研究波浪中浮体结构的动力响应。然而，这些方法不能准确地考虑到粘性和非线性效应。随着计算机技术的飞速发展，人们发展了许多基于非线性Navier-Stokes方程的CFD方法，这些方法由于固有地考虑了流体的黏度，因而更加精确，能够提供丰富的流场信息。因此，CFD方法越来越广泛地应用于浮式系统的运动响应研究。系泊系统对浮式平台的稳定性和安全性起着至关重要的作用。系泊索在使用寿命期间，如果受到极端风浪荷载的作用而突然断裂，会导致浮式平台剧烈运动，严重威胁浮式平台的运行稳定性和安全性。因此，准确预测系泊索的动力响应是至关重要的。常用的系泊分析模型可分为静态、准静态和动态三种。通过考虑系泊线运动方程中的动力效应，动力分析模型能够比静力和准静力模型更准确地预测系泊载荷，从而预测浮式结构的响应。近年来，一些研究者尝试将动力系泊分析模型与CFD求解器相结合，用于系泊浮式结构的响应预测。

但是研究人员建立并验证的系泊浮式结构的耦合分析模型，大多是针对正常海况下的结构响应。然而，系泊浮式结构在使用过程中，当共振现象发生或受到极端波浪荷载时，会激发出较大的非线性运动响应。在这种危险海况下，系泊浮式结构的生存能力研究涉及到强烈的非线性系泊动力学和显著的结构平动和旋转响应，需要高保真耦合CFD-系泊模型。为了适应浮式结构在CFD模拟中的响应，必须对计算网格进行更新。最常用的网格更新技术是动态网格变形方法，该方法采用网格变形方法，不改变网格拓扑结构，但当涉及大结构响应时，可能导致网格质量下降甚至模拟失效。然而重叠网格技术可以克服这个困难，进而才能准确的模拟并预测浮式结构的大幅运动响应。

发明内容

为解决现有技术存在的不足，本发明针对系泊浮式结构的大幅运动，首次提出了一种浮式结构的全耦合CFD-系泊非线性分析方法。首先在开源CFD框架OpenFOAM下，开发了一种基于允许大轴向延伸的细长杆理论的非线性有限元系泊动力分析方法，并创造性的将其与基于重叠网格的两相流RANS求解器和刚体六自由度运动求解器相结合，以实现对恶劣海况下浮式结构的显著运动响应准确的模拟和预测。

本发明的技术方案为：

本发明提供了一种浮式结构的全耦合CFD-系泊非线性分析方法，该方法包括如下步骤：

S1、构建非线性有限元系泊动力学模块：基于细长杆理论的动力系泊分析方法在OpenFOAM中以六自由度刚体约束的形式实现，并以动态库的形式编译成单独的模块，以便于与不同的流体求解器进行耦合；

S2、系泊模块与重叠网格CFD求解器耦合：将具有OpenFOAM重叠网格功能的两相流求解器overInterDyMFoam与开源波浪生成库waVes2Foam相结合，并耦合已开发的基于细长杆理论的动力系泊分析方法，建立新的波-结构-相互作用耦合求解器，并命名为overWaveDyMFoam。

进一步的，步骤S1还包括：

在三维直角坐标系下建立了允许大轴向拉伸的细长杆模型的运动方程，变形杆的中心线用空间曲线r(s，t)表示，空间曲线r(s，t)是弧长s和时间t的函数，空间曲线中的t、n、b分别为切线方向、法线方向和副法线方向上的单位矢量，e

允许大轴向拉伸的细长杆的运动方程和约束方程分别定义为：

式(1)、(2)中，变量上方的点表示对时间t的导数，撇号表示对弧长s的导数；q为作用在系泊索上的外力；M为质量矩阵；

以上变量的表达式分别为：

M＝ρ

其中，I是单位矩阵；

T和N是变换矩阵；

变换矩阵T和N分别定义为：

N＝I-T

作用在系泊索上的外力q包括重力、静水压力和水动力；水动力部分包括惯性力、拖曳力和Froude-Krylov力海水，其中拖曳力由莫里森方程计算；

外力q的表达式为：

式中，C

ρ

D

A

ρ

A

v

a

g为为重力加速度；

将Galerkin数值分析方法用于求解细长杆理论的控制方程(即式(1)的运动方程和式(2)的约束方程)；引入Hermite三次形状函数a

三次Hermite形函数a

二次Hermite形函数p

其中，ξ是无量纲量ξ＝s/L，L是变形前的单元长度；

因此，坐标系中单元位移矢量r、拉格朗日乘子

其中，n＝1～3，i＝1～4和m＝1～3，方程中用形函数描述的主要参数分别为：

其中，u指单元节点位移；

为了空间上离散系泊系统的运动控制方程，将式(1)的运动方程右边的项移到左边得到：

应用Galerkin数值分析方法将式(14)运动方程简化为一组常微分方程：第一，方程两端乘以形函数a

第二，再对上式(15)进行简化并分部积分可得：

式(16)方程右端是单元端点的力，该力可以用广义力f

将式(8)、式(6)和式(7)代入式(16)，可得细长杆单元的一组常微分方程：

其中i，k＝1～4，j，l，m，n＝1～3；

同理，将式(8)、式(6)和式(7)代入式(2)，可得约束方程的表达式：

其中i，k＝1～4，j，l，m，n＝1～3；

其中f

L

对于第一个和最后一个单元，应施加结构的边界条件。

进一步的，步骤S2还包括：

用静力分析方法初始化系泊系统后，将浮式结构的位置和速度从六自由度运动解算器传递到导缆孔；

动态系泊分析求解器更新系泊线的状态，包括张力和节点位置；

将系泊力作用于浮式结构，并将其返回到六自由度运动求解器中，以更新浮式结构的运动响应。

本发明所达到的有益效果为：

第一、本发明在OpenFOAM框架下基于细长杆理论开发非线性有限元系泊动力求解模块，该系泊动态求解方法的显著优势在于：能够处理系泊轴向大变形，能考虑系泊的非线性效应以及计算时能快速收敛，率先解决了常规系泊动力求解方法无法处理的极端环境下系泊轴向大变形的问题。

第二、本发明首次将自主开发并利用有限元方法求解的系泊系统动力求解模块与重叠网格耦合，集成waVes2F

第三、本发明中开发的求解器能准确预测系泊浮式结构物与极端波相互作用的运动响应以及作用在系泊缆和浮式结构物上的荷载，有利于恶劣海况下系泊浮式结构的生存能力研究。

附图说明

图1是本发明的三维直角坐标系下细长杆模型。

图2是本发明的耦合求解器overWaveDyMFoam流程图。

图3是本发明的实验与模拟纵荡的比较RAO。

图4是不规则波中运动响应能谱的比较。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1～4所示，本发明提供了一种浮式结构的全耦合CFD-系泊非线性分析方法，包括如下步骤：

步骤1、非线性有限元系泊动力学模块建立：

基于细长杆理论的动力系泊分析方法在OpenFOAM中以六自由度刚体约束的形式实现，并以动态库的形式编译成单独的模块，便于与不同的流体求解器进行耦合。

步骤2、系泊模块与重叠网格CFD求解器耦合：

将具有OpenFOAM重叠网格功能的两相流求解器overInterDyMFoam与开源波浪生成库waves2Foam相结合，该库采用松弛区技术有效避免了出口边界处的波浪反射，并耦合已开发的基于细长杆理论的动力系泊分析方法，建立新的波-结构-相互作用耦合求解器，并命名为overWaveDyMFoam。

步骤1具体还包括：

在三维直角坐标系下建立了允许大轴向拉伸的细长杆模型的运动方程。如图1所示，变形杆的中心线可以用空间曲线r(s，t)表示，空间曲线r(s，t)是弧长s和时间t的函数。图1中，t、n、b分别为切线方向、法线方向和副法线方向的单位矢量，e

允许大轴向拉伸的细长杆的运动方程和约束方程分别定义为：

式(1)、(2)中，变量上方的点表示对时间t的导数，撇号表示对弧长s的导数；q为作用在系泊索上的外力；M为质量矩阵；

以上变量的表达式分别为：

M＝ρ

其中，I是单位矩阵；

T和N是变换矩阵；

变换矩阵T和N分别定义为：

N＝I-T

作用在系泊索上的外力q包括重力、静水压力和水动力；水动力部分包括惯性力、拖曳力和Froude-Krylov力海水，其中拖曳力由莫里森方程计算；

外力q的表达式为：

式中，C

ρ

D

A

ρ

A

v

a

g为重力加速度；

将Galerkin数值分析方法用于求解细长杆理论的控制方程(即式(1)的运动方程和式(2)的约束方程)；引入Hermite三次形状函数a

三次Hermite形函数a

二次Hermite形函数p

其中，ξ是无量纲量ξ＝s/L，L是变形前的单元长度；

因此，坐标系中单元位移矢量r、拉格朗日乘子

其中n＝1～3，i＝1～4和m＝1～3，方程中用形函数描述的主要参数分别为：

其中，u指单元节点位移；

为了空间上离散系泊系统的运动控制方程，将式(1)的运动方程右边的项移到左边得到：

应用Galerkin数值分析方法将式(14)的运动方程简化为一组常微分方程：

(1)方程两端乘以形函数a

(2)再对上式(15)进行简化并分部积分可得：

式(16)方程右端是单元端点的力，该力可以用广义力f

将式(8)、式(6)和式(7)代入式(16)，可得细长杆单元的一组常微分方程：

其中i，k＝1～4，j，l，m，n＝1～3；

同理，将式(8)、式(6)和式(7)代入式(2)，可得约束方程的表达式：

其中i，k＝1～4，j，l，m，n＝1～3；

其中f

L

对于第一个和最后一个单元，应施加结构的边界条件；如果三维空间的一根系泊线有N个单元，则计算模块中应有15+8(N-1)个独立方程和系数。

步骤1的优点在于：首次在OpenFOAM中开发了非线性有限元系泊动力求解模块，弥补了原生OpenFOAM没有系泊求解模块的缺陷，为后续建立系泊浮式结构耦合模型提供基础。

步骤2具体还包括：

首先，系泊模块与浮式结构的耦合具体实现如下：

(1)用静力分析方法初始化系泊系统后，将浮式结构的位置和速度从六自由度运动解算器传递到导缆孔；

(2)然后，动态系泊分析求解器更新系泊线的状态，包括张力和节点位置；

(3)随后，将系泊力作用于浮式结构，并将其返回到六自由度运动求解器中，以更新浮式结构的运动响应。

其次，对系泊浮式结构耦合模型进行实验验证：

平台的自由衰减响应可以反映一定自由度下的自然周期。在本发明中，对系泊浮式结构最重要的三个自由度，即纵荡、垂荡和纵摇，进行了实验和数值上的自由衰减试验。三个自由度的数值预测自然周期与实验结果的相对差异均小于3％。另外，还模拟了规则波与半潜式平台的相互作用，并将数值预测与实验数据进行了比较。计算域：波的产生/传播/松弛区长度分别等于波长的1/0.5/2倍。对比五种不同工况下系泊平台的预测和实测运动响应，发现系泊平台在规则波作用下表现出周期或准稳态的等幅运动响应，用该耦合模型得到的平台的三个响应与实验测量结果吻合较好。

模型实验验证结果表明：数值预测的结果与实验结果的相对差异均小于5％。

最后，将验证过的系泊浮式结构耦合模型用于系泊浮式结构大振幅运动响应预测：

采用JONSWAP谱产生不规则波，形状参数设为γ＝3.3。显著波高0.12m，谱峰周期1.6s。为了保证数值模型能够准确地产生不规则波，首先进行了不考虑任何结构的单波模拟。模拟时间为谱峰周期的160s或100倍，距离波发生区0.5m处设置测波仪。结果表明数值模拟得到的有效波高与实验测量结果相当。

然后进行能谱分析，在频域中比较数值和实验数据。测量的波表面高程计算的能谱呈钟形，大部分能量分布在0.625Hz附近或规定的谱峰周期为1.6s附近，表明模拟能谱与实验能谱吻合较好。进一步提取出定义频谱的几个重要特征参数，包括显著波高、谱峰周期和谱峰。实验数据与数值数据的比较表明，实测值与模拟值的光谱特征参数差异小于3％。对0～1.4Hz频率范围内的能谱密度积分得到了不规则波的总能量，实验与模拟的差异不超过5％。因此，目前的数值波槽能够以较好的精度产生所需的不规则波。

将验证的不规则波应用于系泊浮式结构，将其模拟运动响应与实验数据进行对比。为排除初始扰动的影响，选取40-200s时间范围内的响应数据进行FFT分析，预测能谱和实验能谱各出现两个不同的谱峰。其中一个峰对应于入射波频率(0.625Hz)，另一个接近三自由度半潜式平台的固有频率。纵荡运动和纵摇运动的共振能量明显大于分布在波(主)频率附近的能量，说明显著的响应幅值出现在共振频率处，不可低估。而对于垂荡运动，两个峰的能量分布是相似的。

提取了纵荡、垂荡和纵摇运动的几个重要特征参数，包括显著响应高度、谱峰周期和谱峰。三个自由度的实验响应和数值响应的特征参数比较，最大误差不超过6.4％。同时，在0～1.4Hz的频率范围内对三个自由度的能谱密度进行积分，得到了每个自由度的总能量，实验与仿真的差异不超过8％。纵荡、垂荡和纵摇运动的共振谱峰值能量与实验相比相差不超过10％，考虑到实验中影响测量精度的因素众多，这是可以接受的。

对于采用悬链线系泊的浮体，系泊系统的刚度对其纵荡响应起着关键作用。当浮式结构进行大幅度的水平运动时，系泊线可能被拉伸到一定程度，使系泊刚度变为非线性。垂荡和纵摇运动的共振频率分别在0.46Hz和0.19Hz左右，对应于其固有频率，且这两个方向的运动受系泊系统的影响较小。

预测结果表明：所建立的耦合数值模型不仅能准确地捕捉到系泊浮体结构在纵荡、垂荡和纵摇运动下的波频响应，而且还能捕捉到系泊浮体结构在三个自由度下的共振响应。

步骤2的优点在于：首次在OpenFOAM中开发了运用重叠网格技术的两相流耦合求解器overWaveDyMFoam，并耦合步骤1中已开发的非线性有限元系泊动力求解模块；开发的数值耦合方法精度相比传统的方法有显著提高(传统方法一般可接受误差为15％以内)；所建立的耦合数值模型解决了传统动网格CFD方法对于系泊浮式结构大幅运动计算发散的问题，同时，验证结果充分展示了全非线性耦合方法在强非线性海洋环境下预测系泊浮式结构运动响应的能力。

模型信息

下面采用一个具体的算例模型来做说明。本算例采用半潜式模型的缩放比设为λ＝1:100，主要参数见表1。波浪槽水深1m，模型吃水0.251m，模型质心距水面0.013m。坐标系原点定义在自由曲面上模型的中心。半潜式模型锚泊系统由8条悬链线组成，每条悬链线由钢链和弹簧组成，锚泊参数见表2。系泊线对称布置，各角处两条系泊线与x轴夹角分别为37.5°和60°。

表1半潜模型主尺度参数

表2系泊参数

首先根据图1展示的坐标系按照步骤1在OpenFOAM中开发有限元系泊动力求解模块；其次是按照图2的流程图完成步骤2：开发系泊浮式结构物耦合求解器；再验证系泊浮式结构耦合数值模型的有效性和准确性；最后对系泊浮式结构大振幅运动响应进行预测，展示本发明的全非线性耦合方法在强非线性海洋环境下预测系泊浮式结构运动响应的能力。本算例仅为本发明的实施示例，不表示本发明的最终实施方案。

计算结果

响应振幅算子(RAO)相比于运动响应能更好的反应系泊浮体与波浪相互作用的运动特性。RAO通过将响应振幅(从最大值到最小值)与入射波高归一化，从时程曲线计算得到。图3为纵摇RAO随入射波周期的变化，从图中可以看出，随着波周期的增加，纵荡RAO几乎呈线性增大，模拟结果与模型实验结果非常接近，相对差异小于5％，显示本发明中提出的耦合分析方法求解精度相比于传统方法(传统方法一般可接受误差为15％以内)精度高。

图4纵摇运动的共振能量明显大于分布在波(主)频率附近的能量，说明显著的响应幅值出现在共振频率处，不可低估。为了进一步研究图4中纵摇运动的共振频率和频谱峰值，还绘制了入射波频率和相应的固有频率。图4中预测和测量的纵摇运动共振频率均接近0.09Hz，而浪涌的固有频率接近0.139Hz。这种频率不匹配的原因可能是系泊线的非线性效应。对于采用悬链线系泊的浮体，系泊系统的刚度对其纵荡响应起着关键作用。当浮式结构进行大幅度的水平运动时，系泊线可能被拉伸到一定程度，使系泊刚度呈现非线性。

可以看出，本发明所提的系泊浮式结构的全耦合-系泊非线性分析方法的计算结果要优于传统的耦合分析方法，且具有更高的精度。同时，所提耦合方法不仅能准确地捕捉到系泊浮体结构在纵荡、垂荡和纵摇运动下的波频响应，而且还能捕捉到系泊浮体结构在三个自由度下的共振响应，解决了传统动网格CFD方法对于系泊浮式结构大幅运动计算发散的问题。同时，本发明开发的非线性有限元系泊求解模块能很好的捕捉到系泊的非线性效应及其对浮式结构物纵摇运动的影响，显著提高数值模拟的精度。以上计算结果充分显示了本发明中提出的系泊浮式结构全耦合分析方法的有效性，以及对于恶劣海况下系泊浮式结构的生存能力研究的工程意义。

以上所述的本发明实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神和原则之内所作的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。