0 引　言

海洋拖曳系统是一种流场和结构变形高度耦合的动力学系统，拖曳缆绕流场中存在的流动特征决定了拖曳系统的动力学响应行为。海洋拖曳缆的雷诺数通常在10⁴～10⁵之间，其尾流性质介于椭圆绕流和圆柱绕流之间。Kundu^[1]指出，在Re范围内，圆柱体或椭圆的尾流没有明显的规则尾流结构，但在圆柱体前缘有一个相对稳定的边界层分离位置。已经通过多相插值和变换来建立流场和电缆动力学之间的联系。Du等^[2]通过计算水下航行器及其螺旋桨流场中的阻力，获得了集总拖曳缆阵列上的阻力，并采用CFD方法进行了计算。Yang等^[3]认为，拖缆上的流体动力学是由连接缆段节点处流场的实时速度决定的。

王飞等^[4]指出缆形取决于法向阻力，拖缆张力与切向阻力平衡。李光明等^[8]采用有限体积法求解拖缆绕流场获得拖曳缆的水动力系数，发现水动力系数的非定常性。章浩燕等^[6]将拖曳缆的阻力近似为柱体在定常流中所受阻力，使用莫里森公式计算拖曳缆的流体阻力。

Zhao^[7]、Guan^[5]、Wu^[9]等通过调整拖曳缆的缆段数量和水流阻力系数等模型参数，更改拖曳缆动力学离散求解格式，结合海试数据修正拖曳缆的动力学模型。将缆段的阻力系数看成常数，造成相同拖曳工况的拖缆张力和缆形与海试数据差异较大。拖曳缆受到水流阻力的计算精度决定了拖曳系统动力学建模计算的精度。Grosenbaugh^[10]、Wang^[11]、Kishore^[12]、Gao^[13]等采用恒定阻力系数计算拖曳缆缆段的水流作用力，对拖曳系统的回转运动进行系统的模拟，不能准确模拟回转操纵中拖曳缆绕流场中存在大范围分离流动形成的流体阻力突降或陡升。

为了改进拖曳缆绕水流阻力预报的精度，更为准确地模拟拖曳缆的水下动力学响应行为，本文将对节点位置有限元法（NPFEM）进行修改，利用计算流体力学方法建模精细计算对应拖曳缆的缆形绕流场，获得拖曳缆的绕流场结构特点和水流阻力，将水流阻力信息传递给节点位置有限元计算模块，提升拖曳系统动力学建模计算的精度。

1 流固耦合方法的构造 1.1 缆的动力学建模

拖曳缆位于水面端点处为坐标原点，水深方向为Z轴负方向，拖曳方向为全局坐标系的X轴方向，采用右手坐标系确定Y轴指向，建立全局坐标系O-XYZ。局部坐标系o-xyz 的x轴沿缆的切线方向，z轴位于缆的法平面内，局部坐标系也服从右手坐标系。局部坐标系和全局坐标系采用欧拉旋转矩阵进行转换。

基于Zhu^[14]和Sun^[15]采用虚功原理推导杆单元模拟拖曳缆段的方法。将拖曳缆离散为一组坐标为(X_i，Y_i，Z_i) (i=1, 2, …, N)节点，节点之间采用杆单元模拟缆段。

根据节点位置有限元法构造的拖曳缆动力学方程为：

$ \underbrace{{\left[\boldsymbol{M} + {\boldsymbol{M}}_{a}\right]}}_{1}{\ddot{\boldsymbol{X}}}_{e} + \underbrace{\boldsymbol{C}}_{2}{\dot{\boldsymbol{X}}}_{e} + \underbrace{\boldsymbol{K}}_{3}{\boldsymbol{X}}_{e} = \underbrace{{\boldsymbol{F}}_{K} + {\boldsymbol{F}}_{a} + {\boldsymbol{F}}_{d} + {\boldsymbol{F}}_{{bg}}}_4。$

(1)

式中：$ {\boldsymbol{X}}_{e} $为节点的全局坐标向量；$ {\dot{\boldsymbol{X}}}_{{e}} $和$ {\ddot{\boldsymbol{X}}}_{{e}} $为时间的一阶和二阶导数。

项1中$ \boldsymbol{M} $为缆元的质量矩阵；$ {\boldsymbol{M}}_{a} $为缆元周围流体产生的附加质量矩阵

$ \boldsymbol{M}=\frac{\rho AL}{6}\left[\begin{array}{cc}2{\boldsymbol{I}}_{3\times 3}& {\boldsymbol{I}}_{3\times 3}\\ {\boldsymbol{I}}_{3\times 3}& 2{\boldsymbol{I}}_{3\times 3}\end{array}\right]。$

(2)

$ {\boldsymbol{M}}_{a}=\frac{1}{6}{C}_{m}{\rho }_{0}AL\left({\boldsymbol{M}}_{a0}-{\boldsymbol{M}}_{a1}\right)，$

(3)

$ {\boldsymbol{M}}_{a0}=2{\mathit{I}}_{6\times 6} \text{，} {\boldsymbol{M}}_{a1}=\left[\begin{array}{cc}2{\mathit{m}}_{0}& {\mathit{m}}_{0}\\ {\mathit{m}}_{0}& 2{\mathit{m}}_{0}\end{array}\right]，$

(4)

$ {\mathit{m}}_{0}=\left[\begin{array}{ccc}{\mathrm{cos}}^{2}{\theta }_{X}& \mathrm{cos}{\theta }_{X}\mathrm{cos}{\theta }_{Y}& \mathrm{cos}{\theta }_{X}\mathrm{cos}{\theta }_{Z}\\ \mathrm{cos}{\theta }_{X}\mathrm{cos}{\theta }_{Y}& {\mathrm{cos}}^{2}{\theta }_{Y}& \mathrm{cos}{\theta }_{Y}\mathrm{cos}{\theta }_{Z}\\ \mathrm{cos}{\theta }_{X}\mathrm{cos}{\theta }_{Z}& \mathrm{cos}{\theta }_{Y}\mathrm{cos}{\theta }_{Z}& {\mathrm{cos}}^{2}{\theta }_{Z}\end{array}\right]。$

(5)

式中：$ \rho $为材料密度；$ A $为缆元的截面面积；$ {\mathit{I}}_{3\times 3}、{\mathit{I}}_{6\times 6} $为单位矩阵；$ {C}_{m} $为拖缆的附加质量系数；$ {\rho }_{0} $为流体密度，$ {\theta }_{X} $、$ {\theta }_{Y} $、$ {\theta }_{Z} $为缆单元在全局坐标系中的空间方位角。

项2为阻尼矩阵，本研究采用瑞利阻尼模型计算阻尼矩阵：

$ \boldsymbol{C}=\beta \left[\mathit{M}+{\mathit{M}}_{a}\right]+\gamma \mathit{K}。$

(6)

式中：$ \beta $和$ \gamma $为瑞利阻尼系数。

项3为缆元刚度矩阵：

$ \boldsymbol{K}={\boldsymbol{Q}}^{\mathrm{T}}{\mathit{K}}_{0}\boldsymbol{Q}，$

(7)

$ {\boldsymbol{K}}_{0}=\frac{EAL}{{L}_{0}^{2}}\left[\begin{array}{cccccc}1& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ -1& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]。$

(8)

式中：$ E $为缆元的杨氏模量；$ L $和$ {L}_{0} $为变形和未变形的缆元长度；$ \boldsymbol{Q} $为局部坐标系向全局坐标系转换的坐标转换矩阵，$ \boldsymbol{Q} $为一个6×6的矩阵。其中非零元素为$ {\boldsymbol{Q}}_{\mathrm{1,1}}= {\boldsymbol{Q}}_{\mathrm{4,4}}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{X} $，$ {\boldsymbol{Q}}_{\mathrm{1,2}}={\boldsymbol{Q}}_{\mathrm{4,5}}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{Y} $，$ {\boldsymbol{Q}}_{\mathrm{1,3}}={\boldsymbol{Q}}_{\mathrm{4,6}}=\mathrm{c}\mathrm{o}\mathrm{s}{\theta }_{Z} $。

项4为缆元受到的各种力向量；$ {\boldsymbol{F}}_{{K}} $为缆元弹性产生的广义节点力向量；$ {\boldsymbol{F}}_{{a}} $为由缆元周围流体的附加质量产生的惯性力向量；$ {\boldsymbol{F}}_{{d}} $为阻力向量；$ {\boldsymbol{F}}_{{bg}} $为浮力和重力向量。

式(1)离散成关于变量(X_i，Y_i，Z_i)的二次方程组，计算规模取决于离散缆段的数量N，数值包IMSL^[16]来求解该大型二次方程组，采用Newmark计算格式求解时间步长。

1.2 直航拖曳缆绕流场计算建模

精确模拟拖曳缆绕流场是准确计算拖曳缆受到的水动力的基础。拖曳缆近似简化为圆形横截面且水力光滑的表面。因拖曳缆的缆径雷诺数在1×10⁴～4×10⁵之间，拖曳缆表面的绕流中普遍存在边界层分离流动现象。应选取能够准确预报拖曳缆分离流动现象的绕流场计算建模技术。

拖曳缆的绕流边界层是流体流过固体表面时形成的一个薄层区域内，速度从零增加到外部流动的速度，速度梯度很大，因此粘性力在边界层内起着重要作用。拖曳缆的绕流边界层可以分为粘性底层、过渡层和完全湍流层等层域，可根据无量纲数y+的大小来确定层域。

为了准确模拟边界层内存在的分离流动，基于雷诺平均方程，采用求解至层流底层的SST$ k-\omega $湍流模型求解包括粘性底层、过渡层和完全湍流层的拖曳缆绕流场。该模型要求近壁面网格高度值满足y+<1，才能准确描述拖曳缆的绕流场特征。为此设计一种确保网格质量的结构化网格生成算法，该算法以绕缆的圆形截面C-H型结构化网格作为网格源面，源面设计中靠近圆形壁面区域确保y+<1的同时确保网格局部密度合理布局，采用节点位置有限元法计算得到的初始缆形作为引线，驱动源面沿引线拉伸，生成三维结构化网格。流场的边界包括速度入口，压力出口和壁面。

使用如式(9)的多项式拟合时刻为t时的拖曳缆空间缆形$ {l}_{t} $，为确保缆形的拟合精度，实际采用N=3的三阶多项式拟合缆形。

$ {l}_{t}\left(x,y,z\right)={C}_{0}+\sum _{n=1}^{N}{a}_{n}{x}^{n}+{b}_{n}{y}^{n}+{c}_{n}{z}^{n}。$

(9)

式中，系数可以使用拟合线性回归方法得到并重写为$ {{\boldsymbol{A}}}=\left[a_{1} \quad a_{2}\quad a_{3}\right]^{{\mathrm{T}}} $，$ {{\boldsymbol{B}}}=\left[{{b}}_{1}\quad {{b}}_{2}\quad {{b}}_{3}\right]^{{\mathrm{T}}} $，$ {{\boldsymbol{C}}}= [{{c}}_{1}\quad{{c}}_{2} {{c}}_{3}]^{{\mathrm{T}}} $，拟合关系式可表述为式：

$ {l}_{t}\left(x,y,z\right)={C}_{0}+\left[\mathit{A},\mathit{B},\mathit{C}\right]\left[\begin{array}{ccc}x& {x}^{2}& {x}^{3}\\ y& {y}^{2}& {y}^{3}\\ z& {z}^{2}& {z}^{3}\end{array}\right] 。$

(10)

使用Gmsh^[17]的二次开发环境设计网格自动生成程序。而后利用Fluent软件脚本读入网格并自动设置边界条件求解雷诺平均方程和SST$ k-\omega $湍流模型，获得每个时刻对应的拖曳缆形对应的阻力$ {\mathit{F}}_{{d}} $，代入式(1)中完成当前时刻的缆形更新求解。直至相邻2个时刻对应的缆形累积误差$ \varepsilon =\left|{l}_{{t}_{i+1}}-{l}_{{t}_{i}}\right|/{l}_{{t}_{i}} $，满足设置的误差要求。

1.3 拖曳缆流固弱耦合模式构造

由定常计算得到全局坐标系下的水流阻力$ {\boldsymbol{F}}_{{d}}= [{F}_{dx},{F}_{dy},{F}_{dz}] $沿着缆长也可拟合为多项式：

$ \left\{\begin{split}&{F}_{dx}={f}_{1}\left({l}_{t}\right)，\\ &{F}_{dy}={f}_{2}\left({l}_{t}\right)，\\ &{F}_{dz}={f}_{3}\left({l}_{t}\right)。\end{split}\right. $

(11)

将拖曳缆的节点坐标代入多项式(11)插值计算出节点位置有限元法模块的离散缆段对应流体阻力值。

$ \left\{\begin{split}&{C}_{x}={F}_{dx}/\left(\frac{1}{2}\rho {V}^{2}{L}_{s}d\right)，\\ &{C}_{y}={F}_{dy}/\left(\frac{1}{2}\rho {V}^{2}{L}_{s}d\right)，\\ &{C}_{z}={F}_{dz}/\left(\frac{1}{2}\rho {V}^{2}{L}_{s}d\right)。\end{split}\right.$

(12)

上述耦合模式的构造过程如图2所示。流场和结构场的数据流向表明只存在流场向拖曳缆的受力数据传递，但没有考虑拖曳缆对流场施加的作用力，因次本研究将该耦合模式称为拖曳缆流固弱耦合模式。

2 拖曳缆动力学流固耦合模式验证 2.1 拖曳系统参数与网格生成

采用Guan^[5]的拖曳实验验证本研究建立的拖曳缆动力学流固耦合方法。拖曳系统的结构参数见表1所示。该拖曳系统采用的拖曳缆长较小，缆径雷诺数为3×10³～1.5×10⁴，拖曳缆为柔性细长体，不考虑缆的弯曲刚度和扭转刚度，为了获得较为准确的初始缆形，表1为一组较为接近真实值的拖曳缆的阻力系数估计值。

图3为拖曳缆周围流场的结构设计以及边界条件，图4为绕缆截面的C-H型拓扑结构。$ \mathrm{\lambda } $为缆的圆形截面，坐标原点位于点A₀。A₀B₁、A₀C₁、A₀C₂和A₀C₃的长度由缆的半径决定。式(13)和式(14)可以用于边界层网格细化。使用y+计算靠近壁面的第一层网格高度。各层网格高度比例为1∶1.05，确保了近壁面区域的网格密度。

$ {d}_{s}=R\sqrt{80}{y}^+{Re}^{-13/14}，$

(13)

$ {d}_{s2}=2{d}_{s}。$

(14)

在图4中，从A₁～B₁的近壁面第一层网格高度用$ {d}_{s} $表示，从B₁～C₁的近壁面第一层网格高度用$ {d}_{s2} $表示。三维网格是通过源面沿缆形曲线拉伸生成，该过程使用Gmsh^[17]二次开发来实现。

2.2 流固耦合模式计算结果

如图5所示，经过流场与缆索结构力学模块的多次迭代计算，缆形趋于稳定，小于累积误差$ \varepsilon $可得到最终的缆形结果。图5显示该缆形的流固计算值与实验值吻合较好。该流固耦合计算模式通常采用较少迭代次数即可达到理想的精度要求，在给定较为接近真实缆形的初始缆形分布时，拖曳缆流固耦合计算收敛速度较快，具有较好的收敛性。

随着拖曳雷诺数的增长，拖曳缆的绕流逐步由层流状态转换为强湍流，拖曳缆段受到水流作用力也趋于显著。图6为不同雷诺数的拖曳缆缆形迭代计算结果，经过数次迭代计算后，缆形较于采用恒定的阻力系数均发生了变化且趋于稳定，随着雷诺数的不断增加，采用恒定的阻力系数模拟的结果与该流固耦合计算模式获得的结果差距在逐渐扩大，该流固耦合计算模式能很好地适应增加的雷诺数。

沿缆的载荷分布如图7～图10所示，其中Z向阻力具有显著的波动，因为椭圆绕流的尾流会导致高升力，而圆柱绕流的尾流具有交替的剪切层脱落，这产生了不显著的升力和高阻力。

图11为拖曳缆截面的动压力分布情况。Z/d用于描述拖缆所处的水深位置，表示沿Z轴的拖缆横截面椭圆绕流向圆柱绕流的过渡。在低速拖曳下，可明显观察到交替的漩涡脱落。在高速拖曳下，拖缆的变形和高频振动对尾流结构的分布有很大影响。

图12为2种拖曳速度下涡度为6 s⁻¹的等值面。在0.3 m/s的拖曳速度下，涡度浓度沿拖缆的倾角分布，尾流漩涡脱落的空间频率表现为近表面端的椭圆尾流和低速拖曳时的圆柱尾流结构。对于1.5 m/s的高速拖曳，拖缆周围的涡度浓度更高。

大的倾斜角决定了椭圆长短轴的比例，当比例达到1时，拖缆周围的流动类似于圆柱绕流。拖缆的尾流可分为2个区域以及其之间的过渡区域，过渡区具有叠加的尾流结构，在0.3 m/s的尾流中产生了涡流融合，在高速拖曳下，剪切层脱落向前推进，并将这些涡流破碎成远尾流。沿缆的不同位置处的流动情况存在较大差异，采用恒定的阻力系数模拟就会带来较大误差。

图13和图14为不同拖曳速度下水深、分离点之间的关系。拖缆周围有2个基本对称的分离点。拖曳缆在不同水深处的绕流边界层影响着分离点的变化，这也说明沿着缆长的拖曳缆绕流场存在着较大差异，采用恒定的阻力系数并不能很好地模拟拖曳缆的绕流场情况。

3 结　语

本文建立一种用于拖曳系统的流固耦合计算模式。利用流体动力学经验系数代入节点位置有限元方程组计算得到初始缆形，围绕缆形建立流体计算域采用求解雷诺平均方程的方法，求解拖曳缆绕流获得拖曳缆流体阻力传递给节点位置有限元方程组，迭代计算得到满足误差要求的拖曳缆的缆形。该方法经过已有拖曳试验验证具有良好的精度，通过分析拖曳缆的绕流场给出了拖曳缆的绕流场中边界层分离流动特征，结论如下：

1）单向耦合拖曳系统的流固耦合计算模式，对直航拖曳具有良好的计算结果。

2）提出一种拖缆绕流场的结构化网格生成方法，该网格能准确描述拖曳缆的绕流场特征。

3）改进选择固定阻力系数进行数值模拟的方法，采用流固耦合计算模式更好地模拟了拖曳系统的真实水下运动情况。