Numerical Comparison for Focused Wave Propagation Between the Fully Nonlinear Potential Flow and the Viscous Fluid Flow Models

Numerical simulations on focused wave propagation are carried out by using three types of numerical models, including the linear potential flow, the nonlinear potential flow and the viscous fluid flow models. The wave—wave interaction of the focused wave group with different frequency bands and input wave amplitudes is examined, by which the influence of free surface nonlinearity and fluid viscosity on the related phenomenon of focused wave is investigated. The significant influence of free surface nonlinearity on the characteristics of focused wave can be observed, including the increased focused wave crest, delayed focused time and downstream shift of focused position with the increase of input amplitude. It can plot the evident difference between the results of the nonlinear potential flow and linear potential flow models. However, only a little discrepancy between the nonlinear potential flow and viscous fluid flow models can be observed, implying the insignificant effect of fluid viscosity on focused wave behavior. Therefore, the nonlinear potential flow model is recommended for simulating the non-breaking focused wave problem in this study.


Introduction
With the offshore industry moving towards deeper water and undergoing harsher environment, increasingly accurate and efficient numerical models for the interaction of extreme waves with the structure are attracting more and more attention. As is well-known, ocean and offshore structures must be designed to survive in very hostile environment, generation of extreme waves in the experimental laboratory or the numerical code is an important aspect for hydrodynamic analysis. Primarily, only regular or random waves can be adopted for simulations. The regular wave cannot represent the actual extreme event, while the random wave takes too long time to achieve that for its quite low probability of occurrence in any random wave series. In order to overcome these disadvantages, the focused wave group is usually used to simulate the extreme hydrodynamic condition. According to adjusted phase shifts, wave components in a spectrum can focus simultaneously at a position in space, and therefore the extreme wave profile can be expressed equivalently with a random process generated with a specified wave energy spectrum.
Experiments and numerical simulations for focused wave groups have been studied by many researchers for generating the transient focused wave event in recent decades. Baldock et al. (1996) conducted an experimental investigation for producing the large transient wave group at the specified point in space and time domain by using the linear superposition of regular wave trains. Johannessen and Swan (1997) generated two-dimensional extreme water waves in an irregular sea state based on an appropriate description of underlying wave spectra. Borthwick et al. (2006) presented the water particle kinematics of focused wave groups at normal and 20° incidence to a 1:20 plane beach. A laboratory study of the focusing directionally spread surface water waves was conducted by Johannessen and Swan (2001), in which a large number of waves with varying frequency and different propagating directions had been used to generate the extreme wave. Li et al. (2012Li et al. ( , 2014 performed an experiment in the wave basin to investigate the interactions between multi-directional focused wave and vertical bottom-mounted cylinder. Chen et al. (2018) measured nonlinear wave forces on the surface-piercing column exerted by the focused wave. According to the extensive experimental investigation, the technique of focused wave group has been adopted to simulate the extreme wave load on marine structures.
With the development of the computational technique, the numerical wave flume has been developed as an effective and efficient tool for modelling various water wave problems. Zang et al. (2006) examined the effect of the second-order nonlinearity on focused wave behavior and the interaction with a fixed ship hull. Bai and Eatock Taylor (2007) simulated a fully nonlinear focused wave propagation and diffraction problem based on a fully nonlinear potential flow model. Ning et al. (2008Ning et al. ( , 2009) simulated the propagation of nonlinear focused wave groups by using a nonlinear boundary integral equation solved by a higher-order boundary element method. Westphalen et al. (2012) conducted the generation to investigate the behavior of extreme focused wave groups in a numerical wave tank, in which the nonlinear effect of the extreme wave was discussed. According to the comparison, the evident discrepancy between the linear potential solution and the laboratory test can be observed, indicating that the accurate description of an extreme wave event must incorporate the influence of free surface nonlinearity. Recently, the Naiver− Stokes based CFD solver is also employed for studying the focused wave problem. Bihs et al. (2017) investigated the focused wave generation, kinematics, and its interaction with a vertical circular cylinder using the open-source CFD model REEF3D. Vyzikas et al. (2018) studied the evolution of unidirectional wave groups by using an experimental measurement and a Reynolds Averaged Navier−Stokes (RANS) simulation based on OpenFOAM package, in which the evolution of free and bound waves during dispersive focusing was investigated. Stagonas et al. (2018) simulated breaking focused waves by using a numerical wave flume with a piston-type wave maker built in the CFD model olaFlow. It can be found that the nonlinear potential flow model cannot simulate the breaking phenomenon of focused wave motion, leading to the inaccurate results.
However, to the best of the authors' knowledge, there are quite few jobs that have been done to make a detailed comparison of numerical results between the nonlinear potential flow and viscous fluid flow models for non-breaking focused wave. A comprehensive investigation on the influence of the free surface nonlinearity and the fluid viscosity is also helpful for further understanding the behavior of fo-cused wave propagation. In this work, the behavior of focused wave group is simulated by using the linear potential flow, nonlinear potential flow and viscous fluid flow models. The wave-wave interaction of the focused wave group is investigated, by which the influence of the free surface nonlinearity and fluid viscosity on the related phenomenon of focused wave is discussed. This is also the primary objective of this study.

Governing equations and boundary conditions
Φ Ω In order to simulate a three-dimensional wave−wave interaction problem, a right-handed Cartesian coordinate system Oxyz with the origin O in the mean water surface and zaxis pointing upwards vertically is defined as shown in Fig. 1. Based on the potential flow theory, the fluid flow can be formulated in terms of a velocity potential (x, y, z, t), which satisfies Laplace's equation within the fluid domain , and various boundary conditions on all surfaces S of the fluid domain.
On the instantaneous free surface S F , the kinematic and dynamic boundary conditions in the Lagrangian description can be expressed as: where D/Dt is the material derivative, X denotes the position of water particles on the free surface, t is the time, and g is the acceleration of gravity. The boundary condition on the wave maker can be given as: where U(t) is the velocity of the wave maker at the instantaneous position, and n = (n x , n y , n z ) is the normal unit vector pointing out of the fluid domain. On solid surface boundaries, including the bottom and sidewalls, the no-flux condition can be imposed as: Since a time domain approach is applied to formulate the present problem, an initial condition is also required. For every case considered here the water surface is initially calm, so that, In addition, a proper far-field condition should be implemented to avoid the wave reflection from the downstream end of the numerical wave flume, which will be defined and discussed in Section 2.3.

Higher-order boundary element formulation
The boundary integral equation for the potential flow problem can be obtained through Green's second identity, where C(x 0 ) is the solid angle at the source point x 0 , and n is measured from the field point x. The solid angle C(x 0 ) can be computed by ξ η A Higher-Order Boundary Element Method (HOBEM) based on quadratic isoparametric elements is adopted to discretize the integral surface. After introducing shape functions N j ( , ) in each surface element, the potential and its derivatives within an element can be expressed in terms of nodal values as follows: Furthermore, the discretized equation of Eq. (6) can be written as: where, N n and N p are the numbers of elements on the Neumann and Dirichlet boundaries, respectively. The terms in the coefficient matrix A and the right-hand side vector B, including the detailed information about the quadratic isoparametric elements are given in Eatock Taylor (2006, 2009).
In order to accelerate the calculation, the domain decomposition technique through implementation of the D/D-N/N iterative scheme is introduced into this model. By using this method, the number of unknowns in each subdomain becomes much smaller than that of equations in the original whole domain. As the memory and the computational time for assembling the coefficient matrix required by the boundary element method (BEM) depend quadratically on the number of nodes, the total required memory and computational time can be greatly reduced. Therefore, the domain decomposition method is an efficient algorithm for large-scale numerical simulations. The detailed information of the domain decomposition method can be referred to Bai and Eatock Taylor (2007) and Jiang et al. (2018a).

Numerical implementation
A number of aspects related to the detailed performance should be considered in applying the numerical simulation, which has been clarified in many literatures, such as Bai et al. (2014aBai et al. ( , 2014b, Zhou et al. (2013Zhou et al. ( , 2015, and so on. In this paper, only a brief summary is to be provided in this subsection.
In the present numerical method, two types of elements in the computational domain are equipped: the structured quadrilateral and unstructured triangular meshes. The structured 8-node quadrilateral meshes are distributed on both the vertical wall surface and the wave maker. The unstructured 6-node triangular meshes are distributed on the free surface, which are generated by using Delaunay triangulation method. A double or triple node is employed on the intersections between the free surface and vertical surfaces, including solid surfaces and interfaces among subdomains. This multiple-node technique can determine intersection lines conveniently without the additional treatment: only these nodes in the free surface should be updated, and the other multiple nodes on vertical surfaces are allowed to move with them. A typical discretization generated by the above method for the wave flume is depicted in Fig. 2, in which three subdomains are involved.
An artificial damping layer on the free surface is adopted to absorb the wave energy on the downstream side of the numerical wave flume. It can avoid the wave reflection from far-field computational boundaries and simulate a sufficiently long duration in a reasonably sized domain. In this numerical beach, the kinematic and the dynamic boundary conditions in Eqs. (2a) and (2b) are respectively modified by introducing a damping term over a finite length within the free surface as follows: υ where (l) is a damping coefficient only considered in the damping layer. X e = (x e , y e , 0) is a reference value specifying the at-rest position of the fluid particle.
The free surface geometry and potential are updated by the standard 4th-order Runge-Kutta scheme, and a cosine ramp function is used to modulate the boundary condition on the wave maker during the initial time steps. Lapack iterative scheme is adopted for solving the full and asymmetric influence matrix arising from the mixed boundary value problem. The matrix factorizations can be computed once for two types of coefficient matrix corresponding to the problems with Neumann and Dirichlet boundary conditions, respectively, which can be used on all the steps of the iteration, and subsequently in the back substitution. In order to realize a long-time simulation, mesh regeneration on the free surface is required to mitigate the sawtooth numerical instability, and this is implemented by adopting Laplacian smoothing technique to obtain new nodes on the free surface. Variables at new nodes are obtained from the interpolation with those of old ones based on their given horizontal coordinates.

Brief description for the viscous fluid flow solver
In addition to the nonlinear potential flow model, a viscous fluid flow solver based on OpenFOAM package is also adopted for numerical simulations in the present study, where Naiver-Stokes equations are employed as governing equations of incompressible two-phase flows. Volume of Fluid (VOF) method (Hirt and Nichols, 1981) is adopted to capture the free surface motion, especially for the largeamplitude focused wave. The relaxation zone by Jacobsen et al. (2012) is utilized to eliminate the transmission wave at the outlet boundary. A piston-like paddle is equipped at the starting of wave flume to generate the incident wave. In numerical implementation, Finite Volume Method (FVM) is adopted to solve the Naiver-Stokes equations and VOF equation. The velocity and pressure are decoupled by Pressure Implicit with Splitting of Operators algorithm (PISO) (Issa, 1986). Euler method is used to discretize the transient term. The convection term and the diffusion term are discretized by Gauss Limited Linear method and Gauss Linear Corrected method, respectively. The basic numerical imple-mentations of CFD simulations in OpenFOAM can be referred to Jasak (1996) and Rusche (2003).

∂u/∂n
The numerical computation always starts from the still state, which means that the hydrostatic pressure and zero velocity are specified as initial conditions. The no-slip boundary condition is satisfied on solid walls including the wave maker and seabed. At the upper boundary of the numerical wave flume, a reference pressure p = 0 and a velocity condition are implemented with the outward unit normal vector n. The interface tension between air and water phases is neglected in the present study since the dynamic effect from the air phase is very small. At the end of the spongy layer, zero velocities are applied considering that the waves are damped out there by the spongy layer. More detailed information of the formulations and numerical schemes described above have been given in the previous work of Jiang et al. (2018bJiang et al. ( , 2019aJiang et al. ( , 2019b.

Focused wave generation
A linear wave solution expressing the required surface elevation generated at the wave maker can be given as: ω n where N is the total number of wave components, A n , k n and are the amplitude, wave number and angular wave frequency of the n-th wave component, respectively, and x p is the position of the wave maker obtained by the linear prediction. Based on this equation, the linear wave amplitude at the focus position can be obtained as: η where (x p , t p ) is the wave crest elevation of focused wave at the focal time t p and the focal point x p . A i is the input wave amplitude based on the linear potential flow solution. Noted that the above x p and t p in Eqs. (12) and (13) are derived by the linear potential flow model, which are dependent of the input wave amplitude A i. ω n A piston-like paddle is adopted as the wave maker for generating the surface wave at the starting of wave flume. The relationship between the paddle motion and the free surface elevation can be described by using a transfer function T n ( ) as: ω n ω n ω n ω n where A n ( ) is the free surface amplitude measured at the focused point, and a n ( ) is the amplitude of paddle undergoing the sinusoidal motion at the corresponding wave frequency . In this work, the transfer function T n ( ) is obtained analytically based on the first order wave theory (Li et al., 2008), that is, where d is the water depth of wave flume.

Numerical analysis for focused wave group
In this section, the fully nonlinear potential flow and viscous fluid flow numerical wave flumes established above are employed to consider the highly nonlinear wave-wave interactions problem. An experimental measurement by Baldock et al. (1996) for focused waves is adopted for validation and investigation, where the ranges of the input period and amplitude are tabulated in Table 1. The total length of the wave flume L = 10 m with water depth d = 0.70 m is equipped, and the wave maker is situated at x = -5 m. Based on the characteristics of symmetry, only half of the wave flume is considered, and the half-width is taken as 0.14 m in fully-nonlinear potential flow simulations. As for the viscous fluid flow model, a two-dimensional wave flume with 1.00 m high is simulated, where the water depth and air thickness are 0.70 and 0.30 m, respectively. With Case B and Case D in Baldock et al. (1996), two kinds of wave groups of broad-band and narrow-band spectra are adopted, in which their range of the wave periods are from 0.6 to 1.4 s and from 0.8 to 1.2 s, respectively. Numerical simulations all begin at t = -15 s. The focus position and focus time should be x = 0 and t = 0, respectively, according to the linear wave solution in Eq. (12).
The convergence verification of the computation is carried out for two kinds of focused waves, Case B and Case D mentioned above for A i = 0.055 m. Each group consists of 29 individual wave components, which are of the same amplitudes equally spaced within the appropriate period range. In the fully nonlinear potential flow model, there are 12 sub-domains applied in numerical simulations and three different mesh schemes in each subdomain are tabulated in Table 2. Six elements on the vertical side walls of wave flume are chosen to perform the calculation for all the nonlinear potential flow cases, and the time step is Δt = 0.02 s. As for the viscous fluid flow model, uniform and non-uniform (1:8) meshes are adopted in the x-axis coordinates of wave propagation and relaxation zone, respectively. In the vertical coordinate, non-uniform meshes with 1:2 ratio of mesh sizes between free surface and wave flume top/bottom are employed. Table 3 shows three types of mesh schemes in the computational domain, in which the time step is determined according to the Courant−Friedrichs− Lewy (CFL) condition, automatically. Fig. 3 shows numerical results obtained by the former two numerical models with three mesh schemes. It can be observed that those results obtained with all three meshes are very close to each other for nonlinear potential flow model; while Mesh 2 can generate the convergent results in the viscous fluid flow results. In addition, those convergent results in the viscous fluid flow model also indicate that the numerical dissipation has been included in the error of numerical model, which is dependent of the present mesh scheme. Therefore, Mesh 1 and Mesh 2 are considered as baselines of the following nonlinear potential flow and viscous fluid flow simulations, respectively. η η Fig. 4 shows the variation of maximal wave crest (x P , t P ) with input wave amplitudes A i by paddle motion, where experimental results by Baldock et al. (1996), two present numerical simulating results, and the linear potential flow solution are also included for the purpose of comparison. It can be seen from the figure that the linear potential flow model under-estimates the maximal wave crest (x P , t P ), especially for the large input wave amplitude A i . The increase of input wave amplitude can induce the increase of the free surface nonlinearity on the focused wave amplitude, which cannot be simulated by the potential flow model correctly. The present nonlinear potential flow and viscous fluid flow results are in agreement with experimental measurements, implying that two present numerical models can be qualified for predicting the maximal wave crests. Furthermore, the variation of focused time and focused position of Cases B and D with input wave amplitudes A i are illustrated in Fig. 5. Based on linear potential flow theory, the occurrence of the focusing event should be at the time instant t P = 0.00 s and the position x P = 0.0 m. With the increase of the input wave amplitude A i , we can observe the increasingly delayed focused time in Fig. 5a. The larger downstream shift of the focused position can be also observed clearly in Fig. 5b for the larger input wave amplitude. The influence of the wave amplitude on the dispersion relation of water   waves can enlarge the wave length and phase velocity, which is the reason of the delayed focused time and the downstream shift of focused position. Although numerical models are somewhat under-predicted the focused time and  focused position, the tendency of the delay and shift can be simulated correctly. In a word, both of the present nonlinear potential flow and viscous fluid flow models are qualified for predicting the focused wave behavior in this study. Finally, all the phenomena in Figs. 4 and 5 indicate the stronger interactions between wave components due to the increase of wave group amplitude.
The influence of wave components on the free surface elevation of focused wave is investigated, in which the input wave amplitude A i = 0.055 m is used. Fig. 6 shows the free surface elevations at positions of x = 0.00, 0.35 and 0.55 m with different numbers of wave components N by using the nonlinear potential flow model. It can be observed from this figure that there is only a small discrep-ancy among different wave components N results. As for viscous fluid flow results in Fig. 7, all the wave evolutions by different wave components N are almost identical with each other. On the whole, the number of wave components N has insignificant effect on focused amplitudes. In the present work, wave components N = 161 are to be adopted for the following numerical simulations. η Fig. 8 shows the frequency band effect on focused crest value (x P , t P ) with the input wave amplitude A i . The increased focused crest value with the increase of input wave amplitude A i can be observed in this figure, implying that the influence of the free surface nonlinearity approaches to more significant for larger input wave amplitude A i . When the input wave amplitude A i is small enough (A i < 0.025 m  JIANG Sheng-chao et al. China Ocean Eng., 2020, Vol. 34, No. 2, P. 279-288 η η in present work), focused crest values of all the cases are nearly the same. The discrepancy between different frequency bands can only be observed for the larger input wave amplitude A i . Numerical simulations show us that the larger focused crest value (x P , t P ) can be obtained following the sequence of Cases A-D, which is the narrower frequency band. However, the focused crest value (x P , t P ) of Case E is smaller than that of Cases B-D, implying the complex phenomenon of frequency band effect. η η η The variations of focused crest value (x P , t P ) with the input wave amplitude A i among the linear potential flow, the nonlinear potential flow and the viscous fluid flow results are compared in Fig. 9. It can be observed that the focused crest value (x P , t P ) obtained by the viscous fluid flow model is smaller than that given by the nonlinear potential flow model. It indicates that the fluid viscosity can reduce the amplitude of focused wave, which can attribute to the influence of energy dissipation. Compared with the results by the linear potential solution, both the results obtained by the nonlinear potential flow and the viscous fluid flow models can predict the larger focused crest value (x P , t P ). The free surface nonlinearity is the most important factor to the fo-cused wave motion during the process of wave-wave interactions in focused wave group. It can be seen from these that the narrow-banded spectrum can generate more largeamplitude waves around the focused time and focused position. For example, only an isolate focused wave appears in Fig. 10a for Case A, while there are four large-amplitude waves around the focused wave for Case E in Fig. 10e. The stronger wave group behavior can be induced by the narrow-banded spectrum, which can be simulated by all the three types of numerical models in the present study. At the position of x = 0.00, the maximal focused crest (t)/A i at the time instant of t = 0.00 can be obtained by the linear potential flow model. However, the nonlinear potential flow and the viscous fluid flow results cannot approach to the maximal focused crest value. The development of focused wave by the nonlinear potential flow and the viscous fluid flow models are slower than that by the linear potential flow model,   The evident discrepancy of the maximal focused crest (t)/A i can be observed between the linear and nonlinear potential flow results, indicating the significant effect of free surface nonlinearity. Furthermore, the comparisons between the nonlinear potential flow and viscous fluid flow results reveal that the fluid viscosity can slightly reduce the maximal focused wave crest (t)/A i . The corresponding time instant of wave crest by the viscous fluid flow model is a little smaller than that of the nonlinear potential flow model. Generally speaking, the free surface nonlinearity has significant effect on the behavior of focused wave; while the influence of fluid viscosity is quite small in the present study.

Conclusion
A series of numerical simulations on focused wave propagation by using three types of numerical models, which are the linear potential flow, nonlinear potential flow and viscous fluid flow models, are carried out in the present study. The calculations provide comprehensive numerical results of the focused wave crest, focused time and focused position, by which the influence of the free surface nonlinearity and fluid viscosity on the focused wave behavior has been investigated. Numerical simulations reveal that the free surface nonlinearity has significant effect on the characteristics of the focused wave, including the increased focused wave crest, delayed focused time and downstream shift of the focused position. It can be obtained by the evident discrepancy between the linear/nonlinear potential flow and the viscous fluid flow results. However, only a little discrepancy of results between the nonlinear potential flow and viscous fluid flow models can be observed, which implies the insignificant influence of fluid viscosity on related phenomena of focused wave motion. According to all of the numerical results, we recommend the nonlinear potential flow model to be used for simulating the focused wave problem. Finally, we need to reaffirm that the conclusion in the present study is only valid for the non-breaking focused wave propagating problem.