A THREE-DIMENSIONAL NAVIER-STOKES MODEL FOR WAVE AND STRUCTURE INTERACTION

This paper describes the capability of a new model, called IH-3VOF to simulate wave-structure interaction problems using a three-dimensional approach. The model is able to deal with physical processes associated with wave interaction with porous structures. The model considers the VARANS equations, a volume-averaged version of the traditional RANS (Reynolds Averaged Navier-Stokes) equations. Turbulence is modeled using a k- approach, not only at the clear fluid region (outside the porous media) but also inside the porous media. The model has been validated using laboratory data of free surface time evolution in a fish tank containing a porous dam. Numerical simulations were calibrated by adjusting the porous flow empirical coefficients for two granular material characteristics. Sensitivity analysis of porous parameters has also been performed. The model is proven to reproduce with a high degree of agreement the free surface evolution during the seeping process. Simulations of a three- dimensional porous dam breaking problem has been studied, showing the excellent performance of the model in reproducing fluid patterns around a porous structure. The model is powerful tool to examine the near-field flow characteristics around porous structures in three dimensional flow conditions.


INTRODUCTION
Design conditions for coastal structures should include a correct assessment of the structure functional response.The most relevant hydraulic processes to be considered in wave-structure interaction encompass wave reflection, wave dissipation, wave transmission resulting from wave overtopping and wave penetration through the porous structures, wave diffraction, run-up and wave breaking, involving several time and length scales, such as gravity waves (~100 m, ~10 s) or turbulent motions (<1 s, ~10 -3 m).
Traditionally wave-structure interaction has been studied through physical tests (two-and threedimensional), especially small-scale model tests.Empirical formulations arisen from physical modeling present several restrictions and a narrow range of applicability.Many other issues related with scale factors or processes such as porous media flow, wave impacts or viscous effects are not correctly represented in the experiments.
A great effort has been made over the last decade in the numerical modeling of wave interaction with coastal structures to overcome these limitations.Several approaches have been followed to study coastal structure induced hydrodynamics.Among other existing approaches, Nonlinear Shallow Water (NSW), Boussinesq-type or Navier-Stokes equations models have traditionally been used.
Although good results in terms of averaged magnitudes have been obtained using NSW equation for example Kobayashi et al. (1989) or more recently Kobayashi et al. (2010), vertical velocity structure can not be resolved using this approach.Moreover, the energy transfer to higher frequencies occurring before wave breaking cannot be reproduced accurately due to the lack of dispersion.Besides porous media flow is not directly solved, and it has to be included in the momentum equation through an empirical friction term.
Boussinesq-type models are able to include frequency dispersion, a depth-dependent velocity profile, and they can be applied to both breaking and non-breaking wave conditions.A great effort has been made in order to relax the original equations by deriving the extended Boussinesq equations (Kirby, 2003).However, this type of models requires setting both the triggering wave breaking mechanism and the subsequent wave energy dissipation due to wave breaking.These models fail to reproduce the strong nonlinear shoaling prior to wave breaking and the free-surface and velocity higher order statistics which are thought to be relevant for structure stability.Wave interaction with porous structures in the absence of breaking using Boussinesq equations has been performed by Lynett et al. (2000).
Navier-Stokes equations models have been incorporated in the last decade to the study of wavestrucure interaction problems.The number of simplifying assumptions assumed in the equations is lower than in other approaches.These models are able to calculate flows in complex geometries and provide very refined information on the velocity, pressure and turbulence field.
Models based on a two dimensional eulerian Navier-Stokes set of equations (Losada et al., 2008;Lara et al., 2008;Guanche et al., 2009) have proven to be powerful to address wave-induced processes.Moreover, wave flow within the porous media can be simulated.Stability (i.e.: waveinduced loads) or functionality (i.e.: wave reflection or overtopping) have been reproduced numerically with a high degree of accuracy, introducing new models to be used as a complementary tool in the design process.However, the range of applicability of these results is limited to normal wave incidence.
Three dimensional wave-induced processes such as wave radiation and diffraction around breakwater trunks and heads or oblique wave incidence can only be studied based on a three dimensional approach.SPH (Smoothed Particle Hydrodynamics) has recently been applied to coastal engineering.This approach solves the flow from a Lagrangian point of view, calculating the kinematics of each particle and its interaction with neighbour particles.The Lagrangian nature of SPH makes it well suited to simulate free surface flows with rapid changes of the flow field.However, SPH models are very expensive from a computational point of view and they cannot be applied to solve the large domains that wave-structure interaction requires.Recently Shao (2010) it has been presented a two dimensional model which takes into account wave interaction with a porous flow structure.
Several approaches based on the use of Eulerian three dimensional Navier-Stokes sets of equations can be found in the literature, for example Li et al. (2004), Liu et al., (2009), Wang et al., (2009) or Christensen (2009).However, none of them are able to solve porous media flow therefore coastal structures are treated as impermeable.In this work, a new three dimensional Navier-Stokes model, called IH-3VOF is presented.A new set of equations has been derived to include porous media flow in the equations.
The paper is organized as follows.In the next section, the numerical model is presented including the mathematical description and the numerical implementation.The porous flow calibration-validation is then presented for a classic dam-break problem.Finally, the hydrodynamics induced by a porous dam-break problem is described.A series of conclusions close the paper.

NUMERICAL MODEL DESCRIPTION
The present study uses the IH-3VOF model, whose main features are described in this section.IH-3VOF is based on TRUCHAS a model developed at Los Alamos National Laboratory, which solves three dimensional Navier--Stokes equations for a multi-phase transient flow.Free surface is tracked with the volume-of-fluid technique (VOF) presented by Rider et al. (1998).IH-3VOF extends the TRUCHAS model applicability by solving flow in porous media with the VARANS equations and a kε turbulence model.Moreover surface gravity waves can be generated.All the features make IH-3VOF a tool especially well suited to deal with coastal engineering issues.In this section mathematical formulations are presented first and numerical implementation details are given next.

Mathematical formulation
The IH-3VOF model solves the three dimensional Reynolds Averaged Navier-Stokes (RANS) equations, based on the decomposition of the instantaneous velocity and pressure fields into mean and turbulent components.Reynolds stresses are closed with a k-ε turbulence turbulence model.
The direct resolution of the flow field inside the porous media is not affordable, given the complex structure of porous materials.In practice, rubble mound coastal structures are made of individual stones or armor units with a random geometry.Moreover, solving the flow characteristics inside the gaps implies a very fine grid increasing the computational cost.To overcome these limitations, flow within the porous media is then volume averaged in a control volume larger than the porous structure but smaller than the characteristic length scale of the flow.The following spatial average operator is introduced and applied to the RANS equation: where the subscript f means intrinsic magnitude.Intrinsic magnitudes can be related with averaged magnitudes (also called Darcy's magnitudes) with: where n is the porosity and is defined as: Applying this operator implies that the Reynolds averaged magnitudes are split into a mean spatial averaged quantity and a spatial fluctuation part as follows: It is important to note that this operation transforms a porous media made of a solid skeleton and gaps into a continuous and homogeneous medium, characterized by material properties such as porosity or nominal diameter.Once this operator is applied small scale spatial fluctuations disappear from the solution, and then the small details of the porous channel geometry are no longer needed in order to obtain solutions.
Introducing now the decomposition given in (2) into the RANS equations, the volume-averaged Reynolds-averaged Navier--Stokes equations (VARANS) are obtained, and also the mass conservation equation for the volume-averaged ensemble-averaged variables: It can be readily seen that mass conservation preserves its expression but changing the instantaneous variable with the ensemble-averaged velocity or with the volume-averaged ensembleaveraged velocity.However, at each averaging operation applied to the momentum conservation equation, new terms appear that need to be modeled in order to obtain accurate solutions (see second line in equation ( 3)).
The first of these terms is the volume-averaging Reynolds stress-tensor while the last two terms correspond with the effects that the porous medium causes to the flow.These two terms need closure equations to model the behavior of the processes by them represented.This closure model is known as the Forchheimer equation and can be expressed as, where c A is the added mass coefficient and α and β two empirical coefficients associated with the linear and nonlinear drag force respectively.The third term of the right-hand side of equation accounts for the inertial effect.The precise descriptions of the c A , α and β coefficients are still not fully understood.They depend a priori on the pore Reynolds number and flow directions.
In order to close the VARANS equations ( 3) and (4), a two equation turbulence model may be used.IH-3VOF makes use of a k-ε model.In order to also consider the volume averaging within the porous media, the k-ε equations need to be transformed, yielding the form presented in the following equations: It can be shown that during volume averaging of the turbulence model new terms appear representing the effect of spatial averaging (last terms in equations ( 6) and ( 7).These terms model the turbulence production and dissipation that occurs inside the averaging volume which is filtered out in the averaging process.Modelling proposed by Nakayama and Kuwahara (1999) are considered as follows: ( ) ( ) IH-3VOF uses the VOF (Volume of Fluid) technique to track the free surface evolution applying Rider et al. (1998).The algorithm is based upon a well-defined, second-order geometric solution of a volume evolution equation.The method utilizes local discrete material volume and velocity data to track interfaces of arbitrarily complex topology.A linearity preserving, piecewise linear interface geometry approximation ensures that generated solutions retain second-order spatial accuracy.This point is crucial in terms of an accurate reproduction of the free surface evolution in complex geometries such as wave breaking.

Numerical implementation
Equations are discretized in the IH-3VOF model following the finite volume technique.A collocated grid is considered in order to define magnitudes.Both structured and unstructured meshes can be considered.To suppress spatial oscillations in the pressure and velocity, the well known Rhie-Chow method is applied.Parallel computation techniques are used to speed up simulations.Chaco software (Sandia National Laboratory) is used to perform domain subdivision.Chaco does not account for the physical characteristics of the various regions of the mesh, only the connectivity of the mesh.
A two step fractional method is considered to solve velocity and pressure field.The resulting Poisson equation is solved with the FGMRES (generalized minimal residual) method.

FLOW IN POROUS MEDIA
The use of the numerical model as a predictive tool for wave-porous structure interaction requires the previous calibration of the numerical parameters governing the hydrodynamics of the system.Although several works solving different flow problems can be found in the literature, calibration is needed for the new set of equations presented here.Experiments carried out by Lin (1998)

Experiments description
In order to validate the model, experiments carried out by Lin (1998) are used.The main reason for choosing these experiments among other existing in the literature is because they provide free surface evolution in both porous media and clear fluid (outside porous media) regions.The fluid interaction with the porous structure and its evolution in time is achieved in detail.Lin (1998) recorded flow passing through a porous dam in a fish tank, using two types of granular material: crushed stones and glass beads.In this work, only crushed stones have been used.Nominal diameter (D 50 ) and porosity (n) were 1.59 cm and 0.49.The fish tank was 89.2 cm long, 44 cm wide and 58 cm high.The porous dam, placed at the middle of the tank, was 29 cm long, 44 cm wide (it occupies the whole fish tank width) and high enough to ensure that water never reached its upper edge.Water was initially confined in a reservoir between one side of the tank and a gate 2 cm away from the porous dam (see figure 1).Although, three water elevations of the initial reservoir were considered for both granular material characteristics by Lin (1998), 35 cm, 25 cm and 15 cm, only 25 will be used here.At the beginning of the experiments, the gate was opened and the water was released interacting with the porous material and percolating into the dry area.Water motion though the porous dam and the clear fluid region was recorded by a camera.Further details about the experiments can be found in Lin (1998).

Computational domain
The numerical wave tank was designed as similar as possible to the experimental set-up.Different discretizations were used in order to investigate the cell size dependence on the numerical solutions and also to evaluate the computational cost.Two uniform grids were considered.The cell size was 1x1x0.5 cm (Δx/Δy/Δz) and 0.5x0.5x0.25 cm (Δx/Δy/Δz) for the coarser and the finer grid respectively.
Figure 2 shows the influence of the cell size in the numerical computations for the medium water elevation at the initial reservoir for a crushed stones dam.Laboratory data are presented in circles.Coarse and fine grid calculations are presented with a dashed black line and a solid gray line respectively.As can be observed in the figures, both numerical simulations show almost the same results.When water rushes to the dam some reflected wave at the left side part of the porous dam are created, generating a small amplitude bore (see t=0.4 s panel) which is better defined by the finer grid size.A similar behavior is observed at t=0.8 s stage when the reflected wave reaches the fish tank wall.No significant differences are observed when water is seeped out of the dam at its right side.
According to the presented results and also because of the larger computational cost observed for the finer grid, the coarser grid was used for this study.

Numerical model calibration
As previously explained, the flow field in the porous media is described in the IH-3VOF model by volume-averaging the RANS equations.The process of volume-averaging leads to the introduction in the momentum equation of two unclosed terms which are modeled collectively using the Forchheimer's relationship as presented by Hsu et al. (2002).
The expression used in the model for the linear and nonlinear drag forces is given by, ( ) ( ) where α and β are the only two empirical parameters used in the calibration of the numerical model, n is the porosity, ν the kinematic viscosity and f i u stands for velocity.An additional parameter, C A affecting the inertia term is included in the formulation to take into account the added mass.
The precise descriptions of α and β coefficients are still not fully understood for oscillatory flows.
They depend on several parameters such as the Reynolds number, the shape of the stones, the grade of the porous material, the permeability and of course the flow characteristics.van Gent (1995) and Burcharth et al. (1995) have suggested different values of porous drag parameters.However, when comparing different formulations attention has to be paid to the fact that discrepancies may exist in the expressions of the terms and that parameter values change depending on the flow conditions.Based on previous works focused on a two dimensional RANS modeling, Garcia et al., (2004), Lara et al., (2006a), Lara et al., (2006b), Losada et al. (2008) or Lara et al., (2008), we have experienced that under oscillatory flow and waves propagating over slopes or breaking, values existing in the literature may not be valid since the experimental conditions for obtaining those formulae did not consider these effects.Since there is not a predictive methodology to determine α and β coefficients in advance, calibration has to be performed.Model calibration is carried out based on measured data from one single test for each porous material (crushed stones and glass beads), comparing computed and measured time evolution of free surface displacement at the whole fish tank including the porous dam.The middle initial water elevation case (25 cm) was considered for both granular material characteristics.

Porous flow parameters Tested values
Linear friction parameter (α) 1000, 5000, 10000 Non-linear friction parameter (β) 0.5, 3, 6 Tested values for the porous flow parameters are presented in table 1. Linear drag parameter α influence is first studied keeping the non-linear drag parameter (β) constant and equal to 6. Figure 2 shows the comparisons between the laboratory data and three numerical simulations performed varying the α parameter for flow passing through the crushed stones.Laboratory data are presented in circles.
The linear drag parameter α is presented by lines.Dashed lines correspond to α=1000, solid lines to α=5000 and dotted lines to α=10000.Numerical simulations are represented at the middle of the fish tank.
The three numerical simulations presented show a similar behavior, especially from t=0.0 s to t= 0.6 s when the water reservoir is released.Water motion within the porous dam is not very well predicted.In our opinion these discrepancies are due to the fact that water is released in a different way.In the numerical modeling, water is at rest at the initial stage and an instantaneous opening of the gate is considered.However in the laboratory experiments, a manually opened gate is move upwards with an estimated time of 0.1 s, as it is reported by Lin (1998).The water is then first accelerated at the bottom yielding a different flow behavior at early stages.Water starts to flow from the bottom while the gate is opened.
For later stages, flow is very well predicted by the model especially for α = 5000 simulation not only within the porous dam but also on both sides.Linear drag influence is noticed when the water on the right side of the dam is calm down and the water on the left side continues to seep into the clear fluid region (see figure 3).A larger α parameter (10000) exhibits a larger drag as expected, delaying water level drop at the right side of the dam.Flow percolates faster for the smaller α parameter (1000).Non linear drag parameter (β) influence is studied next.The linear drag parameter (α) is kept constant and equal to 5000. Figure 4 shows the comparisons between the laboratory data and the numerical simulations for flow passing through the crushed stones.Laboratory data are presented in circles.Non-linear drag parameter (β) is presented by lines.Dotted lines correspond to β = 0.5, dashed lines to β= 3 and solid lines to β= 6.In this case, larger discrepancies among the three performed simulations are observed.Poor correlations can be found for β=0.5 and 3 for time stages larger than 0.6 s.Early time stages are not considered in this analysis because of the previously commented reason.Both simulations exhibit a lower drag than expected from the experiments.Water is seeping faster into the porous dam at the right side of the tank.Water oscillations at both sides of the dam also appear different from the experiments for those parameter values.However, β=6 simulation shows the best agreement.Larger β values (not shown here) yield larger drag and consequently a worst agreement.
As a conclusion, the computed and measured free surface time evolution agree very well for the studied case.Porous flow parameters have been found to be different as expected to other example in literature because of the different flow characteristics.

THREE-DIMENSIONAL POROUS DAM BREAKING PROBLEM
As a preliminary step to simulate wave interaction problems with porous structures, a three dimensional porous dam breaking problem is presented here.A classical dam breaking problem was modified considering a porous media in order to induce both a three dimensional flow structure and fluid percolation through the dam.A 1.2 m long, 0.6 m wide and 0.6 m high domain was considered.The porous dam was 0.3 m long, 0.3 m wide and 0.6 m high, and it was placed 0.3 m from one side of the domain.Nominal diameter (D 50 ) and porosity (n) were 1.59 cm and 0.49.Same values for the linear and nonlinear flow parameter obtained in the validation were used in the simulation.Water was initially defined as a 0.30 m long, 0.6 wide and 0.4 m high box and was released to freely move according to gravity force.The rest of the domain was filled with water 2.5 cm high.The whole numerical flume was discretized by means of a regular orthogonal mesh with a 1x1x0.5 cm (Δx/Δy/Δz) cell size.
Figures 5 and 6 shows several snapshots of how the fluid is moving inside the domain.Each presented panel contains a couple of plots corresponding to two side views of the same time instant, which is displayed on top of each figure.Moreover, the porous dam is plotted in a gray transparent color in order to visualize the free surface evolution inside the dam.Initial conditions are presented in the upper left panel in figure 5.
Figure 4 shows the initial stages of fluid motion until the water reaches the opposite boundary.From the initial fluid motion, two different fluid patterns are observed.First, water in front of the porous dam percolated into the porous media.Second, the other portion of fluid was accelerated at the bottom developing a bore while propagating, as can be seen in the snapshots corresponding to instants between 0.12 s and 0.35 s.
Water level was observed to decrease faster in the region which was not protected by the porous dam due to the rapid movement of the fluid.This feature created a water level gradient which affects the percolation process.This was reduced because part of the water confined between the dam and the boundary start to flow laterally.Small percolation was also observed while the bore passed by side of the dam.Most of the fluid was moving along the gap existing between the dam and the side wall.
Once the fluid reached the end of the dam, part of the water moved towards the left filling the space between the dam and the rear end boundary.Stages between t=0.4 s and 0.5 s show it.Although three dimensional flow patterns were observed from the initial stage due to the problem geometry, a fully three dimensional free surface is clearly identified in the whole domain at this point.Due to the rapid decrease of the water level in front of the dam, not as much water percolated into the dam, as can be seen in stages between t=0.25 s and 0.5 s.
Fluid first reached the rear wall then the side wall behind the dam, generating a violent splash.The water then moved upwards, developing a three dimensional pattern, as can be seen in the lower snapshots in figure 5.At that time, the water level inside the dam decreased compared to the initial stages because of the water level gradient between the inner and the outer part of the dam.The forcing mechanism induced a flow discharge away from the dam.
Figure 6 shows the flow behavior once the water reaches the back side walls.Water piled up in the rear region, showing a fully three dimensional structure.Due to the lateral movements followed by the water occupying the space behind the dam, more water was accumulated at the region.Water was reflected at the wall and moved backwards.A fraction of the water was seeped into the dam (see t=1.25 s to t=1.5 s).Due to the friction created at the porous media, velocity decreased at the interface and water rose along it (see t=1.25 s and t=1.37 s).A backwards broken wave was developed generating a curly free surface due to the continuous impacts of the splashing drops.
Water continued moving towards the front wall with a fully three-dimensional pattern not only along the gap between the dam and the side wall but also through the porous media, as can be identify from t = 1.5 s to t=2 s.
At the final stage shown here (t=2 s), water continued oscillating and three-dimensional patterns were still clearly seen along the whole domain.

CONCLUSIONS
In this paper, a three dimensional Navier-Stokes model, called IH-3VOF is presented.The model is able to deal with physical processes associated with wave-interaction with porous structures.A new set of porous flow equations are presented.A k-ε model was also considered in the equations, not only for the clear fluid region but also for the porous media flow.
Laboratory experiments of dam-break flow were used in order to calibrate and validate the model.Model calibration is carried out based on measured data, comparing computed and measured time evolution of free surface displacement at fish tank including the porous dam.Numerical simulations were calibrated by adjusting the porous flow empirical coefficients for a crushed stone material.The model is proven to reproduce with a high degree of agreement the free surface evolution during the seeping process.
In order to prove the ability of the model in dealing with three-dimensional flow induced flows around porous structures a porous dam breaking problem was studied.Three dimensional flow partners were observed form the initial stages.Flow percolation inside the porous dam was identified during the simulation.
IH-3VOF is proven to be a powerful tool in examining the near-field flow characteristics around porous structures.Numerical information can be used to study physical magnitudes difficult or impossible to measure in the laboratory, such as particle velocity above or around the structures in wave breaking conditions or flow inside the porous media.
The numerical modelling based on RANS or Large Eddy Simulation (LES) is a breakthrough in our future way of approaching coastal structures design and it is currently in a relative early but solid stage.Additional work is still to be done in order to carry out a reliable use of these models for engineering design.
are presented first.Numerical set-up is discussed next.Calibration and validation work results are shown in the following sections.

Figure 1 .
Figure 1.Initial stage for the porous dam breaking.