CFD-CSD NUMERICAL MODELLING OF WAVE-INDUCED PRESSURES IN OPEN-PORED PBA-REVETMENTS

The highly porous Polyurethane Bonded Aggregates (PBA) revetments represent a novel ecologically friendly solution for the protection of shorelines and vulnerable coastal areas against erosion. Advantages of the open-pored PBA-revetments over conventional smooth impermeable revetments are among others, the reduction of: wave runup/run-down, wave reflection and wave-induced loads on the sand core beneath the revetment. However, the hydrogeotechnical processes involved in the interaction of waves with such PBA-revetments and their foundation are still not sufficiently understood. Therefore, a new 3-dimensional one-way coupled CFD-CSD model system "wavePoreGeoFoam" was developed at the Leichtweiß-Institute (LWI) within the OpenFOAM framework for the analysis of the response of open-pored PBA-revetments due to wave-induced loads. In this way, this paper firstly describes the new CFD-CSD model system. Second, validation of the model "wavePoreGeoFoam" is shown considering large-scale laboratory tests performed in the Large Wave Flume (GWK) at the Coastal Research Center (FZK) in Hanover, Germany (Oumeraci et al. 2010). Third, the relevance of implementing the CFD-CSD coupling for modelling wave-induced pressures on and beneath PBA-revetments is discussed. Fourth, a sensitivity analysis related to the effect of the empirically defined parameters for the numerical model is described. Finally, recommendations and implications of the use of CFD-CSD model for further research will be addressed.


INTRODUCTION
Polyurethane Bonded Aggregates (PBA) revetments are highly porous structures that have emerged as a novel alternative to prevent shoreline erosion and to protect the slope of an embankment subject to wave attack (Fig. 1).The open pored PBA-revetments provide advantages over standard impermeable revetments by reducing the wave run-up/run-down as well as the magnitude of waveinduced loads (Oumeraci et al. 2010;Liebisch et al. 2012;Foyer 2013).The hydro-geotechnical processes and the wave-structure subsoil interaction are still not sufficiently understood despite the comprehensive studies with open pored PBA-revetments conducted in large scale facilities such as the Large Wave Flume (GWK) at the Coastal Research Center (FZK) in Hanover, Germany (Oumeraci et al. 2010).Furthermore, research focused on numerical modelling of these type of revetments (Foyer 2013) was conducted in order to extend the results from the GWK laboratory tests but limitations were found for the accurate and appropriate estimation of the pressures and pore pressures developed on the revetment and in the layers underneath.For this purpose, a new 3-dimensional CFD-CSD model system "wavePoreGeoFoam" was developed and validated with the large scale laboratory tests performed in the GWK (Oumeraci et al. 2010) so that, the knowledge regarding the hydro-geotechnical processes occurred on and beneath PBA-revetments due to wave-induced loads can be enhanced.The CFD and the CSD modules developed at the Leichtweiß-Insittute (LWI) within the free open-source OpenFOAM ® framework were weakly coupled (one-way coupling) for the accurate analysis of wave-induced pressures and pore pressures.
Therefore, this paper will firstly provide a description of the CFD-CSD model system system "wavePoreGeoFoam".Also, it is shown that the CFD-CSD model system is able to reproduce the GWK large-scale tests so that, selected results of the modelling and their comparison with the experimental data are described for the validation of the model.Moreover, the paper will discuss the high relevance of the CFD-CSD coupling for the simulation of PBA-revetments and for the description of the wave-structure-subsoil interaction.A brief sensitivity analysis of PBA-revetments performance due to changes in geotechnical properties (i.e.porosity, permeability) is also addressed and further research steps of the modelling of wave-induced pressures in PBA-revetments will be discussed.

Numerical model description
The CFD-CSD model system is a 3-dimensional model was developed at the Leichtweiß-Institute (LWI) within the free open-source OpenFOAM ® framework.The Volume-Averaged RANS equations were implemented into the CFD solver with inclusion of the VOF method for tracking the free surface elevation.The CFD solver is able to handle two-fluid phases (e.g.water and air) and to explicit define the location of several porous regions with different properties (i.e.d50, porosity).The previous allows the definition of the different layers included in a PBA-revetment (cover layer, filter layer, sand core).On the other hand, the fully dynamic Biot's equations coupled with Darcy's law are implemented in the CSD solver which allows the determination of soil stresses as well as pore pressures developed inside the porous regions.The use of this solver, initially developed for the numerical simulation of the soil foundation underneath caisson breakwaters (El Safti et al. 2013), has been extended for the numerical simulation of PBA-revetments.
The CFD and the CSD solvers were weakly coupled in the numerical model system (one-way coupling).For this purpose, the pressure measured on the PBA-revetment and estimated through the CFD solver is taken as input for the coupling with the CSD solver which calculates the induced pore pressures and soil stresses developed underneath PBA-revetments in the filter layer and in the embankment subsoil (Fig. 2).

Governing equations of the CFD hydrodynamic module
The VARANS equations were implemented in the CFD module considering the approach provided by Hsu et al. (2002).Therefore, the continuity and the momentum balance equations implemented in the model are those described in eqs.( 1) and (2) for unsteady, incompressible, viscous and immiscible fluid flow: where p is the pressure; g the gravitational acceleration; n the porosity; x the Cartesian coordinate axis; U the velocity and ρ the weighted average density of the fluids according to the volume fraction γ occupied by each phase in a cell.The effective dynamic viscosity µ eff is calculated as µ eff = µ t +µ with µ t -Velocities -Pressures -Water surface elevation

CFD-solver (hydrodynamic model)
Results: -Pore pressures -Soil streesses ONE-WAY COUPLING at porous-nonporous boundaries as the dynamic turbulence viscosity (estimated by the turbulence model) and µ as the weighted average dynamic viscosity.Furthermore, the extra term σκ∇γ in the momentum balance equation (already implemented in OpenFOAM ® but not originally considered in Hsu's approach) accounts for the surface tension between the two phases (σ=surface tension and κ is the curvature of the interface κ=∇ • /| | with = ∇γ).
The term [CT] in eq. ( 2) stands for the closure term which is modelled through the Darcy-Forchheimer equation described in eq. ( 3), where the Polubarinova-Kochina terms (PK-term) is included to account for the added mass effect.

[ ] (
) Coefficients a and b in eq. ( 3) could be defined considering different approaches (Sidiropoulou et al. 2006).Therefore, these coefficients are user-defined in "wavePoreGeoFoam".On the other hand, coefficient c A can be determined through the approach described in van Gent (1993) which considered stationary and oscillatory flow through coarse porous media.
In order to track the two-phase flow interface, the volume of fluid (VOF) method is used including the effect of porosity to consider that, e.g. if the porosity of a cell is 0.35 then only 35% of the cell volume is available to be filled by the fluids.The latest is a very important consideration especially for those porous regions located into a dry-wet zone where the exchange of fluids between the free stream and the porous region is a highly predominant process (i.e. the dry-wet zone of a revetment or a breakwater where wave run-up and run-down take part).

Governing equations of the CSD hydro-geotechnical module
The discretization and implementation of the governing equations in the CSD module are described in detail in El Safti and Oumeraci (2012), El Safti et al. (2012) and El Safti and Oumeraci (2013).Therefore, only a brief summary of the governing equations is provided below.
The mass conservation of the fluid is considered in the Biot's theory as described by eq. ( 4) (Zienkiewicz et al. 1999).The mass conservation takes into account the degree of saturation of the porous medium, and therefore, the equation is applicable also to partially saturated porous medium.
In eq. ( 4), ε v is the volumetric strain of the solid matrix described by tr(ε) which is the trace of the strain tensor ε; U is the average Darcy's velocity vector of the percolating fluid and Q=K f /n (with n as the porosity and [K f = (S w /K w +(1-S w )/p 0 ) -1 ] as the bulk modulus of the pore fluid).K w is the bulk modulus for pure water and S w is the degree of saturation rate (S w =V w /V v , with V w and V v as the volume of water and voids, respectively) and p 0 as the absolute zero pore pressure (≈10 5 Pa under atmospheric pressure).
The fully dynamic formulation for the solid-fluid mixture is described in Zienkiewicz et al. (1999) and is written as in eq. ( 5): where σ is the total stress tensor, ℓ is the displacement vector, b is the body force per unit mass tensor (i.e.gravity), n is the porosity and ρ f and ρ s are the densities of the fluid and the solid particles, respectively.
The total stress can be correlated to displacements by the use of the 'constitutive models' and thus, eq. ( 4) and eq.( 5) can be solved simultaneously.Moreover, the calculation of the pore pressure (u) and the effective stresses (σ') is performed considering Terzaghi's principle (eq.( 6)).The effective stresses are calculated based on the relationship between strains and stresses: dσ'=E:dε' and the straindisplacement relationship based on the assumption of small strains as ε=0.5(∇ℓ +(∇ ℓ) T ) where T is the transposition operator.

'
u I σ = σ − , with I as the identity matrix (6) Therefore, the momentum balance equation of the fluid phase considering the solid phase as reference can be described as in eq.(( 7)), where the sink term R is defined as R=U ρ f g /K and represents the viscous drag force according to Darcy's seepage law, with K as the isotropic hydraulic conductivity in m/s.

Validation of the numerical model
Large-scale model tests (Oumeraci et al. 2010) were used for the validation of the CFD-CSD model system.Selected results of the validation are described and further details are provided in Alcérreca Huerta (2014).Thus, two regular wave tests from the GWK study with PBA-revetments (Oumeraci et al. 2010) were selected for the validation of "wavePoreGeoFoam" considering the following wave conditions: • Test case 1 (non-impact wave load condition): wave height H=0.6 m, wave period T=5.0 s, water depth h=3.7 m and surf similarity parameter ξ 0 =2.67.• Test case 2 (impact wave load condition): wave height H=1.0 m, wave period T=3.0 s, water depth h=3.6 m with a surf similarity parameter ξ 0 =1.25.
For the assessment of the Darcy-Forchheimer equation (closure term used by the VARANS equations), the approach provided by Engelund (1953) was used.The porosity n of the revetment, filter and sand core layers as well as the coefficients α and β in Engelund's approach (Engelund 1953) were defined prior the calibration of the CFD-CSD model system and the final values were found within recommendations commonly provided in the literature (e.g.Morris and Johnson 1967).For the particular calibration of the hydraulic permeability it was found that for the sand core of the PBArevetment the permeability is depth-dependent and thus it was set to K=2.9x10 -6 m/s for a depth 0<d<0.40mand K=4.5x10 -6 m/s for a depth d>0.40 m.A summary of wave conditions and the userdefined parameters for the CFD-CSD model system are described in Table 1.ⱡ Parameter for which calibration was necessary.
Additionally, the numerical model setup consisted on a setup similar than that of the GWK-tests which is exemplified in Fig. 3a (further details are described in Oumeraci et al. 2010 andAlcérreca Huerta andOumeraci 2013).Moreover, the location for the extraction of the pressure and porepressures was the same as that described for the pressure transducers used in the GWK-tests (Fig. 3c), so direct comparison can be made between the laboratory tests and the numerical simulations.
The comparison of pressure and pore pressure time series given by the numerical modelling and the GWK-tests is shown in Fig. 4. The comparison includes only the results for the column C3 of pressure transducers normal to the revetment slope (see Fig. 4).However, the systematic validation and comparison are described in Alcérreca Huerta and Oumeraci (2013), where similar results are shown for columns C1, C2 and C4 of pressure transducers (s.Fig. 3c).
The variation between numerical and laboratory results is herein reported through the relative error as in eq.( 8) with P exp and P num as the peak pressures from experimental and numerical results, respectively.For all pressure transducers at a certain depth placed in the GWK tests, an average relative error is herein reported.For the pressure time series on the revetment, the difference between the maximum peak pressures measured in GWK and those calculated by the numerical model provided a relative error RE≈7.9% for both impact and non-impact wave load conditions.Particularly for impact loads, the magnitude of the peak pressure is well captured, however, some difference in the shape of pressure time series is observed as expected since wave-wave interactions occur as consequence of the wave breaker impacting the revetment slope

Relative error [ ]
The pressure time series at PT12 and PT16 (0.20m and 1.00m deep in the sand embankment) have a similar performance for impact and non-impact loads in terms of shape and magnitude.Furthermore, the impact component of the wave loading is not observed for PT12 and PT16 in both numerical and GWK-tests which shows the capability of the CFD-CSD model system to capture the damping of this component through the layers in the sand embankment.The differences observed in the peak pressure for those locations beneath the PBA-revetment provide a RE lower than 25% (RE<25%).Despite that the RE is larger than that on the revetment, the differences are expected specially for layers deeper in the sand foundation where a variation of the permeability may highly affect the magnitude of the pore pressure measured by the pressure transducer in the location.Moreover, calibration of the empiricaldefined parameters for the numerical model was conducted and thus uncertainties on their definition may induce differences with experimental data.

Outcomes of the CFD-CSD coupling consideration
The CFD model was found to correctly estimate the pressures developed on top of PBArevetments and in the free stream region (s.Fig. 6a), as well as to account for the hydrodynamic processes (i.e.wave run-up/down, wave reflection).However, the CFD-CSD coupling is highly important for the correct estimation of the pore pressures, soil stresses and displacements underneath the PBA-revetment.In Fig. 6b, it is demonstrated that the use of the CFD model with the inclusion of the VARANS eqs.overestimates the magnitude of the pore pressures in the porous media by a huge percentage.The later was noticed to occur in CFD models with VARANS eqs.because this set of equations considers the solid skeleton and gaps in the porous region as an unique continuous and homogeneous medium.In fact, the velocity field (U) and the surface elevation (η) at the porous region can be correctly and accurately calculated by considering the aforementioned assumption, but it is not valid for the estimation of the pore pressures.In this sense, the CFD model provides a solution of the pressures but considering the solid skeleton-gaps as an unique media and thus, the CSD model should be used to differentiate between the soil effective stresses taken by the solid skeleton and the pore pressures developed in the water.On the other hand, the CFD-CSD coupling is not needed in the freestream region where no porous media is found (Fig 6a)  The pressures calculated by the VARANS equations (p) and the pore pressures described by the Biot's equations in a porous region (u) are related according to the Terzaghi's principle, where the sum of effective stresses and pore pressures leads to the total stress (σ) developed in the porous media (σ = σ' + u).With the VARANS equations the calculated pressure corresponds to the total stress σ ≈ p since no distinction is made between the pore pressures (u) and the effective stresses (σ', therefore, the CSDmodel is required to generate correct results, as shown in Fig. 6b.
The use of Biot's equations by the CSD module becomes more important for those cases where drainage in the soil is allowed.In the particular case of revetments, even that the degree of saturation of the soil below the water table is almost 100%, the water is not fully confined since the water table offers an open boundary which actually becomes a drainage boundary which results in the development of both effective stress and the pore pressures.Additionally, the coupling between VARANS and Biot's equations becomes relevant since the description of both the hydrodynamic and hydro-geotechnical processes are closely related in the description of the wave-structure interaction.

SENSITIVITY ANALYSIS: EFFECT OF EMPIRICAL-DEFINED PARAMETERS
For the GWK tests used for the validation of the numerical model, prior calibration of the empirically-defined parameters was conducted which led to the sensitivity analysis herein described.For this purpose, the results of the numerical modelling of non-impact loads was considered and the peak pressure on the revetment and the peak pore pressures in the sand core were used for the comparison between numerical and experimental results.The relative error according to eq. ( 8) is applied to quantify the errors at each pressure transducer location.
For the validation with the data from the GWK tests, an average of the relative error of 7.63% was obtained for the pressures on the revetment, and 17.43% for the pore pressures in the sand core (Test ID 14 in Table 2).The relative error for other test cases performed during the calibration became larger and a summary is provided in Table 2.The CFD solver based on VARANS equations can estimate accurately the pressure on the revetment and in the free stream region (non-porous region).

b) Pore pressures underneath the revetment
Coupling between the CFD and CSD solvers into a model system is required beneath the PBA-revetment for the appropriate calculation of the pore pressure (porous region).Considering the calibration test conditions described in Table 2, the effect of the empiricallydefined parameters (e.g.porosity (n), "Darcy-Forchheimer coefficients" (α f and β f ) and degree of saturation (S w ) on the results of numerical simulations was examined.The detailed sensitivity analysis is shown in Alcérreca Huerta and Oumeraci (2013), and a summary is briefly described considering the results shown in Table 2. a) Effects of the variation of porosity and coefficients α f and β f used in the definition of Darcy-Forchheimer equation: • These parameters are all related to the permeability of the soil and it was found that the gravel is less sensitive to a variation of these parameters than the sand.Moreover, for the numerical simulations the effective porosity of the soil should be considered, since only the connected pores contribute to the fluid motion.• A variation of β f (e.g.Test 14 vs. 19) induces larger effects on the flow motion and pressure development in the porous media than the variation of α f (e.g.Test 5 vs. 6).Similar results have been also described by Van Gent (1993).• A large increment of the parameter α f for the revetment cover layer (Test 3 and 6) slightly increases the pressure on the revetment, as expected since the revetment becomes less permeable.
b) Effects of the variation of the degree of saturation S w in the sand core: • An unsaturated sand core will lead to lower pore pressure than a fully saturated sand core.This is because for partially saturated soils, a part of the load is supported by the water and the rest by the soil skeleton.• Considering a degree of saturation S w =95% (Test 9) in the sand core and the porous cover layer, it was observed that the negative pore pressures (uplift) are reduced if compared with a case with S w =100% (e.g.Tests 8 & 9).This effect was also described by De Groot et al. (2006) and it is caused due to the air compressibility which allows the water to occupy a part of the air volume, thus resulting in a slight reduction of the pore pressure.
The performed sensitivity analysis has resulted in a rather qualitative assessment of the effects that the empirically-defined parameters may induce on the results of the numerical simulations.Among these parameters, α f and β f from Darcy-Forchheimer equation together with the porosity are the most relevant for the whole wave-structure-subsoil interaction and they are all related to the permeability of the porous media.Therefore, it can be concluded that the permeability is the key parameter to be determined, since it has a large influence on the magnitude of the pore pressures.However, further research should be conducted in order to describe the aforementioned effects in more detail.

CONCLUSIONS
The new developed weakly coupled CFD-CSD model system "wavePoreGeoFoam" was presented for the description of the hydrodynamic and hydro-geotechnical processes related to waves interacting with PBA-revetments and the soil underneath.Validation of the model was presented considering large-scale tests with PBA-revetments.For this purpose, the pressure time series at different locations were used and numerical and laboratory results were compared.Within the comparison, an averaged relative error lower than 8.0% was found for those pressures on the revetment and less than 25 % for locations beneath the revetment considering both impact and non-impact wave load conditions.
The importance of the CFD-CSD coupling was pointed out for the numerical modelling of waveinduced pore pressures beneath PBA-revetments.The use of the VARANS equations for PBArevetment can be used for the assessment of velocities, water surface elevation in both inside and outside the porous media as well as the pressure development on the PBA-revetment.However, its coupling with the Biot's equations is needed for the correct estimation of the pore pressures and soil stresses in the soil underneath the revetment.Further research may be required to enhance the results of the numerical simulations with a development of a fully coupled CFD-CSD model system for coastal engineering, especially for the solution of 3-dimensional problems for a more accurate calculation of the pore pressures and soil stresses.
During the calibration process for validation of the model, a sensitivity analysis was conducted.As a result, the permeability (related to the porosity and the Darcy-Forchheimer coefficients) has been recognized as the key parameter to be determined for the numerical modelling of a PBA-revetment and its soil foundation.Moreover, it was found that a finer porous media is more sensitive to the variation of the porosity and Darcy-Forchheimer parameters than a coarser porous media.In order to quantitatively define the effects of the input parameters for the numerical model (e.g., porosity, permeability, coefficients α f and β f ), a systematic and comprehensive analysis should be conducted considering numerical results for each of the input variables.
Finally, the new validated CFD-CSD model constitutes a valuable tool for further investigations with the aim of providing a better understanding of the hydro-geotechnical processes which may affect the stability of PBA-revetments and their foundation and to extend the results observed in laboratory tests.However, it should be considered that numerical modelling should always be employed together with laboratory tests for a systematic calibration and validation.Only under such condition, numerical modelling can be considered as a valuable tool to reliably analyse the hydro-geotechnical processes involved in the wave-structure-foundation interaction beyond the tested conditions in the laboratory.

Figure 3 .
Figure 3. Model set-up in GWK tests considered in the numerical simulations: a) model setup in GWK tests and numerical simulations, b) revetment and sand core features and c) location of pressure transducers.

Figure 6 :
Figure 6: Pressure time series comparison at different locations of a PBA revetment: a) on the revetment and b) underneath the revetment.The CFD-CSD coupling is needed to correctly reproduce the pore pressures in the soil underneath the revetment.

Table 2 . Parameters and conditions tests during the calibration of the model for comparison with GWK tests.
= degree of saturation of the porous media; α f and β f tare parameters for calculation of a and b coefficients in Darcy-Forchheimer equation. W