基于Fluent的内孤立波数值模拟及其作用下顶张力立管动力响应分析

Numerical Simulation of Internal Solitary Waves Based on Fluent and Dynamic Response Analysis of Top Tension Risers

摘要: 本文基于Fluent软件，建立具有多种设计振幅的内孤立波数值水槽，通过质量源造波方法实现对内孤立波的模拟，与物理实验结果及理论结果进行对比分析，利用软件自带的监测功能在对大振幅内波进行模拟时获取波致流场数据，结合改进的Morison方程计算内孤立波对顶张力立管的作用，基于欧拉–伯努利梁模型利用有限元法分析其动力响应。结果表明，本文采用的质量源造波方法下的模拟结果，与实验及理论吻合较好，在大振幅内孤立波作用下，立管位于上层流体中的部分所受到的来自内孤立波的水平作用力相比于位于下层流体中的部分大了许多，而在内孤立波作用下的顶张力立管发生的变形和产生的应力都很大，在波谷到达立管时应力与变形达到最大。

Abstract: In this paper, a numerical flume of internal solitary wave with various design amplitudes is estab-lished based on the Fluent software. The internal solitary wave is simulated by mass source wave-generating method and the simulation results are compared with the physical experiment and theoretical results. The data of wave-induced flow field are obtained by monitoring function of the software. Combined with the modified Morison equation, the effect of internal solitary wave on the top tension riser is calculated, and its dynamic response is analyzed by finite element method. The results indicate that the mass source wave-generating method adopted in this paper is in good agreement with the experimental and theoretical results. Under the large amplitude internal solitary wave, the horizontal force of internal solitary wave on the riser in the upper fluid is obviously greater than that on the riser in the lower fluid. The riser under the action of internal solitary wave will undergo large deformation and produce large stress. When the valley reaches the riser, the stress and deformation reach the maximum.

1. 引言

在实际海况中，海水温跃层及盐跃层的存在导致了海水密度垂向分布的不均匀，即垂向上的密度跃层，而内孤立波就是发生在密度稳定层化的海水内部的波动，它是一种海洋中的普遍现象 [1] 。内孤立波的振幅远大于表面波，对海洋工程结构物的影响极为显著，因此，对内孤立波进行相关研究具有重要的工程意义与经济价值。

随着计算流体动力学(CFD)这一领域的蓬勃发展，将CFD方法运用到内孤立波相关模拟中，研究海上结构物与内孤立波之间的相互影响作用问题，已然成为了当下的热门研究方向之一。研究表明，CFD方法不仅可以模拟内孤立波与海上结构物的作用过程中发生的分裂演化现象，还可以获取作用在结构物上的内孤立波荷载等水动力特性，这无疑为内孤立波与海洋结构物相互作用问题的研究又提供了一个新的途径 [2] 。

内孤立波理论随着时间变迁而日趋完善，当下广泛使用的主要有KdV理论，eKdV理论，mKdV理论和MCC理论等。蔡树群等结合KdV理论引入了Morison公式，来计算内孤立波作用下小尺度圆柱桩柱所受的作用力与力矩 [3] 。尤云祥等基于mKdV理论，结合Morison公式研究了内孤立波与张力腿平台的相互作用 [4] 。王旭等结合eKdV理论发展了内孤立波质量源造波方法，所得结果与理论及实验较为吻合 [5] 。Wang等基于内孤立波的KdV、eKdV和MCC理论适用性条件，将内孤立波诱导上下层深度平均水平速度作为入口条件，建立了两层流体中内孤立波对直立圆柱体强非线性作用的数值模拟方法，与相应理论和实验结果一致 [6] 。

目前的研究大部分基于速度入口法进行CFD数值模拟，与理论结果对比分析，而与物理实验对比验证，并将数值模拟下的内孤立波作用于实际尺度的立管，分析其动力响应的则相对较少。鉴于此，本文的研究立足于流体软件Fluent，采用两个点源形式的质量源造波方法，从粘性不可压缩流体的N-S方程出发，结合内孤立波理论中的KdV理论和eKdV理论，建立有限深两层流体中的内孤立波数值模拟方法，在将所得数值模拟结果与理论、实验相对比，验证该方法可靠性的前提下，将质量源造波方法应用到大振幅内孤立波的模拟中，并将模拟所得的内孤立波数据结合改进的Morison方程，求解实际尺度的海洋顶张力立管的动力响应并加以分析。

2. 造波理论与方法

2.1. 理论模型

相关实验研究表明了内孤立波理论的适用范围，发现KdV理论仅适用于振幅水深比小于0.1的情况，eKdV理论适用于振幅水深比超过0.1的广泛振幅范围 [7] 。因此在选取内孤立波理论时，根据后文所要研究的工况，主要选用KdV理论与eKdV理论，采用有限水深情况下两层流体的简化模型，上下层流体深度分别为 $h_1,h_2$ ，各自的密度分别为 $$ ，总水深 $h=h_1+h_2$ ，如图1所示。 $l$ 为内波特征波长， $c$ 为内波相速度，以 $\eta \left(x,t\right)$ 表示内孤立波的波面函数，为振幅大小，当波形为下凹波时为负值，上凸波时为正值。

Figure 1. Sketch of theoretical models for internal solitary waves

图1. 内孤立波理论模型示意图

波面方程的KdV理论解为 [8] ：

$\eta \left(x,t\right)={\eta }_0\sec h^2\left(\frac{x-ct}{l}\right)$

其中：

$\begin{array}{l}c=c_0+\frac{c_1}{3}{\eta }_0\\ l=\sqrt{\frac{12c_2}{c_1{\eta }_0}}\\ c_0=\sqrt{\frac{g{h}_1{h}_2\left({\rho }_2-{\rho }_1\right)}{{\rho }_1{h}_2+{\rho }_2{h}_1}}\\ c_1=-\frac{3c_0}{2}\frac{{\rho }_1{h}_2^2-{\rho }_2{h}_1^2}{{\rho }_1{h}_1{h}_2^2+{\rho }_2{h}_2{h}_1^2}\\ c_2=\frac{c_0}{6}\frac{{\rho }_1{h}_2{h}_1^2+{\rho }_2{h}_1{h}_2^2}{{\rho }_1{h}_2+{\rho }_2{h}_1}\end{array}$

式中： $c$ 和 $l$ 分别为KdV理论下的内孤立波的相速度和特征长度。

波面方程的eKdV理论解 [9] 为：

$\eta \left(x,t\right)=\frac{{\eta }_0}{B+\left(1-B\right){\cosh }^2\left(\frac{x-c_{e}t}{l_e}\right)}$

其中：

$\begin{array}{l}{c}_e=c_0+\frac{{\eta }_0}{3}\left(c_1+\frac{1}{2}{c}_3{\eta }_0\right)\\ l_e=\sqrt{\frac{12c_2}{\left(c_1+\frac{1}{2}{c}_3{\eta }_0\right){\eta }_0}}\\ B=\frac{-c_3{\eta }_0}{2c_1+c_3{\eta }_0}\end{array}$

式中： $c_e$ 和 $l_e$ 分别为eKdV理论下的内孤立波的相速度和特征长度。

2.2. 造波方法

本文采用质量源造波法建立二维内孤立波数值水槽，造波结构示意图如图2所示。

Figure 2. Sketch of the mass source wave-generating method

图2. 质量源造波法示意图

该法是在连续性方程中加入附加质量源项，对其进行特定设置，以此来产生目标流场。上层流体中源S₁的强度为 $S_1\left(x,z,t\right)$ ，下层流体中源S₂的强度为 $S_2\left(x,z,t\right)$ ，通过上层质量源不断释放流体，下层质量源不断吸收流体，从而产生下凹型的内孤立波。

流体控制方程在加入质量源项后，形式变为下式：

$\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=S\left(x,z,t\right)\\ \frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\mu \frac{\partial u_i}{\partial x_j}\right]-\frac{\partial p_i}{\partial x}+\rho g_i\end{array}$

式中， $$ 为流体密度， $\mu$ 为动力粘度系数，g为重力加速度。

假定由质量源所引起的流体质量的增减与孤立波波面的变化相对应，则对于源S₁、S₂有：

${\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle {\int }_{S1}{S}_1\left(x,z,t\right)\text{d}{S}_1\text{d}t=-{\rho }_1{\displaystyle \underset{0}{\overset{t}{\int }}c\eta \left(t\right)\text{d}t}}}$

${\displaystyle \underset{0}{\overset{t}{\int }}{\displaystyle {\int }_{S2}{S}_2\left(x,z,t\right)\text{d}{S}_2\text{d}t={\rho }_2{\displaystyle \underset{0}{\overset{t}{\int }}c\eta \left(t\right)\text{d}t}}}$

质量源为点源，在质量源区域内均匀分布，则可以认为质量源强度仅为时间t的函数，则两源的源强表达式如下：

$S_1\left(t\right)=-\frac{{\rho }_{1}c\eta \left(t\right)}{A_1}$

$S_2\left(t\right)=\frac{{\rho }_{2}c\eta \left(t\right)}{A_2}$

其中A₁、A₂分别为源S₁、S₂的面积。

3. 模型建立与验证

3.1. 数值水槽的建立

为了便于与实验结果进行对比，本文在建立数值水槽时，参照所在课题组在中国海洋大学物理海洋实验室进行的内孤立波物理实验。本实验采用Oster双缸法制取密度分层的流体，采用重力塌陷法制造不同振幅的内波 [10] 。水槽长15 m，宽0.35 m，高0.7 m，额定水深0.6 m，上下层流体高度分别为9.5 cm、48.5 cm，密度分别为1035 kg/m³、1054 kg/m³。

3.2. Fluent设置

边界条件设置：数值水槽的左侧边界设置为对称边界，右侧边界设置为压力出口边界，水槽顶部与底部都设为固壁边界。

UDF设置：利用Fluent中的DEFINE_SOURCE宏编写理论UDF，将源项添加到造波源区域，嵌入到造波区控制方程中。

液面追踪方法：采用Volume of Fluid方法即VOF法来追踪流体中上下层分界处的内孤立波波面。VOF定义为目标流体的体积与网格体积的比值，只要知道这个函数在每个网格上的值，就可以实现对运动界面的追踪。当某一单元体被目标流体占满，其F值为1；当它被另一相非目标流体占满，其F值为0；当它同时含有目标流体和非目标流体时，其F值在0到1之间，表明该单元是位于两流体交界处的单元。

求解设置：压力速度耦合方式采用PISO算法，压力插值选取体积力加权法，动量方程、湍动量方程和湍动耗散率方程均采用二阶迎风格式，多次模拟得知这样的设置可以提高解的精确性。二相流体界面构造方法选用几何重构法，时间步长设置为0.01 s，每个时间步最大迭代次数为20次。

3.3. 数值造波模拟结果及分析

为便于与实验结果进行对比，在数值模拟时选取的工况设置如表1所示。根据振幅水深比的不同，各工况选取适合的内波理论进行模拟。

Table 1. Setting of numerical simulation condition

表1. 数值模拟工况设置

在进行模拟时，首先通过对x轴上不同位置处波高变化的监测，找出达到设计振幅时内波所在的位置。之后通过对该位置进行监测，得到随时间变化的z轴方向上的波高。图3表示了三种工况下波高随时间的变化，截取了内波一个周期中的主要部分来展示其波形。其中工况1和工况2振幅达到设计振幅时的位置在x = 150 cm处，工况3则在x = 180 cm处。

由图3可以看出，工况1的数值造波波形与实验结果和理论解都吻合得较好，工况2和工况3的数值造波波形在前半段与实验结果、理论结果都吻合得较好，在后半段相较于实验结果，与理论解更为一致。这一结果表明本文所采用的质量源造波方法对振幅水深比小于0.1和大于0.1的内孤立波均有较好的适用性，根据其各自对应的合适的内孤立波理论进行数值模拟所得波形及振幅与实验解和理论解都较吻合，能有效实现对内孤立波的数值模拟。

Figure 3. Comparison of numerical waveform with experimental and theoretical results

图3. 数值造波波形与理论及实验结果对比

4. 大振幅内波下立管的动力响应

4.1. 模型验证

参照实际海况资料，本文进行计算的内孤立波参数如表2所示 [11] 。数值水槽总长度为5000米，其中工作区长度3000米，消波区长度1000米。

Table 2. Parameters of internal solitary wave

表2. 内孤立波参数

利用前文所述方法，基于eKdV理论，在Fluent中利用两个质量点源进行造波，运用Fluent的监测功能将所得波高随时间变化的相关数据导出，并将CFD数据下的数值造波波形图与理论情况下的波形进行对比，对比结果如图4所示。

Figure 4. Comparison of numerical waveform with theoretical results

图4. 振幅40 m时数值模拟结果与理论解波形对比

在该模拟中，设计振幅为40 m，通过监测对比知，在x = 130 m处达到该模拟振幅最接近设计振幅的值，为39.7 m。图4展示了这一位置处的波面高度z随时间t的变化情况，由图可以看出，本文所采用的数值模拟方法下的波形与理论波形基本一致，在振幅大小、波高变化等方面吻合度较高，说明利用该方法对大振幅内波的模拟是可行的。

4.2. 结果与分析

参照相关文献，顶部张紧式立管(TTR)的参数选取如表3所示 [12] 。

Table 3. Parameters of top tension riser

表3. 顶张力立管参数

在计算小尺度结构物所受波浪力时通常选用Morison公式，故本文对立管所受水平作用力的计算采用改进的Morison公式，其表达式如下 [13] ：

$f_x=\frac{1}{2}{C}_D\rho D\left(U-\stackrel{\dot{}}{x}\right)\left|U-\stackrel{\dot{}}{x}\right|+C_M\rho \frac{\text{π}{D}^2}{4}\left(\stackrel{\dot{}}{U}-\ddot{x}\right)$

式中： $\stackrel{\dot{}}{x}$ 与 $\ddot{x}$ 分别为立管顺流向运动的速度和加速度， $\rho$ 为外部流体密度，D为立管外径， $C_D$ 与 $C_M$ 分别为拖曳力系数和惯性力系数， $U$ 与分别为内孤立波引起的水质点的速度和水平加速度。

在进行计算时，立管顶张力取立管浮重的1.2倍。立管总长240米，全部位于水下。将立管按长度方向划分成长3米的单元，共划分出80个单元，81个节点。Fluent计算中将内波时间步长设置为0.1 s，为保持一致，在进行立管响应计算时时间步也设置为0.1 s。

由Fluent数据结果知，在第482 s时内孤立波波谷到达立管轴线，内波计算时共历时1600 s。将监测所得的内孤立波流场数据与Morison方程、立管控制方程结合，采用有限单元法在时域中对大振幅内孤立波作用下顶张力立管的动力响应进行数值模拟，最终利用MATLAB编写程序进行结果展示。

通过上述方法，得出了不同时刻立管整体位移和应力沿垂向的分布图，分别如图5、图6所示。

图5和图6中不同的线条代表不同的时刻，每条线则代表在该时刻，在内波的作用下，立管整体的位移情况与应力情况。在初始时刻，内波振幅较小，对立管的作用力微弱，无论是立管的位移还是应力，都处在一个基本没有变化的状态。后期随着波谷靠近，作用力变大，立管的位移与应力都有了明显的变化。

Figure 5. Displacements of riser at different moments

图5. 不同时刻立管整体位移图

Figure 6. Stresses of the riser at different moments

图6. 不同时刻立管应力图

在两图列出的五个时刻中，波谷到达时刻为482 s，120 s和300 s为波谷到达立管之前的两个时刻，680 s和1000 s为波谷到达立管之后的两个时刻。从图中可以看出处于上层流体部分的立管，其顺流向位移及应力都明显大于处于下层流体的部分，原因是内孤立波波致流场中上层的水平流速大于下层水平流速，使得立管上半段产生较大的水平作用力。同时可以看出，在内孤立波的波谷靠近立管的过程中，立管的位移与应力都在逐渐增大。在内孤立波的波谷到达立管时，立管的位移与应力均达到最大，而后随着内孤立波的离去，立管的位移和应力又都逐渐减小。

图7为沿立管垂向不同高度处的节点，其位移随时间变化的图像。选择的节点分别是距离底端42 m、72 m的位置，两层流体中间120 m的位置，内孤立波最大振幅所在的170 m的位置，以及上下层流体分界处的210 m的位置。从图中可以看出，不同深度的节点，其位移均随时间发生不同程度的变化，位移最大值均出现在波谷到达时刻，这是因为波谷到达时振幅最大，对立管的作用效果最显著。而在t相同的同一时刻，几处不同高度的节点，位移最大的都是z = 210 m处的节点，这是因为该节点位于上下两层流体的交界面处，受到内波的作用最大。

Figure 7. Time-histories of displacement for nodes in different depth

图7. 不同深度处节点位移时程图

图8和图9为内波的波谷到达立管时刻，即对立管的作用最强时立管的位移图与应力图，由图8可知在波谷到达时上下层流体分界面处立管位移最大，最大位移达到立管直径的10倍。由图9可知内孤立波波谷到达时上下层流体分界面处立管的应力最大，最大应力约为93 Mpa。

Figure 8. Displacement of riser with trough arrives

图8. 波谷到达时刻立管位移图

Figure 9. Stresses of riser with trough arrives

图9. 波谷到达时刻立管应力图

5. 结论

本文基于Fluent软件，建立具有多种设计振幅的内孤立波数值水槽，通过质量源造波方法实现了对内孤立波的模拟，与物理实验结果及理论结果进行对比分析，验证了造波方法的可行性与准确性。

利用软件的监测功能在对大振幅内波进行模拟时获取波致流场数据，结合改进的Morison方程计算内孤立波对顶张力立管的作用，通过传统的有限元方法得到立管的动力响应，结果表明在大振幅内孤立波作用下，上层流体中的立管所受内孤立波水平作用力明显大于下层流体中的立管，内孤立波作用下的立管会发生较大变形，产生较大应力，尤其在波谷到达立管时应力与变形达到最大，严重影响其安全性能与正常使用性能，在进行立管设计安装以及海洋工程实际中，考虑内孤立波的影响十分重要。

基金项目

国家重点研发计划项目(No. 2016YFC0802301)；山东省自然科学基金博士基金(ZR2017BEE041)；山东科技大学人才引进科研启动基金项目(2016RCJJ038)。

NOTES

^*通讯作者。