半潜平台(semi-submersible)是开采3 000 m水深石油资源的重要装备，其结构设计必须保证其在整个使用期的绝对安全^[1-2]。上述严格的使用要求使得在设计半潜钻井平台时必须同时考虑以下3个方面：1)高精度的三维建模方法；2)考虑实体模型在不同波浪方向和频率作用下的水动力响应和波浪力计算方法^[3-9]；3)基于实体模型的全船有限元分析方法^[11-13]。

进行全船结构有限元分析的前提是准确预报作用在海工结构物湿表面下的流体压力分布包括静水压力和运动响应引起的波浪压力，这些分布力需要准确投影到进行结构分析的有限元模型上。

半潜平台是一个复杂的空间结构，在波浪的作用下其六个方向的刚体运动都很明显，在计算波浪压力分布时，需要同时考虑半潜平台的辐射运动和绕射运动的影响，是一个极为复杂的大规模数值分析系统^[14]。基于波浪辐射和绕射理论的边界元方法^[15]是目前求解海洋工程波浪与结构物相互作用的重要方法，为了获得良好的考虑不同波向和频率作用下半潜平台波浪压力分布的稳定数值解，研究中采用了高阶边界元方法并考虑了物体的几何对称性^[16-17]。

目前较先进的有限元建模软件包括GAMBIT, MSC Pastran和ANASYS等，这些软件可以根据不同的流体计算和结构分析的要求构建复杂的有限元模型；本研究中首次将半潜平台有限元建模、运动参数计算以及边界元数据输出统一成一个完整体系，同时发展了一个重要的新手段就是在半潜平台的几何奇异点处(角点或者物体表面相交的直角边处)按照几何连续性和物理连续性的概念布置双重和多重节点，在该位置处的节点具有不同的法向导数、属于不同的单元，但具有共同的位置坐标和速度势。数值实践表明，该方法极大的提高了波浪压力分析的准确性和数值稳定性，上述方法是实现半潜平台这种复杂空间结构波浪压力高精度结果的重要保证^[18]。本文依据波浪辐射、绕射理论，通过有限元模型发展了新的边界元分析几何数据模型，同时引入高阶边界元方法对半潜平台湿表面的动波压力进行了研究。

1 边界值问题 1.1 波浪辐射绕射基本理论

定义如图 1所示笛卡尔坐标系，该坐标系相对于无扰动的自由表面和物体处于静止状态，原点取在自由表面上，z轴垂直向上为正。假定流体为无粘性、不可压缩并且运动是无旋的，则流体速度可表示为速度势Φ的梯度。在波浪作用下，半潜平台的运动可假定为一个简谐振荡系统，流体由于半潜平台的摇荡运动产生的速度势可设为Φ=Re(Φe^－iωt)，式中Re表示实部，ω为入射波浪频率，t为时间，Φ为与时间变量无关的空间速度势，一般也为复数。这样所有的边界值问题都可以用复速度势Φ来表示，最后的计算结果为复数变量与时间因子e^－iωt的乘积。

速度势Φ的求解为定常问题，在整个流体域内也满足拉普拉斯方程。

$ {{\nabla }^{2}}\mathit{\Phi} =0 $

(1)

线性化的自由表面条件为

$ {{\mathit{\Phi} }_{Z}}-K\mathit{\Phi} =0,z=0 $

(2)

式中：K=ω²/g为深水波数, g为重力加速度。波浪入射速度势可表示为

$ {{\mathit{\Phi} }_{I}}=\frac{igA}{\omega }\frac{\cosh \left[ k\left( z+H \right) \right]}{\cosh kh}{{\text{e}}^{-\text{i}k\left( x\cos \beta +y\sin \beta \right)}} $

(3)

式中：A为入射波浪振幅，H为水深，β为入射波浪传播方向与x轴正向的夹角; 波数k为散射方程$\frac{{{\omega ^2}}}{g}$=ktanhkH的实根。根据波浪运动的线性假定，流体总速度势Φ可以分解为辐射势(Φ_R)和散射势(Φ_S)两部分，分别表示如下：

$ \left\{ \begin{align} & \mathit{\Phi} ={\mathit{{\Phi }}_{R}}+{\mathit{{\Phi }}_{S}} \\ & {{\mathit{\Phi} }_{R}}=-\text{i}\omega \sum\limits_{j=1}^{6}{{{\xi }_{j}}{{\mathit{\Phi} }_{j}}} \\ & {{\mathit{\Phi} }_{S}}={{\mathit{\Phi} }_{I}}+{{\mathit{\Phi} }_{D}} \\ \end{align} \right. $

(4)

式中：散射势(Φ_S)包含入射势(Φ_I)和绕射势(Φ_D)；ξ_j表示物体做辐射运动时在6个刚体自由度方向上的复数振幅；Φ_j为单位速度辐射势，Φ_D为入射势作用在固定物体上时产生的绕射势。在静平衡位置上做简谐摇荡运动的物体，其边界条件一般可由总速度势的各个分量表示，其中辐射势和散射势应满足以下边界条件：

$ \left\{ \begin{align} & {{\mathit{\Phi} }_{jn}}={{n}_{j}} \\ & {{\mathit{\Phi} }_{Sn}}=0 \\ \end{align} \right. $

(5)

(n₁, n₂, n₃)=n, (n₄, n₅, n₆)=x×n, x=(x, y, z)，单位向量n沿法线方向指向物面为正，x为物面坐标。同时辐射速度势Φ_j(j=1, …, 6)和绕射势Φ_D还必须满足波浪在无限远处传播的辐射条件：

$ \left\{ \begin{align} & \underset{R\to \infty }{\mathop{\lim }}\,\sqrt{R}\left( \frac{\partial {{\mathit{\Phi} }_{j}}}{\partial R}-\text{i}k{{\mathit{\Phi} }_{j}} \right)=0 \\ & \underset{R\to \infty }{\mathop{\lim }}\,\sqrt{R}\left( \frac{\partial {{\mathit{\Phi} }_{D}}}{\partial R}-\text{i}k{{\mathit{\Phi} }_{D}} \right)=0 \\ \end{align} \right. $

(6)

总波压力$\mathop P\limits^ -$为：

$ \begin{matrix} \overline{p}=-\rho gz-\rho \frac{\partial \left( \mathit{\Phi} {{\text{e}}^{-\text{i}\omega t}} \right)}{\partial t}=-\rho gz+\text{i}\rho \omega \mathit{\Phi} {{\text{e}}^{-\text{i}\omega t}}= \\ \ \ {{P}_{0}}+P{{\text{e}}^{-\text{i}\omega t}}={{P}_{0}}+\left( {{P}_{S}}+{{P}_{R}} \right){{\text{e}}^{-\text{i}\omega t}} \\ \end{matrix} $

(7)

式中：P₀为静压，P为动压，而散射波压力P_S为

$ {{P}_{S}}=\text{i}\rho \omega \left( {{\mathit{\Phi} }_{I}}+{{\mathit{\Phi} }_{D}} \right) $

(8)

辐射波浪力P_R为

$ {{P}_{R}}=\text{i}\rho \omega \left( -\text{i}\omega \sum\limits_{j=1}^{6}{{{\xi }_{j}}{{\mathit{\Phi} }_{j}}} \right)=\rho {{\omega }^{2}}\sum\limits_{j=1}^{6}{{{\xi }_{j}}{{\mathit{\Phi} }_{j}}} $

(9)

1.2 速度势积分方程

上述关于物体表面的辐射和散射速度势的边界值问题(1)～(6)可通过格林定理在物面边界上建立积分方程式的方法求解。其中物面上的辐射势满足的积分方程为:

$ \begin{matrix} {{\mathit{\Phi} }_{j}}\left( {{r}_{0}} \right)\left[ 1-K\int\limits_{{{S}_{{{f}'}}}}{G\text{d}s} \right]+ \\ \iint_{{{S}_{b}}}{\left[ {{\mathit{\Phi} }_{j}}\left( r \right)-{{\mathit{\Phi} }_{j}}\left( {{r}_{0}} \right) \right]}\frac{\partial G\left( r;{{r}_{0}} \right)}{\partial n}\text{d}s\text{=}\iint_{{{S}_{b}}}{{{n}_{j}}G\left( r;{{r}_{0}} \right)\text{d}s} \\ \end{matrix} $

(10)

而散射势满足的积分方程为

$ \begin{matrix} {{\mathit{\Phi} }_{S}}\left( {{r}_{0}} \right)\left[ 1-K\int\limits_{{{S}_{{{f}'}}}}{G\text{d}s} \right]+ \\ \iint_{{{S}_{b}}}{\left[ {{\mathit{\Phi} }_{S}}\left( r \right)-{{\mathit{\Phi} }_{S}}\left( {{r}_{0}} \right) \right]}\frac{\partial G\left( r;{{r}_{0}} \right)}{\partial n}\text{d}s\text{=}\iint_{{{S}_{b}}}{{{\Phi }_{I}}\left( {{r}_{0}} \right)} \\ \end{matrix} $

(11)

式中：r=(x, y, z)表示域点, r₀(ξ, η, ζ)表示源点, S_f'为物体内部水线面，S_b表示物体在静水中的湿表面，K为深水波数，格林函数G(r; r₀)指的是波浪脉动源，它表示位于r₀点的单位强度点源在r点处产生的速度势。格林函数满足自由表面和远场的波浪辐射条件，有限水深时格林函数可定义为：

$ \begin{matrix} \ G\left( r;{{r}_{0}} \right)=-\frac{1}{4\text{ }\!\!\pi\!\!\text{ }}\left[ \frac{1}{r}+\frac{1}{{{r}''}}+ \right. \\ 2\int_{0}^{\infty }{\left. \frac{\left( k+K \right)\cosh k\left( z+H \right)\cos k\left( \zeta +H \right)}{k\sinh kH-K\cosh kH}{{\text{e}}^{-kH}}{{\text{J}}_{0}}\left( kR \right)\text{d}k \right]} \\ {{\left( {{r}''} \right)}^{2}}={{\left( x-\xi \right)}^{2}}+{{\left( y-\eta \right)}^{2}}+{{\left( z+\zeta +2H \right)}^{2}} \\ \end{matrix} $

(12)

式中：J₀为零阶贝塞尔函数，k的积分路径取k值在实轴上的实根上半部分，从而满足波浪的远场辐射条件。格林函数的具体求解方法见文献[19]。

1.3 节点重置方法及方程离散求解

半潜平台的空间结构是非常复杂的，当采用边界元方法求解时，针对半潜平台在面-面相交形成的大范围的直角边和角点的情况，研究中采用在该几何奇异位置布置双重和多重节点的方法，例如在原来直角边的一个控制点处再增加一个控制点，变成直角边同一位置具有两个控制点，对于三个面构成的角点位置可以增加两个控制点，变成三个控制点。这些新的节点具有相同的空间坐标和速度势，但是具有不同的法向导数并且属于不同的单元。同时研究中采用6节点三角形单元和8节点四边形单元对半潜平台湿表面进行剖分。与重节点方法相匹配，在进行高阶单元网格剖分时，单元不能跨越奇异边界和节点，必须逐面剖分^[20]。由于半潜平台具有两个对称面，研究中引入了几何对称性的方法，提高了计算精度和矩阵的组装时间。空间复杂结构的波浪漂流力计算在短波时会有数值不稳定现象出现，研究中采用了Gauss直接方法求解高频波浪问题，而在其它波浪频率内采用LU分解方法。

2 半潜平台波浪压力分布研究波浪压力分布研究

半潜平台是由立柱和浮箱组成的，长宽比适中具有透空间隔的极为复杂的三维大型结构物，它通常具有两个对称面(关于x和y轴对称)。研究中根据其几何对称性的特征，在迎浪β=180°至横浪β=90°的第一象限波浪入射区域内，沿15°角的间隔共假定了7个方向的波浪入射，在每个波浪入射方向上考虑了21个波浪频率，频率间隔为Δω=0.05 rad/s，这样的计算体系能够完全涵盖半潜平台的所有方向。表 1为研究中采用的ISSC型半潜平台的基本设计和工作参数，其中半潜平台各立柱的位置坐标如表 2所示。

本研究首先利用Gambit软件对半潜平台进行三维实体建模和网格剖分，并通过自主开发的fortran算法对Gambit所输出的三维体网格进行处理，以形成边界元计算所需的经双重及多重节点重置得到的结构湿表面网格，之后通过自有的边界元计算算法^[21-23]对半潜平台波动压力进行计算。

研究中选取了对运动响应和波浪力影响明显的波浪频率ω=0.2、0.9 rad/s以及波向β=135°、180°进行了比较分析。图 2～5为ω=0.2 rad/s, β=135°、图 6～9为ω=0.2、0.9 rad/s, β=180°时入射势(Φ_I)、绕射势(Φ_D)、散射势(Φ_S)及总速度势(Φ)产生的动压分布图。由图 2～5可以发现，在长波ω=0.2rad/s及首斜浪β=135°情况下，总速度势产生的动压主要由入射势决定。总速度势(Φ)及入射势(Φ_I)产生的动压均为正值，在立柱A₂, A₃及B₃, B₄处最大，其中入射势波压在立柱近自由表面处波压为(9 800～10 400 Pa)，而总速度势产生的波压在立柱近自由表面处为(9 400～9 800 Pa)，入射势和总速度势产生的波压在A₅, B₁处最小，沿浮箱底部至立柱，波压力变化范围为(7 600～8 800 Pa)。绕射势(Φ_D)产生的波压最小并且有负压区产生，负压区主要分布在浮箱的底部、侧面以及立柱(A₁, A₄, A₅)以及(B₁, B₂, B₅)的背浪面。绕射势产生的最大正压区位于浮箱顶部，数值区间为(400～600 Pa)，其余正压区位于立柱A₃及B₃表面及以其为中心的外传波浪与其他圆柱的迎浪面，数值区间为(100～300 Pa)。

在迎浪β=180°情况下由图 6～9可以看到，在长波ω=0.2 rad/s时，总速度势产生的动压也是主要由入射势决定。总速度势(Φ)及入射势(Φ_I)和散射势(Φ_I+Φ_D)产生的动压均为正值并且由于半潜平台的几何对称性，各速度势产生的压力分布也是对称的。波动压力在立柱A₃及B₃处最大，其中入射势波压在立柱A₃及B₃处为(9 800～10 400 Pa)，而总速度势产生的波压在A₃及B₃处为(9 400～9 800 Pa)。入射势(Φ_I)、散射势(Φ_I+Φ_D)和总速度势(Φ)产生的波压在(A₁, A₅), (B₁, B₅)处最小，沿浮箱底部至立柱，波压力变化范围为(7 000～8 400 Pa)。绕射势(Φ_D)产生的波压最小并且有负压区产生，负压区主要分布在浮箱的底部、侧面以及立柱(A₁, A₅), (B₁, B₅)处，数值区间为(－500～0 Pa)。绕射势产生的最大正压区位于浮箱顶部，数值区间为(400～600 Pa)，其余正压区位于立柱A₃及B₃表面及以其为中心的外传波浪与其余圆柱的迎浪面，数值变化区间为(200～400 Pa)。

由图 10～13可以看到在首斜浪ω=0.9 rad/s, β=135°时，总速度势产生的动压主要由散射势(Φ_I+Φ_D)决定。总速度势(Φ)及入射势(Φ_I)和散射势(Φ_I+Φ_D)产生的动压在半潜平台不同位置出现量级相当的正负压区域。入射势(Φ_I)产生的压力变化区间(－18 000～18 000 Pa), 最大正压区位于(A₂, A₄), (B₂, B₄)，而最小负压区位于(A₁, A₃, A₅), (B₁, B₃)。绕射势(Φ_D)产生的压力变化区间为(－12 000～14 000 Pa), 最大正压区位于(A₅, A₃)，而最小负压区位于(B₂, B₄)。而总速度势(Φ)和散射势(Φ_I+Φ_D)产生的波压分布基本相同，在A₄处均为正压，而在(A₃, A₁)处为负压区，在其余立柱表面会同时存在正负压分布区。总速度势(Φ)和散射势(Φ_I+Φ_D)在浮箱表面产生的波压力变化范围为(－400～400 Pa)。

由图 14～17可以看到在迎浪ω=0.2 rad/s, β=180°时，总速度势产生的动压主要由散射势(Φ_I+Φ_D)决定。总速度势(Φ)及入射势Φ_I和散射势(Φ_I+Φ_D)产生的动压在半潜平台A区及B区对称分布，正负压同时出现在不同位置但量级相当。入射势(Φ_I)产生的压力变化区间(－11 000～11 000 Pa), 最大正压区位于(A₁, B₁)，而最小负压区位于(A₁, A₅), (B₁, B₅)。绕射势(Φ_D)产生的压力变化区间为(－8 000～11 000 Pa), 最大正压区和最小负压区同时位于立柱(A₂, A₅)和(B₂, B₅)处。总速度势(Φ)和散射势(Φ_I+Φ_D)产生的波压分布基本相同，压力变化区间为(－18 000～18 000 Pa)，正压主要分布在立柱(A₂, A₃), (B₂, B₃)迎浪面，而负压区分布在(A₁, B₁)背浪面及(A₅, B₅)迎浪面。在其余立柱表面会同时存在正负压分布区。总速度势(Φ)和散射势(Φ_I+Φ_D)在浮箱表面靠近立柱(A₅, B₅)处有最小负压存在，波压变化范围为(－800～200 Pa)。

3 结论

为了提高半潜平台波压力的计算精度，本文研究中考虑了半潜平台的几何对称性并应用了双重和多重节点的重置方法去除速度势求解过程中的几何奇异性。研究发现：

1)半潜平台的波压力主要由散射势和辐射势两部分组成，其中散射势产生散射波浪力，辐射速度势所产生的动压与运动响应相关。

2)ω=0.2 rad/s时无论迎浪β=180°还是首斜浪β=135°产生的波压都是正值，而且在起控制作用的散射速度势中入射波浪速度势起决定作用，绕射势产生的波压幅值很小。

3)在ω=0.9 rad/s时迎浪β=180°和首斜浪β=135°产生的波压在半潜平台湿表面特别是立柱表面，同时出现正负压分布区域，而且量值相当。整体而言，对于双浮筒立柱结构的半潜平台，控制波浪频率在ω=0.9 rad/s附近。