0 引　言

我国近岛礁自然地理条件非常复杂^[1]，海底地质大多为沙质，水下地形相对平缓，吃水较浅，且水深变化剧烈，船舶吃水上靠岸困难^[2]，而浮桥能很好地解决这一难题。按照我国交通运输卸载装备发展的规划，未来海上输送卸载系列装备将共同形成一个海岸卸载系统，主要用于近岸及岸滩区域实施大规模装备、物资卸载和人员换乘，完成多种船舶的卸载和越滩上岸任务，浮桥在海岸卸载系统中将发挥重要作用。当波浪从深海向岛礁附近传播时，由于复杂的地形条件，波浪会发生一系列变化^{[3 − 4]}，对浮桥的安全性造成了严重威胁。除了波浪荷载外，还应考虑波浪引起的运动响应和结构形变，所以研究浮桥的运动响应和结构形变对浮桥的稳定运作具有重要意义。

目前，国内外对浮桥的水动力性能主要研究方法有理论分析、模型试验、数值模拟等方法。Seif等^[5]基于三维分布源法，计算了不同波高、周期及入射角度的波浪作用下浮桥的动力响应问题。Ikegami等^[6]在考虑浮桥弹性变形的基础上，将波浪作用直接转化为恒定的波浪力施加在浮桥结构上，然后分析了波浪周期、水深等因素对浮桥位移、加速度及变形的影响。Chen等^[7]对一箱型桥节进行了水池试验和数值计算分析，研究了波浪参数等对桥节运动特性和系泊力的影响。郭飞等^[8]基于水弹性方法，对靠岸工况、旁靠工况和不规则海底工况浮桥的位移响应进行数值预报，得到了海岸、运输船和海底的存在对浮桥的影响规律。何良德等^[9]基于三维势流理论和刚体运动方程建立了浮箱水动力响应数值模型，分析了面元及莫里森单元模拟减摇板影响的适用性。

可见，在处理外部波浪载荷问题时，一般将波浪载荷简化为作用在浮桥表面上的激励力，很少考虑流体与浮桥结构的耦合作用。另外，在少数考虑流体和浮桥结构的耦合作用的文献中，采用的也多是势流理论，没有考虑流体的粘性。基于此，本文采用双向流固耦合的CFD-FEA数值方法，在粘性流体-柔性结构耦合作用下，考虑浮桥结构形变的影响，将浮桥看作梁单元，对规则波工况中的气囊式承重浮桥模型的运动响应和垂直弯矩进行预报。

1 数值方法和模型设置

对气囊式承重浮桥进行建模、仿真和数据分析。采用双向流固耦合的方法，将CFD求解器STAR-CCM+和FEA求解器Abaqus以迭代交错耦合方式相结合。

1.1 浮桥模型描述

以气囊式承重浮桥1∶1模型^[10]为研究对象，浮桥模块由上部结构和气囊组成。浮桥模型的三维视图如图1所示。浮桥模型的主要参数见表1。浮桥的运作稳定性和工作强度按照南海复杂的海况来设计，环境载荷模拟也满足浮桥使用和规范要求。

1.2 流体动力学求解器

采用有限体积法（FVM）模拟流体区域，利用流体体积法（VOF）捕捉自由表面。通过求解非定常雷诺平均Navier-Stokes（URANS）方程来表示流体流动。根据半隐式压力关联方程（SIMPLE）型算法获得整体求解过程，以实现压力和速度之间的隐式耦合。由于Realizable k-ε湍流模型在浮体耐波性研究中得到了广泛应用，本文采用Realizable k-ε湍流模型。

在RANS-VOF求解器中，粘性流体被简为不可压缩流动，因此，连续性和动量方程如下：

$ \frac{\partial_{\rho}}{\partial_{t}}+\frac{\partial}{x_{i}}\left(p u_{i}\right)=0 ，$

(1)

$ \frac{{\partial \left( {\rho {u_i}} \right)}}{{\partial t}} + \frac{\partial }{{{x_j}}}\left( {\rho {u_i}{u_j}} \right) = - \frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{{x_j}}}\left( {\mu \frac{{\partial {u_i}}}{{\partial {x_j}}} - \rho \overline {u_i^\prime u_j^\prime } } \right) + {S_i} 。$

(2)

式中：$ i $和$ j $的取值范围是（1, 2, 3）；$ \rho $为流体密度；$ u $为动力粘性系数；$ {u_i} $和$ {u_j} $为速度分量时均值；$ p $为压力时均值；$ {S_i} $为动量方程广义源项。

1.2.1 流体域和边界条件

在STAR-CCM+软件的帮助下，建立了数值波浪水槽。坐标系的原点位于浮桥迎接波浪面与静水自由面交界的中心点处。浮桥的纵向方向与x轴对齐，侧壁在y轴上延伸，并且空气在指向天空的z轴上延伸。一个计算流体域与浮桥模型和符号选定边界条件的一般视图描绘如图2所示。采用重叠网格技术模拟浮桥在波浪中的大幅运动。流体域的边界条件：速度入口应用于流体域的垂直边界（包括波浪入口和侧壁）以及底部和顶部，以便根据需要定义这些边界条件下的流体速度；在边界出口处应用压力出口。

1.2.2 网格生成

使用STAR-CCM+中的自动网格划分功能执行网格生成。采用切割体网格法生产高质量的网格。网格细化是围绕浮桥周围、浮桥表面、自由液面进行的，采用面控制法和体控制法来实现这些区域中的细化网格密度。图3为充气式承重浮桥周围的精细网格。图4为充气式承重浮桥浮桥表面的精细网格。浮桥表面周围网格上的压力和壁剪切应力在每个求解时间步长提取，并将其应用于FEA模型。图5显示了CFD求解器中浮桥运动的垂直位移幅度示例。

1.2.3 波浪模拟

准确模拟入射波浪是研究浮桥波浪响应的基础。本研究使用五阶Stokes波，与线性Airy波相比，它更接近于真实的的非线性波。为了防止波在传播过程中的衰减和波从边界反射，选择使用EOM来解决这个问题，并获得了良好的结果。

为了监测模拟的入射波，在波浪水槽中布置了5个波浪探头的阵列。G₁探头的位置为原点沿z轴方向与水面的交点，每隔1 m再设置一个探头。图6给出了在波浪工况（ T = 1.75 s ）下，在G₁的监测值和目标值之间的波高比较示例。结果表明，基于EOM的造波方法能够较好地再现入射波，误差小于2%。

1.3 结构求解器

浮桥的运动和接头的结构变形是在Abaqus软件的帮助下通过FEA求解器计算的。Abaqus是一种基于位移的方法，有限元解基于虚功或虚位移原理^[11]。使用隐式动力分析求解器来求解控制方程。

1.3.1 结构模型

为了与CFD求解器中的流体动力学模型一致，还在相同的模型比例下进行FEA计算。CFD和FEA模型在其坐标系中具有相同的尺寸和相同的位置。对于FSI建模的结构部件，目前文献中2种最常用的方法是梁模型和FEA（有限元分析）模型。梁模型在计算上是高效的，并且通常能给予合理的结果。因此浮桥上部结构的有限元模型包括无质量表面壳体和骨干梁，气囊采用刚体单元构建。通过运动耦合约束连接到梁节点的3D壳单元将外部流体载荷从CFD解决方案传输到1D梁单元，以进行结构动态分析。梁的振动变形也将反馈到表面壳体，用于后续的水动力分析。浮桥的上部结构每隔0.55 m设置一个梁节点，共11个节点，作为11个监测站，第一个节点为0号监测站，第11个节点为10号监测站。在相邻节点之间进行监视，获得各节梁中心处的截面力和弯矩。图7所示为在FEA解算器中某一时刻充气式承重浮桥运动时梁上分布的结构应力图。

1.3.2 约束和荷载

对有限元模型进行约束，以实现y轴（无摇摆）平移以及绕x轴和z轴的旋转（无滚转和偏航）。此外，有限元模型还应限制在x轴方向，以防止波浪引起的纵向漂移。

1.4 双向耦合方法

本文采用双向耦合法，Abaqus中的位移数据通过网格映射技术传递到STAR-CCM+里面，STAR-CCM+中压力场发生变化，通过迭代，压力达到新的平衡，每个体积单元中心的压力和壁剪切应力传回到Abaqus里面，Abaqus在新的压力载荷作用下产生新的位移，新的位移数据再传递到STAR-CCM+里面去，实现了一个双向的、循序渐进相互耦合的过程，以获得更精确的结果。在有限元求解器中求解整体浮桥运动和结构变形。

双向耦合可以进一步分为显式耦合和隐式耦合。在本研究中，CFD和FEA之间的双向交互模拟通过隐式耦合算法进行。在隐式耦合中，流体载荷和结构变形在每个时间步交换3次，每个时间步长内的迭代次数为12，能够准确再现在波浪作用下充气式承重浮桥的运动响应和垂直弯矩。

1.5 计算条件

表2总结了本研究中涉及的调查条件，波浪工况满足浮桥实际作业的使用要求。选择0.001 s的恒定时间步长进行双向协同模拟。协同模拟由具有32个核心运行内存128 GB的工作站以并行模式运行，在一般情况下，需要大约15天的时间才能获得浮桥模型响应的18 s物理时程。

2 网格无关性验证

众所周知，由CFD求解器对流体动力学建模和模拟引起的不确定性通常比由FEA求解器与结构响应相关联的不确定性大得多。因此，目前的不确定性研究只涉及浮桥的运动响应，以第一个工况下（T = 1.75 s）的浮桥运动响应用于不确定性分析，有限元网格尺寸不变，根据不同的流体网格尺寸，共设计了3套方案。粗糙网格模型网格总数为83万，中等网格模型网格总数为97万，精细网格模型网格总数为129万。在不同网格密度的数值模型中运行联合仿真计算，获得波浪工况下浮桥运动响应的时程曲线如图8所示。可以看出，粗糙网格模型的计算结果与中等、精细网格模型差异较大，而中等、精细网格模型的计算结果十分接近。故综合考虑数值模拟结果准确性以及计算成本，本文选择中等网格划分流体域模型。

3 全局运动

充气式承重浮桥在不同波浪周期规则波工况中的纵摇和垂荡运动时间序列的综合比较如图9所示。在不同波浪周期的规则波工况中，纵摇、垂荡的峰值和幅值都是从时间序列中识别出来的，总结在表3中。可知，在规则波工况中，浮桥的纵摇运动的幅值随着波浪周期的递增表现先增加后减小的趋势；在$ T = 1.75\;{\mathrm{s}} $时，浮桥的纵摇运动幅值最大，最大值为10.1°；在规则波工况中，浮桥的垂荡运动幅值随着波浪周期的递增表现增大的趋势，在$ T = 2.78\; {\mathrm{s}} $时，浮桥的垂荡运动幅值最大，最大值为0.32 m。

4 垂直弯矩 4.1 不同波浪工况下的垂直弯矩

充气式承重浮桥在不同波浪周期规则波工况中的中剖面垂直弯矩时间序列的综合比较如图10所示。在不同波浪周期的规则波工况中，垂直弯矩的最大值如图11所示。可知，在规则波工况中，充气式承重浮桥的垂直弯矩随着波浪周期的递增表现先增加后减小的趋势；在T = 1.75时，充气式承重浮桥的中剖面垂直弯矩最大，最大值为48.8 kN·m。

4.2 沿桥长不同位置处的垂直弯矩

表4为浮桥在沿桥长不同位置处最大的垂直弯矩。可知，充气式承重浮桥在同一波浪工况下，浮桥中间所受垂直弯矩最大，最易受到损坏，应受到重视。

5 结　语

采用CFD-FEA双向耦合流固耦合方法，对充气式承重浮桥的运动响应和垂直弯矩进行了数值模拟。系统分析了充气式承重浮桥在不同波浪周期波浪条件下的计算结果。得到主要结论如下:

1）在规则波工况中，充气式承重浮桥的纵摇运动幅值随着波浪周期的递增表现先增加后减小的趋势；

2）在规则波工况中，充气式承重浮桥的垂荡运动幅值随着波浪周期的递增表现增大的趋势；

3）充气式承重浮桥中剖面的垂直弯矩随波浪周期的增加呈先增大后减小的趋势；

4）通过与充气式承重浮桥浮桥的联合仿真，验证了所开发的耦合方法能够模拟波浪载荷下浮桥运动响应和结构变形范围内的大部分物理特性。