Numerical Modeling of Tsunami Waves Interaction with Porous and Impermeable Vertical Barriers

Abstract

Tsunami wave interaction with coastal regions is responsible for very important human and economic losses. In order to properly design coastal defenses against these natural catastrophes, new numerical models need to be developed that complement existing laboratory measurements and field data. The use of numerical models based on the Navier-Stokes equations appears as a reasonable approach due to their ability to evaluate complex flow patterns around coastal structures without the inherent limitations of the classical depth-averaged models. In the present study, a Navier-Stokes-based model, IH-3VOF, is applied to study the interaction of tsunami waves with porous and impermeable structures. IH-3VOF is able to simulate wave flow within the porous structures by means of the volume-averaged Reynolds-averaged Navier-Stokes (VARANS) equations. The equations solved by the model and their numerical implementation are presented here. A numerical analysis of the interaction of a tsunami wave with both an impermeable and porous vertical breakwater is carried out. The wave-induced three-dimensional wave pattern is analysed from the simulations. The role paid by the porous media is also investigated. Finally, flow around the breakwater is analyzed identifying different flow behaviors in the vicinity of the breakwater and in the far field of the structure.

1. Introduction

The interaction of tsunami waves with the coast has become of great interest in the last years due to the devastating effects observed in the last large events: The Indian Ocean and Japan tsunamis. Tsunami wave effects leaded to huge human loses and turned out to be expensive in natural resources. Moreover, the economic impact of derived effects is huge, and coastal areas affected by the tsunami wave attack need much time and a lot of economical resources to be recovered. Most of the existing coastal defenses are designed to deal with wave storms but not with tsunami wave attack. Tsunami waves differ significantly in nature from storm waves, because they are characterized by longer wave lengths and they show a clear transient effect. Traditional designs for coastal structures are not valid to provide the required protection, and new defenses are studied in order to protect the coast.

Several approaches have been followed in the literature to study tsunami wave action on coastal structures, see [1] for a review of the state-of-the-art. The lack of knowledge in the interpretation of many of the factors present in the tsunami wave interaction has motivated their study by means of physical tests, especially small-scale model tests. Formulations derived from these experimentations are, in most of the cases, semiempirical in nature with their form based on physical considerations. Moreover, scale factors are present in formulations, and their role for dissipation mechanisms due to wave breaking, turbulence and generation of eddies in the fluid region as well as turbulence and friction within the porous material, is still unsolved.

The use of numerical models in the study of tsunami wave propagation is quite popular by means of the depth averaged nonlinear shallow water (NSW) equations (i.e., COMCOD, [2–4]). NSW set of equations has revealed, as a powerful tool, to be used in large domain areas, where the long wave approach is still valid. In the vicinity of coastal structures, the inherent assumptions of NSW equations are violated because of the relevance of the vertical flow component and the nonhydrostatic behavior of the wave induced pressure.

In the last decade, the use of the Navier-Stokes (NS) equations applied to coastal engineering processes has become more popular. The increment of the computational resources and the improvement of the numerical aspects mainly related with the boundary conditions has made possible to overcome the inherent limitations of classical depth averaged models. NS models are free of the simplifications behind wave theories, and they are able to deal with the flow complexity derived from the wave interaction with coastal structures. Two-dimensional Reynolds-averaged Navier-Stokes (RANS) models [5–7] have revealed that structural functionality and stability can be studied with a high degree of accuracy, even in the presence of granular material layers. Volume-averaged Reynolds-averaged Navier-Stokes (VARANS) equations have been solved to characterize wave induced flows within porous structures. Moreover, tsunami wave transformation has been studied [8] showing a high degree of accuracy in reproducing nonlinear flow characteristics.

Recently, VARANS equations have been extended to three-dimensional problems [8, 9] showing a high degree of accuracy in predicting magnitudes related to the functionality and the stability of coastal porous structures. Wave transformation processes around coastal structures, such as wave reflection, wave penetration through porous structures, wave diffraction, run-up, and wave breaking can be now analyzed in three-dimensional problems. The model presented in [9], called IH3-VOF, appears to be suitable to be used in the analysis of tsunami wave induced flow around coastal structures, even in the presence of porous media and a two-phase flow. In the present paper, an analysis of tsunami wave interacting with coastal structures is done. The specific case of vertical walls is chosen as a first attempt to deal with more complex breakwater configurations. Both impermeable and porous breakwaters are studied, in order to identify tsunami-induced hydrodynamic.

The work is organized as follows. The mathematical model followed by IH-3VOF is presented first. A detailed description of the numerical implementation and the tsunami wave generation by the numerical model is described in Section 3. Next, numerical simulations are shown and described, studying the more relevant tsunami wave induced processes around vertical porous and impermeable walls. Finally, conclusions about the work are drawn.

2. Mathematical Modeling

IH-3VOF makes use of the volume-averaged Reynolds-averaged Navier–Stokes (VARANS) equations to solve for the pressure and velocity fields of the flow. These equations derive from the Reynolds-averaged Navier-Stokes (RANS) equations, by application of a volume-averaging procedure. This additional averaging of the equations eliminates the requirement of a complete geometrical description of the flow domain and therefore allows to use the equations to solve porous media flow.

In the present section, the main results of the mathematical model are presented. For a complete description see [10] or [9].

2.1. Volume-Averaging

Volume-averaging is a filtering operation, mathematically defined as ⟨𝑎⟩^𝑓=1𝑉_𝑓𝑉𝑓𝑎𝑑𝑉, (2.1) where ⟨⟩^𝑓 is the volume-averaging operator, 𝑎 is the variable to be averaged, and 𝑉_𝑓 is the averaging volume, the volume within which the operation is carried out. This operator is applied to the original RANS equations to obtain the VARANS equations.

It must be noted that the 𝑓 superscript indicates that the volume average is calculated only on the fluid part of the averaging volume. The so-calculated average receives the name of intrinsic average. The total averaging volume 𝑉 presents the same size at every point. However, the fluid part of the averaging volume 𝑉_𝑓 presents a variable size, depending on the amount of solids within the averaging volume. Calculating the intrinsic volume average, a coherent metric for the momentum is kept.

The drawback of the intrinsic averaged variables is that they contain information on two different variations, the variations of the variable of interest itself and the variation of the size of 𝑉_𝑓. This way, a gradient of an intrinsic averaged variable may be due only to variations of 𝑉_𝑓 and not to the variable itself. To overcome this issue extended averaged magnitudes are introduced.

Extended averages magnitudes are defined to verify the following equation: ⟨𝑎⟩^𝑓𝑉_𝑓=_𝑉𝑓𝑎𝑑𝑉=⟨𝑎⟩𝑉, (2.2) where 𝑉 is the complete averaging volume and ⟨𝑎⟩ is the extended averaged variable. As 𝑉 is constant all over the domain, all the variations are accounted by the extended averaged variable. The relation between the intrinsic and extended averages can be done by means of the porosity. The porosity is the ratio between the volume of fluid contained in the porous solid volume and the total volume of the porous body itself. Mathematically, 𝑉𝜙=^𝑓𝑉. (2.3) Introducing this last relation in (2.2), ⟨𝑎⟩=𝜙⟨𝑎⟩^𝑓, (2.4) the equation relating intrinsic and extended averages is obtained. This relation is used to transform the equations in the following pages from their intrinsic averaged form to the extended averaged one.

An important theorem is to be presented for the reader to completely understand the volume averaging process. It is the theorem for the local volume average of a gradient. A complete derivation and proof can be found in [11]. It states that the local volume averaged of a gradient is equal to the gradient of the volume averaged variable plus the integral of the flux of the variable divided by the averaging volume. Mathematically, 𝜕𝜕𝑥_𝑖𝑎^𝑓=𝜕𝜕𝑥_𝑖⟨𝑎⟩^𝑓+1𝑉_{𝑓𝜕𝑉𝑓}𝑎𝐝𝐒, (2.5) where 𝑎 is the variable which gradient is to be volume averaged, 𝜕𝑉_𝑓 is the solid part of the surface enclosing the averaging volume, and 𝐝𝐒 is a surface differential element. In the next section, the theorem is applied.

2.2. Model Equations

In the present section, the model equations are presented. They consists of the volume-averaged mass and momentum conservation equations as well as the volume-averaged turbulence equations and the porous media closure models for all of them.

2.2.1. Mass Conservation

The volume-averaged expression of mass conservation is shown in the following equation: 𝜕𝜕𝑥_𝑖⟨𝑢_𝑖⟩𝜙=0. (2.6) This equation expresses that when dealing with volume-averaged magnitudes, the divergence-free field is not the velocity, but instead the field obtained by dividing the volume-averaged velocity by the porosity field.

2.2.2. Momentum Conservation

The final expression momentum conservation equations are presented as follows: 𝜕𝜕𝑡𝑢_𝑖+𝑢_𝑗𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖𝜙=−𝜌𝜕𝜕𝑥_𝑖1𝜙𝑝+𝜙𝑔_𝑖𝜕+𝜙𝜕𝑥_𝑖𝜈𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖𝜕−𝜙𝜕𝑥_𝑗𝑢_𝑖′𝑢_𝑗′^𝑓𝜕−𝜙𝜕𝑥_𝑗𝑢_𝑖^′′𝑢_𝑗^′′𝑓−1𝑉_𝜕𝑉𝑠1𝜌𝑝^𝑓+𝑝^′′+1𝐧𝑑𝐴𝑉_𝜕𝑉𝑠𝜈𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖+𝜕𝜕𝑥_𝑗𝑢_𝑖^′′𝐧𝑑𝐴. (2.7)

The first line of (2.7) is the expression of the Navier-Stokes with the averaged variables. The second one shows the volume-averaged Reynolds stresses that appear due to turbulence. Finally the third and fourth lines contain the terms that appear due to the volume average process. It must be noted that in the case of volume-averaging a solid-free domain, both integrals vanish, remaining only the effect of the nonresolved subgrid velocities. These integrals introduce the forces exerted by the solid objects against the fluid in the form of pressure and tangential forces.

2.2.3. Turbulence Modeling

Volume-averaging affects all the model equations, including the turbulence model. In this case a 𝑘-𝜖 model is volume-averaged in order to provide a closure for the turbulent terms in the simulations.

The volume-averaged 𝑘 equations is presented as follows: 𝜌𝜕⟨𝑘⟩+𝜕𝑡𝑢_𝑗𝜕𝜕𝑥_𝑗1𝜙⟨𝑘⟩=𝜙𝜏_𝑖𝑗^𝑓𝑆_𝑖𝑗^𝑓⟨𝜖⟩𝜕−𝜌+𝜙𝜕𝑥_𝑗𝜇𝜇+^𝑇𝜙𝜎_𝑘𝜕𝜕𝑥_𝑗1𝜙⟨𝑘⟩+𝜙[𝐶𝑇]_𝑘, (2.8) where the production term can be developed as ⟨𝜏_𝑖𝑗⟩^𝑓𝑆_𝑖𝑗^𝑓=2𝜇_𝑇^𝑓𝑆_𝑖𝑗^𝑓−23𝜌⟨𝑘⟩^𝑓𝛿_𝑖𝑗𝑆_𝑖𝑗^𝑓=2𝜙𝜇_𝑇𝑆_𝑖𝑗^𝑓−23𝜌𝜙⟨𝑘⟩𝛿_𝑖𝑗𝑆_𝑖𝑗^𝑓, (2.9) being the intrinsic average of the rate of strain tensor: 𝑆_𝑖𝑗^𝑓=12𝜕𝑢_𝑖𝜕𝑥_𝑗+𝜕𝑢_𝑗𝜕𝑥_𝑖^𝑓=12𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖+𝜕𝜕𝑥_𝑖1𝜙𝑢_𝑗, (2.10) the extended average of the turbulent viscosity is expressed as 𝜇_𝑇=𝜌𝐶_𝜇⟨𝑘⟩²⟨𝜖⟩, (2.11) and [𝐶𝑇]_𝑘 are the terms of (2.8) that will be modeled with a closure model.

Volume-averaging the 𝜖 equation𝜌𝜕⟨𝜖⟩+𝜕𝑡𝑢_𝑗𝜕𝜕𝑥_𝑗1𝜙⟨𝜖⟩=𝜙𝐶_𝜖1⟨𝜖⟩⟨𝑘⟩𝜏_𝑖𝑗^𝑓𝑆_𝑖𝑗^𝑓−𝐶_𝜖2𝜌⟨𝜖⟩²⟨𝑘⟩𝜕+𝜙𝜕𝑥_𝑗𝜇𝜇+^𝑇𝜙𝜎_𝜖𝜕𝜕𝑥_𝑗1𝜙⟨𝜖⟩+𝜙[𝐶𝑇]_𝜖, (2.12) where ⟨𝜏_𝑖𝑗⟩^𝑓,⟨𝑆_𝑖𝑗⟩^𝑓, and ⟨𝜇_𝑇⟩ are expressed as per (2.9), (2.10), and (2.11), and [𝐶𝑇]_𝜖 are the terms that will be modeled with the closure model.

2.2.4. Porous Media Closure

Volume-averaging the momentum conservation equations eliminates the need of a detailed description of the porous medium geometry in order to carry out simulations. This simplification on the initial requirements, however, introduces new terms in the equations that need to be modeled. These new terms conform the last line of (2.7).

These terms that represent the effect introduced in the flow by the solid matrix are modeled with the Forchheimer model. This model consists of three terms: a linear term introducing viscous effects: a nonlinear term, proportional to the square of the velocity, taking into account turbulent effects and other non linear interactions; an inertia term that accounts for the added mass effect. Mathematically, it can be written as follows: [𝐶𝑇]=𝑎𝑢_𝑖+𝑏𝑢_𝑖𝑢_𝑖𝜕+𝑐𝜕𝑡𝑢_𝑖. (2.13)

[𝐶𝑇] makes reference to all these closure terms. 𝑎,𝑏, and 𝑐 are constants that depend on flow and porous medium characteristic and must be determined empirically. In the state-of-the-art chapter the relation of this coefficients with porous media characteristic is given.

Engelund's [12] formulas are used for the Forchheimer model. The closure terms for the 𝑘-𝜖 equations are modeled by means of the closure model presented in [13] [𝐶𝑇]_𝑘=𝜖_∞[𝐶𝑇]_𝜖=𝐶_𝜖2𝜖²_∞𝑘_∞, (2.14)𝜖_∞ and 𝑘_∞ being 𝑘_∞=3.71−𝜙𝜙^3/2_𝑖𝑢_𝑖²𝜖_∞=39.01−𝜙^5/2𝜙²_𝑖𝑢_𝑖^23/21𝐷₅₀, (2.15) where 𝜙 is the porosity and 𝐷₅₀ the mean pore size.

3. Numerical Implementation

In coastal engineering problems, free surface gravity waves are present in almost every application. All the interactions and generated dynamics are conditioned by the free surface evolution; hence, it is a key aspect of a coastal engineering model. In order to deal with free surface flows in IH3VOF, a free surface tracking algorithm is needed. Among all the possible methods, the volume of fluid one [14] was chosen. It is used together with the volume tracking algorithm presented in [15].

The main reason for choosing the volume of fluid method is that it allows the representation of very complex free surface geometries in a very neat form. It can also be used to consider all the interfaces between an arbitrary number of fluids in the computational domain by adding little complexity to the algorithm. Moreover, the equation describing the evolution of the free surface is a simple transport equation that shares operators with the Navier–Stokes equations, simplifying its implementation. Its performance highly depends on the volume tracking algorithm used. The one presented in [15] represents the interface as a piecewise linear polynomial in every cell, allowing a very accurate capture of the free surface geometry. This accurate representation is also used to quantify the mass and momentum fluxes across cell faces, what drastically reduces the numerical diffusion of the free surface. The choice of this volume tracking method configures a very robust and accurate representation of the free surface. The main drawback being the computational cost that is a bit high in comparison with other tracking methods.

In order to maintain the performance of the model within the limits that engineering practice demands, an explicit method in time for the integration of the momentum equations is required. For this reason, a fractional two-step method [16] was chosen. This method consists in solving the equations in two different steps. First, a prediction of the new time step velocity is obtained. The prediction obtained is nondivergence free velocity field that must be corrected. This step is called predictor. Then, the pressure field is solved imposingthatthe new time step velocity field must be solenoidal. This second step is called projection because the equation linking the pressures and predicted velocities is obtained by a projection procedure. Obviously, as the partial differential equation representing the coupling of velocity and pressure is elliptical, it is implicit in space. Consequently, a system of equations must be solved.

Turbulence is modeled with a first-order accurate explicit scheme in time and space. This makes the solution for the turbulent variables less accurate of what it could be achieved with implicit or semiimplicit schemes. This choice is made mainly for three reasons. The first one is that due to the size of the domains where coastal engineering problems are solved, boundary layers and small scale details of the flow cannot be solved for. This already implies a reduction on the quality of the solutions, independently of the numerical schemes used. Second reason is that turbulence is not the key dynamics of coastal engineering problems. It has been already said that free surface evolution is the most important dynamics and therefore is the one that must concentrate the biggest efforts. Final reason is that explicit schemes are computationally cheaper than implicit ones and having chosen a slower free surface tracking method, turbulence must be solved with a faster method.

Assuming the solution for time 𝑛 is known, all the methods enumerated in previous paragraphs are implemented into a general algorithm that is described as follows.(1) The flow domain for time 𝑛+1 is calculated by means of the VOF method. It makes use of time 𝑛 velocities. The new flow domain is represented by means of time 𝑛+1 cell densities, 𝜌^𝑛+1. (2) Turbulence equations are solved over the new domain with time 𝑛 velocities. Turbulent viscosity is calculated previous to predictor step. (3) Velocity is calculated at an intermediate time ∗ (𝑢^∗) as a prediction of the new time velocity (𝑢^𝑛+1). The predictor step accumulates the momentum contribution of all the terms in the Navier-Stokes equations but the new time pressure gradient term. (4) The increment of pressure gradient at time 𝑛+1 is calculated in order to convert 𝑢^∗ into a solenoidal velocity field. (5) With the time 𝑛+1 pressure gradients, the pressure for time 𝑛+1 is calculated. The new pressure gradients are also applied to the predicted velocities (𝑢^∗) to obtain the new time velocities (𝑢^𝑛+1). (6) The process is carried out until the required final time is reached.

3.1. Free Surface Modeling

The main characteristic of free surface flows is that one of the boundaries, that receives the name of free surface, is not fixed. It moves with the flow. This moving boundary appears as a result of having more than one fluid in the solution domain. Free surface modeling makes part of multiphase flow modeling. Therefore, in addition to the problem of integrating the Navier-Stokes equations, free surface flows present the problem of determining the real domain occupied by every one of the fluids present in the problem.

In a general case, no mathematical expression can be assumed to represent the free surface. Moreover, in complex flows, the free surface may be split in different multiconnected surfaces. Therefore, its numerical treatment as an analytical surface by means of an equation ruling its behavior is very complicated. To avoid this complication, different methods have been developed by different authors. Among them the volume of fluid [14] method is chosen to be implemented in IH3VOF. The volume of fluid method is a numerical technique for tracking the free surface.

It is based on tracking mass fluxes trough volume fractions across the mesh cells. In order to do this, a volume fraction function is defined for every material (phase) at every mesh cell. This function is defined as the volume of fluid 𝑘 contained in the cell divided by the cell volume: 𝜉_𝑘=𝑉_𝑘𝑉_𝑐. (3.1) By means of this volume fraction, the density of a cell can be directly calculated as 𝜌_𝑐=𝜉_𝑘𝜌_𝑘, (3.2) where 𝜌_𝑐 is the cell density and 𝜌_𝑘 is the density of fluid 𝑘. Einstein notation of summation over repeated indices is used.

That is, the cell density is calculated as a weighted mean, where the contribution of every phase is directly proportional to its volume fraction at that particular cell. This procedure is used for every other fluid property, for instance viscosity: 𝜈_𝑐=𝜉_𝑘𝜈_𝑘. (3.3) But it can be applied to any other fluid property just by replacing it by viscosity in (3.3).

To obtain now the evolution equation for the volume fraction, (3.2) is introduced in the mass conservation equation to obtain 𝜕𝜉𝜕𝑡^𝑘𝜌_𝑘+𝜕𝜕𝑥_𝑖𝜉_𝑘𝜌_𝑘𝑢_𝑖^𝑓=0, (3.4) a system of 𝑘 equations, one per fluid. It has been said that in coastal engineering application fluids are supposed to be incompressible, then 𝜌_𝑘 is a constant that can be taken out from last equation. Expressing the equation on extended volume average magnitudes 𝜕𝜉𝜕𝑡^𝑘+𝜕𝜕𝑥_𝑖𝜉_𝑘𝑢_𝑖𝜙=0, (3.5) the expression to implement numerically is obtained.

All the equations and methods explained to treat the free surface in the flow are also used to study two-phase flows. The VOF technique can track several fluids, reconstructing all the interfaces between them. Assigning afterwards to every cell of the domain the fluid characteristics obtained by combining the values of the different fluids in the cell, using as the weight the relative volume occupied by each phase, the flow equations can be solved. Then, every fluid phase is moved and tracked independently, and cell properties are reconstructed for a new time step.

3.2. Time Discretization

It has already been pointed out that an explicit in time numerical scheme is required. A backward step in time derivative will be used. Then, time discretization can be written as 𝜕𝜕𝑡𝑢_𝑖^𝑛+1+𝑢_𝑗𝜕1+𝛾𝜕𝑥_𝑗1𝜙𝑢_𝑖^𝑛𝜙=−𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙𝑝^𝑛+1+𝜙𝑔1+𝛾^𝑖𝑛+𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖𝜇𝜇+^𝑡𝜙𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖^𝑛−23𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖𝜌𝜙⟨𝑘⟩^𝑛−𝛼1+𝛾𝑢_𝑖^𝑛−𝛽1+𝛾𝑢_𝑖𝑢_𝑖^𝑛. (3.6)

The pressure gradient term can be split into two contributions, one at time 𝑛 and another at time 𝑛+1: 𝜕𝜕𝑥_𝑖1𝜙𝑝^𝑛+1=𝜕𝜕𝑥_𝑖1𝜙𝑝^𝑛+𝜕𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1. (3.7) This decomposition expresses the pressure gradient at final time as the pressure gradient at initial time plus a pressure variation. This pressure increment is generated by the forces acting at time 𝑛. Hence, the variation is the flow response to the external actions. This variation is the variable resolved and not the total pressure.

Finally, the time derivative can be developed as a function of three time instants: the current time (superscript 𝑛), the intermediate time used for the two-step projection methods (superscript *) and the new time (superscript 𝑛+1). The expression then results in 𝜕𝜕𝑡𝑢_𝑖^𝑛+1≈𝑢_𝑖^𝑛+1−𝑢_𝑖^𝑛=Δ𝑡𝑢_𝑖^𝑛+1−𝑢_𝑖^∗+𝑢_𝑖^∗−𝑢_𝑖^𝑛=Δ𝑡𝑢_𝑖^𝑛+1−𝑢_𝑖^∗+Δ𝑡𝑢_𝑖^∗−𝑢_𝑖^𝑛.Δ𝑡 (3.8)

Introducing expressions (3.7) and (3.8) into (3.6), the complete numerical scheme is obtained: 𝑢_𝑖^𝑛+1−𝑢_𝑖^∗+Δ𝑡𝑢_𝑖^∗−𝑢_𝑖^𝑛+Δ𝑡𝑢_𝑗𝜕1+𝛾𝜕𝑥_𝑗1𝜙𝑢_𝑖^𝑛𝜙=−𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙𝑝^𝑛−𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1+𝜙𝑔1+𝛾^𝑖𝑛+𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖𝜇𝜇+^𝑡𝜙𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖^𝑛−23𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖𝜌𝜙⟨𝑘⟩^𝑛−𝛼1+𝛾𝑢_𝑖^𝑛−𝛽1+𝛾𝑢_𝑖𝑢_𝑖^𝑛. (3.9)

This complete scheme is to be solved with a fractional method where the complete scheme is split in different steps that are solved in sequence instead of simultaneously. In this case, there are two different steps, predictor and projector. Predictor step is completely explicit and is expressed as 𝑢_𝑖^∗−𝑢_𝑖^𝑛+Δ𝑡𝑢_𝑗𝜕1+𝛾𝜕𝑥_𝑗1𝜙𝑢_𝑖^𝑛𝜙=−𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙𝑝^𝑛+𝜙𝑔1+𝛾^𝑖𝑛+𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖𝜇𝜇+^𝑡𝜙𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖^𝑛−23𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖𝜌𝜙⟨𝑘⟩^𝑛−𝛼1+𝛾𝑢_𝑖^𝑛−𝛽1+𝛾𝑢_𝑖𝑢_𝑖^𝑛. (3.10) It can readily be seen that ⟨𝑢_𝑖⟩|^∗ is an explicit function of time 𝑛 terms only.

The rest of the terms are collected for the projector step 𝑢_𝑖^𝑛+1−𝑢_𝑖^∗𝜙Δ𝑡=−𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1, (3.11) where time ∗ values are known and an implicit relation is expressed between time 𝑛+1 velocities and pressures. This expression must be projected onto a solenoidal field by taking the divergence of the expression, resulting in the equation used for the projection step, yielding a well-known Poisson-type equation 𝜕𝜕𝑥_𝑖𝑢_𝑖^∗=𝜕Δ𝑡𝜕𝑥_𝑖𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1. (3.12) Equation (3.12) states the implicit relation that exists between pressure and velocity in the Navier-Stokes equations. This implicitness cannot be avoided because both variables are strongly coupled. Time 𝑛+1 pressure turns the 𝑢^∗ velocity field into the solenoidal time 𝑛+1 velocity field. This relation is obtained by rearranging terms in (3.11), obtaining 𝑢_𝑖^𝑛+1=𝑢_𝑖^∗−Δ𝑡𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1. (3.13)

3.3. Predictor Step

Making use of the finite volume method, the spatial discretization of the predictor steps results in 𝑢_𝑖^∗=𝑢_𝑖^𝑛−Δ𝑡𝑉_𝑐𝑢_𝑗1+𝛾_𝑓𝑢_𝑖𝑓𝜙_𝑓𝑛_𝑗𝐴_𝑓−Δ𝑡𝜙𝜌11+𝛾𝑉𝑐_𝑓1𝜙_𝑓𝑝_𝑓𝑛_𝑗𝐴_𝑓+Δ𝑡𝜙𝑔1+𝛾^𝑖+Δ𝑡𝜙𝜌11+𝛾𝑉_𝑐𝑓𝜇_𝑓+𝜇_𝑡𝑓𝜙_𝑓𝜕𝜕𝑥_𝑗1𝜙𝑢_𝑖𝑓𝑛_𝑗𝐴_𝑓−23Δ𝑡𝜙𝜌11+𝛾𝑉𝑐_𝑓𝜌_𝑓𝜙_𝑓⟨𝑘⟩_𝑓𝑛_𝑗𝐴_𝑓−Δ𝑡𝛼1+𝛾𝑢_𝑖−Δ𝑡𝛽1+𝛾𝑢_𝑖𝑢_𝑖. (3.14)

3.4. Projection Step

The same procedure is now carried out on the projection step equations to obtain its finite volume spatial discretization. Integrating the left hand side of (3.11) in first place _𝑉𝜕𝜕𝑥_𝑖𝑢_𝑖^∗Δ𝑡𝑑𝑉=_𝑓𝑢_𝑖^∗_𝑓𝑛Δ𝑡^𝑖𝐴_𝑓, (3.15) and then the right hand side _𝑉𝜕𝜕𝑥_𝑖𝜙𝜌𝜕1+𝛾𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1𝑑𝑉=_𝑓𝜙_𝑓1+𝛾_𝑓𝜌_𝑓𝜕𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1_𝑓𝑛_𝑖𝐴_𝑓, (3.16) the final expression can be constructed _𝑓𝑢_𝑖^∗_𝑓𝑛Δ𝑡^𝑖𝐴_𝑓=_𝑓𝜙_𝑓1+𝛾_𝑓𝜌_𝑓𝜕𝜕𝑥_𝑖1𝜙Δ𝑝^𝑛+1_𝑓𝑛_𝑖𝐴_𝑓, (3.17) where the pressure gradient term is calculated the same way shear stresses were, face velocities are interpolated with a Rhie-Chow scheme, and the rest of face values are linearly interpolated. Once (3.17) is solved, the pressure increments at the new time-step are obtained. These increments are added to the previous time-step pressure solution obtaining the pressure solution for time 𝑛+1.

The divergence theorem can also be used to approximate the value of a divergence or gradient in the center of a cell, because a divergence term can be treated also in the general way, so _𝑉()∇𝛼𝑑𝑉=∇𝛼_𝐶𝑉_𝑐∬_𝜕𝑉𝛼_𝑓^→𝑛𝑑𝑆, (3.18) and then the gradient at the centroid can be expressed as (∇𝛼)_𝐶=∬_𝜕𝑉𝛼_𝑓^→𝑛𝑑𝑆𝑉_𝑐. (3.19)

With the updated pressures, the new pressure gradients at cell centers are calculated by means of (3.19), interpolating time 𝑛+1 pressure values to cell faces. At this point, all the terms needed to obtain the new time-step velocities have been calculated. New time-step velocities are then obtained by substitution in (3.13).

3.5. Turbulence Model

A backward step in time discretization will be used as in all the previous equations, providing an explicit scheme for both the 𝑘 and the 𝜖 equations. The spatial discretization is performed by means of the finite volume technique. The numerical scheme for the turbulent kinetic energy is then ⟨𝑘⟩^𝑛+1=𝜌^𝑛𝜌^𝑛+1⟨𝑘⟩^𝑛−𝜌^𝑛𝜌^𝑛+1𝜙𝑉_𝑐𝑓𝑢_𝑗^𝑛_𝑓⟨𝑘⟩^𝑛_𝑓𝜙²_𝑓Δ𝑡𝑛_𝑗𝐴_𝑓+12𝜌^𝑛+1𝜇_𝑇^𝑛_𝑖𝑗𝜕𝜕𝑥_𝑗𝑢_𝑖^𝑛𝜙+𝜕𝜕𝑥_𝑖𝑢_𝑗^𝑛𝜙²−𝜌^𝑛𝜌^𝑛+1⟨𝜖⟩^𝑛+𝜙𝜌^𝑛+1_𝑓𝜇𝜇+^𝑇𝑛𝜙𝜎_𝑘𝑓𝜕𝜕𝑥_𝑗⟨𝑘⟩^𝑛𝜙_𝑓𝑛_𝑗𝐴_𝑓+𝜙𝜌^𝑛+1𝜖_∞, (3.20)

and for the 𝜖 equation ⟨𝜖⟩^𝑛+1=𝜌^𝑛𝜌^𝑛+1⟨𝜖⟩^𝑛−𝜌^𝑛𝜌^𝑛+1𝜙𝑉_𝑐𝑓𝑢_𝑗^𝑛_𝑓⟨𝜖⟩^𝑛_𝑓𝜙²_𝑓Δ𝑡𝑛_𝑗𝐴_𝑓+𝐶_𝜖1⟨𝜖⟩^𝑛⟨𝑘⟩^𝑛12𝜌^𝑛+1𝜇_𝑇^𝑛_𝑖𝑗𝜕𝜕𝑥_𝑗𝑢_𝑖^𝑛𝜙+𝜕𝜕𝑥_𝑖𝑢_𝑗^𝑛𝜙²−𝐶_𝜖2𝜌^𝑛𝜌^𝑛+1⟨𝜖⟩^𝑛2⟨𝑘⟩^𝑛+𝜙𝜌^𝑛+1_𝑓𝜇𝜇+^𝑇𝑛𝜙𝜎_𝜖𝑓𝜕𝜕𝑥_𝑗⟨𝜖⟩^𝑛𝜙_𝑓𝑛_𝑗𝐴_𝑓+𝐶_𝜖2𝜙𝜌^𝑛+1𝜖_∞²𝑘_∞. (3.21)

3.6. Boundary Conditions

The model makes use of the no-slip boundary condition at solid boundaries, imposing a zero velocity at the fluid. The only exception is the inflow and outflow boundaries in which the velocity profile is specified. In Section 3.6.1, the wave generation inflow condition is presented.

At solid boundaries where the no-slip boundary condition for the velocity holds, the flow solution is forced to match the law of the wall: 𝑈=𝑢_𝜏1𝜅𝑢ln^𝜏𝑦𝜈+𝐵. (3.22) In the law of the wall equation, 𝑈 is the velocity in the first node over the surface, 𝑢_𝜏 is the friction velocity, 𝜅 is the Kármán constant, 𝑦 is the distance between the node and the surface, 𝜈 is the kinematic viscosity, and 𝐵 is the constant of the law of the wall that takes an exact value of 5.0.

Once the friction velocity has been obtained by solving this transcendental equation, the values of 𝑘 and 𝜖 to be applied at the surface are given by expression (3.23): 𝑢𝑘=^2𝜏𝐶_𝜇𝑢𝜖=^3𝜏.𝜅𝑦 (3.23)

3.6.1. Wave Generation

Essentially, surface gravity waves generation is a boundary condition that introduces a time-dependent mass and momentum flux through one of the domain faces. The flux varies in time because the velocities and the area through which the flow enters the domain vary in time. Therefore, the surface gravity wave generation boundary condition does not only affect velocities but also requires to impose conditions on the free surface elevation, represented in the model by means of the VOF function.

In the present study where only tsunami waves are considered, solitary wave is the only wave theory used. Solitary waves are modelled by means of Boussinesq equations (obtained from [17]) that are presented as follows: 𝜂=𝐻sech34𝐻ℎ𝑋ℎ²inwhich𝑋=𝑥−𝐶𝑡𝐶=𝐻𝑔ℎ1+ℎ𝑢=𝜂𝑔ℎℎ𝜂1−+ℎ4ℎ3𝜂²1−3𝑧²2ℎ²𝑑²𝜂𝑑𝑋²𝑣=𝑔ℎ−𝑧ℎ𝜂1−2ℎ𝑑𝜂+ℎ𝑑𝑋3²𝑧1−²2ℎ²𝑑³𝜂𝑑𝑋³, (3.24)(3.25)(3.26)(3.27) where 𝜂 is the free surface elevation respect to mean level, 𝐻 is the wave height, ℎ is the water depth, 𝐶 is wale celerity, 𝑔 is the gravity acceleration, 𝑢 is the horizontal component of velocity, 𝑣 is the vertical component, and 𝑧 is the depth of the point in which the speed is being calculated with respect to the bottom.

The variation of the area through which the flow enters the domain is represented in the model by a variation of the VOF function values at the faces of the wave inflow boundary. Figure 1 shows in a simplified bidimensional case two different time instants of the wave generation process. The generated wave travels from left to right. This fact is simulated in the figure moving the generation boundary (red line) from right to left.

Figure 1 

VOF boundary condition.

At time 𝑡₀, every cell position is compared with the free surface elevation. Free surface elevation is obtained from (3.24). Those cell faces that lie completely below the free surface elevation are given a VOF value of 1.0 (blue dots). If the cells lie completely above the free surface level, VOF is set to 0.0 (yellow dots). Cell faces that are cut by the free surface receive a VOF value equal to the relation of the wet surface of the cell face and the total surface of that face (purple dots).

At a later time instant 𝑡₁, the free surface position has changed, and, therefore, VOF function values at the faces of the inflow boundary must be updated to reflect the current time conditions. It can be observed in Figure 1 that as the free surface elevation raises at the inflow boundary, there are more full cells (blue dots), because the free surface elevation raises over more cells and less void cells (yellow dots). The position of the partially filled cell (purple dot) has changed, because the intersection point between the free surface and the inflow boundary has changed.

Having set the VOF values at the boundary, velocity values are calculated according to (3.26) and (3.27) and applied to the boundary faces.

Once the VOF function values and the velocity values have been calculated, the initial problem has been transformed into the imposition of Dirichlet conditions on the VOF and the velocity variables, which is a direct task.

4. Results

The model IH3-VOF has been used to study tsunami wave interaction with vertical impermeable and porous coastal structures. The model results have been previously validated using laboratory experiments as presented in [8, 10]. In this section, only numerical results are presented.

4.1. Numerical Simulations

The numerical simulations are carried out in a numerical wave basin (see Figure 2) 17.8m long, 8.6m wide, and 0.8m high. The structure is located 11.0m away from the wavemaker. The structure is 0.5m long, 4.0m wide, and 0.5m high. The dimensions of the structure are given considering the length is measured in the wave propagation direction, the longest direction of the flume, the width is measured in the horizontal orthogonal direction, and height in the vertical direction. The structure is located against the left wall, following the direction of wave propagation, of the wave basin.

Figure 2 

Bodies to define mesh geometry for the impervious structure case.

An orthogonal mesh is used for both the impervious and porous structure investigations. The mesh discretization in the longitudinal direction (𝑋 direction) is variable. At the wavemaker (in Figure 2 left end of the fuchsia body), Δ𝑥=4cm. From that point, Δ𝑥 is reduced until reaching a point 2.0m away from the structure (left limit of the cyan body), where Δ𝑥=2.0cm. The area around the structure (cyan body) is meshed with a Δ𝑥=2cm resolution until a point that is 2.0m downstream the structure (right limit of the cyan body). From that point, it grows linearly again until it reaches Δ𝑥=4cm at the end of the flume (right end of the green body). The most interesting processes occur around the structure and, therefore, in order to study them into more details, the finer discretization is used in this area. In the spanwise (𝑌 direction) direction, the discretization is homogeneous all along the flume and is Δ𝑦=2cm. In the vertical direction, very important to model the free surface processes correctly, Δ𝑧=0.75cm. Mesh dimensions have been calculated to have at least 6 cells to represent the wave height.

The mesh built this way contains about 12.000.000 elements, which takes 950MB of disk space. In order to run the simulations there is a need to make use of the parallel capabilities of IH-3VOF. Depending on the supercomputer used to run the simulations, the optimal number of cells per processing task ranges from 15.000cells/core to 50.000cells/core. The wave basin simulations have been run making use of 200 processor cores. The total simulated time for all the different cases has been 20s. The simulations take around 72h to complete and generate an output data of 500GB.

Numerical simulations have been carried out to reproduce and compare the results with the laboratory experiments. A 𝑘-𝜖 turbulence model is used for all the numerical simulations. The pore-based Reynolds number expected in the simulations is high, very similar to the narrow flume case. Therefore, the flow inside the porous medium is expected to be fully turbulent. For this reason and also because of the good agreement found in the narrow flume experiments, the porous medium parameters are taken equal to the crused rock values obtained in the previous chapter (𝛼=10.000, 𝛽=3.0 and 𝑐=0.34). The no-slip boundary condition is used for all the solid boundaries at which a wall function for the turbulence model is also imposed. An atmospheric boundary condition (𝑝=0) is used for the top face of the mesh. Water properties used are 𝜌=1000kg/m³ and 𝜇=1⋅10⁻³m²/s. Air properties used are 𝜌=1.20kg/m³ and 𝜇=1⋅10⁻⁵m²/s.

4.2. Flow Description

The left-hand column of Figures 3 and 4 shows different snapshots of the solution obtained with the numerical model, where the red color indicates higher elevation than the blue one. A solitary wave of height 𝐻=9cm has been generated in a water depth of ℎ=40cm. At time 𝑡=6.0s, the wave approaches the structure, but no disturbance due to reflection is observed so far. At the next time step, 𝑡=7.0s, the wave that impacts against the structure induces a free surface elevation raise over the seaward face of the structure, while the rest of the wave continues to propagate. Maximum height reached by the wave against the structure face is located at the wall.

(a)

(b)

Figure 3 

Snapshots at different times of the solitary wave interacting with the impermeable structure (a) and the porous one (b).

(a)

(b)

Figure 4 

Snapshots at different times of the solitary wave interacting with the impermeable structure (a) and the porous one (b).

At time 𝑡=8.0s, it can be clearly seen that two waves have been generated in the interaction. The first one is formed by the original wave (near the left wall) and the diffracted wave that propagates to the shadow area in the leeward side of the structure. The second one is the reflected wave that propagates back to the wavemaker. It can be observed that the reflected wave height is higher near the wall where the confinement conditions configure a bidimensional-like reflection. The three-dimensional structure of the reflected wave generates also the diffraction of these waves.

Time 𝑡=9.0s panel shows the moment in which the diffracted wave generated at the structure head corner reaches the side wall of the basin. It can be seen that the red intensity of both waves is reduced due to wave diffraction effects. It is also very important to note the wave radiation produced by the structure head. Time 𝑡=10.0s illustrates the propagation of the waves radiated from the structure head. Moreover, it shows the impact of the wave against the basin end wall. At that point, it can be observed, the high three-dimensional structure of the wave. Indeed the spanwise propagating component of the wave generates a high point at the right end corner of the basin on time 𝑡=11.0s. At that time instant, it can be observed that wave diffraction has generated a spreading of the original reflected wave that makes it occupy almost all the basin width.

The last two time instants 𝑡=13.0s and 𝑡=15.0s show the propagation of the waves reflected back from the end wall and the wavemaker. It can be seen that free surface patterns are very complicated and that waves are propagating in every direction inside the basin. At this point, it is very difficult to separate the effect of the different processes taking place at the basin.

The right hand column of Figures 3 and 4 presents different snapshots showing the interaction of the solitary wave of height 𝐻=9cm and water depth ℎ=40cm, with the porous structure. There are very small differences between this case and the impermeable structure. At time 𝑡=6.0s, the solitary wave reaches the position of the structure. At time 𝑡=7.0s, the free surface elevation rises at the seaward face of the structure. The elevation presents a smoother free surface gradient than the one found for the impermeable structure.

The reflected wave at time 𝑡=8.0s has a smaller wave height in the porous structure case than in the solid structure one. The diffracted wave crest is slightly larger for the solid structure than for the porous one. This can be observed at time 𝑡=9.0s at which the diffracted wave reaches the right wall of the wave basin. The most important difference at the diffracted wave appears in the last time instants, 𝑡=13.0s and 𝑡=15.0s, where the dissipation produced by the porous medium results in a smaller free surface elevation for the porous body experiments than for the solid structure case. Wave radiation can be observed at different time instants of the experiment (𝑡=9.0s, 𝑡=10.0s, and 𝑡=11.0s).

5. Discussion

In order to provide new insights about the interaction of a tsunami wave with a vertical structure, new plots are presented here. Both impermeable and porous vertical walls are compared.

Figures 5 to 7 present free surface elevation isolines comparing the different wave patterns developed in the interaction of a solitary wave with a porous or an impermeable structure. The upper part of the panels shows the impermeable structure isolines, and the lower part shows the porous ones. Free surface elevation is normalized by the incident wave height to simplify the interpretation and comparison of the results. Isolines are drawn at 0.25 intervals.

(a)

(b)

Figure 5 

Free surface elevation isolines for the impermeable and porous structure under the action of the solitary wave of height 𝐻=9cm and water depth ℎ=40cm, part I.

Time 𝑡=7.0s panel (see Figure 5) shows the moment in which the incident wave splits into a reflected and an incident wave. The free surface elevation pattern around the structure is quite different in both cases. In the solid body case, the free surface elevation in the seaward face of the structure reaches twice the incident wave height. In addition, the wave height in front of the structure head is slightly higher than the incident wave height.

In contrast, in the porous structure experiment, the wave height in front of the structure is slightly smaller than the incident wave height. The free surface elevation at the seaward face of the structure only reaches one and a half times the incident wave height. Free surface gradients are smoother in the porous structure case except inside of the porous structure, where the gradient is very high (15cm of elevation difference in 50cm). This free surface elevation gradient is responsible for the flow within the porous structure.

At time 𝑡=7.5s (see Figure 6), the differences in the reflection and diffraction processes for both structures can be clearly seen. In the solid structure case, the free surface elevation gradient between the transmitted wave and the shadow area leeward of the structure is more important than in the porous medium case. This gradient is responsible for the trough that appears around the structure head in the impermeable case. The higher gradient also generates a more important diffraction as it can be observed at times 𝑡=8.0s and 𝑡=8.5s.

(a)

(b)

Figure 6 

Free surface elevation isolines for the impermeable and porous structure under the action of the solitary wave case of height 𝐻=9cm and water depth ℎ=40cm, part II.

(a)

(b)

Figure 7 

Free surface elevation isolines for the impermeable and porous structure under the action of the solitary wave case of height 𝐻=9cm and water depth ℎ=40cm, part III.

Wave reflection patterns are also different. While the reflected wave in the impermeable structure case presents a wave height equal to the incident wave height, in the porous medium case it is only slightly higher than half the incident wave height. It is quite important to note also the pronounced three-dimensional structure of the wave reflected in the impermeable structure. Free surface gradients in the spanwise direction are more important in this case. Within the porous medium, free surface gradients have been reduced and are similar to the gradients outside the porous medium.

At time 𝑡=8.0s (see Figure 7), the depression located around the impermeable structure head has been entrained by the reflected wave and has moved to the seaward face of the structure. This depression does not exist in the impermeable structure simulation, in which the free surface gradients within the porous medium have almost completely dissipated. The shape of the diffracted wave crest is slightly different in both cases. In the solid structure case, the wave front is orthogonal to the structure, while in the porous case, it is not. This effect is produced by the porous medium drag that forces the wave front to bend towards the porous medium (by means of wave refraction).

The difference in wave diffraction induced by the porous medium presence can be clearly seen at time 𝑡=8.5s (see Figure 6). In the impermeable structure case, the 0.25 isoline is sticked to the leeward side of the structure, while a depression is present in the seaward side. In the porous simulation, the 0.25 isoline is not in touch with the structure anymore. Free surface elevation, however, is lower in the leeward face of the porous structure, but no depression is found in the seaward face. This clearly shows that the differences in wave diffraction are due to the link between the leeward and seaward faces that is established by the porous structure.

Time 𝑡=9.0s shows how the depression in the seaward face of the structure grows as the reflected wave travels back to the wavemaker position. In the porous structure case, this “suction” deforms the diffraction tail towards the structure. At 𝑡=9.5s, it is shown how the initial highest free surface gradient present in the impermeable structure case produces a higher wave height at the right wall of the basin. It is also important to note how the wave height along the left wall of the wave basin has decreased from its initial value. Initially, the incident wave height extended for half of the basin width. However, at the last time step, there is only a small portion of the wave that keeps the incident wave height.

One of the main advantages of using the IH-3VOF model is that the NS equations resolve the vertical structure of the flow velocity. The importance of the structure analysis relies on the fact that the associated flow processes induced by a tsunami wave, such as potential bed erosion around the breakwater head, can be detected. As a preliminary step, an analysis of the wave induced flow around the breakwater head is carried out.

Numerical results are presented in Figures 8 and 9. In both, horizontal velocity module is represented using isocontours for the impermeable and the porous structure simulations. In Figures 8 and 9, horizontal velocity is presented at two levels, the upper part at 𝑧=0.02m and the lower one at 𝑧=0.2m, in order to study vertical variations of the flow. Only snapshots corresponding to 𝑡=7.0s and 7.5s are presented, because the vertical flow structure was more relevant for those time instants, which correspond to the wave crest passing around the breakwater head.

(a)

(b)

Figure 8 

Velocity module isolines for the impermeable and porous structure under the action of the solitary wave case of height 𝐻=9cm and water depth ℎ=40cm. Data in m/s. Upper figure slice at plane 𝑧=0.02m. Lower figure slice at plane 𝑧=0.2m, part I.

(a)

(b)

Figure 9 

Velocity module isolines for the impermeable and porous structure under the action of the solitary wave case of height 𝐻=9cm and water depth ℎ=40cm. Data in m/s. Upper figure slice at plane 𝑧=0.02m. Lower figure slice at plane 𝑧=0.2m, part II.

As shown in Figure 8, the tsunami wave induces an important increment of the velocities around the structure head at both breakwaters. Velocities are observed to be around 30% higher at the head than the ones predicted by the model at the wave crest. The increment is similar in both impermeable and porous cases. However, the velocity increment affects to a slightly smaller area in the case of the porous structure, as a result of the mass flow percolating the porous vertical wall. It is also important to note the vertical structure of the flow, as it can be identified analyzing differences between horizontal velocity near the bottom and at middepth (see Figure 8). The vertical structure of the flow can be clearly observed on the upper wall of the domain, where the velocity near the bottom is almost a 20% smaller than at middepth.

At the breakwater head, as a result of the disturbance introduced in the flow by the structure, the numerical results at both levels show similar flow patterns, as shown in the lower and upper figure panels in Figure 9, revealing a quasi-uniform flow in depth. This increase in near-bottom velocity around the breakwater head is responsible for the initiation of erosion processes that reduce structural stability.

The wave reflection is also visible in the velocity analysis. Due to the wave mass and energy transmission through the porous wall, the model predicts lower velocity values at the seaward side of the porous breakwater. At time 𝑡=7.5s (see Figure 9), differences can be found for the flow horizontal velocity in front of the breakwater. The reflected wave at the wall is lower in the porous wall, due to transmission, inducing lower velocities.

6. Conclusions

In this work, a numerical analysis of the interaction of a tsunami wave with a vertical structure has been performed. Both an impermeable and porous structure has been considered in order to analyse the difference in the breakwater typology on the wave induced hydrodynamics. The use of a Navier-Stokes model, called IH-3VOF [9], has been considered in order to describe the nonlinear interaction of the tsunami wave with a vertical breakwater. The model is also able to calculate the three-dimensional wave induced patterns in the vicinity of the breakwater considering also the turbulent magnitudes. Moreover, the ability of the model to evaluate wave flow within the porous media by means of the volume-averaged Reynolds-averaged Navier-Stokes equations (VARANS) is used here, to describe different tsunami wave transformation processes in the near field. The model, previously validated by [8, 9], has been described in detail, paying special attention to the numerical scheme used in both time and spatial discretization.

Numerical experiments have been presented in terms of free surface evolution and wave induced velocities. First, free surface evolution has been investigated showing larger differences between the two studied typologies. Reflected wave at the breakwater has been detected to be lower in the porous breakwater, due to the wave transmission and the energy dissipation through the porous media. Diffracted wave has been also analyzed, identifying clearly three-dimensional wave patterns. However, differences between the two studied typologies are not as large as the observed for the reflected wave. It has been also noted that high three-dimensional wave features have appeared at the breakwater head, with lower run-up than in the breakwater trunk.

Similar conclusions can be drawn from the analysis of the velocity field. The wave induced pattern around the breakwater is highly affected by the presence of the porous media. At the breakwater head, similar flow characteristic has been detected. The flow clearly increases at the breakwater head at both typologies due to the existence of a flow separation and a fully three-dimensional wave pattern. However, due to the flow percolation through the porous breakwater, lower velocities are observed for the porous vertical wall. The flow has been detected to be close to uniform in depth at the breakwater head. This difference in the vertical structure of the flow, departing from the theoretical profile of the wave, is of great importance when considering the stability of the foundation of the structure. The accurate determination of the erosion velocities allows a proper design of the structure's foot protection and therefore improves the overall stability of the structure.

A more evident three-dimensional flow structure has been observed in the flow far from the head. The influence of the porous structure is detected as a decrement in the velocity magnitude. Residual motions after the tsunami wave have also been observed at both breakwaters, being lower in magnitudes the ones for the porous structure.

The IH-3VOF model presented here is revealed as a very promising tool to study the interaction of a tsunami wave with coastal structures. In this work, a simplified configuration of a vertical breakwater has been performed. Results shown here present IH-3VOF model as a promising tool to be used in the future with more complex configurations, such as traditional rubble-mound breakwaters.

Acknowledgments

M. del Jesus is indebted to the Spanish “Ministerio de Educación y Ciencia” (M. E. C) for the funding provided in the FPU Program (AP2007-02225). J. L. Lara is indebted to the M.E.C. for the funding provided in the “Ramon y Cajal” Program (RYC-2007-00690). The work is funded by Projects BIA2008-05462, BIA2008-06044, and BIA2011-26076 of the “Ministerio de Ciencia e Innovación” (Spain).