0 引　言

随着我国深远海养殖产业的快速发展，国内海工业界对深远海养殖平台进行了大量探索实践，研发了船式、半潜式、框架式等多种养殖平台形式。框架式平台因其用钢量少、建造工艺简单的特点，成为我国现阶段主流深远海养殖平台形式。目前，业界关于深远海养殖平台的研究更多地集中在渔网水动力载荷试验和数值模拟^{[1 − 4]}、平台主体水动力分析和定位系泊系统分析^{[5 − 8]}以及平台结构强度分析^{[9 − 12]}的研究；气隙性能研究更多地针对于传统海洋油气钻井、生产平台^{[13 − 18]}；对于深远海养殖平台的气隙性能研究较少。海洋平台气隙是指波浪表面到平台最下层甲板底部的垂直距离，设计时要求最下层甲板具有足够的高度，以保证平台气隙值为正，避免波浪砰击到甲板结构，引起结构损坏和平台工作人员伤亡^[13]。框架式养殖平台相对于油气平台干舷设计值更小，最下层工作甲板设置离水面更近，方便于进行养殖操作，但会导致平台初始静气隙过小，在环境载荷作用下容易出现负气隙，造成结构砰击损坏和作业人员安全风险。

框架式养殖平台结构型式与传统海洋油气钻井、生产平台差别较大，配备大型网衣附属结构，平台结构构件尺度更小，较小干舷的框架式养殖平台在波浪中水线面和湿表面会发生很大变化，因此其受到的环境载荷与传统油气平台存在较大差别。深远海养殖平台作为新兴海洋结构物，其载荷的复杂性导致业界对其在风浪流环境条件中的运动规律尚无统一、清晰的认识。采用传统频域线性气隙分析方法可能会忽略对养殖平台的气隙性能影响的关键因素，导致低估了最危险气隙情况。

为了提高深远海框架式养殖平台气隙计算结果准确性和可靠性，采用时域动态方法对某型深远海框架式养殖平台进行数值分析，对考虑瞬时湿表面、二阶波浪载荷、定位系泊系统和风浪载荷因素的气隙结果进行对比分析，确定框架式养殖平台气隙计算的关键非线性影响因素。

1 计算方法 1.1 方法概述

目标框架式养殖平台主体由板梁结构和细长杆件组成，配备尼龙网衣。本文假设网衣处于张紧状态，网衣变形可忽略。框架式养殖平台数值模型由湿表面面元网格单元、结构杆件莫里森单元和网衣等效杆件莫里森单元组成。首先，采用基于势流理论的面元法，对初始平均湿表面模型进行频域分析，得到平台平均湿表面对应的傅汝德-克雷洛夫力（FK力）、静水恢复力、绕射力和附加质量、辐射阻尼的频域解。再采用时域方法模拟平台在波浪中的运动，考虑瞬时入射波面并计算对应的非线性FK力和静水恢复力，每个时间步根据瞬时入射波面情况对平台湿表面进行修正并构建流体动力学模型，据此模型计算瞬态的静水恢复力、FK力、绕射力以及辐射力；根据莫里森方程计算结构杆件的粘性拖曳力和等效网衣杆件的拖曳力，将这些载荷施加在平台结构上，通过求解刚体动力学方程得到平台运动。

1.2 平台在波浪中运动频域分析

基于三维势流理论来求解框架式养殖平台的运动，采用速度势函数来描述流场。求解平台运动的本质就是求解流场中平台存在时的速度势^[17]。流场中总的速度势可表示为$ \mathrm{\varnothing }\left(x,y,z,t\right)=Re\left\{\varphi \left(x,y,z,t\right){e}^{i\omega t}\right\} $，总速度势可线性分解为：

$ {\varphi }_{T}\left(x,y,z\right)={\varphi }_{I}\left(x,y,z\right)+{\varphi }_{D}\left(x,y,z\right)+{\varphi }_{R}\left(x,y,z\right) 。$

(1)

式中：$ {\varphi }_{I}\left(x,y,z\right)={\zeta }_{a}{\varphi }_{0}\left(x,y,z\right) $为入射势；$ {\varphi }_{R}\left(x,y,z\right)= \sum _{j=1}^{6} \left[i\omega {\eta }_{j}{\varphi }_{j}\left(x,y,z\right)\right] $为绕射势；$ {\varphi }_{R}\left(x,y,z\right)=\sum _{j=1}^{6} $ $ \Bigr[i\omega {\eta }_{j} {\varphi }_{j} \left(x,y,z\right)\Bigr] $为辐射势。

通过求解这3个部分速度势，叠加得到总速度势。使用拉普拉斯方程$ {\nabla }^{2}\mathrm{\varnothing }\left(x,y,z,t\right)=0 $来描述流场中的总速度势。

入射波势为：

$ {\varphi }_{I}\left(x,y,z,t\right)=\frac{g{\zeta }_{a}}{{\omega }_{0}}{e}^{{{k}_{0}}z}\sin\left({k}_{0}x\cos\beta +{k}_{0}y\sin\beta -{\omega }_{\varepsilon }t\right)。$

(2)

式中：$ {\omega }_{0} $为波浪自然频率；$ {k}_{0} $为波浪自然波数；$ {\zeta }_{a} $为入射波波幅；$ g $为重力加速度；$ \beta $为浪向角；$ {\omega }_{\varepsilon } $为遭遇频率。

通过求解定解条件下的$ {\nabla }^{2}\mathrm{\varnothing }\left(x,y,z,t\right)=0 $可得到辐射势和绕射势，需要求解的方程组如下：

$ \left\{\begin{aligned} &{\nabla }^{2}\mathrm{\varnothing }\left(x,y,z,t\right)=0，\\ &\frac{{\partial }^{2}\varphi }{\partial {t}^{2}}+g\frac{\partial \varphi }{\partial z}=0，\\& \frac{\partial \varphi }{\partial \mathrm{n}}=v\cdot n，\\ &\nabla \varphi =0，\\ &\underset{R\rightarrow \mathrm{\infty }}{\lim } \sqrt{R}\left(\frac{\partial \varphi }{\partial R}-i\nu \varphi \right)=0。\\ \end{aligned} \right.$

(3)

波浪中平台刚体运动方程频域表达如下：

$ \left[\boldsymbol {M}+\mu \left(\omega \right)\right]\ddot{\eta }+\lambda \left(\omega \right)\dot{\eta }+C\eta =F\left(\omega \right)。$

(4)

式中：$ \boldsymbol{M} $为平台广义质量矩阵；$ C $为静水回复力系数；$ \mu \left(\omega \right) $为附加质量；$ \lambda \left(\omega \right) $为阻尼；$ F\left(\omega \right) $为波浪激励力。根据前述定解条件和方程组，可以求得流场中的各速度势，积分后得到运动方程两端的各项系数，再对方程进行求解，得到平台的频域运动响应。

1.3 平台在波浪中运动时域分析

采用时域方法求解养殖平台的运动，时域运动方程可以表示为^[19]：

${ \left[M+m\right]\ddot{\eta }\left(t\right)+\int _{-\mathrm{\infty }}^{t}\left[K\left(t-\tau \right)\right]\dot{\eta }\left(t\right){\mathrm{d}}\tau +C\eta \left(t\right)=F\left(t\right)+{F}_{ex}\left(t\right) 。}$

(5)

式中：$ M $和$ m $分别为平台质量和频率趋向无穷时$ \left(\omega \rightarrow \mathrm{\infty }\right) $对应的附加质量；$ K\left(t-\text{π} \right) $为延迟函数；$ C $为平台刚度；$ F\left(t\right) $为时域波浪激励力；$ {F}_{ex}\left(t\right) $为除波浪激励力以外的载荷，包括杆件拖曳力、网衣拖曳力、风力、流力和系泊力等。

时域运动方程的附加质量$ m $可用无穷频率处的频域附加质量表示：

$ {\mu }_{ij}\left(w\right)-{\mu }_{ij}\left(\mathrm{\infty }\right)=-\frac{1}{\omega }\int _{0}^{\mathrm{\infty }}{K}_{ij}\left(t\right)\sin\omega t{\mathrm{d}}t。$

(6)

通过傅里叶逆变换得延迟函数$ {K}_{ij}\left(t\right) $：

$ {K}_{ij}\left(t\right)= -\frac{2}{\text{π} }\int _{0}^{\mathrm{\infty }}\omega \left[{\mu }_{ij}\left(w\right)-{\mu }_{ij}\left(\mathrm{\infty }\right)\right]\sin\omega t{\mathrm{d}}\omega 。$

(7)

根据Cummins脉冲理论，依据式(8)直接对频域波浪力进行转换得到时域波浪力：

$ \left\{\begin{aligned} & F_{i}^{\left(1\right)}= \int _{0}^{t}h_{i}^{\left(1\right)}\left(t-\tau \right)\eta \left(\tau \right){\mathrm{d}}\tau，\\& F_{i}^{\left(2\right)}= \int _{0}^{t}\int _{0}^{t}h_{i}^{\left(2\right)}\left(t-{\tau }_{1},t-{\tau }_{2}\right)\eta \left({\tau }_{1}\right)\eta \left({\tau }_{2}\right){\mathrm{d}}{\tau }_{1}{\mathrm{d}}{\tau }_{2}。\\ \end{aligned} \right.$

(8)

其中：

${ \left\{\begin{aligned} &h_{i}^{\left(1\right)}\left(t\right)= \frac{1}{\text{π} }\int _{0}^{\mathrm{\infty }}H_{i}^{\left(1\right)}\left(\omega \right){e}^{iwt}{\mathrm{d}}\omega，\\ &h_{i}^{\left(2\right)}\left({t}_{1},{t}_{2}\right)= \frac{1}{2{\text{π} }^{2}}\int _{0}^{\mathrm{\infty }}\int _{0}^{\mathrm{\infty }}H_{i}^{\left(2\right)}\left({\omega }_{1},{\omega }_{2}\right){e}^{i\left({\omega }_{1}{t}_{1}+{\omega }_{2}{t}_{2}\right)}{\mathrm{d}}{\omega }_{1}{\mathrm{d}}{\omega }_{2}。\\ \end{aligned}\right. }$

(9)

$ H_{i}^{\left(1\right)}\left(\omega \right) $和$ H_{i}^{\left(2\right)}\left({\omega }_{1},{\omega }_{2}\right) $分别为平台受到的一阶和二阶波浪激励力的响应传递函数。

时域分析中每个时间步的绕射力可通过辐射力计算得到^[20]：

$ {\begin{split}{F}_{d}\left(t\right)=&\boldsymbol{E}\left(t\right)\sum _{{N}_{d}}^{m=1}\sum _{{N}_{m}}^{j=1}{a}_{jm}{F}_{d}\left({\omega }_{jm},{\beta }_{m}\right)\times\\&{e}^{i\left({k}_{jm}{X}_{g}\cos {\chi }_{m}+{k}_{jm}{Y}_{g}\sin {\chi }_{m}-{\omega }_{m}t+{\alpha }_{jm}\right)}。\end{split}}$

(10)

式中：$ \boldsymbol{E}\left(t\right) $为时间步$ t $处的欧拉旋转矩阵；$ {\beta }_{m}= {\chi }_{m}-{\theta }_{3}\left(t\right) $为波浪与平台间的相对角度；$ {\chi }_{m} $为浪向角；$ {\theta }_{3} $为平台艏摇角；$ {F}_{d}\left({\omega }_{jm},{\beta }_{m}\right) $为频率$ {\omega }_{jm} $、相对浪向角$ {\beta }_{m} $对应的单位波幅入射波引起的辐射力，辐射力可由基于求解拉普拉斯方程的频域方法计算得到；$ {N}_{d} $为频域分析时的浪方数量；$ {N}_{m} $为第$ m $个浪方的组成随机波的规则波数量。

1.4 气隙计算方法

平台气隙观测点处气隙最小值等于静气隙减去观测点垂向相对运动最大值^[21]，有

$ AG\left(x,y,t\right)={AG}_{0}\left(x,y\right)-\gamma \left(x,y,t\right)。$

(11)

式中：$ \left(x,y\right) $为气隙观测点坐标；$ AG\left(x,y,t\right) $为观测点气隙瞬时值；$ {AG}_{0}\left(x,y\right) $为观测点初始静气隙；$ \gamma \left(x,y,t\right) $为观测点处垂向相对运动。

观测点垂向相对运动与观测点处运动和观测点处波面升高有关：

$ \gamma \left(x,y,t\right)=\alpha \cdot \eta \left(x,y,t\right)-\delta \left(x,y,t\right) 。$

(12)

式中：$ \eta \left(x,y,t\right) $为观测点处波面升高；$ \delta \left(x,y,t\right) $为观测点处垂向运动；$ \alpha $为波浪不对称系数。

波面升高包含入射波面升高、扰动波面升高以及二阶波面升高成分，平台气隙与波面升高关系密切，业界相关研究已非常成熟^{[13 − 16]}。为了更清晰地分析瞬时湿表面、二阶载荷、定位系泊系统和风、流载荷等因素对气隙的影响，减少其他变量的影响，波面升高只考虑入射波面成分。

2 计算模型

目标框架式养殖平台坐标系原点为中纵剖面、中横剖面在基线处的交点，Y轴正向指向左舷，X轴正向指向艏部，Z轴正向垂直向上，如图1所示。平台总长为72 m、型宽为38 m、型深为15 m，设计吃水为13.80 m，排水量为1800 t，重心位置为X=0.00 m、Y=0.00 m、Z=4.90 m，回转半径为$ {R}_{xx}=17.50\;\mathrm{m} $、$ {R}_{yy}=23.60\;\mathrm{m} $、$ {R}_{zz}=28.20\;\mathrm{m} $，作业水深30 m。平台定位系泊系统由4角4根长度为190 m、链径58 mm的锚链组成，锚链与X轴夹角为60°，锚链单位长度空气中质量为58.03 kg/m、拉伸刚度为249000 kN。

养殖平台营运海域为珠江口外海，计算平台自存工况下的气隙，对应环境条件波浪有义波高为1.96 m，谱峰周期为5 s，采用Jonswap谱模拟，谱峰因子为1.9；流速为2 m/s；10 m高度处1 min平均风速41.6 m/s；短期预报模拟时长3 h。选取P1～P12气隙观测点见图1。利用商业软件AQWA进行框架式养殖平台频域、时域模拟分析，数值模型如图2所示。

3 计算结果 3.1 瞬时湿表面影响分析

为研究瞬时湿表面对平台气隙的影响，提取0°、90°浪向考虑平均湿表面和瞬时湿表面的气隙时域分析结果统计值对比如图3(a)所示。0°浪向最小气隙分布趋势一致，考虑瞬时湿表面后最小值显著减小（约0.4 m）；90°浪向最小气隙分布趋势发生了明显变化，考虑瞬时湿表面后，平台左舷背浪侧气隙显著减小（约0.8 m），平台右舷迎浪侧气隙变化较小（约0.1 m）。选取90°浪向P3、P9点气隙时历结果对比见图4，气隙时历变化趋势基本一致，背浪测P3变化趋势较迎浪侧P9波动更大，在时厉曲线极值处波动更为明显，且极值出现的时刻也发生了变化。由此可知，瞬时湿表面对平台气隙结果影响显著，考虑瞬时湿表面变化会明显减小气隙值。

3.2 二阶波浪载荷影响分析

为研究二阶波浪载荷对平台气隙的影响，提取0°、90°浪向考虑二阶波浪载荷和只考虑一阶波浪载荷的气隙时域分析结果统计值对比见图3(b)。0°浪向平台最小气隙分布变化很小，考虑二阶波浪载荷气隙减小约0.1 m；90°浪向平台最小气隙分布变化较大，考虑二阶波浪载荷平台左舷背浪侧气隙减小值（约0.4 m）大于右舷迎浪侧气隙减小值（约0.2 m）。选取0°浪向P6点和90°浪向P9点的气隙时历结果对比见图5(a)和图5(b)，0°浪向气隙时历变化很小，趋势和极值基本重合；90°浪向气隙时历变化趋势一致，但在极值处会出现波动，考虑二阶波浪载荷气隙最小值减小。由此可知，不同浪向二阶波浪载荷对气隙结果影响存在差别，0°浪向影响较小，90°浪向使得气隙极值显著减小。

3.3 定位系泊系统影响分析

为研究定位系泊系统对平台气隙影响，将定位系泊系统预张力分别设置为1T（锚链松弛状态）、5T（锚链悬链线状态）和40T（锚链张紧状态）。提取0°、90°浪向气隙时域分析结果统计值对比见图3(c)：随着预张力增大，气隙呈现减小趋势；0°浪向气隙减小值约为0.1 m，90°浪向气隙减小范围为0.1～0.2 m。分别选取0°浪向P6点和90°浪向P9点的气隙时历结果对比见图5(c)和图5(d)：锚链松弛状态和悬链线状态时历曲线基本重合，张紧状态气隙时历曲线的趋势和极值都发生较大变化。定位系泊系统不同预张力导致平台整体刚度发生变化，从而引起平台运动和气隙时历结果的改变；松弛和悬链状态，定位系泊系统对平台垂向运动限制较小，因此气隙结果相差不大；张紧状态定位系泊系统对平台垂向运动限制较大，引起平台气隙时历发生较大变化。

3.4 风载荷影响分析

为研究风载荷对平台气隙影响，考虑平台在无风载荷、1倍风载荷和2倍风载荷作用的3类工况进行时域分析。0°、90°环境载荷方向气隙时域分析结果统计值对比见图6(a)：考虑风载荷作用，平台不同区域气隙变化的情况不同；再结合图1气隙观测点坐标分布发现，0°风载作用，迎风面观测点P4/P5/P6/P7/P8气隙增大，背风面观测点P1/P2/P10/P11/P12气隙减小，中横轴上观测点P3/P4基本无变化；90°风载作用，迎风面观测点P7/P8/P9/P10/P11气隙增大，背风面观测点P1/P2/P3/P4P5气隙减小，中纵轴上观测点P6/P12基本无变化。

为进一步分析风载荷影响：剔除波动成分干扰，提取平台垂荡、横摇和纵摇运动时历统计平均值，见表1。将无风载荷的气隙平均值$ {AG}_{0} $叠加风载荷引起的平台垂荡$ H $和倾角$ R/P $的变化，根据几何关系计算出气隙值$ {AG}_{C1} $/$ {AG}_{C2} $：0°浪向$ {AG}_{C1}={AG}_{0}+\Delta H-X\times {\tan }^{-1} (P) $；90°浪向$ {AG}_{C1}={AG}_{0}+\Delta H+Y\times{\tan }^{-1} (R) $。计算气隙值$ {AG}_{C1} $/$ {AG}_{C2} $与时域模拟统计结果$ {AG}_{1} $/$ {AG}_{2} $相等，见图7(a)和图7(b)。垂向运动引起的气隙变化均值$ {DV}_{Z1/Z2}=\Delta H-X\times{\tan }^{-1} (P) $或$ {DV}_{Z1/Z2}=\Delta H+Y\times {\tan }^{-1} (R) $与时历统计结果气隙差值$ {DV}_{1/2} $（由考虑风载荷结果减去无风载荷结果得到）对比见图8(a)和图8(b)。除了中横轴和中纵轴上的气隙观测点，风载荷导致平台垂向运动变化引起的气隙变化均值占到平台气隙极值变化约80%（0°方向）、约100%（90°方向）。由此可知，风载荷对平台气隙的影响主要原因是风载荷引起平台垂向自由度运动的变化，平台下沉和倾斜引起了气隙观测点几何位置的改变，导致气隙最终结果变化。

3.5 流载荷影响分析

为研究流载荷对平台气隙影响，考虑平台在无流载荷、1倍流载荷和2倍流载荷作用的3类工况进行时域分析。0°、90°环境载荷方向气隙时域分析结果统计值对比见图6(b)，再结合图1气隙观测点坐标分布发现：0°流载作用，平台整体气隙都减小，迎流面观测点P4/P5/P6/P7/P8气隙相较于背风面观测点P1/P2/P10/P11P12气隙减小更多，与风载影响规律有区别；90°流载作用，迎流面观测点P7/P8/P9/P10/P11气隙增大，背流面观测点P1/P2/P3/P4/P5气隙减小，与90°风载影响一致。

进一步分析其中原因，对比表1中风载荷和流载荷的结果，流载荷引起平台垂荡运动变化更大，垂荡运动减小约0.34~0.61 m，风载荷影响约0.10~0.19 m；流载荷引起平台倾斜（横摇/纵摇）变化更小，0°方向纵倾角约为0.20°，风载荷影响约0.49～1.00 °；90°方向横倾角约为1.05～1.95 °，风载荷影响约2.17～4.15 °。采用3.4节风载荷气隙平均值和运动几何关系计算方法，得到计算气隙值$ {AG}_{C1} $/$ {AG}_{C2} $与时域模拟统计结果$ {AG}_{1} $/$ {AG}_{2} $，二者也是相等，与风载荷影响规律一致，见图7(c)和图7(d)；由图8(c)和图8(d)可知，流载荷导致平台垂向运动变化引起的气隙变化均值占到平台气隙极值变化约70%～95%（0°方向）、约25%、88%（90°方向）。由此可知，流载荷对平台气隙的影响原因与风载荷一样，引起平台垂向自由度运动的变化，导致气隙的改变；但流载荷使得平台下沉运动更大，倾斜角度更小。

4 结　语

深远海框架式养殖平台气隙对海上养殖作业安全和平台结构完整非常重要，框架式养殖平台因其结构特点，需要采用时域方法进行气隙分析考虑其中非线性因素的影响，上述研究表明：

1）相对于平均湿表面假设，采用时域方法考虑平台瞬时湿表面变化，平台气隙时历结果峰值/谷值出现波动且对应时刻发生了变化；瞬时湿表面变化对平台气隙结果影响显著，会极大减小平台气隙值。

2）不同浪向平台二阶波浪载荷对气隙结果影响存在差别：0°浪向对平台气隙结果影响较小，90°浪向平台使得气隙显著减小。定位系泊系统不同锚链预张力引起平台整体刚度变化，导致平台气隙时历发生改变，对平台气隙极值结果影响较小，但气隙极值出现的时刻变化较大。

3）风、流载荷导致平台垂向相对运动（垂荡和横摇/纵摇）变化，平台下沉、倾斜引起几何位置改变，导致气隙极值结果变化。相较于风载荷，流载荷使得平台下沉运动更大、倾斜角度更小。