悬链线管体三维流场数值模拟研究

2021-10-27何建勇高洋洋王立忠沃恩海张卓先

海洋工程 2021年5期

何建勇，高洋洋，王立忠，沃恩海，张卓先

(1. 浙江大学海洋学院，浙江舟山 316021； 2. 浙江省海港投资运营集团有限公司，浙江宁波 315040)

悬链线式细长结构在深海油气悬链线立管、锚泊系统及海上风电管缆等海洋工程中广泛应用。超细长悬链线管体结构由于自重和环境荷载作用，展向曲率沿程变化，相比于竖直管体，循环往复的海流诱发变曲率管体产生涡激振动的过程更为复杂。上部平台大幅运动下会引起悬链线管体结构触底点位置不断变化，导致管体倾斜角度和曲率剧烈变化。在涡激振动与平台运动作用下极易引起悬链线管体产生疲劳破坏，导致立管油气泄漏、锚链走锚等事故发生。

目前关于管体绕流流场特性研究主要针对竖直圆柱和倾斜圆柱，开展了大量试验研究和数值模拟。研究表明圆柱流场特性主要和雷诺数密切相关，随着雷诺数的增大，圆柱尾流会从层流转变为湍流，转变过程中伴随着尾流过渡、剪切层过渡和边界层过渡，最明显的特征是不稳定性。Williamson[1-2]研究发现圆柱尾流旋涡形成区域的不稳定性导致了尾流的紊乱。该不稳定性包括小规模流向涡的形成，大规模的旋涡变形以及剪切层的不稳定。随着雷诺数的增大，旋涡脱落会呈现不连续性、斑状大尺度结构等特征，后驻点压力对尾流的不稳定性特征更为敏感[3]。Karniadakis和Triantafyllou[4]研究发现在雷诺数Re=200附近尾流特征由二维过渡为三维，随着雷诺数增大，靠近圆柱展向涡呈波浪状，下游涡结构则呈肋条状。Mittal和Balachandar[5]研究表明二次涡的存在对圆柱表面压力分布没有显著的影响。Behara和Mittal[6]研究在150≤Re≤350区间内圆柱尾流的转捩特性时发现，波浪状的展向涡充分发展必然引起旋涡错位，导致涡脱频率随时间变化，升阻力系数均值相应减小。Shlrakashi等[7]试验研究了倾斜角度变化对圆柱尾流中卡门涡街的干扰作用，当倾角增大时，圆柱轴向出现的二次流会导致涡脱频率减小，旋涡的规则性降低。Matsuzaki等[8]通过风洞试验发现减小长径比或增大倾斜角度，圆柱底端脱落的旋涡对整体展向涡的干扰作用增强，沿圆柱轴向涡脱强度增大。Razali等[9]试验研究了倾斜角度对圆柱尾流速度和涡量三维特征的影响，当倾斜角度为45°时，圆柱展向出现轴向流，尾流中存在大量流向涡，轴向流干扰卡门涡街的形成，导致涡脱变形，尾流的三维特性增强。Zhou等[10]研究发现当圆柱的倾斜角度大于15°时，尾流中出现二次轴向涡，当倾斜角度增大到45°时，速度频谱图中峰值减小，谱峰值区域变宽。Zhao等[11]通过三维数值模拟发现随着倾斜角度的增大，展向涡形状和圆柱轴线保持平行，圆周时均压力系数和无量纲旋涡脱落频率减小。Wang等[12]研究发现随着倾斜角度增大，展向涡的一致性明显降低而流向涡一致性增加，在向下游发展的过程中展向涡之间的距离增大，导致涡脱之间相互作用减弱。

对于悬链线管体，沿展向曲率变化较大，其尾流特征比竖直和倾斜管体更为复杂。Miliou等[13]数值研究了雷诺数Re为100和500条件下弯管尾流中涡形成及脱落特征，结果表明，弯管凹形布置时，由于弯管展向强烈的轴向流存在，导致尾流中涡脱被完全抑制。Gallardo等[14]采用直接数值模拟方法研究了Re=3 900下弯管凸形布置时曲率对弯管尾流的影响，结果表明弯管竖直段的涡管形态与弯管轴线保持平行，而靠近水平段由于局部倾角较大，导致涡脱强度减弱，尾流的三维特性增强。Assi等[15]通过粒子图像测速技术研究了Re=1 000条件下弯曲管体的尾流形态，研究表明凹形布置时水平段的存在会使来流出现分离，同时产生的旋涡冲击弯曲段，造成管体整体尾流湍动性增强，而凸形弯曲管体展向涡脱相关性强，整体涡脱更加规律。

综上所述，已有研究主要集中于雷诺数变化对竖直和倾斜圆柱流场形态及受力特性的影响研究，对于悬链线管体弯曲段尾流流场特征及受力特性研究较少。由于局部倾角沿弯管展向变化，产生的轴向流不仅使尾流中展向涡发生倾斜，同时会对涡脱起到抑制作用，引起时均压力系数和回流区的变化。采用计算流体力学开源代码OpenFOAM开展了不同雷诺数条件下悬链线弯曲管体绕流三维流场特性研究，主要分析了雷诺数对凸形布置弯管瞬时、时均流场特征，表面压力分布和涡脱频率变化特性的影响规律。

1 数学方法

1.1 计算模型

对于不可压缩黏性流体，在笛卡尔坐标系中，三维Navier-Stokes(N-S)方程和连续性方程为：

(1)

(2)

式中:xi为x,y,z三个方向上的坐标；ui为xi方向的速度；t，ρ，p和ν分别表示时间、流体密度、压力和流体运动黏度。

大涡模拟采用滤波函数将流场中的旋涡分解为大小尺度两部分。大尺度结构可以通过直接求解N-S方程得到，小尺度结构用亚格子模型求解。经过滤波处理的N-S方程和连续性方程为：

(3)

(4)

采用有限体积法进行离散。在100≤Re≤500和1 000≤Re≤3 900范围内分别采用直接数值模拟和大涡模拟方法开展弯管绕流流场三维数值模拟。在低雷诺数100≤Re≤500范围内，时间项采用二阶Crank-Nicolson隐式离散，对流项、压力项和拉普拉斯项均采用二阶高斯线性格式离散。在高雷诺数1 000≤Re≤3 900区间，时间项离散采用backward二阶隐式，对流项采用高斯线性格式离散，压力项和拉普拉斯项采用高斯线性limited格式离散，且均为二阶精度。

1.2 边界条件

将悬链线管体简化成由长度6.0D的竖直段，半径为12.5D的四分之一圆环以及10.0D的水平段构成的弯管模型，其中D=0.01 m为弯管竖直段截面直径，单弯管绕流的计算域如图1所示，计算域尺寸为32.5D×20.0D×28.5D(长×宽×高)。入口边界至弯管竖直段轴线的距离为10.0D，底部边界至水平段轴线的距离为10.0D，两侧边界距离弯管竖直段轴线为10.0D，弯管的曲率半径设置为R=12.5D，与Miliou等[13]、Gallardo等[14，16]曲率保持一致，竖直段和水平段的长度分别为6.0D和10.0D，以消除顶部和出口边界对弯曲段的影响[17-18]。s为弯曲段的弧长，即s=RΦ，其中Φ为从竖直段末端开始沿弯曲段逆时针方向的角度。定义竖直段与弯曲段交接处s/D=0，水平段与弯曲段交接处s/D=19.6。

图1 单弯管绕流计算域示意Fig. 1 Computational domain of flow past single curved cylinder

入口边界采用狄利克雷边界条件，流向速度u分别设置为0.01 m/s、0.03 m/s、0.05 m/s、0.10 m/s、0.15 m/s 和0.39 m/s，对应雷诺数Re=100、300、500、1 000、1 500和3 900，横向速度与展向速度设置为v=w=0。出口边界采用纽曼边界条件，流速沿出口方向梯度和压力均设置为0。弯管表面采用无滑移边界条件。两侧边界和上下边界设置为自由滑移边界条件。当雷诺数处于100≤Re≤500和1 000≤Re≤3 900时，无量纲时间步长ΔtU/D分别取0.01和0.001，计算过程中最大库朗数CFL维持在0.55左右，以保证数值模拟的稳定性和收敛性。

1.3 网格敏感性验证

首先在Re=500和3 900条件下，采用三套不同密度的网格对竖直圆柱和弯管竖直段开展了网格敏感性研究，以确保后续模拟结果的准确性和可靠性。竖直圆柱展向高度与弯管竖直段保持一致，即z=6.0D。计算域整体网格和弯管局部网格如图2所示。

图2 计算网格示意Fig. 2 Schematic view of the computational mesh

阻力系数CD、升力系数CL、压力系数Cp和斯特劳哈尔数St的定义分别为：

(5)

(6)

(7)

(8)

式中:FD和FL分别为圆柱受到的阻力和升力；U为来流速度；A为圆柱沿流向的投影面积，对于竖直圆柱而言，A=LD，L为圆柱展向高度，对于弯管而言，A=(R+D/2+Lv)D，Lv是弯管竖直段的长度；p∞为入口处的压力大小；fs为旋涡脱落频率。

表1 Re=500和3 900时网格敏感性研究结果Tab. 1 Results of mesh independence test for Re=500 and 3 900

表2 数值计算结果与文献结果对比Tab. 2 Comparisons between the numerical results and previous literature

文中弯管的几何尺寸与Miliou等[13]一致，将模拟得到的瞬时涡量图、时均流线图和时均压力系数与文献[13，29]结果(如图3～4所示)对比，进一步验证文中数值模型的可靠性。如图3所示，当Re=100时，靠近弯管的涡管轴线垂直，向下游发展的过程中，涡管发生弯曲变形，文中模拟结果和Miliou等[13]结果非常吻合。时均流线图和Canabes[28]的结果基本一致，从图3(b)中可以发现Re=100时的回流区主要集中在弯管轴线后方2D左右的范围内，沿展向回流区面积逐渐减小。如图4所示，文中得到的前后驻点线时均压力系数变化趋势和文献结果一致，即时均压力系数绝对值沿展向逐渐减小，该值在弯曲段末端约为-0.1，且Re=500时后驻点线时均压力系数绝对值大于Re=100的时均压力系数。

图3 Re=100时弯管尾流特征Fig. 3 Wake flow for a curved cylinder at Re=100

图4 展向时均压力系数Fig. 4 Comparisons of the time averaged pressure coefficients along the spanwise

2 结果与讨论

2.1 瞬时流场特征

图5 不同雷诺数下弯管的涡量等值面三维图和俯视图Fig. 5 Instantaneous iso-surface of vorticity contours of curved cylinder at different Reynolds numbers

图6 不同雷诺数下竖直圆柱的涡量等值面三维图和俯视图Fig. 6 Instantaneous iso-surface of vorticity contours of straight cylinder at different Reynolds numbers

图7 Re=100时的瞬时涡量Fig. 7 Instantaneous spanwise vorticity along the span at Re=100

图8 Re=3 900时的瞬时涡量Fig. 8 Instantaneous spanwise vorticity along the span at Re=3 900

图9 竖直圆柱瞬时涡量Fig. 9 Instantaneous spanwise vorticity along the span

图10 不同雷诺数下不同水平位置处的流向瞬时涡量Fig. 10 Instantaneous streamwise vorticity along the x-direction at different Reynolds numbers

2.2 时均流场特征

图11为Re=100、500和3 900时沿展向不同位置处的时均流线拓扑，如图所示，尾流中形成的回流区关于弯管中轴线对称，随着雷诺数的增大，同一展向位置处的回流区逐渐减小。当Re=100时，在z/D=27.5，22.5和17.5处的回流区距弯管圆心2D左右，而Re=500和3 900时，相同位置处的回流区向下游分别发展了1.5D和1.0D。由于z/D=12.0位置处局部倾角(Φ=56.44°)较大，回流区消失。对于弯管而言，沿展向局部倾角增大，靠近弯曲段末端的过程中回流区逐渐消失。

如图12所示为不同雷诺数下弯管中心剖面位置处(y/D=10)的时均流线拓扑。在不同雷诺数条件下，上游区域的时均流线变化较为一致，即流线基本平行分布，方向和来流保持一致，流线形状在接近弯曲段的过程中发生变化，逐渐与弯管凸面相切。上述变化说明在弯管前驻点线附近产生了轴向流。竖直段后方流线仍然保持水平，但是方向却与来流相反，表明该范围内出现了回流区，竖直段附近的回流区大小随着雷诺数的增大逐渐减小。弯曲段流线与来流方向之间存在一定夹角，当100≤Re≤500时，在0°<Φ<42.5°范围内流线与x轴的正方向夹角大于90°，流线存在沿x轴负方向和z轴负方向的分量，说明既有回流现象又有轴向流。当42.5°<Φ<90°时，流线与x轴的正方向夹角小于90°，该区间只有轴向流，不会出现回流区。当1 000≤Re≤3 900时，回流区出现的范围为0°<Φ<44.5°，相比低雷诺数有所增大。高洋洋等[32]通过数值模拟研究了倾斜角度变化对圆柱流场特征的影响，结果表明雷诺数较小时，若倾斜角度大于15°，圆柱尾流中不会出现回流区，而在高雷诺数下，倾斜角度超过30°时回流区消失。Gallardo等[30]在Re=3 900时观察到回流区在Φ≥45°后基本不会出现。这和文中得到的高雷诺数下弯管尾流中回流区的分布规律相似。

图12 不同雷诺数下在xoz平面时均流线拓扑Fig. 12 Time-averaged streamline topology at xoz plane for different Reynolds numbers

为了深入探究尾流的展向分布特征，图13给出了不同z/D处〈u′u′〉/U2和〈v′v′〉/U2沿管体中心剖面位置(y/D=10)的分布。坐标轴中xb定义为22.5-(R-D/2)cosΦ，即对应z在基点线上的横坐标。从基点到〈u′u′〉/U2最大值位置之间的流向距离即为旋涡形成长度，随着z/D减小，〈u′u′〉/U2和〈v′v′〉/U2峰值减小，即展向涡脱强度降低，同时向下游移动，表明剪切层的发展更加稳定，旋涡形成长度变大，图7和8的瞬时展向涡量反映了相同的变化规律。当Re=100时，在z/D=11.0位置处没有出现〈u′u′〉/U2峰值，且〈v′v′〉/U2基本为0，由图7(e)可知管体两侧仅有剪切层存在，无旋涡脱落现象。如图8(e)所示，当Re=3 900时，在z/D=11.0位置处剪切层卷起产生旋涡脱落，〈u′u′〉/U2出现峰值，〈v′v′〉/U2大于0并出现峰值，但峰值范围相比于其他位置扩大，能量耗散严重，即使旋涡脱落，但涡脱强度仍然较低。随着雷诺数的增大，可以发现相同z/D位置处的〈v′v′〉/U2值增大，即雷诺数增大，涡脱强度相应增大，和图5中瞬时三维涡量的变化趋势相似。

图13 不同雷诺数下不同展向位置处的雷诺正应力Fig. 13 The normal Reynolds stress of the streamwise along the centreline of the wake at different depths

2.3 时均压力系数

图14～15为时均压力系数沿圆周和展向的分布情况，其中定义前驻点所在的位置处θ=0°，沿圆周顺时针方向为正。

图14 不同雷诺数下的圆周时均压力系数Fig. 14 Time-averaged pressure coefficients along cylinder surface for different Reynolds numbers

如图14所示，不同展向位置的Cp分布关于θ=180°(后驻点)对称。在z/D=27.5和22.5位置处的Cp值大小非常接近。在弯曲段上端z/D=17.5处，Cp值相比竖直段有所减小，但变化趋势仍与竖直段保持一致。在z/D=12.0处局部曲率较大，该处的Cp在θ=0°时明显小于其他展向位置的值，沿圆周的变化幅值也相应减小。Zhao等[11]对倾斜圆柱时均压力系数的研究结果表明，圆周Cp会随倾斜角度的增大而减小。研究发现分离点对雷诺数的变化较为敏感，当100≤Re≤500时，分离点θ≈85°，当1 000≤Re≤3 900时，分离点减小为θ≈75°。

从图15中可以看到，在不同雷诺数下，竖直段的前驻点时均压力系数Cps基本不会随展向位置而改变。Cps沿弯曲段逐渐减小，产生的压力差导致流向发生变化，因此时均流线在上游区域靠近弯管时形状发生弯曲。相比较前驻点时均压力系数，雷诺数变化对后驻点时均压力系数Cpb的影响更为明显，Cpb的绝对值随着雷诺数的增大而增加。另外，弯曲段上部(0°<Φ<45°)的Cpb小于下部(45°<Φ<90°)，因此沿弯管展向回流区内外的压力差减小，时均流线逆时针方向旋转，最终与来流方向基本平行。前后驻点时均压力系数在弯曲段和水平段交接处(s/D=19.6)二者相等，该位置处时均压力系数值约为-0.1，与Miliou等[13]、Gallardo等[30]的研究结果基本一致。

图15 不同雷诺数下时均压力系数沿展向分布情况Fig. 15 Variation of time-averaged pressure coefficients along the span of the curved cylinder for different Reynolds numbers

2.4 旋涡脱落频率

为获得弯管在不同展向位置处的旋涡脱落频率，对距离弯管轴线后方5D位置处的瞬时横流向速度进行快速傅里叶变换，开展频谱分析。图16表示100≤Re≤3 900的速度频谱，横坐标为无量纲的涡脱频率，纵坐标根据相应的最大谱峰值进行无因次处理。尽管在不同展向位置处都会出现明显的频谱峰值，但是随着z/D的减小，对应的峰值能量逐渐降低。当Re=300、500和1 000时，频谱峰值频率带随z/D的减小而变宽，这表明沿展向局部倾角增大导致涡脱强度衰减，影响了脱落旋涡的完整性，能量出现耗散，旋涡的紊乱无序导致谱峰频率范围扩大。当Re=1 500和3 900时由于尾流三维特性，频谱图中存在较多的谐波。

图16 不同雷诺数下不同展向位置处的速度频谱Fig. 16 Power spectrum of velocity at different depths for different Reynolds numbers

图17对比了不同雷诺数条件下弯管在不同展向位置处的无量纲涡脱频率St，可以看到弯管沿展向的St并非保持一致。当Re=100时，在z/D=27.5、22.5和17.5处的St相等，均为0.167，在z/D=12.0处St略有减小。Miliou[13]在该雷诺数下观察到了类似现象，即频谱能量沿展向衰减，但涡脱频率保持不变。当300≤Re≤3 900 时，展向St的变化规律较为相似，即从z/D=12.0到17.5变化逐渐增大，从z/D=17.5到27.5位置处St逐渐减小，在z/D=27.5处的St和竖直管体St值基本保持一致，弯曲段和水平段的存在对竖直段涡脱频率的影响较小[30-31]。Norberg[33]通过总结竖直圆柱St与Re的关系发现，当Re≤300时，Re的变化对St的影响较大，而当Re>300时，在亚临界范围内，Re对St的影响减弱。如图17所示，当100≤Re≤300时，弯管St的范围为0.166～0.203，当1 000≤Re≤3 900时，弯管的St基本在0.21左右，说明在高雷诺数范围内，St对Re的变化并不敏感。由于Re=300正好处在St和Re关系的临界点，其展向St的变化相对其他雷诺数变化波动更大。