利用平滑粒子流体动力学方法研究波浪作用下水下浮式隧道的数字化孪生技术

《Applied Ocean Research》：Digital twinning of submerged floating tunnel under wave actions using smoothed particle hydrodynamics

　　邵家琪|胡浩|梁东方
剑桥大学工程系，英国剑桥CB2 1PZ

**摘要**
本文探讨了无网格平滑粒子流体动力学（SPH）方法在模拟受波浪作用影响的潜浮式隧道（SFT）水动力响应方面的工程应用。与基于网格的数值方法相比，SPH方法在处理复杂的自由表面和流固耦合问题时具有独特的优势。首先，该模型被应用于一个受锚固约束的二维圆形隧道，并通过改变SFT的浮力-重量比来研究其三自由度（3-DOF）运动。随后，SPH模型被进一步应用于更具挑战性的三维场景，其中SFT由单个或双圆形截面组成。在规则波和孤立波的影响下分析了结构响应。本研究的一个显著特点是，通过一个附加接口将SPH数据与商业软件Blender连接起来，增强了SPH模拟的可视化解释能力。通过完全实现数值模拟与复杂计算机可视化的集成，实现了SFT的逼真数字孪生，以展示SPH建模方法在实际工程应用中的潜力。

**1. 引言**
潜浮式隧道（SFT）是一种创新的水下基础设施，旨在实现安全高效的水上通行，例如穿越山谷、海峡或深水通道。与通常建在海床下的传统地下隧道不同，SFT悬挂在水下指定深度。这种设计利用浮力使隧道保持在水面以下，通常通过固定在海床上的缆绳或浮筒来支撑。SFT在极端水下环境中具有许多独特优势，特别是在传统桥梁或隧道难以建造的情况下。由于位于水面以下，SFT可以避免风浪、船只和极端天气的影响，同时保持足够的灵活性以跨越相当的距离。图1展示了根据Chen等人（2020年）和Fan等人（2022年）的研究结果制作的一个SFT原型示意图。

**下载：**下载高分辨率图片（273KB）
**下载：**下载全尺寸图片

**图1. 潜浮式隧道示意图（AI生成）**

到目前为止，现场还没有完全建造的潜浮式隧道。它们大多仍处于理论研究阶段，尽管一些国家已经提出了一些计划。其中一个代表性项目是挪威的E39高速公路。这条沿海道路横跨1000多公里，包含七个渡轮通道，总行驶时间为21小时。挪威公共道路管理局计划用各种类型的桥梁替换所有七个渡轮通道，以将行驶时间缩短至11小时。其中，Sognefjord是一个特别具有挑战性的深海峡，水深可达1308米，使得许多传统结构不适用。因此，他们考虑使用由水面上的浮筒支撑的双管道潜浮式隧道（Dunham等人，2016年）。

潜浮式隧道（SFT）的工程重要性引起了研究人员和从业者的极大兴趣。尽管在这个领域缺乏全面的知识，但仍有一些已建立的研究。理论上，Di Pilato等人（2008年）对SFT在地震和水动力激励下的三维动态响应进行了研究。原型模型是穿越意大利墨西拿海峡的隧道，总长度为4680米，最大深度为285米，假设其受到地震、波浪和流速的联合载荷作用。他们将锚索建模为细长杆，并对整个结构使用了特定的有限元方法。Chen等人（2018年）提出了一个简化的模型，用于分析在波浪-流速联合作用下的SFT管缆系统的振动和共振问题，其中控制方程是根据汉密尔顿原理求解的。Jin和Kim（2020年）研究了在波浪激励下的动态隧道-锚固-列车SFT系统，开发了一个时域耦合的流体-弹性模型。隧道和锚固线的运动方程基于有限元杆理论，列车则采用多刚体方法进行建模。此外，Jeong和Kim（2024年）提出了一种基于一些实际假设的分析方法来研究SFT的响应，假设自由端边界条件；而Wu等人（2025年）使用超弹性本构关系研究了SFT中半刚性和半柔性接头在各种波浪和载荷条件下的力行为。

此外，还有许多实验研究报道。例如，Chen等人（2023年）对双圆形SFT在波浪作用下的运动进行了实验室测试，测量了结构响应和锚固线力，并使用Mathieu不稳定性理论分析了SFT的动态行为。Yuan等人（2023年）研究了缆绳布局对SFT在波浪作用下的运动影响。测量是在一个缩比原型模型上进行的，收集了三种不同缆绳布局下的隧道运动、缆绳力和波高等数据。最近，Yang等人（2025年）对在复杂波浪作用下的长跨度SFT的三维流体-弹性响应进行了大规模实验研究。与之前的文献不同，本次实验系统研究了SFT端部的约束、沿隧道的弹性变形以及管缆系统之间的耦合。Chen等人（2026年）使用受限截断模型进行了实验室实验，以研究SFT在波浪作用下的动态运动，并比较了受限边界条件和自由边界条件下的结构响应和锚固力。

此外，文献中还有许多数值研究报道。Chen等人（2020年）使用重叠网格解决了SFT在表面波作用下的二维Navier-Stokes方程，以动态配置系统。该模型应用于中国广东和海南省之间的琼州海峡原型SFT的波浪-结构相互作用。Jeong等人（2022年）研究了使用组合垂直和倾斜缆绳来稳定隧道的不同SFT配置。为了评估，使用ABAQUS-AQUA有限元程序进行了水动力模拟，并检查了不同缆绳间距和倾斜角度的组合。Jin等人（2023年）使用完全耦合的时域离散模块梁方法研究了固定SFT的流体-弹性行为，其中流体在频域的流体动力学参数来自三维势理论，结构在时域使用多体Cummins方程和Euler-Bernoulli梁及Saint-Venant扭转的耦合刚度矩阵。通过将结果与OrcaFlex和ABAQUS两个商业软件进行比较，验证了集成模型的准确性。如上所述，尽管存在许多关于各种波浪-结构相互作用的成熟数值解，但大多数都是基于固定网格方法的。对于涉及自由表面大变形和移动边界的剧烈流动，这些基于网格的方法需要使用非常复杂的网格处理来保持计算精度。另一方面，SPH方法的无网格特性利用粒子的动态运动避免了这一问题，提供了自然准确的解决方案。

平滑粒子流体动力学（SPH）方法的最初应用记录在溃坝流和波浪破碎研究中。SPH非常适合这种类型的流体动力学，因为它能够自然处理自由表面、大变形和复杂的边界相互作用。Monaghan（1994年）首次提供了使用SPH方法模拟自由表面流动的基础见解，即溃坝洪水的传播。Colagrossi和Landrini（2003年）开发了一个更先进的SPH模型，该模型包括两种不同的流体相，如水-空气相和水-油相。Gómez-Gesteira等人（2010年）对SPH方法在各种自由表面流动应用中的进行了全面回顾。Jian等人（2015年）研究了大坝崩塌后上下游水界的复杂演变，通过SPH流体动力学方法评估了实际自然灾害。近年来，SPH在波浪-结构相互作用中的应用作为CFD建模工具的替代方案受到了越来越多的关注。例如，He等人（2023年）使用开源代码SPHysics（http://www.sphysics.org）研究了由多个浮体组成的复杂浮动结构的波浪诱导响应。他们通过自己的实验室实验验证了数值锚固和碰撞力。此外，还进行了关于波浪和结构参数的全面研究，以评估浮动平台的敏感性和性能。

然而，关于SPH在潜浮式隧道领域的应用的研究非常有限，近年来只有少数案例被记录下来。例如，Luo等人（2019年）研究了SFT在潜艇冲击载荷下的损坏情况。在他们的模型中，水使用SPH模拟，而锚索、SFT和潜艇则使用有限元方法（FEM）建模。Shao等人（2022年）结合了SPH和机器学习（ML）算法来评估缆绳型SFT支撑系统的损坏情况。在他们的模型中，SPH用于预测结构响应和水动力力，以分析SFT的振动，而基于ML的向量自回归移动平均模型（VARMA）和Mahalanobis距离（MD）算法用于识别损坏的严重程度和位置。此外，Yang等人（2023年）对悬挂在固定水平平台上的潜浮式隧道进行了密集的SPH计算。他们确定了几个关键参数，如多波条件、锚固配置和隧道的浸没深度，这些参数对SFT结构的水动力响应有影响。他们还发现，隧道-平台系统的自然周期在波浪冲击引起的隧道运动中起着关键作用。

本文旨在充分探索SPH建模方法在这种创新工程应用中的潜力。该研究将解决3D SPH计算和数值结果可视化演示方面的知识空白，以实现真实世界现象的完全数字化。为了实现上述目标，我们建议使用基于GPU的SPH求解器来管理计算负载，并使用高性能渲染工具Blender来可视化数值数据。

**2. 方法论和算法**
本节提供了数字孪生工程潜浮式隧道的关键方法论和算法，从而将数值计算与计算机可视化相结合。

**2.1. SPH的控制方程**
在基于SPH的方法中，控制方程通过分布在计算域中的大量粒子离散化，使用粒子相互作用模型来处理流体动力学方程中的每一项。SPH方法具有许多优势，例如有效且准确地跟踪自由表面和流固界面。

流体的质量和动量守恒方程构成了该模型的核心组成部分，表示为：
\[1\] \(\rho \frac{d\rho}{dt} + \nabla \cdot u = 0\)
\[2\] \(\frac{du}{dt} = -\frac{1}{\rho}\nabla P + g + \nu_0\nabla^2u + \frac{1}{\rho}\nabla \cdot \tau\)
其中，\(\rho\) = 流体密度；\(t =\) 时间；\(u =\) 速度；\(P =\) 压力；\(g =\) 重力加速度；\(\nu_0 =\) 运动粘度；\(\tau =\) 额外剪切应力，包括非牛顿或湍流效应。在SPH框架中，\(\tau\) 通常通过Gotoh等人（2001年）提出的涡粘性公式来建模，表示为：
\[3\] \(\tau_{ij}/\rho = 2\nu_T S_{ij} - 2\frac{3}{k}\delta_{ij}\)
其中，\(\nu_T =\) 湍流或非牛顿流的粘度；\(S_{ij} =\) 平均流动的应变率；\(k =\) 动能；\(\delta_{ij} =\) 克罗内克德尔塔。在湍流中，\(\nu_T\) 通常通过众所周知的Smagorinsky模型计算，表示为：
\[4\] \(\nu_T = (C_s\Delta X)^2|S|\)
其中，\(C_s =\) Smagorinsky常数；\(\Delta X =\) 表示空间分辨率的粒子间距；\(|S| = (\frac{2S_{ij}}{S_{ij})^{1/2}\) 是综合应变率。方程（3）和（4）应与SPH主方程（1）和（2）耦合，称为亚粒子尺度（SPS）湍流模型（Gotoh等人，2001年），并在SPH社区中得到广泛应用。

为了计算方程（2）中的流体粒子压力，我们采用了弱可压缩SPH（WCSPH）方法的概念。它通过状态方程明确解决了压力问题，因此即使在工程模拟中涉及大量粒子，也能实现高度的计算效率。在这种解决方法中，压力是通过粒子的密度变化来计算的，使用以下公式（5）：\[P = B\left(\frac{\rho}{\rho_0}\right)^{\gamma - 1}\]，其中 \(\rho_0\) 是流体的参考密度（例如水的密度为 1000.0 kg/m3）；\(\rho\) 是中间时间步的密度；\(\gamma\) 是一个常数（通常值为 7）。此外，\(B = \frac{c_0^2}{\rho_0/\gamma}\)，其中 \(c_0\) 是该流体在参考密度下的声速。

2.2. SPH 的数值方法
不同的 SPH 相互作用模型需要用来离散化控制方程（1）和（2）中的每一项。例如，粒子的密度通过以下公式计算：\[\rho_a = \sum_{b \in \mathcal{N}} m_b W\left(|r_a - r_b|, h\right)\]，其中 \(a\) 和 \(b\) 分别是参考粒子和其邻居粒子；\(m_b\) 是粒子质量；\(r_a\) 和 \(r_b\) 是粒子位置；\(W\) 是插值核（例如三次样条或 Wendland 核）；\(h\) 是平滑距离（即粒子相互作用的核范围）。

2.3. 沉浸式浮动隧道的模型
为了计算 SFT 的运动，我们首先需要通过每个边界粒子 \(\mathcal{N}_k\) 整合作用在结构边界壁上的流体力：\[f_k = \sum_{a \in \mathcal{Fluid}} f_{ka}\]，其中 \(f_{ka}\) 是从流体粒子 \(a\) 作用到边界粒子 \(\mathcal{N}_k\) 的单位质量力。然后，通过以下方程计算整个 SFT 的线速度和角速度：\[M\frac{dV}{dt} = \sum_{k \in \mathcal{SFT}} m_k f_k\]，\[I\frac{d\Omega}{dt} = \sum_{k \in \mathcal{SFT}} m_k (r_k - R_0) \times f_k\]，其中 \(M\) 是 SFT 的总质量；\(I\) 是惯性矩；\(V\) 是线速度；\(\Omega\) 是角速度；\(R_0\) 是整个 SFT 的质心。通过积分方程（11）和（12），可以更新围绕 SFT 表面的每个边界粒子 \(\mathcal{N}_k\) 的速度 \(u\)：\[u_k = V + \Omega \times (r_k - R_0)\]。

2.4. 锚泊线的模型
锚泊线限制了 SFT 在六个自由度上的运动。结构通过位于结构下方的锚泊点通过锚泊线固定在海床上。锚泊线的动力学通过集中质量方法进行建模，该方法将整个锚泊链划分为 \(N\) 个等长的段，每个段有 \(N+1\) 个连接节点。作用在每个节点上的所有相关力都是根据 Hall 和 Goupee（2015）的方法计算的。每个段内的拉力来自材料刚度 \(T\) 和内部阻尼 \(C\)，分别通过以下公式计算：\[T_{i+1/2} = \frac{E\pi d^2}{4\epsilon_{i+1/2}\]，\[C_{i+1/2} = \frac{C_{int}\pi d^2}{4\partial\epsilon_{i+1/2}\partial t}\)，其中 \(d\) 是锚泊线的直径；\(E\) 是杨氏模量；\(C_{int}\) 是内部阻尼系数；\(\epsilon\) 是段的应变。然后，每个锚泊节点的动态方程表示为：\[m_{i+1/2}\frac{d^2 r_i}{dt^2} = T_{i+1/2} - T_{i-1/2} + C_{i+1/2} - C_{i-1/2} + W_i + S_i + D_{tra,i}\]，其中 \(T_{i+1/2}\) 和 \(T_{i-1/2}\) 是内部轴向拉力；\(C_{i+1/2}\) 和 \(C_{i-1/2}\) 是阻尼力；\(W_i\) 是浸没重量；\(S_i\) 是海底接触力；\(D_{tra,i}\) 是横向阻力；\(m_i\) 是段的质量；\(m_{a,i}\) 是段的附加质量。

2.5. 模型实现
由于 SPH 是一种无网格数值方法，因此可以通过精确使用粒子自动跟踪流动的自由表面，而无需数值扩散。Crespo 等人（2007）提出的动态边界条件用于处理波浪-结构相互作用的流固界面。数值解采用简单的预测-校正方案进行，计算稳定性遵循 Monaghan（1994）提出的标准 Courant 和粘性约束。SPH 模拟在 3D 大规模工程场景中计算成本较高。Crespo 等人（2015）探索了图形处理单元（GPUs）的巨大潜力，并将其纳入他们的 SPH 代码中，通过在多 GPU 上优化代码架构来执行 SPH 算法的计算，从而开发了 DualSPHysics（可在 https://dual.sphysics.org/ 获取）。此外，MoorDyn 是 Hall（2020）开发的另一个开源库，用于建模锚泊线的动力学。DualSPHysics 和 MoorDyn 之间的耦合通过松散的策略实现。也就是说，首先将 SPH 计算出的 SFT 运动学传递给 MoorDyn，后者计算锚泊链的动力学，并将信息传回 DualSPHPhysics。最终，DualSPHPhysics 整合了锚泊线和流体力，并更新了沉浸式浮动隧道的位置。

3. SPH 可视化
在工程计算中，尽管有许多用于可视化 SPH 计算结果的渲染工具，但总是缺乏真实世界的逼真度。例如，Huang 等人（2011）开发了他们的 VisualSPH 软件，具有交互式的用户界面，用于后处理滑坡流体化运动。Wang 等人（2017）结合了基于视频的重建和基于 SPH 的模拟来可视化动态水面变形和流体-结构相互作用。此外，Ihmsen 等人（2012）提出了各种算法来渲染流体行为，包括喷雾、泡沫和气泡，特别是针对基于粒子的模拟。然而，所有这些之前的可视化技术都依赖于不同研究小组的内部代码，在渲染过程中缺乏统一性。另一方面，Blender 是一个非常强大的渲染软件，已在各种可视化场景中取得了成功应用。将先进的 SPH 模拟与专业的 Blender 结合起来，应该是一种将数值模拟数字化为工程现实的有效方法。在这项工作中，我们将通过使用高质量的计算机可视化来填补这一知识空白。

3.1. VisualSPHPhysics
DualSPHPhysics 的主要结果以 vtk 文件的形式提供。这些文件记录了模拟每个时间步所有粒子的位置、速度、压力等参数。VisualSPHPhysics（Garcia-Feal 等人，2022）使用 Python 和 C++ 开发，便于在商业软件 Blender 中后处理模拟结果。VisualSPHPhysics 作为接口，能够将科学上准确的结果作为原生 Blender 对象真实地表示出来，即它支持将动画导入 Blender。这是通过按给定时间帧顺序加载相应的模拟输出文件来实现的。在加载过程中，这个 VisualSPHPhysics 插件会传输上一个时间步使用的所有材料和纹理属性，以便在后续时间步的新导入文件中保留所有这些视觉属性。此外，VisualSPHPhysics 还可以模拟扩散粒子和水花，以进一步增强真实感。所提出的方法（Ihmsen 等人，2012）在计算精度和 CPU 成本之间保持了适当的平衡，考虑了关键控制因素，如被困空气的循环和波峰检测。各种扩散粒子的分类取决于它们相对于宿主流体的相对位置。例如，流体外的粒子被分类为喷雾，内部的为气泡，表面上的为泡沫。

3.2. Blender
Blender（https://www.blender.org）是一个功能强大的开源 3D 创建工具，用于建模、动画、渲染、视频编辑等。由于其多功能性和免费访问性，它在艺术家、开发人员和爱好者中非常受欢迎。Blender 支持从为游戏和电影创建 3D 资源到制作视觉效果、模拟甚至 2D 动画的各种任务。其强大的社区和丰富的教程库使其成为数字内容创作的初学者和专业人士的绝佳工具。一旦 SPH 模拟数据作为 Blender 对象导入 Blender，所有后续工作都在内部进行。创建最终图像渲染的典型工作流程如下图 2（a）-（d）所示，基于波浪场中自由浮动结构的 2D SPH 模拟。图 2（a）显示了 SPH 数据导入后的模拟结果等值面。目前渲染结果是灰色的，缺乏任何材质效果。图 2（b）显示了在每个渲染元素上应用材质效果后的结果。创建的材质反映了准确的视觉效果，包括水反射和折射、金属光泽等细节；以及木材上的凹凸和漩涡等纹理效果。图 2（c）显示了使用更重但更真实的 Cycles 引擎进行渲染后的场景，与之前的 Eevee 引擎相比。它还展示了完整的场景设置和全局照明。照明根据场景设置定义进行了科学控制，例如散射、亮度，甚至基于材料类型的反射和折射。最后，图 2（d）显示了经过进一步润色的最终产品。这里使用的技术包括泡沫和水花模拟，以及额外的运动模糊和图像合成。

3.3. 集成 SPH 模型的流程图
总之，数字孪生的总体工作流程通过 DualSPHPhysics、VisualSPHPhysics（即 Blender 插件接口）和 Blender 本身来实现。它们之间的相互关系可以在图 3 中示意性地说明，左侧链遵循 Garcia-Feal 等人（2022）的开创性工作。另一方面，右侧链展示了数值数据的标准后处理，例如图表和图形绘制，通过众所周知的软件 Tecplot 和 Origin 进行。侧视图快照显示在图5（a）-（d）中，代表了波浪到达和经过的不同时间点，以及浮动结构的最高和最低位置。下载：下载高分辨率图像（758KB）下载：下载全尺寸图像下载：下载高分辨率图像（248KB）下载：下载全尺寸图像图5. 不同时间和位置的波浪、浮动结构和系泊绳的侧视图快照：(a) 波浪接近；(b) 浮筒的顶部位；(c) 波浪经过；(d) 浮筒的底部位置。数值结果使用之前提到的Blender渲染，以获得逼真的视觉效果。然而，由于MoorDyn库中使用的段数不足，导致系泊绳的渲染效果不够理想，该库与DualSPHysics耦合使用。数字快照展示了在生成谐波后结构的6自由度运动。这表明所提出的模型在研究更复杂的浮动设施（如浮动风力涡轮机和波浪能转换器）方面具有巨大潜力，并有助于设计能够在极端海况下稳定海上结构的有效系泊系统。本次的大规模SPH模拟是在配备24GB内存的Nvidia GeForce RTX 4090 GPU卡上实现的，该GPU卡部署在云服务器上，同时使用的是AMD EPYC 7402 24核处理器（2.79 GHz，10个核心，32GB内存）。模拟50秒的波浪传播（5个波周期）大约花费了30分钟的计算时间，使用了100万个粒子和110K个时间步长。执行过程使用了88MB的CPU内存和190MB的GPU内存。然而，后处理数据的渲染时间略长于SPH主计算时间，尤其是在需要高分辨率图像时。进一步升级显卡应该可以通过商业软件Blender提高渲染效果的质量。

5. 2D水下浮动隧道的应用
在本节中，所提出的SPH模型被应用于一个实际海洋环境中的2D SFT（水下浮动隧道），该隧道受到原型水深和波浪条件的影响。该设计最初是为日本北海道的Funka Bay的交通路线而设计的。概念模型如图6所示，基于Chen等人（2021年）的研究。在实际应用中，它横跨Funka Bay 30公里的范围。SFT内部有一个由两条铁路线和四条车道组成的结构。隧道受到规则波浪载荷的影响，我们使用SPH模拟来研究3自由度的结构运动，以评估SFT的性能。下载：下载高分辨率图像（130KB）下载：下载全尺寸图像图6. 日本Funka Bay水下浮动隧道的概念模型。

5.1. 模型设置和计算参数
如图6所示，浸没深度d是一个关键参数，较大的值表示较小的环境影响和更安全的航行，但建设成本也会更高，反之亦然。Chen等人（2021年）在他们的CFD（计算流体动力学）模拟中研究了三种不同的浸没深度，因为他们还考察了三种不同的水深h。在当前的SPH计算中，我们的主要目的是展示SFT的数字化孪生特性，因此我们将水深固定为h = 100米，浸没深度d = 30米，符合常见的工程实践。隧道的直径根据原型设计值固定为D = 23米。
浮力与自重比（BWR）是另一个反映SFT稳定性的关键指标，也决定了系泊系统的初始张力。为了全面了解BWR的影响，我们在SPH计算中将这个比率从1.0变化到2.0。需要注意的是，Chen等人（2021年）的CFD模拟中设计了总共49个案例研究，其中波周期范围从8.0秒到13.0秒，波高从2.0米到9.3米。在当前的SPH模拟中，我们只选择了波周期固定为T = 13秒、波高H = 9.3米的设置。由于粒子分辨率的限制，没有考虑波高H = 4.0米和2.0米的案例。也就是说，如果使用波高作为控制变量来设置粒子大小，那么较大的水深h = 100米将产生过多的粒子，从而影响计算资源。总共进行了11次不同的SPH测试，BWR的值在1.0到2.0的范围内以0.1的增量进行设置。
数值水槽长度为800米，高度为120米。波浪挡板位于SFT中心外海侧300米处。在水槽末端放置了一个200米长的阻尼区以防止波浪反射。水槽内充满深度为100米的水，隧道顶部边缘与水面的距离为30米。图6中的系泊绳角度β = 60°，遵循Chen等人（2021年）的设计。由于Chen等人（2021年）建议使用刚性模型，因此在SPH计算中假设系泊绳直径为5-7毫米，弹性模量为5.0 [Math Processing Error]× 10^11帕。两条系泊绳的长度各为58.16米。SPH粒子大小被选为0.5米，以平衡计算精度和效率，总共有282,720个粒子参与了计算。

5.2. SFT运动的数字化孪生
图7（a）-（c）显示了一个波浪周期内SFT运动和波浪速度场的数字化图像。渲染过程是在Paraview上实现的，因为Blender在3D可视化方面更为直观有效。如图7（a）所示，在时间t = t0时，入射波峰直接经过水下结构。在波浪力和系泊绳的约束下，SFT开始逆时针旋转并向水槽的岸侧移动。由于上方水动力压力的作用，它还会下沉到最低位置。图7（b）显示，在时间t = t0 + 0.5T时，SFT表现出完全相反的行为，同时岸侧系泊绳受到的张力大于图7（a）中的海侧系泊绳。相比之下，在两个波零交叉点t = t0 + 0.25T或t = t0 + 0.75T时，SFT几乎处于垂直位置，如图7（c）所示，这意味着两条系泊绳承受的张力相似。下载：下载高分辨率图像（609KB）下载：下载全尺寸图像图7. 一个波浪周期内的数字化波浪速度场和SFT运动：(a) t = t0；(b) t = t0 + 0.5T；(c) t = t0 + 0.25T或t0 + 0.75T。
关于速度大小，图7显示，较高速度区域总是出现在水面附近、波前和浮动结构周围的边界区域。然而，应该注意的是，当零交叉波经过结构时，SFT的存在显著改变了流场，如图7（c）所示，较大的速度区域沿着深度扩展到水面。总体而言，SPH模拟很好地捕捉了非线性波浪变形，包括SFT周围的波浪反射和传输，以及结构和系泊绳的运动。

5.3. 模型验证和分析
浮力与重量比BWR与3自由度结构响应（即纵荡、升沉和俯仰运动）之间的定量关系分别显示在图8（a）-（c）中。此外，SPH结果还与Chen等人（2021年）提出的基于有限体积方法的CFD（FLUENT）结果进行了比较，以进行验证。结果显示SPH数据与FLUENT计算结果吻合良好。相对误差分别为4.5%、6.3%和6.5%，分别对应纵荡、升沉和俯仰运动。下载：下载高分辨率图像（154KB）下载：下载全尺寸图像下载：下载高分辨率图像（278KB）下载：下载全尺寸图像图8. BWR与结构运动之间的关系：(a) 纵荡运动；(b) 升沉运动；(c) 俯仰运动。
结果还清楚地揭示了BWR对SFT运动响应的影响。从图8（a）-（c）可以看出，BWR范围内存在一个转折点值，这个值大约在BWR = 1.6左右，适用于所有3自由度运动。这表明1.6应该是一个阈值，反映了当前波浪和结构条件下的SFT稳定性。结构运动中出现这个转折点可以从系泊绳力的组成来解释。当BWR较小时，主要的恢复力是静止水水平面上的系泊绳初始力。它无法有效约束较大的结构运动，随着BWR的增加，SFT的运动也随之增加。然而，当BWR达到1.6时，主要的恢复力变成了动态波浪条件下的额外系泊绳力。因此，这种不断增加的稳定力约束了SFT的运动，导致纵荡、升沉和俯仰运动减小。换句话说，对于具有软系泊系统的SFT，较高的BWR会减少结构的有效质量并增加其对波浪载荷的敏感性。另一方面，对于那些具有刚性系泊绳或高预张力的SFT配置，增加BWR可能会增强以张力为主的行为，最终抑制结构运动。
此外，图8（c）显示，在最不利的条件下，SFT的倾斜角度可达到30°，这可能会对内部交通的安全造成风险。相比之下，如图8（b）所示，升沉运动的幅度远小于其他两种运动模式，而图8（a）中的纵荡运动在计算上更为准确。在另一个阈值条件BWR = 1.0下，所有三个运动分量都表现异常，这意味着系泊绳有松弛的风险，结构可能会不稳定。总之，这些发现表明，无论是过大的还是过小的BWR都不利于SFT的实际设计。前者可能导致严重的结构运动，而后者可能导致系泊绳松弛和结构潜在的倒塌。因此，应综合考虑波浪条件、结构特性和系泊绳强度等因素来确定BWR的最佳值。从本研究来看，BWR = 1.1 - 1.3的值似乎为实际水下浮动隧道的初步设计提供了良好的参考。
图9显示了不同BWR值下的系泊绳力，其中展示了SPH结果和文献数据中的最大岸侧和海侧力。可以看出，随着BWR的增加，这两种力都会增加。这是可以理解的，因为在波浪-隧道相互作用过程中，当浮力效应增强时，会产生更大的系泊力。图9还表明，在当前波高H = 9.3米和水深d = 30米的条件下，岸侧系泊力始终大于海侧系泊力。此外，当BWR非常小（如约1.0 - 1.1）时，系泊力的幅度迅速减小到零，表明系泊绳可能有松弛的风险，隧道结构也可能不稳定。进一步比较SPH计算和Chen等人（2021年）的FLUENT数据，发现岸侧和海侧系泊力的相对误差分别为4.5%和3.9%。下载：下载高分辨率图像（206KB）下载：下载全尺寸图像图9. BWR与最大系泊绳力之间的关系。3D水下漂浮隧道的应用
DualSPHysics在3D工程应用中的强大性能已通过许多成熟的研究得到验证。例如，Suzuki等人（2022年）使用DualSPHPhysics和SWASH的单向耦合模型来可视化洪水事件中室内家具与波浪-物体-结构之间的相互作用。Ng等人（2022年）开发了一个基于3D粒子的流体-结构相互作用模型，用于具有自由表面的水弹性问题，其中采用了体积补偿粒子方法来考虑可变形的固体。此外，在最近的一项研究中，Chatzoglou等人（2025年）提出了一种基于动态Smagorinsky模型的先进湍流闭合方案，与DualSPHPhysics结合使用，以模拟具有河床植被的复杂明渠流动。因此，以下不再提供额外的3D验证案例，因为我们的主要关注点在于模型应用和水下漂浮隧道的数字孪生技术。

大多数关于SFT的研究都是基于传统的外部载荷，如波浪力或波浪-水流力。然而，在实际的海况下，例如地震或山体滑坡引发的海啸波中，SFT的设计还需要考虑极端事件的影响。此外，由于计算资源的限制，大多数大规模的数值SFT模拟都是在2D环境中进行的。在本节中，我们使用所提出的SPH模型和可视化工具Blender来进行更实际场景的3D数字孪生。考虑了两种比较案例，包括单管和双管水下漂浮隧道。首先展示了两种波浪-结构相互作用场景的数字可视化结果，然后进一步研究了在规则波和孤立波作用下的结构运动，旨在检验系统在恶劣海洋环境中的稳定性和安全性，这些内容将在第7.6.1节中提供。

模型设置和计算参数
这里考虑了两种示意性设计的SFT场景，包括单管和双管配置，以研究规则波作用下的流体动力学和结构运动。SFT的流体动力学特性与其横截面配置密切相关。单管SFT概念已应用于日本Funka湾的实际设计中，而双管SFT则应用于挪威的E39高速公路项目中。

所提出的数值波浪水槽长度为70米，宽度为16米，高度为15米。水槽最初充满水，水深为8米。单管和双管SFT分别通过四根系泊绳固定在水槽底部，系泊绳的倾斜角度为60度。每个圆柱体的长度为14米，直径为1.0米，浮力与重量比为1.25，浸入水中的深度为2米。两个SFT位于水槽横向的中心位置。在纵向方向上，波浪挡板与结构系统的重心之间的距离为30米。双管SFT由两个大小相同的圆柱体组成，它们之间相隔1米，通过三个均匀分布的接头连接。

在SPH模拟中，使用设计波高H=2米、波周期T=5秒的规则波。总模拟时间为100秒，覆盖了20个波周期的传播过程。SPH粒子的分辨率为0.2米，数值模拟中使用了大约110万个粒子。每根系泊绳的长度为8.15米，线刚度为5×10^6牛顿，直径为0.1米。每根系泊绳被划分为十个段，以便与Moordyn库进行耦合。在水槽末端放置了一个长度为20米的数值波浪阻尼层，以防止波浪反射。

6.2 SFT运动的数字孪生
数值SPH数据通过插件VisualSPHysics导入Blender，图10（a）-（c）展示了双管SFT配置在一个波周期内的运动情况，展示了波浪-结构-系泊绳相互作用过程中的典型运动特征。图10（a）对应于波浪到达的时刻。此时，过大的流体动力压力将SFT推向岸边，并使其向顺时针方向严重倾斜。图10（b）显示当波峰刚刚经过SFT上方时，SFT已经开始逆时针旋转。然后图10（c）显示当波浪经过结构时，观察到完全相反的现象，即SFT主要向海侧倾斜。

为了对比，图11（a）-（c）展示了单管SFT的数字孪生结果，发现波浪和结构运动非常相似。然而，单管SFT的旋转运动似乎更快，因为在波峰刚刚经过结构上方时，它已经达到了最大的逆时针位置。当波浪过去后，SFT大多恢复到直立平衡状态。

6.3 单管和双管SFT运动的比较与分析
图12（a）-（c）分别展示了单管和双管SFT的涌浪、纵荡和横摇运动的计算时间历史。需要注意的是，时间原点实际上是从数值数据中的t=50秒开始的。在SPH计算中发现，在第10到第20个波周期的时间窗口内，波浪场的重复性得到了很好的保持，而在前10个波之前观察到了不真实的SFT运动振荡。图12（a）和（b）显示两种SFT在两个不同方向上的过渡运动幅度大致相同。然而，图12（c）显示双管SFT的旋转运动（即横摇）显著减少。与单管SFT配置相比，顺时针和逆时针运动分别减少了约45%和60%。这表明双管SFT在运行中可以实现更好的稳定性和安全性。还可以注意到，两种SFT设置下的纵荡运动都表现出明显的更高频率成分。

图12中的两种SFT运动可以这样解释：涌浪和纵荡运动主要取决于结构的投影面积和位移体积。在单管和双管SFT的情况下，圆柱体的几何形状和材料属性相同，因此作用的流体动力对其过渡运动产生类似的影响。另一方面，旋转运动（即横摇）强烈依赖于作用力相对于旋转中心的分布。双管SFT的形状显著改变了力臂和旋转惯性，即使流体动力相同，产生的旋转力矩也不同，从而导致较小的角运动。此外，在双管SFT配置下，两个圆柱体间隙之间的相互作用流可以改变波浪速度和压力分布，起到一定的缓冲作用，减轻了原始波浪场的旋转冲击。

7. 单管和孤立波作用下的SFT运动
在本节中，我们使用SPH模拟来揭示这种差异。考虑到SPH实现的技术可行性，仅生成了一个孤立波。这个波的高度与第6节中的规则波相同。所有其他计算参数，包括流体动力学和结构属性，都与规则波时使用的参数保持一致。由于单管和双管SFT之间的比较没有本质区别，因此下面以双管SFT为例进行说明。

在SPH设置中，数值水槽的长度设置为规则波水槽的三倍。这是因为孤立波是通过在岸边推动波浪挡板产生的，这可能会有效减小水槽的体积，并在末端形成人工水面堆积。随后的回流涌浪会降低模拟精度。较长的数值水槽通过提供额外的空间来延迟回流波的生成，从而有效解决了这个问题。因此，在当前的SPH模拟中使用了大约400万个粒子。

7.1 两种波作用下的结构运动比较
图13（a）-（c）分别展示了双管SFT在规则波和孤立波作用下的涌浪、纵荡和横摇运动的计算时间序列。需要注意的是，时间原点实际上是从数值数据中的t=50秒开始的。在SPH计算中发现，在第10到第20个波周期的时间窗口内，波浪场的重复性得到了很好的保持，而在前10个波之前观察到了不真实的SFT运动振荡。图13（a）和（b）显示两种SFT在两个不同方向上的过渡运动幅度大致相同。然而，图13（c）显示双管SFT的旋转运动（即横摇）显著减少。计算结果显示，与单管SFT配置相比，顺时针和逆时针运动分别减少了约45%和60%。这表明双管SFT在运行中可以实现更好的稳定性和安全性。还可以注意到，两种SFT设置下的纵荡运动都表现出更明显的高频成分。

图13中的两种SFT运动可以这样解释：涌浪和纵荡运动主要取决于结构的投影面积和位移体积。在单管和双管SFT的情况下，圆柱体的几何形状和材料属性相同，因此作用的流体动力对其过渡运动产生类似的影响。另一方面，旋转运动（即横摇）强烈依赖于作用力相对于旋转中心的分布。双管SFT的形状显著改变了力臂和旋转惯性，即使流体动力相同，产生的旋转力矩也不同，从而导致较小的角运动。此外，在双管SFT配置下，两个圆柱体间隙之间的相互作用流可以改变波浪速度和压力分布，起到一定的缓冲作用，减轻了原始波浪场的旋转冲击。

7.2 规则波和孤立波作用下的SFT运动可视化
为了揭示上述两种波条件下SFT运动差异的机制，图14（a）和（b）分别绘制了规则波和孤立波作用时波浪接近结构时的流速场。显然，前者的流速能量主要集中在水面附近，分布较广，而后者则更集中在流动深度方向。因此，孤立波能够提供更强的能量和动量，使水下漂浮隧道产生更大的运动。

8. 结论
本文结合了数值模拟和计算机可视化技术，实现了复杂海洋环境中SFT运动的数字孪生。该模型基于开源代码DualSPHPhysics，可视化通过商业软件Blender实现。该模型首先被应用于规则波作用下的二维SFT（Semi-Particle Method）运动中，并通过现有的CFD（Computational Fluid Dynamics）数据对SPH（Smoothed Particle Hydrodynamics）计算结果进行了验证。随后，该模型被用来研究由单圆管和双圆管组成的更复杂的3D SFT结构，这些圆管都固定在海床上。通过比较这两种结构的运动特性，可以揭示波浪与结构之间的复杂相互作用。为了提供更实际的见解，还进行了另一项对比研究，以探讨规则波与孤立波对SFT的影响差异，旨在考察极端海洋条件对固定浮动结构的影响。研究发现，由于孤立波具有更高的能量和动量集中度，因此可能更具危险性。总体而言，所提出的SPH模型能够通过全面预测实际设计和运行中的流体动力响应和结构响应，为海洋工程领域提供一种有前景的工程解决方案。此外，通过与渲染软件Blender的结合，该研究实现了数值结果的高度真实可视化，从而能够在许多工程场景中对SPH模拟进行数字孪生处理。

**作者贡献声明：**
- 邵家琪（Jiaqi Shao）：撰写、审阅与编辑、初稿撰写、可视化处理、软件开发、方法论设计。
- 胡浩（Hao Hu）：撰写、审阅与编辑、资源收集、数据整理、概念构思。
- 梁东方（Dongfang Liang）：撰写、审阅与编辑、初稿撰写、项目监督、资金筹集。