Chinese Journal of Ship Research
多液舱晃荡与养殖工船时域耦合运动的数值模拟
引用格式：肖凯隆,陈作钢.多液舱晃荡与养殖工船时域耦合运动的数值模拟[J]. 中国舰船研究, 2020, 15(1): 136–144. XIAO K L, CHEN Z G. Numerical simulation of aquaculture ship motions coupled with tanks sloshing in time domain[J]. Chinese Journal of Ship Research, 2020, 15(1): 136–144.肖凯隆1,3，陈作钢*1,2,3 1上海交通大学海洋工程国家重点实验室，上海 200240 2高新船舶与深海开发装备协同创新中心，上海200240 3上海交通大学船舶海洋与建筑工程学院，上海 200240
摘 要:［目的］为研究多液舱养殖工船在横向波浪激励下，其液舱晃荡与船体运动的时域耦合问题，对其进行数值模拟。［方法］采用脉冲响应函数法预报船体运动，基于CFD理论模拟计算液舱晃荡问题，并讨论壁面剪切力对液舱晃荡的影响。为验证数值方法的可靠性，进行载液工船模型横摇自由衰减试验，并将试验结果与数值结果进行对比。同时，在此基础上，提出多液舱模型的二维简化方法。［结果］计算结果表明：在横向波浪激励下，波浪频率为0.5 rad/s时，液舱晃荡使养殖工船横摇幅值响应算子峰值增大了43%。液舱晃荡可能会给载液船舶的安全性带来风险，应该在设计阶段充分考虑。［结论］研究成果可为多液舱船的优化设计提供参考。关键词：多液舱晃荡；养殖工船；脉冲响应函数；耦合运动；时域模拟中图分类号: U661.32 文献标志码:A DOI：10.19693/j.issn.1673-3185.01709
Numerical simulation of aquaculture ship motions coupled with tanks sloshing in time domain
XIAO Kailong1,3, CHEN Zuogang*1,2,3
1 State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China 2 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, Shanghai 200240, China 3 School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
Abstract: ［Objectives］ For multi-tank aquaculture ships, the numerical method is used to study the coupling problem between tank sloshing and hull motion under transverse wave excitation in the time domain. ［Methods］ Hull motion is predicted using the impulse-response function method, the sloshing phenomenon is simulated based on the CFD theory, and the influence of the wall shear force on sloshing is discussed. In order to verify the effectiveness of the numerical method, a free-rolling experiment on a ship model with liquid-loaded tanks is carried out, and the experimental results are compared with the numerical results, showing a good match. Furthermore, a two-dimensional simplified model of the multi-tank model is proposed ［Results］ The calculation results show that the sloshing phenomenon increases the maximum value of the rolling amplitude response operator by 43% under transverse wave excitation with wave frequency at 0.5 rad/s. ［Conclusions］ The sloshing phenomenon may bring risks to the safety of liquid-loaded ships, and should be fully considered during the design process. Key words: tank sloshing； aquaculture ship； impulse-response function； coupling motion； time domain simulation
收稿日期: 2019–08–19 修回日期: 2019–12–16 网络首发时间: 2020–03–30 14:11基金项目:中国水产科学研究院渔业机械仪器研究所资助项目（15DZ1202100）作者简介: 肖凯隆，男，1995 年生，硕士生。研究方向：船舶水动力学。E-mail: 530570352@sjtu.edu.cn陈作钢，男，1967年生，博士，研究员，博士生导师。研究方向：计算流体力学与船舶水动力学。E-mail: zgchen@sjtu.edu.cn *通信作者:陈作钢
0 引 言
我国是海水养殖发达国家之一，养殖面积和产量均居世界首位[1]。现阶段，我国海水养殖方式还比较粗放，且内陆和沿海的养殖空间受到挤压，导致养殖密度过大、病害频发、水产品品质下降和环境恶化等问题日益突出[2]。因此，走向深远海、发展绿色养殖生产新方式已成必然趋势。深远海养殖工船作为深远海养殖的一个发展方向，集成了繁育、养殖、加工、冷冻冷藏等鱼货物供给的一条龙功能，可为阶段性的分舱养殖、自动化投喂、自动排污、机械化起捕等提供先进的养殖和生产手段。依托养殖工船进行深远海养殖，可以实现高度集约化、生态化、规模化的健康养殖，同时将鱼类养殖区域推向深远海，能拓展养殖空间，有效地推进海洋渔业的转型发展[3]。养殖工船通常系泊工作于水质和风浪条件合适的深远海海域，并汲取海水至养殖液舱。但液舱晃荡产生的晃荡力作用于船体，可能会加剧船体运动，使船舶安全性面临更大的挑战。因此，在设计阶段对液舱晃荡与船体运动的耦合作用进行充分论证，是确保养殖工船安全性的关键。目前，国内外许多学者针对液舱晃荡与船体时域耦合运动展开了研究。其中，液舱晃荡的模拟方法主要有2类：非线性势流方法[4] 和黏流方法[5]。船体运动预报主要基于势流理论，采用脉冲响应函数（IRF）法。船体运动预报和液舱晃荡模拟分别采用势流和黏流的方法，简称为势流—黏流耦合法；如果都采用势流方法，则简称为全势流法。全势流法假设舱内液体无黏无旋，但当晃荡剧烈时，该假设就会出现局限性。势流—黏流耦合法是基于CFD理论模拟液舱晃荡，以牺牲计算效率为代价来提高计算精度。势流—黏流耦合法兼顾了势流方法计算速度快以及CFD 方法相对准确的优点，基于此，李裕龙等[6] 和操戈等[7]将计算结果与试验结果进行对比，验证了该方法的可靠性。本文将采用势流—黏流耦合法计算液舱晃荡与船体运动耦合的问题。首先，通过SESAM 软件计算船体水动力系数及波浪力，基于IRF 法计算时延函数，建立船体运动预报方程。接着，通过 FLUENT软件数值模拟液舱晃荡，采用流体体积（VOF）法捕捉自由面，并将液舱晃荡力实时添加至船体运动方程，让船体运动方程求解与液舱晃荡模拟交替进行，以随着时间的步进，重复解耦工作，然后在自由横摇衰减试验的基础上验证该方法的有效性。最后，计算得到船体运动幅值响应算子（RAO），对比分析考虑液舱晃荡力与不考虑液舱晃荡力下的计算结果，探究液舱晃荡力对船体运动的影响。同时，还将研究壁面剪切力对数值模拟液舱晃荡准确度的影响，并针对多液舱模型，探索一种二维计算方法，以极大地减小计算规模。
1 载液船舶在波浪上的运动预报理论1.1 速度势叠加理论
势流理论假设流体无黏无旋，并引入速度势φ描述流场。φ满足式中，∇为拉普拉斯算子。根据叠加原理，可以把速度势分解为入射势φI、绕射势 φD、辐射势φR，如果只保留一阶项，则
得到
各项速度势沿船体表面积分可得到流体对船体的各项作用力。其中，入射力FI是假设船体不存在时入射波在船体本该存在的位置产生的力；绕射力FD是假设船体固定时入射波绕射过船体所产生的力；辐射力FR是假设水面静止时船体运动所产生的力。因此，一阶波浪力F(1)为H
1.2 船舶频域运动方程
根据牛顿第二定理，刚性船体的运动方程为
式中：M 为船体质量矩阵； ξ¨为船体加速度；FS= −Kξ ，为船体回复力，其中K为回复系数矩阵，ξ为船体位移；FL为船体受到的其他外力，如液舱晃
荡力。
式中： Ma为船体附加质量矩阵；C为船体阻尼系数矩阵；ξ˙为船体速度。将式（3）和式（5）代入式（4），得到船体频域运动响应方程
式中： ω为波浪频率；FI,D 为 Froude-Kriloff 力，为入射力与绕射力之和，表示船体在规则波中受到的波浪激励力[8]。
2 时域耦合方程的数值计算方法
本文将基于水动力学软件SESAM 计算船体质量、附加质量、阻尼系数和波浪载荷等水动力系数，基于CFD 软件 FLUENT 数值模拟液舱晃荡，并在用户自定义函数（UDF）中实时将晃荡力添加至船体运动方程。具体运算步骤为： 1） 采用 SESAM 软件计算船体水动力系数，直接输入或等待UDF 读取； 2） 给液舱设置微小的初始速度作为第1个时间步的边界条件； 3） 采用软件 FLUENT 计算当前时间步的液舱晃荡力，并添加至船体运动方程； 4） 求解船体运动方程，给出下一时间步的船体运动速度，并将其作为液舱新时间步下的运动边界条件； 5） 步骤3）和步骤4）交替进行，随着时间的步进，重复解耦工作，直至计算结束。UDF的编制及计算过程需考虑以下细节： 1） UDF 在软件 FLUENT 中运行时，新时间步下自定义物理量会被重新定义，而求解船体运动方程需要参考历史物理量，因此要将有关物理量定义成静态变量，以避免在新时间步被清空。2） 为避免每一时间步都要重新计算时延函数，可编写外部程序预先计算，并在UDF 运行的初始时刻读取。3） 船体水动力是基于船体运动重心位置确定的，转换液舱晃荡力至船体运动方程或将船体运动响应转换至液舱晃荡边界条件时，需要根据实际情况确定。
3 不规则多液舱简谐晃荡的数值模拟3.1 养殖工船物理模型
养殖工船垂线间长 241 m、型宽 48 m、型深18.5 m，经快速性多目标数值优化[11]，得到型线图如图1所示。在主甲板下设置了关于船舯对称的 6对养殖水舱，其中艏、艉2对为不规则液舱，船舯4对为尺寸一致的矩形规则液舱，各液舱高度均为 16.1 m ，载液率为75%。图 2 所示为液舱简化俯视图，图中仅绘制了1对规则液舱来代替4对规则液舱。
横浪工况下，船舶以横摇和垂荡运动为主[12]，当只计算这2 个自由度时，4对规则液舱的流动
情况一致。因此，可用1 对规则液舱晃荡力的4倍代替4对规则液舱的晃荡力，以减小计算规模。
3.2 壁面剪切力对液舱晃荡数值模拟的影响
正确处理壁面黏性效应，是CFD数值模拟的关键。壁面函数法是常用的壁面处理方法，其本质是在黏性子层及过渡层使用称为壁面函数的半经验公式计算壁面与充分发展湍流区域之间的黏性影响区域。壁面函数的使用依赖于壁单位（Wall unit）y+的正确选取，y+表示第一层网格与壁面间的无量纲距离，其计算公式可参考文献[13]。当壁面函数选取增强壁面处理（ Enheanced wall treatment）时，要求 y+保持在 5以下。壁面函数法对y+的要求导致壁面附近的网格量剧增、计算时间变长。若壁面黏性效应对液舱晃荡影响极小，不妨将壁面设置为滑移壁面，此时壁面附近没有强剪切，离壁第1层网格的高度也不再受壁面函数的约束。本文设置了3组案例（Case1, Case2, Case3)来对比验证壁面剪切力对液舱晃荡数值模拟的影响。3组案例均采用30 m×16 m的二维矩形液舱，载液率均为50% ，液舱绕其中的点作频率为0.085cos0.6t rad/s 的横摇运动。其壁面函数的选择与网格划分情况如表1所示。
Case1 与 Case2 均采用壁面加密结构网格，如图 3所示。要满足y+小于 5，Case1 采用增强壁面处理的壁面函数法；Case2 选择滑移壁面作为边界条件，忽略壁面剪切力；Case3采用均匀的结构网格，同样忽略壁面剪切力。图 4 为 3组案例相对回转中心的晃荡力矩系数 d的计算结果。由图可见，3条曲线基本重合，但由于网格单元数减少了90% 以上，导致 Case3的计算时间减少了约80%。因此，可以认为壁面剪切力对液舱晃荡的数值模拟影响不大，可以选择滑移壁面作为壁面条件，同时采用均匀的结构网格，降低网格数量，提高计算效率。
4 数值求解方法的有效性验证
为了验证本文数值方法的有效性，本节对养殖工船模型在载液率为62.5% 时的静水横摇自由衰减试验进行了数值模拟。
4.1 横摇自由衰减运动的模型试验
试验在上海交通大学循环水槽中进行，该水槽试验段长 8m，宽 3m。根据水槽尺寸，试验模型缩尺比选为1∶100。试验开始前，利用惯量架测量船模转动惯量及重心位置。安装、调整压载铁，使吃水位于设计水线处，并记录压载铁的数量和位置。试验时，人工施加恒力于一侧船舷，使船体产生稳定的初始横倾角，当船体内、外水静止后释放船体，船体开始做自由横摇衰减运动，同时用陀螺仪记录横摇角时历。图5所示为载液及压载铁安装示意图。
4.2 数值模拟结果与试验结果对比
在已知船模、压载铁、液舱水体的重心位置和质量的情况下，根据初稳性理论可以求出横摇回复力系数K。假设船体绕着重心所在轴线做自由横摇，可以求出转动惯量。采用SESAM 软件计算横摇阻尼系数及附加转动惯量。FLUENT数值模拟液舱晃荡中，经网格无关性验证，设置液舱模型网格尺寸为6 mm，时间步长为 0.005 s，迭代步数为 60。图 6所示为试验与数值模拟得到的横摇角时历对比图，图中同时还绘制了不考虑液舱晃荡的模拟结果。可以发现，考虑液舱晃荡的数值模拟结果与试验结果吻合良好，证明本文数值方法对液舱晃荡与船体运动耦合方程的求解有效。
5 实船计算结果与分析5.1 船体水动力系数、时延函数及CFD 参数设置
横摇临界阻尼系数λ取3% 时，SESAM 软件计算得到的横摇阻尼系数C44 如图 7(a) 所示；由式（10）计算出的时延函数 R44如图 7(b)所示；波高2m时的一阶波浪横摇力矩F44 如图 7(c) 所示。液舱晃荡的数值模拟中，经网格无关性验证，网格模型尺寸取为 0.6 m，时间步长为 0.005 s，每一时间步的迭代步数为20。初始时刻液舱自由面示意图如图8所示。
5.2 RAO时域计算结果
图 9所示为横向规则波激励下，考虑液舱晃荡力和不考虑液舱晃荡力情况下计算得到的横摇RAO4及升沉 RAO3。如图 9(a) 所示，当波浪频率较低时，液舱晃荡力会加大船舶横摇运动幅值；而在
波浪频率较高时，横摇幅值有所降低，液舱起到了减摇水舱的作用。同时，考虑液舱晃荡力时， RAO4的峰值较不考虑液舱晃荡力时的情况增大了 43% ，且峰值对应的波浪频率有所降低。如图 9(b) 所示，当波浪频率较低时，液舱晃荡力会加大船舶升沉运动幅值；而在波浪频率较高时，船舶运动幅值有所降低。考虑液舱晃荡力后， RAO3在频率更低的位置达到了峰值，且峰值略大于不考虑液舱晃荡力的情况。
5.3 液舱晃荡力对养殖工船运动的影响
如图 9(a) 所示，在波浪频率ω=0.4 ， 0.6 rad/s的激励情况下，液舱晃荡力会分别增大、减小船体横摇运动幅度。本节将对比分析横摇运动时历及有关力矩，以说明液舱晃荡力对船体运动的影响机制。图 10(a) 所示为在波浪频率ω=0.4 rad/s 的横向规则波激励下，考虑液舱晃荡力和不考虑液舱晃荡力的船体横摇运动时历ξ4(t) 对比，图 10(b)所示为晃荡横摇力矩与波浪横摇力矩时历F4(t)的对比。从中可以看出，运动稳定后，晃荡横摇力矩与波浪横摇力矩的相位十分接近，此时，液舱晃荡力将增大船体横摇运动幅度。图 11(a) 所示为在波浪频率ω=0.6 rad/s 的横向规则波激励下，考虑液舱晃荡力和不考虑液舱晃荡力的船体横摇运动时历对比，图11(b) 所示为晃荡横摇力矩与波浪横摇力矩的时历对比。从中可以看出，运动稳定后，晃荡横摇力矩与波浪横摇力矩相位相差较大，此时，液舱晃荡力将降低船体横摇运动幅度。
6 多液舱晃荡的二维简化计算
本文虽然用1对规则液舱代替了4对规则液舱，但受限于计算规模，仍有必要进一步简化计算模型。图12所示为波浪频率 ω=0.5 rad/s 时，艏、艉液舱横摇力矩与4对规则液舱横摇力矩的对比，从中可以看出，前者数值较小，因而对船体运动的影响有限。图 13所示为液舱总横摇力矩与规则液舱横摇力矩的比值k随频率ω的变化示意图。如图所示，k的变化范围很小，可近似认为总力矩与规则液舱力矩间的比值恒定。因此，可选取某一波浪频率下的k作为所有频率下的定值，并将k 倍的规则液舱力矩近似看作总力矩，以避免计算艏、艉不规则舱。波浪频率ω=0.5 rad/s 对应的k 值为1.131，选其为本文定值 k。
矩形液舱长、宽比大于1 时，三维液舱晃荡与二维液舱晃荡的数值模拟结果经过长度换算后吻合1.6，可数值模拟其二维横截面液舱晃荡，并将晃荡力乘以液舱长度代替三维液舱的晃荡力，以进一步简化模型。通过引入参数k替代艏、艉不规则舱的晃荡模拟以及采用二维模型代替三维液舱模型的策略，可将计算效率提升90% 以上。图 14所示为在波浪频率 ω=0.3，0.7 rad/s 的激励下，通过二维方法计算的与包含不规则液舱的三维方法计算得到的船体横摇运动时历对比。虽然 k 取波浪频率 ω=0.5 rad/s 对应的值，但在波浪频率 ω=0.3 ， 0.7 rad/s情况下的横摇运动模拟中，二维计算结果与三维计算结果吻合良好，表明用单一波浪频率下的k近似代替其他频率下的k是可行的。同时也表明，用二维液舱模型代替三维规则液舱模型来计算液舱晃荡是可行的。
如图 9（a）所示，在波浪频率 ω=0.5 rad/s 时 RAO4达到峰值，该频率下，波幅为1m的横向波浪可使船体横倾角最大达到11.5°。图 15所示为在该波浪激励下，船体运动稳定后左倾、右倾达到最大位置以及在中间位置处（从左到右）的液舱晃荡7 结 论本文求解了多液舱养殖工船在横向规则波激励下的时域响应，得出如下主要结论： 1） 考虑壁面剪切力与否对液舱晃荡的数值模拟影响不大，壁面条件可选择滑移壁面，同时采用均匀结构网格代替边界加密网格，提高计算效率。2） 载液量较大时，液舱晃荡力会对船体运动产生明显影响。液舱晃荡力与波浪力相位接近时，将增大船体运动幅度；反之，将减小。对于养殖工船，液舱晃荡力对船体横摇运动影响较大。不考虑液舱晃荡力时，横摇RAO 在波浪频率ω=0.6 rad/s 达到峰值，为 0.14 rad/m ；考虑液舱晃荡力时，在波浪频率ω =0.5 rad/s 达到最大，为0.20 rad/m，增大了 43%。3） 养殖工船艏、艉不规则舱液量较少，对船体运动影响有限，可将某一波浪频率下，液舱总晃荡力与规则液舱晃荡力的比值视作该比值在所有波浪频率下的值，从而避免计算艏、艉不规则液舱的晃荡。对于规则舱，可用其二维横截面模型代替三维液舱模型，进一步简化计算。参考文献： [1] 蔡计强, 张宇雷, 李建宇, 等. 10万吨级深远海养殖平
台总体技术研究[J]. 船舶工程, 2017, 39(增刊 1): 198–203. CAI J Q, ZHANG Y L, LI J Y, et al. General technology research of 100 thousand ton deep sea aquaculture platform[J]. Ship Engineering, 2017, 39(Supp 1): 198–203 (in Chinese). [2] 徐皓, 谌志新, 蔡计强, 等.我国深远海养殖工程装备发展研究 [J]. 渔业现代化, 2016, 43(3): 1–6. XU H, CHEN Z X, CAI J Q, et al. Research on the de
云图。从中可以看出，由于周期较大，横摇角度较小，即使在 RAO4峰值对应的波浪条件下，舱内水体的晃荡依然比较平稳，自由液面也未发生明显形变。velopment of deep sea aquaculture engineering equipment in China[J]. Fishery Modernization, 2016, 43(3): 1–6 (in Chinese). [3] 徐皓.水产养殖设施与深水养殖平台工程发展战略[J].中国工程科学, 2016, 18(3): 37–42. XU H. Development strategy for aquaculture facility and deepwater aquaculture platform[J]. Engineering Science, 2016, 18(3): 37–42 (in Chinese). [4] KIM Y, NAM B W, KIM D W, et al. Study on coupling effects of ship motion and sloshing[J]. Ocean Engineering, 2007, 34(16): 2176–2187. [5] LEE S J, KIM M H, LEE D H, et al. The effects of LNG-tank sloshing on the global motions of LNG carriers[J]. Ocean Engineering, 2007, 34(1): 10–20. [6] 李裕龙, 朱仁传, 缪国平, 等. 基于 OpenFOAM 的船舶
与液舱流体晃荡在波浪中时域耦合运动的数值模拟[J].船舶力学, 2012, 16(7): 750–758. LI Y L, ZHU R C, MIAO G P, et al. Simulation of ship motions coupled with tank sloshing in time domain based on OpenFOAM[J]. Journal of Ship Mechanics, 2012, 16(7): 750–758 (in Chinese). [7] 操戈, 李旭, 张咏鸥, 等. FPSO液舱晃荡与船舶时域耦
合运动数值模拟[J]. 中国舰船研究, 2015, 10(1): 88–96. CAO G, LI X, ZHANG Y O, et al. Numerical simulation of ship motion coupled with tank sloshing for FPSO[J]. Chinese Journal of Ship Research, 2015, 10(1): 88–96 (in Chinese). [8] 朱仁传, 缪国平. 船舶在波浪上的运动理论 [M]. 上海:上海交通大学出版社, 2019: 67-79. ZHU R C, MIAO G P. Theory of ship mtion in waves[M]. Shanghai: Shanghai Jiao Tong University Press, 2019: 67-79 (in Chinese). [9] 周立师.排水型船横摇阻尼的理论与试验研究 [D]. 哈尔滨:哈尔滨工程大学, 2012. ZHOU L S. Theoretical and experimental research on
roll damping of the displacement type vessel[D]. Harbin: Harbin Engineering University, 2012 (in Chinese). [10] FALTINSEN O M. A numerical nonlinear method of sloshing in tanks with two-dimensional flow[J]. Journal of Ship Research, 1978, 22(3): 193–202. [11] FENG Y K, CHEN Z G, DAI Y, et al. Multidisciplinary optimization of an offshore aquaculture vessel hull form based on the support vector regression surrogate model[J]. Ocean Engineering, 2018, 166: 145–158.
[12] 洪亮, 朱仁传, 缪国平, 等.波浪中船体与液舱晃荡耦合运动的时域数值计算 [J]. 哈尔滨工程大学学报, 2012, 33(5): 635–641, 647. HONG L, ZHU R C, MIAO G P, et al. Numerical calculation of ship motions coupled with tank sloshing in time domain based on potential flow theory[J]. Journal of Harbin Engineering University, 2012, 33(5): 635–641, 647 (in Chinese). [13] 章梓雄, 董曾南. 粘性流体力学 [M]. 北京: 清华大学出版社, 1998: 291−297. ZHANG Z X, DONG Z N. Viscous fluid mechanics[M]. Beijing: Tsinghua University Press, 1998: 291−297 (in Chinese). [14] 文攀.船舶与船舶液舱晃荡耦合运动分析 [D]. 武汉:华中科技大学, 2013. WEN P. Analysis on coupling effects of ship motion and tank sloshing[D]. Wuhan: Huazhong University of Science and Technology, 2013 (in Chinese).