一种用于浮式海洋结构物静水平衡状态计算的有限元模型

摘要

一种用于浮式海洋结构物静水平衡状态计算的有限元模型，属于海洋工程结构物设计方法领域。其特点在于：将浮式海洋结构物有限元模型中湿表面单元的浮力作为单元的固有属性，建立单元浮力修正矩阵，通过浮力修正矩阵修正结构刚度矩阵。模型的垂向位置由重力与浮力平衡关系确定，避免了施加垂向自由度边界约束。通过求解总体平衡方程，计算各节点的广义位移向量，得到计入浮态与变形耦合作用的静水平衡状态。本发明适用于结构相对刚度较低的浮式海洋结构物，如下水工作船、超大型浮体、及玻璃钢等复合材料船舶的浮态与结构变形耦合作用下的平衡状态分析，可用于精确计算上述浮式海洋结构物的浮态、总纵强度、结构变形和结构应力。

说明书

技术领域

本发明涉及一种用于浮式海洋结构物静水平衡状态计算的有限元模型，其属于海洋工程结构物设计方法领域。

背景技术

船舶、浮船坞、浮箱和超大型浮体等浮式海洋结构物的船体结构具有一定的柔性，在水中漂浮时受不均衡的浮力与重力作用，将发生一定的弹性变形。结构弹性变形量与浮体的浮态相关，同时弹性变形的发生将进一步引起浮体所受静水浮力分布的变化，引起浮态的进一步变化。浮式海洋结构物的浮态与结构弹性变形相互影响达到平衡，所以其全浮状态下静水平衡状态，本质上是一种流固耦合现象。

现有的浮式海洋结构物浮态计算方法中，忽略浮态与结构变形的耦合问题，即在计算浮态时，将船体(广义船体，包括船舶和各类浮式海洋结构物)视为刚体，不考虑弹性变形引起的浮力变化。现有的船体结构变形分析方法中，将浮力作为固定载荷施加，不考虑结构变形进一步引起的浮力载荷的变化。对于常规运输船舶，船体结构刚度较大，结构变形对浮体分布的影响及浮力变化引起的船体结构变形均相对较小，所以采用上述简化原则加以简化，并对设计结果留有适当的安全余量，仍能够满足安全性、经济性等要求。但对于超大型浮体、浮箱、下水工作船、浮船坞等型长型深比相对较大的特殊浮式海洋结构物，以及采用玻璃钢等高强度、低刚度复合材料的船舶，船体的相对刚度显著低于常规船舶。同样载荷作用下，这类浮式海洋结构物的船体变形量显著大于钢制运输船舶，此时如果忽略船体结构变形与浮态的耦合作用，对设计中的浮态计算、总纵强度计算、结构局部强度分析等计算精度均会造成较大的影响。

发明内容

现有的浮式海洋结构物设计方法不能准确求解计入结构变形与浮态耦合作用的静水平衡状态。针对该问题，本发明提出一种用于浮式海洋结构物静水平衡状态计算的有限元模型。该有限元模型将浮式海洋结构物有限元模型中湿表面单元的浮力作为单元的固有属性，建立单元的浮力修正矩阵，通过浮力修正矩阵修正单元的结构刚度矩阵。单元的结构刚度矩阵经过浮力修正矩阵修正后，单元受到的浮力载荷与其变形量满足流体压力原理，可准确模拟浮式海洋结构物在静水中漂浮时，结构湿表面所受压力随浮态与变形的变化而改变的问题。有限元模型的垂向位置由重力与浮力平衡关系确定，不需要施加垂向自由度边界约束。通过求解总体平衡方程，得到各节点的广义位移向量，其中包括浮态变化引起的刚性位移及结构弹性变形，实现浮态与变形耦合求解，得到浮式海洋结构物自由漂浮时的静水平衡状态。

一种用于浮式海洋结构物静水平衡状态计算的有限元模型，其发明内容如下。所述的用于浮式海洋结构物静水平衡状态计算的有限元模型，基本原理如图1所示，其单元刚度矩阵、浮力单元的浮力修正矩阵、总体刚度矩阵、载荷及边界条件、模型的求解方法如下。

1)单元刚度矩阵。

对于浮式海洋结构物湿表面的壳单元或者膜单元(简称浮力单元)，采用浮力修正矩阵修正单元的结构刚度矩阵

$K^e=K_S^e+K_B^e$

其中，K^e为经浮力修正矩阵修正后的单元刚度矩阵，为单元的结构刚度矩阵，为浮力修正矩阵。对于模型中的非浮力单元，不对单元的结构刚度矩阵做修正。

2)浮力单元的浮力修正矩阵。

如图2所示，对于三角形、四边形浮力单元，C为单元的坐标中心点。将C与各边的中点M_i(1≤i≤N_E)连接得到N_E条直线段，N_E为单元的边数。N_E条直线段将单元的几何面划分为N_E个面片S_i。通过曲面积分算法，计算S_i的面积A_i和法向

对于4节点四边形壳单元、4节点四边形膜单元、8节点四边形壳单元及8节点四边形膜单元，所述浮力修正矩阵为

其中，ρ为水密度，g为重力加速度。v_x＝(1,0,0)，v_y＝(0,1,0)，v_z＝(0,0,1)分别为X轴、Y轴和Z轴的矢量方向，其中X方向、Y方向和Z方向分别为浮式海洋结构物长度方向、宽度方向及型深方向。

对于3节点三角形壳单元、3节点三角形膜单元、6节点三角形壳单元及6节点三角形膜单元，所述浮力修正矩阵为

所述的浮力修正矩阵，对于4节点四边形壳单元、3节点三角形壳单元，R_f取值为2，R_m取值为5，R_e取值为3，C_m取值为3，C_e取值为3，对于4节点四边形膜单元、3节点三角形膜单元，R_f取值为2，R_m取值为2，R_e取值为0，C_m取值为0，C_e取值为0。

所述的浮力修正矩阵，对于8节点四边形壳单元、6节点三角形壳单元，R_f取值为2，R_m取值为5，R_e取值为6N_E+3，C_m取值为3，C_e取值为6N_E+3；对于8节点四边形膜单元、6节点三角形膜单元，R_f取值为2，R_m取值为2，R_e取值为3N_E，C_m取值为0，C_e取值为3N_E，其中，N_E为单元的边数。

3)总体刚度矩阵。

所述用于浮式海洋结构物静水平衡状态计算的有限元模型，采用常规有限元法，将单元刚度矩阵组装，得到总体刚度矩阵。

4)载荷与边界约束。

所述用于浮式海洋结构物静水平衡状态计算的有限元模型，将浮力载荷作为单元的固有属性，体现于单元刚度矩阵中，所以载荷向量中不包含浮力分量。

所述用于浮式海洋结构物静水平衡状态计算的有限元模型，模型的Z方向位置由隐含于总体平衡方程中的重力浮力平衡关系确定，所以边界约束中不包括Z方向线位移约束以及绕X轴、绕Y轴的转动约束。

5)总体平衡方程求解。

所述用于浮式海洋结构物静水平衡状态计算的有限元模型，通过求解总体平衡方程，计算各节点的广义位移向量，完成浮体平衡状态求解。广义位移向量包括两部分，一是因浮式海洋结构物浮态变化引起的位移向量，二是浮式海洋结构物结构弹性变形引起的位移向量。

附图说明

图1为一种用于浮式海洋结构物静水平衡状态计算的有限元模型基本原理图。

图2为单元曲面片分割算法示意图，其中，图2(a)为三角形单元示意图，图2(b)为四边形单元示意图。

图3为实施例中浮体示意图。

图4为实施例中有限元网格划分示意图。

图5为实施例中单元划分及曲面片编号示意图。

具体实施方式

以实例说明所述一种用于浮式海洋结构物静水平衡状态计算的有限元模型的具体实施方式。

如图3所示组合浮体由一块钢板与一块泡沫构成。钢板尺寸为2000mm×1000mm，厚度为10mm，其上方粘连的硬泡沫尺寸为2000mm×1000mm，厚度400mm。将该组合浮体置于淡水中，钢板中间位置处，左右两边各施加1000N集中力。钢的密度为7860kg/m³,杨氏模量为2.06×10⁶Mpa,泊松比为0.3。水密度为1000kg/m³。重力加速度取9.8。忽略泡沫的质量，忽略泡沫因水压力产生的变形，求组合体在静水中的平衡状态。

采用所述用于浮式海洋结构物静水平衡状态计算的有限元模型，计算上述浮体平衡状态具体实施方式如下。

该算例中，泡沫仅起提供浮力作用，所以只建立钢板的有限元模型，把泡沫提供的浮力作为钢板单元的属性，将钢板单元定义为浮力单元。

1)划分网格，建立单元结构刚度矩阵。

如图4所示，将板划分为两个边长为1000mm的正方形4节点四边形壳单元e₁和e₂。e₁依次由节点1,2,3,4构成，e₂依次由节点2,5,6,3构成。

对于壳单元e₁和e₂，结构刚度矩阵为24×24对称矩阵，且两个单元的结构刚度矩阵相同。采用常规有限元法，计算单元结构刚度矩阵

$K_S^e=\left(\begin{array}{cccc}{k}_{11}& k_{12}& k_{13}& k_{14}\\ k_{21}& k_{22}& k_{23}& k_{24}\\ k_{31}& k_{32}& k_{33}& k_{34}\\ k_{41}& k_{42}& k_{43}& k_{44}\end{array}\right)$

其中，

$k_{12}=\left(\begin{array}{cccccc}-539029& -27802& 0& 0& 0& -495192\\ 27802& 195695& 0& 0& 0& 495192\\ 0& 0& -95& 8430& -42622& 0\\ 0& 0& 8430& 3124428& 766369& 0\\ 0& 0& 42622& -766369& 9491186& 0\\ -495192& -495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{13}=\left(\begin{array}{cccccc}-592840& -367362& 0& 0& 0& -495192\\ -367362& -592840& 0& 0& 0& 495192\\ 0& 0& -9& 9373& -9373& 0\\ 0& 0& -9373& 7074176& 1591690& 0\\ 0& 0& 9373& 1591690& 7074176& 0\\ 495192& -495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{14}=\left(\begin{array}{cccccc}195695& 27802& 0& 0& 0& -495192\\ -27802& -539029& 0& 0& 0& 495192\\ 0& 0& -95& 42622& -8430& 0\\ 0& 0& -42622& 9491186& -766369& 0\\ 0& 0& -8430& 766369& 3124428& 0\\ 495192& 495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{21}=\left(\begin{array}{cccccc}-539029& 27802& 0& 0& 0& -495192\\ -27802& 195695& 0& 0& 0& -495192\\ 0& 0& -95& 8430& 42622& 0\\ 0& 0& 8430& 3124428& -766369& 0\\ 0& 0& -42622& 766369& 9491186& 0\\ -495192& 495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{23}=\left(\begin{array}{cccccc}195695& -27802& 0& 0& 0& -495192\\ 27802& -539029& 0& 0& 0& -495192\\ 0& 0& -95& 42622& 8430& 0\\ 0& 0& -42622& 9491186& 766369& 0\\ 0& 0& 8430& -766369& 3124428& 0\\ 495192& -495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{24}=\left(\begin{array}{cccccc}-592840& 367362& 0& 0& 0& -495192\\ 367362& -592840& 0& 0& 0& -495192\\ 0& 0& -9& 9373& 9373& 0\\ 0& 0& -9373& 7074176& -1591690& 0\\ 0& 0& -9373& -1591690& 7074176& 0\\ 495192& 495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{31}=\left(\begin{array}{cccccc}-592840& -367362& 0& 0& 0& 495192\\ -367362& -592840& 0& 0& 0& -495192\\ 0& 0& -9& -9373& 9373& 0\\ 0& 0& 9373& 7074176& 1591690& 0\\ 0& 0& -9373& 1591690& 7074176& 0\\ -495192& 495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{32}=\left(\begin{array}{cccccc}195695& 27802& 0& 0& 0& 495192\\ -27802& -539029& 0& 0& 0& -495192\\ 0& 0& -95& -42622& 8430& 0\\ 0& 0& 42622& 9491186& -766369& 0\\ 0& 0& 8430& 766369& 3124428& 0\\ -495192& -495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{34}=\left(\begin{array}{cccccc}-539029& -27802& 0& 0& 0& 495192\\ 27802& 195695& 0& 0& 0& -495192\\ 0& 0& -95& -8430& 42622& 0\\ 0& 0& -8430& 3124428& 766369& 0\\ 0& 0& -42622& -766369& 9491186& 0\\ 495192& 495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{41}=\left(\begin{array}{cccccc}195695& -27802& 0& 0& 0& 495192\\ 27802& -539029& 0& 0& 0& 495192\\ 0& 0& -95& -42622& -8430& 0\\ 0& 0& 42622& 9491186& 766369& 0\\ 0& 0& -8430& -766369& 3124428& 0\\ -495192& 495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{42}=\left(\begin{array}{cccccc}-592840& 367362& 0& 0& 0& 495192\\ 367362& -592840& 0& 0& 0& 495192\\ 0& 0& -9& -9373& -9373& 0\\ 0& 0& 9373& 7074176& -1591690& 0\\ 0& 0& 9373& -1591690& 7074176& 0\\ -495192& -495192& 0& 0& 0& 495192308\end{array}\right)$

$k_{43}=\left(\begin{array}{cccccc}-539029& 27802& 0& 0& 0& 495192\\ -27802& 195695& 0& 0& 0& 495192\\ 0& 0& -95& -8430& -42622& 0\\ 0& 0& -8430& 3124428& -766369& 0\\ 0& 0& 42622& 766369& 9491186& 0\\ 495192& -495192& 0& 0& 0& 495192308\end{array}\right)$

2)建立单元的浮力修正矩阵。

如图5所示，采用本发明提出的单元曲面片分割算法，取4个边的中点M₁,M₂,M₃和M₄，连接单元的坐标中点C与四个中点，得到四条直线段。四条直线段将单元划分为S₁,S₂,S₃和S₄四个曲面片。四个曲面片的面积均为250000mm²，法向均为(0,0,1)。建立单元的浮力修正矩阵

$K_B^e=\left(\begin{array}{cccc}{b}_{11}& b_{12}& b_{13}& b_{14}\\ b_{21}& b_{22}& b_{23}& b_{24}\\ b_{31}& b_{32}& b_{33}& b_{34}\\ b_{41}& b_{42}& b_{43}& b_{44}\end{array}\right)$

其中，

$b_{11}=\left(\begin{array}{cccccc}0& 0& 1.378& 0& 0& 0\\ 0& 0& 1.378& 0& 0& 0\\ 0& 0& 1.378& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),b_{12}=\left(\begin{array}{cccccc}0& 0& 0.459& 0& 0& 0\\ 0& 0& 0.459& 0& 0& 0\\ 0& 0& 0.459& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),b_{13}=\left(\begin{array}{cccccc}0& 0& 0.153& 0& 0& 0\\ 0& 0& 0.153& 0& 0& 0\\ 0& 0& 0.153& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right)$

b₁₄＝b₂₁＝b₂₃＝b₃₂＝b₃₄＝b₄₁＝b₄₃＝b₁₂

b₂₄＝b₃₁＝b₄₂＝b₁₃

b₂₂＝b₃₃＝b₄₄＝b₁₁

3)建立单元刚度矩阵。

计入浮力修正的单元刚度矩阵为单元的结构刚度矩阵和浮力修正矩阵之和，即浮力修正后的单元刚度矩阵K^e为

$K^e=\left(\begin{array}{cccc}{K}_{11}& K_{12}& K_{13}& K_{14}\\ K_{21}& K_{22}& K_{23}& K_{24}\\ K_{31}& K_{32}& K_{33}& K_{34}\\ K_{41}& K_{42}& K_{43}& K_{44}\end{array}\right)$

其中，K_ij＝k_ij+b_ij(1≤i≤4,1≤j≤4)

4)建立总体刚度矩阵。

根据模型中两个单元的拓扑关系，将两个单元的单元刚度矩阵组装得到总体刚度矩阵。总体刚度矩阵为36×36对称矩阵，表达为

$K^G=\left(\begin{array}{cccccc}{K}_{11}& K_{12}& K_{13}& K_{14}& 0& 0\\ K_{21}& K_{22}+K_{11}& K_{23}+K_{14}& K_{24}& K_{12}& K_{13}\\ K_{31}& K_{32}+K_{41}& K_{33}+K_{44}& K_{34}& K_{42}& K_{43}\\ K_{41}& K_{42}& K_{43}& K_{44}& 0& 0\\ 0& K_{21}& K_{24}& 0& K_{22}& K_{23}\\ 0& K_{31}& K_{34}& 0& K_{32}& K_{33}\end{array}\right)$

5)建立总体载荷向量。

所述用于浮式海洋结构物静水平衡状态计算的有限元模型，将浮力作为结构的固有属性，所以不将单元的浮力作为载荷考虑，等效节点载荷向量中没有浮力分量。

各节点承担一部分钢板的重力载荷，节点2与节点3受两个1000N的集中力作用。总体载荷向量为

F^G＝{0,0,-193,0,0,0,0,0,-1385,0,0,0,0,0,-1385,0,0,0,0,0,-193,0,0,0,0,0,-193,0,0,0,0,0,-193,0,0,0}^T

6)施加边界约束。

对节点1和节点4，固定X方向线位移和Y方向线位移自由度。对节点5和节点6，固定Y方向线位移自由度。其余自由度均自由。所述用于浮式海洋结构物静水平衡状态计算的有限元模型，模型的垂向位置由重力浮力平衡条件确定，所以不约束Z方向的线位移、绕X轴的转角位移和绕Y轴的转角位移。

采用置大数法处理边界约束。将上述6个固定约束对应的自由度在总体刚度矩阵中相应的对角线元素设置为10²⁰。

7)总体平衡方程求解。

施加约束后，消除总体刚度矩阵的奇异性，总体平衡方程有解。求解总体平衡方程，得到各节点广义位移向量

δ^G＝{δ₁,δ₂,δ₃,δ₄,δ₅,δ₆}^T

其中

δ₁＝{0,0,-175.91,0,0.01,0}^Tδ₂＝{0,0,-185.37,0,0,0}^T

δ₃＝{0,0,-185.37,0,0,0}^Tδ₄＝{0,0,-175.91,0,0.01,0}^T

δ₅＝{0,0,-175.91,0,-0.01,0}^Tδ₆＝{0,0,-175.91,0,-0.01,0}^T

将所述用于浮式海洋结构物静水平衡状态计算的有限元模型应用于浮式海洋结构物的静水平衡状态分析时，将浮式海洋结构物的湿表面壳单元或者膜单元定义为浮力单元并对其单元刚度矩阵进行修正，其余过程与本例相同。