# Relaxation zones

2.2.1 Background and development

The relaxation zone method for generating and absorbing waves is an abstract method that is based on the concept of dedicating a part of the CFD model domain to perform wave generation and/or absorption. The concept can be very broadly summarised in the following equation, which, without loss of generality, is written in ID form:

where *F* is a model variable or parameter in broad terms (e.g., velocity, free-surface elevation, pressure, mass/momentum sink scaling parameter, numerical parameters, mesh size, among others), _{r} is the coordinate along the wave propagation axis, with *x _{r}* = 0 at the interface of the numerical domain and the relaxation zone, and

*L*the length of the relaxation zone along the ay-axis. Variables or terms with indices

_{r}*b*(boundary) and

*cl*(domain), indicate that these are calculated or imposed at the boundary or the internal domain, respectively. The application of this concept in a NWT is shown in Figure 1.

In the context of wave absorption and generation problems, it has been demonstrated that an appropriate form of Equation 4 can be introduced to the flow equations (physics of the model) or the numerical solution procedure (numerics of the model) to relax the solution towards the boundary condition, within the relaxation zone domain. It has been additionally shown that given enough length, this transition becomes smooth and reflections of outgoing waves are minimised. Techniques based on Equation 4 may be referred to as “sponge layer”,

Figure 1: Illustration of relaxation zones in a NWT, according to (Jacobsen et al., 2012).

“numerical beach”, “generation/absorption layers” among others. In this paper, all these techniques are considered to fall under the term “relaxation zone met hods/techniques”. It is worth noting that to classify as such, the relaxation zone must share an interface with the domain boundary, otherwise it falls under the internal wave maker techniques, which will be discussed in the following section.

Early work on passive absorption

Relaxation zone methods for NWTs were first applied to absorb outgoing waves at the inshore/landwards outlet, in combination with a Dirichlet boundary condition or an internal wave maker, rather than used to fully equip a NWT. In this case, the boundary condition at the outlet can be a no-slip wall (as in an actual wave tank) or a free-flow outlet (constant sea level). The purpose of the relaxation zone is to entirely dissipate the flow velocities (passive absorption) before the waves reach the outlet boundary.

The application of the relaxation zone presents a physical analogy with layouts used in the laboratory for absorbing waves such as absorbing foam, screens and beaches. In this case, the absorbing layout usually takes over few or several meters of the experimental flume or tank (as opposed to a paddle, which is the equivalent of a boundary condition in numerical models) and the absorbing layout properties (primarily the resistance it offers to the flow) are gradually increased towards the outlet wall, in a manner consistent with Equation 4. Examples of absorbing layouts in experiments can be found in Tiedeman et al. (2012) and Delafontaine (2016), among others.

An early investigation of water wave absorption using a relaxation zone was conducted by Le Mehaute (1972), by introducing an “internal friction term” in the dynamic (Bernoulli) equation of motion for the pressure, in the context of potential flow equations. The friction term is proportional to the flow potential and is broadly equivalent to a Darcy resistance term in the Navier-Stokes equations. This technique can be considered as “physical forcing” of the relaxation zone. Indeed, the physical analogy presented in Le Mehaute (1972) is parallel plates with varying density aligned with the wave direction. This may be now considered primitive, but it is consistent with the concept of increasing viscous resistance towards the outlet boundary, which is used in modern absorption techniques in the laboratory (e.g., screens, foam or beach). The concept of adding an internal friction term in a relaxation zone to facilitate wave absorption has nevertheless proven timeless, as many modern depth-integrated and CFD models still utilize this idea (Peric and Abdel-Maksoud, 2015; Dimakopoulos et al., 2019).

A generic formulation for using a sponge layer for passive absorption of waves was established by Israeli and Orszag (1981), where it was demonstrated that sponge layers, as opposed to simple radiation conditions, were capable of absorbing multiple wave frequencies.

They also demonstrated that a combination of radiation conditions and relaxation zones can optimise the absorption efficiency. Early applications of this concept can be found in depth integrated or potential flow models which matured earlier than CFD models for ocean and coastal engineering applications. In the Boussinesq model proposed by Larsen and Dancy (1983), a sponge layer method is applied by “dividing, at each timestep, the surface elevation and the flow on a few grid lines next to the boundary by a set of numbers which increase towards the boundary”. In essence, they use a weighted denominator to divide the pressure and the free-surface at each timestep, right after these values are calculated by the numerical solution, and it is demonstrated that this approach efficiently absorbs outgoing waves. Their scheme is coupled with an internal source generator and it is argued that active absorption capabilities are not needed, as long as the sponge layer is introduced in all domain boundaries. The concept of introducing intermediate steps in the numerical algorithm to absorb waves, survives in many modern NWT approaches, based on CFD (e.g., in Jacobsen et al. (2012)). This concept will be hereafter mentioned as “numerical forcing”, as opposed to “physical forcing” which is based on introducing appropriate source terms in the momentum equations, while preserving the original numerical solution algorithm.

Cao et al. (1993) argued that the formulation of Le Mehaute (1972) is not bounded at positive values and could sometimes result in negative dissipation, which is not desirable. Grilli et al. (1997) proposed an improved implementation of the relaxation zone technique for potential flow equations, within the context of NWTs using the Boundary Element Method (BEM). Their method included a modified dynamic free surface boundary condition for the pressure, resistant to the motion of the free surface. The influence of the modified boundary condition was increased approaching to the boundary by using polynomial weighing function with an order of 2-3. This technique was shown to achieve good absorption characteristics for short waves (< 3% reflection) achieved with a 2 wavelengths long absorption zone, but relatively poor performance for longer waves (up to 40%). This is because the modified boundary condition affected the wave kinematics close to the free-surface, thus being more appropriate for deep water waves, where the wave energy is concentrated in the upper water layers. To improve performance for longer waves, Grilli et al. (1997) combined the relaxation zone method with an absorbing piston at the outlet based on the Sommerfeld equation, which is known to yield good absorption performance for longer, shallower waves. The relaxation zone approach proposed by Grilli et al. (1997), without radiation conditions, was adopted for a NWT based on a single phase Navier-Stokes model with a moving free-surface boundary in Dimas and Dimakopoulos (2009), where good performance was achieved, for intermediate waves, with a 4 wavelengths long relaxation zone.

The application of the relaxation zone method for passive absorption was also further elaborated in Boussinesq modelling by adding linear resistance terms in the momentum equations (Kirby et al., 1998). A step further was taken in Nwogu and Demirbilek (2001), by adding a linear mass sink term. Whilst this term interferes with mass conservation, Nwogu and Demirbilek (2001) recommend its inclusion, arguing that it further enhances absorption performance. This is the physical equivalent of, e.g., adding appropriate overflow arrangements in a numerical wave tank, with the overflow level reaching the mean water level as the onshore boundary is approached.

A relaxation zone method for wave absorption in CFD models was proposed by various researchers, starting from the late 90's. Mayer et al. (1998) used a “numerical forcing” approach which relaxes the free surface elevation and velocities to prescribed values at each time-step, using a 3rd order polynomial equation for *a.* Note that a fractional time-step scheme with pressure correction is employed for solving the Navier-Stokes equations. Whilst Mayer et al. (1998) do not report the absorption efficiency, they demonstrate a good overall comparison with analytical solutions and experimental data, including structures interacting with waves and currents. The COBRAS model (Lin and Liu, 1999) also employed the formulation of the Israeli and Orszag (1981) approach, coupled with an internal source wave maker, and has since been used for numerous applications during the late 90’s and 00’s (see section on internal wave makers below).

Early work primarily comprised passive absorption schemes, while wave generation was mostly achieved by using internal wave makers. More recently, CFD-based NWTs were developed by fully utilising the relaxation zone method, once this concept was demonstrated to be capable of simultaneously absorbing and generating waves at the boundary Jacobsen et al. (2012). This process has the same result as the concept of the active wave/generation absorption applied to absorbing boundary conditions and numerical wave paddles, but strictly speaking, it cannot be described as active, as the relaxation zone method does not adapt to the local wave properties. In order to distinguish between the two, the term “simultaneous wave absorption and generation” (SWAG) is hereafter used in the context of using relaxation zone method for both wave generation and absorption.

Simultaneous wave absorption and generation

The conversion of the passive absorption zone techniques presented above to SWAG techniques are -in hindsight -relatively straightforward, particularly for the ones that rely on the notion of dissipating waves by forcing the solution to constant or zero target velocity and free surface elevation values. The generation can be achieved by introducing a target velocity or free surface that corresponds to a wave theory. An early demonstration of this idea was performed by Afshar (2010) who, in terms of classification, stands between the relaxation zone and the internal wavemaker concept (see Figure 2). A layout with three relaxation zones is proposed, with the first two located at the offshore/generation boundary (Zone I and Zone II) for “ramping up” and “ramping down” the wave field. A third relaxation zone (Zone III) is introduced at the outlet for wave absorption. Generation and outlet boundaries were considered as solid walls, and in this sense, the overall scheme approaches the idea of an internal wave maker. Pressure, velocities and volume fractions were relaxed with 3rd and 6th order functions, using “numerical forcing”. Zone I and II were one wavelength long while Zone III was two wavelengths. The method was tested for linear, standing and nonlinear waves and was demonstrated to be successful at a conceptual level. Several areas of improvements were identified, including the interpolation of forcing terms between cell centres to avoid the appearance of a saw-tooth profile in generation zones and the removal of high air velocities in the relaxation zone.

Jacobsen et al. (2012) devised the waves2Foam library which was written for the OpenFOAM® platform and was fully based on the relaxation zone method, as the relaxation zone concept was combined with a Dirichlet condition for wave generation at the offshore boundary, as in Figure 1. A substantial improvement was achieved over the early concept

Figure 2: Illustration of relaxation zone from Afshar (2010).

of Afshar (2010). The proposed technique addressed the interpolation issues by projecting the location of the free-surface in a boundary face or a computational cell and calculating the water fraction based on the interpolated volumes. Air velocities were forced to zero through the relaxation zone technique, and pressure forcing was made redundant. The relaxation function had an exponential form, following recommendations from Mayer et al. (1998). The method was shown to be able to efficiently generate and absorb waves of linear and nonlinear regimes with high efficiency, as wave reflections were < 3% for linear and nonlinear waves and relaxation zone length of two wavelengths.

This idea was largely successful and was picked up and further developed by the research community. Dimakopoulos et al. (2016) verified a similar level of reflection efficiency by testing the waves2Foam library against regular wave conditions proposed by Higuera et al. (2013). In addition, they further demonstrated the capability of the toolkit as a SWAG scheme for generating regular, random, directional waves, with or without currents. The capability of the scheme to simultaneously absorb regular and random waves was demonstrated through numerical tests of wave reflection against solid walls. The evolution of wave height of standing regular and random waves was compared against analytical solutions and agreement was excellent for low-steepness waves. It was also further demonstrated that the method is capable of performing well within the context of a 3D numerical basin, showing a good performance for generating short-crested waves, as well as wave and current interaction in oblique directions.

The relaxation zone technique of Jacobsen et al. (2012), is also adapted to the Reef3D model (Miquel et al., 2018), and includes numerical forcing of the dynamic pressures. Their implementation proved to be highly efficient for random waves absorption (< 2%) and less efficient for regular waves (< 10%), something that seems counter-intuitive and contradicts findings from Jacobsen et al. (2012) and Dimakopoulos et al. (2016). It should nevertheless be noted that the implementation of the method in Reef3D includes relaxat ion of the pressure terms and the reflection analysis that they use is more likely to be suitable for random rather than regular waves, so this may justify the different findings. Additional models that use the relaxation zone approach are the OpenFOAM®-based toolkit Naval Hydro pack (Jasak et al., 2015) or Star-CCM CFD platform (Peric et al., 2018).

As the relaxation zone technique became increasingly widespread, further analysis was performed regarding the parameters relating to the relaxation zone. Peric and Abdel- Maksoud (2016) used mass and momentum sources by setting *Fb = F _{s}cr(x)(f* —

*f*

_{t}) and

*Fj =*0 in Equation 4, where /,

*ft*are the calculated and target flow variables (volume of fluid or velocity) and

*F*is a scaling coefficient. Their analysis showed that

_{s}*F*should be proportional to the wave frequency while the length of the relaxation zone should be proportional to the wave length for the absorption zone to achieve relaxing capabilities. In addition, they demonstrated that the optimal values of these parameters are

_{s}*F*= 7Г

*ш*and

*x*= 2Л for an exponential type relaxation function. Similar values were consistently used by researchers using trial and error approaches, e.g., see Nwogu and Demirbilek (2001) or Jacobsen et al. (2012). In Peric and Abdel-Maksoud (2018), a theory was developed to predict reflection coefficients for monochromatic waves, given the parameters of the relaxation zone. The theory was found to agree with numerical tests and the researchers or engineers can use it for selecting appropriate parameters for the relaxation zone conditions, given a target reflection coefficient. The authors nevertheless recommend further work on extending the theory to cover irregular or random waves.

In Dimakopoulos et al. (2019), a relaxation zone technique is proposed, which is very similar to the one presented in Peric and Abdel-Maksoud (2016). Dimensional analysis confirms findings from Peric and Abdel-Maksoud (2016) that the scaling factor *F _{s}* must be proportional to the wave frequency to achieve good performance and the value of the scaling coefficient is set to 5 times the angular wave frequency, following recommendations from Nwogu and Demirbilek (2001). The method is implemented in the CFD toolkit Proteus (Kees et ah, 2011), showing good performance for generating and absorbing random waves, although the approach can be possibly improved by modifying the mass conservation equation, as the method only considers momentum sources.

2.2.2 Examples and applications

The relaxation zone method for CFD models became widespread after the release of the waves2Foam library and it is usually combined with the OpenFOAM® CFD toolkit. The method proposed by Jacobsen et al. (2012) has been used for a wide range of applications, including modelling wave structure interaction processes in seawalls and breakwaters (Richardson et al., 2014; Jensen et al., 2014; Medina-Lopez et al., 2015; Hu et al., 2016; Jensen et al., 2017), marine and offshore structures (Palm et ah, 2013; Fuhrman et ah, 2014; Eskilsson et ah, 2015; Jasak et ah, 2015; Moradi et ah, 2016; Vvzikas et ah, 2017; Chen and Christensen, 2018), naval and ship hydrodynamics (Piro, 2013; Winden et ah, 2014; Shen and Korpus, 2015; Vukcevic et ah, 2017; Kim and Lee, 2017; Liu et ah, 2018) and soil, sediment and scour processes (Stahlmann and Schlurmann, 2012; Jacobsen and Fredsoe, 2014; Karagiannis et ah, 2016; Cheng et ah, 2017; Elsafti and Oumeraci, 2017). The list is not exhaustive as many more researchers and engineers have used the relaxation zone method in the context of the OpenFOAM® CFD toolkit.

The relaxation zone method has been also used for modelling engineering applications in the coastal and marine environments in other models as well, such as the Reef3D model (Bihs et ah, 2015; Chella et ah, 2015; Kamath et ah, 2016; Ong et ah, 2017)) and the Proteus CFD toolkit (Cozzuto et ah, 2019; de Lataillade et ah, 2017; Mattis et ah, 2018), covering a broad range of applications, similar to the one listed for OpenFOAM®.

2.2.3 Advantages and disadvantages

The main advantage of the relaxation zone method is primarily the high fidelity of wave generation along with high efficiency in absorption. This is because the method is not susceptible to boundary instabilities due to the sudden change of flow equations, while forcing analytical solutions or linear approximations meaning that evanescent modes and parasitic waves (typical in moving paddles or internal wave makers) are not present or very quickly dissipated.

Another advantage is the ease of implementation. It is relatively straightforward to add source terms to existing equations or to use “numerical forcing”, without particular considerations for numerical stability and accuracy, as opposed to, e.g., radiation boundary conditions. For example, the method used in (Dimakopoulos et ah, 2019) performed a simple adaptation of the Darcy term in the porous media equation to implement SWAG. This method is therefore suitable to implement in newly and/or rapidly developed models, where for example an adaptation of an existing model component may allow the model to be used as a NWT. The inclusion of the mass source term through the volume of fluid scheme may require some treatment to, e.g., avoid non-physical oscillations at the generation zone (Afshar, 2010; Jacobsen et ah, 2012). Whilst this may require some additional efforts, including mass sources is recommended, as it has been reported that it enhances absorption (Nwogu and Demirbilek, 2001) and it is also capable of maintaining the water level in the domain constant, thus not needing corrections for influx due to a mass imbalance between generation of wave troughs and crests.

The relaxation zone also needs a suitable selection of parameters such as relaxat ion function, length and scaling of the coefficient , thus being subject to uncertainties. Recent research and engineering practice have addressed these issues. Based on the literature consensus, the authors recommend the following:

- • Length of generat ion zone:
*L*~ Л (corresponding to peak frequency for random waves)_{r} - • Length of absorption zone:
*L*~ 2A (corresponding to peak frequency for random_{r}

waves)

- • Scaling parameter: F
_{s}=3-5*ш* - • Relaxation function:
*a*= '^{x|},’_{x[)}i—-j—- ■

Some of these values are theoretically proven (see Peric and Abdel-Maksoud (2015) for regular waves. For the remaining parameters or wave regimes, these recommendations are demonstrated to have a reasonable performance throughout multiple studies in literature. Note that frequency and wavelength parameters correspond to peak frequency and wavelength for random wave regimes.

Undoubtedly the greatest disadvantage of the method is the additional computational resources required for running engineering applications, compared to other methods presented in this article. Regarding wave absorption, the domain must be extended by about two wavelengths and this can significantly increase the length of the domain. Peric et al. (2018) also demonstrated that a careful selection of the relaxation zone parameters may allow for a reduction of the size of the absorption zone, e.g., using 0.5-1 wavelenghts, rather than two. (Dimakopoulos et ah, 2016) have demonstrated that increasing the mesh size with a mesh progression ratio of 10%, the mesh cells in the absorption zone can be reduced by > 95%, whilst maintaining equal or better absorption efficiency. A similar technique can be applied in unstructured meshes (de Lataillade, 2019). The caveat of this approach is that the free surface tracking scheme should be robust enough in order to cope with changes in mesh size along the free surface without introducing spurious oscillations.

Dedicating up to a wavelength for SWAG cannot essentially be considered as an increase of the computational domain, as for all generation methods, it is considered good practice to leave a buffer zone after the generation boundary to allow the waves to develop before entering the main domain (e.g., approaching a slope or a structure). In this case, the role of the buffer zone can be played by the generation zone. It was nevertheless demonstrated by Dimakopoulos et al. (2016) that the computational cost associated with generating nonrepeating random wave sequences for, e.g., simulating a storm event may be disproportionate to the overall cost of the simulation. For example, in Dimakopoulos et al. (2016), it was demonstrated that using 500 wave components for generating a random wave series which corresponds to 250 non-repeating waves doubles the simulation time in a NWT of 240,000 mesh cells, when compared to simulations that use 50 components. This is primarily due to the fact that trigonometric and hyperbolic functions are rather expensive to calculate for each cell and each time step in the relaxation zone.

Several approaches have been proposed to address this problem. Dimakopoulos et al. (2016) showed that by replacing native С++ and FORTRAN trigonometric functions with 4th order Taylor approximations results in a reduction of the computational cost by an order of magnitude, whilst the approximation errors remained < 1%. This reduced the overall cost but it did not completely eliminate it, e.g., simulations with 500 components in the NWT mentioned before need 1.4 more times, rather than double with the native functions. In Jacobsen (2017) a stretched distribution of the spectrum is proposed for improving nonrepeatability of the series, thus needing less wave components for describing a full storm, but this approach does not offer a consistent discretisation of the spectra in the frequency range away from the peak. Similarly the computational time in Jacobsen (2017) may be further reduced by manipulating the trigonometric and hyperbolic functions to split spatial and temporal terms, storing spatial terms in memory and only calculating temporal terms during time advancement. This has demonstrated to speed up the calculation time needed to compute the signal by 50 times for 1000 wave components. More recently, Dimakopoulos et al. (2019) developed a method based on pre-processing the random wave signal using spectral window decomposition methods. Their technique was demonstrated to be able to generate non-repeating signals of 0(1000) waves, size using 0(10) frequency components, thus achieving similar reductions in computational cost. The technique was tested in the CFD toolkit Proteus.

Overall, the relaxation zone is an accurate and easy technique to implement for generating and absorbing all wave regimes, but the method is indeed computationally expensive. For its merits, thhe relaxation zone is one of the most preferred methods to use in NWT and researchers have been making significant advances in addressing its drawbacks relating to speed. Including these advances in future works will further increase the attractiveness of the technique.