杨静宇，李 靖

（上海交通大学 船舶海洋与建筑工程学院 海洋工程国家重点实验室，上海200240）

0 引言

随着经济的发展，人们的生活水平大幅提升，对水产品的需求与日俱增[1]，渔业的健康发展已成为全世界众多国家所共同关心的问题。我国从1970年代开始发展浅水网箱，到2018年全国传统近岸网箱数量已达到149.3万只，2019年深水网箱的养殖体积已经达到19 358 969 m3[2]。相比于试验，网衣的数值模拟要更加简便和经济[3]。近年来，随着科技的迅速发展，计算机的计算能力也突飞猛进，有限元方法、集中质量法和质量加弹簧的模拟方法被广泛地应用到网衣数学建模和动态模拟上。NIEDZWIED和HOPP[4]以网目为基本单位预测了网衣的形状和载荷，并计算出中层拖网的三维形状；PRIOUR等[5]利用有限元的方法计算出网衣下垂的形状，并采用有限元方法研究刚度对网衣的影响；LEE等[6]将网衣离散成弹簧和质量点的结合体，模拟出拖网在运动时的状态；TAKAGI等[7]考虑了渔网的弹性对网衣变形的影响，建立渔网的数学模型，对渔网的运动状态进行动态模拟，在此基础上开发出软件NALA来计算网衣的变形；AARASNES等[8]把柔性的网衣简化为一些平面网片铰接在一起，将这些小网片视为由纯二维平面组成，忽略了这些网片的三维特性。在国内，李玉成等[9]建立了各集中质量点的运动方程，通过Runge-Kutta-Verner方法求解运动方程，模拟结果与试验相比吻合较好；袁军亭[10]利用有限元理论，建立网衣的三维力学数学模型（3DNNM），计算出平面网衣和圆柱形网片的张力分布和总阻力，与试验结果相比误差较小；赵云鹏等[3，11]将网衣离散成质量点和弹簧的组合体，求解得到网衣在波浪中各个时刻和重力式网箱在来流中的受力和变形；孙宵峰等[12]采用隐式的Newmarkβ方法对网衣的力学模型进行求解，并建立网衣的可视化工具。

本研究将柔性网衣离散为弹簧和质量点的组合体，通过Morison公式计算网线所受水动力，并建立质量点的运动方程，基于迭代方法编程求解网衣质量点的坐标，最终得到网衣的变形和所受阻力。

1 运动方程的建立

1.1 网绳物理模型

网衣的基本单位为网绳，因此，首先建立网绳的运动方程，并将其拓展以模拟计算三维网衣。在实际应用中，组成渔网的网绳大部分都是柔性体，在水流的冲击下形状会发生改变，形状的改变又会引起作用在网绳上的受力发生变化，形状和受力之间密切相关，相互耦合。为了解决这一耦合问题，将网绳离散成无质量的圆柱形弹簧和无体积的质量点的结合体。圆柱体在水中的受力不仅与水流的相对速度也与圆柱体的方位有关，这里对该模型做出如下假设：（1）网绳只能承受拉力，不能承受压力；（2）张力在整个横截面上为常数；（3）网绳在拉伸过程中横截面积保持不变；（4）两个质量点之间的圆柱呈直线，不能发生弯曲。

完全浸没在水中的网绳受到的重力和浮力的合力FB可以表示为：

其中，L为网线长度，dw为网线的直径，ρnet和ρ分别指网线的密度和水的密度。

网线的直径通常不超过1 cm，相比常见海浪的波长是小量，因此可以利用Morison公式计算网线的水动力。Morison公式为半经验半理论公式，将静止在水中的构件的水动力分为拖曳力和惯性力两部分：

其中，FD为拖曳力，FI为惯性力，CD为拖曳力系数，A为沿来流方向的投影面积，CM为惯性力系数，V为来流的速度。

当构建在水中运动时，根据Brebbia和Waler[13]的建议，Morison公式可以写为：

其中，Cm=CM-1为附加质量力系数，VR为水质点相对于构件的移动速度，R为构件的位移。

对于聚乙烯网衣，利用Wilson[14]推导的关于网绳伸长量和张力之间的关系，网绳的拉力为：

其中，T为网线所受到的拉力，L0为网绳的原始长度，L为网绳变形后的长度，弹性系数C1和C2的选取参考Gerhard Klust[15]在《纤维绳索》中弹力与变形量之间的实验数据，对于聚乙烯材料网线弹性系数C1=345.37×106，C2=1.012 1。

对于尼龙网衣，采用H.Le Dret[16]在EZ 20 Lloyd拉力机上对直径为4 mm的尼龙网线进行实验并根据实验结果拟合得到的拉力与相对伸长量之间的关系式：

1.2 运动方程的建立

将网绳的质量和所受的重力、浮力、水动力以及拉力分配到两端的质量集中点上，如图1和图2所示，在每个质量点上建立运动方程。图2中Bi、Wi分别表示分配在第i个质量点的重力和浮力，Ti1和Ti2分别表示质量点i所受质量点i-1和i+1的拉力，Fi表示分配在第i个质量点的水动力。

图1 网绳模型示意图

图2 模型质量点受力分析示意图

根据牛顿第二定律，每一段网绳在水中的受力和加速度可以表示为：

教师要根据高职学生的实际学情，指导学生参与选择实用性强、富有趣味的传统文化内容，因材施教，灵活运用情境教学法、讨论法、朗读法、质疑探究法、多媒体教学法等多种教学方法，营造出传统文化学习的氛围，带领高职学生走进传统文化情境中，增进学生对有关传统文化的理解，促进学生人文素养的提高。例如在教学《鸿门宴》这一文章时，可以用多媒体播放作者的生平介绍及课文中出现的刘邦、项羽、张良、范增等人物生平经历、古代座次文化等知识，激发学生参与语文学习活动的兴趣，让学生更深刻地理解和掌握有关文学、文化知识，促进了高职学生人文素养的不断提升。

其中，M为网绳的质量，将式（3）代入（6）中，可以得到：

若将离散了的每一段网绳的水动力、重力和浮力平均分配到两端的质量点上在加上拉力，如图2所示，得到的第i个质量点的运动方程为：

其中，Mi为分配在第i个质量点的质量，ΔMi为分配在第i个质量点的附加质量，ai为第i个质量点的加速度，下标x、y、z分别表示在大地坐标系下x、y、z三方向的分量。

2 编程计算过程

Matlab软件是一款数学软件，有着强大的矩阵运算功能和简单的语法结构，使用Matlab软件编程来模拟绳子的受力变形过程，过程如下：

（1）输入绳子的参数：绳总长la、绳子直径dw，分段数N，网线材料弹性系数C1，C2，网线材料密度ρnet，流速V、密度ρ、阻力系数CD和附加质量力系数Cm，计算步长δt，最终迭代进度tend。

（2）将整个绳子划分为单独的点，并将这些点的位置、速率、加速度进行初始化。同时，建立时间迭代步的位置、速度和加速度矩阵，并进行初始化。

（3）依照网绳质量点i当前的位置，计算出第i个点在第j次迭代中，在大地坐标系下的受力大小。

（4）由步骤（3）得到的受力大小进一步计算出节点加速度aix、aiy、aiz；根据加速度和速度计算出在该时间段内绳子的质点的运动速度Uix、Uiy、Uiz；根据加速度和速度计算出网绳质量点的位移Rix、Riy、Riz，并更新质量点i的位置。

（5）网衣的计算结果收敛时质量点受力达到平衡状态，加速度趋于0，中间迭代过程的加速度并不影响最终的计算结果，所以两次迭代计算的加速度的相对误差小于0.1时，可认为当前位置已经足够准确，进入步骤（6）进行下一节点的迭代计算，否则，令j=j+1，进入（3）继续循环迭代。

（6）若i是最后一个点质量点则进入下一整体迭代步（t=t+δ）的计算，否则令i=i+1，j=0进入（c）中继续进行下一节点的迭代计算。

（7）重复（3）至（6）的步骤直至计算到预先设定的迭代进度tend。

在步骤（4）中采用NewMark-β方法以增加计算的稳定性。NewMark-β方法是对物体运动的加速度做出线性变化的假定，采用该方法后速度Ui和位移Ri的基本方程分别为：

其中，Ri、Ui，ai分别表示质量点i的位移、速度和加速度。令α=0.5，β=0.25这相当于加速度在δt内为常量，其值为两个迭代步（t，t+δt）两端加速度的平均值。

网衣是由纵横的网线构成，三维网衣的模型与网绳离散不同的是，质量点（i，j）周围有4个与之相连的圆柱杆件（i，j-1），（i，j+1），（i+1，j），（i-1，j）。将每个圆柱杆件的质量平均分配到两端的质量点上，并且质量点承担了圆柱杆件所受的水动力、浮力、重力以及拉力，质量点所在位置即为网衣目脚交点的所在位置，如图3所示。网衣质量点的运动方程与式（7）相同，编程求解过程与网线一致。

图3 网衣离散示意图

当流速大于0.5 m/s时，网绳轴线方向的阻力系数CDτ和垂直于轴线方向的阻力系数CDn的选取参考Casarella[17]的实验结果，分别取0.8和1.8。沿网线轴线方向的附加质量力系数Cmτ和垂直于轴线方向的阻力系数Cmn选取，参考赵云鹏[3，18]的研究成果，分别取0.0和1.0。

3 网衣的模拟

3.1 四周固定网衣的模拟

建立单位面积（1 m×1 m）无结尼龙网衣，命名为S0，S0四边固定，网目边长为0.0625 m，网目数为16×16，网线直径为0.002 m，密度为1.14 g/cm3。袁军亭[11]在其博士论文中将网衣离散成集中质量和弹簧的组合体，建立平面网衣数学力学模型3DNNM，采用Runge-Kutta方法求解运动方程，对网衣S0模拟计算。本小节将采用迭代求解的方法，在不同的流速下，计算出网衣S0的变形和受到的水动力以及网线中的拉力，并对比两种计算方法的结果。

设置计算的时间步长δt=10-4，t=48后结果收敛，这里的时间t并不是指物理时间，而是代表网衣质量点的迭代进度。结果如表1所示，其中，r为相对误差，T为网线中的最大拉力，F为网片所受的总水动力，X为节点的最大偏移量，下标ite和3D分别代表本文迭代法和3DNNM方法的计算结果。

表1 不同流速下两种计算结果的比较

网线之间的最大拉力，网衣所受的水动力以及网衣中的节点的最大偏移量随着流速的增加不断增大，不同流速下网衣的变形如图4所示，网衣的最大变形量出现在中间节点。两种计算方法相比，最大拉力之间的相对误差在10%以内，总水动力之间的相对误差在5%以内，最大偏移量之间的相对误差在10%以内，当流速大于0.8 m/s后相对误差在5%以内，表明通过迭代法求解网衣的变形的结果是可靠的。

3.2 四角固定网衣的模拟

建立模拟网衣S1，S1四角固定，方形网目大小为4 cm（网目目脚长度为4 cm），网线直径为1.18mm，网衣边长为80 cm×32 cm，网目数量为20×8，材料为聚乙烯。对S1进行模拟，设置流速分别为1.0 m/s、1.5 m/s和2.0 m/s。计算中迭代步长δt=2×10-5，t=50后质量点的坐标收敛，停止计算，计算结果如表2所示。结果表明，流速越大，网衣中网线的最大拉力、网衣所受的总水动力以及网衣的最大偏移量越大，不同流速下网衣的变形如图5所示。不论是网衣S0还是网衣S1，流速越大，网衣的变形量越大，但表1和表2中Fite/V2的数值均越小。根据网衣阻力系数CDn的计算公式（11），网衣的面积An、水流的密度ρ和网线的阻力系数在模拟中保持恒定，表明影响网衣阻力系数的是形变量，且网衣的形变量越大其阻力系数越小。

图5 不同流速下S1的变形

表2 不同流速下四角固定网衣S1计算结果

4 结语

通过将网绳离散成无质量的弹簧和质量点的模型，利用Morison公式求解网线在来流中的水动力，将网线的受力分配集中到质量点上，在质量点上建立运动方程。采用迭代方法求解各个的运动方程，并利用NewMark-β算法使迭代过程更加稳定。将网线的求解模型拓展到三维网衣，在四周固定网衣S0模拟计算中，与通过Runge-Kutta求解运动方程的3DNNM模型的结果相比，网线之间的拉力和网衣中最大的变形量的相对误差均在10%以内，网衣的水动力的相对误差在5%以内，表明本文通过迭代法求解网衣的变形的方法是可靠的，具有工程应用价值。在四角固定网衣S1的计算结果中发现，网衣的网绳最大拉力、总水动力和最大形变量随着流速的提高都逐渐增大。模拟结果表明，不论是四周固定网衣S0还是四角固定网衣S1的阻力系数均随着网衣变形量的增大而逐渐减小。