筏式海水养殖设施在波浪作用下的水弹性分析

张光明，储 备，张 尧，张新曙

(1. 上海交通大学 海洋工程国家重点实验室，上海 200240；2. 高新船舶与深海开发装备协同创新中心，上海 200240；3. 自然资源部海洋减灾中心，北京 100194)

近年来，水产品逐渐成为人类获取优质蛋白质的重要来源。据联合国粮农组织2020 年报告：世界鱼类产量在2018年已经达到1.8亿t，其中0.3亿t来自海洋水产养殖，而中国贡献了35%的鱼类产量[1]。海水养殖主要分为筏式养殖和深水网箱养殖，筏式养殖设施广泛分布在我国福建和浙江沿海地区，筏式养殖的产量占海水养殖总产量的33%。

根据相关文献描述[2]，柬埔寨很早就采用了渔业网箱养殖方式，至今已经经过160多年的发展。挪威的深海养殖技术在世界上处于领先地位，其深水网箱的结构型式多样，主要包括可移动框构型、环型、张力腿型[3]深水网箱。而美国海上养殖设施中应用更多的是刚性更强的网箱，例如锚拉式或者崇碟型的网箱以及碟型样式的养殖网箱[4]。大概在20世纪70年代我国才开始使用近海类型网箱养殖，这种类型的构造比较简易，抵抗风浪的能力比较弱。到20世纪80年代开始进行抗风浪类型的网箱研究，但是成果不多，自从90年代末期引进了挪威产的养殖网箱，我国在该领域的研究才开始进入迅速发展时期[5]。不同于深水网箱，筏式养殖设施主要分布在我国近海，在我国海水养殖业中占据着重要的位置。如图1所示，在福建省连江县的木制筏式设施在台风作用下中遭受严重破坏。然而目前对该设施的破坏机理分析较少，其破坏产生的主要原因尚不明确。

图1 损坏的传统木制筏式设施(2019年8月拍摄于福建省连江县)Fig. 1 The damaged traditional wooden raft facilities (The photo is taken in Lianjiang, Fujian Province in August 2019)

对浮筏结构和深海网箱结构的研究有许多相似性，譬如结构都涉及水弹性问题且都有锚绳系统。Fredriksson等[6]开发利用有限元方法来研究高密度聚乙烯(high density polyethylene，简称HDPE)材料构成的网箱结构承载能力。吴常文等[7]对深水网箱抗风浪流性能进行验证试验，发现风力主要影响扶栏系统,流速主要影响网具系统,而波浪对网箱整体结构均有明显影响。曾启东等[8]建立了网衣与其在波浪中的三维数值模型，来研究波浪经过网衣后的变化规律，通过与物理模型试验结果的对比验证了该方法的可行性。对于浮筏养殖设施来说，水弹性问题和系泊系统是整个结构系统安全性的重要影响因素。

目前，对筏式养殖设施动力响应的研究并不充分，而涉及该结构水弹性问题的研究则更少。筏式养殖设施与浮桥结构有一定的相似性，对于一座浮桥的研究通常需要包括对水弹性进行有关的分析。对简化的浮桥结构模型，Viuff等[9]提出了其在非规则波中计算分析的方法。Fu等[10]在三维势流理论和有限元方法的基础上，提出了浮桥在非规则波中的水弹性响应数值模型。Viuff等[11]先分别用AQWA和WADAM进行浮筒的水动力系数求解，然后分别施加到SIMO-RIFLEX和OrcaFlex两个海洋工程应用软件建立的模型中，来对比研究浮桥的水弹性动力响应，分析了不同商业软件对浮桥问题研究的误差和可行性。文中所研究的筏式养殖设施如图2所示，纵横交错的筏式养殖设施相对于长线型的浮桥结构更加复杂。Chu 等[12]对与文中相同的筏式养殖设施在规则波与非规则波下的水弹性响应以及锚泊系统非线性的影响进行了研究。而文中则基于实际应用中的两种常见筏式养殖设施(木制浮筏与HDPE型浮筏)，研究走板材料(浮筏刚度)对其抗风浪性能的影响。

由于走板材料的不同，浮筏在波浪中的破坏程度也会不同。首先利用势流面元软件WAMIT对浮子的水动力系数进行求解，基于海洋工程动力分析软件OrcaFlex对使用不同走板材料的浮筏结构在波浪中的水弹性响应进行数值模拟和分析，比较了两类不同走板材料结构上的水动力载荷。

1 数值计算模型

1.1 浮筏结构动力响应方程

筏式养殖设施可以看作一个有一定弹性的梁单元系统，而这个梁单元系统是由多个提供浮力的刚性浮子连接在线单元的走板上组成的。利用OrcaFlex软件[13]对筏式养殖设施进行建模计算，其计算模型如图2所示，浮子采用的是6自由度浮体(6D buoy)，可看成一个刚体，且可提供所需浮力。通过在浮子上设置附加质量系数和阻力系数来分别考虑惯性力和黏性力的影响。该浮筏结构不考虑波浪二阶力的影响。水深为10 m，锚绳在4个方向分别布置2根，总共8根，每个方向的2根锚绳一端连接最靠近结构角点处的浮子，另一端固定在海底，锚绳连接点的水平距离为50 m，锚绳总长为52 m。

图2 筏式养殖设施模型示意Fig. 2 Sketch of the whole structure modelled in OrcaFlex

理论上求解筏式养殖设施结构的时历运动响应需要利用卡明斯方程，如式(1)所示：

(1)

由于浮子尺寸很小，而且根据Chu等[12]的结果，浮子的附加质量系数在所关心的波浪周期内变化不大，波浪阻尼系数相对于附加质量系数很小。因此，忽略式(1)中的延迟函数项，采用Morison公式求解浮子在波浪中水动力。

(2)

其中，ρ表示流体密度；Ad表示浮子截面积；U，af分别表示浮子中心处流体速度与加速度；CD表示浮子在波浪中的拖曳力系数。从而，整个浮筏结构运动方程可以简化为:

(3)

1.2 附加质量计算

浮子的附加质量系数通过势流软件WAMIT[14]求解得到。WAMIT是MIT的Newman教授团队开发的软件，它通过基于频域的面元模型来求解浮体的辐射绕射问题。图3(a)表示单个浮子模型和主尺度，浮子长度为0.9 m，直径为0.65 m，设计吃水d为0.2 m，oxyz为浮子局部坐标系，x方向为轴向，y方向为法向，z方向垂直向上；图3(b)表示单个浮子在WAMIT中的面元模型；图3(c)中定义了波浪的入射方向，全局坐标系OXYZ，浪向角β，该坐标系与图2(b)中的坐标系重合。

图3 浮子的主尺度、面元模型及坐标系Fig. 3 The principal dimensions, panel model of the buoy and the incident wave directions

单个浮子可以看成浮于海面的水平圆柱，在WAMIT中采用高阶面元法来求解浮子的水动力系数。基于势流理论，速度势Φ需满足拉普拉斯方程：

∇2Φ=0

(4)

总的速度势Φ可以写为：

Φ=Re(φeiωt)ω

(5)

其中，φ为空间速度势，Re表示取实部，ω是物体运动频率，t表示时间。空间速度势的自由面条件为：

φz-Kφ=0，z=0

(6)

其中，K=ω2/g表示无限水深波数。入射波速度势定义为：

(7)

其中，k表示有限水深入射波的波数，为ω2/g=ktanhkH的实根，A表示入射波的波幅，H为水深。

将空间速度势φ分解成辐射势φR和绕射势φD两部分：

φ=φR+φD

(8)

(9)

φD=φ0+φs

(10)

其中，ξj是浮体6自由度运动的幅值，φj表示相应的辐射势，j=1,2，……，6。速度势φs表示当浮体固定时由入射波引起的散射速度势，而入射波速度势和散射波速度势合成了绕射势φD。辐射势需在无穷远处满足辐射条件。辐射势和绕射势的物面分别如下：

φjn=nj

(11)

φDn=0

(12)

其中，向量n为物面的法向量，n=(n1,n2,n3)，x=(x,y,z)，x×n=(n4,n5,n6)。

浮体附加质量和波浪阻尼系数可通过式(13)计算：

(13)

其中，Sb表示平均湿表面。附加质量和波浪阻尼系数的无量纲形式为：

(14)

其中，p=3对应i,j=1,2,3;p=4对应i(j)=1,2,3,j(i)=4,5,6;p=5对应i,j=4,5,6。

在对浮子水动力系数进行网格收敛性分析时，设置了3组不同尺度的面元模型，面元尺寸分别为0.20 m，0.10 m和0.05 m。图4表示浮子无因次横荡附加质量随周期的变化，如图所示网格大小为0.10 m时结果已收敛，后续采用0.10 m的网格计算附加质量。

图4 浮子横向附加质量系数随周期的变化Fig. 4 The added mass of the buoy in sway with respect to wave periods

1.3 阻力系数计算

浮子所受的黏性力对其动力响应有较大影响。计算浮子在波浪中的阻力系数时需要注意以下两点：流体运动的周期性和自由面的影响。由于自由面的影响，浮子的阻力系数可能与Sarpkaya等[15]的浸没圆柱在振荡流中的阻力系数不同。依据挪威船级社规范DNV-RP-C205[16]里建议的方法，当浮子的运动幅值比结构的特征尺度大时，结构在静水中的摇荡阻尼系数可以用作结构在波浪中的阻力系数。因此，浮筏设施的浮子阻力系数可以通过浮子在纵荡方向的阻尼系数来估算。令浮子在纵荡方向进行强迫振动x(t)=X0cos(2πt/T)，X0表示振荡运动的幅值，T是振荡的周期。使用商业软件STAR-CCM+[17]计算浮子的阻力系数。选择求解雷诺平均N-S方程(RANS)并采用k-ω SST湍流模型，用流体体积方法(volume of fluid，简称VOF)来捕捉自由液面。

对浮子在横荡方向和纵荡方向的受力时历曲线进行傅里叶变换可以得到其所在对应方向的阻力系数。利用莫里森方程拟合得到浮子的阻力系数，作用在浮子上的纵向力可以表示为：

(15)

式中：V是浮子的体积，Ax是浮子在x方向的湿表面投影面积。在计算模拟中经过3个振动周期后，浮子的附加质量系数Ca和阻力系数CD趋于平稳。当各系数达到稳定后，对其取平均值，得到浮子在轴向和径向的阻力系数分别为0.58和0.23。这两个系数是后续在OrcaFlex中计算的设置系数。

2 结果分析

2.1 浮筏计算参数设置

浮筏设施各组成部分的类型和参数如表1所示，为了对比不同材料刚度对浮筏结构响应的影响，对连接走板设置了3种不同的材料刚度，HDPE材料和木材走板(WOOD)，以及设置了一半木板刚度的走板(HWOOD)，其刚度参数见表2。走板的截面形状如图5所示。表1中走板的直径为软件OrcaFlex中视图的直径，计算载荷是利用截面形状的参数来后处理的。连接浮子可以传递相交走板的载荷，其体积和质量可以忽略。每两个浮子间的连接走板长度多数为0.78 m，只有浮筏内部Y方向为1.08 m。

表1 浮筏设施组成部分的建模类型和参数Tab. 1 The parameters of components of the floating raft

表2 不同走板的材料及其对应参数Tab. 2 The parameters of different plates

比较不同走板材料对浮筏结构动力响应的影响。设置了不同的非规则波条件如表3所示。非规则波采用Jonswap波浪谱，其中参数参见Isherwood等[18]研究中的取值，波浪过零周期Tz从5～9 s每间隔0.5 s取值，有效波高2 m，浪向角0°。每组计算模拟筏式养殖设施在非规则波中1 h的结构动力响应时历。

图5 走板横截面模型示意Fig. 5 The sketch of the cross section of the plate

表3 非规则波参数Tab. 3 The parameters of the irregular waves

2.2 不同走板刚度的浮筏结构动力响应

锚绳所受到的张力过大导致的锚绳断裂和走板上受力超过其承受极限而导致的破坏是浮筏结构在波浪中损坏失效的两种主要形式。计算不仅需要考虑锚绳所受到的张力影响，也要考虑连接走板上所受的力和力矩的影响。由软件可以计算得到拉力，剪力，弯矩和扭矩等值，然后根据经验公式[19]来估算走板截面上的正应力和剪应力：

(16)

式中：σpl和τpl分别为走板横截面上的正应力和剪应力；FTpl，FSpl，Mpl和TSpl分别为其横截面上的拉力，剪力，弯矩和扭矩；h0为走板高度，b0为宽度，d0为其壁厚，如图5中所示；Ap和Iyy分别表示其横截面积和横截面在y轴方向的惯性矩。走板横截面正应力是由其横截面上所受到的拉力和弯矩合成的最大应力，剪应力是由其所受到的剪力和扭矩合成的最大应力。

图6表示不同材料刚度的浮筏结构其锚绳和走板上载荷的均方差(STD)随波浪过零周期的变化，均方差能够反映结构上的载荷在时历曲线中的总体变化情况。从图6中可以很直观地看到两种不同材料对浮筏结构响应所带来的影响，HDPE走板材料下的锚绳张力和走板上的拉力、剪力、弯矩以及扭矩的均方差明显小于另外两种木材走板的计算结果。随着波浪周期的增大，各走板材料下锚绳的张力、走板截面上的剪力和扭矩先增大随后减小，这些载荷的均方差在波浪周期6 s左右达到最大值。而走板截面上拉力和弯矩随着波浪周期的增大逐渐变小，在短波条件下(Tz=5 s)拉力和弯矩有最大值。如图6(d)所示，从载荷均方差的相对大小可以看出，材料刚度对结构响应中走板截面上的弯矩影响最显著。而从材料的特性可知，木板的刚度最大，而木板制成的浮筏在实际海况中最容易破坏，这说明在实际海况中，所用的材料对走板截面上的弯矩承受能力有很大的影响，从而很可能对结构的安全性能产生很大的影响。

图6 锚绳和走板上载荷的均方差随波浪周期的变化Fig. 6 The STD of the loads of the structure in irregular waves with respect to Tz

根据式(16)可得到走板截面的正应力与剪应力。图7中展示了不同非规则波下走板截面的正应力和剪应力的均方差。图7(a)和(b)中，木质走板下的应力明显比HDPE材料走板下的应力要大。图7(a)中所示同一波浪条件下，走板刚度越大其截面所受正应力的均方差越大。另外，走板截面所受的正应力明显大于其所受剪应力。所以走板材料刚度越小时，走板上所受到的弯矩越小，则其正应力也越小，结构越不容易破坏。由于HDPE材料具有很好的柔性而木板很硬几乎是刚体，所以当浮筏受到波浪载荷时，柔性材料可以很好缓解结构受力不均所带来的弯矩过大问题，而木板则产生很大的弯矩，结果表明木制结构浮筏可能更容易遭受破坏。

图7 走板截面的应力均方差随波浪周期的变化Fig. 7 The STD of the stress of the plates in irregular waves with respect to Tz

图8展示了非规则波下浮筏结构中心处的水平位移和垂直位移均方差随波浪周期的变化。从图8中可以看出，走板材料的不同，结构的位移也会有差异，木制走板材料的水平位移比HDPE大，前者垂向位移比后者小。随着波浪周期的增大，结构中心的水平位移先增大后减小，在波浪周期6.5 s左右，水平位移均方差有最大值，垂向位移几乎随着波浪周期的增大而增大。结构的水平位移总体要比垂向位移大许多，几乎是2倍关系。结构在波浪中的水平位移越大导致锚绳所受的张力也越大，浮筏结构走板截面上的各受力和力矩以及应力也越大，结构越容易破坏。HDPE材料制成的走板浮筏结构不仅在水平位移上要比木板小得多，其锚绳所受的张力、木板上所受的力和力矩以及应力也更小。在同等的破坏强度下，综合来看采用刚度更小、更具有柔性的HDPE材料作为浮筏走板材料更加安全可靠。从以往台风影响下近海浮筏渔排受损情况来看，HDPE材料制成的浮筏比传统木板制成的浮筏渔排更能抵抗风浪，这也和文中的结果一致。

图8 浮筏中心的位移均方差随波浪周期的变化Fig. 8 The STD of the motion of the structure in irregular waves with respect to Tz

2.3 浮筏动力响应的极值估计

根据Faltinsen[20]的研究，海洋平台等结构物在随机海况中的响应极值可以通过经验估算得到大概值，前提是这个值的数学统计结果有足够长的时历曲线结果，且这个统计结果要满足瑞利分布，最可能出现的最大值估计(MPM)可以通过式(17)估算：

(17)

式中：σm(X)表示变量X的均方差，td表示监测的持续时间，Tz表示波浪的过零周期。

由前面的分析可知，浮筏结构在波浪周期为6 s时所受载荷比较大，因此对波高为2 m，波浪周期为6 s，浪向角为0°的浮筏位移进行最大值估计分析。为了方便表示，文中定义M1=|MAX′-MEAN|，M2为式(17)所求，Errorm=(M2-M1)/M1。其中MEAN表示时历曲线的平均值，MAX′表示离平均值最远的值，M1表示时历曲线中离平均值最远距离的绝对值，M2表示变量的最大值估计MPM，Errorm表示最大值估计M2与M1误差百分比。

图9为浮筏中心处位移的时历曲线，不同材料刚度的走板位移时历曲线没有表现出很大的差异。图10表示浮筏中心位移的功率谱，水平位移功率谱的最大值约为垂向的5倍，说明水平方向的能量更大，而木板WOOD材料的功率谱值比HDPE材料的值更大，反映出在同样海况下木板结构的水平位移会比HDPE的大。图11是浮筏中心处位移的直方图分布情况，直方图中使用木板情况下浮筏中心位移绝对值大于2 m的N(X)数量比使用HDPE情况下的更多，这些结果反映了和功率谱一致的规律，即木板材料刚度下结构出现较大水平位移的可能性更大。根据公式(17)的最大值估计与时历曲线的最大值进行统计分析得到表4中的结果。表4中结果表明HDPE材料下的时历曲线中位移最大值小于木板的情况。垂向位移的最大值估计MPM结果与实际最大值的误差约为1%，水平位移的最大值估计结果误差在16%左右，这是由锚绳系统的非线性特性导致的。

图9 非规则波下浮筏中心的位移时历曲线(Tz=6 s)Fig. 9 The time history curves of the structural motion (Tz=6 s)

图10 非规则波下浮筏中心的位移功率谱(Tz=6 s)Fig. 10 The power spectrum of the structural motion (Tz=6 s)

图11 非规则波下浮筏中心的位移分布直方图(Tz=6 s)Fig. 11 The histogram of the structural motion(Tz=6 s)

表4 非规则波中不同走板材料的结构位移统计分析Tab. 4 The structural motion in the irregular waves with different plate materials

3 结 语

对非规则波下采用不同走板材料浮筏设施的水弹性响应进行了计算和比较。利用WAMIT计算了浮子的附加质量系数，利用CFD软件计算了浮子的阻力系数，然后使用OrcaFlex对采用HDPE和木板两种不同走板材料的3种不同刚度浮筏结构进行了动力响应计算，得出如下结论：

1) 随着波浪周期的增大，不同材料刚度下锚绳的张力、走板截面上的剪力和扭矩先增大随后减小，波浪周期Tz在6 s左右这些载荷的均方差达到最大值。而走板截面上拉力和弯矩随着波浪周期的增大逐渐变小，在短波条件下拉力和弯矩有最大值；材料刚度对走板截面上的弯矩影响最大。

2) 同一波浪条件下，走板刚度越大其截面所受正应力的均方差越大，且正应力明显大于剪应力，走板材料刚度越小时走板上受到的弯矩越小，则其正应力也越小，这表明当浮筏在波浪中运动时，柔性材料可以很好缓解结构受力不均所带来的弯矩过大问题。

3) 相同波浪条件下，当走板采用木板材料时浮筏结构的水平运动最大幅值比使用HDPE材料的大，结构在波浪中的水平运动最大幅值越大则锚绳所受的张力也越大，结构走板截面上的各受力和力矩以及应力也越大，结构也越容易破坏。

4) 垂向位移的最大值估计MPM结果与直接模拟最大值的误差约为1%，水平位移的最大值估计结果误差约为16%。综合来看，在同等的材料破坏强度下，采用刚度更小、更具有柔性的HDPE作为浮筏的走板材料可能更加安全可靠。