A Novel FDTD–PML Scheme for Noise Propagation Generated by Biomimetic Flapping Thrusters in the Ocean Environment

: Biomimetic ﬂapping-foil thrusters can operate efﬁciently while offering desirable levels of thrust required for the propulsion of a small vessel or an Autonomous Underwater Vehicle (AUV). These systems have been studied both as main propulsion devices and for augmenting ship propulsion in waves. In this work, the unsteady hydrofoil loads are used to calculate the source terms of the Ffowcs Williams–Hawkings (FW-H) equation which is applied to model noise propagation in the underwater ocean acoustic environment. The solution provided by a simpliﬁed version of the Farassat formulation in free space is extended to account for a bounded domain and an inhomogeneous medium, characterizing the sea acoustic waveguide. Assuming the simplicity azimuthal symmetry of the environmental parameters, a numerical model is developed based on a Finite Difference Time Domain (FDTD) scheme, incorporating free-surface and seabed effects, in the presence of a variable sound speed proﬁle. For the treatment of the outgoing radiating ﬁeld, a Perfectly Matched Layer (PML) technique is implemented. Numerical results are presented illustrating the applicability of the method. problem, the hydroacoustic model is extended to treat the bounded waveguide featuring the inhomogeneity of parameters. To address the latter complexity, an axisymmetric cylindrical environment is used to account for the geometrical spreading. The developed numerical method is based on a ﬁnite difference scheme in the time domain (FDTD) able to incorporate the effects imposed by the free surface, seabed boundaries, as well as a variable sound speed proﬁle. For the open radiation boundaries, a PML reformulation of the acoustic wave equation for the truncated inhomogeneous ocean waveguide is presented, and subsequently employed in the numerical scheme.


Introduction
Ship traffic is responsible for the significant rise in ambient noise at low frequencies (10-100 Hz) in many ocean regions, whose levels have increased by 15-20 dB over the last century [1,2]. Underwater radiated noise from anthropogenic sources has a negative impact on the ocean acoustic environment, disrupting navigation and communication of marine mammals. In order to control the ambient ocean noise levels due to human activities the International Maritime Organization [3] have issued non-mandatory guidelines to minimize the introduction of incidental noise from commercial shipping operations into the marine environment to reduce potential adverse impacts on marine life. Ships generate noise primarily by (a) propeller action, (b) propulsion machinery, and (c) hydraulic flow over the hull. Significant contributions of propeller noise are associated with cavitation, especially volume fluctuation of the sheet and tip vortex cavities and bubbles collapsing in the slipstream, as well as with the periodic blade passage and other sources such as propeller singing if present, producing an overall broadband noise spectrum that ranges from a few Hz to over 100 kHz [3,4]. In the lower frequency band, noise excitation from marine propellers, especially under partial cavitating conditions, matches the first octave bands, having a negative impact on marine mammals. This might interfere with animal communication and hinder acoustic signal detection. Under such framework, the hydroacoustic performance of ships and marine vehicles has become a critical factor to be evaluated and controlled to maintain a certain level of radiated noise. In recent years, flapping-foil thrusters have been investigated both as main and auxiliary propulsion devices for ship propulsion in waves [5][6][7] contributing to the reduction of power feeding to propellers and mode of combined/complex oscillatory motion of a biomimetic propulsion system with significant contribution to propulsion energy savings; see [6,19,20].
The flapping thruster operates as an unsteady hydrofoil combining heaving and pitching motion with a phase difference about 90 degrees, as depicted in Figure 1. Th most important non-dimensional parameters are the Strouhal number is the frequency (with ω denoting the angular frequency), 0 h the heav ing amplitude and U the forward travelling speed, the heave to chord ratio is 0 hc, th pitching amplitude is 0  , and the phase difference ψ between the foil heave and pitc oscillatory motions. The case of the flapping wing of a large aspect ratio is considered in this study. Fo simplicity, using the unsteady hydrofoil theory (see, e.g., [21]) the lift force L F can b estimated as follows: The case of the flapping wing of a large aspect ratio is considered in this study. For simplicity, using the unsteady hydrofoil theory (see, e.g., [21]) the lift force F L can be estimated as follows: where p/c is pitching axis location, k = ωc/(2U) is the reduced frequency, C(k) is the Theodorsen function (lift deficiency factor), A is the area of the foil, and H = AR(AR + 2) −1 is a 3D correction from lifting-line theory (elliptic wing). The instantaneous angle of attack α(t) is given by where the foil pitching and heaving oscillatory motions are defined as follows: θ(t) = θ 0 sin(ωt + 0.5π) and h(t) = h 0 sin(ω t) The horizontal and vertical foil forces can be calculated as follows: F z = F L cos(θ) + F D sin(θ) (5) where F D = 1 2 ρU 2 A C D (Re, α) denotes a drag force that can be approximately considered by using a drag coefficient based on the Reynolds number and corrected for the angle of attack; see also [22].
The above simplified method provides us with a first estimation of the hydrodynamic loads that can be exploited in the sequel for estimating the noise generation and propagation by the oscillating foil. More accurate predictions for the hydrodynamic pressure field and loads associated with the operation of the flapping thruster can be derived by the implementation of unsteady BEM and CFD, as, e.g., presented by [23] taking also into account the effects of the free surface and incident waves. Moreover, systematic results and technical diagrams are provided by BEM [24].
As an example, a flapping wing of AR = 6 with symmetrical NACA0012 sections is considered and results concerning the time history of the dynamic angle of attack and the flapping thruster loads against non-dimensional time t/T, are presented in Figure 2. In this case, the flapping thruster operates at a Strouhal number Str = 2 f h 0 /U = 0.23, h 0 /c = 0.75, θ 0 = 23 deg, where f = ω/2π is the frequency and T = 1/f the period, respectively, U is the foil steady travelling speed, h 0 the heaving amplitude, c is the chord in the centerplane of the wing and θ 0 is the self-pitching amplitude about a pivot axis located at a distance of c/3 from the wing's leading edge; for more details see Ref. [25].
denotes a drag force that can be approximately considered by using a drag coefficient based on the Reynolds number and corrected for the angle of attack; see also [22]. The above simplified method provides us with a first estimation of the hydrodynamic loads that can be exploited in the sequel for estimating the noise generation and propagation by the oscillating foil. More accurate predictions for the hydrodynamic pressure field and loads associated with the operation of the flapping thruster can be derived by the implementation of unsteady BEM and CFD, as, e.g., presented by [23] taking also into account the effects of the free surface and incident waves. Moreover, systematic results and technical diagrams are provided by BEM [24].
As an example, a flapping wing of 6 AR = with symmetrical NACA0012 sections is considered and results concerning the time history of the dynamic angle of attack and the flapping thruster loads against non-dimensional time t/T, are presented in Figure 2. In this case, the flapping thruster operates at a Strouhal number is the frequency and T = 1/f the period, respectively, U is the foil steady travelling speed, 0 h the heaving amplitude, c is the chord in the centerplane of the wing and 0  is the self-pitching amplitude about a pivot axis located at a distance of /3 c from the wing's leading edge; for more details see Ref. [25]. In Figure 2a, the dynamic angle of attack is plotted and in Figure 2b the vertical (lift) force coefficient is shown by using a black line and the horizontal (thrust) force by using a cyan line, as calculated using the unsteady BEM method [26] and they are compared against the measured data [25] shown using symbols and the unsteady hydrofoil theory results using red lines. It is observed that the present method provides compatible predictions concerning the integrated quantities with unsteady hydrofoil theory and the experiment satisfactorily approximating the maxima of both the lift and thrust forces.

Noise Propagation in Free Space
The first integral approach for acoustic propagation is the Lighthill Acoustic Analogy [27]. The Ffowcs Williams and Hawkings (FW-H) [28] equation represents an extension of Lighthill's Analogy, derived by rearranging the Navier-Stokes equations into a wave equation [29]. In recent years, this numerical approach is being used increasingly for applications to marine and maritime problems [30][31][32][33][34]. This method propagates the acoustic shown by black and red lines and experimental data shown by using circles, and C X = F X / 0.5ρU 2 A shown using cyan lines and experimental data shown using crosses.
In Figure 2a, the dynamic angle of attack is plotted and in Figure 2b the vertical (lift) force coefficient is shown by using a black line and the horizontal (thrust) force by using a cyan line, as calculated using the unsteady BEM method [26] and they are compared against the measured data [25] shown using symbols and the unsteady hydrofoil theory results using red lines. It is observed that the present method provides compatible predictions concerning the integrated quantities with unsteady hydrofoil theory and the experiment satisfactorily approximating the maxima of both the lift and thrust forces.

Noise Propagation in Free Space
The first integral approach for acoustic propagation is the Lighthill Acoustic Analogy [27]. The Ffowcs Williams and Hawkings (FW-H) [28] equation represents an extension of Lighthill's Analogy, derived by rearranging the Navier-Stokes equations into a wave equation [29]. In recent years, this numerical approach is being used increasingly for applications to marine and maritime problems [30][31][32][33][34]. This method propagates the acoustic pressure by means of the analytical solution of the D' Alembert's equation in free space. Therefore, it cannot incorporate inhomogeneities as well as boundaries and interactions with surfaces. The pressure field and the forces produced by the flapping thruster motion, obtained using the BEM model described in the previous section, are used to calculate the dipole and monopole forcing terms associated with the initial boundary value problem governed by the FW-H equation. The acoustic model in free space is derived based on the Farassat formulation [35].
As the flapping thruster operates, it is subjected mainly to unsteady pressure loads. Low-frequency noise is caused by the fluctuations of foil pressure and volume flow distur-bance due to oscillatory motion. The formulation for the acoustic pressure generated from rotating machinery based on the FW-H equation is described as follows: where c is the speed of sound in the medium (1500 m/s ≤ c ≤ 1550 m/s for seawater), while the various terms on the right-hand side correspond to the acoustic monopole, dipole, and quadrupole source terms [36], respectively, and are defined as follows: where u n is the normal velocity on the foil, l i represents the corresponding loads, and T ij denotes the stresses.
The quadrupole term becomes more important for high frequency and strongly transonic flow phenomena. Simplified models based on the above assumption have been presented by various authors, e.g., by [37] in the case of cavitating marine propellers. Considering that the speed of sound in water is much greater than the flow velocities and focusing on the low-frequency part of the generated noise spectrum, the contributions by the latter term are neglected in the present work. The Farassat formulation is employed offering an integral representation of the solution of Equation (6) forced by the monopole and dipole terms. This is a rearrangement of the basic conservation laws of mass and momentum into an inhomogeneous wave equation, where the different noise generation mechanisms are identified and expressed by separate source terms. Using the free-space Green's function for the wave equation allows to write the two linear terms in an integral form. In the infinite homogeneous domain, the solution is given by thickness and loading components as follows: The loading term is given by where f = 0 indicates the moving surfaces, dp denotes the pressure difference between the lower and upper sides of the foil surface, M r denotes the Mach number in the r− direction, and the integrand is calculated at the retarded time. For relatively large distances of the observation point from the wing, we use the approximation where x T denotes the center of lift and thrust force on the flapping wing. Using the fact that the Mach number is very small, the following simplification is derived from Equation (9): where F(t r ) denotes the fluctuating unsteady part of the foil force, mainly composed from vertical (lift) and horizontal (thrust) forces, and t r = r/c denotes the retarded time between the observation point x 0 and the disturbance generating point x T . Similarly, for the thickness effect we have: which is approximated by where x Q (t r ) is the center of volume Q c displaced by the foil. In the case of an unsteady cavitating foil thruster, the latter term will also include the volume fluctuation of the sheet and tip vortex cavities and bubbles which in several cases, e.g., flapping foils and marine propellers, has a periodic or nearly periodic character. Indicative results concerning the calculated acoustic field in free space generated by the flapping thruster operating in water (c = 1500 m/s), as obtained using the above equations, are presented in Figure 3. In particular, the case of the foil of Figure 1 is considered, flapping at Str = 0.23, h 0 /c = 0.75 and θ 0 = 23 deg, and the acoustic field is calculated using Equations (11) and (13) using the hydrodynamic loads of Figure 2b. The result in the near field is plotted at two instants, when the nose of the flapping thruster is directed upwards and downwards, respectively, and the field is normalized with respect to its value at the distance of 1 m from the source point. (1 ) is the center of volume c Q displaced by the foil. In the case of an uns cavitating foil thruster, the latter term will also include the volume fluctuation of the and tip vortex cavities and bubbles which in several cases, e.g., flapping foils and m propellers, has a periodic or nearly periodic character.
Indicative results concerning the calculated acoustic field in free space generat the flapping thruster operating in water ( 1500 c = m/s), as obtained using the above tions, are presented in Figure 3. In particular, the case of the foil of Figure 1

 =
, and the acoustic field culated using Equations (11) and (13) using the hydrodynamic loads of Figure 2b result in the near field is plotted at two instants, when the nose of the flapping thrus directed upwards and downwards, respectively, and the field is normalized with re to its value at the distance of 1 m from the source point.  In the case of periodic signals, the acoustic force F(t) and cavity volume data Q(t) can be represented by the Fourier series and using the results in conjunction with the Green's function of the Helmholtz equation in free space: and its derivative: we find that the expressions provided by Equations (11) and (13) are the time-domain equivalent of the solution of the Helmholtz equation with complex monopole and dipole source intensities provided by: that is The above equations provide us with the expression of the complex monopole and dipole source terms at the basic frequency and its multiples, as provided by the hydrodynamic flapping foil responses described above. The above are sufficient for the study of acoustic noise generation and propagation from the flapping thruster in free space (unbounded domain). However, in the case of a bounded domain as it happens in underwater acoustics where the propagation is affected by the presence of the free surface and seabed boundary, the above results can be used to determine the source intensity. Subsequently, the propagation of noise emitted by the flapping foil in the ocean acoustic waveguide is obtained by formulating and numerically solving the forced acoustic equation, considering refraction, reflection, and scattering effects, as described in the next section.

Noise Propagation in the 3D Ocean Acoustic Waveguide
We consider the ocean acoustic waveguide limited above by the sea surface and below by the seabed; see Figure 4. Both Cartesian x = (x, y, z) and cylindrical coordinates , z are used in the formulation with the origin at a point located above the acoustic source on the free surface. For the considered problem, the hydroacoustic model is extended to treat the bounded waveguide featuring the inhomogeneity of parameters. To address the latter complexity, an axisymmetric cylindrical environment is used to account for the geometrical spreading. The developed numerical method is based on a finite difference scheme in the time domain (FDTD) able to incorporate the effects imposed by the free surface, seabed boundaries, as well as a variable sound speed profile. For the open radiation boundaries, a PML reformulation of the acoustic wave equation for the truncated inhomogeneous ocean waveguide is presented, and subsequently employed in the numerical scheme.
The above equations provide us with the expression of the complex monop dipole source terms at the basic frequency and its multiples, as provided by the hy namic flapping foil responses described above. The above are sufficient for the s acoustic noise generation and propagation from the flapping thruster in free spa bounded domain). However, in the case of a bounded domain as it happens in und acoustics where the propagation is affected by the presence of the free surface and boundary, the above results can be used to determine the source intensity. Subseq the propagation of noise emitted by the flapping foil in the ocean acoustic waveg obtained by formulating and numerically solving the forced acoustic equation, co ing refraction, reflection, and scattering effects, as described in the next section.

Noise Propagation in the 3D Ocean Acoustic Waveguide
We consider the ocean acoustic waveguide limited above by the sea surface low by the seabed; see

The 3D Acoustic Problem in a Constant Depth
The forcing terms associated with the monopole and vertical dipole terms are clearly azimuthally symmetric and in an environment where the parameters do not depend on θ the generated fields are also azimuthally symmetric. In the case of excitation by the horizontal dipole source, the forcing term is explicitly dependent on θ, and we approximately assume that the same holds also for the generated field in the above environments. Using the above facts, assuming the azimuthal symmetry the hyperbolic acoustic equation expressed in cylindrical coordinates x = (r, θ, z) takes the following form: The homogeneous Dirichlet boundary condition is used for modeling the free surface Equation (20), and for simplicity the seabed is modeled as a hard boundary using the Neumann boundary condition, Equation (21). In Equation (19), F(x; t) is the forcing term that includes both dipole and monopole sources, and c(x) denotes the variable sound speed. Since the hydroacoustic problem is inherently defined in the infinite ocean waveguide, the numerical methods need to incorporate appropriate termination conditions that ensure field radiation at the numerically truncated boundary.
The solution of problem Equations (19)-(21) will be based on the time integration technique using a finite-difference numerical scheme. Expanding the idea regarding the determination of the monopoles and dipoles in the previous section, we need to formulate suitable expressions for the acoustic sources into the time domain produced by the Helmholtz equation, as well. Thus, the F(x; t) will be determined in the source near field by the Green's function of the Helmholtz equation in the bounded waveguide and its derivatives, which will be presented in Section 4.2. Additionally, the outward radiation condition will be satisfied by a Perfectly Matched Layer (PML) technique and will be discussed in detail in Section 4.3. Finally, the FDTD-PML numerical scheme and the numerical results and the corresponding verifications with analytical solutions will be presented Section 4.4.

Determination of the Source Near-Field Data in the Acoustic Waveguide
The singular nature of the cylindrical Laplacian operator, along with the presence of point singularities (monopole and dipole sources) pose difficulties in numerical implementation. In the present work, these challenges are overcome by the enforcement of the field data in the near-field region. In this procedure, the original domain is divided into two subdomains, the computational domain Ω and the residual domain D (or relaxation zone) as shown in Figure 5. The problem is solved analytically in D = {r < r REL , −h < z < 0, 0 ≤ θ < 2π} for monopole and dipole force terms, and the data concerning the acoustic field are used to propagate the solution numerically in  For this purpose, the Green's function G(x|x s ) of the Helmholtz equation, satisfying the free-surface and hard seabed boundary conditions, is implemented in the near-source subregion of the waveguide: where δ(x − x S ) is the Dirac delta function. The corresponding expressions modeling the dipole terms associated with the lift and thrust foil forces are defined as solutions of: and are obtained as: where G is the Green's function and x s denotes the position of the singularity; see [38]. Thus, the present monopole and dipole (horizontal and vertical) sources data in the near-field region are obtained from the corresponding Green's function of the Helmholtz equation by using the harmonic time dependence, i.e., G(x|x s ; t ) = Re G(x|x s ; ω )e −iωt . Next, we consider the source region to be range independent and thus the monopole source denoted by G M (x|x s ; ω ) is obtained using separation of the variables as follows [12]: where the normalized vertical eigenfunctions are: where δ nm is the Kronecker delta and solutions to the corresponding Vertical Eigenvalue Problem (VEP) are obtained as: with boundary conditions: where k rm represents the eigenvalues. In the case of a homogeneous environment, the VEP is solved analytically, whereas in the case of a vertical sound speed profile the Sturm-Liouville problem is solved numerically by, e.g., the Finite Difference Method (see [12]).
In the case of a sound field generated by a horizontal dipole (associated with the foil thrust forcing) and a vertical dipole (associated with the lift force) the corresponding expressions are as follows: where ∂ x s r = (x − x s )/r = cos θ and θ is the azimuthal angle. The above analytical expressions obtained using the normal-mode expression of the Green's function of the Helmholtz equation and its derivatives in D are converted to the time domain as follows F(x, t) = Re ∑ n (Q n G(x|x s ; ω n ) + X n G H (x|x s ; ω ) + Z n G H (x|x s ; ω )) e −iω n t , r < r REL , (32) in order to be used to enforce near-source field data in Equation (19) by applying a relaxation scheme. Subsequently, the problem is numerically solved in cylindrical coordinates by using a Finite Difference Time Domain solver, in conjunction with a Perfectly Matched Layer technique for the satisfaction of the radiation condition (FDTD-PML), which is described in detail in the next subsections.

Incorporation of the Perfectly Matched Layer into the Acoustic Equation
To address the need of modeling general boundaries and inhomogeneities in a large domain, an absorbing boundary technique dissipating outgoing waves with minimum backscattering and without excessive computational cost is considered; see, e.g., [39]. In the present work, a Perfectly Matched Layer (PML) method is used based on Berenger's [39] original PML formulation. In general, PML formulations in the time domain reorganize the second-order wave equation into a system of first-order equations, by introducing pseudovariables [40][41][42][43][44][45]. Alternatively, the direct treatment of the second-order wave problem has been proposed by [46,47]. Finally, optimal PML formulations for the Helmholtz equation are proposed in [48,49].
In this work, the derivation of the PML method for a 3D cylindrical waveguide is based on Berenger's model. The basic step consists of applying a complex change in variables as follows (coordinate stretching): However, in our case only r-stretching is required where L denotes the PML thickness. In the present case, the variable absorption coefficient within the layer is denoted as σ r and is defined as a function of r, as it will be provided in the sequel. Equation (19) can be expressed as a system of equations by introducing the pseudovariable ∂ t w = ∇p and the following system is derived, where B(x) = σ r 0 0 σ z , A(x) = w r 0 0 w z and w = (w r , w z ).
The damping terms σ r , σ z are assumed to be non-negative functions. Within the absorbing layer, where σ r , σ z > 0, the wavelike solutions are transformed into exponentially decaying ones. In the computational domain of interest, the above terms vanish (σ x = σ z = 0) and the system reduces to the 3D cylindrical acoustic equation. The numerical solution is obtained using a finite-difference discretization and time marching integration. Noticing the excitation foil forces in Figure 2 and the different main frequency and harmonics involved in the time series of lift and thrust, the time-domain scheme is considered more efficient for the propagation of the source data.

The FDTD-PML Model
The aim of this section is to introduce the Finite Difference Time Domain-Perfectly Matched Layer (FDTD-PML) solver used as a propagator, along with the relaxation scheme described above, for discretizing the system of Equation (35) and enforcing the source data, respectively. All spatial derivatives are discretized with centered finite differences, which guarantees second-order approximation in space. For time discretization, centered finite differences for the first-and second-order time derivatives on a uniform mesh is also employed.
Let ∆t > 0 denote the time step size and ∆r, ∆z > 0 denote the spatial mesh sizes in the in the r, z directions, respectively. Denoting the nth time step as t n = n∆t and the spatial nodes r i = i∆r, z i = j∆z, the Laplacian is discretized as follows: The first-and second-order time derivatives are provided, respectively, as, Next, by substituting Equations (36) and (37) into the modified wave equation, Equation (35), we obtain the following discrete scheme: The above scheme is used explicitly together with the system of equations which includes the pseudovariable w = (w r , w z ), where ∇ i,j = (∂ r , ∂ z ) i,j denotes the finite difference approximation of the gradient on the i,j grid point of the spatial mesh.
In the sequel, we restrict our attention to problems formulated at the given azimuthal angle θ and the PML subregions correspond to large r > r PML . Although the PML has very good properties in its continuous form, a number of factors including the spatial discretization, the PML thickness L, the mesh density, and the absorption function need to be optimized in the discrete version. These aspects have been studied for classical wave problems in various works [50][51][52]. Within the context of the present work, the absorbing function takes Berenger's [40] polynomial form, where σ 0 is a positive parameter, L is the thickness of the layer, and m is the polynomial degree. Suitable choices for σ 0 , L, and m are required to ensure minimal numerical solution contamination within the region of interest through interface reflection as well as back propagation from the fictitious termination boundary. In the present work, L is selected to be of the order of characteristic acoustic wavelength and m = 3.
As already discussed in the previous subsection, the source, the free-surface, and the hard seabed data are enforced by a relaxation scheme at each time step providing the field data in the notes p n i,j = G M,H,V r for (i,j) in the source field region, where h denotes the waveguide depth.

Numerical Results and Discussion
A series of numerical results, obtained by means of the developed FDTD-PML scheme, including the free surface and seabed effects, are presented in this section. Acoustic propagation is considered in both homogeneous environments, featuring constant velocity c = 1500m/s, as well as in more realistic cases where the velocity exhibits vertical variation. The employed variable velocity profile is shown in Figure 6 and described as c(z) = 5 · 10 −5 z 2 + 0.05z + 1.5125 · 10 3

Numerical Results and Discussion
A series of numerical results, obtained by means of t scheme, including the free surface and seabed effects, are prese tic propagation is considered in both homogeneous environm locity 1500 c m s = , as well as in more realistic cases where th variation. The employed variable velocity profile is shown in F ( )  First, a homogenous (isovelocity) environment is considered, with the constant depth set as h = 3000 m and the calculated acoustic field excited by monopole, horizontal, and vertical dipole sources of constant intensity is presented in Figures 7-10. In order to illustrate the time evolution of the solution, in Figure 7 the acoustic field generated by the monopole source is plotted at time instants after 1, 3, 7, and 14 periods from the beginning of the calculation. Subsequently, in Figures 8-10 the periodic results obtained after 15 periods are compared against the solutions of the Helmholtz equation, for the monopole and dipole sources, respectively. In order to avoid numerical dissipation, the horizontal length of the relaxation source field zone is selected in the order of wavelengths. The results presented here and in the sequel correspond to delta-type forcing in the right-hand side of Equation (22), and thus, are normalized with respect to the monopole/dipole source amplitude.
In the case of the isovelocity environment, a comparison of the computed solution using the present FDTD-PML scheme and the Helmholtz equation is presented in Figures 8-10. Numerical results concerning the acoustic field excited by the monopole/ dipole sources at f = 5 Hz are obtained using a r − z grid consisting of 300 by 150 points, and ∆t = 0.005 s corresponding to 40 time steps in a period, which provides a Courant number of C = c/(∆x/∆t) = 0.75, that is sufficient for convergence and stability of the present FDTD method.
trate the time evolution of the solution, in Figure 7 the acoustic field generated by monopole source is plotted at time instants after 1, 3, 7, and 14 periods from the beginn of the calculation. Subsequently, in Figures 8-10 the periodic results obtained after 15 riods are compared against the solutions of the Helmholtz equation, for the monopole dipole sources, respectively. In order to avoid numerical dissipation, the horizontal len of the relaxation source field zone is selected in the order of wavelengths. The results sented here and in the sequel correspond to delta-type forcing in the right-hand sid Equation (22), and thus, are normalized with respect to the monopole/dipole source plitude. trate the time evolution of the solution, in Figure 7 the acoustic field generated by the monopole source is plotted at time instants after 1, 3, 7, and 14 periods from the beginning of the calculation. Subsequently, in Figures 8-10 the periodic results obtained after 15 periods are compared against the solutions of the Helmholtz equation, for the monopole and dipole sources, respectively. In order to avoid numerical dissipation, the horizontal length of the relaxation source field zone is selected in the order of wavelengths. The results presented here and in the sequel correspond to delta-type forcing in the right-hand side of Equation (22), and thus, are normalized with respect to the monopole/dipole source amplitude.       Similar results obtained in the case of underwater acoustic propagation in the sea acoustic waveguide characterized the vertical sound-speed profile, shown in Figure 8, are presented for the monopole, horizontal and vertical dipole source in Figures 11-13, respectively, along with comparison with the results from the frequency domain where the VEP is solved with a finite difference scheme. It is seen in the aforementioned figures that the propagation pattern is modified by the presence of additional boundaries and refraction phenomena that manifest at higher frequencies. In particular, the effect of the free surface and hard bottom modifies the noise pattern significantly.
Similar results obtained in the case of underwater acoustic propagation in the sea acoustic waveguide characterized the vertical sound-speed profile, shown in Figure 8, are presented for the monopole, horizontal and vertical dipole source in Figures 11-13, respectively, along with comparison with the results from the frequency domain where the VEP is solved with a finite difference scheme. It is seen in the aforementioned figures that the propagation pattern is modified by the presence of additional boundaries and refraction phenomena that manifest at higher frequencies. In particular, the effect of the free surface and hard bottom modifies the noise pattern significantly.  acoustic waveguide characterized the vertical sound-speed profile, shown in Figure 8, are presented for the monopole, horizontal and vertical dipole source in Figures 11-13, respectively, along with comparison with the results from the frequency domain where the VEP is solved with a finite difference scheme. It is seen in the aforementioned figures that the propagation pattern is modified by the presence of additional boundaries and refraction phenomena that manifest at higher frequencies. In particular, the effect of the free surface and hard bottom modifies the noise pattern significantly.

Extension to the Case of the Impedance-Type Seabed Model
Next, the present acoustic model is extended in the case of a more realistic sea bottom which is represented as an impedance boundary. In order to quantify the seabed acoustic parameter, the Pekeris waveguide consisting of a homogeneous layer overlying a homogeneous half-space is considered. In the latter two-layer medium, the sound velocity, attenuation, and density can be discontinuous from layer to layer. In this work, we will proceed to the formulation of an impedance-type seabed model using the Pekeris waveguide analysis; see [12]. In this case, an isovelocity water layer lies over an isovelocity half-

Extension to the Case of the Impedance-Type Seabed Model
Next, the present acoustic model is extended in the case of a more realistic sea bottom which is represented as an impedance boundary. In order to quantify the seabed acoustic parameter, the Pekeris waveguide consisting of a homogeneous layer overlying a homogeneous half-space is considered. In the latter two-layer medium, the sound velocity, attenuation, and density can be discontinuous from layer to layer. In this work, we will proceed to the formulation of an impedance-type seabed model using the Pekeris waveguide analysis; see [12]. In this case, an isovelocity water layer lies over an isovelocity half-space bottom and a boundary condition is applied in the interface between the layers. The general solution in the half space is given by: where and c b represent the sound speed in the seabed. Requiring a bounded solution, the above equation reduces to: horizontal, and vertical dipole sources, respectively, over an impedance seabed boundary. These expressions will be used to enforce near-source field data by applying, again, a relaxation scheme, as already discussed in Section 4.2. Next, the numerical results obtained using the present FDTD-PML numerical scheme for the case of an impedance-type BC for the modeling of the seabed effects in the acoustic propagation are presented for the monopole, horizontal and vertical dipole source in Figures 15-17 Thus, after deriving the eigenvalues I zn k , the above analytical expression obtained using the normal-mode expression of the Green's function and its derivatives in the , s s xz direction of the location of the source could be similarly exploited to model the monopole, horizontal, and vertical dipole sources, respectively, over an impedance seabed boundary. These expressions will be used to enforce near-source field data by applying, again, a relaxation scheme, as already discussed in Section 4.2.
Next, the numerical results obtained using the present FDTD-PML numerical scheme for the case of an impedance-type BC for the modeling of the seabed effects in the acoustic propagation are presented for the monopole, horizontal and vertical dipole source in Figures 15-17, respectively, and compared and verified with results from the solution of the Helmholtz equation. The effect of the bottom scattering field is clearly observed in contrast with the hard bottom case.    Finally, having obtained normal mode expressions of the Green's function for the impedance sea bottom G I and for the hard seabed G N , respectively, the results can be used to derive an alternative representation of the impedance sea bottom singularity in the form: G I = G N + p scat , where p scat is a regular scattering field excited by the sea bottom data.
As an example, the scattering field is presented in Figure 18 as obtained using p scat = G I − G N in the case of a monopole source located at z 0 = −300 m and with parametersÂ = 0.001, h = 3000 m, and c = 1500 m/s. Comparing this result with Figure 11, we observe that the scattering field due to the effect of the impedance seabed is an order of magnitude less than the source field and is maximized on the sea bottom. is a regular scattering field excited by the sea bottom data.
As an example, the scattering field is presented in Figure 18 as obtained using  Figure 11, we observe that the scattering field due to the effect of the impedance seabed is an order of magnitude less than the source field and is maximized on the sea bottom.
The above alternative representation will be exploited in future work for formulating and solving scattering problems in the inhomogeneous ocean waveguide characterized by an impedance seabed with the additional presence of obstacles modeling the effects of ocean vehicles or AUVs modifying the directivity pattern of the acoustic field generated by the monopole and dipole sources associated with the hydrodynamic flapping thruster loads. For this purpose, appropriate 3D BEM methods have been developed; see, e.g., [53].

Conclusions
The effects of noise generation and propagation by flapping thrusters used for small ships and AUV propulsion were studied. A new and efficient, low-cost numerical method based on finite differences and absorbing layer techniques was presented and validated for treating noise propagation in large domains and lengthy time simulations. The present The above alternative representation will be exploited in future work for formulating and solving scattering problems in the inhomogeneous ocean waveguide characterized by an impedance seabed with the additional presence of obstacles modeling the effects of ocean vehicles or AUVs modifying the directivity pattern of the acoustic field generated by the monopole and dipole sources associated with the hydrodynamic flapping thruster loads. For this purpose, appropriate 3D BEM methods have been developed; see, e.g., [53].

Conclusions
The effects of noise generation and propagation by flapping thrusters used for small ships and AUV propulsion were studied. A new and efficient, low-cost numerical method based on finite differences and absorbing layer techniques was presented and validated for treating noise propagation in large domains and lengthy time simulations. The present FDTD-PML numerical scheme extends the standard acoustic model from free space to the inhomogeneous ocean waveguide in a 3D cylindrical environment including freesurface and bottom-boundary effects, as well as refraction/diffraction phenomena due to sound-speed variations. For simplicity, the axisymmetric sea acoustic waveguide was considered. A Perfectly Matched Layer technique was employed for the treatment of the outgoing radiating field behavior in the open boundaries. For the verification of the present method, the numerical results in the case of fields excited by monopole and dipole sources were presented with a comparison against a frequency domain solution both for an isovelocity environment and variable sound-speed profile in the case of a hard seabed boundary. The acoustic model is further extended by employing an impedance-type seabed model, allowing for a more realistic approach of the sea bottom. The FDTD-PML numerical scheme is in good agreement in comparison with analytical solutions and exhibits low computational cost, which makes it appropriate for large domains and lengthy time simulations. Future work will be directed towards the modeling of more complex range-dependent waveguide configurations, as in the case of variable seabed profiles using a curvilinear FDTD-PML scheme on a structured grid derived using conformal mapping and differential system techniques. The Green's function for an impedance-type seabed will used for the implementation of 3D BEM treating noise propagation by flapping propulsors of ocean and underwater vehicles, providing useful information concerning the directivity pattern of the acoustic field generated by the hydrodynamic flapping thruster loads.
Author Contributions: This work was supervised by K.B. The numerical scheme was developed by I.M. The draft of the text was prepared by both authors and the numerical simulations were handled by I.M. All authors have read and agreed to the published version of the manuscript.

Funding:
The present work has been supported by the Seatech H2020 project. The APC have been offered by JMSE.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable.
Acknowledgments: The present work has been supported by Seatech H2020 project received funding from the European Union's Horizon 2020 research and innovation program under the grant agreement No 857840. The opinions expressed in this document reflect only the author's view and in no way reflect the European Commission's opinions. The European Commission is not responsible for any use that may be made of the information it contains.

Conflicts of Interest:
The authors declare no conflict of interest.