CN111737897B - 一种深海网箱高密度养殖鱼群的数值模拟方法 - Google Patents

Abstract

本发明公开了一种深海网箱高密度养殖鱼群的数值模拟方法，该方法采用刚体鱼模型对养殖鱼群进行模拟，刚体鱼包括鱼身和鱼尾部分，在鱼身和鱼尾分别施加拖曳力和升力，模拟高密度鱼群对网箱周围流场和变形的影响。采用k‑ωSST湍流模型模拟网箱及高密度养殖鱼群周围的流场，采用非线性结构有限元模型计算网箱的结构变形，最终实现了高密度养殖鱼群对网箱周围流场和变形影响的数值模拟。采用本方法可以计算得到不同的养殖鱼群分布方式、养殖密度、游动速度和加速度条件下网箱周围流场和结构变形。该模型解决了深海网箱内部高密度养殖鱼群的数值模拟的问题，可以很好地用于实际网箱的动力响应分析，为深海网箱的结构设计提供支撑。

Description

一种深海网箱高密度养殖鱼群的数值模拟方法

技术领域

本发明属于海洋工程和水产养殖工程技术领域，特别涉及到一种波流作用下网箱养殖设施内部高密度养殖鱼群对水流影响的三维数值模拟方法。

背景技术

随着海洋养殖从近岸走向深水的过程中，大型网箱是未来海洋养殖的发展趋势。利用网箱进行水产养殖是一种高密度、集约化的水产养殖方式，网箱内部高密度的养殖鱼群的活动将对网箱周围的流场产生显著的影响，势必也将影响到网箱的荷载以及结构变形；利用数值模拟的方法研究鱼群对于网箱周围流场及其所受荷载的影响机制是建立网箱养殖设施结构动力响应分析的精细化模型关键。

目前，国内外关于波流作用下网箱动力响应的数值模拟研究，并没有考虑高密度养殖鱼群对网箱内部流场及其变形的影响，这将影响到网箱内部流场和结构变形的准确计算。传统的关于鱼的游动的数值模型，是采用计算流体动力学的方法对鱼的游动进行模拟，但是这一类方法在进行高密度养殖鱼群的模拟时，需要大量的细致的网格，会严重影响计算效率，使得计算无法进行。因此，提出一种行之有效的高密度养殖鱼群的数值模拟方法是十分必要的，对于深水网箱的设计具有重要意义。

发明内容

为解决现有技术存在的上述问题，本发明的目的是提供一种波流作用下养殖网箱内部高密度鱼群的数值模拟技术。

本发明的技术方案：

一种深海网箱高密度养殖鱼群的数值模拟方法，步骤如下：

通过流场模型计算获得网衣周围的流场和压力场，然后对网衣表面的压力进行积分获得网衣受到的水动力荷载，进一步通过网衣的非线性有限元模型计算得到网衣的应力和变形；具体方法如下：

A、流场计算模型

将网衣简化成一系列的圆柱结构，模拟网箱对周围流场的影响；采用有限体积法求解网箱周围流场的控制方程，计算获得网箱周围的流场；

B、网衣结构有限元模型

波流作用下柔性网衣将产生大变形，采用大变形非线性结构模型建立网衣的有限元模型，计算网衣的应力、应变和位移；网衣简化成一系列杆单元和与之相连的球形铰单元；杆单元和与之相连的球形铰单元之间无相对平动位移，可产生相对转动；

C、高密度养殖鱼群的模拟

采用刚体鱼模型对网箱中的养殖鱼群进行模拟，分析高密度养殖鱼群对网箱周围流场及其变形的影响；刚体鱼模型包括两部分：鱼身和鱼尾，忽略鱼鳍和尾鳍；鱼身部分受到水流产生的拖曳力(阻力)，而鱼尾部分产生推进力。

本发明的有益效果：

(1)本发明采用刚体鱼模型对高密度养殖鱼群进行模拟，分别在鱼身和鱼尾施加拖曳力和推进力，可以对高密度养殖鱼群进行模拟，高效地实现养殖鱼群对网箱周围流场和变形影响的数值模拟，为深海网箱的设计奠定重要基础。

(2)本发明提出的刚体鱼模拟，可以很好地模拟养殖鱼群的分布方式、游动速度和加速度等参数对网箱流场和变形的影响，具有很好的实用性。

(3)传统方法对鱼的运动的模拟采用计算流体动力学方法实现，需要大量的计算网格，计算效率低，无法准确控制养殖鱼群的分布方式和游动速度，本发明提出的养殖鱼群的数值模拟方法，能够减少计算网格数量，实现养殖鱼群分布方式、游动速度等模型输入参数的精确控制。

(4)相对于物理模型试验而言，该发明可以直接对原型网箱内部的养殖鱼群进行模拟，解决了模型试验无法满足的鱼群相似性的问题。

附图说明

图1是无节网结构的几何模型，图1a是网衣结构，图1b是网箱结构。

图2是数值模型计算中的流体网格。

图3是刚体鱼周围的网格。

图4是网线表面的结构单元。

图5是网衣结构的杆单元和球形铰单元，F_wi＝F_hi+F_bi，其中：F_hi是单元i水动力，F_bi是单元i浮力，F_g是单元i的自重，C_i是单元i的几何中心。

图6是网衣结构表面的流体单元节点和结构单元节点的映射关系示意图。

图7是网衣结构水动力荷载传递示意图。

图8是刚体鱼模型示意图，其中V_F是水流流速矢量，V_a是鱼游泳速度矢量，θ是鱼身与流向夹角，F_D和F_L分别是鱼所受拖曳力和升力。

图9是刚体鱼模型验证计算布置图，其中V_C是鱼在水流中的相对速度矢量。

图10是不同运动速度条件下刚体鱼模型和经验公式计算的拖曳力系数对比。

图11是不同回转角条件下刚体鱼模型和经验公式计算的拖曳力系数对比。

图12是不同回转角条件下刚体鱼模型和经验公式计算的升力系数的对比。

图13是数值模拟和模型试验的网衣变形对比。

图14是数值模拟和模型试验网衣拖曳力对比。

图15是数值模拟和模型试验的网箱变形对比。

图16是不同流速条件下网箱内部固定测点位置处的数值模拟和试验测得的流速对比。

图17是不同位置处数值模拟与试验测得的流速对比。

图18是不同流速条件下数值模拟和试验测得的网箱拖曳力对比。

图19是夜间分布方式时网箱内部的流线图。

图20是白天分布方式时网箱内部的流线图。

图21是不同的网箱养殖密度时网箱内部流速分布。

图22是不同的网箱养殖密度时网箱的变形。

具体实施方式

下面结合附图对本发明进行进一步地描述。如图1-22所示，采用本发明模拟养殖鱼群对网箱周围流场和变形的影响，具体过程如下：

一种深海网箱高密度养殖鱼群的数值模拟方法，步骤如下：

通过流场模型计算获得网衣周围的流场和压力场，然后对网衣表面的压力进行积分获得网衣受到的水动力荷载，进一步通过网衣的非线性有限元模型计算得到网衣的应力和变形；具体方法如下：

A、流场计算模型

将网衣简化成一系列的圆柱结构，模拟网箱对周围流场的影响；采用有限体积法求解网箱周围流场的控制方程，计算获得网箱周围的流场；

采用k—ωSST湍流模型模拟网箱周围的流场特性，具体控制方程如下：

连续性方程：

式中，ρ是水的密度，t表示时间，x，y，z表示坐标，u，v，w表示速度；

雷诺平均Navier-Stokes方程：

式中，μ是流体粘性系数；P是压力；

和

是沿三个坐标轴的时均速度分量，u^′，v^′和w^′是沿三个坐标轴方向的速度脉动分量；

湍动能方程：

式中，x_i为i方向坐标分量，U_i为i方向速度分量，k为湍流动能，σ_k为湍流动能普朗特数，μ_t为湍流粘度，

x_j为j方向坐标分量，U_j为j方向速度分量，ω是湍动能耗散率；

湍动能耗散率方程：

式中，α为旋度生成系数，β为湍流耗散系数，v_t是湍流粘性系数，σ_ω1和σ_ω2是常数，混合函数F₁如下所示：

其中，

y′是距壁面的距离，v为流体运动粘性系数；

湍流粘性系数如下：

其中，S为应变率张量，混合函数F₂由下式确定：

k—ωSST湍流模型常数如下：α₁＝5/9，α₂＝0.44，β₁＝3/40，β₂＝0.0828，β^*＝9/100，σ_k1＝0.85，σ_k2＝1，σ_ω1＝0.5，σ_ω2＝0.856，a₁＝0.31；

B、网衣结构有限元模型

波流作用下柔性网衣将产生大变形，采用大变形非线性结构模型建立网衣的有限元模型，计算网衣的应力、应变和位移；网衣简化成一系列杆单元和与之相连的球形铰单元；杆单元和与之相连的球形铰单元之间无相对平动位移，可产生相对转动；球形铰单元和杆单元之间的约束关系如下：

其中，局部坐标系C₀-x₀-y₀-z₀的原点位于球形铰的中心C₀，

是局部坐标系沿各坐标轴的单位向量，E_x,E_y,E_z是整体坐标系沿各坐标轴的单位向量。

网衣非线性结构模型控制方程如下：

[K(x)]x＝F_h+F_b+F_g＝Q (11)

其中，F_h是水动力，F_b是浮力，F_g是重力，x是网衣的位移向量，K(x)是网衣非线性刚度矩阵；基于虚功原理，网衣结构的有限元模型控制方程如下：

其中，δ是节点位移增量向量，^t+ΔtQ是节点荷载向量，

是等效单元应力的节点力向量，

就是小位移情况下的单元刚度矩阵，

是单元初位移刚度矩阵，

是单元的初应力刚度矩阵；

和

表示为：

其中，₀D是材料本构矩阵，

和

是第二类Piola-Kirchhoff应力矩阵和向量，

和

分别是线性应变和非线性应变和位移的转换矩阵：

其中，₀N_k，j为k节点在j方向上初始位形的差值偏导，j＝1,2,3；N_k为节点k关联的差值函数，⁰x_j为单元节点初始位形的j方向坐标分量，j＝1,2,3；

为k节点初始位形的差值函数，

为k节点在j＝1方向上的总体坐标分量，

其中，₀N_k,1和₀N_k,2分别为k节点在j＝1和j＝2方向上参考初始位形的差值偏导，

和

分别为节点k在时间t位形的i＝1和i＝2方向上的位移分量，

其中，

为t时刻参考初始位形的单元应变张量，

为t时刻参考初始位形的单元材料本构张量具体表达为：

其中，G为材料的剪切模量，γ为材料泊松比，δ_ik、δ_jl、δ_ij和δ_kl均为克罗内克函数(即

)，Q为结构受到的外载荷，包括重力载荷Q₁、浮力Q₂和水动力Q₃；Q＝Q₁+Q₂+Q₃；重力Q₁＝ρ₁gV，浮力Q₂＝ρgV，水动力力Q₃＝PA，其中，ρ₁是网衣的网线密度，ρ是水的密度，g是重力加速度，V是网衣结构单元的体积，P是作用于网线外表面的流体压力，由k-ωSST湍流模型计算得到，A是网线单元的表面积；

结构计算时，同样采用有限体积法离散控制方程，离散格式采用Newton-Raphson迭代算法，平衡方程按二阶格式进行计算，

其中，

^t+Δtδ^(l+1)＝^t+Δtδ^(l)+Δδ^(l)

其中，l为迭代次数，^t+Δtδ^(l+1)为t+Δt位形下第l+1次迭代的节点总体位移增量，^t+Δtδ^(l)为t+Δt位形下第l次迭代的节点总体位移增量，Δδ^(l)是第l次迭代的节点位移增量，^t+ΔtQ是t+Δt位形下的节点荷载，

是t+Δt位形下第l次迭代的等效节点力，

是小位移情况下的单元刚度矩阵，

是单元初位移刚度矩阵，

是单元的初应力刚度矩阵；

C、高密度养殖鱼群的模拟

采用刚体鱼模型对网箱中的养殖鱼群进行模拟，分析高密度养殖鱼群对网箱周围流场及其变形的影响；刚体鱼模型包括两部分：鱼身和鱼尾，忽略鱼鳍和尾鳍；鱼身部分受到水流产生的拖曳力(阻力)，而鱼尾部分产生推进力；

刚体鱼模型的合力包括鱼身受到的拖曳力和鱼尾产生的推进力，作用于刚体鱼模型上的合力沿x和y坐标轴的分量如下所示：

其中，F_D、F_L和F_h分别是鱼身的拖曳力、升力和鱼尾的推进力，τ_d和τ_h分别是鱼身和鱼尾的壁面剪切应力；鱼身受到的拖曳力和升力与鱼身的形状阻力和粘性阻力有关，鱼身的拖曳力和升力采用下式计算：

C_L＝C′_Lβα₀＝0.935α₀ (22)

其中，S_fish是鱼身的表面积；V_c是鱼的相对运动速度；θ是鱼身的冲击角；C_D是拖曳力系数，取为0.1936；C_L是升力系数，取为0.1412 rad^-2；α₀是鱼尾的摆动角；鱼身冲击角和鱼尾摆动角之间的关系为θ＝0.2315α₀。

图1是无结节网衣的几何模型，网目形状为方形，目脚尺寸为20mm，网线直径为2.6mm，网线材料为聚乙烯，网箱结构的浮架系统为高强度聚乙烯材料。

图2是网衣的数值模型计算的网格布置，为了保证计算结果的可靠性，在网衣周围对网格进行了加密。

图3是刚体鱼周围的网格设置，在刚体鱼附近对网格进行了加密处理。

图4是网线结构表面的网格设置，在网线结构的表面采用O形网格。

图5是网线的结构模型设置，网线由杆单元和与之相连的球单元相连，杆单元与球单元可以自由转动，杆单元与球单元之间不传递弯矩荷载。

图6是网衣的流场模型与结构模型的节点之间的数据传递关系，通过在流场单元上计算获得的压力进行插值，得到结构单元上的流体压力值，之后通过积分，计算得到网线的水动力荷载。

图7是网衣结构的水动力荷载，主要包括杆单元和球单元受到的水动力荷载。

图8是刚体鱼的数值模型示意图，刚体鱼数值模型包括两部分，鱼身和鱼尾，分别在鱼身和鱼尾施加拖曳力和推进力，模拟刚体鱼对周围水流的作用。

图9是对刚体鱼数值模型的验证，考虑刚体鱼不同游动速度V_c、不同艏摇角ε的计算工况。

图10是数值模拟和经验公式计算得到的不同游动速度条件下刚体鱼的拖曳力系数，结果表明不同的游动速度条件下，刚体鱼的拖曳力系数数值模拟结果与经验公式计算结果吻合良好。

图11是数值模拟和经验公式计算得到的不同艏摇角条件下刚体鱼的拖曳力系数，结果表明不同的艏摇角条件下，刚体鱼拖曳力系数的数值模拟结果与经验公式计算结果吻合良好。

图12是数值模拟和经验公式计算得到的不同艏摇角条件下刚体鱼的升力系数，结果表明不同的艏摇角条件下，刚体鱼升力系数的数值模拟结果与经验公式计算结果吻合良好。

图13给出了数值模拟和模型试验得到的网衣的变形，结果表明该模型可以很好地模拟网衣的变形。

图14表明该数值模型可以很好地模拟网衣受到的拖曳力。

图15表明该数值模型可以很好地模拟不同流速条件下网箱的变形。

图16-17给出了不同流速条件下、不同位置处的网箱周围流速的对比，结果表明该数值模型可以很好地模拟网箱周围的流场。

图18给出了不同流速条件下网箱的拖曳力的对比，表明该数值模型可以很好地模拟网箱受到的拖曳力荷载。

图19-20给出了鱼群的夜间分布和白天分布时网箱周围的流线图，表明该模型可以很好地分析鱼群的运动引起的网箱内部的水流流动。

图21给出了不同的养殖密度条件下网箱内部的流速分布，结果表明该数值模型可以很好地模拟养殖鱼群对网箱内部流场的影响。

图22给出了不同养殖密度条件下网箱的变形，结果表明该数值模型可以很好地模拟养殖鱼群对网箱变形的影响。

Claims (1)

1.一种深海网箱高密度养殖鱼群的数值模拟方法，其特征在于，步骤如下：

通过流场模型计算获得网衣周围的流场和压力场，然后对网衣表面的压力进行积分获得网衣受到的水动力荷载，进一步通过网衣的非线性有限元模型计算得到网衣的应力和变形；

具体方法如下：

A、流场计算模型

将网衣简化成一系列的圆柱结构，模拟网箱对周围流场的影响；采用有限体积法求解网箱周围流场的控制方程，计算获得网箱周围的流场；

采用k—ωSST湍流模型模拟网箱周围的流场特性，具体控制方程如下：

连续性方程：

式中，ρ是水的密度，t表示时间，x，y，z表示坐标，u，v，w表示速度；

雷诺平均Navier-Stokes方程：

式中，μ是流体粘性系数；P是压力；

和

是沿三个坐标轴的时均速度分量，u′，v′和w′是沿三个坐标轴方向的速度脉动分量；

湍动能方程：

式中，x_i为i方向坐标分量，U_i为i方向速度分量，k为湍流动能，σ_k为湍流动能普朗特数，μ_t为湍流粘度，

x_j为j方向坐标分量，U_j为j方向速度分量，ω是湍动能耗散率；

湍动能耗散率方程：

式中，α为旋度生成系数，β为湍流耗散系数，υ_t是湍流粘性系数，σ_ω1和σ_ω2是常数，混合函数F₁如下所示：

其中，

y′是距壁面的距离，υ为流体运动粘性系数；

湍流粘性系数如下：

其中，S为应变率张量，混合函数F₂由下式确定：

k—ωSST湍流模型常数如下：α₁＝5/9，α₂＝0.44，β₁＝3/40，β₂＝0.0828，β^*＝9/100，σ_k1＝0.85，σ_k2＝1，σ_ω1＝0.5，σ_ω2＝0.856，a₁＝0.31；

B、网衣结构有限元模型

波流作用下柔性网衣将产生大变形，采用大变形非线性结构模型建立网衣的有限元模型，计算网衣的应力、应变和位移；网衣简化成一系列杆单元和与之相连的球形铰单元；杆单元和与之相连的球形铰单元之间无相对平动位移，可产生相对转动；球形铰单元和杆单元之间的约束关系如下：

其中，局部坐标系C₀-x₀-y₀-z₀的原点位于球形铰的中心C₀，

是局部坐标系沿各坐标轴的单位向量，E_x,E_y,E_z是整体坐标系沿各坐标轴的单位向量；

网衣非线性结构模型控制方程如下：

[K(x)]x＝F_h+F_b+F_g＝Q (11)

其中，F_h是水动力，F_b是浮力，F_g是重力，x是网衣的位移向量，K(x)是网衣非线性刚度矩阵；基于虚功原理，网衣结构的有限元模型控制方程如下：

其中，δ是节点位移增量向量，^t+ΔtQ是节点荷载向量，

是等效单元应力的节点力向量，

就是小位移情况下的单元刚度矩阵，

是单元初位移刚度矩阵，

是单元的初应力刚度矩阵；

和

表示为：

其中，₀D是材料本构矩阵，

和

是第二类Piola-Kirchhoff应力矩阵和向量，

和

分别是线性应变和非线性应变和位移的转换矩阵：

其中，₀N_k,j为k节点在j方向上初始位形的差值偏导，j＝1,2,3；N_k为节点k关联的差值函数，⁰x_j为单元节点初始位形的j方向坐标分量，j＝1,2,3；

为k节点初始位形的差值函数，

为k节点在j＝1方向上的总体坐标分量，

其中，₀N_k,1和₀N_k,2分别为k节点在j＝1和j＝2方向上参考初始位形的差值偏导，

和

分别为节点k在时间t位形的i＝1和i＝2方向上的位移分量，

其中，

为t时刻参考初始位形的单元应变张量，

为t时刻参考初始位形的单元材料本构张量具体表达为：

其中，G为材料的剪切模量，γ为材料泊松比，δ_ik、δ_jl、δ_ij和δ_kl均为克罗内克函数(即

)，Q为结构受到的外载荷，包括重力载荷Q₁、浮力Q₂和水动力Q₃；Q＝Q₁+Q₂+Q₃；重力Q₁＝ρ₁gV，浮力Q₂＝ρgV，水动力力Q₃＝PA，其中，ρ₁是网衣的网线密度，ρ是水的密度，g是重力加速度，V是网衣结构单元的体积，P是作用于网线外表面的流体压力，由k-ωSST湍流模型计算得到，A是网线单元的表面积；

结构计算时，同样采用有限体积法离散控制方程，离散格式采用Newton-Raphson迭代算法，平衡方程按二阶格式进行计算，

其中，

^t+Δtδ^(l+1)＝^t+Δtδ^(l)+Δδ^(l)

其中，l为迭代次数，^t+Δtδ^(l+1)为tΔt位形下第l+1次迭代的节点总体位移增量，^t+Δtδ^(l)为tΔt位形下第l次迭代的节点总体位移增量，Δδ^(l)是第l次迭代的节点位移增量，^t+ΔtQ是t+Δt位形下的节点荷载，

是tΔt位形下第l次迭代的等效节点力，

是小位移情况下的单元刚度矩阵，

是单元初位移刚度矩阵，

是单元的初应力刚度矩阵；

C、高密度养殖鱼群的模拟

采用刚体鱼模型对网箱中的养殖鱼群进行模拟，分析高密度养殖鱼群对网箱周围流场及其变形的影响；刚体鱼模型包括两部分：鱼身和鱼尾，忽略鱼鳍和尾鳍；鱼身部分受到水流产生的拖曳力，而鱼尾部分产生推进力；

刚体鱼模型的合力包括鱼身受到的拖曳力和鱼尾产生的推进力，作用于刚体鱼模型上的合力沿x和y坐标轴的分量如下所示：

其中，F_D、F_L和F_h分别是鱼身的拖曳力、升力和鱼尾的推进力，τ_d和τ_h分别是鱼身和鱼尾的壁面剪切应力；鱼身受到的拖曳力和升力与鱼身的形状阻力和粘性阻力有关，鱼身的拖曳力和升力采用下式计算：

C_L＝C′_Lβα₀＝0.935α₀ (22)

其中，S_fish是鱼身的表面积；V_c是鱼的相对运动速度；θ是鱼身的冲击角；C_D是拖曳力系数，取为0.1936；C_L是升力系数，取为0.1412rad^-2；α₀是鱼尾的摆动角；鱼身冲击角和鱼尾摆动角之间的关系为θ＝0.2315α₀。

Images