深水域钢管桩平台与波浪流固耦合作用数值模拟

Numerical Simulation of Fluid-Structure Interaction between Steel Pipe Pile Platform and Wave in Deep Water

摘要: 本文采用k-ε湍流模型和VOF运动界面追踪法进行波浪数值模拟，通过位移–压力格式的有限单元法实现钢管桩平台与波浪的双向流固耦合作用，以长江口某桥墩施工平台为例，分析了钢管桩平台在波浪力作用下的受力和变形特征，为平台设计提供参考依据。

Abstract: In this paper, the k-ε turbulence model and VOF motion interface tracking method are used to simulate the numerical wave, and the two-way fluid-structure interaction between the steel pipe pile platform and the wave is implemented by using displacement-pressure format finite element method. The stress and deformation characteristics of the steel pipe pile platform under the action of wave force are analyzed by taking the construction platform at the Yangtze River Estuary as an example, which provides a reference for platform design.

1. 引言

近年来，人们在深水域环境中建造了大量结构物，如石油钻井平台、海上风电、港口码头及跨海、跨江大桥等。波浪对深水域建(构)筑物的造价、安全度及寿命起着重要作用 [1] 。实际工程中，小直径桩柱大多以群桩的形式存在，群桩间的相互作用使得结构物附近的波浪场更加复杂。因此，随着深水域工程施工的发展，与之相适应的波浪与结构物耦合作用的研究得到了广泛的关注，许多学者在波浪数值模拟以及波浪与桩柱的耦合作用方面开展研究工作。王安庆 [2] 采用Adina-STRUCTURE分析了流固耦合条件下波浪对风机桩基础受力性能的影响。周远洲等 [3] 采用CFD方法计算比较了桥梁基础在规则波作用下的各类波浪荷载计算方法的差异。朱晓玲 [4] 运用ANSYS CFX对单桩柱和波浪场进行了双向、单向流固耦合模拟计算，并对串列双桩柱和波浪场进行了双向流固耦合模拟。雷欣欣 [5] 结合运用物理模型试验和数值模拟的研究手段，开展五桩结构在波浪作用下的水动力特性研究。李志阳等 [6] 通过实验研究分析了波浪力作用下非对称群桩结构的群桩系数。何一宽等 [7] 采用半解析方法推导了承台–群桩结构波浪力的理论公式。

目前关于水域结构平台与波浪双向耦合作用的数值模拟研究较少，本文以带有空间夹角和不同桩间距的九桩平台结构为研究对象，建立了双向流固耦合模型，分析钢管桩平台与波浪之间的耦合作用，为工程实践提供参考。

2. 理论模型与控制方程

2.1. 流固耦合系统

假设流体为无粘性，可压缩和小扰动的，流体自由液面为小波动；固体为线弹性材料。线弹性结构与理想可压缩流体组成的流固耦合系统的有限元分析模型主要有两种：位移–位移格式和位移–势格式。由于前一种格式中流体剪切刚度为零，计算会出现零能模式，目前工程中主要采用后一种格式，其中固体域的方程采用位移u作为基本未知量，流体域的方程以流场压力p作为基本未知量；虽然有限元方程中的系数矩阵是非对称的，但没有零能模态的困难，相应的有限元表达格式称为流固耦合分析的位移–压力(u, p)格式。

2.1.1. 流体域

流体场方程为

$p_{,ii}-\frac{1}{c_0^2}\ddot{p}=0$ (1)

式中，p为流体压力，c₀为流体中的声速。

流体壁面边界条件和自由面边界条件分别为

壁面： $\frac{\partial p}{\partial n^{\left(f\right)}}=0$ (2)

自由面： $\frac{\partial p}{\partial z}+\frac{1}{g}\ddot{p}=0$ (3)

式中， $n^{\left(f\right)}$ 为流体壁面法线方向，g为重力加速度。

2.1.2. 固体域

固体场方程为

${\sigma }_{ij,j}+f_i={\rho }^{\left(s\right)}{\ddot{u}}_i$ (4)

式中， ${\sigma }_{ij}$ 为固体应力分量，u是固体位移分量， $f_i$ 为固体体积力分量， ${\rho }^{\left(s\right)}$ 为固体质量密度。

固体边界条件为

力边界： ${\sigma }_{ij}{n}_j^{\left(s\right)}={\bar{T}}_i$ (5)

位移边界： $u_i={\bar{u}}_i$ (6)

式中， ${\bar{T}}_i$ 和 ${\bar{u}}_i$ 分别是固体上的已知面力分量和已知位移分量。

2.1.3. 交界面协调条件

钢管平台与波浪相互作用的特征是耦合作用仅仅发生在流固交界面上，方程耦合是由流固耦合面的动力平衡及运动学协调关系引入的。

1) 运动学条件：流固交界面上法向速度保持连续，即

$v_n^{\left(f\right)}=v^{\left(f\right)}\cdot n^{\left(f\right)}=v^{\left(s\right)}\cdot n^{\left(f\right)}=-v^{\left(s\right)}\cdot n^{\left(s\right)}=v_n^{\left(s\right)}$ (7)

利用流体运动方程上式可改写为

$\frac{\partial p}{\partial n^{\left(f\right)}}+{\rho }^{\left(f\right)}\ddot{u}\cdot n^{\left(f\right)}=0$ (8)

式中， ${\rho }^{\left(f\right)}$ 为流体质量密度。

2) 动力学条件：流固交界面上法向力保持一致，即

${\sigma }_{ij}{n}_j^{\left(s\right)}={\tau }_{ij}{n}_j^{\left(f\right)}=-{\tau }_{ij}{n}_j^{\left(s\right)}$ (9)

式中， ${\tau }_{ij}$ 是流体应力张量的分量，对于无粘流体，

${\tau }_{ij}=-p{\delta }_{ij}$ (10)

上式代入式(9)，得到

${\sigma }_{ij}{n}_j^{\left(s\right)}=p{n}_j^{\left(s\right)}$ (11)

2.2. 湍流模型

考虑流体为湍流流体，需要额外补充湍流方程，k-ε湍流模型的脉动动能 满足方程：

$\frac{\partial \left(\rho k\right)}{\partial t}+\nabla \cdot \left(\rho uk\right)=\nabla \cdot \left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\nabla k\right]+P_k-\rho \epsilon +S_k$ (12)

式中，

$P_k={\mu }_t\nabla u\cdot \left(\nabla u+\nabla u^{\text{T}}\right)-\frac{2}{3}\nabla \cdot u\left(3{\mu }_t\nabla \cdot u+\rho k\right)$ (13)

${\mu }_t$ 是涡粘度，

${\mu }_t=C_{\mu }\rho \frac{k^2}{\epsilon }$ (14)

湍流耗散率 $\epsilon$ 满足方程：

$\frac{\partial \left(\rho \epsilon \right)}{\partial t}+\nabla \cdot \left(\rho u\epsilon \right)=\nabla \cdot \left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\nabla \epsilon \right]+\frac{\epsilon }{k}\left(C_{\epsilon 1}{P}_k-C_{\epsilon 2}\rho \epsilon \right)+S_k$ (15)

k-ε湍流模型中， $S_k$ 和 $S_k$ 为自定义源项， $C_{\mu }$ 、 ${\sigma }_k$ 、 ${\sigma }_{\epsilon }$ 、 $C_{\epsilon 1}$ 、 $C_{\epsilon 2}$ 均为经验常数，取值可参考文献 [8] 。

2.3. 波浪理论

采用线性波理论模拟波浪，波面方程为

$\eta =\frac{H}{2}\cos \left(kx-\omega t\right)$ (16)

波浪运动的速度场方程为

$\{\begin{array}{l}u=\frac{kHg}{2\omega }\frac{\text{ch}k\left(h+z\right)}{\text{ch}kh}\cos \left(kx-\omega t\right)\\ w=\frac{kHg}{2\omega }\frac{\text{sh}k\left(h+z\right)}{\text{ch}kh}\sin \left(kx-\omega t\right)\end{array}$ (17)

压强分布如下

$p=-\rho gz-\rho g\frac{H}{2}\frac{\text{ch}k\left(h+z\right)}{\text{ch}kz}\cos \left(kx-\omega t\right)$ (18)

以上各式中， $\eta$ 是波面离平均水面的垂直距离， $k=\frac{2\pi }{\lambda }$ 为波数， $\lambda$ 为波长； $\omega$ 为波浪圆频率；H为波高；h为固定边界到静水面的距离。

2.4. VOF界面追踪法

VOF (Volume of Fluid)方法既能够保证质量的守恒，又能够在不引起数值耗散的情况下处理大变形自由液面，特别适用于三维计算，是发展较为成熟的波面追踪方法。

引入相体积分数函数 ${\varphi }_q$ ， ${\varphi }_q$ 定义为单元内第q相体积与单元体积的比值。 ${\varphi }_q=1$ 表示该单元内全部为 相物质， ${\varphi }_q=0$ 表示q相物质为空， $0<{\varphi }_q<1$ 则表示该单元为交界面单元。通过相体积分数 ${\varphi }_q$ 的计算，求出各网格单元的 ${\varphi }_q$ ，进而在整个计算域内构建出界面，实现对相间界面的追踪。

对于深水域的波浪问题，流体域只有气液两相物质，相控方程如下

$\frac{\partial {\varphi }_q}{\partial t}+\frac{\partial \left(u{\varphi }_q\right)}{\partial x}+\frac{\partial \left(v{\varphi }_q\right)}{\partial y}=0$ (19)

${\displaystyle \underset{q=1}{\overset{2}{\sum }}{\varphi }_q}=1$ (20)

3. 数值求解方法

上述流固耦合(FSI)系统的偏微分方程复杂程度很高，再加上流体域和固体域边界条件的复杂性，通常很难得到精确解，工程上常采用数值方法求解。在流固耦合的位移–压力格式有限元模型中，流体域接受固体域传递的节点位移，并以节点位移为边界条件，求解出以流场压力为未知变量的方程，然后将边界压力传递给固体域，固体域以节点位移为基本方程未知量，求解出节点位移，再将节点位移传递给流体域作为其求解的边界条件，依次循环，固体域和流体域分别根据更新后的边界条件进行分区迭代计算，直至收敛。

本文利用ANSYS Mechanical采用有限单元法求解钢管桩平台的CSD (Computational Solid Dynamics)方程，利用CFX采用有限体积法求解波浪的CFD (Computational Fluid Dynamics)方程，借助ANSYSMFX多场求解器交换钢管桩平台和波浪耦合面的数据，具体实现思路如图1所示。

Figure 1. Realization idea of two-way fluid-structure interaction

图1. 流固双向耦合实现思路

4. 波浪数值模拟验证

参考长江口水文水资源勘测局徐六泾水文站的观测数据设计波浪，波高0.7 m，波长16 m，水深25 m，最大流速3.86 m/s。设计水槽尺寸为100 m (长) × 14 m (宽) × 30 m (深)，边界设置如下：顶面为自由面，为避免回流，出口和入口均使用开敞边界条件(Opening)，前后侧面为对称面，水槽底部设为壁面。

将数值计算得到的波面与设计波面对比，进行波浪验证。图2是t = 2.4 s时的模拟波面和设计波面对比，二者基本吻合，最大相对误差为0.28%，平均误差为0.009%，满足计算要求。

Figure 2. Comparison between simulated wave surface and design wave surface

图2. 模拟波面与设计波面对比

5. 实例分析

长江口某特大型桥梁桥墩的施工平台主要用于放置变电站及材料，并作为钻孔灌注桩施工完毕后塔吊的退位支撑基础，平台沿波向纵长24~27 m，横宽12 m，共有9根φ1400 (δ = 16 mm)钢管桩桩基，其中有5根桩与波向形成空间夹角，且各桩的桩间距不同，属于复杂不规则平台。钢管桩通过上下两道平联联成整体，上道平联位于桩顶，上下平联高程间距3 m，桩顶(平台)比静水面高5 m。

波浪水槽尺寸为100 m (长) × 50 m (宽) × 25 m (深)，平台离水槽入口距离不小于一个波长，离水槽出口距离不小于三个波长。波浪水槽和平台模型见图3所示。

Figure 3. DM modeling diagram of platform and wave tank

图3. 平台和水槽DM建模图

波高受群桩影响较小，波面变化主要发生在钢管桩桩周，如图4所示，入射波和反射波叠加后，桩周的波面随着钢管桩所处的波浪位置不同发生规律性的变化，说明平台桩间距较大时，单根钢管桩所形成的反射波无相互干涉，平台群桩对波面的影响规律与单桩影响规律类似，波面影响的群桩效应不明显。

与波面不同，波速受群桩的影响较大，截取t = 5.0 s时的波面流速，见图5，每根钢管桩周围的流速变化大致规律类似，桩前为减速区，两侧为加速区，桩后为速度较小的尾涡区，但是，减速区、加速区和尾涡区不再像孤立单桩一样呈规则半圆形或椭圆形，而且，群桩对波浪速度的影响范围比单桩大，在桩附近和桩间都发生变化，影响区连成一片，甚至会影响平台前后远处的波速。波浪对波面以下流场速度的影响随着深度增加而快速降低，低于波面3 m的流体受波浪影响较小。

跟波面速度类似，群桩对波面压强的影响范围也较大，见图6，不仅涉及整个平台区域范围，而且在平台前后形成了明显的压强变化台阶，说明群桩效应显著。由于群桩的各桩间距及其与波浪的空间夹角不规则，平台区域内的波面压强在群桩间的变化规律较复杂，这也意味着平台在波浪作用下的受力特征比较复杂。

Figure 4. Wave surface near the platform (t = 5.0 s)

图4. 平台附近波面图(t = 5.0 s)

Figure 5. Flow velocity nephogram of wave surface (t = 5.0 s)

图5. 波面流速云图(t = 5.0 s)

Figure 6. Pressure nephogram of wave surface (t = 5.0 s)

图6. 波面压强云图(t = 5.0 s)

图7~14是每T/8时间间隔的平台应力云图，提取整理一个波浪周期内平台最大位移及其构件最大应力，见表1。

Figure 7. Stress nephogram of platform (t = 0)

图7. 平台应力云图(t = 0)

Figure 8. Stress nephogram of platform (t = T/8)

图8. 平台应力云图(t = T/8)

Figure 9. Stress nephogram of platform (t = T/4)

图9. 平台应力云图(t = T/4)

Figure 10. Stress nephogram of platform (t = 3T/8)

图10. 平台应力云图(t = 3T/8)

Figure 11. Stress nephogram of platform (t = T/2)

图11. 平台应力云图(t = T/2)

Figure 12. Stress nephogram of platform (t = 5T/8)

图12. 平台应力云图(t = 5T/8)

Figure 13. Stress nephogram of platform (t = 3T/4)

图13. 平台应力云图(t = 3T/4)

Figure 14. Stress nephogram of platform (t = 7T/8)

图14. 平台应力云图(t = 7T/8)

Table 1. Maximum displacement of platform and maximum stress of components

表1. 平台最大位移和构件最大应力

在波浪作用下，平台的最大位移、构件最大应力变化幅度较大，平台强度和刚度满足要求。平联构件最大应力发生在与钢管桩的连接端，钢管桩最大应力发生在底部嵌固端，计算结果表明，平台平联与钢管桩的连接点应力集中明显。

6. 结论

1) 采用VOF界面追踪法和k-ε湍流模型进行数值模拟，得到的数值波与设计波吻合较好。

2) 借助ANSYSMFX多场求解器可以较好地通过位移–压力格式实现双向流固耦合数值模拟，但对复杂平台，双向耦合计算成本较高。

3) 桩间距较大时，平台群桩对波面影响不大，但对波速影响较大，这会导致钢管桩上的波浪作用发生变化，不规则平台波浪力影响的群桩效应复杂。本文建立的流固耦合数值模型可以得到不规则平台的波浪作用效应，这对于分析和评价实际工程平台结构的工作性能非常重要。

4) 长江口深水区某九桩工作平台在波浪作用下，最大应力集中在钢管桩与平联的连接节点处，且以第一排桩和斜桩与平联连接处最为集中，平台设计需重视并加强节点连接构造。

基金项目

国家重点研发计划项目(2022YFC3803403)资助。