王旭，林忠义，尤云祥

(1.上海交通大学海洋工程国家重点实验室，上海200240;2.嘉兴南洋职业技术学院，浙江嘉兴314003)

Spar平台是当今深海资源开发中的主力深海平台之一，由于通常永久系泊于特定海域进行作业，海洋环境对其安全性有很大影响，因此必须对其在各种海洋环境中下的水动力性能进行研究［1］。南海油气资源丰富，已成为我国深海资源开发的主战场，但南海内孤立波活动频繁。1990年，在流花油田就曾发生过因内孤立波导致缆绳拉断、船体碰撞，甚至拉断和挤破漂浮软管的事故［2］。同年，在南海陆丰油田也发生过因内孤立波导致半潜钻井船与锚定油轮在连接输油管道时发生困难等问题［3］。因此，内孤立波已成为南海深海资源开发中需要考虑的海洋环境因素之一。内孤立波是一种最大振幅发生在密度稳定层化海洋内部的波动，由于非线性和色散效应在一定尺度上的平衡，在其传播过程中可以保持波形和传播速度不变，一般地可以用 KdV(Korteweg-de Vries)、eKdV(extended KdV)［4］和MCC(Miyata-Choi-Camassa)［5］等理论模型来描述。这3类理论中弱非线性和弱色散这2个条件仅仅为定性描述，为此，黄文昊等以系列实验为依据给出了这两个条件的定量表征方法［6］。

直立圆柱型结构是Spar等各种深海平台的主体结构形式，因此研究内孤立波作用下直立圆柱体的载荷特性，具有现实的理论和工程意义。程友良等［7］和蔡树群等［8-9］将Morison公式与KdV理论结合，而Xie等［10］则将Morison公式与MCC理论结合，研究了内孤立波作用下直立圆柱体的载荷特性问题。尤云祥等［11-12］将Morison公式与eKdV理论结合，研究了内孤立波作用下张力腿和半潜式平台的载荷与动力响应问题，而宋志军等［13］则将Morison公式与KdV理论结合，研究了内孤立波作用下Spar平台的载荷与动力响应问题。需要指出的是，在这些文献中，关于Morison公式中惯性力和拖曳力系数都是参照表面波的方法选取的，但这种选取方法缺乏理论和实验依据。为此，黄文昊等以系列实验为依据对圆柱型结构及张力腿平台给出了这两个系数的选取方法［14-15］。计算流体力学(computational fluid dynamics，CFD)应用方面，关辉等［16-18］分别利用 KdV、MCC、eKdV 理论解作为初始条件，采用CFD方法分别研究了内孤立波作用下水下潜体和海洋立管的载荷特性等问题。但需要指出的是，在这些文献中，由于没有考虑KdV、eKdV和MCC理论的适用性条件等问题，致使其数值模拟结果均不同程度地出现了内孤立波振幅及其波形不可控等问题。

鉴于此，本文采用Navier-Stokes方程，结合文献［6］中提出的KdV、eKdV和MCC理论的适用性条件，建立振幅及其波形可控的内孤立波CFD数值模拟方法。对内孤立波与直立圆柱体的相互作用特性进行数值模拟，进而分析内孤立波作用下直立圆柱体各种载荷成分的形成机理及其影响程度，直立圆柱体对内孤立波的波形及其流场特性的影响，以及Morison公式的适用性等问题。

1 数值方法

考虑两层流体中内孤立波与直立圆柱体的相互作用问题。设各层均为不可压流体，上层流体深度与密度分别为h1和ρ1，下层流体深度与密度分别为h2和ρ2。内孤立波为平面前进波，界面位移为ζ，沿ox轴正方向传播，圆柱体直径为D，吃水为d。建立直角坐标系oxyz，其中oxy平面位于流体静止时两层流体的界面上，oz轴与圆柱体中心轴重合且以竖直向上为正。

采用求解Navier-Stokes方程的方法对内孤立波诱导流场进行数值模拟，其中流场控制方程为

式中:(u1，u2，u3)为速度，p为动压力，t为时间，ρ为流体密度，当 ζ＜z＜h1时，ρ=ρ1;当-h2＜z＜ζ时，ρ=ρ2。

将圆柱体直立置于两层流体中，圆柱体壁面取为无滑移不可穿透边界条件，而水面和水底要求满足如下壁面条件:

流场计算的控制区域如图1所示，包括内孤立波生成传播区和消波区2个区域。采用速度入口方法生成内孤立波，当造波区中形成稳定的内孤立波后，对所生成内孤立波的传播特性进行监测分析，并对直立圆柱体的内孤立波载荷进行计算。设相速度为c，则其诱导上下层流体中的层深度平均水平速度分别为［5］

速度入口处速度分布为:在上层流体入口处，速度取为u-

1;在下层流体入口处，速度取为u-2。在内孤立波生成与传播过程中，两层流体的界面会发生变化，采用VOF(volume of fluid)方法追踪两层流体界面的变化。利用海绵层消波方法对水槽尾部的内孤立波进行消波处理，该方法在消波区中通过在动量方程右端添加源项-μ(x)ui的方式实现，其中μ(x)为海绵层衰减系数。

图1 数值水槽示意图Fig.1 Schematic diagram of numerical flume

直立圆柱体的内孤立波水平力由摩擦力和压差力组成:

直立圆柱体的内孤立波垂向力同样由表面摩擦力和压差力组成，如下:

式中:S为直立圆柱体浸湿侧表面积;式中第一项为直立圆柱体侧表面和底部的载荷摩擦力;第二项为直立圆柱体侧表面的载荷压差力;Ut为表面切向速度，(nx，ny，nz)为直立圆柱体表面法向分量，正方向指向直立圆柱体内部。

在数值模拟中，如何选择合适的内孤立波理论来计算入口速度是一个需要解决的关键问题。为此，定义非线性参数 ε=a/h和色散参数 μ=(h/λ)2，其中λ为内孤立波特征宽度。根据上下层流体深度、密度及内孤立波振幅等条件，利用KdV、eKdV和MCC理论解，分别计算相应的内孤立波非线性参数ε和色散参数μ，根据文献［6］中3类内孤立波的适用性条件，入口速度(3)的计算方法如下。

利用商业软件Fluent进行数值模拟与分析。采用有限体积法离散动量和连续性方程，对流项离散采用QUICK(quadratic upstream interpolation for convective kinetics)格式，压力插值格式采用彻体力加权(body force weighted)方法，压力速度耦合迭代采用PISO(pressure implicit with splitting of operators)算法，两层流体界面的构造方法选用几何重构法。初始时间步长为Δt=0.005 s，计算过程中根据每个时间步长的收敛情况逐渐增加时间步长以缩短计算时间。

2 结果与分析

文献［14］利用大型密度分层水槽，对内孤立波作用下直立圆柱体的载荷特性进行了系列实验。其中，圆柱体直径D=0.15 m，吃水d=0.535 m，实验水槽长为30 m，水深为1 m，上下层流体密度分别为ρ1=998 kg/m3和ρ2=1 025 kg/m3，上下层流体深度比分别为h1∶h2=1 ∶9、2 ∶8、3 ∶7，本文将结合该文献中的相关实验结果进行数值模拟与分析。为此，数值水槽主尺度、上下层流体密度及其深度比均与该文一致。在数值模拟中，内孤立波生成传播区的长度为18 m，而消波区的长度为12 m，圆柱体中心轴距速度入口端9 m。计算区域采用六面体结构化网格进行离散，沿圆柱体周向的单元数为80个，总的单元数量为1 468 472个。

2.1 内孤立波数值模拟结果

对内孤立波生成传播特性进行数值模拟与分析。在利用CFD方法对内孤立波进行数值模拟时，流体粘性对其生成与传播特性的影响是一个需要重点考虑的重要问题。为此，在数值模拟中考虑了如下两种情况:一种为考虑流体粘性的情况，运动粘性系数取为ν=1.0×10-6m2/s，称为 N-S模拟;另一种为不考虑流体粘性的情况，即运动粘性系数取为ν=0，称为Euler模拟。

图2中给出了当h1∶h2=3∶7和ad/h=0.101时，采用N-S方程和Euler方程对内孤立波生成与传播特性的数值模拟结果，其中ad为内孤立波设计振幅。由图可知，在入口处上层流体向右而下层流体向左运动过程中，为保持流体质量的守恒，上层流体产生向下塌陷现象，形成一个下凹型鼓包;在约化重力作用下凹型鼓包逐渐形成一个向水槽右方传播的内孤立波，在有粘和无粘两种情况下，数值模拟所得内孤立波在其传播过程中均保持波形稳定、振幅衰减很小，没有明显的尾波现象。结果表明，在有粘和无粘两种情况下，数值模拟所得内孤立波振幅与其设计振幅的相对误差均在5%以内，因此对内孤立波的CFD数值模拟，采用基于N-S和Euler方程的两种方法均是合适与可靠的。为此，在下文中，如无特别声明，所有数值模拟均是在有粘情况下进行的。

图2 两种不同模拟工况下内孤立波数值模拟结果Fig.2 The numerical results for internal solitary wave waveforms between two cases

图3 内孤立波波形数值模拟结果与理论和实验结果比较Fig.3 Comparisons of the numerical results for internal solitary wave waveforms with theoretical and experimental ones

图3给出了在3种工况下，内孤立波波形的数值模拟结果，并与相应理论和实验结果进行了比较，实验结果取自文献［14］。其中，工况A为h1∶h2=3∶7和ad/h=0.101，此时内孤立波为中等非线性和弱色散的，选择eKdV理论计算入口速度;工况B为h1∶h2=2∶8和ad/h=0.052，此时内孤立波为弱非线性和弱色散的，选择KdV理论计算入口速度;工况C为h1∶h2=1∶9和ad/h=0.086，此时内孤立波为强非线性和弱色散的，选择MCC理论计算入口速度。结果表明，在各工况下，数值模拟所得内孤立波的波形，不仅与初始内孤立波理论解波形一致，而且与实验所得波形吻合，这表明依据3类内孤立波理论的适用性条件，采用本文所述数值模拟方法所得内孤立波的波形是准确可控的。

2.2 内孤立波载荷特性

为后文陈述方便，记Fx、Fz和My分别为内孤立波作用下直立圆柱体的水平力、垂向力及其力矩，其中力矩转动中心至圆柱体底部距离为d-=d+0.15。定义和分别为无因次内孤立波水平力、垂向力及其力矩。

图4 内孤立波无因次水平力、垂向力及其力矩幅值数值与实验结果比较Fig.4 Results of numerical and experimental amplitudes for dimensionless loads(horizontal forces，vertical forces，torques)due to internal solitary waves

图5给出了在A工况下，内孤立波无因次水平力、垂向力及其力矩时历的数值模拟结果，并与文献［14］中的实验结果进行了比较。由图可知，直立圆柱体内孤立波载荷时历数值结果与实验时历结果吻合，表明采用本文所述直立圆柱体内孤立波载荷的计算方法合理可行。结果同时表明，对水平力、垂向力绝对值及其力矩，在内孤立波波谷到达直立圆柱体中心轴之前的某个时刻达到其最大值，在该时刻之前随时间增加而增大，在该时刻之后则随时间增大而减小;在A工况下，直立圆柱体底部始终位于内孤立波的波面下方，由于在波面下方内孤立波诱导水平速度方向与其传播方向相反，因此波面下方的流体动压力是负的，从而导致直立圆柱体的垂向力是负值。

图5 Case A工况内孤立波无因次水平力、垂向力及力矩时历特性Fig.5 The time-variant characteristics for dimensionless loads(horizontal forces，vertical forces，torques)due to internal solitary waves for Case A

下面考虑直立圆柱体内孤立波水平和垂向力中压差力和摩擦力两种成分的特性。在图6中，给出了A工况下，内孤立波无因次水平压差力Fpx及其摩擦力fx，无因次垂向压差力Fpz及其摩擦力fz时历特性的数值模拟。结果表明，无论是直立圆柱体内孤立波水平力还是垂向力，摩擦力与压差力相比都是一个小量，可以忽略。这意味着直立圆柱体内孤立波水平和垂向力的主要成分为压差力。

图6 内孤立波无因次压差力及摩擦力时历特性Fig.6 The time-variant characteristics for pressure forces and friction forces due to internal solitary waves

对直立圆柱体的内孤立波水平和垂向压差力，可以进一步将其分为波浪压差力和粘性压差力两个部分。其中，波浪压差力与内孤立波诱导的水质点波动有关，可以采用Euler模拟方法得到;粘性压差力为由于流体粘性效应导致的压差力增量，将N-S模拟所得压差力减去波浪压差力即得粘性压差力。在图7中，给出了在A工况下，内孤立波无因次水平波浪压差力及其粘性压差力，无因次垂向波浪压差力及其粘性压差力时历特性的数值模拟。结果表明，对直立圆柱体内孤立波水平压差力，其波浪压差力和粘性压差力量级相当;对直立圆柱体内孤立波垂向压差力，其粘性压差力与波浪压差力相比是一个小量，可以忽略。这意味着对直立圆柱体的内孤立波水平力必须考虑流体的粘性效应，但对其垂向力则可以忽略流体的粘性效应。

图7 内孤立波无因次波浪压差力及粘性压差力时历特性Fig.7 The time-variant characteristics for wave pressure forces and viscous pressure forces due to internal solitary waves

3 结论

基于两层流体KdV、eKdV和MCC理论解的适用性条件，以Navier-Stokes方程为流场控制方程，以内孤立波在上下层流体中诱导的深度平均水平速度作为入口条件，建立了内孤立波与直立圆柱体强非线性作用的数值模拟方法。结果表明:

1)该方法对直立圆柱体内孤立波水平力、垂向力和力矩幅值及其时历变化特性的数值模拟结果与相应实验结果一致，可以用于内孤立波与海洋浮式结构强非线性作用的数值模拟。

2)直立圆柱体内孤立波水平和垂向力由波浪压差力、粘性压差力和摩擦力组成。其中，水平摩擦力、垂向摩擦力及其粘性压差力均很小，可以忽略。

3)水平力的主要成分为波浪压差力和粘性压差力，而垂向力的主要成分为波浪压差力。对于直立圆柱体内孤立波水平力必须考虑流体的粘性效应，但对于其垂向力则可以忽略流体的粘性效应。

［1］黄维平，白兴兰，孙传栋，等.国外Spar平台研究现状及中国南海应用前景分析［J］.中国海洋大学学报，2008，38(4):675-680.HUANG Weiping，BAI Xinglan，SUN Chuandong，et al.State of the art of studies on spar platforms abroad and outlook of their application in the South China Sea［J］.Periodical of Ocean University of China，2008，38(4):675-680.

［2］陈景辉.南海流花11-1深海油田开发工程［J］.中国海洋平台，1996，11(1):43-45.CHEN Jinghui.Liuhua 11-1 deepwater fields development project in the South China Sea［J］.China offshore Platform，1996，11(1):43-45.

［3］BOLE J B，EBBESMEYER C C，ROMEA R D.Soliton currents in the South China Sea:measurements and theoretical modeling［C］//The 16th Offshore Technology Conference.Houston，America，1994:367-376.

［4］HELFRICH K R，MELVILLE W K.Long nonlinear internal waves［J］.Ann Rev Fluid Mech，2006，38:395-425.

［5］CHOI W，CAMASSA R.Fully nonlinear internal waves in a two-fluid system ［J］.J Fluid Mech.，1999，396:1-36.

［6］黄文昊，尤云祥，王旭，等.有限深两层流体中内孤立波造实验及其理论模型［J］.物理学报，2013，62(8):084705:1-14.HUANG Wenhao， YOU Yunxiang， WANG Xu， et al.Wave-making experiments and theoretical models for internal solitary waves in a two-layer fluid of finite depth［J］.Acta Phys.Sin.，2013，62(8):08475:1-14.

［7］CHENG Y L，LI J C，LIU Y F.Recent Advances in fluid Mechanics［M］.Berlin:Springer，2004.

［8］CAI S Q，LONG X M，GAN Z J.A method to estimate the forces exerted by internal solitons on cylinder piles［J］.O-cean Engineering，2003，30(5):673-689.

［9］CAI S Q，LONG X M，WANG S G.Force and torques exerted by internal solitions in shear flows on cylinder piles［J］.Applied Ocean Research.2008，30(1):72-77.

［10］XIE J S，JIANG Y J，YANG L G，et al.Strongly nonlinear internal solution load on a small vertical circular cylinder in two-layer fluids［J］.Applied Mathematical Modelling，2010，34(8):2089-2101.

［11］尤云祥，李巍，时忠民，等.海洋内孤立波中张力腿平台的水动力特性［J］.上海交通大学学报，2010，44(1):12-17.YOU Yunxiang，LI Wei，XIANG Zhongmin，et al.Hydrodynamic characteristics of tension leg platforms in ocean internal solitary waves［J］.Journal of Shanghai Jiaotong University，2010，44(1):12-17.

［12］尤云祥，李巍，胡天群，等.内孤立波中半潜平台动力响应特性［J］.海洋工程，2012，30(2):1-7.YOU Yunxiang，LI Wei，HU Tianqun，et al.Dynamic responses of a semi-submersible platform in internal solitary waves［J］.The Ocean Engineering，2012，30(2):1-7.

［13］宋志军，勾莹，滕斌，等.内孤立波作用下Spar平台的运动响应［J］.海洋工程，2010，32(2):12-19.SONG Zhijun，GOU Ying，TENG Bin，et al.The motion responses of a Spar platform under internal solitary wave［J］.The Ocean Engineering，2010，32(2):12-19.

［14］黄文昊，尤云祥，王旭，等.圆柱型结构内孤立波载荷实验及其理论模型［J］.力学学报，2013，45(5):1-13.HUANG Wenhao，YOU Yunxiang，WANG Xu，et al.Internal solitary wave loads experiments and its theoretical model for a cylindrical structure［J］.Chinese Journal of Theoretical and Applied Mechanics，2013，45(5):1-13.

［15］黄文昊，尤云祥，王竟宇，等.张力腿平台内孤立波载荷及其理论模型.上海交通大学学报，2013，47(10):1494-1502.HUANG Wenhao，YOU Yunxiang，WANG Dingyu.Internal solitary wave loads and its theoretical model for a tension leg platform［J］.Journal of Shanghai Jiaotong University，2013，47(10):1494-1502.

［16］关晖，苏晓冰，田俊杰.分层流体中内孤立波数值造波方法及其应用［J］.力学季刊，2011，32(2):218-223.GUAN Hui，SU Xiaobing，TIAN Junjie.Numerical method of solitary wave generation in a stratified flow and its applications［J］.chinese quarterly of mechanics，2011，32(2):218-223.

［17］高原雪，尤云祥，王旭，等.基于MCC理论的内孤立波数值模拟［J］.海洋工程，2012，30(4):29-36.GAO Yuanxue，YOU Yunxiang，WANG Xu，et al.Numerical simulation for the internal solitary wave based on MCC theory［J］.The Ocean Engineering，2012，30(4):29-36.

［18］陈杰，尤云祥，刘晓东，等.内孤立波与有航速潜体相互作用数值模拟［J］.水动力学研究与进展，2010，25(3):344-350.CHEN Jie，YOU Yunxiang，LIU Xiaodong，et al.Numerical simulation of interaction of internal solitary waves with a moving submarine［J］.Chinese Journal of Hydrodynamics，2010，25(3):344-350.