波浪入射角及地形对浮体水动力学特性影响的有限元分析

摘要

利用有限元方法分析了二维情况下波浪入射角及地形对浮体水动力学特性的影响，计算了不同情况下浮体的水动力学系数；通过与简单地形情况下解析解的比较，验证了数值方法的正确性；并在分析中采用了BiCGStab算法，高效地求解离散形成的线性方程组，大大地提高了分析的效率。研究结果为海上结构物的设计提供了关键的水动力学数据。

正文

《1 前言》

1 前言

随着人类对海洋开发的日益重视, 波浪和浮体相互作用的研究已引起人们越来越多的关注。这方面的研究不仅有重要的理论意义, 还有广泛的工程应用价值。如, 海洋结构物的设计必须考虑波浪引起的动载荷;船舶设计必须考虑兴波阻力和耐波性;海洋波浪能利用必须考虑波浪与波能装置的相互作用等。

关于波浪—浮体相互作用的研究方法大致上可分为2类:解析方法^[1,2,3,4,5]和数值方法^{[6,7,8,9,10,11,12,13]}。解析方法虽然精度和效率一般都要优于数值方法, 但它只适用一些简单的情形。对复杂问题 (如非线性、复杂地形等) 基本上无法使用解析方法, 只能采用数值方法。

目前国内外对波浪—浮体系统的研究大多限于正向波和简单地形情况, 而对实际的斜向波和复杂地形联合作用下浮体水动力学特性的研究^[10,12,13]还很少。随着海上运输、海洋资源开发向着更具挑战性的海域推进, 波浪入射角及地形对浮体水动力学特性的影响成为急需解决的重要课题。作者采用基于三角形网格的有限元方法进行分析, 这种方法能很好地适应地形。在验证了数值方法的正确性后, 计算了不同情况下浮体的水动力学系数, 分析了波浪入射角和地形对其水动力学特性的影响。这些结果对设计更加安全可靠、经济合理的海洋结构物有重要参考价值。在设计分析所用的有限元程序时, 形成一个有效的任意平面区域三角形网格的全自动生成器, 有效地扩展了分析的范围。在求解离散形成的线性方程组这一分析中最耗时的部分时, 采用了目前被认为是最快算法之一的BiCGStab算法^[14,15], 大大地提高分析的效率。

《2 数学模型》

2 数学模型

研究复杂海底地形上, 一截面形状任意的长柱体在斜向入射波作用下的水动力学特性。入射波与x轴的夹角为θ, 问题的几何描述和坐标布局如图1所示。

《图1》

图1斜向波与浮体相互作用示意图

Fig.1 Interaction of oblique wave and buoy

在分析波浪入射角θ和地形对浮体水动力学特性的影响时, 为使问题简化, 建立数学模型时做了如下假设:

1) 流体为无粘、无旋的理想流体, 因此可以应用势流理论。控制方程为:

$Δ^{2} Φ (x, y, z; t) = 0 (1)$ $Δ^{2} Φ (x, y, z; t) = 0 (1)$

2) 波倾系数ε=A/λ为小量, 因此能用线性理论处理。

3) 现象是周期性的, 因此各物理量的时间项可按 (2) 式分离 $(i = \sqrt{- 1})$ $(i = \sqrt{- 1})$ 。

$F (x, y, z; t) = f (x, y, z) e^{- i ω t} (2)$ $F (x, y, z; t) = f (x, y, z) e^{- i ω t} (2)$

4) 浮体在y方向为无限长, 因此速度势在该方向可假定为周期性变化, 这样就可按以下方法处理成二维问题。

在上述假定下, 图1所示的流场速度势可表示为

$Φ (x, y, z) = ϕ (x, z) e^{i k y s i n θ} (3)$ $Φ (x, y, z) = ϕ (x, z) e^{i k y s i n θ} (3)$

由于图1中等深线为平行直线, 根据Snell定律

$k s i n θ = c o n s t (4)$ $k s i n θ = c o n s t (4)$

利用式 (1) ～ (4) , 原三维拉普拉斯方程变成下面的变形海姆霍兹方程^[10,12,13]

$Δ^{2} ϕ (x, z) - (k s i n θ)^{2} ϕ (x, z) = 0 (5)$ $Δ^{2} ϕ (x, z) - (k s i n θ)^{2} ϕ (x, z) = 0 (5)$

《2.1速度势的控制方程及边界条件》

2.1速度势的控制方程及边界条件

为了求解方便, 对于线性波浪问题, 通常利用Haskind方法, 将ϕ分解为

$ϕ = ϕ^{D} + \sum_{α = 1}^{6} V_{α} ϕ_{α} (6)$ $ϕ = ϕ^{D} + \sum_{α = 1}^{6} V_{α} ϕ_{α} (6)$

式中, ϕ^D表示入射波在浮体固定不动时的绕射势;ϕ_α表示浮体速度为1时强迫运动的广义模式, 即辐射势。α对应的运动模式分别为:1.横荡;2.纵荡;3.垂荡;4.横摇;5.纵摇;6.首摇。对二维问题, α为1, 3, 5。

经以上处理后, ϕ^D可由下列方程及边界条件求出

$\begin{array}{l} Δ^{2} ? ϕ^{D} - (k s i n θ)^{2} ϕ^{D} = 0 (在流体中) (7) \\ \frac{\partial ϕ^{D}}{\partial n} - \frac{ω^{2}}{g} ϕ^{D} = 0 (在自由表面上) (8) \\ \frac{\partial ϕ^{D}}{\partial n} = 0 (在水底) (9) \\ \frac{\partial ϕ^{D}}{\partial n} = 0 (在浮体表面) (10) \\ ϕ^{D} - ϕ^{Ι} (在无穷远处 ‚ 波以有限值向外传播) (11) \end{array}$ $\begin{array}{l} Δ^{2} ? ϕ^{D} - (k s i n θ)^{2} ϕ^{D} = 0 (在流体中) (7) \\ \frac{\partial ϕ^{D}}{\partial n} - \frac{ω^{2}}{g} ϕ^{D} = 0 (在自由表面上) (8) \\ \frac{\partial ϕ^{D}}{\partial n} = 0 (在水底) (9) \\ \frac{\partial ϕ^{D}}{\partial n} = 0 (在浮体表面) (10) \\ ϕ^{D} - ϕ^{Ι} (在无穷远处 ‚ 波以有限值向外传播) (11) \end{array}$

以上各式中, n均为流体外法向。

$ϕ^{Ι} = \frac{- i g A}{ω} \frac{c o s h k (h + z)}{c o s h k h} e^{i k x c o s θ}$ $ϕ^{Ι} = \frac{- i g A}{ω} \frac{c o s h k (h + z)}{c o s h k h} e^{i k x c o s θ}$

表示入射势。ω和k满足色散关系ω²=gktanhkh。

ϕ_α满足下列方程及边界条件:

$\begin{array}{l} Δ^{2} ϕ_{α} - (k s i n θ)^{2} ϕ_{α} = 0 (在流体中) (12) \\ \frac{\partial ϕ_{α}}{\partial n} - \frac{ω^{2}}{g} ϕ_{α} = 0 (在自由表面上) (13) \\ \frac{\partial ϕ_{α}}{\partial n} = 0 (在水底) (14) \\ \frac{\partial ϕ_{α}}{\partial n} = n_{α} (在浮体表面) (15) \\ ϕ_{α} (在无穷远处 ‚ 波以有限值向外传播) (16) \end{array}$ $\begin{array}{l} Δ^{2} ϕ_{α} - (k s i n θ)^{2} ϕ_{α} = 0 (在流体中) (12) \\ \frac{\partial ϕ_{α}}{\partial n} - \frac{ω^{2}}{g} ϕ_{α} = 0 (在自由表面上) (13) \\ \frac{\partial ϕ_{α}}{\partial n} = 0 (在水底) (14) \\ \frac{\partial ϕ_{α}}{\partial n} = n_{α} (在浮体表面) (15) \\ ϕ_{α} (在无穷远处 ‚ 波以有限值向外传播) (16) \end{array}$

n_α为流体的广义外法向量, 定义为:

$[n_{α}] = [\overset{⇀}{n} (\overset{⇀}{x} - \overset{⇀}{X}_{0}) \times \overset{⇀}{n}]^{Τ}$ $[n_{α}] = [\overset{⇀}{n} (\overset{⇀}{x} - \overset{⇀}{X}_{0}) \times \overset{⇀}{n}]^{Τ}$

式中 $\overset{⇀}{X}_{0}$ $\overset{⇀}{X}_{0}$ 为转动中心。

由于计算时边界不可能取在无穷远处, 因此边界条件 (11) 、 (16) 需进行近似处理。边界条件 (11) 可处理为以下形式:

$\begin{array}{l} \frac{\partial ϕ^{D}}{\partial n} - i k c o s θ ϕ^{D} + 2 i k c o s θ ϕ^{Ι} = 0 \\ (在边界 Ⅰ 处) (17) \\ \frac{\partial ϕ^{D}}{\partial n} - i k c o s θ ϕ^{D} = 0 \\ (在边界 Ⅱ 处) (18) \end{array}$ $\begin{array}{l} \frac{\partial ϕ^{D}}{\partial n} - i k c o s θ ϕ^{D} + 2 i k c o s θ ϕ^{Ι} = 0 \\ (在边界 Ⅰ 处) (17) \\ \frac{\partial ϕ^{D}}{\partial n} - i k c o s θ ϕ^{D} = 0 \\ (在边界 Ⅱ 处) (18) \end{array}$

边界条件 (16) 可近似为:

$\begin{array}{l} \frac{\partial ϕ_{α}}{\partial n} - i k c o s θ ϕ_{α} = 0 \\ (在边界 Ⅰ ‚ Ⅱ 处) (19) \end{array}$ $\begin{array}{l} \frac{\partial ϕ_{α}}{\partial n} - i k c o s θ ϕ_{α} = 0 \\ (在边界 Ⅰ ‚ Ⅱ 处) (19) \end{array}$

《2.2浮体的水动力学系数》

2.2浮体的水动力学系数

浮体的水动力学特性可由一组系数来描述。这组系数可经如下方法得到:通过2.1节中的控制方程和边界条件求出绕射势和辐射势, 然后利用式 (20) ～式 (23) 可得到作用在浮体上的波浪激励力、附加质量和辐射阻尼。

$\begin{array}{l} [F^{D}] ∶ F_{α}^{D} = i ρ ω \int_{S_{B}} \int ϕ^{D} n_{α} d S \\ (激励力向量) (20) \\ [f] ∶ f_{β α} = i ρ ω \int_{S_{B}} \int ϕ_{β} n_{α} d S \\ (辐射作用力矩阵) (21) \end{array}$ $\begin{array}{l} [F^{D}] ∶ F_{α}^{D} = i ρ ω \int_{S_{B}} \int ϕ^{D} n_{α} d S \\ (激励力向量) (20) \\ [f] ∶ f_{β α} = i ρ ω \int_{S_{B}} \int ϕ_{β} n_{α} d S \\ (辐射作用力矩阵) (21) \end{array}$

式中, f_βα表示简振模式β引起的在α模式下的水力学反作用。

$\begin{array}{l} [μ] : μ_{β α} = \frac{1}{ω} Ι m [f_{β α}] \\ (附加质量矩阵) (22) \\ [λ] : λ_{β α} = - Re [f_{β α}] \\ (辐射阻尼矩阵) (23) \end{array}$ $\begin{array}{l} [μ] : μ_{β α} = \frac{1}{ω} Ι m [f_{β α}] \\ (附加质量矩阵) (22) \\ [λ] : λ_{β α} = - Re [f_{β α}] \\ (辐射阻尼矩阵) (23) \end{array}$

将这些参数无因次化后即得到浮体的水动力学系数。

《3 数值模型》

3 数值模型

为适应复杂的分析域, 并保证合理的计算精度, 作者选择有限元方法作为分析工具。

《3.1分析域的空间离散》

3.1分析域的空间离散

生成合适的分析网格是有限元分析的前提。三角形单元以其良好的区域适应性成为二维有限元分析的常用网格单元。目前比较流行的三角形网格生成方法^[16]有:网格模板法、Delaunay三角化方法 (Delaunay triangulation——DT) ^[17] 、推进波前法 (advancing front method——AFM) ^[18]。作者采用基于推进波前法的网格生成算法。该方法生成的单元质量高, 种类多, 网格密度控制方便, 自动化程度高, 自适应能力强。为进一步提高单元质量, 在网格生成后采用拉普拉斯光滑技术^[19] 来加以优化。

《3.2控制方程的离散以及边界条件的引入》

3.2控制方程的离散以及边界条件的引入

为便于处理, 在确定2.1中定解问题的泛函时, 将各边界条件统一为如下形式:

$[\frac{\partial ϕ}{\partial n} + p ϕ + q] |_{Γ} = 0 (24)$ $[\frac{\partial ϕ}{\partial n} + p ϕ + q] |_{Γ} = 0 (24)$

这样满足控制方程 (5) 和边界条件 (24) 的泛函为:

$\begin{array}{l} J = \int_{D} \int (ϕ_{x}^{2} + ϕ_{z}^{2} + (k s i n θ)^{2} ϕ^{2}) d x d y + \\ \int_{Γ} (p ϕ^{2} + 2 q ϕ) d s (25) \end{array}$ $\begin{array}{l} J = \int_{D} \int (ϕ_{x}^{2} + ϕ_{z}^{2} + (k s i n θ)^{2} ϕ^{2}) d x d y + \\ \int_{Γ} (p ϕ^{2} + 2 q ϕ) d s (25) \end{array}$

内部点的泛函只需去掉边界曲线积分即可。

采用三角形单纯形单元离散后的控制方程和边界条件如下:

$\begin{array}{l} {(\frac{1}{2 A})^{2} (\begin{matrix} b_{i} & c_{i} \\ b_{j} & c_{j} \\ b_{m} & c_{m} \end{matrix}) (\begin{matrix} b_{i} & b_{j} & b_{m} \\ c_{i} & c_{j} & c_{m} \end{matrix}) \int_{D} \int d x d z + \\ \int_{D} \int (k s i n θ)^{2} [\begin{array}{l} Ν_{i} \\ Ν_{j} \\ Ν_{m} \end{array}] [Ν] d x d z + \int_{Γ} p [\begin{array}{l} Ν_{i} \\ Ν_{j} \\ Ν_{m} \end{array}] [Ν] d s} \\ [ϕ]^{e} = - \int_{Γ} q [\begin{array}{l} Ν_{i} \\ Ν_{j} \\ Ν_{m} \end{array}] d s (26) \end{array}$ $\begin{array}{l} {(\frac{1}{2 A})^{2} (\begin{matrix} b_{i} & c_{i} \\ b_{j} & c_{j} \\ b_{m} & c_{m} \end{matrix}) (\begin{matrix} b_{i} & b_{j} & b_{m} \\ c_{i} & c_{j} & c_{m} \end{matrix}) \int_{D} \int d x d z + \\ \int_{D} \int (k s i n θ)^{2} [\begin{array}{l} Ν_{i} \\ Ν_{j} \\ Ν_{m} \end{array}] [Ν] d x d z + \int_{Γ} p [\begin{array}{l} Ν_{i} \\ Ν_{j} \\ Ν_{m} \end{array}] [Ν] d s} \\ [ϕ]^{e} = - \int_{Γ} q [\begin{array}{l} Ν_{i} \\ Ν_{j} \\ Ν_{m} \end{array}] d s (26) \end{array}$

式中 $Ν_{i} = \frac{1}{2 A} (a_{i} + b_{i} x + c_{i} z)$ $Ν_{i} = \frac{1}{2 A} (a_{i} + b_{i} x + c_{i} z)$ 为形函数, A为单元面积。式 (26) 即为离散后的控制方程, 等式左边[ϕ]^e的系数即为单元刚度矩阵。

《3.3离散后形成的线性方程组的求解》

3.3离散后形成的线性方程组的求解

将单元刚度矩阵组装形成总体刚度后得到一系数为大型稀疏矩阵的线性方程组。有效地求解该方程组是整个数值分析的基础。在计算中必须把求解的时间复杂度和空间复杂度都降到一个合理的水平, 以便现有的软、硬件能够承受。这就要求选择合适的算法、设计有效的数据结构。作者采用基于Krylov子空间^[14]的BiCGStab算法 (H.A Van Der Vorst^[15], 1992) 。该算法的效率远高于基本迭代法且有良好的收敛性^[14]。为节约稀疏系数矩阵的存储空间, 作者采用了行压缩存储格式 (compressed sparse row, CSR^[14]) 。这种存储格式具有运算简单、操作方便、通用性好等优点。

《图2》

图2分析3种地形离散后的情况

Fig.2 Discretizations of analysis areas

《4 模型验证》

4 模型验证

作者所分析的3种地形如图2所示。为验证数值模型的正确性, 在分析复杂地形情况之前使用简单地形 (Layout 1) 情况下的解析解来验证数值解。图3为简单地形情况下数值解与解析解的比较。该图中各物理量为无因次化的水动力学系数, 意义分别为: $F_{α} = \frac{F_{α}^{D}}{2 ρ g r A} — α$ $F_{α} = \frac{F_{α}^{D}}{2 ρ g r A} — α$ 模式的广义波浪激励力系数;θ_α—α模式的广义波浪激励力系数的辐角; $m_{β α} = \frac{μ_{β α}}{ρ D}$ $m_{β α} = \frac{μ_{β α}}{ρ D}$ —简振模式β引起的α模式的附加质量系数; $l_{β α} = \frac{λ_{β α}}{ρ ω D}$ $l_{β α} = \frac{λ_{β α}}{ρ ω D}$ —简振模式β引起的α模式的辐射阻尼系数。其中r为浮体半径, A为入射波振幅, D为浮体的浸润面积。

二维情况下的运动模式为1 (x方向运动—横荡) 、3 (z方向运动—垂荡) 、5 (绕y轴转动—纵摇) 。因此在进行上述无因次化时, 若对应模式中有一个为5时, 需再除以r, 若对应模式中 (α, β) 均为5时, 需除以2 r², 以达到无因次化的目的。

《4.1误差分析》

4.1误差分析

从图3来看, 数值解与解析解基本一致。计算误差可归结于2.1节中开边界I、II的处理, 以及计算区域的网格分布不尽合理。

在数学上要求开边界Ⅰ、Ⅱ的位置要尽可能远离扰动区 (关心区域) , 最好是选在无穷远处, 以便局部扰动完全衰减。但在数值处理中, 如果开边界取的太远使得分析所需的计算时间和存储空间都大为增加, 同时由于格式的耗散或频散, 分析精度会随着分析区域的增大而大大降低。在实际分析时必须折衷考虑这两个因素。

在网格生成中, 自动化程度越高, 质量越难以控制。因此, 也不可避免地会对计算产生不良影响, 带来误差。

《5 波浪入射角及地形对浮体水动力学特性的影响》

5 波浪入射角及地形对浮体水动力学特性的影响

为更好地研究波浪入射角和地形的影响, 分析了地形占水深的比例分别为0.29和0.46两种情况 (浮体所占水深为0.4) 。实际的海底地形非常复杂, 为便于研究作者将地形加以简化 (见图2的Layout 2—地形2 和 Layout 3—地形3) 。这里比较了入射波角度分别为0°、60°时地形对浮体的水动力学特性的影响, 结果示于图4 (图中各物理量的意义同图3) 。

《5.1波浪入射角的影响》

5.1波浪入射角的影响

图4显示, 在kh>π即深水情况时, F₁和F₅对波浪入射角不再敏感, 而此前斜向波的这两个波浪激励力系数明显减小。斜向波的F₃在1<kh<10时也有所减小。相反, 各波浪力的辐角在kh<π时变化不大, 而此后随着kh的增加即进入深水, 斜向波的辐角比正向波的大得越多。

波浪入射角对附加质量和辐射阻尼的影响更为显著。斜向波情况下的附加质量不仅峰值有所增加, 而且几乎都有先增加后减少的趋势, 并且在kh增大到一定程度时与正向波的附加质量曲线相交。辐射阻尼在出现比正向波大得多的峰值后与其差别逐渐减小。

《5.2地形的影响》

5.2地形的影响

从图4中可以看出, 地形2和3对浮体的各水动力学系数都有影响, 地形3的浅化影响尤为显著。

《图3》

图3简单地形下数值解与解析解的比较

Fig.3 Comparisons between the numerical and the analytical methods for simple seabed

图4显示, 地形的影响在使得F₁和F₅达到更高峰值后迅速地减小, 各波浪激励力辐角变化与此相似。而F₃在地形的影响下出现先减小后增大的趋势, 在kh达到一定值后变化已不明显。

附加质量在地形的影响下先增大后减小, 此后随着kh增加影响变小。而辐射阻尼出现峰值明显增大并且前移的现象。

出现这些现象的原因可归结为:kh>6时所分析的问题已进入深水情况, 此时地形对波浪的影响变小。而此前, 由于地形引起的绕射和反射使得浮体的各水动力学系数都有一定程度的变化。

《6 结论》

6 结论

在验证了数值方法的正确性后, 计算了不同情况下浮体的水动力学系数, 分析了波浪入射角和地形对其水动力学特性的影响。数值分析不仅得到了与理论分析一致的趋势, 而且给出了影响的大小。这些结果为设计安全可靠、经济合理的海洋结构物提供了关键的水动力学数据, 有非常重要的价值。

《图4》

图4波浪入射角及地形对浮体水动力学特性的影响

Fig.4 Effects of incident wave angle and seabed on hydrodynamic characteristics of buoy

展示更多