A Case Study of Nearshore Wave Transformation Processes along the Coast of Mexico near the Laguna Verde Nuclear Power Plant Using a Fast Simulation Method

Abstract

A numerical model based on the mild-slope equation of water wave propagation over complicated bathymetry, taking into account the combined effects of refraction, diffraction, and reflection due to breakwater, is presented. The numerical method was developed using a split proposed version of the mild-slope equation and solved by an implicit method in a finite volume grid; this technique easily allows model the wave effects caused by the breakwater building in coastal waters, where industrial and other economic activities take place. Controlled case studies have been made and the results match very well with the reference solution. The capability and utility of the model for real coastal areas are illustrated by application to the breakwater of the Laguna Verde Nuclear Power Plant (LVNPP).

1. Introduction

Recently a great deal of public concern and interest are being directed towards the coastal waters, where industrial activities take place. As a result engineers are faced with and increasing demand for solutions to coastal problems. Among these is the understanding of the water wave transformation over the complicated bathymetry, from deep water to shallow sloping beaches, which dominates various aspects of coastal hydrodinamic processes, especially the generation of alongshore currents as well as sediment transport.

The mild-slope equation- (MSE-) based wave models emerged in late 1960s and have since reached total maturity by early 1990. A number of free and commercially MSE models are available today like CGWAVE and CREDIZ models, but in the field it is not possible to use them because the engineers do not have the enough time to create the meshes and wait for the results in field. The key of the model presented is the previous calculus of some important nearshore wave processes like wave current and wave-wave interactions made with a research hydrodinamic model early developed [1].

The mild-slope equation which describes the propagation of periodic, small amplitude surface gravity waves over a seabed of mild slope was first established by Berkhoff [2]. This original equation takes into account the combined effects of refraction, diffraction and reflection, while the influence of wave breaking, bottom friction, currents and wind is neglected; for this reason we utilize a hydrodinamic model [1] to calculate these effects.

The solution of the mild-slope equation was obtained by using the finite-element method [3, 4] and finite-difference schemes [5, 6]. However, this mild-slope equation is essentially of elliptic type and therefore its solution is very difficult due to the necessity of giving boundary conditions along the closed curve and the requirement of solving a large number of simultaneous equations. The latter restriction makes it computationally costly to find the solution for short waves over a larger area. To overcome this problem, in this research the mild-slope equation is split in tree equations to perform the calculations of the phase velocity and wave height.

The proposed numerical scheme requires the solution of a tridiagonal and pentadiagonal systems to solve the phases and the wave height, respectively. The type of grid used is the staggered cell with the possibility of a grid refinement in zones of interest. The numerical code developed was coupled with a hydrodynamic code [7] which provides the hydrodinamic in the study area. The numerical code was proved with examples where the analytical solution exists and has constant conditions obtaining acceptable values with an error less than 1%.

The study area corresponds to the breakwater of the Laguna Verde Nuclear Power Station located on the coast of the Gulf of Mexico, in Alto Lucero, Veracruz (LVNPP). The real meteorological parameters and the waves direction information of the area were acquired in order to obtain more accurate results.

2. Mathematical Model

The propagation of periodic, small amplitude, linear surface waves over a seabed of mild slope can be described by the following equation [8]: ∇𝐶𝐶_𝑔−𝐶∇𝜂^𝑔𝐶𝜕²𝜂𝜕𝑡²=0, (1) where 𝜂 is the wave height, 𝐶 is the phase velocity, and 𝐶_𝑔 is the group velocity.

Equation (1) accounts for the combined effects of refraction, diffraction, and reflection and can be used for the computation of wave propagation in coastal regions and for the computation of wave transformation with a reflecting boundary such as harbour or breakwater [9]. However, it is essentially for elliptic form, and its solution needs a boundary condition along the closed curve and a great amount of computing time and storage. Therefore for short waves over an open coast area, it is infeasible.

The aim of the numerical solution proposed for (1) is to split the equation, in its phase velocities components 𝑥 and 𝑦 and the wave height 𝜂. The equations for the phase velocities 𝑄_𝑥 and 𝑄_𝑦 in the directions 𝑥 and 𝑦, respectively, are expressed as 𝐶²𝜕𝜂𝜕𝑥=−𝜕𝑄_𝑥𝐶𝜕𝑡₂,𝜕𝜂𝜕𝑦=−𝜕𝑄_𝑦𝜕𝑡. (2)(3)

The wave height is determined by 𝜕𝜂1𝜕𝑡=−𝑛𝜕𝜕𝑥𝑛𝑄_𝑥+𝜕𝜕𝑦𝑛𝑄_𝑦, (4) where 𝑛 (group factor) is defined as follows: 1𝑛=21+2𝑘ℎ(),senℎ2𝑘ℎ (5) where 𝑘 denotes the wave number and ℎ the depth.

The advantage of dividing the general (1) in (2), (3), and (4) is that it models the presence of breakwaters in a simple manner.

3. Numerical Scheme

The numerical solution of (2), (3), and (4) is performed using a staggered grid illustrated in Figure 1, in the grid wave height is calculated at the solid grid points, labeled (𝑖−1,𝑗), (𝑖,𝑗), (𝑖+1,𝑗), (𝑖,𝑗−1), and (𝑖,𝑗−1), and so forth, and the velocities are calculated at the open grid points, labeled (𝑖−1/2,𝑗), (𝑖+1/2,𝑗), (𝑗+1/2,𝑖), and (𝑗−1/2,𝑖). Specifically, 𝑥 is calculated at points (𝑖−1/2,𝑗) and (𝑖+1/2,𝑗) and 𝑦 is calculated at different points (𝑖,𝑗+1/2) and (𝑖,𝑗−1/2). The key feature here is that wave height is calculated at different grid points. In Figure 1, the open grid points are shown equidistant between the solid grid points, but this is not a necessity.

Figure 1 

Staggered grid.

The discretization of (2) proceeds as follows: 1𝐶²_{𝑖−1,𝑗+1/2}𝑄𝑥^𝑡+1_{𝑖−1,𝑗+1/2}+22𝐶²_{𝑖,𝑗+1/2}𝑄𝑥^𝑡+1_{𝑖,𝑗+1/2}+1𝐶²_{𝑖+1,𝑗+1/2}𝑄𝑥^𝑡+1_{𝑖+1,𝑗+1/2}=1𝐶²_{𝑖−1,𝑗+1/2}𝑄𝑥^𝑡_{𝑖−1,𝑗+1/2}+22𝐶²_{𝑖,𝑗+1/2}𝑄𝑥^𝑡_{𝑖,𝑗+1/2}+1𝐶²_{𝑖,𝑗+1/2}𝑄𝑥^𝑡_{𝑖+1,𝑗+1/2}−24𝜂Δ𝑥^{𝑡𝑖+1/2,𝑗+1/2}−𝜂^𝑡_{𝑖−1/2,𝑗+1/2}; (6) analogous for (3) the discretization is expressed as 1𝐶²_{𝑖+1/2,𝑗−1}𝑄𝑦^𝑡+1_{𝑖+1/2,𝑗−1}+22𝐶²_{𝑖+1/2,𝑗}𝑄𝑦^𝑡+1_{𝑖+1/2,𝑗}+1𝐶²_{𝑖+1/2,𝑗+1}𝑄𝑦^𝑡+1_{𝑖+1/2,𝑗+1}=1𝐶²_{𝑖+1/2,𝑗−1}Q𝑦^𝑡_{𝑖+1/2,𝑗−1}+22𝐶²_{𝑖+1/2,𝑗}𝑄𝑦^𝑡_{𝑖+1/2,𝑗}+1𝐶²_{𝑖+1/2,𝑗+1}𝑄𝑦^𝑡_{𝑖+1/2,𝑗+1}−24𝜂Δ𝑦^{𝑡𝑖+1/2,𝑗+1/2}−𝜂^𝑡_{𝑖+1/2,𝑗−1/2}. (7)

The solution of (6) and (7) produces a tridiagonal system for every row 𝑖 and 𝑗 column, this implies that we have to solve 𝑛𝑖×𝑛𝑗 lineal systems, for this reason it is important to have an efficient method, and in this research we use the Thomas algorithm. Figure 2 shows the computational molecule for the calculation of 𝑄𝑥 and 𝑄𝑦.

Figure 2 

Computational module for the calculation of 𝑄𝑥 and 𝑄𝑦.

The approximation for the wave height is given by 𝑛_{𝑖−1/2,𝑗+1/2}𝜂^𝑡+1_{𝑖−1/2,𝑗+1/2}+𝑛_{𝑖+3/2,𝑗+1/2}𝜂^𝑡+1_{𝑖+3/2,𝑗+1/2}+𝑛_{𝑖+1/2,𝑗−1/2}𝜂^𝑡+1_{𝑖+1/2,𝑗−1/2}+𝑛_{𝑖+1/2,𝑗+3/2}𝜂^𝑡+1_{𝑖+1/2,𝑗+3/2}+20𝑛_{𝑖+1/2,𝑗+1/2}𝜂^𝑡+1_{𝑖+1/2,𝑗+1/2}=𝑛_{𝑖−1/2,𝑗+1/2}𝜂^𝑡_{𝑖−1/2,𝑗+1/2}+𝑛_{𝑖+3/2,𝑗+1/2}𝜂^𝑡_{𝑖+3/2,𝑗+1/2}+𝑛_{𝑖+1/2,𝑗−1/2}𝜂^𝑡_{𝑖+1/2,𝑗−1/2}+𝑛_{𝑖+1/2,𝑗+3/2}𝜂^𝑡_{𝑖+1/2,𝑗+3/2}+20𝑛_{𝑖+1/2,𝑗+1/2}𝜂^𝑡_{𝑖+1/2,𝑗+1/2}−Δ𝑡𝑛Δ𝑥^{𝑖+1,𝑗−1/2}𝑄𝑥^𝑡+1_{𝑖+1,𝑗−1/2}−𝑛_{𝑖,𝑗−1/2}𝑄𝑥^𝑡+1_{𝑖,𝑗−1/2}+𝑛_{𝑖+1,𝑗+3/2}𝑄𝑥^𝑡+1_{𝑖+1,𝑗+3/2}−𝑛_{𝑖,𝑗+3/2}𝑄𝑥^𝑡+1_{𝑖,𝑗+3/2}−Δ𝑡𝑛Δ𝑦^{𝑖−1/2,𝑗+1}𝑄𝑦^𝑡+1_{𝑖−1/2,𝑗+1}−𝑛_{𝑖−1/2,𝑗}𝑄𝑦^𝑡+1_{𝑖−1/2,𝑗}+𝑛_{𝑖+3/2,𝑗+1}𝑄𝑦^𝑡+1_{𝑖+3/2,𝑗+1}−𝑛_{𝑖+3/2,𝑗}𝑄𝑦^𝑡+1_{𝑖+3/2,𝑗}−22𝑛_{𝑖+1,𝑗+1/2}𝑄𝑥^𝑡+1_{𝑖+1,𝑗+1/2}−𝑛_{𝑖,𝑗+1/2}𝑄𝑥^t+1_{𝑖,𝑗+1/2}+𝑛_{𝑖+1/2,𝑗+1}𝑄𝑦^𝑡+1_{𝑖+1/2,𝑗+1}−𝑛_{𝑖+1/2,𝑗}𝑄𝑦^𝑡+1_{𝑖+1/2,𝑗}. (8) The last (8) is fourth-order-accurate and produces a pentadiagonal banded lineal system. Such system requires the previous calculated phases (𝑄𝑥;𝑄𝑦) and the wave height in the previous time. We utilize a sparse matrix iterative method to solve the system. The values of 𝜂, 𝑄𝑥, and 𝑄𝑦 could be zero at time 𝑡 except at the boundaries where the waves are incidents; this condition is necessary to determine the new values at time (𝑡+Δ𝑡). Figure 3 shows the computational molecule for the calculation of the wave height.

Figure 3 

Computational module for the calculation of the wave heights.

4. Boundary Conditions

To obtain the wave propagation at open boundaries, we utilize Snell's law, which permits that the waves leave the domain without modifying their direction nor magnitude. The obstacles and buildings inside the domain are considered wall boundaries with a reflection factor (𝐾_𝐿), ranging from 0.0 (no reflection) to 1.0 (completely reflected).

To calculate the wave height at the borders of any breakwater, we formulate two expressions: 𝜂₀^𝑛+1=1+𝐾_𝐿∗𝜂^𝑛_{()𝐶Δ𝑡/Δ𝑥}−𝐾_𝐿𝜂₀^𝑛−1,𝜂_𝐿^𝑛+1=1+𝐾_𝐿∗𝜂^𝑛_{()𝐿−𝐶Δ𝑡/Δ𝑥}−𝐾_𝐿𝜂₀^𝑛−1. (9)(10)

Equation (9) is used at the start of the breakwater of any length and (10) is used at the end. Figure 4 shows how the boundary conditions are applied.

Figure 4 

Boundary conditions used in the open domain and near of the breakwater.

4.1. Energy Dissipation

For the wave height of the quasioscillating waves produced by the shock of the incident waves of height (𝐻_𝑖) versus the reflecting waves of height (𝐻_𝑟), the new wave height could be estimated as 𝜂^𝑛+1_𝑖,𝑗=𝐻_𝑖,𝑗2cos𝐾𝑥_𝑖,𝑗𝑋_𝑖cos𝜃_𝑖,𝑗+𝐾𝑦_𝑖,𝑗𝑌_𝑗sin𝜃_𝑖,𝑗−2𝜋𝑇∗𝑡, (11) where 𝐾𝑥 is the factor of group in 𝑥, 𝐾𝑦 factor of group in 𝑦, 𝑋 is the distance with respect to 𝑥-axis, and 𝑌 the distance with respect to 𝑦-axis.

At the wall, the proper boundary condition for (6) at 𝑖=1 and 𝑖=𝑛𝑖 along the line (𝑦=𝑗×Δ𝑦+1/2) is determined as 𝑄𝑥^𝑡+1_1,𝑗+1/2=𝑄𝑥^𝑡_1,𝑗+1/2−𝐶²_1,𝑗+1/2Δ𝑡𝜂Δ𝑥^{𝑡5/2,𝑗+1/2}−𝜂^𝑡_3/2,𝑗+1/2,𝑄𝑥^𝑡+1_{𝑛𝑖,𝑗+1/2}=𝑄𝑥^𝑡_{𝑛𝑖,𝑗+1/2}−𝐶²_{𝑛𝑖,𝑗+1/2}Δ𝑡𝜂Δ𝑥^{𝑡𝑛𝑖,𝑗+1/2}−𝜂^𝑡_{𝑛𝑖−1,𝑗+1/2}. (12)

Analogous for (7) at 𝑗=1 and 𝑗=𝑛𝑗 along the line (𝑥=𝑖×Δ𝑥+1/2) are used the next equations: 𝑄𝑦^𝑡+1_𝑖+1/2,1=𝑄𝑦^𝑡_𝑖+1/2,1−𝐶²_𝑖+1/2,1Δ𝑡𝜂Δ𝑦^{𝑡𝑖+1/2,5/2}−𝜂^𝑡_𝑖+1/2,3/2,𝑄𝑦^𝑡+1_{𝑖+1/2,𝑛𝑗}=𝑄𝑦^𝑡_{𝑖+1/2,𝑛𝑗}−𝐶²_{𝑖+1/2,𝑛𝑗}Δ𝑡𝜂Δ𝑦^{𝑡𝑖+1/2,𝑛𝑗}−𝜂^𝑡_{𝑖+1/2,𝑛𝑗−1}. (13)

Finally, the angle frequency could be estimated by the relation between the phases 𝑄𝑥 and 𝑄𝑦 as 𝜔=tan⁻¹𝑄𝑦^𝑛+1𝑄𝑥^𝑛+1. (14)

5. Controlled Numerical Experiments

To carry out the validation of the numerical code, theoretical examples were used in order to compare the results produced by the effects of refraction, diffraction, and reflection. In our first example we want to reproduce the diffraction and reflection caused by the presence of two symmetrical breakwaters with an opening in the middle.

Figure 5 depicts the geometry of the obstacles and the expected waves diagram at the exit of the opening.

Figure 5 

Diffraction-Reflection diagram produced by two breakwaters.

In Figure 6 we can see the effects of diffraction and reflection, where the behaviour of the waves agrees with the expected solution depicted in Figure 5. The size of each grid cell is Δ𝑥=Δ𝑦=10 m and the time step Δ𝑡 is of 1  seconds, the grid size is of 110×110 and the breakwaters are located in the middle of the grid with an opening of 10 m.

(a)

(b)

(c)

(d)

Figure 6 

Snapshots of the normalized wave height at time 16  seconds (a), 34  seconds (b), 50  seconds (c), and 64 seconds (d).

In our next example we model the wave transformation produced by the breakwaters of a small rectangular harbour, depicted in Figure 7; the experimental solution (reference solution) for this problem was presented by Unluate and Mei [10] and a numerical solution by Lee [11] and J. Maa et al. [12].

Figure 7 

Small harbour geometry.

The miniexperimental harbour geometry is of 0.3212 m in length × 0.0605 m of wide and 0.2576 m of depth, the angle of incident of the waves is zero degrees, the initial wave height is 0.01 m, and the initial parameters summary of the numerical test is depicted in Table 1.

Table 1 

Parameters used for the small harbour.

Figure 8 shows the contour of the wave height obtained with the numerical code developed when the waves are stable; the boundary condition at the wall was considered with a reflection coefficient of 1.0 that means a total reflecting wall.

Figure 8 

Relative wave height distribution into the small harbour.

Later a quantitative comparison between the solutions obtained by Unluate and Mei [10] and J. Maa et al. [12] was made. Figure 9 shows normalized wave height obtained at the receiver compared with the reference solution and the numerical solution obtained by J. Maa et al. with a max relative error of 0.75%; the technique used to approximate the error could be seen in [13].

Figure 9 

Time evolution of the wave height at the numerical receiver.

6. Numerical Experiments in the Breakwater of the LVNPP

The numerical simulation was carried out in the coastal area of the LVNPP which is located on the coast of the Gulf of Mexico, in Alto Lucero, Veracruz. It is the only electric power generating nuclear plant in Mexico and produces about 4.5% of the country's electrical energy. LVNPP has an original installed capacity of 1,365 Megawatts (Mw). It consists of two units General Electric Boiling Water Reactors (BWR-5) using Uranium (U235 Isotope 3% enriched) as fuel. Unit-1 (U-1) started its operation on July 29, 1990. Unit-2 (U-2) started its operation on April 10, 1995. The plant is owned and operated by the national electric company owned by the Mexican government (see Figure 10).

Figure 10 

Location map of the study area and the location of the weather stations.

The generation of electric power at the plant is based on the technology of nuclear fission of uranium atoms, which takes place in the reactor. The energy released by the nuclear fission is transferred as heat from the fuel to the cooling water, which boils into steam. The water discharge of the cooling system is carried out in a trough channel with a flow rate of 63 m³/s and a mean speed of 1.4 m/s.

The type of the cooling system is considered open, actually the intake of the water is superficial, the breakwater which protects the intake has two piers, the main problem is the blockage, and the recirculation of hot water has causing deficiencies in the cooling system.

The grid of the study domain is shown in Figure 11, it is conformed by 5,896 cells with a constant size of Δ𝑥=Δ𝑦=250 m.

(a)

(b)

Figure 11 

Grid map (a) and the bottom topography contour (b).

The initial parameters for this simulation are defined in Table 2.

Table 2 

Parameters used for the coastal area of the LVNPP.

In Figure 12 is shown the result for the coastal area, we can see the refraction effects when the waves are near the coastal line produced by the bottom topography.

Figure 12 

Wave height distribution in the coastal zone.

Studying in detail the behaviour of waves and their circulation inside the breakwater is very important to study the sedimentation effects, because the sediments transports have the potential to increase the maritime life producing the growing of different animals like mollusc, shells, and clams. This issue could cause the siltation and obstruction of the intake area reducing the performance of the cooling system.

Due to the importance of the breakwater area it is necessary to make a refinement of the grid into the breakwater; this refinement is made with 600 cells with a size of Δ𝑥=Δ𝑦=5 m. In Figure 13 is shown the grid and the bathymetry of the breakwater area.

(a)

(b)

Figure 13 

Grid map (a) and the bottom topography contour (b) of the breakwater area.

To probe the model, as first step, is proposed a breakwater geometrically idealized (see Figure 14); this means a controlled scenario to represent the wave transformations into the breakwater. This idealized breakwater is very close to the way and localization of the real breakwater to observe the behavior of wave propagation.

Figure 14 

Geometry of the idealized breakwater.

The results of the wave propagation are shown in Figure 15 for different simulation times; we can see that the diffraction and reflection phenomena are diminished by the breakwater structure.

(a)

(b)

(c)

(d)

Figure 15 

Snapshots of the wave height distribution into the breakwater idealized at 16 seconds (a), 34  seconds (b), 50 seconds (c) and 64 seconds (d).

Once the simulations over the idealized breakwater were done and show results with well-defined, patterns we carry out numerical experiment with the real bathymetry and the real breakwater form. The simulation result of the stable waves is depicted in Figure 16.

Figure 16 

Stable wave height distribution into the breakwater.

In Figure 17 we can see that the waves height into the breakwater is minimal; we find just a small disturbance in the free surface ranging from 0.05 to 0.1 m; this pattern is consistent with the average data acquired in the meteorological stations.

Figure 17 

Stable wave heights contour distribution into the breakwater.

7. Final Remarks

In this research we develop a very fast numerical code for the study of the wave propagation over irregular bottom topographies. The model has produced consistent results with controlled cases and the result for the coastal area of the LVNPP is very acceptable compared with the data acquired in situ.

The aim of the model is to split the governing equation, in its phase velocities components and wave height, reducing significantly the computing time and the complexity of the implementation using sparse matrices. We can solve the system over a grid of 5896 cells in just 120 s using sparse matrix algorithms in a standard workstation.

The model is capable of reproducing the effects of the reflection and diffraction with accuracy and the refraction could be reproduced only for domains with cells greater or equal to 10² m. In future work we will present the sediment transport and particle transport coupled with the wave transformations.

Nomenclature

𝐶:	Phase velocity (m/s)
𝑄𝑥:	Phase velocity in the 𝑥 direction (m/s)
𝑄𝑦:	Phase velocity in the 𝑦 direction (m/s)
𝜂:	Wave heights (m)
𝜔:	Angle of frequency (angle in radians)
𝜑:	Wave potential function
𝑛:	Group factor
𝐾𝐿:	Reflection factor.

Acknowledgments

The authors would like to thank engineer Ivan Campos for the information he provided to calibrate the numerical code. The authors would like to thank the reviewers for their helpful comments to improve the manuscript.