SPH Simulation of Solitary Wave Interaction with Cylinder

1. Introduction

Tsunami is a kind of disastrous wave that can cause severe damage to marine structures. Various marine structures have cylindrical components, such as the columns of wind turbines and oil/gas platforms, piers of sea-crossing bridge and subsea pipelines. Studying the mechanism of tsunami impacting on cylindrical structures is of great significance for the design and disaster prevention of marine structures. In this context, the solitary wave (which represents tsunami due to their resemblance) interaction with horizontal and vertical cylinders is numerically studied in this work.

Computational Fluid Dynamics (CFD) has become a numerical tool with great potential based on the continuous development of computer functions. Among them, the Smoothed Particle Hydrodynamics (SPH), as a Lagrangian meshless particle method, can naturally capture rapidly moving interface and handle large deformation of free surface. Because of its own advantages, it has been widely developed and utilized in the field of ocean engineering (Luo et al. [2021]; Zhang et al. [2022, 2023]; He et al. [2023]; Ju et al. [2023]). In the last two decades, the SPH theories and open-source codes have been getting significant developments, e.g. GPUSPH (Hérault et al. [2010]; Shi et al. [2017]), DualSPHysics (Crespo et al. [2015]), SPHinXsys (Zhang et al. [2021]), SPHERA (Salis et al. [2022]) and SPHydro (Lyu et al. [2023]).

The SPH method has been used for simulating solitary wave interaction with cylindrical structures. For example, Wei et al. [2015] used GPUSPH to study the impact of tsunami on bridge piers, in which the velocity field around cylindrical piers and wakes were studied and the numerical accuracy was verified. Based on SPH, Pan et al. [2016] studied the motion effect of a solitary wave impacting upon a floating tension-leg platform (TLP) whose tension-leg consists of four vertical cylinders. Aristodemo et al. [2017] conducted the SPH numerical and experimental study of the horizontal and vertical hydrodynamic forces induced by solitary waves on a horizontal cylinder. Fathi & Ketabdari [2018] studied the climbing and overtop phenomena of solitary waves on semicircular breakwater using the SPH method, and the validation is addressed by comparing with experimental results. Domínguez et al. [2019] studied the impact of double solitary waves on two cylindrical reservoirs based on the DualSPHysics model, and the research results show that DualSPHysics can reflect the main hydrodynamic forces when solitary waves interact with coastal structures. Han and Dong [2020] studied the interaction of solitary waves with underwater breakwaters consisting of rectangular and semicircular sections based on the SPH method, and verified the validity of the proposed numerical model. By using the SPH open-source code DualSPHysics, Pringgana et al. [2021] studied the impact of tsunami on coastal structures in different layouts. Reis et al. [2022] studied the horizontal and vertical pressure distributions imposed by tsunami-like waves on elevated structures, and compared the experimental data to calibrate the numerical model. Francis et al. [2023] simulated the interaction between solitary waves with a thin and upright porous barrier using the SPH method, and studied the energy dissipation characteristics of the porous barrier at different water depth ratios.

In recent years, the weakly compressible SPH code SPHinXsys (Zhang et al. [2021]) has obtained significant developments and has been demonstrated to be a powerful model in simulating fluid–structure interaction problems. This study extends the application of SPHinXsys to the solitary wave and its interaction with a submerged horizontal cylinder and a surface-piercing vertical cylinder. The key hydrodynamic parameters including wave elevation, wave profile and impact load on cylinder will be elaborated in comparison with the results produced by the experiment. The physics behind the wave–structure interaction process is then discussed, as well as the accuracy of SPHinXsys is verified. In what follows, Sec. 2 presents the key mathematical formulations of SPHinXsys. Section 3 introduces the numerical examples. Section 4 concludes.

2. Mathematical Formulations

2.1. Governing equations

The governing equations are the mass and momentum conservation (i.e. Navier–Stokes equation) laws in the Lagrangian form

where $\rho$ is the density of a fluid, v is the particle velocity vector, p is the fluid pressure, $\nu$ the kinematic viscosity, g is the gravity acceleration and t is the time. The computational domain is discretized by nonconnecting particles which carry a fixed mass and move under the joint action of the gravity, pressure gradient and viscous force as depicted in Eq. (1).

This work adopts the weakly compressible SPH method. The fluid pressure is obtained through an equation of state (Zhang et al. [2021]) :

where ${\rho }_0$ is the reference density, $c_0$ is the artificial speed of sound adopted to be 10 times of the projected maximum particle velocity during the simulation.

2.2. Numerical formulations of SPHinXsys

The governing equations are discretized in the particle interpolation form (Zhang et al. [2021])

where ${\bar{v}}_{ij}=(v_i+v_j)/2$ and ${\bar{p}}_{ij}=(p_i+p_j)/2$ are the particle average velocity and pressure, $m_j$ the mass of particle j which is constant in simulation and W the kernel function. A fifth-order Wendland kernel function (Wendland [1995]) is used in particle interpolation in this study : where ${\alpha }_d$ is a parameter equal to $7/(4\pi h^2)$ for two-dimensional (2D) simulations and $21/16\pi h^3$ for 3D simulations, and $q=|r_i-r_j|/h=r_{ij}/h$ with h being the smoothing length.

A key feature of the SPHinXsys for fluid dynamics is adoption of the Riemann solver (Zhang et al. [2017]) to mitigate numerical noises (e.g. the spurious pressure fluctuation). In the Riemann solver, a fictitious interface is founded between two particles, and equations of state are constructed at the left and right sides of the interface :

where L and R indicate left and right states, respectively.

A low-dissipation Riemann solver :

where $\beta =min(3max(U_L-U_R,0)/c,1)$ is the low-dissipation confinement and $c=({\rho }_L{c}_L+{\rho }_R{c}_R)/({\rho }_L+{\rho }_R)$ is applied. The variables of $p^{\ast }$ and $v^{\ast }$ are then used to substitute the particle average terms of ${\bar{p}}_{ij}$ and ${\bar{v}}_{ij}$ in Eq. (3). This results in the particle-interpolated governing equations.

For the time marching, the time step size is concluded by the Courant–Friedrichs–Lewy (CFL) condition :

where the coefficient CFL is adopted to be 0.6 and h is the smoothing length.

In SPHinXsys, the wall boundaries of the cylinders and the computational domain are modeled by dummy particles, while the fluid–solid particle interaction is resolved by one-side Riemann solver (Zhang et al. [2017]). Through these, the nonslip and impermeable conditions are imposed for velocity, as well as the Neumann boundary condition is applied for pressure. Specifically, the pressure values on wall particles are solved by $P_R=P_L+{\rho }_{f}g\cdot r_{fw}$, where $r_{fw}=r_w-r_f$ with $r_w$ and $r_f$ being the positions of a wall particle and its relating fluid particle (Zhang et al. [2021]).

3. Numerical Examples

3.1. Solitary wave generation and propagation

The free surface formula of Boussinesq-type solitary wave reads (Goring [1978])

where H is the wave height and $c_{\mathrm{\text{w}}}$ is the wave speed. The wave celerity with initial water depth being d can be assessed as and the parameter k can be solved as

The particle method can model the moving wave maker with ease due to its Lagrangian nature. Therefore, this study adopts the moving wave maker approach. For this, a simplified time series function of the wave maker movement is adopted as (Farhadi et al. [2016]; Wang et al. [2020])

where S is the stroke of the wave maker and $\tau$ is the duration of the wave maker movement, calculated, respectively, as

At the downstream side of the domain, a sponge layer is set to dissipate wave energy and minimize the wave reflection. The fluid particles in the sponge layer are influenced by excessive viscous force and their velocities are governed by $v=v_0(1.0-\mathrm{\Delta }t\alpha (\frac{r-r_0}{r_1-r_0}))$, where $v_0$ is the initial velocity of a fluid particle just entering the damping zone, $\mathrm{\Delta }t$ the time step, $r_0$ and $r_1$ the starting and ending locations of the damping zone, respectively, and $v$ the velocity at the location $r$. The coefficient $\alpha$ controls the damping intensity and is adopted to be 5.0 in this study.

The 2D solitary wave generation and propagation by the SPHinXsys model and the convergence with respect to the particle size (dp) and smoothing length (h) have been demonstrated (Cai et al. [2022]). In this paper, the 2D case is based on the verified model, so the particle size $\mathrm{\text{dp}}=0.005$m and smoothing length $h=1.7$dp are used. A new research on convergence of particle size with $h=1.7$dp is studied in a 3D empty water tank. The verified model is shown in Fig. 1, the parameter settings of model are described in Sec. 3.3 (only one vertical cylinder is added).

The simulated time histories of wave height are presented in Fig. 2, with the relative errors listed in Table 1. Four particle sizes are tested, i.e. 0.04, 0.03, 0.02 and 0.018m. The minimum particle size is determined to make the simulation time on a workstation equipped with an AMD Ryzen Threadripper PRO 5975WX CPU (32-Cores 3.60GHz) is manageable. As can be seen, the wave crest approaches the analytical solution with a decrease in particle size, with the relative error reducing from 13.5% to around 5%. Note that the error in the case of particle size $\mathrm{\text{dp}}=0.018$m is slightly larger than that of $\mathrm{\text{dp}}=0.02$m, which is due to the numerical noises of the simulation (albeit significantly mitigated) (Ren et al. [2023]). This shows the particle convergence of the present SPH model in simulating the solitary wave case. Balancing computational efficiency and accuracy, the particle size of $\mathrm{\text{dp}}=0.02$m is adopted in the following simulations.

3.2. Solitary wave interaction with a submerged horizontal cylinder

The case of solitary wave interaction with a submerged horizontal cylinder is studied with experimental validation (Aristodemo et al. [2017]). As sketched in Fig. 3, the effective length of the computational domain is adopted to be $L_{\mathrm{\text{w}}}=6.5$m and a piston-type wave maker is placed on the left side of the domain. A fixed cylinder of diameter 0.127m is placed such that its center is 2.55m from the initial location of the wave maker and 0.2m from the bottom of the tank. Note that the wave flume length is reduced in SPHinXsys simulation as compared to the experimental wave flume for saving computational cost. Therefore, the numerical arrival time of the solitary wave crest has been shifted by the reduced flume length divided by the wave celerity, and the time instant when the wave crest reaches right above the cylinder is set $t=0$s. Other parameters considered are: water depth $d=0.4$m and solitary wave height $H=0.071$m. Studied include the wave elevation at $x=2.55$m (WG in Fig. 3), as well as the horizontal and vertical components of the wave force applied on the cylinder. In simulations, a fixed particle size of $\mathrm{\text{dp}}=0.005$m is adopted. With this setup, the numerical domain contains a total of 102591 particles and the computational time is about 11min for simulating 1s physical process based on a personal computer equipped with an Intel i7-10700 CPU (8-Cores 2.90GHz).

The wave profiles and velocity fields at typical time instants of the wave–cylinder interaction process are presented in Figs. 4 and 5. With the arrival of the wave crest, the velocity on the upper and lower sides of the cylinder gradually increases, forming the flow around the cylinder (see Figs. 4(a) and 4(b)). As the wave crest moves down, the flow behind the cylinder forms two vortices that rotate in opposite directions, but the overall flow velocity reduces gradually (see Figs. 4(c) and 4(d)). In this process, the solitary wave does not break, but the internal flow field is disturbed to generate a large flow field velocity, which will directly impact and damage the pipeline. At the same time, the stream will generate horizontal thrust and alternating vertical force (see Fig. 8), which can cause the dynamic response of the structure. All of the above are the hidden hazard for submerged structures. The present SPH model captures the flow field in the process satisfactorily with regular particle and smooth pressure distributions.

The wave elevations predicted by the present SPH model are compared with published results (Aristodemo et al. [2017]) in Fig. 6. It is noted that results obtained by SPHinXsys are close to the experimental results, and the relative errors of wave crest with experimental results are −0.289%. Generally, the present SPH results are in a closer agreement with the experimental data.

Wave forces exerted on the cylinder are analyzed through integrating the pressures on the boundary particles of the cylinder. The data produced in this way contain noises, as shown in Fig. 7. These noises are related to the spurious fluctuations contained in the pressure results. For brevity, the force results are filtered by a low-pass filter with the cut-off frequency of 20Hz. The filtered force components are then normalized by dividing the maximum horizontal force ($F_{\mathrm{\text{Hmax}}}=20$N) for comparison with the experimental data, as shown in Fig. 8.

For the horizontal force $F_{\mathrm{\text{H}}}$, the overall trend and magnitude predicted by the SPH model are in good agreement with the published experimental data. For the vertical force $F_{\mathrm{\text{v}}}$, the positive peak predicted by the SPH model is slightly larger than that of the experimental result, and the negative forces show some discrepancies. Regarding the negative forces, that captured by SPH is short-lived, and it is only around $t=0$s when the wave crest is directly above the cylinder, while that in experiment lasts longer. Due to the obstruction of the cylinder, the stream bypasses the cylinder. The lower side of the cylinder has a small watershed, which produces a large flow field velocity, resulting in a decrease in pressure. But the upper side of the cylinder has a large watershed where being a small flow field velocity and a strong pressure. The above may result in a negative pressure. The wave force gap captured by the SPH model is the focus of future research.

3.3. Solitary wave interaction with a surface-piercing vertical cylinder

The benchmark case of solitary interaction with a surface-piercing vertical cylinder (Yang et al. [2021]) is studied. As sketched in Fig. 9, the origin of the coordinate system is at the lower left corner of the tank at the initial instant, the effective length of the computational domain is $D_L=14$m and a piston-type wave maker is placed on the left side of the domain. A fixed cylinder of radius $R=0.125$m is placed such that its center is 8m from the initial location of the wave maker and its height is 1.2m from the bottom of the tank. Note that the length of the flume has been reduced in SPHinXsys simulation for saving computational cost, hence the arrival time of the solitary wave crest of SPH has been shifted for corresponding to the experiment time. Other parameters of the simulation are: water depth $d=0.6$m and solitary wave height $H=0.2$m. The wave elevation at $x=6$m and $z=1$m (WG in Fig. 9), as well as the horizontal wave force applied on the cylinder, are studied.

In simulations, a fixed particle size is adopted to be $\mathrm{\text{dp}}=0.02$m with the smoothing length of 1.7dp based on the convergence analysis in Sec. 3.1. With this, the numerical domain contains a total of 3100272 particles. The computational time of simulating 1s physical process is 2.55h based on a workstation equipped with an AMD Ryzen Threadripper PRO 5975WX CPU (32-Cores 3.60GHz).

As shown in Fig. 10, the maximum $\eta /d$ obtained by SPHinXsys is 0.348 and the experimental $\eta /d$ is 0.338, with the relative error being 3.05%. The time history of the solitary wave obtained by SPHinXsys agrees well with the experiment result. The perspective view of the wave profiles and the top view of the velocity fields at typical time instants are presented in Figs. 11 and 12.

It can be observed from Fig. 11(b) that the wave climbs up at the upstream side of the cylinder and a wake region is formed at the downstream side at $t=4.4$s. Due to the obstruction of the cylinder, the fluid velocities upstream and downstream the cylinder are very small. After the wave crest passes through the cylinder, at $t=4.8$s, it can be noted from Fig. 11(c) that two concave shapes of the free surface behind the cylinder are symmetrical about the x-axis of the cylinder, and the concave shape of the free surface gradually disappears and becomes gentle at $t=5.2$s. It can also be seen from Fig. 12 that the change of the free surface flow field is uniform, forming a typical state of flow around the cylinder.

As shown in Fig. 13, the maximum positive force $F/\rho g{R}^{2}d$ obtained by SPHinXsys is 1.147 and the magnitude of the negative force is 0.537, while those in experiment are 0.744 and 0.5007, respectively. These correspond to relative discrepancies of 54.17% and 7.25%, respectively. The discrepancy of the positive force is relatively large, which may be attributed to that the used wall boundary condition introduces numerical errors. The extreme value of the positive force occurs at $t=4.4$s and the pressure contours of this time instant are shown in Fig. 14. Specifically, Fig. 14(a) shows the pressure distribution of the 3D view, as well as Figs. 14(b) and 14(c) depict the pressure contours in Y–Z and X–Y planes which are obtained by extrapolating particle pressures onto the grids of a particular section. In general, the particle distribution and pressure distribution are smooth in the whole domain.

4. Conclusions

This paper utilizes the weakly compressible SPH code SPHinXsys to simulate the solitary wave interaction with a submerged horizontal cylinder that can represent the subsea pipeline and a surface-piercing vertical cylinder that mimics the pier of coastal bridges and the column of oil/gas platform. Through the benchmarking by published experimental data, the capability of SPHinXsys in predicting large-amplitude wave elevations and wave profiles, the velocity distribution and the wave-induced force, is demonstrated.

As for the physics behind the wave–cylinder interaction process, in the case one of horizontal cylinder, two vortexes rotating in opposite directions are generated behind the cylinder, which can cause wave forces acting in opposite directions on subsea pipeline. In the case of vertical cylinder, a wake region is formed and two concave shapes of the free surface behind the cylinder are also captured. Regarding the wave forces applied on the cylinder, those induced by steady flow fields are captured satisfactorily by SPHinXsys, while when the flow velocity changes significantly, pressure fluctuations exist and the predicted force show discrepancies.

Acknowledgments

The authors would like to express their sincere gratitude to the anonymous reviewers for the insightful and constructive comments that definitely helped in improving the paper. This research was supported by the Science Foundation of Donghai Laboratory (No. DH-2022KF0311) and the start-up funding provided by Zhejiang University to Min Luo.

ORCID

Guozhen Cai  https://orcid.org/0009-0001-0847-0977

Chi Zhang  https://orcid.org/0000-0003-4072-580X

Yafei Yang  https://orcid.org/0009-0005-4477-7844

Min Luo  https://orcid.org/0000-0002-6688-9127