1
TOWARDS A PRACTICAL APPLICATION OF NUMERICAL MODELS TO PREDICT WAVE-STRUCTURE INTERACTION: AN INITIAL VALIDATION
J. van den Bos12, H.J. Verhagen2, M. Zijlema2 and B. Mellink1
A more detailed understanding of porous flow inside a rubble-mound structure may have potential benefits in breakwater design. Numerical models are expected to be a useful additional research tool in this field, provided that their calculation results can be validated against measurements. This paper presents the results of a systematic effort to validate a set of different numerical models for a strongly simplified model set-up. The results of this effort will be used to select models for a next stage in which the complexity of the model set-up is increased. By means of this step-by-step approach a good insight in the capabilities and limitations of the various models is achieved.
Keywords: rubble-mound structures; porous flow; numerical modeling; VOF modelling
INTRODUCTION Background
The study of fluid motion inside a porous structure, such as a rubble-mound breakwater or revetment, is an area of considerable research interest. The present design formulas for this type of structures tend to focus on ‘external’ behavior such as armor layer stability or wave overtopping. The interaction with the internal part of the structure is commonly included in these methods by empirical parameters. Perhaps the best-known example of such a parameter is the Notional Permeability (P) in the rock armor stability formulas of Van der Meer (1988); other examples include roughness parameters r for wave run-up and overtopping, which have a different values for ‘permeable’ or ‘impermeable’ structures. These concepts are not further quantified, i.e. it is not known at present how exactly to quantify the permeability of a structure and how the values of these empirical parameters could be calculated in detail as a function of this value.
Having a better insight into the details of the fluid motion inside the structure would be beneficial to the practical design of these structures, for a multitude of reasons. For instance (1) a better quantification of the empirical parameters such as P or r can be obtained, or – on a more ambitious level – design methods may be developed that do not rely on these parameters at all anymore; (2) problems related to interface stability (‘open filters’) on internal interfaces, which require knowledge of pore velocities and/or pore pressure gradients, can be studied better; and (3) it may open the way for the study of the ecological potential of rubble-mound structures, for instance in terms of providing shelter and/or inducing colonization of species which is also a function of pore velocities, shear stresses, turbulence et cetera.
The classical research methods in coastal engineering, which include physical model tests, have potential limitations when studying porous flow. These limitations include practical difficulties in measuring velocities or pressures inside pores, and scale effects affecting the flow regime in the porous medium. On the other hand, numerical computational methods may have fewer limitations of this kind and may be more suitable to study these phenomena. However, these numerical methods require a careful validation against measurements. This is the topic of the present paper, which describes a series of systematic validation measurements conducted at Delft University of Technology, comparing the prediction of a set of different numerical models against measurements in a physical wave flume.
The experiments were conducted for a strongly simplified structure consisting of (1) a single-layered, vertical structure, (2) with an open rear-side boundary, (3) with known and controlled porous flow parameters, (4) loaded by regular sinusoidal waves (5) on a flat horizontal seabed, and (6) measuring only the ‘macroscopic’ phenomena related to energy balance and dissipation, being reflection and transmission. This is thought to be the simplest setup possible and therefore a reasonable starting point for validation. If the numerical models pass this first validation, the complexity of the structures can be slowly increased, including multiple layers, closed rear-side boundaries, irregular waves or sloping structures. By means of this step-by-step approach a good insight in the capabilities and limitations of the various models is achieved.
1 Royal Boskalis Westminster NV, PO Box 43, 3350 AA Papendrecht, The Netherlands
2 Delft University of Technology, Delft University of Technology, Faculty of Civil Engineering and Geosciences, Department of Hydraulic Engineering, PO Box 5048, 2600 GA Delft, The Netherlands
Additional tests for the purpose of further validation have been carried out at Delft University in the meantime, but this paper only describes the initial validation on the simplified structure.
Selection of models
In this present validation, three different numerical models were studied:
(1) An analytical model based on the work of Madsen and White (1976), which uses an explicit formula to calculate reflection and transmission for exactly the type of vertical, single-layered structures studied in this paper. The model is based on the simplified physics of the Shallow Water Equations and is therefore expected to have a limited range of validity. Madsen and White (1976) report L/h0 > 10.
(2) SWASH (Zijlema et al 2011), developed at Delft University, which solves the Shallow Water Equations on a numerical grid in the horizontal direction. The model is intended to be used for shallow coastal zones to model surf zone dynamics (wave breaking, wave-current interaction, wave-vegetation interaction, transport phenomena et cetera) but is used ‘out of context’ in this paper to represent a numerical wave flume.
(3) IH2VOF (Lara 2005), developed at the Environmental Hydraulics Institute "IH Cantabria”, which uses the full VARANS equations and solves these on a numerical grid in both horizontal and vertical direction. The model is intended to be used as a numerical flume.
These three models were selected to represent a wide range of flexibility, complexity and computational effort. This corresponds to be general purpose of this paper which is to show an initial validation of numerical models and to explore the possibilities.
This selection is not intended to be complete or exhaustive. For instance, more complex models using smooth particle hydraulics (SPH), finite element modelling (FEM) or similar techniques may in the future also be able to model porous flow behavior, and potentially even couple this with movement or deformation of the armor layer. At present, these models require a very large amount of computational effort and were considered to be too impractical yet to be included here. However, as the development of these types of models is progressing, a similar validation program may be performed in the future.
PHYSICAL MODELLING Test setup
Physical model tests were carried out in the Environmental Fluid Mechanics laboratory at Delft University of Technology by Mellink (2012). The wave flume measures 38 m x 0.8 m x 1.0 m and is equipped with an Automatic Reflection Compensator (ARC) at the wave maker. The simplified ‘vertical rubble mound’ structures described above were reproduced by blocks of ElastocoastTM with varying thickness. Six of these blocks were prepared earlier for this purpose (Zeelenberg and Koote 2012) and their hydraulic properties (porosity, permeability, grain diameter) were measured. Since the stones in the ElastocoastTM blocks are bound together they form a rigid block and the properties measured by Zeelenberg and Koote remain fixed throughout the tests. This makes the use of these blocks very practical for model validation purposes.
Zeelenberg and Koote produced and measured six blocks, four of which (#2, 3, 5 and 6) were used in the tests by Mellink.
Table 1. Hydraulic properties of ElastocoastTM blocks Block # Thickness Laminar friction
coefficient
(Stationary) turb. friction coefficient
Porosity Grain size
l [mm] α0 [-] β0 [-] n [-] dn50 [m] 1 39 700 1.10 0.386 0.007 2 88 1200 1.25 0.405 0.020 3 132 1200 1.25 0.423 0.020 4 80 1900 1.70 0.410 0.039 5 160 1150 1.60 0.466 0.039 6 240 1020 1.45 0.460 0.039
The blocks were placed in the flume at a distance x = 28.0 m from the wave maker. Two arrays of three wave gauges were placed, one before the structure and one behind the structure in order to measure reflection and transmission. The wave conditions consisted of regular waves with wave heights H = 0.075 m, H = 0.10 m, H = 0.125 m and H = 0.15 m, in combination with periods T = 1.0 s,
T = 1.5 s, T = 2.0 s and T = 3.0 s. All sixteen combinations of H and T were used, for each of the four blocks, resulting in 64 tests in total. The water depth was h0 = 0.65 m in all tests.
Figure 1. Setup physical model tests (Mellink 2012)
Table 2. Wave conditions
Period Water depth
Local wave length
Ratio Wave steepness s = H/L [-]
T [s] h0 [m] L [m] L/h0 [-] H = 0.075 m H = 0.100 m H = 0.125 m H = 0.150 m
1.0 0.65 1.54 2.4 0.049 0.065 0.081 0.097
1.5 0.65 3.06 4.7 0.025 0.033 0.041 0.049
2.0 0.65 4.50 6.9 0.017 0.022 0.028 0.033
3.0 0.65 7.20 11.1 0.010 0.014 0.017 0.021
Analysis and results
The recorded time series at the six wave gauges were analyzed by decomposing them into an incoming and reflected signal using the two-probe method of Goda and Suzuki (Goda 1985), using the first two probes of each array only.
During the analysis of these tests, it became clear that the tests with T = 1.0 s were less reliable because (1) the very short waves developed a lateral standing wave pattern (across the width of the flume) which could not be compensated by the ARC which works only in the length direction of the flume, and (2) the spacing of the wave gauges (Δx12 = 0.82 m = 0.53.L) was not ideal for these short waves, leading to less accurate reflection and transmission values. For these reasons, the tests with T = 1.0 s are not included in the validation of this paper (although they were modelled in the various numerical models and were included in the grid convergence analyses). Finally, in the tests some unwanted reflection off the far side of the wave flume was observed, for which the measurements were corrected (Mellink 2012).
Figure 2. Impression of physical and numerical modelling.
NUMERICAL MODELLING Analytical model
An analytical model was proposed by Madsen and White (1976) based on the Shallow Water Equations (SWE) for surface elevation and horizontal flow velocity U:
0
0
0x
g
t
U
x
U
h
t
(1)in which h0 is the constant water depth. Inside a porous structure these equations change to
0
0
0n
f
x
g
t
U
n
S
x
U
h
t
n
(2)in which n is the porosity of the material, = 2/T is the wave frequency, S = 1+(1-n) is a coefficient describing the effect of unsteady flow, is an added mass coefficient (Madsen and White report values = 0 – 0.5) and f is a linearized friction factor:
u
n
f
3
8
(3) with
2 50 2 3 01
nd
n
n
(4) 50 3 01
1
nd
n
n
(5)in which α0 and β0 are the laminar and turbulent Forchheimer coefficients of the porous medium, respectively.
Madsen and White propose harmonic solutions of the type =Re{(x)eit} and U=Re{u(x)eit}. The solution is the well-known combination of two waves travelling in opposite direction. In the domain before the structure (x < 0):
x ik r x ik i
e
a
e
a
0
0
(6)
ikx
r x ik ie
a
e
a
h
g
u
0 0 0
(7)In the domain after the structure (x > l) there is no wave travelling in negative direction (only the transmitted wave): ) ( 0 x l ik t
e
a
(8) ) ( 0 0 x l ik te
a
h
g
u
(9)In these equations k0 = /√(gh0) is the wave number in the absence of friction, and ai, ar and at are the amplitudes of the incoming, reflected and transmitted waves respectively (ar and at are complex numbers). Finally, inside the structure the solution is influenced by the friction in the porous medium:
) (x l ik ikx
e
a
e
a
(10)
( )
0 l x ik ikxe
a
e
a
h
g
u
(11)In which k = k0√(S-if) is the complex wave number under the influence of friction, ε = n/√(S-if) and a+ and a- are the (complex) wave amplitudes inside the structure. The system can be solved by setting continuity boundary conditions for and u at (x = 0) and (x = l), giving the explicit solution (Madsen and White 1976):
ikl
ikl i te
e
a
a
2 21
1
4
(12)
ikl
ikl ikl ikl i re
e
e
e
a
a
2 2 21
1
1
(13)And finally the real reflection and transmission coefficient follow as the amplitude of these complex expressions: R = |ar/ai| and T = |at/ai|. The solution of the system now depends on the value of the friction coefficient f and the added mass coefficient S. Madsen and White note that S ≈ 1 for practical purposes and propose the simple relationship S = (n/0.45)2. The friction coefficient can be found as part of the solution from the definition (equation 3).The complete solution requires an iteration on f and |u|.
In this paper two different versions of the model are used. One is the original version as described in the original publication (Madsen and White 1976). The second, ‘modified’ version deviates in three ways:
(1) Madsen and White state that the iteration on f and |u| is “computationally expensive”, and an approximation method is proposed. On modern computers such an iteration is however no longer an issue. In the version used in this paper, this iteration is implemented in Matlab; (2) Van Gent (1995) reports that the value of the turbulent friction coefficient β0 is not constant but
KC
stat5
.
7
1
, 0 0
(14) 50 nnd
T
u
KC
(15)in which KC is the wave Keulegan-Carpenter number and β0,stat is the stationary turbulent friction coefficient (as measured by a simple permeability test on a soil sample, i.e. the parameter given in table 1). This can simply be implemented as another step in the iteration on f and |u|;
(3) The explicit solution for at/ai and ar/ai given by Madsen and White is valid for a single porous layer. If more layers are allowed then the system can be expanded by setting continuity boundary conditions on the additional interfaces (and solving the iteration on f, β0 and |u| for each layer). In that case no explicit solution can be formulated, but instead the model must be solved by linear algebra. The version of the model used in this paper uses this notation, with an eye on possible future expansion to multiple layers, even though in the present tests only single layers were used.
SWASH
The program SWASH is also based on the Shallow Water Equations. It solves these numerically on a Eulerian grid with fixed intervals Δx. Being based on the SWE, there is no solution in the vertical direction, although the vertical domain can be divided into K layers for increased calculation accuracy. Some of the grid cells can be given a nonzero porosity in order to represent porous structures. The porosity can be different for each cell, however the associated friction parameters (α0, β0) can only have a single value for the entire domain. The value for β0 is also treated as stationary, i.e. the correction β0 = β0,stat(1-7.5/KC) is not implemented. Moreover, β0 is limited to the range 1.8 < β0 < 3.6.
Figure 4. Grid convergence numerical models. Top left: SWASH. Top right: IH2VOF (after modification). Bottom: IH2VOF (effect of number of grid cells per wave height and wavelength)
Following the standard recommendations in the SWASH manual a grid size Δx = L/(50 to 100) should be chosen. Given the range of wave lengths present in the tests (see table 2), Δx = 0.03 m was selected. In the vertical, K = 2 layers were chosen. The initial time step was Δt = 0.005 s giving a Courant number Cr = Δt.c/ Δx = 0.42 (using the shallow water approximation c = √gh0), which is below the recommended threshold Cr < 0.5. In order to verify whether this grid size was small enough (i.e. not leading to any inaccuracies in the calculation), the calculations were also made on a smaller grid Δx = 0.02 m (with Δt = 0.0035 s giving Cr = 0.44) and the results were compared, see figure 4. No significant deviations were observed and it was concluded that the calculation grid Δx = 0.03 is adequate.
The model was run for all 64 tests described above, using the same ‘numerical flume’ length as the physical flume length, the same location and parameters of the porous block, the same wave conditions and the same water depth. The waves run in the model from x = 0 (left hand side) to x = 38.0 m (right hand side). Beyond x = 38 m, a sponge layer with a width of 36 m was implemented in order to absorb the incoming wave and avoid reflections back into the calculation domain. The water surface elevations calculated by SWASH were recorded as a function of time, at six locations corresponding exactly to the locations of the wave gauges in the physical model. The thus obtained time series were then analyzed to find the reflection and transmission coefficients in much the same way as was done for the physical model. The reflection coefficient is then calculated as R = ai1/ai1 and the transmission coefficient T = ai2/ai1. The effectiveness of the sponge layer was checked by observing that ar2 was sufficiently low. It was indeed found that the maximum value was ar2,max = 0.0012 m for all tests.
The incoming and reflected waves were decomposed using the two-probe method of Goda and Suzuki, which is included in the standard analysis software in use at the DUT laboratory. The method was applied pair-wise (i.e. gauges 1-2, 2-3, 1-3 and 4-5, 5-6, 4-6). Results were then averaged. If the mutual distance of a pair of gauges was leading to inaccurate results (as indicated by the DUT software), then the results from that pair were not included in the average. The first and last part of the measured signal were removed from the analysis in order to avoid model startup and shutdown effects (if any).
Figure 5. Example of recorded time series showing start-up effects. Left: SWASH. Right: IH2VOF.
IH2VOF
The IH2VOF model solves the full Volume Averaged Reynolds Averaged Navier-Stokes (VARANS) equations, in both x and y direction, on a Eulerian numerical grid in two dimensions. The location of the free surface is tracked using a Volume of Fluid (VOF) routine. A porous structure can be implemented in the model and can in principle consist of multiple layers of any orientation, not just vertically. For each layer the porous parameters (n, α0, β0 and dn50) can be defined. The model then calculates the velocities based on the Forchheimer equation, including the correction β0 = β0,stat(1-7.5/KC). Note that the model is volume-averaged, in other words it gives the average (filter) velocity over a porous grid element Δx.Δy. The real pore velocities deviate from this value.
In this model a grid size must be selected in two directions: horizontal (Δx) and vertical (Δy). The time step is selected automatically by the software in order to satisfy the Courant condition. The software manual recommends to use approximately 10 vertical grid cells per wave height (H/Δy = 10) in combination with an aspect ratio Δx/ Δy = 2.5. On this basis an initial grid size Δx = 0.03 m and Δy = 0.012 m was selected. Like before, the calculations were then repeated for a smaller grid size and compared to verify the grid convergence. This time, deviations of up to 20% were observed (especially for the shorter waves), so a yet smaller grid size was used for further comparison. It seemed that in
particular the number of horizontal grid points per wave length L/Δx was important. This procedure was repeated until a satisfactory convergence of results was obtained. This analysis was done for Block 6 (tests 49-64) only because of the large computational effort involved.
The results of this analysis are presented below as functions of both H/Δy and L/Δx. It follows that the results are indeed sensitive for L/Δx, and that a recommended minimum number of grid points is around L/Δx = 150. With fewer points the results become inaccurate. There is less sensitivity for the number of vertical grid points per wave height H/Δy. In the remainder of the calculations, all grids were based on these conclusions.
Table 3. Grid size IH2VOF modelling
Period Local wave length
Grid size Ratio
T [s] L [m] Δx [m] Δy [m] L/ Δx [-] Δx / Δy [-]
1.0 1.54 0.012 0.008 128 1.5
1.5 3.06 0.020 0.008 153 2.5
2.0 4.50 0.030 0.012 150 2.5
3.0 7.20 0.030 0.012 240 2.5
The model results were then analyzed in the same way as in the SWASH modelling, i.e. a time series of calculated surface elevations was recorded at the six locations of the wave gauges in the physical model, and the reflection and transmission coefficients were obtained using pair-wise application of the 2-probe method of Goda and Suzuki. As before, only the part of the signal that was not influenced by model start-up or shutdown effects was used in the analysis. The effectiveness of the absorbing right-hand side boundary of the model was also verified by observing that the reflected wave behind the structure was sufficiently low (it was found that ar2 < 0.0018 m for all tests).
The length of the modelled time series was always 100 s. The calculation time required for the IH2VOF model varied between approximately 20 minutes for the longer waves to roughly 1 hour for the T = 1.5 s tests and over 3 hours for the very short waves (T = 1.0 s). This performance was achieved on a standard laptop with Windows Vista 32-bit, Intel i7-2640M CPU @2.80 GHz and 4 GB RAM. No attempt was made to use dual cores or parallel computing.
The IH2VOF model tends to generate a lot of data, but since only time series of water surface elevations were saved, the total amount of output data for these present tests was limited. In total 7.4 GB of data was produced for all tests combined. By comparison, if all data for the entire computational domain would have been saved (as would be required to create e.g. the plots shown in figure 2 a much larger volume of data would have resulted, amounting to approx. 1.4 GB for a single test.
RESULTS
The reflection and transmission coefficients calculated by the various models can now be compared directly with the measurements from the physical model tests. This is done in two ways: a visual comparison on a scatter plot showing calculated vs measured data pairs, and a quantitative assessment of the mean prediction error:
measured calculated measured tests RR
R
R
N
Err
1
(1)
measured calculated measured tests TT
T
T
N
Err
1
(1)Since some of the models are based on the Shallow Water Equations it is expected that the predictive power of the numerical models is a function of the validity of that approximation and therefore of the dimensionless water depth (ratio L/h0). In order to test this hypothesis, the mean prediction errors were calculated for (1) the entire population containing all tests, (2) omitting the tests with T = 1.5 s, and (3) further omitting the tests with T = 2.0 s. Note that the tests with T = 1.0 s are omitted from this comparison at all because of the presumed inaccuracies in the physical model.
(1) For reflection, the modified Madsen and White model, SWASH and IH2VOF perform almost equally well, and are able to reproduce the measured values with order 5% accuracy. There is no discernable effect of the L/h0 ratio, although the results for the lower wave period (L/h0 = 4.7) are slightly worse. The original Madsen and White model performs less well, with accuracy in the order 10%.
(2) Transmission is predicted less well than reflection, and the effect of the L/h0 ratio is more pronounced. Overall the models have an accuracy of order 20% for low L/h0 ratios, improving to order 5-10% for higher ratios. The performance of the models improves with their complexity: the IH2VOF model performs better overall than the modified Madsen and White model, and SWASH sits in between. Again, the original version of the Madsen and White model performs less well, this time to a quite considerable degree with an accuracy of 30-50%.
Table 4. Validation results – mean prediction errors
Population Validity Mean prediction error ErrR Mean prediction error ErrT M&W
(1976)
Modified M&W
SWASH IH2VOF M&W (1976) Modified M&W SWASH IH2VOF T 1.5 s L/h04.7 0.110 0.079 0.054 0.076 0.481 0.228 0.237 0.128 T 2.0 s L/h06.9 0.122 0.049 0.044 0.068 0.339 0.160 0.118 0.087 T 3.0 s L/h011.1 0.134 0.052 0.048 0.062 0.281 0.151 0.067 0.063
Figure 6. Validation results – mean prediction errors
CONCLUSIONS
It is concluded that the numerical models studied in this paper are able to give a rather fair reproduction of reflection and transmission for this very simplified set-up. This provides a good basis for further validation on more complex structures.
Regarding the analytical model, it is concluded that the modifications introduced in this paper (full iteration on friction and velocity and introduction of the non-steady correction for β0 proposed by Van Gent 1995) do improve the performance of the model considerably. Also, the validity of the model appears to extend to lower values of dimensionless water depth (tested up to L/h0 > 4.7) than the use of the Shallow Water Equations would theoretically allow. For future use on more complex structures, it is expected that the matrix notation would allow for application of the model for multilayered (though still vertical) structures and closed rear-side boundaries. The model has an added advantage over the other, purely numerical models in the sense that it provides insight into relevant physical parameters such as friction coefficients.
Figure 7. Validation results – scatter plots for T 1.5 s
Comparing SWASH and IH2VOF, it is concluded that both models perform equally well in their prediction of reflection, and that IH2VOF performs slightly better in its prediction of transmission. SWASH on the other hand is more efficient in terms of computational effort, and the increased accuracy of the IH2VOF model does not seem to outbalance that. However, SWASH is expected to be less versatile in modelling the more complex structures: sloping structures cannot be modelled altogether, but also multilayered structures cannot be modelled at present because the model only allows a single value for α0 and β0 (although the porosity n can have multiple values). SWASH is in that aspect even less versatile than the modified analytical model, while its prediction performance on these simplified structures is not significantly better. IH2VOF, by contrast, is expected to be able to model other structures no matter their complexity, including multilayered or sloping structures.
On the basis of these conclusions it is decided that future work will use the IH2VOF model for the modelling of more complex structures, while the modified analytical model may also be used to provide a theoretical framework (pending the validation results on the more complex structures).
REFERENCES
Goda, Y (1985) Random seas and the design of maritime structures, Univ. Tokyo Press, Tokyo, Japan Lara, J (2005) A numerical wave flume to study the functionality and stability of coastal structures,
PIANC magazine, nr 121
Mellink, B (2012) Numerical and experimental research of wave interaction with a porous breakwater, MSc thesis Delft University of Technology, Delft, The Netherlands
Madsen O, White S (1976) Reflection and transmission characteristics of porous rubble-mound breakwaters, Technical Report 76-5, US Army Corps of Engineers, Coastal Engineering Research Center, Fort Balvoir, USA.
Van der Meer, J (1988) Rock slopes and gravel beaches under wave attack, Ph.D. thesis, Delft University of Technology, Delft, The Netherlands
Van Gent, M (1995) Wave Interaction with Permeable Coastal Structures, Ph.D. thesis, Delft University of Technology, Delft, The Netherlands
Zeelenberg, W and Koote, M (2012) The use of Elastocoast in breakwater research, Technical report, Delft University of Technology, Delft, The Netherlands
Zijlema M, Stelling G, Smit P (2011) SWASH: an operational public domain code for simulating wave fields and rapidly varied flows in coastal waters, Coastal Engineering (58) pp 992-1012.
!["Wacława Potockiego Moralia (1688)", wyd. Tadeusz Grabowski, Jan Łoś, [Kraków] 1915-1918 : [recenzja]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACH5BAEAAAAALAAAAAABAAEAAAICRAEAOw==)