Application of Non-Reflective Boundary Conditions in Three-Dimensional Numerical Simulations of Free-Surface Flow Problems

Abstract

This paper considers the issue of using non-reflective boundaries for surface wave simulations within the framework of three-dimensional Navier–Stokes equations. We distinguish a wave damping approach among the known implementations of non-reflective boundary conditions in surface wave simulations. The approach employs a sponge boundary layer to dampen incident waves. In this paper, we describe implementations of this approach on unstructured meshes. We also present the calibration of free parameters, the values of which control wave damping efficiency and the amplitude of reflected waves. Comparison of the results obtained at different types of distribution for the free parameter was conducted. The implemented wave damping approach was tested by simulating a solitary wave propagating in a water tank. We demonstrate the use of damping non-reflective boundary conditions for the case of a wave traveling across the surface of a real body of water near Sakhalin Island while considering its bathymetry.

1. Introduction

Simulations of waves traveling across open-boundary water areas, where the waves are allowed to freely escape from the computational domain, often require the use of non-reflective boundary conditions when modeling the infinite space beyond the domain.

There are several methods for organizing non-reflective boundary conditions for waves. The absence of a single method is explained by the fact that none of them ensures the free exit of the waves from the simulation area with no wave reflections. The first group of methods is generalized by the name “non-reflecting boundary conditions”, which implies the use of decomposition of the solution of the equation in a Fourier series, extrapolation of the equations to asymptotic solutions, or implementation of the method of characteristics [1,2,3,4,5,6] for setting non-reflecting boundary conditions. The last-mentioned method of characteristics is based on the analytical derivation of non-reflecting boundary conditions with Riemann invariants [3,4,6] and is often used in numerical simulations of waves with the help of shallow water equations [7]. In the case of 3D simulations, the method is often used in gas dynamics problems [3,4,5,6]. Implementation of the method of characteristics for 3D simulations of surface waves is rare in the literature. One of the works where an example realization of a nonreflecting boundary condition algorithm for 3D Navier–Stokes equations and the VOF method visibly demonstrated the efficiency of this procedure for structured computational meshes.

A general drawback of the first group of methods is their limited applicability to linear waves; the existence of non-linear waves in the tanks being simulated can result in the growth of the general mass of water in the computation domain with time [1]. The expressions for boundary conditions in the first group of methods can be not “local” and include integral parameters that result in problems during program realization in parallel computation.

The second group of methods, generalized as the wave damping methods, are considered to be relatively simple methods for implementation. They are based on the idea of wave energy forced rejection.The method is known in the literature as a ”sponge layer”, an “absorbing layer”, or a “damping zone” [2,8]. It is popular in the field of hydrodynamic computations due to its usability and because it does not depend on the computational mesh used. Suppression of waves is realized through introduction of an additional source of forces in the momentum conservation equation distributed along the marked area of wave suppression. A disadvantage of the method is the fact that a part of the computation domain near a non-reflected boundary condition is used for damping, which increases the sizes of the simulation domains and the number of computational elements.

Research on this method, especially in terms of the influence of model parameters on the simulation results, has been described in some recent works [2,8,9,10]. The author of [2] showed that, when damping force has the same value, preference should be given to a wide zone of damping, as a bigger zone results in a smoother attenuation of the wave burst. The authors of [11] substantiated the fact that the value of the damping force coefficient is directly proportional to the depth of the water in the computation tank. They practically note that it is possible to produce the optimal values of mathematical parameters for this method with only test computations [2,8,9,10,11]. The authors of [12] studied the reasonability of accounting for the sponge layer source term in the volume fraction transport equation and compared the results obtained at different forms of the distribution of damping force in the vertical direction. It was shown that afore-mentioned modifications do not provide significant advantages in problems involving a free surface. Complex research on the method employed by the authors of [10] shows that the greatest contribution to the damping of waves is stipulated by the width of the damping zone as compared to the computational mesh refinement.

It is also important to note a known technique for organizing non-reflecting boundary conditions that does not require mathematical manipulation or modifications of the initial system of equations. This technique implies artificial dissipation of the perturbation due to exponential “elongation” of the cells of the mesh when approximating outer boundaries of the computation domain [2,4,10,13]. The technique works well in the case of flat flows though causes problems in mesh generation in the case of 3D configurations. There are also some examples where a combination of several methods makes it possible to obtain productive results [2,9].

This paper describes the implementation of a wave damping method for the mathematical model used in the LOGOS software package to simulate wave propagation problems. A formula of the source term for the momentum conservation equation is provided; the procedure for calibrating setup variables to obtain the least amplitude of the reflected waves is also described in detail. Additionally, the paper presents a mathematical model based on a system of three-dimensional Navier–Stokes equations and the VOF method. The implemented algorithm is used to calculate the problem of a solitary wave traveling across a water tank. The use of damping non-reflective boundary conditions is demonstrated in the case of a water wave traveling across the surface of a real body near Sakhalin Island while considering its bathymetry.

2. Governing Equations

2.1. Mathematical Model

We considered a mathematical model for traveling wave simulations based on a system of Navier–Stokes equations in combination with the VOF method [14] implemented in the domestic software package LOGOS (version 5.4) [15,16,17].

This software package underwent large-scale verification on problems such as those involving a free surface [18] and tsunami-related problems [19].

We assume that the flow is isothermal and that all phases have the same velocity field. Considering such assumptions, we obtained a system of equations consisting of equations of conservation of mass and momentum and a volume fraction transfer equation. In Cartesian coordinates, the system has the following form:

(1)$\left\{\begin{array}{l}\frac{\partial u_i}{\partial x_i}=0,\\ \mathsf{\rho }\frac{\partial u_i}{\partial t}+\mathsf{\rho }\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}{\mathsf{\tau }}_{ij}+\mathsf{\rho }{g}_i\\ \frac{\partial {\alpha }_w}{\partial t}+\frac{\partial }{\partial x_i}\left(u_i{\alpha }_w\right)=0,\end{array}\right.,$

where i, j are the subscripts indicating that the vector components belong to the Cartesian coordinates, i, j = {x, y, z}, ρ is the mean density calculated as ρ = (ρ_wα_w + ρ_aα_a), w (water) is the subscript indicating the quantities belonging to the “water” phase, a (air) is the subscript indicating the quantities belonging to the “air” phase, α_w is the volume fraction of water, u_i is the component of the velocity vector i = {x, y, z}, t is time, p is pressure, x_i is the component of the vector of the Cartesian coordinates i = {x, y, z}, and τ_ij is the tensor of viscous stresses which, according to the Boussinesq hypothesis, takes the form

(2)${\mathsf{\tau }}_{ij}=\mathsf{\mu }\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\mathsf{\delta }}_{ij}\right),$

where μ is the dynamic viscosity, δ_ij is the Kronecker delta, and g_i is the component of the gravitational acceleration vector.

Gravity is included using an algorithm based on bulk force correction [20] to prevent spurious oscillations related to the non-colocated position of unknowns on grids with an arbitrary type of cells.

The system of equations is solved by numerical integration on a finite-volume grid. The equations are discretized using an original fully implicit method for solving the Navier–Stokes equations in simulations of free-surface multi-phase flows [17,21].

Before discretizing the equations in (1), it makes sense to use transformations to improve the accuracy and stability of the solution. The momentum equation is written in its half-divergent form because, as shown in [22], this representation compensates for the approximation errors associated with the imperfect fulfillment of the mass balance condition in the cell and resolves the shape of the free surface more accurately:

(3)$\mathsf{\rho }\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_i{u}_j\mathsf{\rho }\right)-u_i\frac{\partial }{\partial x_j}\left(u_j\mathsf{\rho }\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}{\mathsf{\tau }}_{ij}+\mathsf{\rho }{g}_i.$

Thus, the system of Equation (1) has the following ultimate form:

(4)$\left\{\begin{array}{l}\frac{\partial u_i}{\partial x_i}=0,\\ \mathsf{\rho }\frac{\partial u_i}{\partial t}+\frac{\partial }{\partial x_j}\left(u_i{u}_j\mathsf{\rho }\right)-u_i\frac{\partial }{\partial x_j}\left(u_j\mathsf{\rho }\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}{\mathsf{\tau }}_{ij}+\mathsf{\rho }{g}_i,\\ \frac{\partial {\alpha }_w}{\partial t}+\frac{\partial }{\partial x_i}\left(u_i{\alpha }_w\right)=0\end{array}\right.$

The system of Equation (4) must be supplemented by boundary conditions. For tsunami wave simulations, one generally uses “wall”-type boundary conditions for the bottom and non-reflective boundary conditions for the outer boundaries of water bodies. There is no need to introduce any additional boundary conditions for runup simulations; one uses the “wall” boundary condition only for the bottom surface. On solid walls, the gradients of pressure and volume fraction are zero:

(5)$\frac{\partial p}{\partial n}=0,\frac{\partial {\alpha }_k}{\partial n}=0,$

the value of velocity is zero:

u = 0, v = 0, w = 0,(6)

Thus, the fluid can neither slide along nor flow normal to the fluid/wall interface. Velocity and shear stresses at the fluid/fluid interface must be continuous. Air at the upper boundary has zero static pressure, and the gradients of velocity and volume fractions are zero:

(7)$\frac{\partial u}{\partial n}=0,\frac{\partial v}{\partial n}=0,\frac{\partial w}{\partial n}=0,\frac{\partial {\alpha }_k}{\partial n}=0$

Equation (4) is discretized by finite volumes on an arbitrary unstructured mesh and solved numerically by a fully implicit method [17,21] based on the known SIMPLE algorithm. Simulations of free-surface flows involve certain modifications of the SIMPLE algorithm. The basic formulas of the SIMPLE algorithm, its boundary conditions, and implementation in the LOGOS software package are described in detail in [17,19,21]. Following is a summary of the equation discretization procedure given in as much detail as required for the presentation below.

Consider the finite-volume discretization of the equations we use in this work. The basic equation in solving system (4) is the transfer equation of the passive scalar φ:

(8)$\frac{\partial \mathsf{\rho }\mathsf{\phi }}{\partial t}+\frac{\partial }{\partial x_i}\left(\mathsf{\rho }\mathsf{\phi }{u}_i\right)=\frac{\partial }{\partial x_i}{\mathsf{\tau }}_i^{}+Q.$

The first term in (8) is an unsteady term, the second one is a convection term, and the latter one is a diffusion term. The equation can also contain sources and sinks represented by the latter term Q. The tensor τ_i contains spatial derivatives of passive scalar φ. For simplicity, assume that τ_i = μ$\frac{\partial \mathsf{\phi }}{\partial x_i}$. Although this assumption is computationally much more efficient, it has no impact on the generality of the methods.

Consider an arbitrary unstructured mesh as shown in Figure 1:

Here, k is a set of faces of cell P consisting of a set of inner faces, k_int, and a set of outer faces, k_s. The cell’s neighbor sharing inner face k_int is denoted by $M_{k_{int}}$. The area vector of face k is S_i,k, where i = 0 ÷ 2 is the vector component. The vector from the center of cell P to the center of cell M along face k_int is denoted by $d_{i,k_{int}}$ = r_i,M − r_i,P, and the vector from the center of P to the face center is denoted by $d_{k_s}=r_{k_s}-r_P$ where r_i is the radius vector.

Time discretization of Equation (8) by the second-order scheme is performed using the Adams–Bashforth scheme [23]:

(9)$\frac{3{\mathsf{\rho }}^{j+1}{\mathsf{\phi }}^{j+1}-4{\mathsf{\rho }}^j{\mathsf{\phi }}^j+{\mathsf{\rho }}^{j-1}{\mathsf{\phi }}^{j-1}}{2\Delta t}+{\left[\frac{\partial }{\partial x_j}\left(\mathsf{\rho }\mathsf{\phi }{u}_j\right)-\frac{\partial }{\partial x_j}{\mathsf{\tau }}_{ij}^{}-Q\right]}^{j+1}=0$

Here, $\Delta t$ is a time step. To discretize Equation (8) in space, let us integrate it over the volume V_P of cell P and pass on to the integration over the area S_P for the convection and diffusion terms:

(10)${\displaystyle \underset{V_P}{\int }\frac{3{\mathsf{\rho }}^{j+1}{\mathsf{\phi }}^{j+1}-4{\mathsf{\rho }}^j{\mathsf{\phi }}^j+{\mathsf{\rho }}^{j-1}{\mathsf{\phi }}^{j-1}}{2\Delta t}dV}+{\displaystyle \underset{S_P}{\oint }\mathsf{\rho }\mathsf{\phi }{u}_{i}d{S}_i}-{\displaystyle \underset{S_P}{\oint }\mathsf{\mu }\frac{\partial \mathsf{\phi }}{\partial x_i}d{S}_i}-{\displaystyle \underset{V_P}{\int }QdV}=0$

The source and the unsteady term are discretized as follows:

(11)${\displaystyle \underset{V_P}{\int }\frac{{\mathsf{\rho }}^{j+1}{\mathsf{\phi }}^{j+1}-{\mathsf{\rho }}^j{\mathsf{\phi }}^j}{\Delta t}dV=}{\left(\frac{{\mathsf{\rho }}^{j+1}{\mathsf{\phi }}^{j+1}-{\mathsf{\rho }}^j{\mathsf{\phi }}^j}{\Delta t}\right)}_P{V}_P,{\displaystyle \underset{V_P}{\int }QdV}=Q_P{V}_P.$

The discrete equivalent of the diffusion term is written in the following form [22]:

(12)${\displaystyle \underset{S_P}{\oint }\mathsf{\mu }\frac{\partial \mathsf{\phi }}{\partial x_i}d{S}_i}\approx {\displaystyle \sum _k{\left(\mathsf{\mu }\frac{\partial \mathsf{\phi }}{\partial x_i}\right)}_k{S}_{i,k}}={\displaystyle \sum _k{\mathsf{\mu }}_k{\left(\frac{\partial \mathsf{\phi }}{\partial n_k}\right)}_k\left|S_k\right|},$

where n_k is the normal of face k. The product under the summation sign on the right side contains a derivative along $\frac{\partial \mathsf{\phi }}{\partial n_k}$, which, for an orthogonal mesh, can be defined as

(13)$\frac{\partial \mathsf{\phi }}{\partial n_k}\left|S_k\right|=\frac{{\mathsf{\phi }}_M-{\mathsf{\phi }}_P}{\left|d_k\right|}\left|S_k\right|.$

For the finite-volume mesh approximation, the convection term is written as

(14)${\displaystyle \underset{S_P}{\oint }\mathsf{\rho }\mathsf{\phi }{u}_{i}d{S}_i}\approx {\displaystyle \sum _k{\mathsf{\rho }}_k{\mathsf{\phi }}_k{u}_{i,k}{S}_{i,k}}\approx {\displaystyle \sum _k{\mathsf{\rho }}_k{\mathsf{\phi }}_k{F}_k},$

where F_k is the volume flux across face k. The face value of the quantity, φ_k, is determined by the current discretization of the convection term. There are a lot of discretization schemes that can be used for arbitrary unstructured meshes [22,24,25,26]. Among them, one can distinguish some schemes having the highest usability rating in applied simulations: Upwind Differences (UD), Linear Upwind Differences (LUD), QUICK, Central Differences (CD), the Normalized Variable Diagram (NVD), and hybrid schemes (schemes listed above in combination with the upwind scheme to improve monotonicity). The schemes differ in their procedure of target quantity mapping onto faces and, as a result, have different dissipation properties.

To discretize the convection term in the volume fraction transport equation, we use special so-called “compressing” schemes [20,27,28]. Their distinguishing feature is that such discretization must keep the interface as thin as possible, that is, prevent interface smearing, and keep the volume fraction within a range between 0 and 1. One of the known “compressing” schemes is the HRIC scheme [27], which is used in this paper. The HRIC scheme is a “compressing” scheme used for preservation of the inter-phase boundary from dissipation, and it has been implemented previously within the numerical method under consideration [19].

Using this discretization, Equation (8) is replaced by a system of linear algebraic equations (SLAE) written for each cell:

(15)$A_P{\mathsf{\phi }}_P+{\displaystyle \sum _{k_{int}}{A}_{k_{int}}{\mathsf{\phi }}_{M_{k_{int}}}=R_{i,P}}.$

Here, A_P represents the diagonal elements of the matrix, A_kint represents the off-diagonal elements of the matrix, and R_i,P represents the elements of the vector for the right part. To solve this SLAE, we use a multigrid solver, AMG [29].

Below, we describe the implementation of the damping method.

2.2. Damping Method

The wave damping method employs a damping layer next to the open boundaries to absorb the kinetic energy of the incident wave. In real experiments, this effect is realized using a porous intermediate layer located at the edges of the experimental pool. To incorporate such absorption in momentum Equation (1), one can add one more source term, i.e., I_i:

(16)$\frac{\partial u_i}{\partial t}+\mathsf{\rho }\frac{\partial }{\partial x_j}\left(u_i{u}_j\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}{\mathsf{\tau }}_{ij}+\mathsf{\rho }{g}_i+I_i,$

The additional source of momentum I_i is a resisting force, counteractive to the flow pulse. This source is proportional and opposite in its sign to the velocity:

I_i = −δk_sεu_i. (17)

Here, k_s is the drag coefficient, the value of which determines the rate of kinetic energy absorption, and ε is a geometric factor, which is non-zero in the damping layer and which grows linearly from the front of the damping zone to its end:

(18)$\mathsf{\epsilon }=\min \left(\max \right(1-\frac{l}{L},1),0),$

where l is the shortest distance to the boundary and L is the width of the damping zone. Such linear distribution reduces the effect of wave reflection from the front of the damping zone.

According to Equation (18), the source, (17), acts only in the damping region, while its value is reduced to zero in the rest of the computational domain. The solution produced in the damping area is not physical due to artificially enforced suppression of the wave pulse. The basic task of this source is to reduce, at maximum, the amplitude of reflected waves.

There are two free parameters in this method: parameter k_s and the width of the damping layer L. These parameters directly control the efficiency of the wave absorption and, consequently, the amplitude of reflected waves. Their optimal values depend on the parameters of incident waves: wave amplitude, wavelength, and channel depth. To find the optimal values of k_s and L, we solve the problem of a solitary wave of height H traveling across a water tank. Suppose a wave with given parameters travels from the left boundary to the free boundary on the right.

Its initial waveform is defined by

η(x, 0) = Hsech²(γ(x − X_s)), (19)

where $\gamma =\sqrt{\frac{3H}{4d}}$, d is channel depth, L = X_s − X₀ = $\frac{\arccos h\left(\sqrt{20}\right)}{\gamma }$, X_s is a horizontal coordinate in point S, X₀ is a horizontal coordinate of the boundary, and h is the height of the wave ridge.

The initial wave velocity is

(20)$u\left(x,0\right)=\sqrt{\frac{g}{d}}\eta \left(x,0\right).$

To determine the optimal value of the damping parameter k_sour, we run simulations for channels of different depths. The cases under consideration are summarized in Table 1.

First, we run wave simulations with a varied channel depth to choose the optimal damping parameter k_sour, ensuring the minimum amplitude of the reflected wave. For the selected constant depth of the channel, we then vary the width of the damping zone L with the same wavelength.

Let us first consider the four cases specified in Table 1. To choose the optimal damping parameter, we compare the amplitude of the reflected wave with the amplitude of the incident wave at the tide gauge near the damping zone.

Figure 2 shows plots that compare the amplitude of the reflected wave with different damping parameters for the varied channel depths.

The plots suggest that the optimal value of k_sour for the channel depth of 0.32 m is 2000 s⁻¹; for 3.2 m, it is 750 s⁻¹; for 32 m, it is 175 s⁻¹; and for 320 m, it is 75 s⁻¹. Based on the resulting data, the following curve of the damping parameter as a function of the channel depth is plotted (Figure 3):

(21)$k_{sour}=\frac{1.1\cdot {10}^3}{\sqrt{d}},$

where d is the channel depth.

Next, we can define the obtained wave shape as a function of the form of changing damping coefficient ε. Three main options are considered, according to which the form is represented as the function from the exponent, the function from a sine, and the function from an antitangent. The damping coefficient was selected in compliance with Equation (21). A function distribution pattern is given in the following plot (Figure 4).

A series of numerical experiments from Table 1 were then carried out. Coefficient k_sour was chosen according to Equation (21). The character of the curve of the coefficient distribution in the area of the damping zone was a variable parameter. The Figure 5 shows the readings of tide meters for the differently scaled channels for four options of multiplier factor ε distribution. Symbols used in the following figure are introduced here: f(x) is a linear function, f(sin) is a function of sine, f(exp) is an exponential function, and f(arctg) is a function of an antitangent.

The results show that the least amplitude for the reflected wave is produced with a linear distribution of ε, and it is thus used in further computations.

The optimal width of the damping zone is determined by varying L at a constant k_sour. The wavelength λ is 40.24 m, and the wave height is 0.64 m. The channel parameters chosen were a length of 13λ and a depth of 3.2 m, and the width of the damping zone was varied (L = 1÷9λ). The damping parameter is defined by (21).

Figure 6 shows tide gauge data measured at a point located 20 m upstream of the front side of the damping zone. In the plots, the amplitudes of the wave incidents on the tide gauge are compared to the amplitudes of the reflected waves.

The plots demonstrate that the amplitude of the reflected wave predictably grows with an increase in the width of the damping zone. Amplitude relationships between reflected and incident waves for each case are given in Table 2.

The table indicates that the amplitude relationship between the reflected wave and the incident wave for a damping zone with a width of 6λ is less than five percent. A further increase in zone width has almost no effect on the amplitude of the reflected wave. As a result of these numerical experiments, we establish the criterion that controls the degree of wave reflection from the boundary.

Criterion. To ensure that the magnitude of the reflected wave does not exceed ten percent, the width of the damping zone must be in the range of two to six incident wavelengths. If the damping layer is wider than six incident wavelengths, the reflected wave does not exceed five percent.

The wave damping method is preferable when applied to real problems because it can be tuned to any specific problem by choosing the optimal damping parameters.

3. Numerical Experiments. Bathymetry-Aware Simulation of Tsunami Waves Traveling across a Real Water Body

To test the method for applying non-reflective boundary conditions, we ran numerical simulations of a decaying single pulse in a water area near Sakhalin Island while considering its real bathymetry. The problem geometry is shown in Figure 7. The size of the area in the X direction is 300,000 m, and in the Y direction it is 350,000 m. The maximum depth of the aquatic area reaches 7600 m. At the center of the water area, we applied a single pulse with a height of 30 m and a width of 40 km.

To simulate the propagation of the generated waves, a mesh model of 3 million cells was generated. The size of the cells in horizontal directions was ensured to be 2000 m. The mesh became more refined towards the surface of the phase division to simulate wave propagation more accurately; minimum cell size in the vertical direction for the region of phase division was 0.2 m. The required computational resources were 80–120 processors and the computation time was 3–5 h.

Its decay created cylindrical waves that traveled over the water surface. When the waves reached the open boundaries of the water area, they had to leave the domain without reflection. To test the performance of the method, we simulated two cases. In the first case, all the side boundaries were non-reflective. In the second case, the wall boundary conditions were applied.

Figure 7 shows the problem schematic. The non-reflective boundary is implemented as a near-boundary damping layer simulating a porous region, which absorbs incident waves (shown with a dashed line). Figure 7 shows the locations of the tide gauges. Tide gauge 1 was placed in the way of the propagating wave, and tide gauge 2 was placed in the damping zone. The tide gauge data from point 1 allowed us to estimate the amplitude of the waves reflected from the side boundary, while tide gauge 2 indicated the degree of attenuation in the damping zone.

Figure 8 shows the time series for the free surface position for both cases. Figure 9 presents the tide gauge data at points 1 and 2. As one can see in Figure 8, at 700 s, when the waves reached the side boundaries the wave patterns in both cases were the same, which is evidenced by the tide gauge data at point 1 (Figure 9) up to 1000 s. In the first case (non-reflective boundary conditions), at 1200 s and 1400 s, one can see that the waves decayed when they reached the side boundaries, while in the second case (wall boundary conditions) one can clearly distinguish that the waves reflected from the side boundaries. The tide gauge data at point 2 (Figure 9) indicate that the amplitude of the first incident wave in the damping zone was 65% lower. At 1200 s, the wave reflected from the side wall came back to tide gauge 1. The pane with its data (Figure 9) shows that the wave amplitude at this point in space and time, considering the reflected wave, was 80 percent higher in the second case than in the first case.

It is important to note that a part of the land located in the left upper corner of the domain is in the shallow water area, which results in the deceleration of the propagating wave and growth of its height. As research on groundswells is not the purpose of this numerical experiment, we ignored the fact that the wave that came in the direction of this piece of land had changed its propagation shape. In case we need to reproduce the groundswell of such waves at the shore without changes in its shape under the effect of damping force, we could reformulate the initial problem’s setup.

The amplitude of the second peak in the pane displaying tide gauge data at point 2 for the non-reflective boundary conditions is less than 10 cm (indicated by the arrow), which is considered to be an acceptable error due to the mesh resolution in the area of free surface. Even if it returns into the domain, this wave has almost no effect on the primary wave.

Thus, the numerical simulations demonstrate that the implemented method for non-reflective boundaries can be used to simulate tsunami waves in water areas within the world’s oceans while considering their real bathymetry.

4. Conclusions

In this paper, we considered a wave absorption method implementation in free-surface simulations. The discussed wave absorption method employs an artificial damping zone. We identified free parameters for the damping method and revealed an empirical relationship between the values of the damping zone parameters and the depth of the simulated channel. We also investigated different functions of distribution for the free geometric parameter of the method under consideration. As a case study, we simulated waves traveling in a water area in the Pacific Ocean near Sakhalin Island.

Author Contributions

Conceptualization, A.K. (Andrey Kozelkov) and A.K. (Andrey Kurkin); data curation, E.T., V.K. and D.U.; formal analysis, A.K. (Andrey Kozelkov) and V.K.; investigation, E.T., A.K. (Andrey Kozelkov), A.K. (Andrey Kurkin), D.S., V.K. and D.U.; methodology, A.K. (Andrey Kozelkov) and A.K. (Andrey Kurkin); software, E.T., V.K. and D.U.; supervision, A.K. (Andrey Kozelkov); validation, E.T., D.S. and D.U.; visualization, E.T. and D.U.; writing—original draft, A.K. (Andrey Kozelkov) and A.K. (Andrey Kurkin); writing—review and editing, A.K. (Andrey Kozelkov) and A.K. (Andrey Kurkin). All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Computational mesh.

Figure 2. Tide gauge data comparing the amplitude of the reflected wave with different damping parameters for the varied channel depths.

Figure 3. Damping parameter ksour as a function of the channel depth.

Figure 4. Distribution curves of the geometric factor ε.

Figure 5. Tide gauge data with different types of distribution of the geometric factor ε.

Figure 6. Incident and reflected wave amplitudes.

Figure 7. Problem schematic.

Figure 8. Time series of free surface position.

Figure 9. Tide gauge data.

Table 1

Problem settings.

Case	Channel Depth, m	Channel Length, mm	Wave Height, m	Wavelength, m	L, m
1	0.32	23	0.064	4.24	4
2	3.2	230	0.64	40.24	40
3	32	2300	6.4	402.4	400
4	320	23,000	64	4024	4000

Table 2

Relationship between the amplitude of the reflected wave and the width of the damping zone.

Width of Damping Zone	Amplitude Relationship between Reflected and Incident Waves, %
1λ	23.8
2λ	9.7
3λ	7.8
4λ	6.5
5λ	5.3
6λ	4.6
7λ	4.5
8λ	4.4
9λ	4.4

References

1. Fürst, J.; Musil, J. Development of non-reflective boundary conditions for free-surface flows. Proceedings of the Topical Problems of Fluid Mechanics; Prague, Czech Republic, 21–23 February 2018.

2. Mani, A. On the reflectivity of sponge zones in compressible flow simulations. Annual Research Brief; Ames Research Center: Mountain View, VA, USA, 2010.

3. Salvesen, H.-C.; Teigland, R. Non-reflecting boundary conditions applicable to general purpose CFD simulators. Int. J. Numer. Meth. Fluids; 1998; 28, pp. 523-540. [DOI: https://dx.doi.org/10.1002/(SICI)1097-0363(19980915)28:3<523::AID-FLD735>3.0.CO;2-5]

4. Giles, M. Non-Reflecting Boundary Conditions for the Euler Equations; Computational Fluid Dynamics Laboratory, Department of Aeronautics and Astronautics, Massachusetts Institute of Technology: Cambridge, MA, USA, 1988.

5. Granet, V.; Vermorel, O.; Leonard, T.; Gicquel, L.; Poinsot, T. Comparison of nonreflecting outlet boundary conditions for compressible solvers on unstructured grids. AIAA J.; 2010; 48, pp. 2348-2364. [DOI: https://dx.doi.org/10.2514/1.J050391]

6. Dorodnitsyn, L.V. Non-reflecting boundary conditions for the gas dynamics equations systems. J. Comput. Math. Math. Phys.; 2002; 42, pp. 522-549. (In Russian)

7. Givoli, D.; Neta, B. High-order nonreflecting boundary conditions for the dispersive shallow water equations. J. Comput. Appl. Math.; 2003; 158, pp. 49-60. [DOI: https://dx.doi.org/10.1016/S0377-0427(03)00462-X]

8. Carmigniani, R.A.; Violeau, D. Optimal sponge layer for water waves numerical models. Ocean Eng.; 2018; 163, pp. 169-182. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2018.05.068]

9. Wei, G. The Sponge Layer Method in FLOW-3D; Flow Science, Inc.: Santa Fe, NM, USA, 2015.

10. Peric, R.; Abdel-Maksoud, M. Reliable Damping of Free Surface Waves in Numerical Simulations. Ship Technol. Res.; 2016; 63, pp. 1-13. [DOI: https://dx.doi.org/10.1080/09377255.2015.1119921]

11. Hsu, T.W.; Ou, S.H.; Yang, B.D.; Tseng, I.F. On the damping coefficients of sponge layer in Boussinesq equations. Wave Motion; 2005; 41, pp. 45-57. [DOI: https://dx.doi.org/10.1016/j.wavemoti.2004.05.003]

12. Wang, D.; Dong, S. A discussion of numerical wave absorption using sponge layer methods. Ocean Eng.; 2022; 247, 110732. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.110732]

13. Torregrosa, A.J.; Fajardo, P.; Gil, A.; Navarro, R. Development of Non-Reflecting Boundary Condition for Application in 3D Computational Fluid Dynamics Codes. Eng. Appl. Comput. Fluid Mech.; 2012; 6, pp. 447-460. [DOI: https://dx.doi.org/10.1080/19942060.2012.11015434]

14. Kar, S.K.; Turco, R.P. Formulation of a Lateral Sponge Layer for Limited-Area Shallow-Water Models and an Extension for the Vertically Stratified Case. Mon. Wea. Rev.; 1994; 123, pp. 1542-1559. [DOI: https://dx.doi.org/10.1175/1520-0493(1995)123<1542:FOALSL>2.0.CO;2]

15. Kozelkov, A.S. The Numerical Technique for the Landslide Tsunami Simulations Based on Navier-Stokes Equations. J. Appl. Mech. Tech. Phys.; 2017; 58, pp. 1192-1210. [DOI: https://dx.doi.org/10.1134/S0021894417070057]

16. Kozelkov, A.S.; Efremov, V.R.; Kurkin, A.A.; Pelinovsky, E.N.; Tarasova, N.V.; Strelets, D.Y. Three dimensional numerical simulation of tsunami waves based on the Navier-Stokes equations. Sci. Tsunami Hazards; 2017; 36, pp. 183-196.

17. Kozelkov, A.S.; Lashkin, S.V.; Efremov, V.R.; Volkov, K.N.; Tsibereva, Y.A.; Tarasova, N.V. An implicit algorithm of solving Navier–Stokes equations to simulate flows in anisotropic porous media. Comput. Fluids; 2018; 160, pp. 164-174. [DOI: https://dx.doi.org/10.1016/j.compfluid.2017.10.029]

18. Kozelkov, A.S.; Kurkin, A.A.; Sharipova, I.L.; Kurulin, V.V.; Pelinovsky, E.N.; Tyatyushkina, E.S.; Meleshkina, D.P.; Lashkin, S.V.; Tarasova, N.V. Minimal basis tasks for validation of methods of calculation of flows with free surfaces. Trans. R.E. Alekseev NSTU; 2015; 2, pp. 49-69.

19. Tyatyushkina, E.S.; Kozelkov, A.S.; Kurkin, A.A.; Pelinovsky, E.N.; Kurulin, V.V.; Plygunova, K.S.; Utkin, D.A. Verification of the LOGOS Software Package for Tsunami Simulations. Geosciences; 2020; 10, 385. [DOI: https://dx.doi.org/10.3390/geosciences10100385]

20. Efremov, V.R.; Kozelkov, A.S.; Kornev, A.V.; Kurkin, A.A.; Kurulin, V.V.; Strelets, D.Y.; Tarasova, N.V. Method for taking into account gravity in free-surface flow simulation. Comput. Math. Math. Phys.; 2017; 57, pp. 1720-1733. [DOI: https://dx.doi.org/10.1134/S0965542517100086]

21. Chen, Z.J.; Przekwas, A.J. A coupled pressure-based computational method for incompressible/compressible flows. J. Comp. Phys.; 2010; 229, pp. 9150-9165. [DOI: https://dx.doi.org/10.1016/j.jcp.2010.08.029]

22. Hrabry, A.I.; Zaitsev, D.K.; Smirnov, Y.M. Numerical simulation of currents with free surface based on VOF method. Trans. Krylov State Res. Cent.; 2013; 78, pp. 53-64. (In Russian)

23. Roache, P. Computational Fluid Dynamics; Mir: Moscow, Russian, 1980; 618p. (Translated into Russian)

24. Jasak, H. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows. Ph.D. Thesis; Department of Mechanical Engineering, Imperial College of Science: London, UK, 1996.

25. Jasak, H.; Weller, H.G.; Gosman, A.D. High resolution NVD differencing scheme for arbitrarily unstructured meshes. Int. J. Numer. Methods Fluids; 1999; 31, pp. 431-449. [DOI: https://dx.doi.org/10.1002/(SICI)1097-0363(19990930)31:2<431::AID-FLD884>3.0.CO;2-T]

26. Gaskell, P.H. Curvature-compensated convective-transport—SMART, A new boundedness-preserving transport algorithm. Int. J. Numer. Methods Fluids; 1988; 8, pp. 617-641. [DOI: https://dx.doi.org/10.1002/fld.1650080602]

27. Muzaferija, S.; Peric, M.; Sames, P.; Schelin, T. A two-fluid Navier-Stokes solver to simulate water entry. Proceedings of the 22nd Symposium of Naval Hydrodynamics; Washington, DC, USA, 9–14 August 1998.

28. Waclawczyk, T.; Koronowicz, T. Remarks on prediction of wave drag using VOF method with interface capturing approach. Arch. Civ. Mech. Eng.; 2008; 8, pp. 5-14. [DOI: https://dx.doi.org/10.1016/S1644-9665(12)60262-3]

29. Volkov, K.N.; Kozelkov, A.S.; Lashkin, S.V.; Tarasova, N.V.; Yalozo, A.V. A Parallel Implementation of the Algebraic Multigrid Method for Solving Problems in Dynamics of Viscous Incompressible Fluid. Comput. Math. Math. Phys.; 2017; 57, pp. 2030-2046. [DOI: https://dx.doi.org/10.1134/S0965542517120119]