CN102368279A - 一种深海养殖网箱群组水动力响应数值模拟方法 - Google Patents

Abstract

本发明公开了一种深海养殖网箱群组水动力分析的数值模拟方法，包括以下步骤：建立浮架模型；建立网衣模型；建立锚绳模型；建立浮球模型；联立运动微分方程求解。本发明采用集中质量法以及刚体运动学原理模拟深水网箱群组结构在波浪作用下的水动力响应，采用本方法可以获得深水网箱群组结构的锚绳张力和浮架运动；根据我们的物理实验验证发现数值模拟得到的浮架运动以及锚绳张力与实际吻合。本发明能够分析多网箱组成的网箱群组的水动力响应，能够对组合式网箱锚绳系统进行优化。

Description

一种深海养殖网箱群组水动力响应数值模拟方法

技术领域

本发明属于水产养殖工程、水利工程和海洋工程技术领域，特别涉及一种在海洋环境荷载作用下深海养殖网箱群组水动力响应数值模拟方法。

背景技术

早在19世纪70年代，中国就开始利用网箱进行水产养殖活动；在上世纪八、九十年代，网箱水产养殖得到了快速的发展。截至2006年，全国有超过100万只网箱，养殖了30多种鱼类。在中国，网箱养殖的鱼每年大约30万吨，海洋养殖网箱已经成为了水产养殖的支柱产业。但是中国的传统网箱大部分都分布在近岸，这会引起一系列的问题：首先，传统网箱不能抵抗强风浪的侵袭，而分布在拥挤的近岸，近岸养殖密度超过了水体环境的承载极限，从而造成水体的污染以及网衣本身的污染；其次，大部分的近岸水体受到来自海岸线上的工业污染，这将导致养殖鱼类的品质下降，甚至死亡，降低了养殖效率，从而限制了水产养殖的发展。深水网箱养殖成为水产养殖可持续发展的必然选择。深海网箱养殖有以下优点：高科技及自动化，使用年限长以及抗污染、抗风浪能力强，大容量，经济效益更好，风险小，低成本，鱼类产品的高品质以及高价格。虽然深海养殖网箱有诸多优势，但是也同样面临很大的挑战，深海海洋环境相对近岸更加恶劣，网箱受到海洋环境荷载作用力也更大，为了设计合适的深海网箱锚绳系统，需要开展网箱水动力特性的研究。为了研究网箱的水动力特性，学者们开展了广泛的研究，这些研究主要集中在网衣以及单个网箱的水动力特性，对于组合式网箱群组的水动力研究很少，而深海网箱主要是以组合式网箱群组的形式存在，因此开展组合式网箱群组的水动力特性研究对于深海网箱系统的设计有重要的科学意义。

发明内容

为了辅助深水网箱群组结构的设计，本发明的目的是开展深水网箱群组结构在波浪作用下的水动力特性，分析深水网箱群组结构在波浪荷载作用下的锚绳受力以及浮架运动，分析网箱在不同布置条件下的锚绳受力情况，对网箱锚绳系统进行优化。

为了实现上述目的，本发明的技术方案如下：一种深海养殖网箱群组水动力分析的数值模拟方法，包括以下步骤：

A、建立浮架模型

浮架包括两个同心圆环形浮管、扶手和连接构件，浮架一般漂浮于水面之上，两个浮管是主要的受力构件；为了分析浮架的受力，浮架被简化成两个浮管；

为了计算作用于浮管的外力，浮管被分成许多的微段；局部坐标系n-τ-v固定于每一个浮管微段，n轴和τ轴分别为微段的法向和切向，v轴垂直于微段所在的平面；整个浮管受到的外力可以通过每一个微段上的外力求和得到；

因为浮管的管径远小于波浪的波长，采用修正的莫里森公式计算作用于浮管微段上的波浪力，该公式考虑了波浪与浮管微段之间的相对运动，表示如下：

$F_n=\frac{1}{2}{C}_{\mathrm{Dn}}\mathrm{\ρ}{A}_n|{\stackrel{\→}{u}}_n-{\stackrel{\stackrel{\→}{\·}}{R}}_n|\·({\stackrel{\→}{u}}_n-{\stackrel{\stackrel{\→}{\·}}{R}}_n)+\mathrm{\ρ}{V}_0{\stackrel{\→}{a}}_n+C_{\mathrm{mm}}\mathrm{\ρ}{V}_0({\stackrel{\→}{a}}_n-{\stackrel{\→}{\stackrel{\·\·}{R}}}_n)$

$F_{\mathrm{\τ}}=\frac{1}{2}{C}_{\mathrm{D\τ}}\mathrm{\ρ}{A}_{\mathrm{\τ}}|{\stackrel{\→}{u}}_{\mathrm{\τ}}-{\stackrel{\stackrel{\→}{\·}}{R}}_{\mathrm{\τ}}|\·({\stackrel{\→}{u}}_{\mathrm{\τ}}-{\stackrel{\stackrel{\→}{\·}}{R}}_{\mathrm{\τ}})+\mathrm{\ρ}{V}_0{\stackrel{\→}{a}}_{\mathrm{\τ}}+C_{\mathrm{m\τ}}\mathrm{\ρ}{V}_0({\stackrel{\→}{a}}_{\mathrm{\τ}}-{\stackrel{\→}{\stackrel{\·\·}{R}}}_{\mathrm{\τ}})---\left(1\right)$

$F_v=\frac{1}{2}{C}_{\mathrm{Dv}}\mathrm{\ρ}{A}_v|{\stackrel{\→}{u}}_v-{\stackrel{\→}{\stackrel{\·}{R}}}_v|\·({\stackrel{\→}{u}}_v-{\stackrel{\stackrel{\→}{\·}}{R}}_v)+\mathrm{\ρ}{V}_0{\stackrel{\→}{a}}_v+C_{\mathrm{mv}}\mathrm{\ρ}{V}_0({\stackrel{\→}{a}}_v-{\stackrel{\→}{\stackrel{\·\·}{R}}}_v)$

其中，F_n、F_τ和F_v分别表示波浪力在n、τ和v方向上的分量；

和

分别表示水质点在n、τ和v方向上的速度分量；

和

分别表示浮管微段在n、τ和v方向上的速度分量；和

分别表示水质点在n、τ和v方向上的加速度分量；

和

分别表示浮管微段在n、τ和v方向上的加速度分量；ρ是水的密度；V₀是浮管微段排开水的体积；A_n、A_τ和A_v分别是浮管微段在n、τ和v方向上的投影面积；C_Dn、C_Dτ和C_Dv分别是n、τ和v方向上的拖曳力系数；C_mn、C_mτ和C_mv分别是n、τ和v方向上的附加质量系数；

浮管被视为刚体，用六个自由度来描述浮管的运动，纵荡、横荡和升沉用来描述三个平动，纵摇、横摇和回转用来描述三个转动；为了建立浮管的运动方程，定义了两套坐标系；整体坐标系Oxyz以及局部坐标系Gabc；根据牛顿第二定律，在整体坐标系下，三个平动方程为：

${\stackrel{\·\·}{x}}_G=\frac{1}{m_G}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_{x_i},$ ${\stackrel{\·\·}{y}}_G=\frac{1}{m_G}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_{y_i},$ ${\stackrel{\·\·}{z}}_G=\frac{1}{m_G}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{F}_{z_i}---\left(2\right)$

其中，

和

是作用于浮架微段上的外力分量，和

是浮架中心的加速度，N是浮管微段数，m_G表示浮管的质量；

用欧拉方程来描述三个转动方程，在局部坐标系下，三个转动方程如下：

$I_a\frac{\∂{\text{\ω}}_a}{\∂t}+(I_c-I_b){\mathrm{\ω}}_c{\mathrm{\ω}}_b=M_a,I_b\frac{\∂{\mathrm{\ω}}_b}{\∂t}+(I_a-I_c){\mathrm{\ω}}_a{\mathrm{\ω}}_c=M_b$ (3)

$I_c\frac{\∂{\mathrm{\ω}}_c}{\∂t}+(I_b-I_a){\mathrm{\ω}}_a{\mathrm{\ω}}_b=M_c$

其中，下标a、b和c分别表示局部坐标系下的坐标轴a、b和c，I_a、I_b和I_c分别表示惯量沿三个坐标轴的分量，ω_a、ω_b和ω_c分别表示转动角速度的三个分量，M_a、M_b和M_c分别表示力矩的三个分量；

B、建立网衣模型

采用集中质点法建立网衣的模型，网衣假定由无质量的弹簧连接的质点组成，质点位于每一个网目目脚的中间点和端点；根据牛顿第二定律，网衣质点的运动方程为：

$M\stackrel{\stackrel{\→}{\·\·}}{R}=M\frac{{\∂}^2\stackrel{\→}{R}}{\∂t^2}={\stackrel{\→}{F}}_D+{\stackrel{\→}{F}}_I+\stackrel{\→}{T}+\stackrel{\→}{B}+\stackrel{\→}{W}---\left(4\right)$

其中，和

分别是拖曳力和惯性力，

是质点的加速度，

是网目目脚张力，

是浮力，

是重力，M是质点的质量；

C、建立锚绳模型

锚绳被简化成一系列的单元和质点，假定锚绳单元为圆柱体，局部坐标系τ-η-ξ定义在每一个锚绳单元上，η轴位于τ轴和水质点速度

所在的平面内，单元中点处的水质点速度分解成τ和η分量，因此，作用于每一个单元上的外力也被分解成τ和η分量；在整体坐标系下，局部坐标系τ-η-ξ轴的单位向量分别表示为e_τ＝(x_τ，y_τ，z_τ)，e_η＝(x_η，y_η，z_η)和e_ξ＝(x_ξ，y_ξ，z_ξ)。

作用于锚绳单元上的外力表示为：

$F_{\mathrm{\τ}}=-\frac{1}{2}\mathrm{\ρ}{C}_{\mathrm{d\τ}}\mathrm{Dl}|\stackrel{\·}{\mathrm{\τ}}-e_{\mathrm{\τ}}\·V|(\stackrel{\·}{\mathrm{\τ}}-e_{\mathrm{\τ}}\·V)+\mathrm{\ρ}\∀{\stackrel{\text{\→}}{a}}_{\mathrm{\τ}}+C_{\mathrm{m\τ}}\mathrm{\ρ}\∀({\stackrel{\→}{a}}_{\mathrm{\τ}}-e_{\mathrm{\τ}}\·\stackrel{\·}{V})$

$F_{\mathrm{\η}}=-\frac{1}{2}\mathrm{\ρ}{C}_{\mathrm{d\η}}\mathrm{Dl}|\stackrel{\·}{\mathrm{\η}}-e_{\mathrm{\η}}\·V|(\stackrel{\·}{\mathrm{\η}}-e_{\mathrm{\η}}\·V)+\mathrm{\ρ}\∀{\stackrel{\→}{a}}_{\mathrm{\η}}+C_{\mathrm{m\η}}\mathrm{\ρ}\∀({\stackrel{\→}{a}}_{\mathrm{\η}}-e_{\mathrm{\η}}\·\stackrel{\·}{V})---\left(5\right)$

$F_{\mathrm{\ξ}}=-\frac{1}{2}\mathrm{\ρ}{C}_{\mathrm{d\ξ}}\mathrm{Dl}|\stackrel{\·}{\mathrm{\ξ}}-e_{\mathrm{\ξ}}\·V|(\stackrel{\·}{\mathrm{\ξ}}-e_{\mathrm{\ξ}}\·V)+\mathrm{\ρ}\∀{\stackrel{\→}{a}}_{\mathrm{\ξ}}+C_{\mathrm{m\ξ}}\mathrm{\ρ}\∀({\stackrel{\→}{a}}_{\mathrm{\ξ}}-e_{\mathrm{\ξ}}\·\stackrel{\·}{V})$

其中，F_τ、F_η和F_ξ表示作用于锚绳单元的波浪力，G_Dτ、C_Dη和C_Dξ表示拖曳力系数，D是锚绳直径，l是锚绳单元的长度，

和

表示质点的速度分量，

和

表示质点的加速度分量；

其中，水动力系数采用下式计算：

$C_{\mathrm{Dn}}=\left\{\begin{array}{cc}0.0& ({\mathrm{Re}}_n\&le;0.1)\\ \frac{0.45+5.93}{{\left({\mathrm{Re}}_n\right)}^{0.33}}& (0.1<{\mathrm{Re}}_n\&le;400)\\ 1.27& (400<{\mathrm{Re}}_n\&le;{10}^5)\\ 0.3& ({\mathrm{Re}}_n>{10}^5)\end{array}\right.---\left(6\right)$

$C_{\mathrm{D\&tau;}}=\left\{\begin{array}{cc}\frac{1.88}{{\left({\mathrm{Re}}_n\right)}^{0.74}}& (0.1<{\mathrm{Re}}_{\mathrm{\&tau;}}\&le;100.55)\\ 0.062& ({\mathrm{Re}}_{\mathrm{\&tau;}}>100.55)\end{array}\right.---\left(7\right)$

其中，Re_n＝ρ|V_Rn|D/μ，μ是水的粘性系数；

计算了锚绳单元上的外力之后，将外力均匀的分布到与之相邻的质点上，锚绳质点的运动方程表示如下：

$m_i{a}_i=\underset{j=1}{\overset{\mathrm{count}}{\mathrm{\&Sigma;}}}({\text{T}}_j+W_j+B_j+F_j)---\left(8\right)$

其中，下标i表示质点编号，下标j表示与质点相邻的单元编号，count表示与质点相邻的单元数；

D、建立浮球模型

浮球漂浮在水面之上，受到波浪力的作用，作用于浮球上的波浪力采用下式计算：

${\stackrel{\&RightArrow;}{F}}_B={\stackrel{\&RightArrow;}{F}}_D+{\stackrel{\&RightArrow;}{F}}_I=\frac{1}{2}\mathrm{\&rho;}{C}_{D}A\frac{{\stackrel{\&RightArrow;}{V}}_{\mathrm{RB}}\left|{\stackrel{\&RightArrow;}{V}}_{\mathrm{RB}}\right|}{2}+\mathrm{\&rho;}{\&ForAll;}_B{C}_m\frac{\&PartialD;{\stackrel{\&RightArrow;}{V}}_{\mathrm{RB}}}{\&PartialD;t}+\mathrm{\&rho;}{\&ForAll;}_B\frac{\&PartialD;{\stackrel{\&RightArrow;}{V}}_B}{\&PartialD;t}---\left(9\right)$

其中，C_D为拖曳力系数，C_m为浮架质量系数，A是浮球在水质点速度方向上的投影面积，

是浮球的入水体积；拖曳力是雷诺数Re的函数，查表1可以得到拖曳力系数取值，表中未列出的数值采用线性插值得到：

表1浮球的拖曳力系数

(x_B，y_B，z_B)为浮球中心的坐标，浮球入水深度为：

Δh＝η(x，y，t)-(z_B(t)-r)    (10)

其中，η是(x_B，y_B，z_B)处的波面，能够通过线性波浪理论获得，r是浮球的半径；

浮球的投影面积计算如下：

$A_x=A_y=\left\{\begin{array}{cc}\frac{\mathrm{\&pi;}}{4}{D}^2& (\mathrm{\&Delta;h}\&GreaterEqual;D)\\ \frac{1}{8}{D}^2(\mathrm{\&theta;}-\sin \mathrm{\&theta;})& (D/2<\mathrm{\&Delta;h}\&le;D)\\ \frac{\mathrm{\&pi;}}{4}{D}^2-\frac{1}{8}{D}^2(\mathrm{\&theta;}-\sin \mathrm{\&theta;})& (0<\mathrm{\&Delta;h}\&le;D/2)\\ 0& (\mathrm{\&Delta;h}<0)\end{array}\right.---\left(11\right)$

$A_z=\left\{\begin{array}{cc}\frac{\mathrm{\&pi;}}{4}{D}^2& (\mathrm{\&Delta;h}\&GreaterEqual;D/2)\\ \mathrm{\&pi;}({\left(\frac{D}{2}\right)}^2-{(\frac{D}{2}-\mathrm{\&Delta;h})}^2)& (0<\mathrm{\&Delta;h}\&le;D/2)\\ 0& (\mathrm{\&Delta;h}<0)\end{array}\right.---\left(12\right)$

其中，D是浮球半径，θ是浮球与水面相交的弦对应的中心角，通过下式可得

$\mathrm{\&theta;}=2{\cos }^{-1}\left(\frac{D/2-\mathrm{\&Delta;h}}{D/2}\right)---\left(13\right)$

E、运动微分方程求解

上述方程(2)、(3)、(4)和(8)构成了个网箱群组的运动微分方程组，采用Runge-Kutta-Vener六阶数值方法求解该运动微分方程组，从而得到整个网箱群组结构在每一时刻的运动响应以及锚绳受力情况。

与现有的技术相比，本发明具有以下有益效果：

1、本发明采用集中质量法以及刚体运动学原理模拟深水网箱群组结构在波浪作用下的水动力响应，采用本方法可以获得深水网箱群组结构的锚绳张力和浮架运动；根据我们的物理实验验证发现数值模拟得到的浮架运动以及锚绳张力与实际吻合。

2、本发明能够分析多网箱组成的网箱群组的水动力响应，能够对组合式网箱锚绳系统进行优化。

附图说明

本发明共有附图12张，其中：

图1是浮架微段示意图。

图2是浮架局部坐标示意图。

图3是网衣模型示意图。

图4是锚绳模型示意图。

图5是浮球模型示意图。

图6是模型实验的俯视图。

图7是数值模拟与模型实验中的锚绳张力。

图8是数值模拟与模型实验中的两个浮架的垂直运动(第一个网箱)。

图9是数值模拟与模型实验中的两个浮架的垂直运动(第二个网箱)。

图10是数值模拟与模型实验中的两个浮架的水平运动(第一个网箱)。

图11是数值模拟与模型实验中的两个浮架的水平运动(第二个网箱)。

图12是采用本数值模型研究不同的网箱布置形式对锚绳张力分布的影响。

图中：1、锚绳，2、浮球，3、浮架，4、网衣，5、沉子。

具体实施方式

下面结合附图对本发明进行进一步地描述。如图1-12所示，采用本发明的方法进行的模拟如下：

对于浮架3的模拟，采用刚体运动学原理建立浮架3的模型。浮管微段模型如图1所示；定义在浮管微段上的局部坐标系以及整体坐标系如图2所示；应用集中质量法模拟网衣4，网衣4模拟为由无质量的弹簧连接的集中质量点见图3；锚绳1离散成一系列的单元和质点见图4，将作用于每一个单元上的外力均匀的分布到与之相连的质点上；浮球2被模拟成一个质点，浮球2模型示意图如图5所示。建立了整个网箱锚绳1结构模型的运动方程之后，采用Runge-Kutta-Vener六阶数值方法求解运动方程，得到整个网箱锚绳1结构的锚绳1受力和浮架3运动。

模型试验在大连理工大学海岸和近海工程国家重点实验室的港池中进行。试验水池长56米，宽34米，试验水深1米。试验中的网箱锚绳1的具体参数如表2所示，模型布置如图6所示。

表2网格式锚锭网箱系统参数

数值模拟与实验结果比较：数值模拟得到的浮架3水平和垂直运动，锚绳1张力与实验结果进行对比，如图7-11所示，数值模拟的结果与物理模型实验的结果吻合的很好，说明该数值模型能够很好的模拟网箱锚绳1结构在波浪作用下的浮架3运动情况，以及锚绳1受力情况。

采用该数值模型分析了两种不同的网箱布置形式见图12，并分析了不同的波浪方向作用下锚绳1张力分布情况。图12的结果表明，当波浪入射方向为0度时，两种布置形式下的锚绳1最大张力相差很小，但是当入射角为45度时，不同的布置形式将产生不同的锚绳1张力分布。采用本模型能够对网箱锚绳1系统进行优化，设计一个合适的锚绳1形式会有利于网箱锚绳1结构的安全。

Claims (1)

1.一种深海养殖网箱群组水动力分析的数值模拟方法，其特征在于：包括以下步骤：

A、建立浮架(3)模型

浮架(3)包括两个同心圆环形浮管、扶手和连接构件，浮架(3)一般漂浮于水面之上，两个浮管是主要的受力构件；为了分析浮架(3)的受力，浮架(3)被简化成两个浮管；

为了计算作用于浮管的外力，浮管被分成许多的微段；局部坐标系n-τ-v固定于每一个浮管微段，n轴和τ轴分别为微段的法向和切向，v轴垂直于微段所在的平面；整个浮管受到的外力可以通过每一个微段上的外力求和得到；

因为浮管的管径远小于波浪的波长，采用修正的莫里森公式计算作用于浮管微段上的波浪力，该公式考虑了波浪与浮管微段之间的相对运动，表示如下：

$F_n=\frac{1}{2}{C}_{\mathrm{Dn}}\mathrm{\&rho;}{A}_n|{\stackrel{\&RightArrow;}{u}}_n-{\stackrel{\stackrel{\&RightArrow;}{\&CenterDot;}}{R}}_n|\&CenterDot;({\stackrel{\&RightArrow;}{u}}_n-{\stackrel{\stackrel{\&RightArrow;}{\&CenterDot;}}{R}}_n)+\mathrm{\&rho;}{V}_0{\stackrel{\&RightArrow;}{a}}_n+C_{\mathrm{mm}}\mathrm{\&rho;}{V}_0({\stackrel{\&RightArrow;}{a}}_n-{\stackrel{\&RightArrow;}{\stackrel{\&CenterDot;\&CenterDot;}{R}}}_n)$

$F_{\mathrm{\&tau;}}=\frac{1}{2}{C}_{\mathrm{D\&tau;}}\mathrm{\&rho;}{A}_{\mathrm{\&tau;}}|{\stackrel{\&RightArrow;}{u}}_{\mathrm{\&tau;}}-{\stackrel{\stackrel{\&RightArrow;}{\&CenterDot;}}{R}}_{\mathrm{\&tau;}}|\&CenterDot;({\stackrel{\&RightArrow;}{u}}_{\mathrm{\&tau;}}-{\stackrel{\stackrel{\&RightArrow;}{\&CenterDot;}}{R}}_{\mathrm{\&tau;}})+\mathrm{\&rho;}{V}_0{\stackrel{\&RightArrow;}{a}}_{\mathrm{\&tau;}}+C_{\mathrm{m\&tau;}}\mathrm{\&rho;}{V}_0({\stackrel{\&RightArrow;}{a}}_{\mathrm{\&tau;}}-{\stackrel{\&RightArrow;}{\stackrel{\&CenterDot;\&CenterDot;}{R}}}_{\mathrm{\&tau;}})---\left(1\right)$

$F_v=\frac{1}{2}{C}_{\mathrm{Dv}}\mathrm{\&rho;}{A}_v|{\stackrel{\&RightArrow;}{u}}_v-{\stackrel{\&RightArrow;}{\stackrel{\&CenterDot;}{R}}}_v|\&CenterDot;({\stackrel{\&RightArrow;}{u}}_v-{\stackrel{\stackrel{\&RightArrow;}{\&CenterDot;}}{R}}_v)+\mathrm{\&rho;}{V}_0{\stackrel{\&RightArrow;}{a}}_v+C_{\mathrm{mv}}\mathrm{\&rho;}{V}_0({\stackrel{\&RightArrow;}{a}}_v-{\stackrel{\&RightArrow;}{\stackrel{\&CenterDot;\&CenterDot;}{R}}}_v)$

其中，F_n、F_τ和F_v分别表示波浪力在n、τ和v方向上的分量；

和

分别表示水质点在n、τ和v方向上的速度分量；

和

分别表示浮管微段在n、τ和v方向上的速度分量；和分别表示水质点在n、τ和v方向上的加速度分量；和

分别表示浮管微段在n、τ和v方向上的加速度分量；ρ是水的密度；V₀是浮管微段排开水的体积；A_n、A_τ和A_v分别是浮管微段在n、τ和v方向上的投影面积；C_Dn、C_Dτ和C_Dv分别是n、τ和v方向上的拖曳力系数；C_mn、C_mτ和C_mv分别是n、τ和v方向上的附加质量系数；

浮管被视为刚体，用六个自由度来描述浮管的运动，纵荡、横荡和升沉用来描述三个平动，纵摇、横摇和回转用来描述三个转动；为了建立浮管的运动方程，定义了两套坐标系；整体坐标系Oxyz以及局部坐标系Gabc；根据牛顿第二定律，在整体坐标系下，三个平动方程为：

${\stackrel{\&CenterDot;\&CenterDot;}{x}}_G=\frac{1}{m_G}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_{x_i},$ ${\stackrel{\&CenterDot;\&CenterDot;}{y}}_G=\frac{1}{m_G}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_{y_i},$ ${\stackrel{\&CenterDot;\&CenterDot;}{z}}_G=\frac{1}{m_G}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{F}_{z_i}---\left(2\right)$

其中，

和

是作用于浮架(3)微段上的外力分量，

和

是浮架(3)中心的加速度，N是浮管微段数，m_G表示浮管的质量；

用欧拉方程来描述三个转动方程，在局部坐标系下，三个转动方程如下：

$I_a\frac{\&PartialD;{\text{\&omega;}}_a}{\&PartialD;t}+(I_c-I_b){\mathrm{\&omega;}}_c{\mathrm{\&omega;}}_b=M_a,I_b\frac{\&PartialD;{\mathrm{\&omega;}}_b}{\&PartialD;t}+(I_a-I_c){\mathrm{\&omega;}}_a{\mathrm{\&omega;}}_c=M_b$ (3)

$I_c\frac{\&PartialD;{\mathrm{\&omega;}}_c}{\&PartialD;t}+(I_b-I_a){\mathrm{\&omega;}}_a{\mathrm{\&omega;}}_b=M_c$

其中，下标a、b和c分别表示局部坐标系下的坐标轴a、b和c，I_a、I_b和I_c分别表示惯量沿三个坐标轴的分量，ω_a、ω_b和ω_c分别表示转动角速度的三个分量，M_a、M_b和M_c分别表示力矩的三个分量；

B、建立网衣(4)模型

采用集中质点法建立网衣(4)的模型，网衣(4)假定由无质量的弹簧连接的质点组成，质点位于每一个网目目脚的中间点和端点；根据牛顿第二定律，网衣(4)质点的运动方程为：

$M\stackrel{\stackrel{\&RightArrow;}{\&CenterDot;\&CenterDot;}}{R}=M\frac{{\&PartialD;}^2\stackrel{\&RightArrow;}{R}}{\&PartialD;t^2}={\stackrel{\&RightArrow;}{F}}_D+{\stackrel{\&RightArrow;}{F}}_I+\stackrel{\&RightArrow;}{T}+\stackrel{\&RightArrow;}{B}+\stackrel{\&RightArrow;}{W}---\left(4\right)$

其中，

和

分别是拖曳力和惯性力，

是质点的加速度，

是网目目脚张力，

是浮力，

是重力，M是质点的质量；

C、建立锚绳(1)模型

锚绳(1)被简化成一系列的单元和质点，假定锚绳(1)单元为圆柱体，局部坐标系τ-η-ξ定义在每一个锚绳(1)单元上，η轴位于τ轴和水质点速度所在的平面内，单元中点处的水质点速度分解成τ和η分量，因此，作用于每一个单元上的外力也被分解成τ和η分量；在整体坐标系下，局部坐标系τ-η-ξ轴的单位向量分别表示为e_τ＝(x_τ，y_τ，z_τ)，e_η＝(x_η，y_η，z_η)和e_ξ＝(x_ξ，y_ξ，z_ξ)；

作用于锚绳(1)单元上的外力表示为：

$F_{\mathrm{\&tau;}}=-\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{d\&tau;}}\mathrm{Dl}|\stackrel{\&CenterDot;}{\mathrm{\&tau;}}-e_{\mathrm{\&tau;}}\&CenterDot;V|(\stackrel{\&CenterDot;}{\mathrm{\&tau;}}-e_{\mathrm{\&tau;}}\&CenterDot;V)+\mathrm{\&rho;}\&ForAll;{\stackrel{\text{\&RightArrow;}}{a}}_{\mathrm{\&tau;}}+C_{\mathrm{m\&tau;}}\mathrm{\&rho;}\&ForAll;({\stackrel{\&RightArrow;}{a}}_{\mathrm{\&tau;}}-e_{\mathrm{\&tau;}}\&CenterDot;\stackrel{\&CenterDot;}{V})$

$F_{\mathrm{\&eta;}}=-\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{d\&eta;}}\mathrm{Dl}|\stackrel{\&CenterDot;}{\mathrm{\&eta;}}-e_{\mathrm{\&eta;}}\&CenterDot;V|(\stackrel{\&CenterDot;}{\mathrm{\&eta;}}-e_{\mathrm{\&eta;}}\&CenterDot;V)+\mathrm{\&rho;}\&ForAll;{\stackrel{\&RightArrow;}{a}}_{\mathrm{\&eta;}}+C_{\mathrm{m\&eta;}}\mathrm{\&rho;}\&ForAll;({\stackrel{\&RightArrow;}{a}}_{\mathrm{\&eta;}}-e_{\mathrm{\&eta;}}\&CenterDot;\stackrel{\&CenterDot;}{V})---\left(5\right)$

$F_{\mathrm{\&xi;}}=-\frac{1}{2}\mathrm{\&rho;}{C}_{\mathrm{d\&xi;}}\mathrm{Dl}|\stackrel{\&CenterDot;}{\mathrm{\&xi;}}-e_{\mathrm{\&xi;}}\&CenterDot;V|(\stackrel{\&CenterDot;}{\mathrm{\&xi;}}-e_{\mathrm{\&xi;}}\&CenterDot;V)+\mathrm{\&rho;}\&ForAll;{\stackrel{\&RightArrow;}{a}}_{\mathrm{\&xi;}}+C_{\mathrm{m\&xi;}}\mathrm{\&rho;}\&ForAll;({\stackrel{\&RightArrow;}{a}}_{\mathrm{\&xi;}}-e_{\mathrm{\&xi;}}\&CenterDot;\stackrel{\&CenterDot;}{V})$

其中，F_τ、F_η和F_ξ表示作用于锚绳(1)单元的波浪力，C_Dτ、C_Dη和C_Dξ表示拖曳力系数，D是锚绳(1)直径，l是锚绳(1)单元的长度，

和表示质点的速度分量，和

表示质点的加速度分量；

其中，水动力系数采用下式计算：

$C_{\mathrm{Dn}}=\left\{\begin{array}{cc}0.0& ({\mathrm{Re}}_n\&le;0.1)\\ \frac{0.45+5.93}{{\left({\mathrm{Re}}_n\right)}^{0.33}}& (0.1<{\mathrm{Re}}_n\&le;400)\\ 1.27& (400<{\mathrm{Re}}_n\&le;{10}^5)\\ 0.3& ({\mathrm{Re}}_n>{10}^5)\end{array}\right.---\left(6\right)$

$C_{\mathrm{D\&tau;}}=\left\{\begin{array}{cc}\frac{1.88}{{\left({\mathrm{Re}}_n\right)}^{0.74}}& (0.1<{\mathrm{Re}}_{\mathrm{\&tau;}}\&le;100.55)\\ 0.062& ({\mathrm{Re}}_{\mathrm{\&tau;}}>100.55)\end{array}\right.---\left(7\right)$

其中，Re_n＝ρ|V_Rn|D/μ，μ是水的粘性系数；

计算了锚绳(1)单元上的外力之后，将外力均匀的分布到与之相邻的质点上，锚绳(1)质点的运动方程表示如下：

$m_i{a}_i=\underset{j=1}{\overset{\mathrm{count}}{\mathrm{\&Sigma;}}}({\text{T}}_j+W_j+B_j+F_j)---\left(8\right)$

其中，下标i表示质点编号，下标j表示与质点相邻的单元编号，count表示与质点相邻的单元数；

D、建立浮球(2)模型

浮球(2)漂浮在水面之上，受到波浪力的作用，作用于浮球(2)上的波浪力采用下式计算：

${\stackrel{\&RightArrow;}{F}}_B={\stackrel{\&RightArrow;}{F}}_D+{\stackrel{\&RightArrow;}{F}}_I=\frac{1}{2}\mathrm{\&rho;}{C}_{D}A\frac{{\stackrel{\&RightArrow;}{V}}_{\mathrm{RB}}\left|{\stackrel{\&RightArrow;}{V}}_{\mathrm{RB}}\right|}{2}+\mathrm{\&rho;}{\&ForAll;}_B{C}_m\frac{\&PartialD;{\stackrel{\&RightArrow;}{V}}_{\mathrm{RB}}}{\&PartialD;t}+\mathrm{\&rho;}{\&ForAll;}_B\frac{\&PartialD;{\stackrel{\&RightArrow;}{V}}_B}{\&PartialD;t}---\left(9\right)$

其中，C_D为拖曳力系数，C_m为浮架(3)质量系数，A是浮球(2)在水质点速度方向上的投影面积，

是浮球(2)的入水体积；拖曳力是雷诺数Re的函数，查表1可以得到拖曳力系数取值，表中未列出的数值采用线性插值得到：

表1浮球(2)的拖曳力系数

(x_B，y_B，z_B)为浮球(2)中心的坐标，浮球(2)入水深度为：

Δh＝η(x，y，t)-(z_B(t)-r)    (10)

其中，η是(x_B，y_B，z_B)处的波面，能够通过线性波浪理论获得，r是浮球(2)的半径；

浮球(2)的投影面积计算如下：

$A_x=A_y=\left\{\begin{array}{cc}\frac{\mathrm{\&pi;}}{4}{D}^2& (\mathrm{\&Delta;h}\&GreaterEqual;D)\\ \frac{1}{8}{D}^2(\mathrm{\&theta;}-\sin \mathrm{\&theta;})& (D/2<\mathrm{\&Delta;h}\&le;D)\\ \frac{\mathrm{\&pi;}}{4}{D}^2-\frac{1}{8}{D}^2(\mathrm{\&theta;}-\sin \mathrm{\&theta;})& (0<\mathrm{\&Delta;h}\&le;D/2)\\ 0& (\mathrm{\&Delta;h}<0)\end{array}\right.---\left(11\right)$

$A_z=\left\{\begin{array}{cc}\frac{\mathrm{\&pi;}}{4}{D}^2& (\mathrm{\&Delta;h}\&GreaterEqual;D/2)\\ \mathrm{\&pi;}({\left(\frac{D}{2}\right)}^2-{(\frac{D}{2}-\mathrm{\&Delta;h})}^2)& (0<\mathrm{\&Delta;h}\&le;D/2)\\ 0& (\mathrm{\&Delta;h}<0)\end{array}\right.---\left(12\right)$

其中，D是浮球(2)半径，θ是浮球(2)与水面相交的弦对应的中心角，通过下式可得

$\mathrm{\&theta;}=2{\cos }^{-1}\left(\frac{D/2-\mathrm{\&Delta;h}}{D/2}\right)---\left(13\right)$

E、运动微分方程求解

上述方程(2)、(3)、(4)和(8)构成了个网箱群组的运动微分方程组，采用Runge-Kutta-Vener六阶数值方法求解该运动微分方程组，从而得到整个网箱群组结构在每一时刻的运动响应以及锚绳(1)受力情况。

Images