Numerical Simulation of Inline Forces on a Bottom-Mounted Circular Cylinder Under the Action of a Specific Freak Wave

Freak waves pose a great threat to the tension-leg platforms (TLPs) and monopile foundations of offshore wind turbines (OWTs), which necessitates comprehensive investigations on the characteristics of freak waves and the wave actions on those offshore renewable energy structures with circular cylinder. The recorded freak wave series “New Year Wave” (NYW) was numerically simulated using the Computational Fluid Dynamics methods. The compensation measure was adopted to effectively improve simulation accuracy. Under the action of the NYW, the inline forces and secondary load cycle (SLC) on a vertical-mounted cylinder, as the classic form for the TLPs and foundation of OWTs, were fully addressed. The simulation results were compared with the empirical formulations and experimental data to reveal the differences and the possible causes. The development of SLC was found to be closely related to the downstream vortex and return flow, which induces the reduction of the wall pressure and thus the inline force. The maximum inline forces vary with the cylinder position relative to the wave peak, and the simulation results reveal that the linear inline forces calculated by Morison formulation may be less than 65% of the total wave forces.


INTRODUCTION
With the utilization of the renewable marine energy and resources, the offshore structures, including the tension-leg platforms (TLPs) and offshore wind turbines (OWTs) have rapidly developed over the past few decades. Cylindrical structures, e.g., the support structures for TLPs and monopile foundations of the OWTs, are the most common structural forms in ocean engineering, which are easily threatened by extreme freak waves. Freak waves may occur on a calm sea and reach very high amplitude. Klinting and Sand (1987) defined freak waves with the following criteria: (1) H i > 2H s ; (2) H i > H i−1 , H i > H i+1 ; (3) η c > 0.65H i , where H i and H s represent the maximum wave height and the significant wave height, H i−1 and H i+1 are the immediate wave height before and after H i , and η c is the crest height of the maximum wave. Given the difficulty in fully satisfying the above conditions in a real sea environment, Kharif and Pelinovsky (2003) proposed regarding H i > 2H s as the criterion of freak waves. One such example, the "New Year Wave" (NYW, with a maximum wave height of 25 m, and the significant wave height is 10.6 m) in the North Sea, recorded at the Draupner platform on January 1st, 1995 (Haver, 2003), is a well-known representative freak wave. Freak wave is believed to be one major culprit to destruct marine structures and ships (Kharif and Pelinovsky, 2003;Bertotti and Cavaleri, 2008;Kharif et al., 2009), and would pose great threats to the wind turbines, offshore platforms, etc. Therefore, the actions of the freak wave on the marine structures should be discreetly studied to guarantee a safe design.
The freak wave has attracted wide attentions, most of previous studies are limited to moderate waves. For the simulation of the freak wave, it is usually treated as a focused wave, and a large amount of experimental and numerical investigations on the generation of the focused waves have been conducted. Tromans et al. (1991) proposed the NewWave theory, by which the target wave can be effectively generated in the experimental wave tank. Chaplin (1996) proposed the phase convergence method to experimentally generate the focused wave. Based on the wave spectrum energy distribution, Kriebel and Alsina (2000) proposed a method to simulate focused waves and its effectiveness was verified in the wave tank. For the numerical simulation methods, Ning et al. (2008) adopted the quadratic higher order boundary element model to numerically generate the focused waves, and the efficiency of this method was verified. Pakozdi et al. (2016) simulated freak waves by using Star-CCM+, in which the fluid was considered to be inviscid and incompressible. Ghadirian et al. (2016) applied the OceanWave3D potential flow wave model and a coupled OceanWave3D-OpenFOAM solver to simulate two-and threedimensional focused waves, and the results agree well with those of the experiments.
Many methods were proposed to calculate the wave forces (Morison et al., 1950;Rainey, 1995), and some attempts were also made to predict the higher-order wave forces based on the third-order perturbation theories (Faltinsen et al., 1995;Molin, 1995). However, there were limitations to apply those theories to calculate the higher-harmonic wave forces (Rainey, 2007;Paulsen et al., 2014). Experimental and numerical explorations have been extensively performed to reveal the higher-harmonic non-linear characteristics of wave forces induced by the extreme waves. Chaplin et al. (1997) and Bredmose and Jacobsen (2010) experimentally studied the inline forces on the cylinder structures under the action of the freak wave, and found that the maximum wave forces are larger than the results calculated by Morison formulation. Chen et al. (2018) experimentally studied the nonlinear harmonic components of hydrodynamic forces on a fixed cylinder and found the linear forces are only about 60% of the total forces for some cases. Deng et al. (2016) experimentally investigated the freak wave forces acting on a vertical cylinder, and the inline force is found to be closely related to the wave steepness. Liu et al. (2019) numerically studied the freak wave loads on a small diameter cylinder using OpenFOAM. A high slamming force is observed under the spilling breaking waves induced by the hitting of the local flow with high horizontal velocity on the structure. Ning et al. (2017) found the fully nonlinear potential flow theory could well predict the maximum pressure on the structure under the action of focused wave. Kim et al. (2012) found that the inline force can be doubled when the wave elevation increased by 25%. Tang et al. (2020) numerically studied the effects of breaking wave loads on a 10-MW largescale monopile offshore wind turbine, and found the breaking wave created larger inline forces than the traveling wave, and the Morison formulation would underestimate the wave load. Besides, as mentioned by Bachynski and Moan (2014), none of the aforementioned theoretical approaches could capture the secondary load cycle (SLC) that was observed in experimental and numerical studies, and Grue et al. (1993) believed that SLC makes a crucial contribution to higher-harmonic component of inline force, which may cause the ringing vibrations of the marine structures.
To the authors' knowledge, the previous simulations for the specific freak wave "NYW" were not accurate enough, which may cause errors in the simulation of wave forces and its higherharmonic components. In this study, the specific freak wave series, "NYW, " is numerically simulated and further modified. The non-linear inline forces acting on a bottom-mounted vertical cylinder (widely used for the support structures for TLPs and monopile foundations of the OWTs) and the characteristics of the SLC are analyzed. The remainder of the paper is organized as follows. The numerical methodology is introduced in section 2. In section 3, the two-dimensional (2-D) numerical models are used to simulate the "NYW, " and the simulation accuracy is verified. In section 4, the three-dimensional (3-D) numerical simulation of the actions of the "NYW" on a vertical cylinder is presented. The non-linear wave forces and the SLC are comprehensively investigated. In section 5, the influences of wave front steepness (kA) on the non-linear wave force, and the effects of the flow fields on the SLC magnitude are investigated. Finally, the conclusions are given in section 6.

Governing Equations of Fluid
The time-filtered Navier-Stokes equations in the Eulerian frame were applied to describe the incompressible, viscous fluid flow. The continuity and momentum conservation equations can be presented as: where t is the time, and x i and x j (i,j = 1, 2, and 3) and ith and jth Cartesian coordinates,ū i and u i denote the time-averaged and fluctuating flow velocity components, respectively,pdenotes the time-averaged pressure, ρ and υ are the fluid density and kinematic viscosity. f i represents the volume force. The pressureimplicit with splitting of operators (PISO) algorithm was used for the solution of pressure-velocity coupling. The second order upwind scheme was used for the momentum term. The leastsquares cell-based scheme was used for gradient terms.

Free Surface Tracking
The volume of fluid (VOF) (Hirt and Nichols, 1981) scheme was adopted to simulate the two-phase flow and the Geo-reconstruct scheme enables to capture the free surface, and the volume ratio relation between two phases and transport equation can be presented as: where α q is the volume fraction, representing the volume fraction coefficient of the fluid. When the cell is full of the water, the value of the volume fraction is 1, and when the cell is empty of water, the value is 0. Besides, the parameter is between 0 and 1 when the cell interface between the two phases cuts the cell, therefore the volume fraction can be evaluated as:

Wave Generation
The dynamic grid methods were applied to simulate the pistontype wavemaker, and the transfer function (Tr) between wave height (H) and wavemaker stroke (S) can be written as: where k and d are the wave number and still water depth, respectively. For the wave absorbing methods, a source term was introduced into the momentum equation to avoid wave reflection. The damping term was added to the momentum equation to construct the following wave elimination equation: where u and w denote the velocity components in the x and z directions, respectively. µ(x) and µ(z) are the corresponding attenuation coefficients, which can be respectively calculated as: where the empirical coefficients a x and a z in the x and z direction were set as 0.5 and 4, respectively, x 1 and x 2 are the boundaries of the wave absorbing zone.
The "NYW" was recorded at a sampling frequency of 2.133 Hz in North Sea in 1995, and the wave series are shown in Figure 1A.
Considering the wave crest occurred at 264 s, the wave series from 120 to 500 s were selected as the target to save the calculation resources, the number of the total data is 812. In addition, the target wave series were simulated at a scale of 1:100 based on the Froude similarity criterion. Therefore, the sampling frequency of the scaled wave series is 21.333 Hz and the frequency intervals f = 0.0263 Hz. The superposition of the simple harmonics (Longuet-Higgins, 1974) was carried out to generate the focused wave at specified location and time, the wave surface elevations (η) can be calculated as: where k i , a i , and ε i are the component wave numbers, amplitudes and the initial phases, respectively. For the scaled wave series, and the FFT decomposition was performed to obtain a i and ε i , which are shown in Figures 1B,C, respectively. Besides, a i with higher frequency (ω ≥ 30 rad/s) components were almost negligible and only 182 wave components were incorporated for the mathematical formulations of the target wave series.

NUMERICAL SIMULATION OF "NYW"
To improve the simulation efficiency, the 2-D numerical model without the cylinder was firstly built up to simulate the "NYW." The numerical wave tank (Figure 2) is 15 m long, with a 2.5 m length of wave absorbing zone. The still water depth is 0.9 m. The wave surface elevations at 11 positions (arranged with the uniform interval of 0.1 m) were monitored. A mesh convergence study was carried out for the numerical simulation of the "NYW." The mesh arrangements for the 2-D numerical models are list in Table 1.
The target wave series were numerically reproduced at x = 7.3 m and the horizontal velocity of the wavemaker can be presented as: The calculating time step ( t) is preliminarily set as 0.001 s. The simulated waves with stable wave from the Coarse, Medium, and Fine mesh cases, and the target wave surface elevation (η) are shown in Figure 3A, where η c is the crest height, and the stable simulated results were selected to quantify the error. Because the control signal of the wavemaker was based on the linear wave theory which takes no consideration of non-linearity effects and the wave-wave non-linear interaction. Even with the fine mesh model, the errors for the initially calculated H and η c are 14.4 and 22.9%, respectively, and the standard deviation for the η 1 (t)is 0.018, where η 1 (t) represents the differences between η 0 (t) and η 1 (t). The time-step interval convergence analyses were also  conducted, 0.002, 0.001, and 0.0005 s, respectively. The results of wave elevations with three time-steps were shown in Figure 3B.
It is clearly found that the errors between the simulation results and the target values cannot be further reduced by improving the grid quality and decreasing the time-step. The wavemaker control signal required to be modified to eliminate the deviations between the initial simulated wave elevations η 1 (t) and the target data η 0 (t). The modification for the horizontal velocity of the wavemaker can be presented as: where U 1 (t) and U 1 (t) represent the wavemaker horizontal velocities. c is the modification coefficient, which is not only affected by the numerically simulated H and η c , but also the simulated mild random waves. The modification quality criteria could set as: the errors for the simulated H and η c need to be less than 4.8 and 7.6%, respectively, which approximate a third of the errors for the initial simulation. Besides, the standard deviation for should be less than 0.020, which is approximately 1.1 times larger than the standard deviation for η 1 (t). The optimal c can be obtained after several attempts. It's worth noting that the optimal values of c may be different in those three mesh cases. Therefore, the Fine mesh model was selected as the standard case to obtain the optimal c. And the same value was used in the other mesh cases to simulate the wave series. After several attempts, 0.18 was selected as the optimal value and the newly-calculated results for the three mesh cases are shown in Figure 3C. Both H and η c are much closer to the target values than those without using the modification technique. The errors for the newly-calculated H and η c are 2.4 and 7.3%, respectively. In order to further evaluate the effectiveness of the modification method, kth (k = 1, 2, 3, and 4) order statistics of η 1 (t)and η 2 (t)were calculated, which denote the mean value, standard deviation, skewness and kurtosis, respectively. The first four orders statistic of η 1 (t)are 0.0032, 0.018, 0.157, and 3.598, respectively. And the first four orders statistic of η 2 (t) are 0.0024, 0.019, 0058, and 3.239, respectively. The results by using the Medium and Fine meshes are fairly close. In addition, it is apparent that (Figures 3B,D) time-step convergence is reached for t = 0.001 s. It is clear that the modification method can improve the simulation accuracies of H and η c without sacrificing the simulation quality of the mild random wave series. So, the Medium mesh and t = 0.001 s were adopted in consideration of both simulation accuracy and computation time.
Moreover, the wave height, respectively by using the phase iteration (Gao et al., 2015) and amplitude-phase iteration (Deng et al., 2015) methods are 0.240 and 0.236 m, and the wave crest are 0.150 and 0.159 m, the modification results in this study are closer to the target values, by which the efficacy of the modification method can be verified.
In addition, the numerical simulation accuracy can be improved by several iterations, and the results are shown in Figure 4, the mean value and standard deviation of the simulation deviation with the target became smaller with the increase of iterations. However, the modified quality would not be further improved in the fifth iterative simulation. The wave elevations of the first-time modification result monitored at eleven positions are shown in Figure 5A, where 0 m corresponds with x = 7.3 m, A represents the difference between η c and the wave front trough, which varies from the traditional definition of the wave height (H). The development process of the freak wave can be observed. η c and A range from 0.145 to 0.177 m, and 0.219 to 0.266 m, where A is the wave front height. The maximum η c and A occur at x = 0.3 and −0.1 m, respectively.
The wave front steepness kA at different locations are shown in Figure 5B, where k is the wave number calculated by the troughto-trough period of the freak wave. The maximum and minimum kA are 0.296 and 0.275, which occur at x = 0.4 and −0.2 m, respectively. Overall, the kA in this region are very close.

ACTION OF THE "NYW" ON A BOTTOM-MOUNTED VERTICAL CYLINDER Empirical Formulations of Inline Forces
Morison formulation is used to calculating the inline force on the fixed unit length vertical cylinder, which can be presented as: where f M is the wave force acting on unit length body, ρ is the water density, u x is the horizontal water particle velocity, C D , C M are the empirical drag and inertial coefficients, which are selected as 1.0 and 2.0, respectively. However, the Morison formulation can be used only for the calculation of the linear harmonic force and it cannot be used for prediction of the high-order non-linear wave forces. In order to refine the Morison inertia term, the formulations for the "axial divergence force" and "surface intersection force" on a slender body were proposed by Rainey (1995), which brings the accuracy up to the second-order in Stokes expansion. The corresponding correction terms to calculate the inline forces can be written as: where the cylinder submerged depth (d s ) and wave surface elevation (η) are the lower and upper integral bounds for the "axial divergence force, " respectively. Eq. 11a is the second order in Stokes expansion and Eq. 11b represents a point force that exists at the instantaneous free-surface elevation due to the entry and exit forces. The summation of F M , F A , and F I can be denoted as F MAI , i.e., F MAI = F M +F A +F I .

3-D Numerical Simulation of the Inline Forces Induced by the "NYW"
The aforementioned wave generation method of the "NYW" was applied to the 3-D numerical wave tank, and modified the wavemaker control signal for one time. The cylinder is located at x = 7.3 m, and the diameter (D) and submerged depth (d s ) of the cylinder are 0.1 and 0.3 m, respectively. For the 3-D numerical model, the grids around the cylinder in the x-y plane are shown in Figure 6, where N C , N L , and N Y are the cell numbers around the cylinder, body-fitted cell layers, and Y-axis grids, respectively. The simulation need to avoid the buffer layer while determining the y + ( = u * y υ , where u * is the friction velocity; y is the wall spacing, and υ is the dynamic viscosity.), which in other words is to guarantee either y + < 5 (viscous sublayer) or 300 > y + > 30 (logarithmic sublayer), and y+ were kept inside this region (300 > y + > 30) for all the simulations in the present study. The maximum dimensionless wall distance y + of the first layer grids was lower than 100.
In order to further study the influences of the mesh arrangement around the cylinder on the evaluate inline forces accuracy, a mesh convergence study on the 3-D numerical model was carried out, and the grid parameters are listed in Table 2. All of the cases in this study were computed on a workstation with 4 × Intel(R) Xeon 6148 CPUs.
The actions of the "NYW" on a bottom-mounted vertical cylinder with different grid numbers are shown in Figure 7A. Generally, the results from Medium and Fine meshes are fairly close, and the Medium (3-D) mesh is adopted for the following analysis in consideration of both simulation accuracy and computation time. Based on the numerically simulated wave  elevations, the FFT decomposition was performed to obtain the related linear wave components, including the component wave numbers, amplitudes and the initial phases. For an irregular wave series, a common practice is to superimpose water velocity and acceleration of each wave component. The F M , F MAI , and experimental results (Deng et al., 2016) are also shown in Figure 7A for comparisons. Some findings can be summarized as follows: (1) Because the Morison formulation cannot accurately calculate the high-order non-linear wave forces, the maximum F M is smaller than those of the numerical and experimental cases. F MAI improves the accuracy to the second-order in Stokes expansion, and the peak result is closer to the numerical value. As shown in Figure 7A, both F M and F MAI deviate from the numerical and experimental results with SLC, it should be stated that the SLC stage can be predicted by neither F M and nor F MAI .
(2) The numerically simulated inline forces histories are fairly close to the experimental results. The differences are mainly due to the different target wave elevations. The target in the experiment was derived from the JONSWAP spectrum, which unavoidably deviates much or less from the actual "NYW" series. H and η c in the experiment are also smaller than those of the numerical simulation.

Analyses of the SLC
The SLC, which is an additional longitudinal short-duration force prior to the minimum wave force, may appear when a steep crest passes the vertical cylinder, and it was firstly identified by Grue et al. (1993). The wave forces that calculated by the theoretical methods (e.g., Morison formulation, slender body theory, and FNV) have some weaknesses in predicting the inline forces, especially for the minimum forces and SLC stage. Compared with the experimental method, the computational fluid dynamics (CFD) methods have some superiorities in revealing the mechanism of SLC of vertical cylinders, because more flow details around the structure can be obtained. For the SLC, F c ( = F c k/ρgπr 2 ) and F c / F p were respectively used by Chaplin et al. (1997) and Liu et al. (2019) to quantify the magnitude of SLC, where F p is the peak-to-peak value of the total inline wave force, F c is shown in Figure 7B. The F c and F c / F p are 0.043 and 0.053, respectively. In addition, the dimensionless parameter T H ( = T H /T) was used to quantify the duration of the SLC, where T H is the duration of SLC, and the T H = 0.152 in this case. Paulsen et al. (2014) numerically studied the SLC of bottom-mounted cylinders under regular waves conditions, the downstream vortex was observed in the SLC cases. However, the detailed relationship between the downstream vortex and the SLC was not unveiled. In this study, six typical moments were selected, which were marked as A-F in Figure 7B, to investigate the effects of the vortices on the SLC. The instantaneous flow fields in the slice of z = 0.75 m are shown in Figure 8A. As the vortex structures in the SLC stage remain symmetric, the wall pressure distributions also keep symmetry along the wave propagating direction, so the side-wall pressure distribution is shown in Figure 8B.
Two symmetric vortices appear at moment A, the wall pressure shows visible changes near the position of the downstream vortices, the non-smoothed pressure contour line indicates that the wall pressures decrease in those areas. Then the downstream vortices become larger and the influences on the side-wall pressure are more notable. The influences of the vortical structures on the wall pressure becomes more significant with decreasing submerged depths. The inline force continues to decrease with the rising of vortex size, and the inline force is slightly upward during the moment A-C. The downstream vortices become largest at moment E and move from the downstream to the lateral sides of the cylinder. At the same time, the nonsmoothed pressure contour line also moves to the side of the cylinder. At moment F, the vortices structures disappeared and the pressure changes become less pronounced. According to the relationships between the development of SLC and the vortex processes, it can be inferred that the downstream vortices lead to the decrease of wall pressures, which may be the main cause of SLC. Figure 9 presents the flow velocity distributions along the z-direction with 0.5r, r, and 1.5r (r is the radius of the cylinder) downstream the cylinder. At Moment A, the return flow was observed at position 0.5r, the velocity changes more obviously with the increase in the z coordinate, which induced the wall pressure changes more notably near the water surface. Then the range of the return flow area becomes larger with the development of SLC, and the return flow was observed at positions 0.5r, r, and 1.5r. As shown in Figure 9, it indicates that the return flow is induced by the downstream vortex, and larger vortices will produce more distinct effects on the flow velocity, which finally results in the SLC (additional inline forces). According to the analysis of the numerical computations, the SLC may be induced by downstream vortex associated with the changes of the wall pressure and return flow.

INLINE FORCES WITH DIFFERENT CYLINDER POSITIONS
In order to study the effects of kA on the maximun inline forces and the magnitude of the SLC. Seven more cases with different cylinder locations, i.e., x = −0.2, −0.1, 0.1, 0.2, 0.3, 0.4, and 0.5 m were calculated. The time histories of inline forces are shown in Figure 10. The related maximum inline forces are listed in Table 3, where the dimensionless parameter kA could be obtained from Figure 5B.
It found that the F CFD /F M and F CFD /F MAI increase with increasing the kA, which indicates the occupation of the nonlinear components correspondingly become larger. According to the results of F CFD /F M , the ratio of linear inline force to the total is less than 65% when kA equals to 0.296. For the inline forces induced by the freak waves, such as the "NYW, " both F M and F MAI would underestimate the maximum inline forces that contain the higher-order components and slamming forces. Besides, the maximum forces occurred before the wave crest and the delta-T among the range from 0.08 to 0.11 s for the eight cases, which are close to the experimental results (Deng et al., 2016). Based on the numerical results, the moment of the force peak is not the instant when the still water passes through the central axis of the cylinder. Therefore, it is inappropriate to predict the moment of the maximum wave force through the regular wave theory.
As shown in Figure 11, F c and F c / F p vary from 3.3 to 4.3%, and 4.3 to 5.4%, respectively. F c decreases with increasing kA when kA is larger than 0.276. The dimensionless parameter T H varies from 0.127 to 0.140, which show fairly satisfactory agreements with the values (T H < 0.150 and <0.173) that were respectively monitored by Grue et al. (1993) and Liu et al. (2019).  Considering the SLC is a complex phenomenon, the kA in the aforementioned 3-D cases were relatively close, more focused wave cases were numerically simulated to further study the influence of the flow fields on the wave forces. The focused waves for C745, C765, and C775 in the experiments of Chaplin et al. (1997) were simulated. The related parameters of the three waves could be obtain from the descriptions for the focused waves (Chaplin et al., 1997). The still water depth was 0.525 m and the duration between troughs was close to 1.1 s.
As mentioned in the experiment, the focused wave series was generated by the superposition of 34 simple harmonic wave components, the frequencies of which ranged from 0.511 to 1.244 Hz, with a step of 0.0222 Hz. Same with the generation of the NYW, the superposition of the simple harmonics (Longuet-Higgins, 1974) was taken as the target of the numerical simulation. The numerically simulated wave series are shown in Figure 12A, and the experimental results are also offered for comparison. The simulated wave height and wave crest coincide well with those of the experimental results.
Based on the experiments conducted by Chaplin et al. (1997), the vertical cylinder was located in the center of the working zone, which the diameter is 100 mm. The numerical and experimental overturning moments are shown in Figure 12B. Overall, the two sets of histories match fairly well, especially for the maximum moments. For the cases with smaller wave steepness, the F M and F MAI could predict the maximum overturning moments accurately, but both the Morison formula and the slender body theory would underestimate the maximum overturning moments when the wave steepness becomes larger. Similar to the results in Figure 7A, the maximum M M is smaller than those of the numerical and experimental cases. M MAI improves the accuracy to the second-order in Stokes expansion, and the peak result is closer to the numerical and experimental values. In addition, the SLC appeared for C765 and C775, whereas it was not observed for C745. Those existing theoretical formula could not predict SLC stage, because those formulas considered the fluid is idealized and assumed to be irrotational and non-viscous, therefore, the suction forced induced by the downstream vortex could not be monitored. Moreover, the CFD method could predict the SLC stage satisfactorily, this also verified the numerical simulation accuracy.
The horizontal velocity streamlines at Moment E are shown in Figure 12C, in which the velocity streamlines are in a horizontal plane perpendicular to the axis of the cylinder for z = 0.425 m. For C765 and C775, two symmetric vortices were generated on the near-wall surfaces of both sides of the cylinder during the SLC stage. In comparison, the vortices in C745 were FIGURE 11 | The relationship between F c , F c / F p , and kA.
Frontiers in Marine Science | www.frontiersin.org much smaller. The corresponding magnitudes of the SLC, F c of C745, C765, and C775 were 0, 0.044, and 0.059, respectively. Thus, the vortex size is one important factor affecting the magnitude of the SLC.

CONCLUDING REMARKS
The "NYW" was numerically simulated based on the 2-D CFD model and the modification method was applied to improve the simulation accuracy. In addition, the action of the "NYW" on a bottom-mounted vertical cylinder was numerically simulated based on the 3-D model. Sensitivity analysis for the grid resolution was conducted to ensure the computational accuracy, and the simulation results were compared with the experimental and empirical formulation results, the differences were analyzed to unveil the possible causes. Finally, the influences of kA on the maximum inline forces were analyzed. Some major conclusions are drawn as follows: (1) For the simulation of the extreme freak wave, e.g., "NYW, " the simulation accuracy can be improved by the compensation for the difference between the target wave series and initial numerical results. (2) Under the action of "NYW, " the cylinder maximum and minimum inline forces including the high-order nonlinear components cannot be accurately predicted by the Morison and the slender body theoretical formulations. In addition, the pressure contour line and distribution of velocities at six typical moments during the SLC stage shows that the appearance of SLC is closely related to the downstream vortex, the vortex size is one important factor affecting the magnitude of the SLC. (3) The non-linear components of the inline force increase with increasing kA. The numerical results reveal that the linear inline force may be less than 65% of the total wave loading. The Morison and slender body theoretical formulations would underestimate the maximum inline forces that contain the higher-order components andslamming forces. This paper provides a useful reference for the design of the support structures for TLPs and monopile foundations of the OWTs under the actions of the freak waves.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
YW: methodology, formal analysis, data curation, and writingoriginal draft. FX: conceptualization, supervision, resources, writing-review and editing, and visualization. ZZ: software and validation. All authors contributed to the article and approved the submitted version.

FUNDING
This work is supported by the National Science Foundation of China (51978130). The authors would like to sincerely thank Sverre Haver from Statoil for providing the New Year Wave series.