基于桁架模型和多孔介质模型的柔性网衣和周围流场单向耦合方法研究

徐子鸣，林志良, 2，谢 彬, 2

(1. 上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海 200240; 2. 上海交通大学 船海工程数值试验中心，上海 200240)

水产养殖一直是世界粮食生产的重要方式，在全球范围内处于稳步扩张阶段。由于人口密度增加、过度捕捞和不合理开发等因素影响，近海水产资源逐渐枯竭。为了寻求更优越的养殖环境，把网箱移至外海的深海养殖应运而生[1-4]。在深水养殖过程中，网衣的减流效应备受关注，不仅决定了整个网箱的水动力学特性，而且对网箱内水质环境的影响也非常显著，因此需要更准确的水动力分析和计算方法。

网衣属于小尺度、大变形的柔性多孔结构体，周围的流场分布极其复杂，许多学者对此提出了不同研究方法。Lφland[5]通过Rudi等[6]的试验数据拟合出网衣阻力和升力系数与来流速度的关系；宋伟华等[7]通过模型试验研究了方形网箱外部的流速分布；Tsukrov等[8]和Tang等[9]通过大量试验发现，网衣水动力特性主要依赖于雷诺数Re和坚固性两个无量纲变量。网衣的数值模拟难度较大，通常将流场导致的网衣变形和网衣对流场的作用分别进行求解。其中网衣变形模型一般使用经验公式，经验公式有两种，一种是Morison模型[10-11]，把每一根网线当作独立的圆柱并将其作用力通过求解Morison方程计算；第二种为Screen模型[12]，将网衣划分为多个平面网板单元并根据密实度和攻角计算网衣的阻力和升力系数。雷江涛等[13]和Cheng等[14]分别通过Morison模型和Screen模型对网衣在定流速情况下的变形和受力情况进行分析。而网衣周围的流场分布可采用计算流体力学(CFD)方法进行求解。这种方法将网衣表示成多孔介质[15]并通过在Navier-Stokes方程中添加源项来表示网衣对流体的影响，其特点是不需要对网衣的详细几何形状进行建模，因此具有计算简便和效率较高的优点。赵云鹏等[16]结合多孔介质模型模拟网衣，探究了在不同流速、不同攻角、不同网衣片数等工况下水流作用下平面网衣周围流场的情况。在网衣和流场的耦合模拟方面，Bi等[17-18]用软件Fluent将多孔介质模型和集中质量模型耦合，分别模拟了平面网衣和圆形网衣在定常流速下的流场情况，和试验对比了网衣受力。Chen和Christensen[19-20]通过在软件OpenFOAM框架下结合多孔介质模型和集中质量模型，开发出一种新的耦合接口用以流体和结构求解器间的数据交换，模拟结果与试验数据基本一致。但这种耦合方式需要在每个时间步内迭代至稳态，计算量很大，在实际应用中存在一定的局限性。Cheng等[21]利用code_aster和OpenFOAM进行网衣的流固耦合分析。新的耦合算法通过消除多孔介质系数的数据拟合过程并通过Screen模型计算网衣的水动力，从而达到简化计算过程提高模拟精度的效果。

下文提出一种更简便的计算方法，通过Screen模型求解网衣水动力，然后运用桁架模型建立其网衣节点位移的偏微分方程，在code_aster中计算网衣节点的位移值。当解得网衣变形后各节点的坐标后，再将其导入OpenFOAM中生成多孔介质区域并通过求解流体的控制方程来模拟流场情况。计算结果表明，网衣变形、流速分布和网阻力等特征参数与试验数据对比良好。研究结果为养殖网箱周围的流场分析提供了参考。

1 数学模型

1.1 桁架模型

1.1.1 概述

柔性网衣在桁架模型中被视为一组杆件和网格点。对于三维的柔性网，如图1(a)所示，单根杆件包含两个节点，并且每个节点都能在三维空间里任意运动。对于变形的柔性网，图1(a)展示柔性网的配置情况，作用在杆件的力被均匀地分布到其所关联的节点中，并且在每个节点上列出运动偏微分方程。这里结构求解器将根据Cheng等[14]提出的原理进行编写。

1.1.2 节点受力

作用在每个节点上的力包括了水动力Fh、相邻的两个节点的张力Fs、节点重力Fg和浮力Fb。Fb和Fg是恒定的，由网衣的材料所决定。在求解中，水动力采用Screen模型公式计算，如图1(b)所示，点A、B、C、D组成三维空间任意的四边形，一般来说A、B、C、D不在同一个平面上，为了保证节点受力容易计算，将四边形ABCD拆分为四个小三角形(ΔABC，ΔABD，ΔACD，ΔBCD)，以任意的一个小三角形ABC为例，水动力Fh=FD+FL，FD和FL的计算公式为：

图1 桁架模型Fig. 1 Sketch of the truss model

(1)

(2)

其中，At是ΔABC网片的面积，ur为网片与来流的相对速度，iD,iL是阻力、升力的单位力方向矢量，分别与相对速度方向相同和垂直，对于节点A的水动力，为Fh/3。在式(1)、式(2)中的CD,CL由式(3)～(5)计算：

CD=0.04+(-0.04+Sn-1.24Sn2+13.7Sn3)cosθ

(3)

CL=(0.57Sn-3.54Sn2+10.1Sn3)sin2θ

(4)

(5)

式中：Sn为网衣密实度，dw为网衣直径，L为目脚长度(如图2所示)，θ为网片相对来流的迎角。

图2 网衣模型Fig. 2 Net model

节点上的重力和浮力如图1(c)所示，将杆的重力和浮力平均分给每个节点，以杆件AB为例，A节点的重力和浮力为:

(6)

式中：dw为网衣直径，L为目脚长度，ρ和ρnet分别为流体密度和网的密度。

1.1.3 尾流效应

尾流效应是分析绕流的重要机制。在水动力模型Screen模型中，后网的来流速度衰减，因此尾流效应必须考虑。一般来说，尾流效应包含三种，如图3所示，网线之间、单个网衣之间、网笼之间，网线之间的尾流效应已经包含在Screen模型之中[14], 对于单片网衣之间，若前网的速度为u，后网速度为ru，r为衰减因子，Lφland[5]给出的经验公式为r=1-0.46CD,θ=0，CD,θ=0为迎角为0的阻力系数，网笼之间的尾流效应一般与单片网衣的经验公式相同。

图3 尾流效应Fig. 3 Wave effect

1.1.4 运动方程

对网衣节点位移建立偏微分方程，与Cheng等[14]研究有一样的方程为：

(7)

式中：M为质量矩阵，K为刚度矩阵，Fh、Fg和Fb分别为节点的水动力、重力和浮力，已经在1.1.2节中叙述，x为瞬态的节点位移。图4为由网衣节点1和2组成的杆件，在局部坐标系OS中可得：

图4 单杆单元Fig. 4 A rod element

(8)

由于节点1和节点2在全局坐标系oxyz都是任意的运动，都具有三个自由度，转换矩阵为：

(9)

式中：cos(x,s)，cos(y,s)，cos(z,s)分别为os与ox、oy、oz夹角的余弦值。因此，在全局坐标系oxyz下，M=M1·T,K=K1·T。

同时，对于单个节点的位移和速度推进运用HHT-α方法[14]:

(10)

(11)

(12)

其中，取α=0.3，Δt=0.01 s，在code_aster中模拟时间为10 s。

1.2 多孔介质模型

1.2.1 控制方程

Realizable k-ε模型是一种具有广泛适用性的湍流模型，适合解决涉及强压力梯度的流动问题。对于纯水流绕流于平面网，由于网的屏蔽作用，存在很强的压力梯度。因此，选择Realizable k-ε湍流模型来解决这个问题。控制方程包括不可压缩流体的连续性方程和动量守恒方程：

∇·u=0

(13)

(14)

其中，t是时间，p是压力，u是流体速度,μ和ρ分别是流体的动力黏度和密度，g是重力场。S是动量方程的源项。在多孔介质理论中，

(15)

其中，C是多孔介质系数矩阵。在局部坐标系中可表示为：

(16)

其中，Cn是法向惯性阻力系数，Ct是切向惯性阻力系数，对于矩阵C来说，可以由试验数据通过最小二乘法[16]得到，也可以通过Chen等[19]推导的公式获得，如果局部坐标系不与全局坐标系对齐，则需要对系数矩阵进行转换[19]。

Realizable k-ε模型方程为：

(17)

(18)

1.2.2 多孔介质区域的生成

根据Chen和Christensen[20]的方法，对于每个小的网片，如图1(b)所示，是由四个节点组成。但是四个节点不一定在一个平面上，因此将方形网格拆分成两个三角形面板，相应的多孔介质区域是由质量点组成三角形延伸到厚度t0的棱柱来生成的。

图5描述了如何通过全局网格生成多孔介质区域，图中小正方体代表一个网格，x1,x2,x3分别代表一个三角形面板的坐标，x代表全局网格体心的坐标，d代表全局网格体心到三角形面板的距离，n是三角形面板的单位法向量，x0为全局网格体心在三角形面板投影的坐标，其中d和x0的计算为：

图5 三维情况下流体模型与结构模型的映射图解Fig. 5 Illustration on mapping of the element between the fluid and structure model in 3D case

(19)

d=(x-x2)·n

(20)

x0=x-d·n

(21)

如果网格坐标点满足下列两个条件：1)d<1/2t0；2)由点x0x1x2,x0x2x3,x0x3x1组成的三个小三角形面积S1,S2,S3之和等于x1x2x3组成的三角形面积S，即投影点x0应该在x1x2x3组成的三角形内部，则将其标记为多孔区域的网格。多孔区域由所有满足判定条件的网格组成。

2 数值模型验证

2.1 圆形网衣模型验证

采用Bi等[18]所得的试验结果进行模型验证，试验在波浪水槽中进行，测试网为圆形网，在周向有40目，在垂直方向有8目，网衣直径为dw=1.2 mm，目脚长度为λ=20 mm，圆形框架底部直径为D=0.254 m，网衣高度为H=0.15 mm，框架为m=8 g的配重，在本算例中，8 g的配重除去浮力的部分将会平均分给40个底部的网衣节点，每个节点的受力为0.001 73 N。物理模型和测量点如图6所示，圆形网衣需要考虑容积损失，通过Zhao等[22]的方法来计算出容积损失率。由于多孔介质的厚度影响不大[15]，选择多孔介质为20 mm的厚度，多孔介质系数Cn=9.16 m-1,Ct=5.67 m-1。

图6 圆形网衣的数值模型和测量点的设置Fig. 6 Physical model of a circular net and general setting of the measurement points

2.2 网格划分和计算条件

数值模型的网格划分如图7所示。x为来流方向，y为水平面上与x垂直的方向，z为水平面上法向方向。水槽左端定义为速度入口边界；右端定义为出流边界；水槽侧壁和水面定义为滑移边界；水槽底面定义为固壁边界。网格类型为非结构化六面体网格，在多孔介质区域向外加密两层网格。采用有限体积法(FVM)离散控制方程，时间项目采用一阶欧拉隐式格式，对流项采用二阶迎风差分格式，压力—速度耦合选择瞬态的PISO算法，计算时间步长取0.005 s，求解时长为40 s。

图7 计算网格示例Fig. 7 Example of computational grids

3 结果与讨论

3.1 收敛性研究

为了保证在OpenFOAM计算准确，需要进行网格收敛性研究，在流速U=0.242 m/s的情况下，使用三种网格来计算网衣的阻力。图8是在运用1.2.2节所述的多孔介质生成方法，在OpenFOAM得到的多孔区域，从多孔区域看出，其中有部分网格缺失，原因是判断条件算法的局限性和静态网格的使用，多孔区域只能近似变形网的几何形状。表1是在三种不同网格下的计算结果，计算结果表明，在不同的全局网格数和多孔介质区域网格数下，计算流体力学(CFD)模拟的网衣阻力随网格数变化很小，与试验值基本吻合，网格收敛性良好，图9为在三种网格下的切面y=-0.05 m的速度切面云图，来流通过前网和后网都有明显的速度降低，两侧的绕流增加，有明显的速度过渡区域。

图8 OpenFOAM与试验结果之间净变形比较Fig. 8 Comparison of the net deformation between OpenFOAM and the experimental result

表1 OpenFOAM中网格收敛性研究Tab. 1 Mesh convergence study in OpenFOAM

图9 三种网格下切面y=-0.05 m的速度云图Fig. 9 Velocity distribution of three kinds of meshes under the slice y=-0.05 m

3.2 不同来流速度通过单片网衣计算结果

3.2.1 在code_aster中生成圆形网衣的变形

在圆形网衣的数值模拟中，三种不同的来流速度(U=0.122、0.178、0.242 m/s)通过底部配重8 g的单片圆形网衣，图10给出了code_aster计算中网衣在初始时刻和稳定状态下的位置和形状，与试验结果吻合较好，随着流速从0.122 m/s增加至0.242 m/s，圆形网衣的变形加剧，网衣的前端变形未能很好地与试验匹配，主要是因为在试验中，配重片是以圆环的形式存在的，而在数值计算中，是将圆环的配重平均分给了网衣节点，这导致了数值模型和试验模型的差异。表2为网衣的容积损失率，随着流速从0.122 m/s增加至0.242 m/s，容积损失率从20.9%增加至48.5%，这与图10的结果是符合的。

表2 容积损失率Tab. 2 Volume loss rate

图10 数值模拟和模型试验在不同流速下圆形网衣的变形对比Fig. 10 Deformation comparison of the circular net under different current velocities by the numerical simulations and the physical model tests

3.2.2 在OpenFOAM中流场情况

图11为速度U=0.242 m/s的水流通过圆形网衣的数值模拟结果。从图11(a)来看，网衣内部有一段区域流速衰减，在网的后方有相当长的流速衰减区域，流速衰减最小达到了0.12 m/s，出现在网两侧的正后方，绕流速度最大值达到了0.26 m/s，出现区域为网的两侧开外区域。 从图11(b)来看，网的内部和后方流速衰减明显，后网的后方是一片很大的流速衰减区域，网的底部绕流明显，速度分层明显。从图11整体来看，网后方的流速降低主要是由于网的阻塞作用，随着x方向距离柔性网的距离增加，两侧和底部均有绕流，网正后方速度有明显衰减。由于使用了Realizeable k-ε湍流模型，在速度的过渡层区域，急剧的速度梯度被减小，不同速度层之间的动量交换，速度分布比较平滑。

图11 速度0.242 m/s数值模拟计算网周围的流速分布Fig. 11 Flow-velocity distribution around a net calculated by numerical simulation for an incoming velocity of 0.242 m/s

3.2.3 OpenFOAM和code_aster受力对比

图12为在两个求解器中计算的阻力与试验数据互相对比，从数值来看吻合较好。在圆形网衣的计算中，结构求解器code_aster的后网来流速度是根据前网的速度衰减[14]估算得到的，将流速衰减因子r看成常数是不合理的，这与Bi等[23]的测量结果不符，因此在code_aster中计算的受力与试验数据相比存在偏差且最大误差为15.3%，出现在来流速度U=0.178 m/s。OpenFOAM将网衣节点的坐标位置导入计算域并采用多孔介质模型计算出流场分布和受力情况，大小与试验值较为符合。

图12 不同流速下数值模拟和试验数据之间的阻力比较Fig. 12 Comparison of the drag force between numerical simulations and experimental data at different current velocities

3.3 固定来流速度通过不同配重的单片网衣计算结果

3.3.1 在code_aster中生成圆形网衣的变形

对具有不同配重m(m=8、45、367 g)的网衣进行数值模拟，来流速度均为U=0.242 m/s，图13表示了code_aster的计算结果和试验结果对比，给出网衣的初始状态和变形状态，数值模拟结果和试验结果吻合较好，随着配重从m=8 g增加至m=367 g，圆形网衣的变形减小，在w=367 g的变形十分小。表3为网衣的容积损失率，随着配重从m=8 g增加至m=367 g，容积损失率从48.5%减小至2.0%，这与图13所得的结果相符合。

表3 容积损失率Tab. 3 Volume loss rate

图13 数值模拟和模型试验在同流速不同配重下圆形网衣的变形对比Fig. 13 Deformation comparison of the circular net under the same current velocity and different weights by the numerical simulations and the physical model tests

3.3.2 在OpenFOAM中流场情况

图14为OpenFOAM中的模拟结果，结果表明，网的内部和后方都有明显的流速衰减，底部有明显的绕流，过渡区域明显。相较于配重大的网衣，配重m=8 g的网衣后方流速衰减更大，这表明，网衣的屏蔽作用随着配重的增加而降低，尤其在网衣的正后方。

图14 不同配重下数值计算切面y=0的流速分布Fig. 14 Velocity distribution of slice y=0 simulated by numerical simulation in different weights

图15为沿着图6监测点的线中平均流速的分布，可以看出，流速与试验数据大小相近，变化趋势一致，最大的相对误差为 5.6%，模拟结果较好。在配重m=8 g的圆形网衣模拟中，在x=1.25 m左右流速趋向稳定，而在配重m=45 g和m=367 g的情况下，在x=0.75 m左右流速趋向稳定，配重m=8 g的稳定流速要比配重m=45 g和m=367 g更小。

图15 数值模拟和试验中测量点的平均流速的比较Fig. 15 Comparison of the average flow-velocity between numerical simulations and experimental measurement points

3.3.3 OpenFOAM和code_aster受力对比

由于在模型试验中并未提供试验数据，本文只进行两个求解器中网衣阻力的对比，图16为两个求解器的阻力对比，在OpenFOAM中的计算结果要比code_aster中大，这与图12的结果较为一致。随着配重的增加，阻力也持续加大。

图16 不同流速下数值模拟和试验数据之间的阻力比较Fig. 16 Comparison of the drag force between numerical simulations and experimental data at different current velocities

4 结 语

利用结构求解器code_aster和流体求解器OpenFOAM提出了一种更简便的柔性网衣和流场耦合计算模型，通过圆形网衣进行了验证。具体做法是，首先采用Screen模型计算网衣的水动力，通过桁架模型code_aster计算网衣的结构变形，模拟结果与试验相符，然后导出网衣节点坐标至OpenFOAM中，根据判断条件在流场的网格区域内生成多孔区域，通过加多孔介质模型源项的方法，利用OpenFOAM在不同来流和不同配重情况下进行了流场分析，网衣后方有明显的流速衰减，两侧和底部绕流明显，流速衰减与试验数据趋势一致，大小相近，说明多孔介质模型的有效性，最后比较了试验所得网衣阻力与code_aster、OpenFOAM计算的网衣阻力，误差均在合理范围之间。上述结果均显示，与Bi等[18]的模拟结果和试验数据吻合，说明了桁架模型和多孔介质模型耦合方法的准确性。

在柔性网衣与流场的耦合计算中，时间开销主要由流体求解器承担，结构求解器占比较小。由于本模型采用的是单向耦合方式，因此不必考虑两个求解器的数据交换和时间同步。其次，Screen水动力模型是由经验公式得到，在流速不高的条件下求解网衣变形是比较可靠的。因此本文提出的基于桁架模型和多孔介质模型建立的柔性网衣和流场的单向耦合方法对于实际问题的研究是经济、有效的。