基于CFD的刚性养殖网衣流场数值模拟及不确定度分析

2023-06-20张晓莹

渔业研究 2023年3期

张为，郭军*，扈喆，张晓莹，丁兰

(1. 集美大学轮机工程学院，福建厦门 361021；2. 福建省水产研究所，福建省海洋生物增养殖与高值化利用重点实验室，福建厦门361013)

渔业过度捕捞使渔业资源量下降，单纯依靠海洋捕捞早已无法满足人类对水产品日益增长的需求，现在迫切需要海水养殖来弥补水产品消费市场的不足。2021年，我国水产品总产量达6 690×104t，养殖产量与捕捞产量之比约为4∶1，网箱养殖为人们提供了丰富的海鲜食品。随着国内外各种大型养殖设施的建造，养殖业逐渐向深远海方向发展[1]，养殖的种类有鲍鱼[2]、大黄鱼[3]和云龙石斑鱼[4]等。在养殖过程中，网衣对海水有阻流作用，这不仅会影响养殖对象生长，还会阻碍水体交换，影响网箱内部的养殖环境，因此网衣内部的流场特性对网箱养殖非常重要。

研究网衣流场特性通常可以采用模型试验和数值模拟两种方法。虽然模型试验可以进行养殖网箱网衣流场的分析，但成本高，且受试验场地的限制。随着计算流体力学(Computational fluid dynamics，CFD)的不断发展，CFD技术已能较为准确地分析一些流体力学问题，不仅能得到养殖网箱整体受力情况，还可以捕捉养殖网箱内外流场的细节。桂福坤[5]对深水重力式网箱的水动力特性开展了模型试验，给出了重力式、碟形及拟碟形网箱的受力及运动特性；Patursson O等[6-7]对平面网衣的流场进行了模型试验和数值模拟，测量了网衣的阻力系数及升力系数，模拟了网衣后面的流场速度分布；赵云鹏等[8-9]借助商业软件FLUENT对单片网衣周围的流场特性进行了二维数值模拟，结果表明：FLUENT中的多孔介质模型可以用于网衣周围流场的数值模拟；刘兴[10]采用多孔介质模型，分别对平面网衣的流场进行了二维和三维数值模拟，结果表明：三维数值模拟的误差更小；刘春宏等[11]对网衣和鱼类共同作用下的网箱周围流场进行了数值模拟，发现需同时考虑鱼类和网衣才能准确地模拟养殖网箱周围的流场。

为研究及分析养殖网衣周围的流场特性，本文基于CFD理论，采用多孔介质模型对单片网衣的流场进行数值模拟，首先进行不确定度分析，再对网衣在不同流速下的流场进行数值计算，并与试验结果进行了比较，最后对网衣在不同攻角下的流场进行数值计算，分析了网衣的流场特性，给出网衣对流场内速度影响的量化区域及范围，旨在能够较好地为实际养殖生产活动提供参考依据。

1 理论模型

1.1 控制方程

流体力学问题要遵守质量守恒、动量守恒和能量守恒定律，对于不可压缩流体，其连续性方程和动量守恒方程分别如下[12]：

(1)

(2)

式(1)～(2)中：ui和uj为速度；xi和xj为坐标分量；t为时间；p为流体的压力；ρ为流体的密度；μ为流体的黏性系数。

控制方程不封闭，因而需引入湍流模型，使方程能够求解。本文计算采用标准k-ε湍流模型[12]，该模型具有较髙的稳定性、经济性和计算精度，应用广泛，其表达形式如下：

(3)

(4)

式(3)～(4)中：k和ε分别为湍动能和耗散率；μt为湍流黏性系数；Gk和Gb分别为平均速度梯度和浮力引起的湍动能产生项；σk和σε分别为k和ε的普朗特数。常数σk=1.0、σε=1.3、Cε1=1.44、Cε2=1.92。

1.2 多孔介质模型

将养殖网衣视为多孔透水板，采用多孔介质模型进行数值模拟。在数值模拟中，以动量方程的右边增加源项来模拟多孔介质，动量方程[13]变为：

(5)

式(5)中：Si为源项，其表达式为：

Si=-(Pνu+Piuu)

(6)

式(6)中：Pν和Pi分别为多孔介质黏性力系数矩阵和惯性力系数矩阵。

水流作用在多孔介质区域的合力(F)[13]为：

F=SitA

(7)

式(7)中：t为多孔介质的厚度；A为多孔介质的面积。

可将F分解为与水流方向相平行的阻力(D)和与水流方向相垂直的升力(L)，当水流方向与网衣平面垂直时(攻角α=90°，如图1所示)，D和L的表达式[13]分别为：

(8)

(9)

网衣密实度(S)的表达式[13]为：

(10)

式(10)中：d为网线直径；λ为网目长度。

网衣密实度与多孔介质孔隙率的关系为：

S+χ=1

(11)

式(11)中：χ为多孔介质的孔隙率。

2 数值计算方法

2.1 计算模型

本文研究对象为Patursson单片网衣试验过程中的平面网衣[6]，网衣尺寸的具体参数见表1。

表1 网衣的尺寸

网衣试验的水槽长37.00 m，宽3.66 m。网衣距离水槽两侧1.33 m，位于水槽液面0.73 m，以减少兴波的影响。网衣通过钢管固定在水槽中心，钢管直径为0.033 m，钢管可绕中心轴转动，从而改变网衣的攻角[6]。网衣布置图如图2所示。

2.2 网格生成

数值计算的计算域和边界条件如图3所示。计算域的长度为5 m，宽度和高度与试验水槽相同，多孔介质模型距离进流面入口1.5 m，多孔区域厚度为0.05 m。计算域的边界条件设置如下：进流面设置为速度入口条件；出流面为压力出口条件；水槽四周为无滑移壁面条件；水槽与多孔区域的交界面为内部界面条件。

数值计算采用的网格类型为非结构切割体网格。在网衣周围区域和网衣后侧区域进行了网格细化，用于捕捉网衣后侧的流场。计算域的总网格数为31×104，其中加密区域1和加密区域2的网格为18×104和8×104。整体计算域网格和中间水平截面网格如图4所示。

2.3 计算方法

在数值计算中，求解器采用基于分离流的黏性求解器，压力速度耦合算法为SIMPLE法，对流项采用二阶离散格式。湍流模型采用标准k-ε湍流模型，壁面函数采用高y+壁面处理。采用三维定常计算，最大迭代步数为1 000步。

对网衣试验结果进行拟合，多孔介质的法向惯性系数和切向惯性系数分别为2 492.5、830.0kg/m4，法向黏性系数和切向黏性系数分别为75.08、38.31kg/(m3·s)[9]，多孔介质的孔隙率为0.802。

3 数值计算结果

3.1 不确定度分析

网衣阻力系数和升力系数的定义[13]为：

(12)

(13)

式(12)～(13)中：Cd为网衣阻力系数；D为网衣所受的阻力；Cl为网衣升力系数；L为网衣所受的升力。

表2 网格的计算结果

可以看出，3套网格在数值迭代计算75步后，网衣阻力都已经完全收敛，将数值模拟的总迭代步数设置为75～100，可节省计算工作量。1号网格、2号网格和3号网格的网衣阻力系数分别为0.263、0.262和0.260。采用国际拖曳水池(ITTC)的不确定计算方法对网衣阻力进行不确定度分析，不确定度分析分为验证与确认两个过程[14]。

1)验证

数值不确定度由迭代不确定度(UI)、网格不确定度(UG)以及其他因素的不确定度(UP)组成。本文数值模拟采用定常计算，网衣阻力曲线几乎无波动，迭代误差可以忽略不计，验证过程主要为计算数值模拟中的网格不确定度UG。

网格的收敛半径RG为：

(14)

式(14)中：ε21为细网格与中网格的差值；ε32为中网格与粗网格的差值。

网格收敛半径为0.281，0

(15)

(16)

(17)

(18)

2)确认

确认是利用试验数据来评估数值模型不确定度USN的过程，将数值模拟值(S)与试验值(D)进行比较。

(19)

(20)

式(19)～(20)中：E为比较误差；UV为数值结果确认过程的不确定度；UD为试验结果的不确定度。

当E>UV，说明数值结果与试验值偏离较远，需要改进数值计算方法。

3.2 不同流速下的计算结果

在来流为0.125、0.250、0.500和0.750 m/s等4个速度下，对攻角为90°的网衣流场进行数值模拟，并与Patursson O的网衣阻力试验结果[6]进行比较，以验证数值模拟方法的准确性及有效性。

网衣阻力系数计算结果和试验结果对比如图6所示。由图可见，当来流为0.125 m/s时，网衣阻力的数值模拟结果偏大，误差为13.4%，主要原因为低速时网衣阻力的数值较小，相对误差较大；当来流分别为0.250、0.500和0.750 m/s时，网衣阻力的误差分别为-1.85%、1.75%和-1.11%。整体上，网衣阻力的数值模拟结果与试验结果吻合良好，说明本文数值模拟方法有效，可用于养殖网箱流场的数值预报。

在4个速度下，网衣中心截面的速度分布如图7所示。由图7可见，网衣边缘端部的流体速度大约会增加2%；由于网衣的阻挡作用，网衣前后的流体速度都会降低，后方的流体速度大约会降低10%；随着来流速度的增加，网衣两端流体速度增加的区域逐渐变小，网衣后方对来流速度的阻挡作用逐渐减弱。

在网衣中心法方向设置50个监测点以监测网衣前后流场的速度变化，在4个速度下，网衣中心前后的速度分布如图8所示。由图8可见，在网衣前后(网前0.5倍网衣长度到网后1.0倍网衣长度的范围内)，流体速度下降较快，网衣对流体有明显的阻挡作用；超过1.0倍网衣长度，流体速度基本不再变化。

3.3 不同攻角下的计算结果

在来流为0.5 m/s的速度下，对攻角分别为90°、60°、45°、30°和0°的网衣流场进行数值模拟，在不同攻角情况下，网衣阻力系数和升力系数的计算结果和试验结果对比分别如图9和图10所示。由图9可见，不同攻角下网衣阻力的数值模拟结果偏大；随着攻角角度的减小，网衣阻力系数不断减小，而误差逐渐增大。由图10可见，不同攻角下网衣升力的数值模拟与试验结果吻合较好；随攻角角度的增加，网衣升力系数先增大后减小；当攻角为45°时，网衣受到的升力最大。

在不同攻角下，网衣中心截面的速度分布如图11所示。由图11可见，当攻角变化后，网衣后面减流区域的宽度也发生变化，宽度基本上为网衣的投影宽度；网衣上端部的流体速度大于下端部的流体速度，网衣下部对流体的阻挡作用大于上部。

在不同攻角下，网衣中心前后的速度分布如图12所示。由图12可见，当攻角在30°至90°变化时，网衣中心前后速度变化的规律基本相同。当攻角为0°(来流与网衣平行)时，受网衣的阻挡，在网后1.0倍网衣长度范围内，流体速度小于其他攻角的速度；在网后大于1.0倍网衣长度时，发生绕射现象，最终流体速度反而大于其他攻角的速度。