Large scale CFD modeling of wave propagation into Mehamn harbor

Large Scale CFD Modeling of Wave Propagation into Mehamn Harbor VII International Conference on Computational Methods in Marine Engineering MARINE 2017 M. Visonneau, P. Queutey and D. Le Touz´e (Eds) LARGE SCALE CFD MODELING OF WAVE PROPAGATION INTO MEHAMN HARBOR Weizhi Wang∗, Hans Bihs†, Arun Kamath†,and Øivind Asgeir Arntsen† ∗Department of Civil and Environmental Engineering Norwegian University of Science and Technology Campus Gløshaugen, 7491 Trondheim, Norway e-mail: [email protected], web page: www.reef3d.com †Department of Civil and Environmental Engineering Norwegian University of Science and Technology Campus Gløshaugen, 7491 Trondheim, Norway e-mail: [email protected], web page: www.reef3d.com Key words: Large Scale, Irregular Bottom, CFD, REEF3D Abstract. Ocean wave propagations into harbours are large scale phenomena with complex wave transformations. Computational fluid dynamics (CFD) has the advantage of capturing most physics with few assumptions. It has also been successfully applied in marine engineering and coastal engineering. Altogether, CFD is considered to be an ideal tool to analyse the wave propagation at harbours. The most prominent limitation of CFD application is the high requirement of computational resources. However, with increasing computational resources, CFD is becoming an attractive alternative. One of the challenges in large scale CFD simulation is to generate the realistic waves. An extended flat-bottom wave generation zone tends to interrupt the continuity of the sub-sea terrain and introduces unrealistic wave transformations. To locate the wave generation zone far from the topography is one remedy, but it demands more computational resources. Therefore, generating waves over an irregular bottom is of particular interest in CFD simulations. In this paper, the three dimensional large scale numerical simulation of wave propagation into Mehamn harbour is performed with waves generated over an irregular bottom using the open source CFD model REEF3D. The relaxation method is used for the wave generation and absorption. A modified wave generation method considering the local water depths and wave numbers is used in this paper. A study case is performed to demonstrate the effect of the irregular bottom wave generation. Then, the large scale wave propagation into Mehamn harbour is simulated with two different waves generated over the real topography. REEF3D simulates the Mehamn harbour wave propagation by solving Navier-Stokes equations on a staggered grid with the finite difference method. The level-set method is applied to capture the free surface. A fifth-order WENO scheme is applied to the convection terms and a thirdorder TVD scheme is applied on the transient terms. The topography of the harbour is modelled using the local inverse distance interpolation method. The Mehamn harbour simulations show good wave transformation results and indicate a successful application of wave generation over irregular bottoms. 1 873 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen 1 INTRODUCTION Mehamn harbour is located in Nordkinnhalvya in the Finnmark municipality, Norway. The harbour is connected to the open sea to the north and relatively well protected from waves coming from the east and the west. The geography of the harbour is shown in Fig. 1. Commercial cargo vessels and large fishing vessels dock in the outer part of the harbour and smaller vessels are moored in the inner part of the harbour. In order to have calmer wave conditions for the mooring of the vessels, Finnmark municipality and the Norwegian Coastal Administration decide to arrange two breakwaters along the harbour. Therefore, an accurate modelling for the wave propagation into the harbour is demanded to optimise the effect the breakwaters. The Norwegian meteorology institute conducted the measurements of the wave data in the offshore region outside Mehamn harbour during 1955 to 2006 [1]. Based on the measurements, SINTEF[1] performed a model test on the wave propagations into the harbour with a scale factor of 1:80. Several numerical wave models have been developed to simulate the coastal waves, such as Simulating Waves Near shore (SWAN)[3] developed by Delft University of Technology. As a third generation wave model, SWAN has been successfully used in predicting significant wave heights and statistical wave properties. However, the phase-averaged results from SWAN are not able to show the details of wave transformations. A phase-resolved model is needed to show the wave propagation patterns and transformation phenomena, such as diffractions and refractions. Computational Fluid Dynamics (CFD) is able to capture the details and complexities in a fluid domain in an accurate way with less assumptions. Therefore, CFD has been used in a wide range of hydrodynamic applications and has shown great potential for further development. One of the challenges of applying CFD on coastal waves is the high demand of the computational resources. However, with the increasing capacity of the computational resources, CFD can be applied on large scale simulations. REEF3D [4] is a CFD model with specialisation on wave hydrodynamics developed at Norwegian University of Science and Technology (NTNU). The model has been applied to a large range of topics within wave hydrodynamics, such as floating body dynamics in waves[5], breaking wavs[6] and breaking wave-structure interactions[7]. The simulations give high-resolution phase-resolved results which present remarkable details of different wave phenomena. In this paper, REEF3D is used to simulate the large scale wave propagation into Mehman harbour. In order to have phase-resolved results, a good wave generation must be ensured. With the intermediate wave condition in Mehamn harbour, the wave properties are influenced by the water depth variations. Extending the wave tank to arrange a flat bottom wave generation zone might cause unnecessary and unphysical wave transformations because of the sudden geometrical change from the flat bottom to the complex topography. Arranging a flat bottom wave generation zone several wave lengths away from the topography change is a remedy but requires more computational resources. Therefore, a method of generating waves over an irregular bottom with intermediate water condition is needed. With intermediate water condition, the analytical wave velocities and wave lengths are functions of the local wave depth (d) and the wave number (k). 2 874 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen Figure 1: Location and geography of Mehamn harbour and surrounding area Thus a modified wave generation method with the consideration of d and k is used in this paper. One study case is performed to demonstrate the effect of the irregular bottom generation. Then the 3D large scale CFD simulations of wave propagation into Mehamn harbour are performed with the waves generated over the real topography. The simulations are performed with two different wave conditions and the wave transformation phenomena are compared and analysed. 2 NUMERICAL MODEL 2.1 Governing Equations Water wave hydrodynamics are well described by the three dimensional incompressible NavierStokes equations, as shown in Eqn. (1) and Eqn. (2). For large scale wave modelling, the turbulence terms are neglected from the equations. ∂ui ∂xi = 0 (1) ∂ui ∂t +uj ∂ui ∂xj =−1 ρ ∂p ∂xi +∂ ∂xjν∂ui ∂xj +∂uj ∂xi+gi(2) where uis the velocity, ρis the fluid density, pis the pressure, νis the kinematic viscosity and gthe acceleration due to gravity. REEF3D solves the Navier-Stokes equations on a staggered grid with the finite difference method and high-order discretisation schemes. The conservative fifth-order weighted essentially non-oscillatory (WENO) [8] scheme is applied on convective terms, resulting in better robustness and smoother flux. The third-order Total-Variation-Diminishing (TVD) Runge-Kutta scheme 3 875 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen [9] is used on time-dependent terms in both Navier-Stokes equations and level-set equations. The Corant-Friedrichs-Lewy (CFL) criterion is applied with adaptive time steps. HYPRE library [11] is used to solve the pressure from the Poisson equation. Furthermore, REEF3D chooses BiCGStab [12] as the iterative solver and the geometric multi-grid PFMG[13] as the pre-conditioner. The combined solver solves the poisson equation very well at large scale simulations. The free surface is modelled with the level-set method[14]. The method applies a signed distance function denoted as the level-set function φ(x , t ) to capture the change of phases at the free surface. As shown in Eqn. (3), the function always equals to zero at the interface of the two phases and shows different signs in each phase. As a result, REEF3D is able to capture complex free surface accurately. φ(x , t )     >0i f x ∈phase 1 =0i f x ∈Γ <0i f x ∈phase 2 (3) The bathymetry of the harbour is described by a large set of scattered points. The scattered points are interpolated using the inverse distance weighting method [17] to form the continuous sub-sea topography. The improved local inverse distance weighting method [18] is adopted in REEF3D considering the cost efficiency in terms of time consumption. Instead of a global interpolation regardless of the distance from the points (xk,yk) to (x, y), the improved method uses the weight function that only considers the contribution of the points within a certain radii, as shown in Eqn. 4. By assigning the radii as a function of k and using a search technique to find the nearest points near (xk,yk), the method is much faster with due accuracy. The irregular bathymetry and coastal line of the harbour are interpolated into a separate level set function. The function is used to distinguish solid from fluid. This approach has the advantage of being very flexible when complex geometries are encountered. Wk(x, y)=(Rw−dk)+ Rwdk2 (Rw−dk)+=Rw−dkif dk<R w 0if dk≥Rw (4) where Rwis the radii of influence and dkdenotes the Euclidean distance between (x, y) and (xk,yk). The relaxation method [15] is applied to generate and absorb the waves. The velocity, pressure and the free surface are ramped up to the analytical values from the wave theories in the wave generation zone before the waves enter the computational domain. In the numerical beach, a reverse process is carried out so that the velocities are reduced to zero, the pressure is damped to the hydrostatic pressure and the free surface is relaxed to the still water level. The wave 4 876 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen generation and absorption process is described in Eqn. (5). The Γ(x) term in Eqn. (5) is called the relaxation factor, which is described in Eqn. (6) by Jacobsen[16]. Φ(x)relaxed = Γ(x)Φanalytical + (1 −Γ(x))Φcomputational (5) where Φ stands for different parameters in the fluid, including horizontal velocity u, vertical velocity w, pressure pand level-set function φ. Γ(x)=1−e(x3.5)−1 e−1for x∈[0; 1] (6) where xis scaled to the length of the relaxation zone. With the intermediate water condition, wave properties change with the water depth d. As an example, the analytical values of the velocity components of the linear wave theory are given in Eqn. (7). As described in the equations, the velocities are functions of local water depth d and local wave number k. Therefore, in the wave generation zone, dand kare not constant numbers but matrices d(i, j) and k(i, j) associated with the x and y coordinates. As a result, the wave properties are ramped up to different analytical values within the wave generation zone. u(x, z, t)=πH T cosh [k(z+d)] sinh (kd)cosθ w(x, z, t)=πH T sinh [k(z+d)] sinh (kd)sinθ (7) where His the wave height, Lis the wavelength, Tis the wave period, ωis the angular wave frequency, dis the water depth and zis the vertical coordinate measured from the still water level z= 0. kis derived from dispersion relation as a function of ωand d. 3 STUDY CASE OF IRREGULAR BOTTOM WAVE GENERATION A study case is performed to investigate the effect of the modified irregular bottom wave generation method. Two wedges along the wave propagation direction (x-axis) with different initial heights are arranged at the bottom of the wave tank to mimic the depth variation in both x direction and y direction. A two-wave-length numerical beach is located at the end of the tank to absorb the waves. Two wave velocity probes are arranged in the middle of the wave generation zone and in the middle of each wedge in each case. Two wave height gauges are arranged to be 50 m away from the beginning of the bottom topography and at the centre of each wedge in the transverse direction. A linear wave is selected to conduct the study. The chosen wave has a wave height of 1.5 m and a wave period of 15 s. The design of the bottom geometry, the locations of the velocity probes and the wave height gauges and the arrangements of the wave generation zones are all shown in Fig. 2. The velocity profiles at the probes in the wave generation zone and the wave heights at the wave gauges are compared between the irregular bottom wave generation case (IBG) and the extended flat bottom wave generation case (FBG). The mesh convergence study is conducted with the FBG case. The average wave heights at gauge 1 between 220 s and 250 s with different mesh sizes are compared in Fig. 3. The results from all mesh sizes share very similar values with only an error of 0.5 % between 5 877 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen 0 10 Y20 80 10 100 P1 (G1) Z P2 (G2) 20 d D (a) (350,10) X Z (350,20) (0,80) (0,100) (350, 0) (2100,0) (0, 0) (2800, 0) B EF GiA’ Pi’ Pi A(B’) d D (b) Figure 2: Numerical tank set-up for the study case: (a) the configurations and arrangements of wave velocity probes Piand wave gauges Giin the cross section (b) the configurations and the arrangements of wave velocity probes Piand wave height gauges Giin the longitudinal direction, i=1,2. EF is the numerical beach, AB is the wave generation zone for the IBG case, A’B’ is the wave generation zone for the FBG case. D is the bottom geometry thickness, d is the wave depth dx=1 m and dx= 0.8 m cases. Therefore, the mesh size in this study is chosen to be 1 m. The velocity profile at the wave probe P1 and the wave height record at the wave gauge G1 from the IBG case are compared with the counterparts from the FBG case in Fig. 4. The local water depth in the IBG case changes the velocity profile at P1 and as a result the wave height and phase are also slightly different in the computational zone at G1. 4 SIMULATION OF MEHAMN HARBOR The measurements described in the SINTEF report [1] show that the highest waves come from 345◦to 45◦. In this paper, it is considered to be the most dangerous situation when the highest waves come directly into the entrance of the harbour, i.e. when the waves come from 345◦. Therefore, the original bathymetry is rotated 15◦to the east to have the right wave direction. It is stated in the SINTEF report [1] that the 100 year return period extreme wave has Hs100 of 4.5m and Tp100 of 15s. NORSOK [2] recommends that the most extreme wave has the wave height of 1.9Hs100. Therefore, in order to simulate the most dangerous wave, the wave height is chosen to be 9 m. With the same wave height, two regular waves are studied with different wave periods. One wave (W1) has the wave period calculated from Eqn. (8) [2], which corresponds to the extreme design wave. Another wave (W2) has the same wave period as Tp100, which describes the waves contributing the most wave energy. The 5th-order Stokes wave theory is used for W1 while the 2nd-order Stokes wave theory is applied to describe W2. 6 878 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen H (m) 1.489 1.490 1.491 1.492 dx (m) 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 Figure 3: Mesh convergence study of the irregular bottom wave generation study case P1 (IBG) P1' (FBG) z (m) 0 10 20 30 40 50 60 70 80 90 |u| (m/s) 0 0.1 0.2 0.3 (a) G1 (IBG) G1' (FBG) z (m) 79.2 79.4 79.6 79.8 80.0 80.2 80.4 80.6 80.8 t (s) 180 190 200 210 220 230 240 (b) Figure 4: (a)Comparison of velocity profiles at probes P1 and P1’ at t=250 s, (b)Comparison of the wave heights at the wav gauges G1 and G1’ The smooth bottom topography obtained from the local inverse distance interpolation method based on a large set of scattered data is shown in Fig. 5. The mesh size is 4m, which gives about 35 and 80 cells per wave length for each wave. The black box in Fig. 5 shows the wave generation zone, which is one wave length long for both wave conditions. There is no numerical beach implemented in the wave thank because the harbour has a closed domain and the region at the end of the tank outside the harbour has shallow water condition, the wave transformation phenomena is better illustrated by including the shallow water part instead of damping the waves out deliberately. The left and right side of the tank have symmetric boundary condition to reduce the reflection from the boundaries. 6.5H100 ≤T≤11H100 (8) The simulation results for W1 and W2 are shown in Fig. 6 and Fig. 7 respectively. The waves start to refract right after they enter the computational zone because of the deep trench at the beginning of the domain. From both velocity graph and surface elevation graph, one is able to observe that the reflected waves from the peninsulas at the sides of the harbour radiate to the 7 879 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen Figure 5: Configuration and bottom topography of the numerical tank of Mehamn simulation, the black box shows the wave generation zone surroundings. The waves pass through the narrow channel at the entrance of the harbour and curve towards the far end of the harbour. The mean wave level rises up in the shallow area of the harbour and waves dissipate energy resulting in smaller wave amplitude inside the harbour. The various wave transformation phenomena are well visualised and represented in the 3D large scale CFD simulation. The details of the wave patterns facilitate the designs of the breakwaters. 5 CONCLUSIONS In this paper, the local inverse distance interpolation method is successfully utilised to generate the complicated bottom topography at Mehamn harbour. The irregular bottom wave generation method is also successfully applied to generate regular waves for a large scale domain. The CFD simulations at Mehamn harbour are performed with two different wave conditions and both show high-resolution phase-resolved results with good visualisation. The wave transformations into the harbour are well represented in the simulation results. CFD method is proved to be able to simulate large scale costal waves and REEF3D is shown to have satisfactory capability of performing such simulations. Further studies of irregular wave CFD simulation at Mehamn harbour and the validation against experimental data will also be performed in the future research. 8 880 Weizhi Wang, Hans Bihs, Arun Kamath,and Øivind Asgeir Arntsen (a) (b) Figure 6: Free surface velocity (a) and free surface level (b) in Mehamn simulation with H=9 m and T=9.5 s (W1) at simulation time of 200 s ACKNOWLEDGMENT This study has been carried out under the E39 fjord crossing project (No. 304624) and the authors are grateful to the grants provided by the Norwegian Public Roads Administration. This study was supported in part with computational resources at the Norwegian University of Science and Technology (NTNU) provided by NOTUR, http://www.notur.no. REFERENCES [1] Vold, S. and Lothe, A.E. Mehamn Modellforsøk, SINTEF, Klæbuveien 153, Trondheim (2009), project number: 3C0337, report number: SBF IN F09203. [2] NORSOK, Actions and action effect, Standard N-003, Norwegian Oil Industry Association and the Federation of Norwegian Industry (2007). [3] Thomas, T. J., and Dwarakish, G. Numerical wave modelling a review, Aquatic Procedia, 4(2015), pp. 443 448. [4] Bihs, H., Kamath, A., Chella, M.A., Aggarwal, A., and Øivind A.A., A new level set numerical wave tank with improved density interpolation for complex wave hydrodynamics, Computers Fluids, (2016) 140, pp. 191 – 208. [5] Bihs, H. and Kamath, A., A combined level set/ghost cell immersed boundary representation for floating body simulations, International Journal for Numerical Methods in Fluids, (2016),fld.4333. 9 881