Analytical Prediction of Reflection Coefficients for Wave Absorbing Layers in Flow Simulations of Regular Free-Surface Waves

TL;DR

This paper develops an analytical theory to predict reflection coefficients in wave absorbing layers, enabling optimal parameter selection for flow simulations of regular free-surface waves, reducing boundary reflections effectively.

Contribution

It introduces a general analytical model to predict wave reflection coefficients, aiding in the optimal tuning of absorbing layer parameters in flow simulations.

Findings

The theory accurately predicts reflection coefficients in simulations.

Validation shows satisfactory accuracy for practical applications.

The model is applicable to various implementations of absorbing layers.

Abstract

Undesired wave reflections, which occur at domain boundaries in flow simulations with free-surface waves, can be minimized by applying source terms in the vicinity of the boundary to damp the waves. Examples of such approaches are absorbing layers, damping zones, forcing zones, relaxation zones and sponge layers. A problem with these approaches is that the effectivity of the wave damping depends on the parameters in the source term functions, which are case-dependent and must be adjusted to the wave. The present paper presents a theory which analytically predicts the reflection coefficients and which can be used to optimally select the source term parameters before running the simulation. The theory is given in a general form so that it is applicable to many existing implementations. It is validated against results from finite-volume-based flow simulations of regular free-surface waves and found to be of satisfactory accuracy for practical purposes.

Figures

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V\rho \mathrm{d}V+{\int }_S\rho (\text{v}-{\text{v}}_g)\cdot \text{n}\mathrm{d}S=0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V\rho u_i\mathrm{d}V+{\int }_S\rho u_i(\text{v}-{\text{v}}_g)\cdot \text{n}\mathrm{d}S=$$

$${\int }_S({\tau }_{ij}{\text{i}}_j-p{\text{i}}_i)\cdot \text{n}\mathrm{d}S+{\int }_V\rho {\text{gi}}_i\mathrm{d}V+{\int }_V\rho q_i\mathrm{d}V,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_V\alpha \mathrm{d}V+{\int }_S\alpha (\text{v}-{\text{v}}_{\mathrm{g}})\cdot \text{n}\mathrm{d}S={\int }_V{q}_{\alpha }\mathrm{d}V.$$

$$q_{\alpha }=\gamma b(x)\left({\alpha }_{\mathrm{ref}}-\alpha \right),$$

$$q_{\mathrm{i}}=\gamma b(x)(u_{i,\mathrm{ref}}-u_i),$$

$$b(x)=1,$$

$$b(x)=\frac{x-x_{\mathrm{sd}}}{x_{\mathrm{ed}}-x_{\mathrm{sd}}},$$

$$b(x)={\left(\frac{x-x_{\mathrm{sd}}}{x_{\mathrm{ed}}-x_{\mathrm{sd}}}\right)}^2,$$

$$b(x)={\cos }^2\left(\frac{\pi }{2}+\frac{\pi }{2}\frac{x-x_{\mathrm{sd}}}{x_{\mathrm{ed}}-x_{\mathrm{sd}}}\right),$$

$$b(x)=\left(\frac{e^{{\left(\frac{x-x_{\mathrm{sd}}}{x_{\mathrm{ed}}-x_{\mathrm{sd}}}\right)}^2}-1}{e^1-1}\right).$$

$$\gamma \propto \omega ,x_{\mathrm{d}}\propto \lambda ,$$

$$x_{\mathrm{d}}=x_{\mathrm{d},\mathrm{ref}}\cdot \frac{\lambda }{{\lambda }_{\mathrm{ref}}},$$

$$\gamma ={\gamma }_{\mathrm{ref}}\cdot \frac{\omega }{{\omega }_{\mathrm{ref}}},$$

$${\psi }_{tt}=c^2{\psi }_{xx}.$$

$$\psi =\frac{H\omega }{2k}\frac{\sinh (k(z+h))}{\sinh (kh)}{\mathrm{e}}^{i(-\omega t+kx)},$$

$$u={\psi }_z,$$

$$w=-{\psi }_x.$$

$$\mathcal{X}=i\frac{H}{2}\frac{\cosh (k(z+h))}{\sinh (kh)}{\mathrm{e}}^{i(-\omega t+kx)},$$

$$\mathcal{Z}=\frac{H}{2}\frac{\sinh (k(z+h))}{\sinh (kh)}{\mathrm{e}}^{i(-\omega t+kx)}.$$

$${\mathcal{X}}_{tt}=c^2{\mathcal{X}}_{xx}.$$

$${\mathcal{X}}_{tt}=c^2{\mathcal{X}}_{xx}+\underset{q}{\underbrace{\gamma b(x)({\mathcal{X}}_{t,\mathrm{ref}}-{\mathcal{X}}_t)}},$$

$${\psi }_{tt}=c^2{\psi }_{xx}+\gamma b(x)({\psi }_{t,\mathrm{ref}}-{\psi }_t),$$

$$k=\sqrt{\frac{{\omega }^2}{c^2}}=\frac{\omega }{c}.$$

$$k=\sqrt{\frac{{\omega }^2}{c^2}+i\frac{\omega \gamma b(x)}{c^2}}.$$

$$k_j=\sqrt{\frac{{\omega }^2+i\omega \gamma b({\sum }_{n=1}^{j-1}{x}_{{\mathrm{d}}_n}+\frac{1}{2}{x}_{{\mathrm{d}}_j})}{c^2}},$$

$${\psi }_0=\frac{H\omega }{2k_0}\frac{\sinh (k_0(z+h))}{\sinh (k_{0}h)}{\mathrm{e}}^{i(-\omega t)},$$

$$\displaystyle\psi_{j}=\psi_{0}\left(\prod_{n=0}^{j-1}C_{\mathrm{T}_{n}}\right)\cdot\biggl{[}\mathrm{e}^{i\left(\sum_{n=1}^{j-1}k_{n}x_{\mathrm{d}_{n}}+k_{j}\left(x-\sum_{n=1}^{j-1}x_{\mathrm{d}_{n}}\right)\right)}$$

$$\displaystyle-C_{\mathrm{R}_{j}}\mathrm{e}^{i\left(\sum_{n=1}^{j-1}k_{n}x_{\mathrm{d}_{n}}+k_{j}2x_{\mathrm{d}_{j}}-k_{j}\left(x-\sum_{n=1}^{j-1}x_{\mathrm{d}_{n}}\right)\right)}\biggr{]}\quad,$$

$$\displaystyle\psi_{x,j}=ik_{j}\psi_{0}\left(\prod_{n=0}^{j-1}C_{\mathrm{T}_{n}}\right)\cdot\biggl{[}\mathrm{e}^{i\left(\sum_{n=1}^{j-1}k_{n}x_{\mathrm{d}_{n}}+k_{j}\left(x-\sum_{n=1}^{j-1}x_{\mathrm{d}_{n}}\right)\right)}$$

$$\displaystyle+C_{\mathrm{R}_{j}}\mathrm{e}^{i\left(\sum_{n=1}^{j-1}k_{n}x_{\mathrm{d}_{n}}+k_{j}2x_{\mathrm{d}_{j}}-k_{j}\left(x-\sum_{n=1}^{j-1}x_{\mathrm{d}_{n}}\right)\right)}\biggr{]}\quad.$$

$${\left[{\psi }_j={\psi }_{j+1}\right]}_{x={\sum }_{n=1}^j{x}_{{\mathrm{d}}_n}}.$$

Full text

††journal: Ocean Engineering

1 Introduction

In flow simulations, it is usually desired to choose the computational domain as small as possible to reduce the computational effort. However, when simulating free-surface wave propagation, undesired wave reflections at the domain boundaries must be minimized. If this is not achieved, it can lead to large errors in the results. For practical purposes, it is desired to be able to estimate the amount of undesired wave reflection before running the simulations. Of the various techniques for reducing undesired reflections (see Perić and Abdel-Maksoud, 2016), this paper is concerned with the ones that apply source terms to the governing equations in a zone adjacent to the corresponding domain boundaries. Such approaches have been presented under many different names, such as absorbing layers (e.g. Wei et al., 1999), damping zones (e.g. Park et al., 1999; Perić and Abdel-Maksoud, 2016; Jose et al., 2017), dissipation zones (Park et al., 1993), numerical beach (e.g. Clément, 1996), sponge layers (e.g. Israeli and Orszag, 1981; Larsen and Dancy, 1983; Brorsen and Helm-Petersen, 1999; Ha et al. 2013; Choi and Yoon, 2009; Zhang et al., 2014; Hu et al., 2015), Euler overlay method (e.g. Kim et al., 2012; Kim et al., 2013), forcing zones (Perić, 2015; Siemens STAR-CCM+ manual version 11.06), and coupling or relaxation zones (e.g. Jacobsen et al., 2012; Wöckner-Kluwe, 2013; Schmitt and Elsaesser, 2015; Jasak et al., 2015; Vukčević et al., 2016).

The general principle behind all these approaches is that they apply source terms to one, to several, or to all of the governing equations, with the intention of gradually forcing the solution towards some reference solution within a zone (layer) attached to the domain boundary. This damps waves which travel into the layer, but it can also be used to generate waves or to couple different flow solvers (e.g. a viscous solver for the near-field and an inviscid solver for the far-field).

A possible distinction could be that terms like forcing zone and relaxation zone are more general, while others are more specific; for example, the Euler overlay method forces the flow towards the analytical solution of an undisturbed wave. Damping zones, absorbing layers and sponge layers often apply source terms only in a single governing equation, with the forcing term formulated so that it can be interpreted as a damping term. Yet in several cases there seems to be no clear distinction and some of the names are used synonymously.

Thus in the following, the term forcing zone will be used to highlight that the results in this work are applicable to all of the above approaches. In Sect. 2, a general formulation for forcing zones is presented, by which the results in this work can be applied to the other approaches.

The main goal with forcing-zone-type approaches is to achieve reliable minimization of undesired wave reflections. However, the source term function contains tunable parameters which must be adjusted for every simulation (Mani, 2012; Perić and Abdel-Maksoud, 2016). It was argued that these ”parameters and profiles […] can only be determined by trail and error” (Colonius, 2004, p. 337), which at the time of writing is still common practice in flow simulations with ocean waves (Perić and Abdel-Maksoud, 2016). So far, no theoretical foundation has been presented which guides how to choose these parameters and predicts reflection coefficients (Colonius, 2004; Hu, 2008; Modave et al., 2010; Mani, 2012). This is especially problematic since in industrial practice reflection coefficients are seldom evaluated; instead, mostly the default coefficients are used, which can lead to significant errors (Perić and Abdel-Maksoud, 2016). Thus prediction and minimization of undesired wave reflections is a key problem of great practical importance.

The present work aims to solve this problem by presenting an analytical approach to predict the reflection coefficients for forcing-zone-type approaches. The theory can be used to tune the forcing zone parameters before running the simulation. Thus undesired wave reflections can be reliably minimized. The theory is implemented in a computer program, which is made publicly available as free software.222The source code and manual can be downloaded from: https://github.com/wave-absorbing-layers/absorbing-layer-for-free-surface-waves

The theory is derived in Sect. 3 and validated in Sect. 5 using results from finite-volume-based flow simulations; the simulation setup is given in Sect. 4. The theory is shown to predict optimum forcing zone parameters and corresponding reflection coefficients for forcing of $x$-momentum, or $z$-momentum, or $x$- and $z$-momentum combined, or both $x$- and $z$-momentum as well as volume fraction $\alpha$, when forcing to reference values for the undisturbed, calm free surface. The results point out the necessity to adjust the tunable parameters according to the wave parameters. Finally, Appendixes A, B and C discuss common forcing zone implementations, the optimum blending-in of the forcing source terms and the convergence of the analytical solution for the discretized problem towards the analytical solution of the continuous problem.

2 Wave damping using forcing zones

This work considers fluid flow governed by the equation for mass conservation, the three equations for momentum conservation and the equation for the volume fraction, which describes the distribution of the phases:

[TABLE]

[TABLE]

[TABLE]

Here $V$ is the control volume (CV) bounded by the closed surface $\mathrm{S}$, v is the velocity vector of the fluid with the Cartesian components $u_{i}$, $\textbf{v}_{g}$ is the grid velocity, n is the unit vector normal to $S$ and pointing outwards, $t$ is time, $p$ is the pressure, $\rho$ are fluid density, $\tau_{ij}$ are the components of the viscous stress tensor, ij is the unit vector in direction $x_{j}$, with volume fraction $\alpha$ of water. Unless severe wave breaking occurs, the propagation of ocean waves is an approximately inviscid phenomenon. Thus the results in this work apply regardless which formulation for $\tau_{ij}$ is chosen or whether it is neglected altogether.

Undesired wave reflections can be minimized by applying source terms for volume fraction, $q_{\alpha}$, and momentum, $q_{i}$, as

[TABLE]

[TABLE]

with reference volume fraction $\alpha_{\mathrm{ref}}$, reference velocity component $u_{i,\mathrm{ref}}$, forcing strength $\gamma$ and blending function $b(x)$.

The unit of $\gamma$ is $\left[\mathrm{rad/s}\right]$. It regulates the magnitude with which the solution at a given cell is forced against the reference solution. The optimum value for $\gamma$ is case-dependent. Apart from Eq. (5), there exist also source terms which are not directly proportional to the forced quantity; these are discussed in Appendix A.

If the reference solution is the hydrostatic solution for the undisturbed free surface (e.g. $u_{i,\mathrm{ref}}=0$), then the forcing can be interpreted as a ’wave damping’.

The blending term $b(x)$ regulates the distribution of the source term over the domain, where $x$ is the wave propagation direction. Many different types of blending functions can be applied. Figure 1 shows common blending functions, such as constant blending

[TABLE]

linear blending

[TABLE]

quadratic blending

[TABLE]

cosine-square blending

[TABLE]

with start coordinate $x_{\rm sd}$ and end coordinate $x_{\rm ed}$ of the forcing zone. The thickness of the forcing zone is $x_{\mathrm{d}}=|x_{\rm ed}-x_{\rm sd}|$. Though at present it is unknown which blending function would be optimal, generally higher order blending functions are preferred, since they proved more effective in several investigations (e.g. Israeli and Orszag (1981)). Unless mentioned otherwise, an exponential blending is used in this work:

[TABLE]

Perić and Abdel-Maksoud (2016) showed that forcing strength $\gamma$ and forcing zone thickness $x_{\mathrm{d}}$ scale as

$\gamma\propto\omega,\quad x_{\mathrm{d}}\propto\lambda\quad,$

(11)

with angular wave frequency $\omega=\frac{2\pi}{T}$ and wavelength $\lambda$. Thus to achieve the same reflection coefficient and similar free-surface elevations within the forcing zones when running two flow simulations, the parameters $\gamma$ and $x_{\mathrm{d}}$ for the second simulation have to be adjusted as

[TABLE]

[TABLE]

where the $x_{{\rm d,ref}}$, $\lambda_{{\rm ref}}$, $\omega_{{\rm ref}}$ and $\gamma_{{\rm ref}}$ are the corresponding parameters from the first simulation.

Perić and Abdel-Maksoud (2016) demonstrated that the optimum forcing strength $\gamma_{\mathrm{opt}}$ and the corresponding reflection coefficient $C_{\mathrm{R}}$ can be determined with comparatively low computational effort for a given wave period $T$, by running 2D flow simulations for different forcing strength $\gamma$. Via Eqs. (12) and (13) $\gamma_{\mathrm{opt}}$ can then be determined for any wave, for which the same blending $b(x)$ and zone thickness $x_{\mathrm{d}}/\lambda$ relative to the wavelength are used. In this case, $\gamma_{\mathrm{opt}}$ can for practical purposes be considered independent of the wave steepness, independent of the discretization (time step, order and choice of discretization scheme, mesh size), and roughly independent of $x_{\mathrm{d}}$ in the interval $1\lambda\leq x_{\mathrm{d}}\leq 2\lambda$; for example for blending according to Eq. (10), Perić and Abdel-Maksoud (2016) obtained $\gamma_{\mathrm{opt}}\approx\pi\omega$, which gives $C_{\mathrm{R}}<1\%$ for $x_{\mathrm{d}}=2\lambda$. A drawback of this approach is that, for a different blending $b(x)$ or zone thickness (i.e. $x_{\mathrm{d}}<1\lambda$ or $x_{\mathrm{d}}>2\lambda$), the above calibration process has to be carried out again.

3 Theory for predicting wave reflection coefficients

This section considers the propagation of long-crested free-surface waves. As illustrated in Fig. 2, the waves are generated at $x=0$, and travel towards boundary $x=L_{x}$, to which a forcing zone with thickness $x_{\mathrm{d}}$ is attached. The coordinate system origin lies at the level of the undisturbed free surface on the wave generating boundary as illustrated in Fig. 2. The $x$-direction points in the wave propagation direction and the $z$-direction is normal to the undisturbed free surface pointing away from the liquid phase.

Outside the forcing zone, the waves fulfill the one-dimensional wave equation

[TABLE]

with location $x$ in wave propagation direction, velocity stream function $\psi$, and phase velocity $c$.

According to linear wave theory in the complex plane, the velocity stream function for an undamped wave is

[TABLE]

with wave height $H$, angular wave frequency $\omega=2\pi/T$, wave period $T$, wave number $k=2\pi/\lambda$, wavelength $\lambda$, and water depth $h$. Horizontal and vertical velocity components $u$ and $w$ are

[TABLE]

[TABLE]

Horizontal and vertical particle displacements $\mathcal{X}$ and $\mathcal{Z}$ are

[TABLE]

[TABLE]

Thus outside the forcing zone, Eq. (14) is equivalent to

[TABLE]

A detailed discussion of linear wave theory can be found in Clauss et al. (1992).

3.1 Forcing of $x$-momentum

Waves can be damped by attaching a forcing zone to the domain, in which the horizontal particle velocity $u=\mathcal{X}_{t}$ is forced to a reference velocity $u=\mathcal{X}_{t,\mathrm{ref}}$

[TABLE]

with forcing strength $\gamma$ and blending function $b(x)$. The last term in Eq. (19), $q$, is equivalent to the forcing source term in Eq. (5). Since in this work $\mathcal{X}_{t,\mathrm{ref}}=0$, Eq. (19) corresponds to a ’classical’ absorbing layer formulation. Further, Eq. (19) is equivalent to

[TABLE]

with temporal derivative of the reference stream function $\psi_{t,\mathrm{ref}}=0$.

Inserting Eq. (15) into Eq. (14) gives the wave number outside the forcing zone

[TABLE]

Inserting Eq. (15) for piecewise-constant blending $b(x)$ into Eq. (20) gives the wave number inside the forcing zone

[TABLE]

Thus inside the forcing zone, the wave number contains an additional imaginary part which damps the wave amplitude but does not change the wavelength.

The analytical solution to the problem of determining the reflection coefficient for a wave entering a forcing zone according to Eqs. (14) and (20) is obtained by an approach somewhat analogous to the way in which forcing source terms are applied in numerical flow simulations on finite grids. The solution domain is discretized into a finite number of cells as illustrated in Fig. 3. Each cell $j$ corresponds to a forcing zone, within which the stream function is given by $\psi_{j}$, and the complex wave number $k_{j}$ has a constant value. For this, $b(x)$ is evaluated at the cell center:

[TABLE]

with thickness $x_{\mathrm{d},j}$ of zone $j$; $x_{\mathrm{d},j}$ is equivalent to the size of the cell in $x$-direction. Thus the damping is constant within every zone. Reflection and transmission may occur at every interface between two cells.

The benefit of this approach is that even non-continuous blending functions and the influence of the discretization can be considered. With increasing resolution, the theoretical results are expected to converge to the solution of the continuous problem. The latter is not derived here, since for practical purposes only the analytical solution to the discretized problem is of interest. In this manner, the problem remains linear and the solution can be derived as follows.

Consider a wave propagating in positive $x$-direction. The wave is generated at $x=0$ following the coordinate system in Fig. 3. Let the stream function at $x=0$ be

[TABLE]

with wave height $H$, angular wave frequency $\omega$, wave number $k_{0}$, vertical coordinate $z$, and time $t$. Set the transmission coefficient $C_{\mathrm{T}_{0}}=1$ and the reflection coefficient $C_{\mathrm{R}_{0}}=0$, thus the ’inlet’ boundary is perfectly transparent, and waves propagating through it in negative $x$-direction will be fully transmitted without reflection.

For illustration, a domain with $4$ zones is depicted in Fig. 3. Let the wave number $k_{1}=2\pi/\lambda_{1}$ within zone $1$ equal the wave number $k_{0}=2\pi/\lambda_{0}$ of the wave generated at $x=0$, where $\lambda_{0}$ and $\lambda_{1}$ are the corresponding wavelengths. Thus within the first zone $0\leq x\leq x_{\mathrm{d}_{1}}$, there is no wave damping, i.e. $q(x)=0$. At the end of the domain, i.e. at $x=\sum_{n=1}^{4}x_{\mathrm{d}_{n}}$, the boundary is perfectly reflecting (a typical ’wall boundary condition’ in computational fluid dynamics), so the transmission coefficient $C_{\mathrm{T}_{4}}=0$ and the reflection coefficient $C_{\mathrm{R}_{4}}=1$. Within each zone in zones $2$ to $4$, the damping is constant, i.e. $q(x)=q_{j}=-\gamma b(\sum_{n=1}^{j-1}x_{\mathrm{d}_{n}}+\frac{1}{2}x_{\mathrm{d}_{j}})\mathcal{X}_{t}$, yet $q_{2}$ to $q_{4}$ may be of different magnitude, which is indicated through the different shading of the zones in Fig. 3.

By requiring that the particle displacements and velocities must be continuous at every interface between two zones, as they should be at the interfaces between two cells in a flow simulation, the periodic solution is obtained.

For a domain with $j$ zones, the general solution for the velocity stream function $\psi(x)$ within zone $j>0$ can be written as a sum of a right-going (incoming) and a left-going (reflected) wave component

[TABLE]

with $\psi_{0}$ according to Eq. (24). The derivative of $\psi_{j}$ with respect to wave propagation direction $x$ is thus

[TABLE]

Requirement 1: At the interface between zones $j$ and $j+1$, the solution for the zone on the left, i.e. $\psi_{j}$, must equal the solution for the zone on the right, i.e. $\psi_{j+1}$:

[TABLE]

Inserting Eq. (25) into Eq. (27), dividing by

[TABLE]

inserting

[TABLE]

and rearranging for $C_{\mathrm{T}_{j}}$ gives

[TABLE]

Requirement 2: At the interface between zones $j$ and $j+1$, the spatial derivative with respect to $x$ of the solution for the zone on the left, i.e. $\partial\psi_{j}/\partial x=\psi_{x,j}$, must equal the spatial derivative with respect to $x$ of the solution for the zone on the right, i.e. $\partial\psi_{j+1}/\partial x=\psi_{x,j+1}$:

[TABLE]

Inserting Eq. (26) into Eq. (29), dividing by

[TABLE]

inserting $C_{\mathrm{T}_{j}}$, and setting

[TABLE]

gives

[TABLE]

For practical purposes, mainly the ’global’ reflection coefficient $C_{\mathrm{R}}$ is of interest, which is the ratio of the amplitude of the wave, which is reflected back into the solution domain, to the amplitude of the wave, which enters the forcing zone. The global reflection coefficient $C_{\mathrm{R}}$ corresponds to the magnitude of $C_{\mathrm{R}_{j}}$ at the interface to the forcing zone, which depends on all $C_{\mathrm{R}_{j}}$ inside the whole forcing zone. So if the forcing zone starts at zone $1$, then

[TABLE]

where $\mathrm{Re}\{X\}$ and $\mathrm{Im}\{X\}$ denote the real and the imaginary part of the complex number $X$.

Since for the complex exponential functions in this case the order of derivatives and integrals can be exchanged, when setting integration constants to zero one obtains

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Thus $\psi_{j}$, $\psi_{x,j}$ and and all particle displacements and velocities are continuous throughout the domain.

3.2 Forcing of volume fraction $\alpha$, $x$- and $z$-momentum

Combinations of forcing of volume fraction $\alpha$, $x$- and $z$-momentum can be described by the above theory, when $\gamma$ is adjusted as outlined in the following.

Assume that the equations for volume fraction, $x$- and $z$-momentum are coupled, so that forcing of one equation acts on the other equations immediately.

According to linear wave theory, the average energy $\bar{E}=\int_{0}^{\lambda}\int_{-h}^{\eta}\int_{0}^{1\,\mathrm{m}}E\,\mathrm{d}y\mathrm{d}z\mathrm{d}x$ in a regular deep-water free-surface wave can be subdivided as

[TABLE]

where

[TABLE]

with location of the free surface $\eta$ above still water level, water depth $h$, average potential and kinetic energy $\bar{E}_{\mathrm{pot}}$ and $\bar{E}_{\mathrm{kin}}$, and the $x$- and $z$-component $\bar{E}_{\mathrm{kin,x}}$ and $\bar{E}_{\mathrm{kin,z}}$ of the average kinetic energy.

Average potential and kinetic energy in the wave have the same magnitude. Therefore applying forcing with same $\gamma$ to both potential and kinetic energy ($q_{\alpha}$, $q_{\mathrm{x}}$ and $q_{\mathrm{z}}$) produces a forcing of twice the strength as when applying forcing with same $\gamma$ only to the kinetic energy (i.e. $q_{\mathrm{x}}$ and $q_{\mathrm{z}}$); thus $\gamma_{(E_{\mathrm{pot}}+E_{\mathrm{kin}})}=\frac{1}{2}\gamma_{(E_{\mathrm{kin}})}$.

Similarly, since the $x$- and $z$-components of the kinetic energy have on average the same magnitude (i.e. $\bar{E}_{\mathrm{kin,x}}=\bar{E}_{\mathrm{kin,z}}$), applying forcing to only one of these (i.e. either $q_{\mathrm{x}}$ or $q_{\mathrm{z}}$) shifts the optimum of $\gamma$ even further to the right. Therefore, it generally holds that

$\gamma_{(E_{\mathrm{kin,z}})}=\gamma_{(E_{\mathrm{kin,x}})}=2\gamma_{(E_{\mathrm{kin}})}=2\gamma_{(E_{\mathrm{pot}})}=4\gamma_{(E_{\mathrm{pot}}+E_{\mathrm{kin}})}\quad,$

(35)

with forcing strengths for the cases of forcing of $w$-velocity, $\gamma_{(E_{\mathrm{kin,z}})}$, forcing of $u$-velocity, $\gamma_{(E_{\mathrm{kin,x}})}$, forcing of both $u$- and $w$-velocities, $\gamma_{(E_{\mathrm{kin}})}$, as well as forcing of $u$- and $w$-velocities and volume fraction $\alpha$, $\gamma_{(E_{\mathrm{pot}}+E_{\mathrm{kin}})}$.

Therefore, the theory from Sect. 3 can be used for the other forcing source terms as well, given that $\gamma$ is adjusted according to Eq. (35); thus when simultaneously forcing $\alpha$, $x$- and $z$-momentum, each with forcing strength $\gamma=20\,\mathrm{rad/s}$, then the theoretical prediction for the reflection coefficient $C_{\mathrm{R}}$ is the same as when using the theory from Sect. 3.1 with $\gamma=80\,\mathrm{rad/s}$.

In the following sections, the theoretically predicted reflection coefficients will be compared with those obtained from flow simulations based on the solution of the Navier-Stokes equations.

4 Simulation setup

In the flow simulations in this work, a regular long-crested free-surface wave train is created and propagates in positive $x$-direction towards a forcing zone as sketched in Fig. 2, where it is partly reflected and partly absorbed. The wave has height $H=0.16\,\mathrm{m}$, period $T=1.6\,\mathrm{s}$, wavelength $\lambda\approx 4\,\mathrm{m}$ and is moderately non-linear (steepness $H/\lambda$ is $\approx 29\%$ of maximum steepness). Deep water conditions apply ($h/\lambda>0.5$).

The commercial flow solver STAR-CCM+ version 11.06.010-R8 from Siemens (formerly CD-adapco) is used for the simulations. The governing Eqs. are (1) to (3). The volume of fluid (VOF) method is used to account for the two fluid phases, liquid water and gaseous air, using the High Resolution Interface Capturing scheme (HRIC) as given in Muzaferija and Perić (1999). The governing equations are applied to each cell and discretized according to the finite volume method. All integrals are approximated by the midpoint rule. The interpolation of variables from cell center to face center and the numerical differentiation are performed using linear shape functions, leading to approximations of second order. The integration in time is based on assumed quadratic variation of variables in time, which is also a second-order approximation. Each algebraic equation contains the unknown value from the cell center and the centers of all neighboring cells with which it shares common faces. The resulting coupled equation system is then linearized and solved by the iterative STAR-CCM+ implicit unsteady segregated solver, using an algebraic multigrid method with Gauss-Seidel relaxation scheme, V-cycles for pressure and volume fraction of water, and flexible cycles for velocity calculation. The under-relaxation factor is $0.9$ for velocities and volume fraction and $0.4$ for pressure. For each time step, eight iterations are performed; one iteration consists of solving the governing equations for the velocity components, the pressure-correction equation (using the SIMPLE method for collocated grids to obtain the pressure values and to correct the velocities) and the transport equation for the volume fraction of water. Further information on the discretization of and solvers for the governing equations can be found in Ferziger and Perić (2002) or the STAR-CCM+ software manual.

The forcing approaches from Sect. 2 are used to minimize wave reflections. Simulations are either run with forcing of $x$-momentum ($q_{\mathrm{x}}$), of $z$-momentum ($q_{\mathrm{z}}$), of both $x$- and $z$-momentum ($q_{\mathrm{x}}$, $q_{\mathrm{z}}$), or of volume fraction $\alpha$ and $x$- and $z$-momentum ($q_{\alpha}$, $q_{\mathrm{x}}$, $q_{\mathrm{z}}$). The forcing zone has exponential blending according to Eq. (10). Forcing zone thicknesses $x_{\mathrm{d}}$ in the range between $0.25\lambda$ and $2\lambda$ are investigated.

The simulations are run quasi-2D, i.e. with only one layer of cells in $y$-direction and symmetry boundary conditions applied to the $y$-normal boundaries. The domain has dimensions $0\leq x\leq L_{x}$ and $-h\leq z\leq 0.5\lambda$, with domain length $L_{x}=6\lambda$ and water depth $h=4.5\lambda$.

The wave is generated by prescribing the volume fraction and velocities according to Fenton’s (1985) $5^{\mathrm{th}}$-order Stokes theory at the inlet boundary ($x=0$). Due to the Stokes drift of the fluid particles in the wave, the inlet produces a net mass flux into the domain, which acts to raise the mean water level; by prescribing the hydrostatic pressure at either the boundary to which the forcing zone is attached ($x=L_{x}$) or the bottom boundary ($z=-h$), such an undesired accumulation of mass can be avoided. Perić and Abdel-Maksoud (2015) found that the first option can lead to low-frequency fluctuations in the amount of mass within the domain, since the forcing zone delays how the pressure boundary regulates the amount of mass within the domain; prescribing the pressure at the domain bottom resolved this problem. In practice however, it is far more common to prescribe the pressure at the boundary to which the forcing zone is attached. For this reason, all cases investigated in this work were simulated with pressure prescribed at boundary $x=L_{x}$.

Additionally, the simulations for $q_{\mathrm{x}}$ and $q_{\mathrm{z}}$ with $x_{\mathrm{d}}=2\lambda$ were rerun with pressure prescribed at boundary $z=-h$ and boundary $x=L_{x}$ set to no-slip wall. The results confirmed that the choice for the pressure boundary does not significantly influence the reflection coefficients (on average $<1\%$ difference). Therefore in practice, the theory presented in this work can be used for both boundary choices. Yet from an academical point of view, there is a difference: as illustrated in Fig. 4, prescribing the pressure at $x=L_{x}$ produces a node at that boundary; when the boundary condition is set to wall instead, there occurs a maximum amplification point. The theory in Sect. 3 is derived for the latter case. Thus in Sect. 5, comparisons of surface elevation with theory (Figs. 8 and 18) are given for the simulations with pressure prescribed at the domain bottom; all other results are given for pressure prescribed at $x=L_{x}$, because this approach is primarily used in practice and the observed difference in results between the two approaches were insignificant.

All remaining boundaries are no-slip walls. For further details on boundary conditions, see Ferziger and Perić (2002).

The volume fraction and velocities in the solution domain are initialized according to Fenton’s (1985) Stokes $5^{\mathrm{th}}$-order theory to reduce the simulation time. The domain is discretized using a rectilinear grid with local mesh refinement around the free surface. The free surface stays at all times within the region of the finest mesh with $100$ cells per wavelength and $16$ cells per wave height. The computational grid, which consists of $\approx 43\,000$ cells, is shown in Fig. 5.

The total simulated time is $18\,\mathrm{s}\approx 11.3T$ with a time step of $\Delta t=T/1000$.

To calculate the reflection coefficient, the free-surface elevation is written to a file at $80$ evenly spaced time intervals per wave period. From the elevation in an interval of $1.5\lambda$ adjacent to the forcing zone, the overall highest and lowest wave heights are obtained as $H_{\mathrm{max}}$ and $H_{\mathrm{min}}$. From these, the reflection coefficient is calculated as in experiments after the procedure by Ursell et al. (1960) as

[TABLE]

This simple approach contains only a small background noise ($\approx 1\%$, see Fig. 5) due to the finite number of time intervals, which was considered accurate enough for the present purposes.

The reflection coefficient gives the ratio of the the wave heights of reflected to incoming wave; the ratio of the energy of the reflected ($E_{\mathrm{refl}}$) wave to the incoming ($E_{\mathrm{in}}$) wave is $E_{\mathrm{refl}}/E_{\mathrm{in}}\approx C_{\mathrm{R}}^{2}$, since according to linear wave theory the wave energy depends on the wave height squared. Thus a reflection coefficient of $C_{\mathrm{R}}=10\%$ means that the forcing zone reflects $1\%$ of the incoming wave energy.

5 Results and discussion

In this section, reflection coefficients obtained in flow simulations using different forcing approaches and parameters are compared to theoretical predictions according to Sect. 3.

5.1 Forcing of $x$-momentum

Flow simulations are carried out with the setup from Sect. 4 and compared to the theory presented in Sect. 3. Forcing of $x$-momentum according to Eq. (5) with exponential blending from Eq. (10) is used to damp waves with period $T=1.6\,\mathrm{s}$ and height $H=0.16\,\mathrm{m}$. Simulations are run for different forcing strengths ($0.625\,\mathrm{rad/s}\leq\gamma\leq 5120\,\mathrm{rad/s}$) and forcing zone thicknesses ($0.25\lambda\leq x_{\mathrm{d}}\leq 2\lambda$).

Comparing simulation results and theory in Fig. 6 shows that the theory can predict reflection coefficient $C_{\mathrm{R}}$ with high accuracy. Although the waves are moderately nonlinear ($\approx 29\%$ of maximum steepness), the deviations are small; this agrees with results from Perić and Maksoud (2016), where deep-water waves up to $\approx 71\%$ of the maximum steepness were investigated, and the influence of wave steepness on the reflection coefficient in this range was found to be negligible for practical purposes ($<2\%$ difference), except for $\gamma$ roughly an order of magnitude below optimum, where the damping improved noticeably for the steeper wave.

The close-up on the optimum forcing regime in Fig. 7 shows that simulation results were not able to reproduce reflection coefficients below $\approx 1\%$, see e.g. the curve for $x_{\mathrm{d}}=2\lambda$ at $\gamma\approx 20\,\mathrm{rad/s}$ in Fig. 7; this was attributed to a background noise in the approach for calculating $C_{\mathrm{R}}$, due to the finite number of time intervals used; possibly also the interpolation of the surface elevation on the finite grid limits the detection of $C_{\mathrm{R}}<1\%$.

Further, while for thinner zones ($x_{\mathrm{d}}<1.5\lambda$) there is a single optimum for forcing strength $\gamma$, for thicker zones ($x_{\mathrm{d}}=2\lambda$) there can be more than one local optimum. This may have implications for the choice of $\gamma$ when damping irregular waves, which future research may show.

The simulated and theoretical surface elevations agree well as shown in Figs. 8 and 9. The peaks in the partial standing wave have different locations depending on $\gamma$, which shows that, with increasing forcing strength $\gamma$, the effective reflection location shifts from the boundary, to which the forcing zone is attached (here $x=24\,\mathrm{m}$), towards the entrance to the forcing zone (here $x=16\,\mathrm{m}$); this underlines the importance of including reflections which occur within the zone in the analysis. Just as the surface elevations, the theory predictions for horizontal and vertical velocities agree well with simulation results; since there are no significant differences, these are not plotted here.

Figure 10 shows that the forcing produces vorticity (here: $w_{x,j}-u_{z,j}$) in the flow field within the forcing zone. Although the generated amount of vorticity is comparatively small, it is clearly visible. Further, small amounts of vorticity are also generated at the free surface and at locations where the mesh size changes, as artifacts of the interface capturing scheme and the discretization; however, such a vorticity generation is discretization dependent, i.e. it disappears on infinitely fine grids, whereas the forcing-generated vorticity turned out to be grid-independent.

As Figs. 11, 12 and 13 show, the theory predicts the forcing-based vorticity generation remarkably well.

While the velocities and particle displacements are continuous everywhere within the domain according to the theory given in Sect. 3, evaluating vorticity $w_{x,j}-u_{z,j}$ for the theory shows that the vorticity is continuous within each zone, as well as at each interface between two zones where the blending function $b(x)$ has the same value. When $b(x)$ changes between two adjacent zones, then the vorticity must be discontinuous at the interface between these zones.

To show that this is realistic and occurs in the flow simulations as well, the simulations are rerun using a strongly discontinuous blending function

[TABLE]

with forcing zone thickness $x_{\mathrm{d}}=24\,\mathrm{m}-16\,\mathrm{m}=8\,\mathrm{m}\approx 2\lambda$; $b(x)$ is illustrated in Fig. 14.

Figures 15 and 16 show that also in the flow simulations the vorticity is discontinuous between cells where $b(x)$ changes, and is continuous otherwise; the agreement between theory and simulation was again satisfactory, as is exemplarily shown for $\gamma=10\,\mathrm{rad/s}$ in Fig. 16.

5.2 Forcing of $z$-momentum

Repeating the simulations from Sect. 5.1 with forcing of $z$-momentum instead of $x$-momentum shows that the theory predicts optimum forcing strength and corresponding reflection coefficient $C_{\mathrm{R}}$ reliably. Comparing Figs. 17 and 18 with Figs. 7 to 9 shows that, for optimum or lower values of forcing strength $\gamma$, the results for forcing of $x$- and $z$-momentum both agree well with theory predictions.

For stronger forcing than optimum (here: $\gamma>40\,\mathrm{rad/s}$) the theory predictions of $C_{\mathrm{R}}$ are conservative. Figure 18 shows that for such $\gamma$-values the mean surface elevation increases within the forcing zone. This becomes more pronounced with increasing $\gamma$. Thus stronger-than-optimum forcing of $z$-momentum leads to a noticeable net mass flux into the forcing zone in the simulation. In the present case, this mass flux resulted in a lower reflection coefficient $C_{\mathrm{R}}$ compared with $x$-momentum forcing for the same $\gamma$.

However, since this effect is negligible for optimum and lower forcing strength, the theoretical prediction is satisfactory for practical purposes.

5.3 Forcing of $x$- and $z$-momentum

Repeating the simulations from Sect. 5.1 with forcing of both $x$- and $z$-momentum shows that the theory, with the extension given in Sect. 3.2, reliably predicts the optimum forcing strength $\gamma$ and the corresponding reflection coefficient $C_{\mathrm{R}}$. As expected from the results in Sects. 5.1 and 5.2, Fig. 19 shows that for stronger than optimum forcing the theory overpredicts the reflection coefficients.

5.4 Forcing of volume fraction $\alpha$ and $x$- and $z$-momentum

Repeating the simulations from Sect. 5.1 with forcing of volume fraction $\alpha$ and both $x$- and $z$-momentum shows that the theory, with the extension given in Sect. 3.2, reliably predicts the optimum forcing strength and the corresponding reflection coefficient. Figure 20 shows that for stronger than optimum damping the theory overpredicts the reflection coefficients, otherwise the results agree well.

6 Conclusion

To reliably damp surface waves using forcing-zone-type approaches (such as absorbing layers, sponge layers, relaxation zones, etc.), the case-dependent parameters in the source functions must be adjusted to the wave. These adjustments are not ’ad-hoc’. Instead, there is a mathematical foundation for choosing the tunable parameters. The parameter forcing strength $\gamma$ scales with the angular wave frequency $\omega$, and the layer thickness $x_{\mathrm{d}}$ scales with wavelength $\lambda$.

In this work, a theory was presented which reliably predicts the reflection coefficients for forcing zones with sufficient accuracy for practical purposes. The theory applies for forcing of horizontal fluid velocities, forcing of vertical fluid velocities, forcing of the volume fraction, as well as for combinations of these approaches. The theory is given for a general forcing source term formulation, so that it can be easily applied to many existing forcing-zone-type approaches from literature and existing implementations in commercial flow solvers. A computer code for evaluating the theory has been made available to the scientific community as free software.

The main benefit of the theory is that it can be used to optimally select the forcing zone parameters before running the simulation. The present results illustrate the importance of adjusting these parameters for every flow simulation, since the common practice to use default settings can produce large errors.

The theory predictions for forcing of horizontal velocities were remarkably accurate both for reflection coefficients and the flow within the layer. This also held for forcing of vertical velocities and volume fraction, except for the case when the forcing strength was larger than optimum; then, the reflection was lower than predicted by the theory. This was attributed to a net mass flux into the forcing layer, which the theory does not predict. However, for practical purposes the theory works satisfactory for all forcing types discussed, since the optimum parameters and the corresponding reflection coefficients are predicted reliably, and otherwise the predictions are conservative.

Apart from predicting the reflection before the simulation is run, the theory can also be used to assess the amount of undesired wave reflections after the simulation is run. This is especially useful for classification societies, when they have to assess the quality of results from flow simulations that they did not perform themselves. Since it is not industry standard to provide reflection coefficients for flow simulations, it is beneficial to be able to quickly take any forcing zone formulation that can be expressed in terms of Eqs. (5) and (4), obtain the reflection coefficients from theory, and judge whether the settings were appropriate.

Future research is necessary to verify how accurate the theory predicts the damping of three-dimensional waves, irregular waves and highly non-linear waves, such as rogue waves, waves close to breaking steepness and even breaking waves; literature results are promising that the theory covers these cases as well (Perić et al., 2015; Perić and Abdel-Maksoud, 2016). Further, research is necessary regarding the optimum choice of blending function. Moreover, applying and possibly extending the theory to the simultaneous generation and absorption of waves via forcing zones, and subsequently the coupling of different flow solvers (e.g. viscous and inviscid codes), seems promising in the light of the present findings.

Appendix A Forcing proportional to the velocity squared

Examples of widely used wave damping implementations are the approaches by Choi and Yoon (2009, implemented in STAR-CCM+ by Siemens) and by Park et al. (1999, implemented in ANSYS Fluent), which are

[TABLE]

[TABLE]

with vertical velocity $w$, forcing strengths $\gamma$ and $\tilde{\gamma}$, vertical coordinate $z$ with domain bottom at $z_{{\rm b}}$, and vertical location $z_{{\rm fs}}$ of the free surface. Further $b(x)$ is exponential (Eq. (10)) in Eq. (38) and quadratic (Eq. (8)) in Eq. (39). In Eq. (38), the first term corresponds to Eq. (5) and is therefore directly proportional to the vertical velocity, while the second term contains an additional factor $|w|$, which renders this forcing term directly proportional to $w|w|$. Equation (39) corresponds to the second term in Eq. (38), except for a factor $0.5$, an additional vertical blending $\left(1-\frac{z-z_{{\rm fs}}}{z_{{\rm b}}-z_{{\rm fs}}}\right)$ and a slightly different $b(x)$ (see Fig. 1).

At the time of writing, the default values in the commercial codes for forcing strengths in Eqs. (38) and (39) are $\gamma=10.0\,\mathrm{rad/s}$ and $\tilde{\gamma}=10.0\,\mathrm{m^{-1}}$. Perić and Abdel-Maksoud (2016) found that for a $w|w|$-proportional forcing, as in the second term in Eq. (38) and in Eq. (39), the optimum value for $\tilde{\gamma}$ is more than one order of magnitude larger than the optimum value for $\gamma$. Thus with default settings in STAR-CCM+, the second term in Eq. (38) has a negligible effect compared to the first term.

Both $w$-proportional and $w|w|$-proportional forcing produced comparable reflection coefficients at optimum settings, so both approaches can be used to damp waves successfully. However, for a fixed forcing strength, directly-proportional forcing as in Eq. (5) has a wider range of wave frequencies which are damped satisfactorily, as illustrated in Fig. 21. Historically, $w|w|$-proportional forcing terms may have been introduced as analogy to porous media flows, where for larger flow rates effects like turbulence lead to nonlinearities which can be expressed as quadratically dependent on the flow velocity, such as the Forchheimer or Brinkman extension to Darcy’s law, see Straughan (2008). However, this analogy is not entirely valid. Even in steep nonlinear ocean waves, turbulence effects are insignificant unless there is wave breaking, which especially with regard to coupling of different flow solvers should not be provoked inside the forcing zone. Moreover, Fourier approximation methods allow to split nonlinear waves into different regular harmonics, so applying a forcing directly proportional to the velocity according to Eq. (5) acts on the higher harmonics as well, so the damping of the higher harmonics is already accounted for in $w$-proportional forcing.

Appendix B Theory implications for tuning forcing strength $\gamma$ and for the choice of blending function $b(x)$

Since the theory from Sect. 3 was validated in Sect. 5 and can be considered sufficiently accurate for practical purposes, the following sections show theoretical predictions without backup from simulation results. This is mainly because, to obtain the following results, a very large number of simulations would be required, the combined computational effort being out of the scope of this study.

For a typical forcing zone setup, Fig. 22 illustrates the necessity of adjusting the forcing strength $\gamma$ for different waves.

So far, the optimum choice for the blending term $b(x)$ for a given zone thickness $x_{\mathrm{d}}$ is not known. However, Fig. 23 shows that common choices for blending functions such as quadratic blending (Eq. (8)), cosine-square blending (Eq. (9)) and exponential blending according to Eq. (10) perform similarly well, with perhaps the exponential blending performing slightly better. Thus these blending functions can all be used successfully.

Constant blending (Eq. (6)) and linear blending (Eq. (7)) are not recommended, since they lead to considerably higher reflection coefficients.

Appendix C Convergence of theory to solution for continuous blending

The theory from Sect. 3 subdivides the forcing zone into $n$ zones with constant blending $b(x)$ as illustrated in Fig. 24. This section demonstrates that, if $n$ is larger than a certain threshold, then the theory results can be considered independent of $n$. Thus also the wave damping in flow simulations is basically grid-independent, if the number $n$ of grid cells, by which the forcing zone is discretized in wave propagation direction, is above the same threshold.

Figure 25 shows for a subdivision into $32$ zones, that the results are barely distinguishable from subdivisions into larger numbers of zones. This agrees well with findings by Perić and Abdel-Maksoud (2016), where the wave damping was observed to be grid-independent for practical flow simulation setups (i.e. grids with at least $30$ cells per wavelength).

Figure 26 shows that for subdivision into less than $16$ zones, the results differ significantly from the results for $>32$ zones. This is relevant when assessing flow simulations, in which a combination of grid stretching and forcing zones is used to damp the waves; in industrial practice, these two wave damping approaches are sometimes combined with the intention to lower the computational effort and to improve the damping. However, Figs. 25 and 26 show that, if due to the grid stretching the number of grid cells per zone thickness drops below a certain threshold, then grid stretching can significantly increase reflection coefficient $C_{\mathrm{R}}$. Based on the present results, it is recommended to have cell sizes of at least $\lambda/10$ when combining grid stretching and forcing zones.

If the forcing zone is subdivided into a sufficient number of zones $n$, then the difference between the theory solutions for different $n$ can be estimated by a Richardson-type extrapolation. Detailed information on Richardson extrapolation can be found e.g. in Richardson (1911), Richardson and Gaunt (1927), and Ferziger and Perić (2002). Say

[TABLE]

where $A$ is the analytical solution for $n=\infty$, $A_{\mathrm{h}}$ is the analytical solution for $n=x_{\mathrm{d}}/h$, and $\varepsilon_{\mathrm{h}}$ is the error. Let all zones have the same thickness $x_{\mathrm{d},j}=h$, with a total forcing zone thickness of $x_{\mathrm{d}}=\sum_{j=1}^{n}x_{\mathrm{d},j}$. Thus when using zones of twice the thickness, i.e. $x_{\mathrm{d},j}=2h$, one obtains

[TABLE]

and similar for further refinement or coarsening.

Taylor-series analysis of truncation errors suggests that the error $\varepsilon_{h}$ is proportional to some power $p$ of the zone thickness $h$, i.e. $\varepsilon_{h}\propto h^{p}$. It follows that the error with a twice coarser spacing is

[TABLE]

where $p$ is the order of convergence. Setting Eqs. (40) and (41) as equal and inserting Eq. (42) leads to

[TABLE]

Insert Eq. (42) into Eq. (43) written for $\varepsilon_{\mathrm{2h}}$ finally gives

[TABLE]

Figures 25 and 26 show that the deviation of $C_{\mathrm{R}}$ can differ depending on $\gamma$. Thus to estimate $\varepsilon_{\mathrm{h}}$ and $p$ in Eqs. (43) and (44), set

[TABLE]

[TABLE]

where $C_{\mathrm{R,h}}(\gamma)$ is the reflection coefficient for zone thickness $x_{\mathrm{d},j}=h$ and forcing strength $\gamma$, and the $\max\{X(\gamma)\}$-function delivers the maximum value of $X$ of all values $\gamma$ in the considered range $0.625\,\mathrm{rad/s}\leq\gamma\leq 10\,000\,\mathrm{rad/s}$.

Figure 27 shows, exemplarily for $b(x)$ according to Eq. (10), that the analytical solution converges with $2^{\mathrm{nd}}$ order to the analytical solution for $n=\infty$. For the blending functions investigated in Appendix B, all curves showed $2^{\mathrm{nd}}$ order convergence (i.e. $p\rightarrow 2$ if $n\rightarrow\infty$), except constant blending ($b(x)=1$), for which the solution naturally must be exact independent of the number of zones. It is out of the scope of this work to rigorously prove that for all possible $b(x)$ the order of convergence will be at least $p=2$, so this is left for future research; however, the present results suggest that this is the case.

Figure 28 shows that the error estimate $\varepsilon_{\mathrm{h}}$ for different $n$ decays accordingly, and the difference between greatest discrepancies of reflection coefficients for zone thicknesses $h$ and $0.5h$ shows the same rate of decay and lies below $\varepsilon_{\mathrm{h}}$. Thus for practical grids in flow simulations, as well as for the theory results plotted in this work, the forcing zone performance can be assumed independent of the number $n$ of zones or grid cells.

Given the simulation setup in Sect. 4, it is expected that, when the grid resolution increases, the wave damping behavior in the flow simulation will converge towards the solution for the specified continuous blending function $b(x)$. Since Sect. 5 showed that the theory from Sect. 3 predicts flow simulation results with great accuracy, it is expected that the results of the theory from Sect. 3, which is based on discontinuous piece-wise constant blending $b(x)$, will converge towards the analytical solution for any given continuous blending function $b(x)$, if the forcing zone is subdivided into a sufficient number of zones $n$.

Acknowledgements

The study was supported by the Deutsche Forschungsgemeinschaft (DFG).