Absorbing Boundary Conditions for Simulating Water Waves Near Solid Bodies ()
1. Introduction
The propagation mechanisms of waves on a body of water have long been a topic of interest for researchers [1]-[7]. Surface water ripples are caused by a balance between gravity forces that maintain the horizontal free surface of the water, surface tension that preserves the consistency of the air-water interface, water inertia, and the difference between air and water pressures. The waves can be classified into gravity waves, capillary waves, or pressure waves based on the force that causes their motion. This work focuses on capillary-gravity surface waves that are generated by the movement of a solid body. These waves propagate around the body and interact with its rigid surface, resulting in a wake in its vicinity [8]. A common method for addressing scattering problems is to use a frequency domain approach for the entire unbounded structure, with analytical radiation boundary conditions applied at infinity [3]-[6] [8]. In addition to these analytical methods, numerical techniques such as the boundary element method [9] or the finite element method with infinite elements have been developed to address the same problem [10]. Otherwise, the domain is truncated by artificial boundaries to use a classical finite element-based approach with High Order Absorbing Boundary Conditions [11]-[14] or the Perfectly Matched Layer method [15] [16]. The main concerns are avoiding, or at least significantly reducing, the spurious boundary reflection waves and ensuring an accurate approximation of the infinite domain solution within the bounded domain. These methods are generally efficient. However, they require a more complex problem formulation and a large number of variables and computations [17].
Furthermore, when dealing with small wavelength waves, such as our case, a fine mesh is required to reduce dispersion errors and obtain numerical solutions with acceptable accuracy, and this requirement quickly makes the numerical solution process prohibitively expensive [18] [19]. Therefore, in the following, we aim to find an approximate low-order absorbing boundary condition to avoid increasing the already existing costs. After formulating the problem and specifying the underlying hypotheses, a linearisation around a steady state is performed, and new boundary conditions are introduced. A time-domain approach is considered without any prior assumptions about the solution form. The variational formulation of the problem is derived. A finite element approximation in space with a centred finite difference scheme in time is used to approximate the solution. The results obtained are presented and discussed in detail.
2. Statement of the Problem
The aim of this study is to investigate the dynamic behaviour of the water surface in the vicinity of a solid body that moves horizontally with a velocity
and may also have an oscillatory displacement. To achieve this, we study the irrotational and inviscid flow of a compressible and barotropic fluid around the structure fixed for a given period
. The structure is immersed in water, and its shape is simplified to a cylinder due to the presence of singularities at the contact points between the solid surface and the water surface, as well as at underwater angular points. The rectangular open domain Ω with a hole of radius
in its center is considered the computational domain. Its boundary
has a unit outward normal vector
.
corresponds to the bottom of the system,
, to the free surface of the water,
and
, to the sides through which the water flow enters and leaves Ω and
, to the rigid body surface (see Figure 1).
denotes the edges of
. A steady flow passes through Ω with horizontal velocity
and a small disturbance is introduced in the fluid as an initial condition.
![]()
Figure 1. Geometry and notations of the problem.
3. Theoretical Modelling
Several papers have addressed the issue of gravity capillary waves [20]-[23]. However, a relatively small number of studies have focused on the steady flow in the background [7]. For simplicity, it is generally assumed that the convective propagation medium is homogeneous and unbounded along the surface, and that the flow is irrotational and incompressible. Then associate velocity potential Φ satisfies the classical Laplace equation, (
), and the vertical displacement of the surface
verifies:
(1)
(2)
(3)
the subscript
stands for curvilinear abscissa along
and indicates that the differential operator is considered locally.
(4)
The dynamic boundary condition Equation (3) comes from the Young-Laplace hydrostatic balance equation introduced in Bernoulli’s equation for an inviscid and irrotational flow since the free surface of the water is a streamline. The surface tension is noted
, the density of the fluid
and the acceleration of gravity
. Under the small displacement theory, a plane harmonic waves analysis can be carried out, leading to the capillary-gravity wave dispersion relation:
(5)
where
is the angular frequency and
the wave number. The well-known capillary-gravity wave dispersion relation in still water is transferred to the reference frame moving with relative velocity
by the Doppler shift in Equation (5) [7]. Consequently, the relative mean flow rate affects the dispersion of waves according to their direction of propagation by influencing their wavelengths,
. The waves situated upstream of the current will exhibit a shorter wavelength than those located downstream. Consequently, these waves will predominantly manifest as capillary waves (
) to the left of the body case and as gravity waves (
) to its right, where
is the capillary constant [24].
Nevertheless, the relatively robust assumptions posited in the aforementioned formulation are designed for an analytical examination of the problem and may, therefore, be open to question. Firstly, the propagation medium is assumed to be homogeneous, even though there are different forces acting on the surface and within the fluid. This formulation assumes incompressibility of the flow, so the effects of pressure waves propagating in the internal fluid at the surface are not taken into account, as the flow velocities are much lower than the sound velocities in water. In the dynamic equilibrium of the free surface used to compute Equation (3), the inertia of the surface should be included for more physical accuracy. In our case, a time domain approach is more natural than a frequency domain approach, but the results are more difficult to interpret. In practice, the radiation boundary conditions are not satisfied at the boundaries of the domain. Therefore, the solutions obtained propagate to infinity without spatial or temporal decay [7]. Therefore, the propagating medium is then considered below as a stratified irrotational compressible fluid waveguide with a convective uniform mean flow, which is subject to a Gaussian pulse. Accordingly, the previous governing partial differential equations and the solution approach must be modified.
3.1. Hypotheses and Formulation of the Global Model
The propagating medium is considered to consist of two homogeneous layers of fluid with very different properties. The upper layer is free surface water, infinitesimally thick as skin and characterised by capillary-gravity wave propagation. The lower layer is the inner water with finite or semi-infinite thickness characterised by acoustic wave propagation. Therefore, two different models must be introduced to account for these features: an internal fluid model and a surface model.
3.1.1. Formulation of the Inner Fluid Model
We assume that the flow is characterized by two variables modelling the mass density
and the velocity potential Φ that satisfy:
(6)
(7)
where
is the barotropic potential,
, the fluid pressure and T, the simulation time.
3.1.2. Formulation of the Surface Model
Applying Newton’s second law of motion to an infinitesimal small surface element of thickness
that vertically moves of a displacement
, leads to the free surface equilibrium equation:
(8)
where
stands for material derivative. In comparison with the dynamic boundary condition Equation (3) that Equation (8) replaces, the inertia of the surface has been taken into account. The
value is chosen in our case so that the characteristic velocity of Equation (8), namely
, is equal to the velocity
which corresponds to the special case of a non-dispersive capillary-gravity wave with
in Equation (5) [20]. However, this value could be chosen arbitrarily to take account of the presence of particles on the surface.
3.1.3. On the Interface
The kinematic boundary condition Equation (2) becomes the continuity of normal velocity at the interface:
(9)
taking account of the rotation of the normal to the surface.
3.2. Linearization of the Governing Equations
The global nonlinear dynamical model obtained is linearized around a main steady state. Therefore, the global solution is split into a steady state obtained and a transient one. The lateral boundary conditions are defined separately according to the nature of the state. For the steady flow, the most realistic condition is to fix the normal velocity. For transient flow, non-reflecting boundary conditions have to be prescribed for the inlet and the outlet of Ω in order to avoid any spurious reflections of the waves reaching the boundaries of the domain.
3.2.1. Main Background Steady State Flow
We introduce
, the velocity potential corresponding to a steady quasi-uniform horizontal steady flow of an incompressible fluid that enters and leaves Ω at constant unit horizontal velocity with normal surface displacement variation along
regarded as negligible and non penetrability condition on
satisfied. This background flow is stationary with respect to the boat, which was chosen as the frame of reference.
is the solution to the following problem:
(10)
To obtain the related displacement on the surface
of this flow for our model, the stationary balance of forces on free surface Equation (8) along with homogeneous Neumann boundary conditions are applied. According to the frame of reference, it leads in stationary case to:
(11)
with
, the fluid density and
, the gravity acceleration. The solution of Equation (10) corresponds to a steady quasi-uniform horizontal flow with a digging effect due to the term
, as shown in Figure 2. The order of magnitude of free surface strain is about 10−3 m and, therefore, negligible in comparison with initial computational domain Ω dimensions.
Figure 2. Free surface vertical displacement
(m) solution of Equation (11) versus
(m) with the “digging” effect.
3.2.2. Transient State
We study the evolution of a small disturbance around the steady state
. The unsteady waves in the fluid are represented by the perturbation functions
of variables
and t. The problem is formulated with them wherein
,
,
. The solution is split into a main background steady-state component and a transient one. The domain Ω is then cropped by
,
,
,
and
, as shown in Figure 3. The new lateral boundaries
and
correspond to equipotential lines of
passing respectively through the left upper domain corner and right upper corner of Ω. As a result, the main steady flow crosses perpendicularly
and
and no reflected flow remains in the new domain Ω. The artificial boundaries are chosen far enough from the rigid body to consider that steady-state flow is uniform in this area, so the corners of the new domain are right-angled. The disturbance is so small that it is then reasonable to neglect the nonlinear terms in the governing equations Equation (6)-Equation (9) to obtain Equation (12)-Equation (15). Convective derivatives with flow velocity
are used to derive linearized equations Equation (13) and Equation (14). Hence,
are assured of satisfying the linearized enhanced Neumann-Kelvin’s model with capillarity:
![]()
Figure 3. Calculated geometry of the new computational domain in the case of a solid body size large compared to the initial computational domain. The digging effect on
and the inwards bowing of
and
are more pronounced than in our case.
(12)
(13)
given that
and
.
(14)
where
denotes the scalar product. Since the domain of study was reshaped,
and
are set to zero on
.
(15)
(16)
The initial condition corresponding to a disturbance taking place in the fluid: the functions
and
are set as Gaussian pulse functions in Ω and the functions
and
are fixed to zero on
.
3.2.3. Lateral Artificial Boundary Conditions
Non-reflecting boundary condition applied on inner fluid lateral edges is:
(17)
where
is a parameter related to the angle of incidence waves with respect to the normal of the boundary surface. This relation is consistent with Sommerfeld-like or zero-order non-reflecting boundary conditions for a wave that propagates at the phase velocity
corrected by the normal to the boundary component of velocity of the main background steady flow
. In the following, the parameter
is set to 1, this implies that the angles of incidence of impinging disturbances are close to the normal of the boundary. For such an order of approximation of absorbing boundary conditions, it is not beneficial to set
[14].
On the surface bounds
, a Sommerfeld-like non-reflecting boundary condition is difficult to apply since the surface is a very dispersive propagating medium, as can be stated from the dispersion relation Equation (5). Therefore, a value of velocity that approaches the apparent wave velocity value at boundaries can hardly be set. The more different this value is, the more spurious reflections you get. Furthermore, in the case of multi-layer models, numerical instabilities due to singularities appear at these points [25].
As the conditions are to be set for the points
that belong to the surface
and to the lateral boundaries
or
, the equations Equation (14), Equation (15) and Equation (17) are considered to devise the new boundary conditions. The guiding idea is to extend the non-reflective boundary condition applied on the lateral boundaries
or
to the surface
intersecting points
by means of the normal velocity continuity condition. Using the simplifying assumptions
,
and
led by choosing
and
away from rigid body, the following simplified equations must be satisfied on
:
(18)
(19)
(20)
(21)
For the right side (resp. the left side), the solution method consists in derivating first Equation (20) (resp. Equation (21)) with respect to
and Equation (19) with respect to
and
, in order to eliminate partial derivatives of
. Introducing the resulting expression into Equation (18) leads, after integrating with respect to time, to the following new boundary condition in cartesian coordinates for each edge,
(22)
with
depending on
. The symbol
denotes that condition is on the left boundary of
and
on the right one. The boundary conditions obtained are said non-local in time as they depend not only on the time
but also on the entire history of
on
.
4. Solution Method
A classical approach to address wave propagation in layered media is hardly applicable due to the complexity of the coupled equations involved. A weak form of the problem and then a finite element formulation are directly considered to obtain the solution.
Variational Formulation and Numerical Approach
Multiplying Equation (13) by
and Equation (14) by
, respectively, together with Green’s formula application leads to the following coupled variational formulation: find functions
such that
(23)
and
(24)
with
depending on
.
For physical field approximation, a finite dimension subspace
made of piecewise linear functions on a fixed mesh, characterized by element length
, is considered.
is spanned by
with
and
finite element shape functions on Ω and on
respectively.
is the coordinate vector of
relative to this basis. For the time domain approximation, a centered finite difference scheme for derivatives and a trapezoidal rule for the integral over time term are applied. The problem becomes finding
,
,
, such that
(25)
where
,
et
are sparse matrices;
a vector depending on
and
accounting for non-local condition term in time Equation (22) with non-zero component corresponding to the surface edges
.
are prescribed by initial conditions.
The computing process is fully automated. All the geometry operations and meshes are generated and updated automatically according to intermediate results by a batch program using Numpy and Scipy Python routines and GMSH. Due to the complexity of weak formulation terms, low-level generic assembly procedures of GETFEM++ is employed to make the assembly of the involved sparse matrices. To compute the solution of the large sparse system Equation (25), a parallel sparse direct solver (MUMPS) is used. For the post-processing handling, Matplotlib Python libraries, PARAVIEW and GMSH are utilized. A mesh convergence study is performed by reducing the characteristic size of elements,
, from
to
. As shown in Figure 4 and Figure 5, results converge upon the same solution as the mesh density increases. A satisfactory compromise between the accuracy of results and computing time can be achieved by choosing the value of
. This result is consistent with the order of magnitude of the wavelength of gravity-capillary waves of interest (
m). Indeed, the solution is curvy or even oscillates over the wavelength; this means using sufficient fine meshes or high-order piecewise polynomials to get a reliable approximation by the finite element method. A classic rule to reduce interpolation errors is to set the resolution of the wave
(here
) for a linear piecewise polynomial, but at a small wavelength, it appears to be insufficient due to numerical pollution identified in Helmholtz problems [26]. For numerical computations, values of parameters of Table 1 are used. The value of
is set to ensure that surface waves progress over deep water. The value of
is chosen to be less than
m·s−1 so that the flow becomes a uniform stream with constant velocity
at infinity [21]. The value of the half-thickness
of the water surface, calculated accordingly
, is equal to 7 × 10−4 m.
Table 1. Numerical values of parameters of the problem.
Parameters |
L (m) |
H (m) |
R (m) |
U (m·s−1) |
σ (N·m−1) |
Values |
1 |
1 |
5 × 10−2 |
0.15 |
0.075 |
Figure 4. Inner fluid velocity potential
(m2·s−1) at point of coordinates
versus time (10−5 s) for different element sizes of the mesh.
Figure 5. Inner fluid velocity potential
(m2·s−1) for different mesh densities along a part of the middle line of the computational domain (m) at
.
5. Main Results and Comments
Wave propagation phenomenon is monitored by the variation of
in the inner fluid and the variation of
on the surface respectively. The value of the time step chosen is
s in order to properly see wave propagation with a velocity of
across the extend of the computational domain Ω. As can be seen in Figure 6, reflecting waves appear on the bottom of the domain as on the surface of the immersed solid body, where homogeneous Neumann conditions are imposed to model non-penetrability of the fluid through them. In addition, no spurious reflecting wave appears to be present on either lateral side of the computational domain. Thanks to the hyperbolicity of the problem, in order to verify whether non-reflecting boundary conditions are satisfactory on inner fluid lateral edges Equation (17), a similar study is carried out according to the same previous calculation criteria on a larger computational domain in
direction, sized so as to avoid lateral side spurious reflecting waves during the all simulation time [14]. The new solution obtained is regarded as a reference solution. Both resulting waves are in phase, but a variable amplitude difference can be noticed. The wave is slightly reflected especially on its peak of amplitude for times when there is not many interference. Indeed, the chosen absorbing boundary conditions are not intended to handle such interfering situations (see Figure 7).
![]()
![]()
![]()
![]()
Figure 6. Propagation of velocity potential disturbance
at times:
,
,
,
(left to right, top to bottom). The order of magnitude of initial perturbation is 10−2.
Figure 7. Comparison of inner fluid velocity potentials
(10−3 m2·s−1) versus time (10−5 s) between extended and main computational domain on the middle of the right artificial lateral edge.
Figure 8. Propagation of disturbance
in Ω and related normal surface displacement
(10−5 m) on
versus
coordinate (m) at time
. The initial order of magnitude of
is 10−2 m2·s−1. Its propagating order of magnitude is 10−3 m2·s−1 and the order of magnitude of the normal displacement is 10−6 m.
Figure 9. Propagation of disturbance
in Ω and related normal surface displacement
(10−5 m) on
versus
(m) coordinate at time
. Its order of magnitude is 10−3 m2·s−1, and the corresponding normal displacement
order of magnitude is 10−6m.
Waves propagation in inner fluid results in deformation of the surface, as shown in Figure 8 and Figure 9. The corresponding normal displacement
propagates along the surface
. On each side of the surface,
, no spurious reflective wave is noticed. In the inner fluid layer, no wave related to any reflective surface on the surface is neither observed (see Figure 9). Then, the lateral boundary conditions introduced by Equation (22) also seem to be adequate for successfully modelling the propagating phenomenon on the surface. Nevertheless, the velocity of the phenomenon is the same as that of the inner fluid layer, which is not in complete agreement with surface layer material properties and wave propagation in stratified media theories. The expected value should be close to the order of the riddle velocity
. Therefore, no surface propagation phenomenon should be observed with time step
. That’s actually what happens when the initial disturbance is located just below or on the surface, as shown in Figure 10. The observed normal displacements
in Figure 8 and Figure 9 are not related directly to surface wave propagation, but rather primarily to the velocity potential acoustic wave propagation in the inner fluid layer and to the interface coupling between the potential
and the normal displacement
on
given by Equation (15). The energy transmitted to the surface layer by the inner layer remains stationary over the time range considered. Therefore, the application of lateral boundary conditions Equation (22) does not significantly affect the propagation phenomenon, and its accuracy can not be estimated with an initial perturbation in the inner fluid layer. Numerical simulations are then carried out in the case of an initial disturbance of the surface with time step
s. The functions
and
are set to zero in Ω and
and
are introduced as Gaussian pulse functions on
.
![]()
Figure 10. Propagation of normal surface displacement
on
versus
coordinate and related disturbance
in Ω at time
. Initial perturbation is located on the surface of the fluid. The order of magnitude of
is 10−6 m. The order of magnitude of the potential
transmitted to the surface of the inner fluid is 10−6 m2·s−1.
Waves propagate and get out of the computational domain without generating significant spurious reflections (see Figure 11). Similarly to the previous case to verify non-reflecting boundary conditions on surface edges Equation (22), a larger computational domain in
direction is chosen to compute a reference solution. Both resulting waves are in phase, but a varying amplitude difference can be noticed due to existing spurious reflections (see Figure 12), which fade away over time (see Figure 11). According to the ratio between the orders of magnitude of the inner fluid potential and the surface displacement noticed in each calculated case, it comes out that the inner fluid wave propagation effect is not significant in the case of initial disturbance near the surface or on the surface itself. Indeed, a velocity potential of the order of magnitude of 10−3 m2·s−1 on the surface leads to a normal displacement response of the order of magnitude of 10−6 m in the inner fluid initial perturbation case. But in the surface initial disturbance case where the order of magnitude of normal displacement is 10−6 m, the velocity potential barely reaches 10−6 m2·s−1 on the surface, and the linearity of the model leads to a normal displacement response of 10−9 m, therefore negligible compared to 10−6 m (see Figures 8-10, Figure 13). Thus, the waves propagate mainly in the surface layer guided in the medium of smaller velocity in totally agreement with wave propagation theories
![]()
![]()
![]()
![]()
Figure 11. Displacement
(10−5 m) of the surface
versus
coordinate (m) at:
,
,
and
. The initial disturbance is located on the surface
. The waves propagate without any instabilities on surface edges, but spurious reflections on the surface are still present.
Figure 12. Comparison between normal displacements
(10−6 m) of surface
in the initial domain (red) and extended domain (blue) cases versus
(m) at
.
Figure 13. Propagation of normal surface displacement
on
versus
coordinate and related disturbance
in Ω at time
. The initial disturbance is located on the surface of the fluid. The order of magnitude of propagating
is 10−6 m. The order of magnitude of the potential
transmitted to the surface of the inner fluid is 10−6 m2·s−1.
in stratified media. Actually, during surface wave propagation, a small amount of energy is steadily transferred from the surface to the inner fluid, which is immediately removed from the computational domain, as in an incompressible fluid. The velocity potential
rendering (see Figure 13) comes from the superposition of all velocity potential waves generated by the propagating surface wave at all times [27] [28]. Therefore, these results can hardly be analyzed and used to clearly draw any possible conclusions on the compliance of the non-reflecting boundary condition chosen Equation (17) in the inner fluid layer.
6. Conclusion and Suggestion
For the modelling of the wake of a solid body moving through a body of water, a wave propagation problem in a convective stratified media has been considered. To numerically solve the problem, the finite elements method is applied due to its versatility and the accuracy of its results in dealing with complex configurations. Nevertheless, since it is an acoustic scattering problem with large wavenumber waves, which is of interest, a fine mesh in the artificially bounded domain has to be used, and appropriate non-reflecting boundary conditions are to be sought while keeping computational costs low. In addition, the significant differences between layer properties make it difficult to address the entire problem with traditional schemes, a non-local in-time boundary condition has been devised by taking into account all the conditions that must be met on the artificial lateral edges of the computational domain. This work has highlighted some of the complex phenomena that involve the coupled propagation of surface and volume waves at different time scales and with very different orders of magnitude, features that cannot be observed under the assumption of fluid incompressibility. To go further with the same model, it is necessary to reduce the size of the linear system by using suitable enriched basis functions in order to decrease the number of elements per wavelength and to be able to increase the order of the absorbing boundary condition to eliminate spurious reflections [29] [30]. Due to the differences in scale between the phenomena occurring in each layer and the weak feedback from the inner fluid to the surface, a more simple one-dimensional model could also be considered by adding a damping term to model the energy dissipation of the surface propagating wave into the fluid.
Acknowledgements
I would like to express my gratitude to Professor Philippe Destuynder, Emeritus Professor at the University of La Rochelle (Laboratory LaSIE), for the valuable insights and discussions we have had.