Yuan Zhuang, De-cheng Wan

Computational Marine Hydrodynamics Lab (CMHL), School of Naval Architecture, Ocean and Civil Engineering,Shanghai Jiao Tong University, Shanghai 200240, China

Abstract: This paper presents a developed new coupled method which combined our in-house CFD solver naoe-FOAM-SJTU and naoe-FOAM-os with a potential theory High Order Spectral method (HOS). A parametric study of nonlinear wave propagation in computational fluid dynamics (CFD) zone is considered. Mesh convergence, time step convergence, time discretization scheme and length of relaxation zone are all carried out. Those parametric studies verify the steady of this new combined method and give better choice for wave propagation. The dissipation in propagation of nonlinear regular wave can be lower than 3% in static mesh, and less than 2% in overset grid mesh. Meanwhile, a LNG FPSO is put into the viscous wave tank to study the suitable size of CFD zone. To achieve a better solution with least calculating resources and best numerical results, the length of CFD zone is discussed. These parametric studies can give reference upon employment of the potential-viscous coupled method and validation of the coupled method.

Key words: High Order Spectral method (HOS) parametric study, naoe-FOAM-SJTU solver, potential-viscous coupled method

Introduction

In recent years, more and more studies focused on freak waves interacting with structures. The needs for studying freak waves induce the improvement on numerical methods. The traditional computational fluid dynamics (CFD) method has a great advantage in dealing with wave structure interactions, especially in nonlinear phenomenon. However, it will take large amount of computational cost when simulating large domain cases or long-time duration cases. To solve this kind of problem, many researchers start to combine potential solver with CFD method. The viscous effect is not essential in far-field, thus the application of the potential solver is reasonable. The solution near the structure using CFD method can consider not only viscous effect but also the large amplitude motion and nonlinear phenomenon. Many potential-viscous studies have been done in recent years. The main idea of potential-viscous flow coupling method is to solve near field around body with viscous flow method while far-field flow using potential theory. Campana et al.[1]applied Navierstokes equation with potential flow code to simulate a flow past a ship hull. Those kinds of coupling can be described as domain decomposition. Colicchio et al.[2]coupled BEM with Navier-Stokes equation using level-set method to capture the free surface. They tested the coupled method through dam break and impact on the body. Sitanggang and Lynett[3]also studied the shoal wave breaking. They combined Boussinesq model with Reynolds-averaged Navier-Stokes (RANS) model to do the simulation. In general,some potential-viscous coupled methods were developed to study freak waves interacting with structures. Some researchers applied quasi arbitrary Lagrangian-Eulerian finite element method (QALEEM) coupled with OpenFoam[4], Lagrangian wave model with olaFlow[5], or SWENSE method[6]to simulate focusing waves acting on the structures.Iafrati and Campana[7]used the potential-viscous domain-decomposed method to simulate wavebreaking phenomena.

Since the high efficiency and accuracy model High-Order Spectral method (HOS)[8-9]has been widely used and developed, the coupling methods with CFD for wave structure interaction become popular[10-13]. Wave signals can be induced by HOS in a larger computational domain and long-time duration,and then they can be transformed to CFD domain.With the help of open source HOS model HOS-ocean[14], HOS-NWT[15]and wrapped program Grid2Grid[16], Zhuang et al.[17]have been developed combined method with waves2Foam[18]and validated the accuracy of wave evolution in empty field.

This paper presents the coupled method of HOS and CFD in our in-house solver naoe-FOAM-SJTU[19].The in-house solver naoe-FOAM-SJTU are well developed with several functions which can solve lots of problems of hydrodynamics, such as ship motion in waves[20], open-water propeller[21], etc.. Meanwhile,the coupled method of HOS and CFD solver with overset module naoe-FOAM-os is considered. A parametric study is conducted with the naoe-FOAMSJTU and naoe-FOAM-os, and the results are given.Besides, the interface of HOS and CFD is based on a relaxation zone, the wave reflection is also discussed.With a moving FPSO in viscous solver zone and regular wave generated by HOS, the validation and discussion of ship motion in different computational size is given.

1. HOS-CFD coupling method

This section introduces HOS method and CFD method as well as the combination between these two methods. As we all know, the time steps and mesh size in potential theory are not match with that in CFD method. Therefore, the combination must deal with the mismatch and build an interface between these two methods to achieve information transform.

1.1 Viscous solver

naoe-FOAM-SJTU[19,22]is an in-house solver based on the open source software OpenFOAM with six DOF module, mooring line module, overset grids and some turbulence models. The governing equations are:

where U is velocity field, Ugis velocity of grid nodes, pd=p-ρg· x is dynamic pressure. The volume of fluid (VOF) method with bounded techniques is applied.

The in-house solver naoe-FOAM-SJTU and naoe-FOAM-os both applied VOF to capture the free surface. To avoid numerical diffusion, the bounded compression techniques are applied. The PIMPLE scheme is implemented to decouple the pressure and velocity, based on pressure-implicit split operator(PISO) algorithm[23].

1.2 HOS method

The HOS method scheme is a pseudo-spectral method. It based on the partial difference equation on dynamic and kinematic free surface boundary condition shown in below[8]:

wheresφ is the surface velocity potential. After we write φ in a perturbation series and further expand each order of φ evaluated on free surface in a Taylor series, then we have

with the initial conditions known for the simulation,Eq. (5) can be solved for unknow φ.

1.3 Relaxation zone

As we have given time and due to Eqs. (3), (4)we can get the unknown φ. The results of HOS can be combined into CFD zone due to the relaxation zone in waves2Foam[18].

The relaxation weightRα can be chosen as three types, and the default choice is using exponential weight. The exponential weight has the following form[17]

The exponent p is set to 3.5 as default. σ is a local coordinate with relaxation zone, the value is from 0 to 1 and depends on the shape of relaxation zone.

1.4 Coupling method

Figure 1 shows the calculation process of the combined solver of naoe-FOAM-SJTU and naoe-OAM-os, respectively. First of all, the calculation is conducted in HOS method, in this paper, we applies the open source software HOS-Ocean[14]and HOS-NWT[15]. The results of HOS model will use as wave signal input into CFD zone with the help of Grid2Grid[16]and relaxation zone in waves2Foam.After that, the information of HOS wave field is transferred into CFD zone through an interface we build in CFD zone. After receiving the wave input signal, the wave propagation and evolution is in progress with body motion and force update. Those process can be adopted as numerical wave tank, the wave force will act on the body, and the body will be excited by wave force and induce motions. The CFD solver with overset module naoe-FOAM-os has an extra process on iteration of mesh information transformation. With the help of sugar++[24], the hull grids and background grids can keep the topology from each other.

2. Parametric study

A periodic wave simulation was carried out to test the steady of current solver and to chase for a better parameter choice for wave elevation. The wave condition is chosen as fully nonlinear waves (Stream Function wave), and the wave propagates with periodic boundary condition (the fluid from outlet will input into inlet). Therefore, the decay of the wave elevation is due to the numerical scheme or other parameters, not the reason of boundary condition. The wave condition is shown in Table 1. The cases are 2-D with length of 8.0820 m (10)λand height of 0.7438 m (-0.6000 m＜z＜0.1438 m).

Table 1 Wave condition of stream function wave

2.1 Parametric study in static mesh

Firstly, we take the parametric study in static mesh. Before we make simulation, the dissipation in interacted mesh scheme must be discussed. If wave spreads in wave tank in a long-time simulation, one needs to know whether wave propagation will keep still. To figure out the problem, we define an empty wave tank without body to see how it affected the wave propagation.

Table 2 gives six cases when simulating the periodic wave simulation. Three different time integrations, mesh generations and time steps are shown in tables. It needs to be mentioned that when comes to time integration scheme, the Crank Nicolson scheme with weighting factor is adopted. When the weighting factorC0=0, the scheme is implicit Euler integration scheme; when the weighting factorC0=1,the scheme is explicit Crank Nicolson scheme. The three different time integration scheme cases have the same mesh generation with the cells are 100 in every wave length (λ/Δx=100), and 20 cells in wave height (H/ Δz= 20), with time step oft/T=400.Wave dissipates quicker in low order time integration scheme, shown in Fig. 2(a) . Therefore,when the time integration scheme chooses to be the Crank Nicolson scheme and the factor adopted withC0=0.95, the wave dissipation can be dismissed most.

Fig. 1 The calculation process of the combined solver in naoe-FOAM-SJTU and naoe-FOAM-os

Table 2 Test cases in static mesh

The three different time step cases have the same mesh generation with the cells are 100 in every wave length (λ/Δx=100), and 10 cells in wave height(H/ Δz= 10). The time integration scheme is chosento be Crank Nicolson with weighting factorC0=0.95. The time step is 200, 400, 800 per period(t/T=200, 400, 800). The time step cases have little difference in first 15 periods, as shown in Fig. 2(b).The small time step decreases the wave damping and according to the other two time step cases, the decreasing is proportional to the time step.

Fig. 2 (Color online) First harmonic wave elevation in different time integration scheme, time step and mesh generation cases

The last three cases are with different mesh generations and time steps. To keep Courant number the same, the mesh cells increase with time step decreases simultaneously. The time integration scheme is chosen to be Crank Nicolson with weighting factorC0=0.95. Figure 2(c) provides the first harmonic wave elevation in 40 periods, it can be seen that comparing to case mesh-10, case mesh-20 and mesh-40 doesn’t have much difference in wave elevation. The case mesh-40 shows little dissipations in wave generation. In three situations, the best choice in time integration case is Crank Nicolson 0.95 with 8% wave dissipation, in time step case ist/T=800 with 17% wave dissipation, and in mesh generation case is mesh-40 with 3% wave dissipation.

Figures 3, 4 and 5 show the spatial distribution of phase difference of these six cases. The subtitle (a) to(c) illustrates the timet=10T, 20T, 30Tand 40T. It can be seen that along with the simulation, the wave propagation has not only wave dissipation, but also phase difference. The Euler scheme has the most serious dissipation, Crank NicolsonC0=0.9 has a little phase difference. The case of Crank NicolsonC0=0.95 is acceptable both in wave dissipation and phase difference.

Fig. 3 (Color online) Spatial distribution of phase difference in different time integration scheme cases

Fig. 4 (Color online) Spatial distribution of phase difference in different time step case

Fig. 5 (Color online) Spatial distribution of phase difference in different mesh generation and time step cases

The phase difference in time step cases is not obvious between each case, the small time step does not reduce the phase difference. The phase difference and wave dissipation can be acceptable in the simulation time of 20 period.

In the timet=40T, the mesh-10 case has obvious phase difference while mesh-20 and mesh-40 shows little phase difference. Although the result of mesh-40 case is closer to the analytic data, the computational time costs much on it. The parameter of wave flow chooses mesh generation with 200 cells in every wave length (λ/ Δ x =200), and 40 cells in wave height (H / Δ z = 40). The time integration scheme is chosen to be Crank Nicolson with weighting factor C0=0.95 and the time step is 800 per period (t/ T =800). This kind of parameter can give the wave propagation better results. However,this parameter choice only suits for 2-D cases and costs large computational resources. As the case of mesh-20 (mesh generation with the cells are 100 in every wave length (λ/Δx =100), and 20 cells in wave height (H / Δ z = 20)), the time integration scheme is chosen to be Crank Nicolson with weighting factor C0=0.95. The time step is 400 per period (t/ T =400) shows little difference in phase difference and the wave dissipation is acceptable, the simulations in discussing the length of relaxation zone choose this kind of parameters.

Figure 6 shows the wave fields in different time integration scheme, mesh generation and time step at time t=40T. The wave dissipation can be observed in each case, and the extent of the dissipation is almost the same with wave elevation mentioned above. The wave fields in Crank-Nicolson time integration scheme with factor of 0.9 and 0.95 are similar. Among those wave fields, mesh-40 case is the best one which can barely find the wave dissipation. In time step cases, with coarse mesh generation, the wave dissipation is still apparent.

Fig. 6 (Color online) Comparison of wave fields in different time integration scheme, mesh generation and time step cases

Figure 7 shows the velocity fields in different time integration scheme, mesh generation and time step at t=40T. It can be seen that with the increase of the factor in Crank-Nicolson time integration scheme, the velocity around the wave due to the simulation decreases, which means a better steady in simulations. In different mesh generation cases, the analysis threshold value of the velocity keeps the same. Therefore, it can be observed that with finer mesh generation, the simulations are steadier. The velocity fields show little differences in different time step cases.

Fig. 7 (Color online) Comparison of velocity fields in different time integration scheme, mesh generation and time step cases

The length of relaxation zone is discussed for efficient wave generation and absorption. The relaxation zone is not only an interface with transforming the wave signal into CFD zone, but also acts as a zone of wave absorption. Therefore, choose a proper length of relaxation zone can reduce the computational resources and avoid large wave reflection. Figure 8 shows three settings of outlet relaxation zone length. The inlet relaxation zone is 0.8082 m (1)λ length and the computational domain is 8.0820 m (10)λ. The length of outlet relaxation zone is 0.4041 m (0.5)λ, 0.8082 m (1)λ and 1.2123 m (1.5)λ), respectively.

Fig. 8 (Color online) The settings of different length of outlet relaxation zone

Figure 9 illustrates the first harmonic wave elevation of three different relaxation zone length. It can be seen that in the time of first 15 periods the curve of the first harmonic of wave elevation is almost the same. In the time of last 25 periods, the outlet with 1λ (marked with outlet-1) and with 1.5λ (marked with outlet-1.5) doesn’t dissipate anymore, while the value of outlet with 0.5λ continues to decrease.

Fig. 9 (Color online) The first harmonic wave elevation of three different relaxation zone length

Figure 10 illustrates the time history of Courant number of three different relaxation zone length.Within 40 periods of simulations, the Courant number is under value of 2 in those three cases. When the simulation time reaches 40 period, the longest length of relaxation zone case shows large courant number which is larger than 2.

Fig. 10 (Color online) The time history of Courant number of three different relaxation zone length

Figure 11 shows the comparison of velocity fields among those three different relaxation zone lengths. The velocity fields show little difference in the domain faraway from relaxation zone, while the velocity accumulates near the relaxation zone.Comparing with courant number, first harmonic wave elevation and velocity fields, it can be seen that long relaxation zone is not the best choice. When the mesh size is average along the x direction, the outlet wave cannot be wiped. Therefore, in the following simulation, the choice of relaxation zone length must follow the conditions of cases and the mesh should be coarse.

Fig. 11 (Color online) Comparison of velocity in different length of relaxation zones

2.2 Parametric study in overset mesh

In the ocean engineering, the offshore structures will face with complex waves or large steep waves.Therefore, the numerical simulation needs to use overset grids to keep the topology between hull grids and background grids. In order to figure out the steady of overset grids in the new coupled method and the dissipation, the parametric study in overset grids is carried out. The size of computational domain and the choice of wave conditions as well as boundary conditions are the same as that in static mesh. In overset grids cases, there exists an empty hull grid.There is no structure in the hull grid, in order to study the numerical dissipation between hull grid and background grid. The length of the hull grid isλ, ten percent of the background length.

To make the hull grid move violently, the hull grid is set to rotate and move around the zero point periodically. Figure 12 shows the motion of the hull grid, the translation motion isx= 0.5cos(4t) and the rotation motion isω= cos4t.

Fig. 12 (Color online) The setup of simulation of overset grids case

Table 3 gives the test cases in overset grids, they are three different time integration scheme cases and three mesh generations, respectively. Contrary to the static mesh cases, the time step in overset grids changes in order to keep the simulation steady. When the order of time integration scheme is large, the time step needs decrease to make the simulation steady.Meanwhile, the fine mesh cases also need small time steps.

Table 3 Test cases of overset grids

Figure 13 illustrates the first harmonic wave elevation in overset grids cases. It can be seen that the trends of the wave elevation in overset grids are similar to that in static mesh grids. The Euler time integration scheme shows large wave dissipation,while the Crank-Nicolson time integration scheme with factor of 0.95 keeps the wave elevation well.Comparing to the static mesh cases, the time steps in overset grids cases are smaller, which give less wave dissipation in simulations. In Fig. 13(b), the error between simulation and analytical results in the finest mesh cases is about 2%.

Fig. 13 (Color online) First harmonic wave elevation in different time integration scheme, time step and mesh generation cases

Figure 14 shows the phase difference in different time integration schemes. It can be seen that the phase difference is obvious in Euler scheme, while Crank-Nicolson with 0.95 factor still keeps well phases in 40 periods.

Fig. 14 (Color online) Spatial distribution of phase difference in different time integration scheme cases

Figure 15 illustrates the phase difference in different mesh generations. The wave dissipation in coarse mesh is obvious although the order of time integration scheme is high. The medium mesh and fine mesh keep well phases in 40 periods.

Fig. 15 (Color online) Spatial distribution of phase difference in different mesh generation cases

As for overset grids, the dissipation not only happens in wave propagation in background mesh, but also between hull girds and background grids.However, the first harmonic wave elevation and phase difference of overset grids shows the similar trends with that in static mesh, which verifies that the dissipation due to motion of the overset grids is not obvious.

In order to show the numerical transformation more clearly, Figs. 16, 17 gives the wave fields and velocity fields of different mesh generations int= 40T. The wave dissipation in coarse mesh is the most obvious, with small velocity around the free surface. The accumulation of velocity around free surfacein fine mesh is small and keeps well shape of free surface. For those three types of meshes in overset grids, the wave fields and velocity fields show smooth transfer from background mesh to hull mesh, without the sudden change between these two meshes.Therefore, the motion of the hull grids will not affect the wave propagation.

Fig. 16 (Color online) Comparison of wave elevation in different mesh generations

Fig. 17 (Color online) Comparison of velocity field in different mesh generations

It needs to be mentioned that when the motion in overset grids is very large and violent, we need to guarantee the steady of the simulation. In this way, the Euler scheme has advantages over Crank-Nicolson scheme. When compared with the first harmonic wave elevation between overset grids and static grids, the Euler scheme decays slowly when the time step is small. Therefore, if the simulation is stable in 15 periods and the simulation needs to be steady, the Euler scheme can be chosen as time integration scheme.

3. Simulation of FPSO in wave

After we find a better way to set the parameters of the simulation and length of relaxation zone, the size of computational domain is considered. The test cases choose a FPSO with 3DOF and 0 filling ratio in a head wave. The detailed numerical setup of ship and tanks parameters can be found in Zhuang and Wan[25].

The incident waves are regular waves with height of 0.025 m. Three different wave frequencies are chosen to test the parameters of viscous zones. The wave frequencies are normalized as, and the values are 2, 2.5 and 3, respectively. The waves are generated by HOS-NWT. The length of CFD zones without relaxation zones are 2.8500 m (1)λ,4.2750 m (1.5)λ and 5.7000 m (2)λ, respectively.The inlet relaxation zone and outlet relaxation zone are both 2.8500 m (1)λ.

The heave motion of FPSO of three different settings and three wave conditions are illustrated in Fig. 18 (in which when normalized wave frequency equals to 2.5, the data of experiment is absent). To compare with the experimental results in Nam et al.[26],the heave motion is normalized to dimensionless data.The normalized heave motion illustrates as ξ/A,which ξ is the maximum value of heave motion, A is wave amplitude.

Fig. 18 (Color online) Comparison of normalized heave motion in different viscous zone and wave conditions

The heave motion of FPSO of three different settings and three wave conditions are illustrated in Fig. 19. The normalized pitch motion illustrates as θL / 2A, which θ is the maximum degree of pitch motion, L is the perpendicular length of the ship. It can be seen that comparing to the pure CFD method,the HOS-CFD shows a similar value both in pitch and heave motion. The values of different viscous zones show little difference, but they are not the same. This illustrates that there exists reflection in the simulations.

Fig. 19 (Color online) Comparison of normalized pitch motion in different viscous zone and wave conditions

The wave contour around the ship in different viscous zones in normalized wave frequency of 2.5 is shown in Fig. 20. There exists wave absorption behind the ship stern, but the wave reflection still happens in relaxation zone scheme. The way to reduce the wave reflection is to make the mesh in absorption zone coarse to do the numerical dissipation. Although the small viscous zone can reduce the execution time and computational resources, the scattering wave around ship cannot be developed. In Figs. 18, 19, it can be seen that the motion of medium viscous zone shows smaller value of that in other two conditions, gives a wave reflection or missing of developed scattering wave field.

Fig. 20 (Color online) The contour of free surface around the ship in different computational size

The execution time of pure CFD, HOS-CFD with 1λ, 1.5λ and 2λ computational size in normalized wave frequency 2.5 case are 914 631 s, 608 082 s,674 987 s and 854 198 s, respectively. It can be seen that the execution time of HOS-CFD is a little smaller than that of CFD method. This may due to the absence in wave propagation equation solving during the simulations.

4. Conclusions

This paper introduces a potential-viscous combined solver HOS-CFD and gives a parametric study based on this solver. The combination is based on our in-house solver naoeFOAM-SJTU and naoeFOAM-os with overset module. The parametric studies include mesh convergence study, time step convergence study,time discretization study and the discussion on relaxation zone and size of computational domain.The parametric studies of mesh, time step and time discretization scheme not only give a better choice on those parameters during the simulation, but also show the steady of the combined solver. The mesh generation, time discretization scheme as well as time step studies are also carried out in naoeFOAM-os with empty moving hull grids. These parametric studies in overset grids illustrates that the information transform between two kinds of grids (hull grids and background grids) affects little on wave propagations.

The reflection of the relaxation zone is also discussed. This paper gives a better choice on the length of relaxation zone due to the default number and settings in waves2Foam. Meanwhile, a moving FPSO in regular wave is validated and discussed in three different computational sizes and three different wave frequencies. Although the wave reflection in outlet zone is the same for these configurations, the small viscous zone still exists reflection waves. Meanwhile, the non-developed wave scattering around the ship also affects the ship motion. For computational cost, the HOS-CFD solver cost little computational resources to simulate the wave-structure interaction,and its steady and accuracy is verified.

Acknowledgements

This work was supported by the Chang Jiang Scholars Program (Grant No. T2014099), the Innovative Special Project of Numerical Tank of Ministry of Industry and Information Technology of China (Grant No. 2016-23/09).