1. Introduction
A Submerged Floating Tunnel (SFT) is generally suspended 20–50 m below the water surface and can cross the two sides of the strait. It is recognized as the most promising form of cross-sea transportation in the 21st century [
1]. According to different supporting systems, submerged floating tunnels are divided into pontoon type, pier type, anchor cable type, and free type [
2,
3]. The SFT is mainly affected by waves and currents, so the study of its hydrodynamic response is the key to the safe design of a submerged floating tunnel structure.
Theoretical analysis methods can quickly predict the response features of an SFT under hydrodynamic loading. Fan et al. [
4] studied the resonance of a single-span SFT undergoing vortex-shedding. They simplified the tunnel structure as a simply supported Euler–Bernoulli beam at both ends, calculated the fluid damping and drag force using the Morison formula, and found that in a specific range of incoming flow velocities, the fluid–structure coupling effect and nonlinear characteristics have a significant impact on the structural response. Ge et al. [
5] solved the wave field and the interaction between the wave and the SFT based on the potential flow theory and the boundary element method, respectively. They used the prototype design parameters of the SFT in Qiandao Lake as a case to carry out systematic analysis. Numerous scholars have conducted research on the dynamic response of the SFT under wave or current loads based on theoretical analysis methods [
6,
7,
8]. Although the theoretical analysis has high computational efficiency, the many simplification conditions limit its application. For example, the potential flow theory neglects the fluid viscosity and cannot accurately calculate the viscous shearing force. The Morison formula can only calculate the force on small-scale components [
9]. In addition, it is difficult to address the nonlinear problem of a free liquid surface through theoretical analysis.
The model test is an important method for investigating the SFT’s dynamic response. Many researchers have studied the SFTs’ response under the action of waves, flows, or wave-flow coupling [
10,
11,
12,
13]. They mainly focused on the motion of the six degrees of freedom of the SFT, tunnel pipe acceleration, anchor cable tension, and tunnel pipe section pressure. The dynamic response of the SFT is mainly controlled by two sets of parameters. One set corresponds to the loads, including wave height, wave period, and incoming flow velocity. The other set pertains to coefficients characterizing the SFT system, including submersion depth, weight-to-buoyancy ratio, and cable inclination angle. These two sets of parameters have been the central focus of previous experimental investigations. Model tests also have deficiencies; for example, it is difficult to simulate large and steep waves, the test cycle is long, and the cost is high.
With the development of computer performance, numerical simulation has gradually become another important means to study the dynamic response of SFTs. Yang et al. [
14] established a finite element model (FEM) of a SFT and anchor cables based on AQWA. They studied the dynamic response of anchor cables under the action of waves and flow and found that regular waves and uniform flow caused periodic motion of the structure and large displacement and low-frequency movement. Zou et al. [
15] numerically studied the dynamic response characteristics of the SFT-anchor cable coupling system and found that the tube motion response increased with increasing wave height and mooring angle. Luo et al. [
16] numerically solved the simplified differential equations for the tunnel pipe and anchor cable system. The response amplitude of the SFT and the anchor cable tension increase with increasing wave height and wave period and decrease with increasing submersion depth. Zou et al. [
17] used the one-way coupling method of CFD and FEM to study the response of SFT in extreme sea conditions. They pointed out that waves with large wave heights and large periods could avoid exciting the first-order dominant mode of SFT. Chen et al. [
18] numerically studied the nonlinear interaction between waves and SFTs. They proposed that an appropriate buoyancy-to-weight ratio could prevent the large movement of SFTs under adverse wave conditions. Peng et al. [
19] established a CFD model of the interaction between waves and SFTs. They found that the CFD model based on solving the Navier–Stokes equations could fully consider processes such as the nonlinearity of the flow field, vortex formation, boundary layer separation, and wave breaking.
In summary, a large number of researchers have numerically studied the dynamic response of a SFT under the action of waves or flows alone, but there are relatively few studies under the action of wave-flow coupling. In addition, many numerical simulations are based on potential flow theory, which cannot consider the force generated by vortex shedding. In this paper, a CFD model of an SFT is established, which can not only calculate the drag force and lift force caused by vortex shedding by considering the fluid viscosity but also capture nonlinear load characteristics such as free liquid surface reversal and breakage. SFTs in the actual marine environment are generally subjected to the combined action of wave and flow loads, and the directions of wave propagation and flow are generally not perfectly aligned. Hence, when investigating the coupling effects of waves and flows, the impacts of factors such as wave height and flow velocity, as well as the directions of wave and flow propagation, bear enhanced practical significance for the dynamic response of SFTs.
2. Numerical Model
The commercial CFD software STARCCM+ V2206 is utilized for numerical simulations. The fluid domain is discretized based on the finite volume method. In the discretized control volume, the governing equations of fluid flow are converted to algebraic equations. The algebraic equations are iteratively solved, combined with the boundary conditions, and initial conditions to obtain the numerical solution of the real flow [
20]. The compressibility of fluid is not considered, and the controlling equations of fluid flow are the time-averaged continuity equation and Reynolds equation (RANS) shown in Equations (1) and (2):
where
is the time-average value of velocity,
is the time-average value of pressure,
μ is the dynamic viscosity of the fluid,
ρ is the density of the fluid,
is the time-average of body force,
gi is the gravitational acceleration vector, and
is the Reynolds stress term and a second-order tensor. The Reynolds stress term appears in Equation (2), and the governing equations of fluid motion are not closed. The
k-
ε turbulence model is chosen as the supplementary equation to establish the relationship between the Reynolds stress and the average velocity, making Equations (1) and (2) closed and solvable. The Eulerian multiphase flow model was used to create the air phase and the water phase, and the volume-of-fluid (VOF) method captured the interface of the two phases. The left and right boundaries of the calculation domain are set as “velocity inlet” boundaries. This boundary condition can simulate the real physical wave-making process by pushing the plate. The wave parameters are specified based on the intermediate water depth wave theory, and the wave propagation direction is from the left side to the right side. With the “wall” boundary at the bottom and the SFT, this is equivalent to a no-slip wall boundary condition. The “pressure outlet” boundary is at the top. The front and back sides were set as “symmetry planes”, and this boundary condition made the velocity gradient normal to the plane zero, thus converting the simulation problem into a 2D problem. At the same time, the “wave force” option is used on the “velocity inlet” boundary, and the “wave force length” is set to twice the wavelength. This method of wave mitigation allows the selection of a smaller computational domain size and improves computational efficiency [
21]. The selection of boundary conditions and computation domain size are shown in
Figure 1a. The fluid calculation domain is subdivided based on the “trimmed mesh”. Due to the large magnitude of rigid body motion in the SFT, the overset mesh method is chosen to discretize the calculation domain around the SFT. To ensure the accuracy of grid interpolation, the mesh sizes of the background region and the overlapping region should be approximately matched. Near the water surface, significant interactions occur between waves and structures. Therefore, the grid in the vicinity of the water surface needs to be refined. The refinement criterion is to divide a unit wavelength into 80 grids and a unit wave height into 20 grids. Due to the viscosity of the fluid, there is a boundary layer with a certain thickness on the surface of the SFT. In the boundary layer, the velocity gradient is relatively large, and the shear stress cannot be ignored. Therefore, the “prismatic layer mesh” was used for the SFT surface to accurately calculate the viscous shear force.
Figure 1b shows the meshing on the cross-section perpendicular to the
y-axis. The “DFBI” module in STARCCM+ can solve the force and movement of the SFT, and the catenary option under this module can be used to establish the numerical calculation model of the SFT and anchor chain system [
22].
3. Model Validation
Three sets of grids with different sizes were established, with dimensions of 0.226 m, 0.16 m, and 0.113 m, corresponding to grid quantities of 280,000, 790,000, and 2,370,000, respectively. The comparison shows that the calculation on the fine mesh shows that the dimensionless sway amplitude (
δx/
Hi), the dimensionless heave amplitude (
δz/Hi), and the dimensionless pitch amplitude (
αB/(2
Hi)) were all larger and closer to the experimental values. The deviations between the calculated values of
δx/
Hi,
δz/
Hi, and
αB/(2
Hi) obtained from fine and moderate grids are 2.6%, 9.5%, and 6.5%, respectively, as shown in
Table 1. Further refining the grid based on the moderate grid has a significant impact only on
δz/
Hi. Moreover, the deviation between the calculated
δz/
Hi obtained from the moderate grid size and the experimental value is already relatively small. Considering both computational accuracy and efficiency, the subsequent calculations will adopt the moderate grid size as the standard for grid partitioning.
The model experiment involves an SFT with a rectangular cross-section. The relevant parameters of the SFT are provided in
Table 2. Here, B, W, and h denote the SFT model’s length, width, and height, respectively. A series of laboratory experiments were conducted within a 2D glass-walled wave flume situated at the Coastal and Ocean Engineering Laboratory, Department of Civil Engineering, Nagoya University, Japan. The wave tank’s dimensions are 30 m in length, 0.7 m in width, and 0.9 m in depth. The wave channel’s bottom is flat, and a piston-type wave paddle is installed at the closed end of the wave flume to generate waves. On the opposite end of the wave flume, a wave absorber constructed of a rubble mound is employed to minimize the impact of reflected waves. Positioned midway between the wave maker and the wave absorber is a model breakwater that spans the entire width of the wave tank. Consequently, the flow can be treated as two-dimensional. Further details about the experiment can be found in Peng et al. [
19]. The model has a mass of
m = 28.6 kg, and the moment of inertia is
I = 0.435 kg m
2. The immersion depth of the model is
ds = 0.102 m. The parameters of the model experiment are outlined in
Table 3. In this table,
Ti,
Hi, and
Li represent the wave period, wave height, and wavelength, respectively. The variable
d stands for the water tank depth, while
β indicates the inclination angle of the anchor cable and
Fi corresponds to the initial tension of the anchor cable. In the model experiment of Peng et al. [
19], a total of four anchor cables were symmetrically arranged. The numerical computational model incorporates two anchor chains along the mid-longitudinal section. Subsequently, the forces calculated for each anchor cable were equivalently distributed onto the two anchor chains, achieving consistency with the model experiment. The side view of the model test is illustrated in
Figure 2.
Figure 3 presents a comparison between numerical results and experimental data. In the figures, “EXP” represents the outcomes from the model experiment, while “CFD” signifies the results obtained through numerical computation. Here, “
ax” and “
az” stand for the lateral and vertical accelerations of the model, while “
Foff” and “
Fon” correspond to the offshore and onshore tensions of the anchor cable. In accordance with
Figure 3, the time series of various response values of the SFT obtained through numerical simulations exhibit good agreement with the model experiment. Specifically, the amplitude values of SFT’s pitch, roll, and lateral acceleration obtained through numerical solutions are slightly lower than those observed in the model test, whereas sway, vertical acceleration, and the offshore and onshore tensions of the anchor cable are observed to be a little larger than the model test values. Perhaps modifying the added mass coefficient during the numerical simulation process could yield results that are more closely aligned with the model experiments. However, due to the lack of clear experimental guidance on how large the added mass coefficient should be, we did not adjust the added mass coefficient. Furthermore, the nonlinear response characteristics can be effectively captured through numerical simulations. Because the CFD method not only considers the fluid’s viscosity but also accounts for nonlinear loads caused by wave breaking or slamming processes, as shown in
Figure 4. These nonlinear loads, such as second-order wave forces and slamming forces, have a significant impact on the response of the SFT. In the experiment of Peng et al. [
19], the SFT is connected via a bottom fixed end and steel chains. The connection between the steel chain end and the SFT is close to a hinge joint, allowing rotation between the two components. Conversely, in the numerical computations, the anchor chain is fixedly connected to the SFT. This difference in connection might contribute to the observed deviations in the corresponding response values.
To further validate the accuracy of the numerical model presented in this paper, the dynamic responses of a circular cross-section SFT model under regular wave conditions were calculated and compared with relevant model tests and numerical simulation results, as shown in
Figure 5. The numerical simulation results by Peng [
19] were also presented in the figure. The relevant parameters for the model test are listed in
Table 3. In comparison to the model test conditions for a rectangular cross-section tunnel, the initial tensile forces of the mooring line slightly increased to
Fi = 35.04 N. The relevant model parameters for the circular cross-section SFT are presented in
Table 4, where
R represents the radius of the circle. As evident from
Figure 5, the non-dimensional roll,
Foff, and
Fon responses of the circular cross-section SFT obtained from our numerical simulations closely match the model test results and Peng’s numerical simulation results. The respective response amplitudes exhibit errors of 14%, 8%, and 3.8% when compared to the model test. The calculated dimensionless sway amplitude is even closer to the model test results than Peng’s calculations, with an error of 2.5%. However, both our numerical calculations and Peng’s results yield significantly lower heave motion amplitudes compared to the model test, and these oscillations exhibit a frequency close to the wave period without displaying the high-frequency oscillations observed in the model test.
5. Conclusions
A numerical computational model was established to simulate the system comprising SFT and anchor chains. The accuracy of the numerical model was validated through comparison with model test results. The influence of parameters such as wave height, flow velocity, and wave and flow propagation direction on the dynamic response of the SFT was analyzed. The main conclusions are as follows:
Under the isolated action of a wave, an increase in wave height leads to a significant increase in the maximum values of ax, az, Foff, Fon, and δz/Hi, accompanied by an augmentation in the nonlinear characteristics. While wave height variation does not significantly affect δx/Hi, the maximum value of αB/(2Hi) initially increases and then decreases with increasing wave height.
When waves and flows propagate in the same direction with a fixed flow velocity, the effects of wave height on the SFT are similar to those observed under isolated wave conditions. An increase in wave height results in larger maximum values for ax, az, Foff, Fon, and δz/Hi, while δx/Hi and αB/(2Hi) decrease. Moreover, under combined wave and flow propagation, the δz/Hi for the SFT increases for all wave heights compared to cases with regular waves acting alone, while δx/Hi and αB/(2Hi) exhibit smaller vibration amplitudes. The presence of flow causes the equilibrium position of the lateral motion of the SFT to shift along the direction of wave and flow propagation, with this displacement influenced by wave height. Notably, higher wave heights lead to smaller lateral vibration amplitudes. Additionally, the SFT rotates a certain angle against the incoming flow and then undergoes periodic rotational motion at the new position. Similarly, as wave height increases, the angle of counterclockwise rotation diminishes. As wave height increases, the non-linear characteristics of the pitch motion of the SFT amplify.
When waves and flows propagate in the same direction with a fixed wave height under the condition of relatively high flow velocity, the SFT displaces significantly in both the direction of flow and the depth of the water, causing it to deviate considerably from its original equilibrium position. At this juncture, the SFT is largely unaffected by wave loads and predominantly experiences periodic lift and drag forces resulting from vortex shedding along the tunnel’s surface. On the side opposite the flow, the anchor chain remains completely slack.
When waves and flows propagate in opposite directions, the maximum values of ax, az, Fon, and Foff for the SFT exhibit significant increases. Structural considerations are required in the design process to prevent anchor chain rupture.