The loading on a vertical cylinder in steep and breaking waves on sheared currents using smoothed particle hydrodynamics

Combined wave–current fields exist in a wide range of nearshore and coastal regions suitable for the deployment of offshore renewable energy (ORE) systems. Wave characteristics, such as the wavelength, wave-induced velocities, and wave heights, are modified in the presence of a current,^1–3 which results in complex wave–current conditions that determine the loads on offshore systems. Focused waves are often selected as representative for reproducing extreme wave conditions to investigate extreme loads on offshore systems. Meanwhile, a current with a non-zero vorticity is more realistic in coastal regions, and a vertical cylinder is widely used as the substructure of typical offshore fixed wind turbines, which are mainly arranged in nearshore and coastal regions. Thus, prediction of loads acting on a vertical cylinder in combined focused waves and sheared currents conditions is of great significance for practical applications.

Interactions of waves and currents with a cylinder have been studied analytically, numerically, and experimentally. Key findings from studies are focused on the wave–current interaction process, load coefficients (drag and inertia), and structural loadings in combined wave–current conditions. The study of Lin and Li⁴ for a bottom-mounted vertical square cylinder subject to combined wave–current loading indicated that the strength and frequency of vortex shedding caused by the uniform current can be reduced in the presence of waves owing to nonlinear wave–current interaction. It is also shown that the vortex shedding plays an important role in the prediction of in-line force. Jian et al.⁵ demonstrated that the water run-up on a bottom-mounted vertical circular cylinder becomes higher with an increase in the current speed and the total wave loads acting on the cylinder in combined short-crested wave–current conditions are larger compared with the wave loads in pure short-crested wave conditions. Li and Lin⁶ found that the drag coefficients for a stationary submerged horizontal circular cylinder are affected significantly by currents, while the vertical inertia coefficients in combined wave–current conditions are similar to those in pure waves. It is demonstrated that a linear sum of the wave-induced force and the current-induced force cannot predict the actual force on the cylinder accurately in combined solitary wave–current conditions for a horizontal cylinder near the free surface in Ref. 7 and for a vertical cylinder mounted into the flume in Ref. 8. Studies concentrated on interactions of waves and currents with other types of structures, such as the truss model,⁹ the jacket model,¹⁰ the tidal turbine model,¹¹ the submerged plate,¹² and the jack-up platform,¹³ also indicated the non-negligible effect of wave–current interactions on the loading compared to pure wave conditions or pure current conditions. It can be seen that the wave–structure interaction will be affected by the presence of a current, and the wave-induced force and current-induced force cannot be considered separately for the prediction of the actual force.

Harmonics of wave loads and surface elevations for structure in waves and in combined wave–current conditions were also investigated in previous studies. Fitzgerald et al.¹⁴ presented a general phase-based harmonic separation method for fully nonlinear free-surface elevation and hydrodynamic force on a fixed structure, which uses time histories obtained from focused wave groups with phase shifts of $0°$⁠, $90°$⁠, $180°$⁠, and $270°$ interacting with the cylinder, respectively. Mohseni et al.¹⁵ studied wave–cylinder interaction using mesh-based Computational Fluid Dynamics (CFD) with investigation on harmonics in wave run-up. It was found that wave run-up associated with high harmonics is involved in the scattered wave field around the fixed vertical cylinder for both short and long wave conditions. Chen et al.¹⁶ studied the interactions of focused waves and sheared currents with a vertical surface piercing cylinder using mesh-based CFD. Two approaches for numerical modeling (the direct method and the coupling method) were used with the application of an iterative wave focusing methodology to achieve wave focusing at the location of the vertical cylinder in the presence of sheared currents. It is demonstrated that the four-phase based decomposition method can be used to obtain the harmonic structure of forces acting on a cylinder in combined waves and sheared currents accurately. Feng et al.¹⁷ investigated the higher harmonic wave loads and moments on a bottom-mounted vertical cylinder using the four-phase based decomposition method, where it was found that the discrepancy between experimental results and numerical results existed in the third-order harmonic. Ghadirian et al.¹⁸ studied the effects of wave–current interactions on wave forces on a vertical cylinder. Based on the four-phase based decomposition method on the experimental measurements, it was found that the surface elevation time histories of waves on the opposing current are the most nonlinear, while the surface elevation time histories of waves on the following current and the inline force time histories of waves with no current are the most linear. Jiang et al.¹⁹ studied the nonlinear resonant behavior induced by higher-order harmonics for two boxes in a side-by-side arrangement. It was found that the second-order and third-order harmonic-induced wave resonances affect the horizontal wave forces on the boxes in this arrangement. Mj et al.²⁰ investigated the harmonic structure of wave-induced forces on a vertical cylinder for a range of unidirectional and directionally spread random sea-states, and a key finding showed that the harmonic coefficients are almost identical if fitted as powers of the linear inline force and its Hilbert transform for spread and unidirectional seas. Chen et al.²¹ investigated the harmonic structure of hydrodynamic forces on a fixed ship-shaped floating production, storage, and offloading (FPSO) vessel hull, where it was found that the harmonic structure of the force is weakly dependent on the wave steepness and the FPSO bow diameter, while the higher-order harmonics increase rapidly with the decrease in the FPSO length.

Smoothed particle hydrodynamics (SPH),^22,23 as a Lagrangian and meshless technique for CFD, is ideally suited to solve fluid problems involving highly nonlinear deformation.^24–29 For numerical modeling of wave–current interaction, regular waves interacting with steady currents were studied in Ref. 30 using a circulating current system and in Ref. 31 using a non-reflective open boundary condition. Numerical modeling of focused waves on variable sheared currents was presented in Ref. 32 using open boundaries and a modified damping zone. The combination of open boundaries and a modified damping zone is also implemented in the present study. For interactions of waves and currents with a structure using the SPH method, they are rare compared with those using mesh-based CFD methods and experiments. Using a circulating current system based on the SPH method, Shi et al.³³ and Liu et al.³⁴ investigated the hydrodynamic response of a flexible floating boom and silt curtain under combined wave–current conditions, respectively. However, for considering loading on the structure in offshore and ocean engineering based on the SPH method, most previous studies concentrated on waves only without the presence of a current. Wen et al.³⁵ investigated the regular wave interacting with a vertical cylinder using weakly compressible SPH (WCSPH). Lind et al.³⁶ studied wave forces on a vertical cylinder including regular wave loading and focused wave loading using incompressible SPH (ISPH). Chow et al.³⁷ developed a numerical wave basin that is capable of simulating focused breaking waves interacting with a vertical cylinder using ISPH. However, a current is known to affect wave characteristics and thus wave loading. Therefore, one of the main scopes of this paper is modeling complex and extreme events of combined waves and sheared currents interacting with a vertical cylinder using the SPH method which is ideally suited to highly nonlinear problems, and this study has received no attention to date using the SPH method. The second scope is using the four-phase based harmonic separation method to validate the model's ability at capturing the harmonic structure of the elevations and loading.

The solver used in the present study is DualSPHysics which is an open-source SPH code accelerated on a graphics processing unit (GPU).³⁸ Numerical modeling of waves within DualSPHysics was presented in Ref. 39 using a wavemaker and a damping zone and in Ref. 40 using open boundaries. The two methods are well tested within DualSPHysics. The detailed implementation of open boundaries in DualSPHysics was given in Ref. 41. A numerical wave–current flume for waves on variable sheared currents was developed based on the combination of open boundaries and a modified damping zone acting on the vertical velocity component within DualSPHysics in Ref. 32. Focused waves interacting with various sheared currents were validated with analytical linear solutions. Capasso et al.⁴² have presented a similar combined use of a damping zone with inlet–outlet boundary conditions for simulating a moored wave buoy in waves and currents in DualSPHysics using a third-order reconstruction the wave profile and velocity field at the boundaries.

Various methods are used to achieve the focusing of waves at the specified point in space and time. Linear wave theory is used to calculate the initial phase of each wave component at the inlet to achieve a wave group focusing at a particular location and time.⁴³ An iterative wave generation approach is used to calculate the corrected phases for the wave generator input until all the wave components come into phase at the focal location as used in Ref. 44, and this approach was modified to include the dispersion relation considering waves in the presence of a current with a constant shear in Ref. 45. The iterative wave generation approach was extended to include both the phase correction and the amplitude correction in Ref. 46, and this type of iterative method was applied in Ref. 47 to generate waves over complex bathymetry, in Ref. 48 to generate waves on sheared currents, and in Ref. 11 to generate waves in the presence of very fast currents (0.8 m/s of following current). An improved iterative focusing methodology considering a linearized amplitude spectrum instead of a full nonlinear spectrum was presented in Ref. 49, and wave records at different positions are used for phase and amplitude iterations. This iterative wave focusing methodology was applied in experiments for the generation of focused waves on sheared currents⁵⁰ and in the numerical modeling for the generation of breaking focused waves.⁵¹ 

In this paper, interactions of focused waves and sheared currents with a vertical cylinder are investigated. The problem addressed is inertia dominated with limited effect of diffraction. The time histories of the surface elevation and flow kinematics used as boundary conditions at the inlet of the SPH wave–current flume are pre-computed using another Lagrangian numerical model,^52,53 hereby referred to as the Buldakov et al. model. An iterative wave focusing methodology⁴⁹ was used within the model to replicate the experiments. This approach was referred to as the coupling method by using the outputs of the Buldakov et al. model as the inputs of the mesh-based CFD model in Ref. 16. Compared with the direct method in Ref. 16 using the iterative wave focusing methodology in the 2D mesh-based CFD model, the coupling method has a higher computational efficiency. This is due to the application of the iterative wave focusing methodology in the faster Buldakov et al. model to calculate the required boundary conditions that a smaller 3D computational domain of the mesh-based CFD model can be used for wave-current-structure interactions. In this paper, the numerical results from the SPH model are initially estimated with the Buldakov et al. model for surface elevations and horizontal velocity profile without a cylinder in place to validate the accuracy of the generated wave–current conditions. Then, the numerical results from the SPH model are validated with experimental measurements of surface elevations and horizontal forces on the vertical cylinder. Four phase repeats are used to understand the harmonic structure of the surface elevation and associated loading. A convergence study for the horizontal force on the cylinder in focused wave condition and combined focused waves-sheared current conditions is presented. Finally, wave amplitudes are increased up to, and beyond, the breaking threshold, to model steep and breaking waves on sheared currents interacting with the vertical cylinder.

This paper is organized as follows: in Sec. II, the SPH model is outlined. In Sec. III, the experimental setup, the numerical setup, and the harmonic separation method are presented. In Sec. IV, the results of numerical modeling are compared with the experimental measurements including surface elevations, cylinder loading and harmonic components for model validation. In Sec. V, the numerical modeling of extreme wave–current conditions interacting with the vertical cylinder is presented. In Sec. VI, the conclusions and future work are given.

Open boundary conditions are implemented based on the buffer zones that the physical quantities, such as velocity, density, and surface elevation, can be imposed on the buffer particles or extrapolated from the fluid domain. The detailed implementation can be found in Ref. 41. The imposed physical quantities can be from theoretical solutions, experiments, and other numerical tools. The physical quantities extrapolated from the fluid domain are based on the use of ghost nodes. The open boundary conditions are applied in the present study for the generation of wave-alone and wave–current conditions and for particles to leave or enter the fluid domain.

The dynamic boundary condition (DBC)⁶⁰ is the default boundary treatment within DualSPHysics. Modified dynamic boundary conditions (mDBC) were introduced in detail in Ref. 61 such that the unphysical gap between fluid and boundary commonly observed with DBC can be avoided. The boundary particles are arranged in the same way as in DBC. The density of each solid particle is obtained from the position of related ghost node within the fluid domain using a linear extrapolation. In mDBC, a no-slip boundary condition is approximated with zero velocity applied to the boundary particles and is applied in the present study.

Periodic boundary conditions (PBCs) can be used to describe an infinitely long domain by using a periodically repeating finite domain within DualSPHyics.³⁸ To achieve this, particles located within $nh$ (⁠ $n=2$ in the present study) near an open lateral boundary on each lateral side of the domain are allowed to interact with the fluid particles of the opposite open lateral boundary on the other side of the domain completing the truncated support in a cyclic manner. The PBC is applied in the present study instead of solid side walls for the numerical wave–current flume. The general sketch of boundary conditions used in the numerical modeling of wave–current–structure interaction in the present study is shown in Fig. 1.

Experimental studies on interactions of focused wave and sheared currents with a vertical cylinder were performed in a wave–current flume at University College London (UCL). The experimental and numerical setup are presented in this section. The experimental measurements are used to validate the results of the numerical modeling with the vertical cylinder in place. In addition, the harmonic separation method applied to both experimental and numerical data is also presented in this section.

The recirculating wave–current flume at UCL is 20 m long, 1.2 m wide, and 1 m depth with all experiments conducted in a water depth of 0.5 m. Two Edinburgh Design Limited force-feedback “piston-type” wavemakers are installed at each end of the flume. One wavemaker is used to generate waves, and the other is used to absorb waves actively. The flow entered the working section of the flume at approximately 1 m in front of the wavemaker.

The time histories of horizontal velocities and surface elevation used as the inputs at the inlet of the SPH numerical wave–current flume are pre-computed by the Buldakov et al. model with the application of the iterative wave focusing methodology to reconstruct experimental surface elevation and kinematics of incoming focused waves and incoming focused waves on sheared currents. A truncated numerical flume compared with the physical flume is allowed owing to a full development of the combined wave–current conditions within this distance before interacting with the cylinder. A general sketch of the SPH numerical model is shown in Fig. 3. The domain length is 10 m, the damping zone length is 3 m, the water depth is 0.5 m, and the domain width is 1.2 m. For wave–current–structure interaction tests, the front face of the vertical cylinder (diameter D = 0.25 m) is located at the FP (x = 0 m), which is 8.7 m from the wavemaker in the physical flume and is 5 m from the inlet in the SPH numerical flume. The time of linear focus is defined as $t f$ = 0 s. The force provided by the post-processing software ComputeForces³⁸ provides values for the total instantaneous force including the combined effects of waves and currents.

Open boundaries are applied for the generation of wave-alone and wave–current conditions and for particles to leave or enter the fluid domain. At the inlet, time histories of horizontal velocities and surface elevation from the Buldakov et al. model are imposed, and density is extrapolated from the fluid domain. At the outlet, horizontal velocities are imposed according to the current profile, surface elevation is imposed according to the water depth, and density is extrapolated from the fluid domain. Only horizontal velocities are imposed at the inlet since no accuracy improvement can be gained by imposing vertical velocities and a negative impact on the particle spacing can occur.^40,62 It is suggested that at least eight layers of buffer particles arranged in buffer zones can give accurate wave propagation simulations,⁴⁰ and this value is used in the present study.

Following the methodology presented in Ref. 32, the buffer zone at the inlet is divided into eight vertical sections in the present study. In each vertical section, the horizontal velocities (time series) are imposed at three different heights. The horizontal velocities at other heights in this section are obtained according to a parabolic fit function of DualSPHysics for each instant. A general sketch of the buffer zone at the inlet with the horizontal velocities imposing at the inlet and one section of the buffer zone is shown in Fig. 4.

Two focused wave groups generating with the phase shifts of $0°$ (crest focused waves) and $180°$ (trough focused waves) can be used to separate even and odd components.^64,65 With additional two focused wave groups generating with the phase shifts of $90°$ and $270°$⁠, these time histories of surface elevation (⁠ $S$⁠) and force (⁠ $F$⁠) can be used to achieve a more effective separation. Figure 5 shows an example of surface elevations at the focal location of focused wave groups generated with the phase shifts of $0°$⁠, $90°$⁠, $180°,$ and $270°$⁠.

The implementation of the SPH-based numerical wave–current flume and modeling of wave–current–structure interactions in the model is validated in this section. The wave–current generation method is first validated with the Buldakov et al. model for surface elevation and velocity profile in Sec. IV A to demonstrate the accuracy of generation of wave–current conditions. The numerical modeling of interactions of focused waves and sheared currents with the cylinder is validated with experiments for surface elevation and horizontal force on the cylinder in Sec. IV B to demonstrate the model's capability for modeling wave–current–structure interaction. The harmonic analysis is performed using the method described in Sec. III C, and the results are validated with experiments in Sec. IV C. The convergence study with computational performance is presented in Sec. IV D.

The sheared current profiles imposed in the SPH numerical model and their comparison with the profiles used in the Buldakov et al. model (referred to as LaNM in the figure legends) are shown in Fig. 6. The current profiles in the Buldakov et al. model were approximated according to the current profiles measured in the experiments. The sheared current profile is also imposed at the outlet of the SPH numerical model.

The surface elevations are measured at the AMP (x = −4.7 m) and the FP (x = 0 m). At the FP, the horizontal velocity profile is measured at the focal time. The measured numerical results for three cases of focused waves on a following sheared current, focused waves only, and focused waves on an opposing sheared current are compared with the Buldakov et al. model in Fig. 7, demonstrating the accuracy of conditions generated in the SPH numerical model. Previous particle refinement studies showing convergence for focused waves and currents can be found in Ref. 32. Figure 8 shows the simulation snapshot of the horizontal velocity field at the focal time for focused waves on a following sheared current, focused waves only, and focused waves on an opposing sheared current.

The measured surface elevations from the SPH numerical model are compared with the surface elevations measured at two of the wave probe locations (the AMP: x = −4.7 m, and at the front face of the cylinder close to the FP: x = −0.02 m) in the experiments. The measured horizontal forces on the cylinder from the SPH numerical model are also compared with the measured ones in the experiments in this section. As suggested in Ref. 37 for the simulation of focused waves impacting a cylinder using ISPH, the particles of the cylinder are arranged in terms of radial arrangement instead of a Cartesian arrangement, and both arrangements are shown in Fig. 9.

The measured numerical results for focused waves on a following sheared current, focused waves only, and focused waves on an opposing sheared current are compared with the experiments (referred to as EXP) in Fig. 10, demonstrating the accuracy of modeling wave–current–structure interactions in the SPH numerical model. Figure 11 shows the simulation snapshot of the horizontal velocity field at the focal time with the cylinder in place for focused waves on a following sheared current, focused waves only, and focused waves on an opposing sheared current. Four focused wave groups are generated with phase shifts of $0°$⁠, $90°$⁠, $180°,$ and $270°$⁠. For focused wave groups generated with phase shifts of $90°$⁠, $180°,$ and $270°$⁠, the measured numerical results of the surface elevations at the AMP and at the front face of the cylinder close to the FP and the horizontal forces on the cylinder are compared with the experiments and given in the  Appendix. The velocity fields are shown to be smooth and in good agreement with the experiments. Some very small noise is noted to exist in the velocity and associated pressure fields. In the present study, the coefficient in the viscosity term in the momentum equation is set to 0.001 to limit the effects of artificial viscosity on wave propagation. In addition to artificial viscosity, there are new treatments emerging, which offer the possibility to reduce this unphysical noise including the velocity-divergence error mitigating (VEM) and volume conservation shifting (VCS) schemes in Ref. 66 or an acoustic damper term in Ref. 67. Alternatively, the pressure values can be filtered during post-processing.⁶⁸ 

Using the harmonic separation method described in Sec. III C, the second-order subharmonic, the linear signal and the first three superharmonics (to fourth-order) are obtained for the surface elevation at the front face of the cylinder and the horizontal force on the cylinder based on the numerical results in Sec. IV B and the  Appendix. The second-order subharmonic and the fourth-order superharmonic are separated by frequency filtering. The frequency range is $0< f / f p < 1$ for the second-order subharmonic as also used in Ref. 16 and $3< f / f p < 5$ for the fourth-order superharmonic according to $N f p− f p<f<N f p+ f p (N=4)$⁠. The harmonic components of the free-surface elevation at the front face of the cylinder and the horizontal force on the cylinder are shown in Figs. 12 and 13, respectively. Time histories of the separated harmonics (second-order subharmonic, linear harmonic, second-order superharmonic, third-order superharmonic, and fourth-order superharmonic) for focused waves on three different flow conditions (following sheared current case, no current case, and opposing sheared current case) are given. An overall good agreement between harmonic components obtained from experimental measurements and numerical results is achieved. However, less good agreement is found for second-order subharmonic and third-order superharmonic as in Ref. 16. The second-order subharmonic discrepancy is likely to do with the difference in the generated subharmonic error waves between the experiments and combined numerical modeling approach.

The RMSE values (time histories of horizontal force based on a 10 s time duration from −5 to 5 s), runtime (for 15 s of physical time), and total number of particles (⁠ $N p$⁠) are listed in Table I. Using a finer resolution (0.0125 m) improves the accuracy of the simulation compared with the other resolutions (0.02, 0.025, and 0.03125 m) with the RMSE value decreasing for all the following sheared current case, the no current case, and the opposing sheared current case. The computational performance of DualSPHysics on GPUs for modeling wave–structure interactions and wave–current–structure interactions is of practical importance for research and engineering. The computational times in the present study are on the order of hours. All simulations in this section were computed on an Nvidia GeForce RTX 3080 Laptop GPU.

The model validation in Sec. IV and the  Appendix demonstrates a high accuracy for modeling wave–current–structure interaction in the SPH numerical model in the present study. In this section, the sensitivity of peak force on the cylinder in steep and breaking focused waves with phase shifts of $0°$⁠, $90°$⁠, $180°,$ and $270°$ on different flow conditions is investigated. Wave amplitudes are increased up to, and beyond, the breaking threshold. The methodology for increasing wave amplitude is presented in Sec. V A. The sensitivity of peak force with increasing wave amplitude for different focused wave phases is presented in Sec. V B. The comparison with the Morison equation is given in Sec. V C.

As shown in Sec. IV D, though a finer resolution gives a closer agreement with experimental data, $d p$ (0.02 m) gives a very similar results compared with $d p$ (0.0125 m) with a considerably shorter runtime. Thus, all the results presented in Sec. V B are using $d p$ (0.02 m). We find significantly different breaking threshold steepness values for the three current profiles. The breaking thresholds are identified visually from the SPH simulations and are found between $∑ k i a i = 0.2$ and $∑ k i a i = 0.25$ for crest focused waves on a following sheared current, between $∑ k i a i = 0.3$ and $∑ k i a i = 0.35$ for crest focused waves only, and between $∑ k i a i = 0.4$ and $∑ k i a i = 0.45$ for crest focused waves on an opposing sheared current.

We believe these differences in breaking threshold steepness can be explained partially by the modified ratios of total surface particle velocity (=wave-induced particle velocity + current surface velocity at the surface) relative to the current-modified phase speed. The increase in the breaking threshold steepness in negative shear can be explained by the reduction of the total surface particle velocity (=a negative surface current + current-modified wave-induced particle velocity) relative to the phase speed. Similarly, the positive shear increases the total surface particle velocity relative to the phase speed, and thus, breaking onset occurs at a lower steepness. This is consistent with the conclusion in Ref. 45 that a reduction in wave steepness at the breaking threshold is found proportional to the magnitude of the velocity gradient. However, the contribution of the approximate scaling methodology, and associated errors in the bound wave structure, could also affect these preliminary results. These findings do warrant a significant future study on the effect of shear profile on breaking steepness, which is out of the scope of this paper.

The maximum target values of $∑ k i a i$ are 0.3, 0.35, and 0.45 for focused waves on a following sheared current, focused waves only, and focused waves on an opposing sheared current, respectively, to obtain breaking waves for all focused wave groups with phase shifts of $0°$⁠, $90°$⁠, $180°,$ and $270°$⁠. The maximum value of the surface elevation at the front face of the cylinder (x = −0.05 m) and the positive peak horizontal force on the cylinder are plotted as a function of $∑ a i$ in Figs. 15 and 16. According to Fig. 16, among four focused wave groups generated with the phase shifts of $0°$⁠, $90°$⁠, $180°,$ and $270°$⁠, it was found that the maximum positive peak horizontal force on the cylinder is obtained in focused waves generated with phase shift $270°$ for two flow conditions (focused waves only and focused waves on an opposing sheared current). For focused waves on a following sheared current, focused waves generated with phase shift $0°$ give a slightly larger positive peak horizontal force on the cylinder than focused waves generated with phase shift $270°$ (200.8 N and 196.5 N). This is contrary to that expected from linear theory for an inertia dominated structure where peak forces would be expected for the $90°$ phase shift condition where components are focused in such a way to provide maximum acceleration. Although this appears to be the case for the more linear conditions (lower $∑ a i$⁠), this is no longer the case for larger waves. The results demonstrate extreme sensitivity of maximum force to the wave phase and hence equivalently the relative location of the cylinder to a large wave event, suggesting different phases/relative locations should routinely be considered in design.

To assess this relationship further, we look at the times where maximum forces are recorded, relative to the focal time, for each condition. The times when the cylinder experiences the positive peak horizontal force are presented in Table II. These are shown for the original and the largest $∑ k i a i$ values normalized according to $( t − t f ) / T p$ (linear focal time $t f$ is 0 s). These values further highlight the complex relationship between the phase of the waves and the peak force as the times where peak force are expected, based on the linear theory, for an inertia dominated structure are not what we observe in the outputs. For example, the maximum value of acceleration and hence expected force in a $0°$ focused wave would be expected at around −0.25  $T p$ whereas this is observed at around −0.1  $T p$ for small amplitude waves and −0.2  $T p$ for steep breaking waves. Similarly, the maximum force for a 90° focused wave would be expected to occur at the focal time (0  $T p$⁠) yet values closer to 0.2  $T p$ are recorded. This problem becomes complex as we increase amplitude due to nonlinear dispersive effects along with the effect of bound harmonics on the kinematics and is further complicated by the presence of wave breaking. This is demonstrated in Fig. 17 showing the force time-history for the $90°$ focused wave cases in opposing current as we increase amplitude. A moderate shift in the time of maximum force is observed with amplitude until the largest case, which is probably due to wave nonlinearity. For the largest wave case where there is breaking, the wave prior to the largest crest predicted by linear theory breaks and introduces a large slamming force, explaining the −1.05  $T p$ time recorded for the maximum force in Table II. However, the positive peak horizontal force occurs during a period of a positive slope in the time histories of the free-surface elevation at the front face of the cylinder for almost all cases, which corresponds to a time associated with large particle accelerations. This is broadly consistent with the NewForce approach⁷⁰ and the findings in the experimental investigation of near-breaking and spilling breaking wave forces on a cylinder,⁷¹ yet the precise times are more difficult to predict due to the highly nonlinear waves and the presence of wave breaking. Wave breaking and the associated slamming, energy redistribution, and dissipation can affect this significantly (as indicated in Fig. 17), and it is noted that for the largest $∑ k i a i$ for the $270°$ focused wave without current, the maximum force occurs closer to the time maximum elevation where wave breaking and slamming forces dominate. The complex relationship observed for which phase introduces the largest force and when this force occurs, particularly in breaking waves, warrants further exploration. It may be required to test a wide range of focused wave phases and relative positions in order to truly identify the peak design loads.

In Fig. 18, the simulation snapshots of the horizontal velocity field at the instant listed in Table II are shown (only the largest $∑ k i a i$⁠). It demonstrates the effectiveness of using the SPH method for numerical modeling of wave–current–structure interaction with high nonlinearity and wave breaking. The results are all for a fixed depth with intermediate-depth waves. The force magnification by breaking over non-breaking waves has been shown experimentally (in regular waves without currents) to be dependent on the depth parameter $kd$ with magnification varying from about $×$ 3 in shallow water to negligible in deep water waves,⁷² and the effect of depth parameter will be studied in further work.

According to Tables III and IV, the values of $C M$ are close to 2, while the values of $C D$ vary from 0.41 to 2.27. The inertia term is dominant in the present study, and the effect of the drag term is relatively small so that the values of $C D$ are less important to determine the total force. It is consistent with the investigation of the variability of the inertia and drag coefficients considering the Morison wave force on a fixed vertical cylinder due to irregular waves in Ref. 75. In their study, irregular wave tests were also inertia dominated that it was found that there is a high variability in the constant drag coefficient compared to the inertia coefficient (close to 2) and that the precise value has a negligible effect on the wave force prediction owing to the predominance of the inertia term. The values of $C M$ are close to the potential flow value of 2, and small differences here are due to viscous effects.⁷⁶ 

The time histories of horizontal force using the Morison equation with the values of $C M$ and $C D$ in Tables III and IV are compared with the time histories of horizontal force computed in the SPH model (subtracting the small mean force) for focused wave groups with the phase shifts of $0°$ at the original $∑ k i a i$ and the largest $∑ k i a i$ (the breaking wave cases) in Fig. 19. It is noted that the Morison equation is not recommended to predict the horizontal force due to breaking waves due to slamming as indicated in Fig. 18. The reason for the observed discrepancy in cases at the largest $∑ k i a i$ in Fig. 19 is likely due to the slamming. The inertia terms and drag terms for these cases are compared in Fig. 20. It is noted that the drag-like term becomes more important in breaking wave cases.

For focused wave groups generated with phase shifts of $90°$⁠, $180°,$ and $270°$ for waves on a following sheared current, waves only, and waves on an opposing sheared current at the original $∑ k i a i$⁠, the inertia terms and drag terms for these cases are compared in Fig. 21. It is indicated that the wave phase plays a less important role in determining the contribution of the inertia term and the drag term in the present study as the effect of drag term is small at the original $∑ k i a i$⁠.

Prediction of loadings on substructures of ORE systems in complex and extreme events of combined waves and sheared currents is of great significance for practical applications in offshore renewable energy. Numerical modeling of interactions of focused waves and sheared currents with a vertical cylinder using SPH-based DualSPHysics code has been introduced in the present study. To generate the conditions in the SPH numerical wave–current flume, open boundaries and a modified damping zone are implemented as in Ref. 32. The time histories of the surface elevation and flow kinematics used as boundary conditions at the inlet of the SPH numerical model are obtained using the Buldakov et al. model with an iterative wave focusing methodology to replicate the wave conditions and wave–current conditions generated in the physical flume. Additionally, a relatively long distance between the wavemaker and the cylinder in the physical flume is not required for the distance between the inlet and the cylinder in the SPH numerical model.

The numerical results from the SPH wave–current flume are validated with the Buldakov et al. model for surface elevation and velocity profile, demonstrating the accuracy of the method used to generate focused waves on different flow conditions (focused waves on a following sheared current, focused waves only, and focused waves on an opposing sheared current). Agreement between the numerical results from the SPH wave–current flume and experimental measurements for surface elevations and horizontal force on the cylinder is achieved, demonstrating the model's capability in the present study for modeling wave–current–structure interactions. Four phase repeats are used in the SPH model to understand the harmonic structure and extract the harmonic components of the surface elevation at the front face of the cylinder and associated loading on the cylinder. Harmonic components and amplitude spectra from harmonic analysis based on the numerical results and the experimental measurements show an overall good agreement. At the same time, a satisfactory computational performance is gained for 3D simulations by only using a laptop GPU, which is of great importance for practical applications in engineering.

The SPH-based numerical flume is well suited for modeling highly nonlinear fluid–structure interaction problems including wave breaking, and the peak forces on the cylinder in steep waves on sheared currents were investigated by increasing the wave amplitudes up to the breaking limit. Among the four focused wave groups with the phase shifts of $0°$⁠, $90°$⁠, $180°,$ and $270°$⁠, it was found that the maximum positive peak horizontal force on the cylinder is obtained in focused waves generated with phase shift $270°$ for two flow conditions (focused waves only and focused waves on an opposing sheared current), and with phase shift $0°$ for one flow condition (focused waves on a following sheared current). This is counter-intuitive to what we would expect from the linear theory, and the relationship between peak force and wave phase, along with when the peak force occurs, is found to be complex and warrants further study. An optimization approach was applied to calculate the inertia coefficient $C M$ and drag coefficient $C D$ in the Morison equation using kinematics from the SPH model. The comparisons indicate that the inertia effect is dominant and the effect of drag effect is relatively small in the present study. However, slamming is shown to be important for breaking wave cases, and the Morison equation is shown not to be able to capture the peak loads for these cases.

The authors would like to acknowledge funding for this project from the School of Engineering at the University of Manchester. S. Draycott acknowledges the Dame Kathleen Ollerenshaw Fellowship.

The authors have no conflicts to disclose.

Yong Yang: Conceptualization (equal); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (equal); Validation (lead); Writing – original draft (lead); Writing – review & editing (equal). Peter Kenneth Stansby: Conceptualization (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Benedict Rogers: Conceptualization (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Eugeny Buldakov: Conceptualization (equal); Methodology (equal); Resources (equal); Writing – review & editing (equal). Dimitris Stagonas: Conceptualization (equal); Methodology (equal); Resources (equal); Writing – review & editing (equal). Samuel Draycott: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

Four focused wave groups generating with the phase shifts of $0°$⁠, $90°$⁠, $180°$⁠, and $270°$ are required for the harmonic separation. Here, the measured numerical results of the surface elevations at the AMP and at the front face of the cylinder close to the FP and the horizontal forces on the cylinder are compared with the experiments for focused waves generating with the phase shift $90°$⁠, focused waves generating with the phase shift $180°$ and focused waves generating with the phase shift $270°$ on different flow conditions (following sheared current case, no current case, and opposing sheared current case) and are given in Figs. 22–24, respectively.