一种波浪作用过程中深水网箱浮架变形的数值模拟方法与流程

本发明涉及一种数值模拟方法，尤其涉及一种用于模拟分析波浪作用过程中深水网箱浮架变形的数值模拟方法。

背景技术：

深水网箱系统主要由浮架、网衣及锚泊系统三大部分组成，其中浮架系统负责提供足够的浮力保持网箱处于漂浮状态，为网箱日常管理提供操作平台，并与锚泊系统相连共同抵抗波浪的冲击，在网箱抗风浪性能中扮演着极为重要的角色。随着深水网箱大型化发展及深远海养殖拓展，占我国深水网箱总数90％以上的HDPE网箱在恶劣海况下受到的海浪冲击作用更大，致使网箱浮架系统有可能无法承受过大的锚绳力而发生结构变形破坏，给养殖户或养殖企业带来巨大损失。如2013年超强台风“尤特”和2014年超强台风“威马逊”造成广东100多个深水网箱破坏，其中多数网箱的浮架因变形严重而发生断裂。

获取波浪作用下网箱浮架受力变形数据是开展深水网箱科学设计和制作安装的关键，水动力模型试验、计算机数值模拟是开展网箱力学研究的主要手段。一方面，在实验室开展网箱水动力模型试验时，因制作的小比例尺网箱模型很难满足浮架刚度相似条件，在实验室造波条件下浮架基本是未发生变形的，只能获得浮架部分受力及运动数据，获取不到重要的浮架变形数据，这与实际中浮架是柔性弹性体的情况不符；另一方面，国内外开展的网箱数值模拟研究因计算理论的缺失和计算条件的限制，绝大多数都假设浮架是刚体，未考虑浮架的变形。而在有限的浮架变形模拟研究中，仅考虑集中荷载情况，不能模拟波浪作用下浮架的实时变形情况。

技术实现要素：

本发明所要解决的技术问题是：提供一种波浪作用过程中深水网箱浮架变形的数值模拟方法。

解决上述技术问题，本发明所采用的技术方案如下：

一种波浪作用过程中深水网箱浮架变形的数值模拟方法，其特征在于：所述的数值模拟方法包括：

步骤S1、建立模型，包括：

步骤S1.1、用有限元法分别将受模拟深水网箱浮架与水体直接接触的浮管以及与所述浮管相连的锚绳分解成多个线微元和节点，建立每个所述节点的运动方程；

步骤S1.2、建立坐标系，获取步骤S1.1分解得到的每个所述节点在初始时刻的速度向量v₁和对应所述坐标系的位置向量p₁；

步骤S2、设置n个测试时刻，将所述初始时刻设定为第一个测试时刻，并设定相邻的第t个测试时刻t与第t+1个测试时刻之间的时间步长为dT，其中，n为正整数，t为取值在1到n之间的正整数；从第一个测试时刻即t取值为1开始，依次循环执行以下步骤S2.1至步骤S2.5，直至t的取值达到预设的最大值或者本数值模拟方法被停止：

步骤S2.1、依据每个所述线微元在第t个测试时刻的速度向量v_t和位置向量p_t，相应计算每个所述线微元在第t个测试时刻所受到的结构荷载和非结构荷载，其中，对于由所述浮管分解得到的线微元，其结构荷载包括张力、弯矩、剪力和扭矩，对于由所述锚绳分解得到的线微元，其结构荷载包括张力，每个所述线微元的非结构荷载均包括水动力、重力和浮力；

步骤S2.2、将每个所述线微元在第t个测试时刻所受到的结构荷载和非结构荷载平均分配到与该线微元相邻的两个所述节点上，并对每个所述节点所分配到的力和力矩分别进行合成，得到每个所述节点在第t个测试时刻的总合力和总力矩；

步骤S2.3、将每个所述节点在第t个测试时刻的总合力和总力矩代入对应节点的运动方程中，并采用半隐式欧拉积分方法求解出每个所述节点在第t个测试时刻的加速度a_t；

步骤S2.4、按以下公式计算出每个所述节点在第t+1个测试时刻的速度向量v_t+1和位置向量p_t+1：

v_t+1＝v_t+a_tdT

p_t+1＝p_t+v_t+1dT；

步骤S2.5、将t的取值加1；

步骤S3、对于每个所述测试时刻，用所述各个节点在该测试时刻所受到的结构荷载，计算出所述各个节点在该测试时刻的von Mises应变和von Mises应力，用以反映受模拟深水网箱浮架对应所述各个节点的位置在所述n个测试时刻的变形情况。

与现有技术相比，本发明具有以下有益效果：

第一，本发明作为一种深水网箱浮架在波浪作用过程中变形的模拟方法，其采用有限元法建立浮架模型，通过求解结构载荷来获得浮架各位置的应力和应变值，以综合反映浮架的变形情况；且可根据实际情况或模拟需求设定有限元法所分解线微元的数量，其计算效率高，计算时间基本在1小时以内，解决了现有浮架有限元模拟中存在的网格数量多、成本高、耗时长的问题。

第二，本发明和网箱水动力模型试验相比，计算精度较高，可以替代物理模型试验，解决现有模型试验过程中浮架不能变形的问题，能够模拟任意波浪条件下原尺寸深水网箱的浮架实时变形情况。

第三，本发明可用于评估浮架极限承载能力，将特定波况下浮架变形的应力值与HDPE管材的屈服强度进行比较，来判断浮架是否失效。

综上所述，本发明能够综合反映浮架的变形情况，具有计算效率高、精度高、用途广泛的优点。

附图说明

下面结合附图和具体实施例对本发明作进一步的详细说明：

图1为本发明的数值模拟方法的流程框图；

图2为简化后的网箱浮架结构模型示意图；

图3为有限元模型计算示意图；

图4为浮管结构张力计算示意图；

图5为规则波条件下不同时刻的浮架整体变形外观图；

图6为规则波和不规则波条件下浮架最大变形位置的应变历时曲线；

图7为各种不同波况下浮架最大变形位置处的应变峰值。

具体实施方式

如图1所示，本发明的波浪作用过程中深水网箱浮架变形的数值模拟方法，包括：

步骤S1、建立模型，包括：

步骤S1.1、用有限元法分别将受模拟深水网箱浮架与水体直接接触的浮管以及与浮管相连的锚绳分解成多个线微元和节点，建立每个节点的运动方程；

其中，上述节点的运动方程形式如下：

M(p，a)+C(p，v)+K(p)＝F(p，v，t) (1)

式中：M(p,a)是惯性载荷，C(p,v)是阻尼载荷，K(p)是刚度载荷，F(p,v,t)是外部载荷，p,v和a分别是代表位置向量、速度向量和加速度向量，t是计算时间。

参见图3，由有限元法分解得到的每个线微元均与两个节点相邻且相邻的两个线微元之间存在一个节点，如线微元segmengt1的一侧相邻的节点node1和另一侧相邻的节点node2，节点node2为存在于线微元segmengt1与线微元segmengt2之间的节点。其中，由于深水网箱浮架系统漂浮于水面上，通常由浮管、扶手、立管、三通及其它连接件组成，而与水体直接接触的浮管及浮管连接件是浮架承受波浪冲击的主要构件，为了方便计算浮架受力，因此，以获取波浪作用过程中深水网箱浮架变形数据为目的时，可以将受模拟深水网箱浮架简化成浮管及浮管连接件；上述步骤S1.1即是基于该思路设计的。

步骤S1.2、建立坐标系，获取步骤S1.1分解得到的每个节点在初始时刻的速度向量v₁和对应坐标系的位置向量p₁；

步骤S2、设置n个测试时刻，将初始时刻设定为第一个测试时刻，并设定相邻的第t个测试时刻t与第t+1个测试时刻之间的时间步长为dT，其中，n为正整数，t为取值在1到n之间的正整数；从第一个测试时刻即t取值为1开始，依次循环执行以下步骤S2.1至步骤S2.5，直至t的取值达到预设的最大值或者本数值模拟方法被停止：

步骤S2.1、依据每个线微元在第t个测试时刻的速度向量v_t和位置向量p_t，相应计算每个线微元在第t个测试时刻所受到的结构荷载和非结构荷载，其中，对于由浮管分解得到的线微元，其结构荷载包括张力、弯矩、剪力和扭矩，对于由锚绳分解得到的线微元，其结构荷载包括张力，每个线微元的非结构荷载均包括水动力、重力和浮力；

其中，上述张力、弯矩、剪力、扭矩和水动力，重力和浮力为常规计算方式，在此不赘述：

第一，张力计算：

浮架和锚绳各线微元的有效张力表达式可以统一写成如下形式：

T_e＝T_w+(P₀*A₀-P_i*A_i) (2)

式中，T_w代表环壁张力，可表示为如下形式：

T_w＝EA*ε-2μ(P₀*A₀-P_i*A_i)+EA*C₁(dL/dt)/L₀ (3)

式中，EA是线微元的轴向刚度，ε是轴向应变，L₀和L分别是线微元在初始时刻的长度和在第t个测试时刻的长度，μ是泊松比，P_i和P₀是内、外压力，A_i和A₀是线微元的内、外横截面积，C₁是阻尼系数。对于锚绳，由于微元的内、外压力及横截面积基本相等，故T_e等于T_w。

第二，弯矩计算：

浮架线微元的弯矩计算表达式为：

M＝EIk+kd(λ_b/100)D_c/dt (4)

式中，EI为弯曲刚度，k为曲率，λ_b为目标弯曲阻尼，D_c是弯曲临界阻尼值，可以表示为L₀(mEIL₀)^1/2。其中，将位于两个相邻线微元之间的节点记为中间节点，与该中间节点相邻的两个线微元的弯矩分别为M₁和M₂。

第三，剪力计算：

由于线微元在弯曲方面是刚性的，弯矩沿线微元长度是线性变化的，如下所示：

Q＝S_Z×(M₂-M₁)/L (5)

式中，S_Z是轴向单位向量，L是线微元在第t个测试时刻的长度。

第四，扭矩计算：

浮架线微元的扭矩由下式计算：

T＝Kτ/L₀+C₂(dτ/dt) (6)

式中，K为扭转刚度，τ为扭转角度，L₀为线微元在初始时刻的长度，C₂是扭转阻尼系数。

第五，水动力载荷计算：

采用扩展的莫里森方程计算浮架和锚绳的每个线微元的波浪力。波浪力由两部分组成，一部分为与水质点速度相关的速度力，另一部分为与水质量加速度相关的惯性力。采用前述莫里森方程计算的波浪力表达式为：

$F_W=\frac{1}{2}{\mathrm{\ρC}}_d{\mathrm{SV}}_r\left|V_r\right|+\mathrm{\ρ}\∀a_f+\mathrm{\ρ}\∀C_a{a}_r---\left(7\right)$

式中，ρ为海水密度，S为阻力面积，V_r为水质点与构件的相对速度，为构件体积，a_f为水质点加速度，a_r为水质点与线微元的相对加速度，C_d和C_a分别为速度力系数和附加质量力系数。

步骤S2.2、将每个线微元在第t个测试时刻所受到的结构荷载和非结构荷载平均分配到与该线微元相邻的两个节点上，并对每个节点所分配到的力和力矩分别进行合成，得到每个节点在第t个测试时刻的总合力和总力矩；

其中，在波浪作用下，浮架的结构载荷包括张力、弯矩、剪力及扭矩，而对于与浮架相连的锚绳，其结构载荷则只有张力。用有限元法分解得到的每个线微元可视为无质量的弹性构件，用于模拟每个微元的结构载荷，对于锚绳而言，由于锚绳是柔软的可任意弯曲，故在计算锚绳力学性能时可不考虑锚绳的弯矩和扭矩特性；上述步骤S2.1和步骤S2.2即是基于该思路设计的。

步骤S2.3、将每个节点在第t个测试时刻的总合力和总力矩代入对应节点的运动方程即上述公式(1)中，并采用半隐式欧拉积分方法求解出每个节点在第t个测试时刻的加速度a_t；

步骤S2.4、按以下公式计算出每个节点在第t+1个测试时刻的速度向量v_t+和位置向量p_t+1：

v_t+1＝v_t+a_tdT (8)

p_t+1＝p_t+v_t+1dT (9)；

步骤S2.5、将t的取值加1；

步骤S3、对于每个测试时刻，用各个节点在该测试时刻所受到的结构荷载，计算出各个节点在该测试时刻的von Mises应变和von Mises应力，用以反映受模拟深水网箱浮架对应各个节点的位置在n个测试时刻的变形情况。

其中，上述von Mises应变和von Mises应力的计算表达式如下：

${\mathrm{\ϵ}}_{vm}=\sqrt{{\mathrm{\ϵ}}_{zz}^2+{\mathrm{\ϵ}}_{cc}^2-{\mathrm{\ϵ}}_{zz}{\mathrm{\ϵ}}_{cc}}---\left(10\right)$

${\mathrm{\σ}}_{vm}=\sqrt{\frac{{({\mathrm{\σ}}_1-{\mathrm{\σ}}_2)}^2+{({\mathrm{\σ}}_2-{\mathrm{\σ}}_3)}^2+{({\mathrm{\σ}}_3-{\mathrm{\σ}}_1)}^2}{2}}---\left(11\right)$

式中，ε_zz是轴向应变，ε_cc是环向应变，σ₁,σ₂,σ₃分别是三个方向的主应力。

下面结合一个应用实例，说明本发明的数值模拟方法的实际应用效果：

针对图2所示的网箱浮架结构(图中示出了深水网箱浮架与水体直接接触的浮管以及该浮管与A’点、B’点、C’点、D’点、E’点和F’点之间连接的锚绳，以及步骤S1.2所建立坐标系OXYZ)，采用有限元法将浮架和锚绳分解成多个线微元及节点，如图3所示(图中举例示出了线微元segmengt1、线微元segmengt2、线微元segmengt3、节点node1、节点node2和节点node3的位置关系，EndA为浮架或锚绳的一端终点，另一端终点EndB的情况未示出)。计算的网箱浮架结构详细参数见表1所示。

表1浮架结构详细参数表

分别设定波浪条件为规则波和不规则波情况，见表2所示。为了比较规则波和不规则波条件下的浮架变形情况，取不规则波的有效波高和平均周期分别为规则波的波高和周期。其中，不规则波浪采用JONSWAP谱模拟，其谱密度函数表达式如下：

$S\left(f\right)=\frac{{\mathrm{\βg}}^2}{16{\mathrm{\π}}^4{f}^5}\exp \[-\frac{5}{4}{\left(\frac{f_p}{f}\right)}^4\]{\mathrm{\γ}}^{\exp \[\frac{-{(f/f_p-1)}^2}{2{\mathrm{\σ}}^2}\]}---\left(12\right)$

其中，f是波浪频率，f_p是谱峰频率，σ是谱形状参数，当f≤f_p时，σ＝0.07；当f≥f_p时，σ＝0.09。γ是峰高因子，正常设为3.3。

表2波浪计算条件

分别计算浮架和锚绳受到的所有载荷如张力、弯矩、剪力、扭矩以及水动力、重力和浮力等，其中浮架和锚绳受到的张力可统一参考图4(图中示出了相邻的节点node n和节点node n+1之间的线微元segmengt mid-point)进行计算。在获得每个节点的合力及力矩后，通过求解运动方程获取节点的加速度、速度及相关位置，从而进一步求解获得浮架各位置的变形情况。

图5为规则波波高6m、周期8s条件下浮架在一个波浪周期内不同时刻的变形形状图。可以看出浮架在波浪作用下发生明显的变形，且浮架在t＝T/2和t＝T时刻相比其它时刻具有较大的变形。

图6为波高6m、周期8s时对应规则波和不规则波条件下的浮架最大变形位置的应变历时曲线图。可以看出，不规则波下浮架最大应变值达到了1.23％，是规则波条件下最大应变值0.49％的2.5倍。

图7为规则波、不规则波分别在不同波高和周期条件下的浮架最大变形比较图。表明最大应变值随着波高的增加而增加。随着波浪周期的增加，规则波条件下的浮架最大应变值变化很小，而不规则波条件下的浮架最大应变值不断递减。对应于各种波况，不规则波条件下的浮架最大应变均要大于规则波情况，且两者之间的差异随着波高的增加而增加，随着波浪周期的增加而减小。

本发明不局限于上述具体实施方式，根据上述内容，按照本领域的普通技术知识和惯用手段，在不脱离本发明上述基本技术思想前提下，本发明还可以做出其它多种形式的等效修改、替换或变更，均落在本发明的保护范围之中。