Free-Surface Effects on the Performance of Flapping-Foil Thruster for Augmenting Ship Propulsion in Waves

Flapping foils located beneath or to the side of the hull of the ship can be used as unsteady thrusters, augmenting ship propulsion in waves. The basic setup is composed of a horizontal wing, which undergoes an induced vertical motion due to the ship’s responses in waves, while the self-pitching motion of the wing is controlled. Flapping foil thrusters can achieve high level of thrust as indicated by measurements and numerical simulations. Due to the relatively small submergence of the above biomimetic ship thrusters, the free-surface effects become significant. In the present work, the effect of the free surface on the performance of flapping foil thruster is assessed by means of two in-house developed computational models. On one hand, a cost-effective time-domain boundary element method (BEM) solver exploiting parallel programming techniques and general purpose programming on graphics processing units (GPGPU) is employed, while on the other hand a higher fidelity RANSE finite volume solver implemented for high performance computing (HPC) is used, and comparative results are presented. BEM and RANSE calculations present quite similar trends with respect to mean submergence depth, presenting 12%, 28%, and 18% of differences concerning the mean values of lift, thrust, and moment coefficients, respectively. The latter differences become very small after enhancement of the BEM model to include viscous corrections. Useful information and data are derived supporting the design of the considered biomimetic thrusters, for moderate submergence depths and conditions characterized by minor flow separation effects.


Introduction
Biomimetic systems such as flapping foils, located on the hull of the ship, are examined for the exploitation of energy in wave induced motions by direct conversion to useful propulsive power. The kinematics of a flapping thruster involve two oscillatory motions heave and pitch with appropriately defined phase difference; see, e.g., [1][2][3]. With appropriate selection of the pivot point, the most energy demanding foil oscillatory motion is the heaving motion. On the other hand, the most frequent situation of ship operation is in a wavy sea condition, causing in many cases responses of significant amplitude. In particular, the vertical motion of the ship stations especially near the bow (and perhaps also at the stern) can be found exploitable to provide the vertical oscillation of an unsteady flapping foil propulsion system free of cost with concurrently enhancing ship stability and reducing unwanted responses. At the first stage of research, passive flapping wings-systems, attached beneath the ship hull, have been considered in order to produce thrust, exploiting energy that is stored in ship motions; see, e.g., [4,5].
More recently, the replacement of passively controlled systems by actively enforced foil self-pitching motion, based on the (irregular) history of wing's vertical oscillation, has been studied by various authors; see, e.g., [6,7]. In order to produce thrust, a flapping foil performs a complex motion that can be decomposed to a linear and a rotational oscillation, while it advances with the ship's forward speed. The vertical foil oscillation is provided by ship's irregular motion in waves and is obtained by the combined effect of vessel's heaving and pitching at the station alongside the ship where the biomimetic thruster is arranged. Moreover, in the context of the project BIOPROPSHIP [8] and in the works [6,9], a method for the coupling of ship dynamics with unsteady flapping wing hydrodynamics is presented as obtained by using linear seakeeping analysis in conjunction with unsteady lifting-line theory and nonlinear 3D panel methods. Numerical results reveal that significant wave energy can be extracted from ship motions by the flapping-foil thruster and could be used for augmenting ship propulsion in waves, with concurrently reducing ship's responses, enhancing the dynamic stability of the combined system with positive effects on safety, passenger and crew comfort, and general behavior in a seaway. In addition, it is shown that the actively controlled self-pitching motion of the thruster, which is tuned on the basis of the ship vertical motion at the station where the foil is arranged, is very efficient and is characterized by very low cost since the power required is found to be only a very small fraction of the developed thrust power.
Focusing on the importance of the effects associated with the interaction of the foil with the sea waves and the free-surface boundary, in [10] the authors developed a two-dimensional time-domain boundary element method (BEM) to study the performance of a hydrofoil that performs unsteady motions in the proximity of the free surface, in harmonic waves. For the treatment of the dynamical system, the BEM is coupled to the dynamical equations and is integrated in the time domain by appropriate schemes. The above method is validated demonstrated through comparison of calculations against linear theory, RANSE solvers, and experiments. In addition, the significance of the effects of the freely moving boundary over the foil is underlined, and it is shown that the system's performance in harmonic-wave conditions is very promising. Furthermore, aiming to investigate the performance of the combined system (ship in waves and flapping-foil thruster) in more realistic conditions representing specific sea states, in [9], the performance of the system is studied in irregular waves obtained through simulation from frequency spectra. The effect of various parameters of the ship-foil system operating in realistic sea conditions have been presented and discussed, and the results indicate that single or multiple arrangements of oscillating foils attached to the ship hull could enhance the performance in a broad range of parameters.
The increased interest in wave propulsion systems by flapping foils is indicated by recent research on the topic. In fact, results from experimental research in ship models equipped with a foil below the hull demonstrated 40-45% reduction of vessels' heaving and pitching responses with simultaneous 60% reduction of the resistance in waves [11]. The study of Eco-Ship with active foils for wave-empowered propulsion including free-running model test experiments [12], showed speed improvements up to 6.2%. Additionally, applications concerning wave propulsion systems for unmanned wave glide vehicle (UWGV) are presented in [13], and the design of a flapping-foil wave powered vessel using numerical and experimental techniques for wave propelled unmanned surface vehicles (USVs) is studied in [14] in order to control forward speed and recover wave energy. Moreover, in [15], a general cargo ship with retractable foils is studied, reporting important fuel savings at full speed. In [16], the propulsion characteristics of flexible biomimetic propulsors in regular waves are studied, demonstrating that flexibility could enhance the propulsive performance and increase the wave energy extraction; see also [17,18]. A more detailed review can be found in [19,20]. In particular, a 3D nonlinear BEM is proposed in [20] for the flapping foil beneath the free surface in oblique incident waves, and a computational code is developed exploiting parallel programming techniques and general purpose programming on graphics processing units (GPGPU). In complementarity to applications concerning wave propulsion, flapping foils are also studied for the development of hydrokinetic energy devices [21][22][23] and hybrid devices [24], exploiting combined wave and tidal energy resources.
In the present work, the effect of the free surface on the performance of flapping foil thruster in irregular wave conditions is assessed by means of two in-house developed computational models.
On one hand, a cost-effective time-domain BEM solver exploiting parallel programming techniques and general purpose programming on graphics processing units (GPGPU) is employed, while on the other hand, a higher fidelity RANSE finite volume solver implemented for high performance computing (HPC) is used. Comparative results are presented providing qualitative and quantitative insight on the effects of viscosity, even in cases where the foil operates beneath the free surface performing general motions. Useful information and data are derived exploited for the enhancement of the BEM model and supporting the design of the considered biomimetic thrusters, for moderate submergence depths and conditions characterized by minor flow separation effects.

Augmentation of Ship Propulsion in Waves by Flapping Thruster-Problem Definition
We consider a ship in head waves advancing at constant forward speed U; see Figure 1. In the present approach, we use the equations of motion of the ship, derived in the body-fixed frame of reference and linearized by assuming small oscillatory amplitudes and small-slope incident waves; see, e.g., [25]. Standard linear seakeeping method in the frequency domain is used, to calculate the coupled heaving (ξ 3 ) and pitching (ξ 5 ) responses of the system (ship and foil) in the vertical plane; see [6,26]. The flapping thruster is located at the bow, at forward station (x wing ), and operates in random wave conditions, gaining its vertical motion from the ship heave ξ 3 and pitch ξ 5 responses. Thus, the vertical heaving motion of the foil is ξ 3 − x wing ξ 5 and the controlled self-pitching motion is denoted by θ. In the present work, the effect of the free surface on the performance of flapping foil thruster in irregular wave conditions is assessed by means of two in-house developed computational models. On one hand, a cost-effective time-domain BEM solver exploiting parallel programming techniques and general purpose programming on graphics processing units (GPGPU) is employed, while on the other hand, a higher fidelity RANSE finite volume solver implemented for high performance computing (HPC) is used. Comparative results are presented providing qualitative and quantitative insight on the effects of viscosity, even in cases where the foil operates beneath the free surface performing general motions. Useful information and data are derived exploited for the enhancement of the BEM model and supporting the design of the considered biomimetic thrusters, for moderate submergence depths and conditions characterized by minor flow separation effects. Figure 1. Ship hull, advancing in head waves, equipped with a flapping thruster located at the bow, in a forward station (x wing ) augmenting ship propulsion in waves. The biomimetic thruster operates in random wave conditions, gaining its vertical motion from ship heave ξ 3 and pitch ξ 5 responses. The vertical heaving motion of the foil is ξ 3 − x wing ξ 5 and the controlled self-pitching motion is denoted by θ.

Ship Dynamics Coupled with Unsteady Flapping Foil Thruster
The system of the equations that describe the dynamics of the system is as follows [6]: where ξ 30 , ξ 50 are the complex amplitudes of ship heaving and pitching oscillations, D 53 = −ω 2 en (a 53 + I 53 ) + iω en b 53 + c 53 , D 55 = −ω 2 en (a 55 + I 55 ) + iω en b 55 + c 55 (3) a jk and b jk , j, k = 3, 5, are the added mass and damping coefficients, m is the total mass of the ship and the wing, and p = −iω en Um is a Coriolis term. Concerning the elements of the hydrostatic matrix, c 33 = ρgA WL , c 35 = c 53 = −ρg x f A WL and c 55 = m g GM L , where A WL is the waterline area, x f is the longitudinal center of flotation, GM L the longitudinal metacentric height, and g the acceleration of gravity. Additionally, I 55 = mR 2 yy and I 35 = I 53 = −m X G are inertia coefficients, where X G is the longitudinal position of the center of gravity and R yy the radius of gyration with respect to the transverse axis associated with the ship pitching motion. In this work, bow waves are considered for the excitation of motions, and the encounter frequency is ω en = ω + kU, where k = ω 2 /g is the wavenumber of the incident waves and ω is the absolute frequency. Additionally, with F j0 , j = 3, 5, we represent the Froude-Krylov and diffraction force and moment complex amplitudes, respectively.
In the above model, the coupling with the flapping thruster dynamics is achieved through the terms X j0 , j = 3, 5, that denote complex amplitudes of excitation due to the wing and are functions of heaving (ξ 3 ) and pitching (ξ 5 ) responses of the ship. The latter vertical force and moment consist of two parts: X j0 = X j0,A (ξ 3 , ξ 5 ) + X j0,B (ϕ INC ), one part dependent on the oscillatory ship amplitudes and one dependent on the incoming wave potential ϕ INC = Λ exp(kz + ikx), where Λ is a constant coefficient involving the incident wave height and frequency.
In the sequel, we will briefly describe the approach that is used in order to analyze the wing forces and moments and more details can be found in [9]. Beginning with the kinematics of the flapping foil, the effective angle of attack is approximately given by where θ(t) denotes the angle of rotation of the foil. In Equation (4) the effect of wing's vertical oscillation on the angle of attack is ε(t) ≈ 1 U dh(t) dt , where h(t) = −ξ 3 (t) + x wing ξ 5 (t) denotes the oscillatory motion of the ship at the longitudinal position x wing of the foil; see Figure 1. Following the works [6,27], wing self-pitching is selected to be a linear function of the angle ε, i.e., θ(t) = w ε(t), and the multiplier w stands for the pitch proportional control parameter. Therefore, we finally obtain the following expression for the angle of attack.
If we assume a large aspect ratio, unsteady lifting line theory can be used to analyze foil excitation. In the present work, under the additional assumption of small Strouhal numbers, the calculation of forces is based on a quasi-steady approximation and spanwise integration of sectional lift forces, resulting in the following expressions [6].
where c R denotes the root chord of the foil, S W is the wing planform area, AR its aspect ratio, ρ the water density, and h 0 and a 0 are the amplitude of heaving motion and the complex amplitude of the angle of attack, respectively. Furthermore, C 3D L , C 3D M , and C 2D L are the 3D and 2D lift and moment coefficients. For more details see, e.g., [6], Equations (14)- (17). Moreover, X 30,B (ϕ INC ) is the contribution of the incoming wave and G = χ πρU 2 S w AR AR+2 . In the latter expression, χ denotes an empirical correction factor calculated by comparing the quasi-steady lifting line estimations with a posteriori predictions by more accurate time domain panel methods for the responses of the oscillatory wing, following the same combined forward and oscillatory motion. Applying the aforementioned analysis in the dynamic system Equations (1)-(3), we introduce the following variations to the system coefficients at the left hand side, due to the operation of the flapping foil: wing Gc/ 4U 2 (8) Furthermore, we obtain the next set of expressions concerning the incident wave-field dependent forces: where the vertical location of the foil beneath the free surface is at z = −d, where d is the foil submergence depth and c the foil root chord length.

Flapping Thruster Performance in Irregular Waves
Using the responses of the system, the kinematics of the flapping wing are obtained, permitting the detailed formulation and study of its performance in random head waves. For this reason, the panel method developed in [10] is applied for the numerical simulation of the performance of hydrofoils at low submergence depth including the free surface effects. We use the assumption of waves with small amplitude in comparison with the wavelength and we work in the context of linear water wave theory. Realistic sea conditions are modeled by parametric frequency spectra and the realization of the free surface elevation and of other wave signals is implemented by using the random phase model. More details can be found in [9].
The Froude number based on the length of the ship is F n = U/ gL and the corresponding one based on the root chord c of the foil is F f oil = U/ √ gc. We denote with S(ω) the wave spectrum in the earth fixed frame of reference and we can obtain an expression for the spectrum according to a moving observer with the foil mean velocity U, as follows: S U (ω e ) = S(ω) dω dω e , with dω = dω en / 1 + 2U gY , and Y = tan h(kH) + kH sec h 2 (kH) (12) Then we use the response amplitude operator RAO = ξ 3 − x wing ξ 5 at the position x wing of the wing, in order to obtain the spectrum of wings vertical motion: The vertical motion of the flapping thruster is simulated by the random-phase model by discretizing the band of encounter frequencies into a large number of consecutive bins of width δω centered at discrete values ω e,n η(t) = n η n cos ω n t + kx wing + ε n , where η n = 2S U f oil (ω e,n )δω e (14) ε n are random phases that uniformly distributes in [0, 2π) ; see, e.g., [28]. In Equation (12), the effect of finite water depth H has been included, and the case of deep water conditions is modeled by using appropriate large values. Additionally, the rotational motion of the flapping wing, about its pivot axis is defined as negative at the clockwise direction. The pivot point is located at a distance x R = c/3 from the leading edge.
The self-pitching motion of the foil is set by a simple proportional control law based on its linear oscillation history: where 0 < w < 1 is the control parameter; see, e.g., [6,27].
As an example, in Figures 2 and 3 we present selected results in the case of a ship model of a fast ferry travelling at F n = U/ gL = 0.25 with a flapping thruster operating at the bow, as schematically shown in Figure 1. The ship length to beam ratio is L/B = 7.7, the beam to draft ratio B/T = 3.71, and the block coefficient C b = 0.45. Moreover, the waterplane area coefficient is C WL = 0.69, the midsection coefficient C m = 0.81, and the flotation-center x f /L = −0.05 (LCF aft midship). Additionally, the metacentric radius BM L /L = 1.8. The vertical position of the buoyancy center is KB/L = 0.021 (from BL), and its position on the horizontal axis is LCB = −0.019 L. of finite water depth H has been included, and the case of deep water conditions is modeled by using appropriate large values.
Additionally, the rotational motion of the flapping wing, about its pivot axis is defined as negative at the clockwise direction. The pivot point is located at a distance / 3 R x c = from the leading edge. The self-pitching motion of the foil is set by a simple proportional control law based on its linear oscillation history: where 0 1 w < < is the control parameter; see, e.g., [6,27].
As an example, in   We assume that the horizontal position of the gravity center is the same as the location of the center of buoyancy, i.e., X G /L = −0.019 (aft midship), Y G = 0, and KG/T = 0.2 (from BL). Furthermore, the longitudinal metacentric height is GM L ≈ BM L . Finally, the radii of gyration about the x-axis and y-axis, respectively, are taken R xx /B = 0.25, R yy /L = 0.2. The foil is located at a distance x wing /L = 0.6, with respect to the midship section, and vertically at submergence d/T = 1.71 below WL. The wing planform shape is orthogonal, its chord is c/B = 0.27, and its span to chord ratio is s/c = 5, and thus, its aspect ratio is AR = 5. The wing planform area is S w /BT = 1.35, and the wing sections are symmetrical NACA0012.
In particular, in Figure 2, the ship heave and pitch responses are presented against nondimensional wavelength λ/L, for Froude number. F n = U/ gL = 0.25 and head waves (β = 180 • ). Additionally, the system responses are compared with (bold line) and without (dashed line) the operation of the flapping wing thruster. It is seen that the operation of the foil at a corresponding Froude number F f oil = U/ √ gc = 1.35 produces a modification of the heave responses due to changes of the resonance and has a more signifficant effect on the pitch responses, especially as frequency becomes lower than the resonant value. The latter is due to the antipitching moment generated by the foil lift and the large arm (x wing /L = 0.6) due to the arrangement of the foil in front of the bow of the ship; see Figure 1. Next, in Figure 3, the frequency spectrum concerning the vertical foil motion is presented against the ones of the incident head waves in the earth-fixed system and in the ship frame of reference. The incident frequency spectrum is modeled by the Bretschneider model and the frequency spectrum for significant wave height H s /L = 0.03 and peak period T p U/L = 0.7, for head seas (β = 180 • ), is shown by dotted line. The corresponding wave spectrum in the ship frame of reference travelling at F n = 0.25 (solid line), is shown by using a solid line, and line and the spectrum of the foil's vertical motion by dashed line, respectively. Moreover, the ship and foil system RAO modulus at the position of the foil (x wing /L = 0.6) is shown by using a thick line. In the examined case, the encounter frequency spectrum is concentrated around the resonant behavior in the vertical ship motion at the bow, and we are able to observe in Figure 3 the magnification of the spectrum density around the peak value due to the large responses. The above effect is illustrated also by means of a realization of the vertical motion of the foil arranged in front of the bow of the ship for the above wave conditions and for a time duration of 25 peak periods as presented in Figure 4. The result is compared with the ship heaving motion shown in the same plot by using a dashed line, making clearer the enlargement of the foil vertical motion. of the foilʹs vertical motion by dashed line, respectively. Moreover, the ship and foil system RAO modulus at the position of the foil (xwing/L = 0.6) is shown by using a thick line. In the examined case, the encounter frequency spectrum is concentrated around the resonant behavior in the vertical ship motion at the bow, and we are able to observe in Figure 3 the magnification of the spectrum density around the peak value due to the large responses. The above effect is illustrated also by means of a realization of the vertical motion of the foil arranged in front of the bow of the ship for the above wave conditions and for a time duration of 25 peak periods as presented in Figure 4. The result is compared with the ship heaving motion shown in the same plot by using a dashed line, making clearer the enlargement of the foil vertical motion.

Flapping Thruster Performance in Irregular Waves
In Section 2, the augmentation of ship propulsion in waves using flapping foils is formulated and a mixed frequency-time domain method based on linearized ship dynamics coupled with unsteady lifting line theory is presented. The latter method is applied to obtain the responses of the ship-foil system in irregular waves corresponding to specific sea states. Using the responses of the system, the kinematics of the flapping thruster are obtained, permitting the detailed study of its performance in random head waves by cost-effective BEM and high fidelity RANSE solvers presented in more detail in this section.

Flapping Thruster Performance in Irregular Waves
In Section 2, the augmentation of ship propulsion in waves using flapping foils is formulated and a mixed frequency-time domain method based on linearized ship dynamics coupled with unsteady lifting line theory is presented. The latter method is applied to obtain the responses of the ship-foil system in irregular waves corresponding to specific sea states. Using the responses of the system, the kinematics of the flapping thruster are obtained, permitting the detailed study of its performance in random head waves by cost-effective BEM and high fidelity RANSE solvers presented in more detail in this section.

BEM Model of the Flapping Thruster in Waves
In this subsection, the BEM model, developed by [10], is applied to obtain the performance of flapping wing thruster operating in random head waves. The problem is treated in the time domain and the flapping hydrofoil is denoted with a closed boundary ∂D B (t) moving with respect to an earth fixed reference frame. The model provides the two-dimensional lifting flow characteristics taking into account interaction with the free surface, and three-dimensional corrections can be approximately taken into account by means of coefficients based on aspect ratio, as discussed in the previous section (see, e.g., Equation (6)).
In the context of linear wave theory, the total wave potential Φ T (x, y; t) is decomposed to the incident irregular wave potential Φ I (x, y; t), and its disturbance is due to the operation of the foil Φ(x, y; t). The problem has been reformulated with respect to the disturbance potential satisfying on the solid boundaries ∂D B the no entrance Neumann condition, and the corresponding bottom boundary condition in the case of finite water depth. Moreover, the linearized free surface boundary conditions are applied on the mean free-surface level. In the above, n denotes the normal to the boundary unit vector, directed to the interior of the domain D.
In order to treat lifting flow problems around bodies with sharp boundaries as the trailing edge of a hydrofoil, a vortex wake model must be used in order to model the generation of lift on the foil and the introduction of vorticity on the wake. In the present work, a curve of potential discontinuity is generated behind the foil and a simplified wake model is used. The trailing vorticity model as well as the implementation of a pressure-type Kutta condition are applied on the foil's trailing edge.
More details are provided in [10]. Following a direct approach based on application of Green's theorem, the following set of boundary integral equations (on body boundary and the free surface) is obtained, with unknowns the potential on the body boundary Φ B and on the free surface Φ F , as well as their normal derivatives where G(x 0 x) denotes the Green's function, which includes the mirror of the source singularity with respect to the bottom surface (y = −H): In the above equations, x 0 = (x 0 , y 0 ) denotes a control point, x 0 = (x, y) is the integration point, and x m = (x, −y − 2H) its image with respect the horizontal bottom.
Applying a low-order boundary element method, the body contour is approximated by a closed polygonal line of N B straight line elements, and also the boundaries of the free surface and of the vortex wake are represented by N F and N W (t) elements, respectively. It is noted here that the problem is treated in the time domain and thus the size vortex sheet emanating from the foil trailing edge increases in time, as it also happens with the number of wake elements. The various hydrodynamic quantities on the panels, such as the potential, its normal derivative, are assumed to be approximated by piecewise constant distributions. Moreover, a collocation scheme is used with control points at the center of each panel, and the following discretized form of Equation (18) is obtained.
where the vectors Φ B , Φ F , b, ∂Φ F ∂n , and µ W consist of hydrodynamic quantities at collocation points of each boundary, µ W is the dipole strength in the foil's wake, and the matrixes D, P, Z are calculated as functions of panel integrals, that in the case of low order panel methods can be obtained analytically. The detailed expressions are provided in [10]. Equation (20) is the discrete Dirichlet to Neumann (DtN) operator that relates the potential with its derivative on the boundary, and also includes the unknown value of the potential jump or the dipole intensity µ W1 on the wake at the vicinity of the trailing edge. Applying the discrete DtN map, Equation (20), in the discrete form of the pressure-type Kutta condition and the free surface conditions, a system of ODEs is finally derived governing the evolution of the dynamical system in time, within the context of our approximation, More details concerning the analytic structure of f(U) and the numerical method for the time integration, as well as for the numerical treatment of the horizontally infinite boundary can be found in [10].
After the solution has been obtained (at each time step), the pressure is calculated through the unsteady Bernoulli's theorem. Furthermore, forces and the moment are directly calculated through the integration of the instantaneous pressure acting on the solid boundary. The forces are corrected to include viscous effects according to the following empirical shear stress coefficient (see also [29]), as a function of Reynolds number Re = Uc/v (where v is the kinematic viscosity) and the angle of attack a, as follows c cor = 0.0858 In the present study, where Re is of the order of 10 6 , the parameter c a is calibrated by RANSE simulation of flapping foil in random motion in unbounded domain as described in Section 4.3, and the corrections are used successfully in all the cases where the foil operates with the effect of the free surface in (Sections 4.4 and 4.5) and also exploited for a corrected prediction in Section 4.2 where also incident waves are considered.

CFD model of the Flapping Thruster in Waves
In this section, the theoretical background of the CFD methodology used in the present work is described. The incompressible equations are solved using the artificial compressibility formulation, while for the two phase flows, an additional equation of the volume fraction (VOF) [30], is added to the coupled system of equations; see also [31]. The governing equations are written in the following integral form: In the above system of equations, Q = [p, u, a] T is the vector of the primitive variables, which contains the pressure, the velocity vector, and the volume fraction. The real and the fictitious time are denoted with t and τ, respectively. The inviscid fluxes are denoted with F c , while F v are the viscous fluxes, and S q are the source terms. The above system of equations is closed with the artificial compressibility parameter β, where it is assumed that there is a relation between density and pressure in pseudo-time β = ∂ρ ∂p τ .
Although the system of equations is casted in primitive form, for time true computations, the conservative variables, U = [0, ρu, a] T , are used to advance the solution in time. For this purpose, the Jacobian transformation matrix Γ e is introduced in Equation (23). Furthermore, in order to achieve better convergence characteristics, the preconditioner matrix Γ is used [32]. The preconditioner multiples the fictitious time derivatives and is also used in the inviscid flux evaluation to remove the dependency of the eigenvalues from the density, when multiphase flows are considered. The two matrices are given by The governing equations are discretized using the finite volume method. The discretized form of the system Equation (23) for a given control volume Ω I , is written as The residual R Ω is computed as a sum of the fluxes over the interfaces of the control volume with the addition of the source terms S q , which are considered to be constant in Ω I .
The inviscid fluxes are evaluated by solving the preconditioned local Riemann problem between the neighbors of the face, using the approximate Riemann solver [33]. For the reconstruction of the velocity field, a piecewise linear interpolation scheme is used. For the pressure field, a similar methodology is followed, but due to the discontinuity of the pressure gradient, a density weighted interpolation scheme is adopted [34]. Furthermore, for the reconstruction of the VOF field, a high order interface capturing scheme is used [35], in order to reduce the numerical diffusion near the free surface.
For the turbulence modeling, the two equations model k-ω SST [36] is used, while a buoyancy term is added to the kinetic energy equation to suppress the turbulent viscosity at the free surface [37].
MaPFlow implements the technique of the deforming grids, which is described in [38], in order to account flows with objects that perform small amplitude motions.
For unsteady simulations, an implicit second order backwards difference scheme (BDF) is used [39], along with a dual time-stepping technique in order to facilitate convergence. Finally, for the inversion of the implicit operator, the Gauss-Seidel iterative method is used combined with the Reverse Cuthill-Mckee reordering scheme.
The above methodology was implemented as an extension in the compressible solver MaPFlow, see [40,41]. MaPFlow solves the unsteady Reynolds averaged Navier-Stokes equations on unstructured grids, in a multiprocessor environment utilizing the MPI protocol and is developed in the National Technical University of Athens.

Grid Generation
In this section, the basic numerical parameters for the CFD methodology are described. Figure 5 illustrates the basic numerical setup for the CFD simulations. All the CFD simulations employ the volume of fluid (VOF) approach to resolve both the air and the water phase. For the mesh generation, the BETA CAE ANSA preprocessor was used. The mesh is unstructured and extends 160 chords in the chordwise and 70 chords in the vertical direction.

Grid Generation
In this section, the basic numerical parameters for the CFD methodology are described. Figure 5 illustrates the basic numerical setup for the CFD simulations. All the CFD simulations employ the volume of fluid (VOF) approach to resolve both the air and the water phase. For the mesh generation, the BETA CAE ANSA preprocessor was used. The mesh is unstructured and extends 160 chords in the chordwise and 70 chords in the vertical direction.  At the inlet and outlet boundaries, far field conditions are imposed while damping zones are employed to prevent numeral reflections from the outer boundaries reaching the domain of interest. Far-field conditions are also imposed on the top boundary (10 chords from the free surface) where only the air-phase is present, and similarly at the bottom is located 60 chords below the free surface to ensure a deep-water approximation. Finally, no-slip condition is employed for the NACA0012 hydrofoil. The body undergoes a combined heaving and pitching motion and consequently properly resolving the generated wake, which is essential to improve the accuracy of the numerical predictions. To that end, two refinement regions are defined in the hydrofoil wake. One refinement zone is at the near hydrofoil region spanning 3 chords (RF1) and another spanning 30 chords to resolve the far wake (RF2), as can be seen in Figure 6. The free surface is modeled though the VOF approach while the STACS scheme [35] is used for the reconstruction of the volume fraction. The computational mesh in the free surface region is structured (while the rest of the mesh remains unstructured) and the resolution was chosen to have 180 points per wavelength and 50 points per wave height (the maximum is considered). Regarding the boundary layer region, the first spacing is 1e-06 chords, ensuring a y+ below 1, while 600 nodes were placed around the hydrofoil. The boundary layer consists of 65 layers with a growth factor of 1.1. A snapshot of the computational mesh in the vicinity of the body is shown in Figure 7. At the inlet and outlet boundaries, far field conditions are imposed while damping zones are employed to prevent numeral reflections from the outer boundaries reaching the domain of interest. Far-field conditions are also imposed on the top boundary (10 chords from the free surface) where only the air-phase is present, and similarly at the bottom is located 60 chords below the free surface to ensure a deep-water approximation. Finally, no-slip condition is employed for the NACA0012 hydrofoil. The body undergoes a combined heaving and pitching motion and consequently properly resolving the generated wake, which is essential to improve the accuracy of the numerical predictions. To that end, two refinement regions are defined in the hydrofoil wake. One refinement zone is at the near hydrofoil region spanning 3 chords (RF1) and another spanning 30 chords to resolve the far wake (RF2), as can be seen in Figure 6. The free surface is modeled though the VOF approach while the STACS scheme [35] is used for the reconstruction of the volume fraction. The computational mesh in the free surface region is structured (while the rest of the mesh remains unstructured) and the resolution was chosen to have 180 points per wavelength and 50 points per wave height (the maximum is considered). Regarding the boundary layer region, the first spacing is 1e-06 chords, ensuring a y+ below 1, while 600 nodes were placed around the hydrofoil. The boundary layer consists of 65 layers with a growth factor of 1.1. A snapshot of the computational mesh in the vicinity of the body is shown in Figure 7.  In order to simulate the hydrofoil motion, the technique of deforming grids is utilized. In particular, the near grid follows the prescribed body motion, while the rest is adapted so that the   In order to simulate the hydrofoil motion, the technique of deforming grids is utilized. In particular, the near grid follows the prescribed body motion, while the rest is adapted so that the displacement of the nodes decays as the distance from the moving boundary increases. At each timestep, the grid velocities are obtained by differences based on three snapshots.
Three successively refined grids are considered. Details are presented in Table 1. For the refinement zones (RF1, RF2) the cell sizes are presented with respect to the chord length. Gridindependence and convergence of numerical results is discussed in Section 4.5. In order to simulate the hydrofoil motion, the technique of deforming grids is utilized. In particular, the near grid follows the prescribed body motion, while the rest is adapted so that the displacement of the nodes decays as the distance from the moving boundary increases. At each timestep, the grid velocities are obtained by differences based on three snapshots.
Three successively refined grids are considered. Details are presented in Table 1. For the refinement zones (RF1, RF2) the cell sizes are presented with respect to the chord length. Grid-independence and convergence of numerical results is discussed in Section 4.5.

Numerical Results and Discussion
In this section, a comparison between the BEM and the RANSE solvers is presented in a case of a steady foil beneath the free surface (Section 4.1). Then, the mixed frequency-time domain method presented in Section 2 is exploited to obtain the kinematics of the ship and wing in waves and study the performance of the flapping thruster using the BEM solver (Section 4.2). Results concerning the flapping foil at great submergence (without free-surface effects) are presented in Section 4.3, and close to the free surface in Section 4.4, as obtained by using the BEM and RANSE solvers. Finally, in Section 4.5 systematic investigation of free-surface, viscosity, and incident-wave effects is presented and discussed including convergence study of both solvers.

BEM and RANSE Calculations in Steady Conditions
Before we proceed to the unsteady calculations, we present a comparison between the BEM and the RANSE solvers in a case of a steady NACA0012 foil (Figure 8a) advancing below the free surface with velocity corresponding to foil Froude number F f oil = U/ √ gc = 1.35, with angle of attack a = 5 and submergence depth d/c = 1.71. The rotation axis is located on the chord at a distance X R /c = 0.33 from the leading edge. This case corresponds to the wave resistance problem of a submerged hydrofoil steadily travelling near the calm free surface. The foil accelerates gradually until steady conditions are obtained, which is achieved after travelling a length of more than 100 wavelengths, which as predicted by linear theory in this case is λ/c = 2πF 2 f oil . Comparison of simulations is presented in Figure 8, and it is observed that the calculated wavelength by both methods is in agreement with the prediction by the dispersion relation of linear theory. Calculated free surface elevation is comparatively shown in Figure 8b,c, where also the position of the foil is indicated at x/c = 0, d/c = 1.71 Furthermore, in Figure 8d,e, the calculated pressure coefficient on the foil c p = p − p atm /0.5ρU 2 is shown, and the results by both methods indicate very good agreement. Calculations have been done by using pitch control parameter 0.5 w = . In particular, in Figure  9, the ship travelling in irregular waves is shown at five instants within a time interval corresponding to one modal period. The foil is located in front of the bow of the ship, at distance xwing /L = 0.6 with respect to the midship section of the ship. In the same figure, the foil wake is plotted, with the calculated dipole intensity (potential jump) on the vortex sheet, which is illustrated by using arrows normal to the wake curve with length proportional to the local dipole strength. The latter result is associated with the memory effect of the generated lifting flow around the flapping foil operating in random incident waves.
Moreover, in the right subplots of Figure 9, the instantaneous distribution of the pressure coefficient on the foil, at the same time instants, as calculated by the present method. From the calculated pressure distributions, lift and thrust components are obtained at each time step by integration. In order to illustrate the relative magnitude of the incident and disturbance flow generated by the flapping motion of the foil in waves, for the same as before wave conditions and hydrofoil data, the calculated free-surface elevation normalized with respect to the significant wave height plotted in Figure 10 over the horizontal domain, at some instances during two modal periods.
Furthermore, we present in Figure 11 the time series of several calculated quantities concerning thrust production by the examined system, as calculated by the present BEM for a time interval from 3 to 18 peak periods. The dynamic system is integrated starting from rest, and the transient effects have been died out in the first three periods. The responses of the examined ship and flapping foil system in waves are calculated by both the present BEM without (solid line) and with viscous corrections.

Foil Motion and BEM Calculations in Waves
Next, unsteady numerical calculations are presented in Figure 9, considering the performance of the flapping thruster in the bow of the ship subjected to vertical oscillations due to ship responses in head waves and at the same time performing self-pitching oscillations about its pivot axis at x R /c = 0.33 with controllable amplitude based on its vertical motion, as described in previous Sections 2 and 3.1. The flapping foil is located at the same mean submergence d/c = 1.71, advancing with the same as before velocity F f oil = 1.35 F ship = 0.25 , in irregular waves (H s /c = 0.86 and Strouhal number ω p,en H s /(2πU) = 0.067, where ω p,en is the frequency of encounter that corresponds to the peak period of the spectrum; see Section 2.
Calculations have been done by using pitch control parameter w = 0.5. In particular, in Figure 9, the ship travelling in irregular waves is shown at five instants within a time interval corresponding to one modal period. The foil is located in front of the bow of the ship, at distance x wing /L = 0.6 with respect to the midship section of the ship. In the same figure, the foil wake is plotted, with the calculated dipole intensity (potential jump) on the vortex sheet, which is illustrated by using arrows normal to the wake curve with length proportional to the local dipole strength. The latter result is associated with the memory effect of the generated lifting flow around the flapping foil operating in random incident waves.  Figure 3). Results plotted at various instances in one modal period. In the right subplots, the distribution of the unsteady pressure coefficient on the flapping foil, located at the bow, is plotted. The unit length of both horizontal and vertical axes is equal to L/10.  Figure 3). Results plotted at various instances in one modal period. In the right subplots, the distribution of the unsteady pressure coefficient on the flapping foil, located at the bow, is plotted. The unit length of both horizontal and vertical axes is equal to L/10. Moreover, in the right subplots of Figure 9, the instantaneous distribution of the pressure coefficient on the foil, at the same time instants, as calculated by the present method. From the calculated pressure distributions, lift and thrust components are obtained at each time step by integration. In order to illustrate the relative magnitude of the incident and disturbance flow generated by the flapping motion of the foil in waves, for the same as before wave conditions and hydrofoil data, the calculated free-surface elevation normalized with respect to the significant wave height plotted in Figure 10 over the horizontal domain, at some instances during two modal periods.  ), and slightly reduce the corresponding rms value from 5.2 to 4.9 . In general, the angle of attack is an indicator of the importance of viscous effects, which remain at moderate levels, rendering cost-effective potential methods capable to provide accurate predictions for the design of flapping foil thrusters as will be illustrated in the sequel through comparison against RANSE calculations in Sections 4.3-4.5.
In Figure 11b, we present the time history of foilʹs heaving motion, as produced by combined oscillatory heaving and pitching motion of the coupled ship-flapping foil system, together with the generated sectional lift coefficient We observe that significant amplitudes of the lift force are produced. Moreover, the phase lag between foilʹs motion and lift force is approximately 180 ; Therefore, the generated lift acts as a restoring force, reducing the responses of a ship equipped with flapping hydrofoils. In Figure 11c, the dynamic evolution of sectional thrust coefficient Furthermore, we present in Figure 11 the time series of several calculated quantities concerning thrust production by the examined system, as calculated by the present BEM for a time interval from 3 to 18 peak periods. The dynamic system is integrated starting from rest, and the transient effects have been died out in the first three periods. The responses of the examined ship and flapping foil system in waves are calculated by both the present BEM without (solid line) and with viscous corrections.
In particular, in Figure 11a, the time history of the angle of attack is presented excluding (a) and including (a w ) the effects of the incident wave field. The angle is calculated using Equation (5) and the effects of the free surface, the disturbance of the foil and the wake are approximately excluded from the calculation. The effect of the wave is either to increase or decrease the temporal value of the angle of attack by a small amount, resulting in a maximum value of a w = 15 instead of a = 12.9 (at t/T p = 9.7), and slightly reduce the corresponding rms value from 5.2 to 4.9. In general, the angle of attack is an indicator of the importance of viscous effects, which remain at moderate levels, rendering cost-effective potential methods capable to provide accurate predictions for the design of flapping foil thrusters as will be illustrated in the sequel through comparison against RANSE calculations in Sections 4.3-4.5.  In Figure 11b, we present the time history of foil's heaving motion, as produced by combined oscillatory heaving and pitching motion of the coupled ship-flapping foil system, together with the generated sectional lift coefficient C L = F Y /0.5ρU 2 c. We observe that significant amplitudes of the lift force are produced. Moreover, the phase lag between foil's motion and lift force is approximately 180; Therefore, the generated lift acts as a restoring force, reducing the responses of a ship equipped with flapping hydrofoils. In Figure 11c, the dynamic evolution of sectional thrust coefficient C T = −F x /0.5ρU 2 c is shown, in the same time interval. We observe in this subplot that the thrust oscillations are in the interval −0.06 ≤ C T ≤ 0.49, with an average value of 0.051 for the purely potential calculation, which is indicated by a horizontal black solid line and 0.033 for the corrected results indicated with a horizontal red dashed line. Although viscous effects cause almost a 35% decrease in mean value, the corresponding decrease is only 3.7% as concerns the peak value. In the examined case, the thrust production compared to the calm water resistance is found to reach levels of the order of 30%, which is considered to be important.
Additionally, the time history of moment required for the active pitching motion of the foil is presented in Figure 11d in terms of sectional moment coefficient C M = M/0.5ρU 2 c 2 . As expected, viscous corrections increase the required torque input for the system to operate. Moreover, in Figure 11e, the power extracted by the examined system from the waves C PT = C T is compared against the corresponding power that is necessary for the self-pitching motion of the foil C Pθ = P θ /0.5ρU 3 c, actually for tuning the instantaneous angle of attack in order to produce positive thrust. We observe that the latter quantity is practically insignificant and that the viscous effects modify slightly the predictions by reducing C PT = C T and increasing C Pθ . In the present section, the wave effects in the performance of the foil have been studied by means of a BEM and the effects of viscosity have been considered using empirical corrections. In the following subsections, the effects of viscosity and the free-surface boundary are also investigated by means of the CFD model and comparison between the BEM and RANSE are presented and discussed. First, in Section 4.3, the effects of the free-surface are neglected in order to derive calibration factors for the BEM model for viscous corrections. Subsequently, in Section 4.4, both models are applied and the predictions including the free-surface effects are compared.

BEM and RANSE Calculations of Flapping Foil in Unbounded Domain-Calibration of Viscous Corrections
In the present subsection, the effects of the wave and the free surface boundary are excluded from the modeling in order to derive estimates of the effects of viscosity for the flapping thruster. BEM and RANS calculations are compared to identify the range of applicability and the limitations of each scheme and to calibrate the correction formula (Equation (22)). To be more specific, in Figures 12-14, simulations with both numerical methods are presented for a flapping foil in heaving motion, pitching motion, and angle of attack as in Figure 11.
In Figure 14, we present the loads of the system in waves in unbounded domain, as calculated by the present BEM including corrections and RANSE. In subplots (a), (c), and (e), the plots of C L , C T , C M histories are presented, respectively, whereas in subplots (b), (d), and (f), details during the 16th and 17th period are provided for comparison between the two methods.
The lift oscillations are in the interval −1.24 ≤ C L ≤ 1.23, with an rms value of 0.45 for the BEM calculation and 0.39 for RANSE. Although viscous effects cause almost a 13% decrease in rms value of lift, the reduction is only 4.9% of the maximum peak value. The thrust oscillations are in the interval −0.023 ≤ C T ≤ 0.51, with a mean value of 0.066 as obtained from BEM and 0.048 from RANSE. Although viscous effects cause almost a 27% decrease in mean value of thrust, the reduction is only 3.5% of the maximum peak value. The moment oscillations are in the interval −0.13 ≤ C M ≤ 0.13 from BEM reaching the maximum value of 0.094, with mean value of 0.032. The corresponding value for RANSE is 0.038. Although viscous effects cause almost a 19% increase in rms value of foil self-pitching moment, that increase is only 6.4% of the maximum peak value. The corrected BEM results with appropriate selection of coefficient c a provide compatible predictions with CFD model.

BEM and RANSE Calculation of Flapping Foil Beneath the Free Surface
Following the previous analysis, in this section the effect of the free surface is now taken into account. The flapping foil undergoes the same heaving and pitching motion as before, advancing with the same forward speed as in Section 4.2. A comparison between BEM and RANS numerical predictions is presented in terms of foil loading as well as the foil-induced wave system.
In Figure 15, three snapshots of the flapping foil numerical predictions can be seen. BEM and RANSE simulations are presented during a time interval of one peak period. As in the previous subsection, BEM results concerning the dipole intensity in the vortex wake are shown, while for the CFD predictions vorticity contours are plotted.

BEM and RANSE Calculation of Flapping Foil Beneath the Free Surface
Following the previous analysis, in this section the effect of the free surface is now taken into account. The flapping foil undergoes the same heaving and pitching motion as before, advancing with the same forward speed F f oil = 1.35. However, now it is considered completely submerged at a depth d/c = 1.71 as in Section 4.2. A comparison between BEM and RANS numerical predictions is presented in terms of foil loading as well as the foil-induced wave system.
In Figure 15, three snapshots of the flapping foil numerical predictions can be seen. BEM and RANSE simulations are presented during a time interval of one peak period. As in the previous subsection, BEM results concerning the dipole intensity in the vortex wake are shown, while for the CFD predictions vorticity contours are plotted.
Comparing the wake evolution between the two methods, it is clear that there is good qualitative agreement between the two solvers. Indeed, as the instantaneous pressure coefficients suggest (Figure 15), results are consistent between the two numerical methods. As stated in the previous section, apart from the absence of viscosity in the BEM predictions, the discrepancies between the two methods are due to the linearization of the wake vortex dynamics. When the free surfaces are considered, another source of error between the two methods emerges since the way the free surface is modeled is completely different. For the RANSE simulations the VOF approach is employed and thus the free surface is not explicitly defined. On the other hand, in the BEM approach, the free surface elevation is calculated as an explicit function of horizontal space variable x, also linearization of free-surface boundary conditions has been applied.
In Figure 16, the free surface elevation as well as the pressure distribution are plotted seven times between the 4th and the 6th period. It is observed that BEM and RANSE are in very good agreement both for the predicted free surface elevation as well as for the instantaneous pressure coefficient. It is noted that, contrary to the RANSE simulations (where the foil is stationary and the flow comes with U ∞ = −U), in the case BEM, the problem is solved in the earth-fixed frame of reference (the foil is travelling with forward speed U). However, the results in Figure 15 suggest that the predicted surface elevation in the near-foil regime is in very good agreement between the two methods. The same can be concluded by comparing the pressure coefficients presented in Figure 16c,d, from which it is evident that both solvers predict the very similar results for c P .

Effects of Waves, Free Surface, Viscosity, and Numerical Convergence
In this section, detailed investigation of viscosity, free-surface, and wave effects in a wide range of mean submergence of the flapping thruster is performed. Moreover, a detailed convergence study of the two numerical schemes is presented. In all cases, a NACA0012 heaving foil is considered in forward motion that corresponds to  Finally, in Figure 17, we present the evolution of the forces and moment acting of the flapping foil operating beneath the free surface for the simulation of Figure 15. Good comparison between the present BEM and RANSE results is obtained. Moreover, calculations with the corrected BEM method are provided in the same figure. The viscous corrections are calibrated using the results of Section 4.3 for time integrated values (average and RMS). Again, simulation is over 18 peak periods, and the time-integrated quantities (averaged and RMS values calculated excluding the first 3 periods) are presented with horizontal lines. In Figure 17a,c,e, the time history of force coefficients is presented, whereas in Figure 17b,d,f, details during the time interval between 16th and 17th period are provided.
The sectional lift coefficient oscillations are found to be in the interval −1.23 ≤ C L ≤ 1.26 that is slightly wider than the corresponding interval calculated in the previous section in unbounded domain (−1.  The self-pitching moment oscillations are in the interval −0.13 ≤ C M ≤ 0.14 from BEM, reaching a maximum value of 0.097 and rms value 0.032, while the prediction from CFD analysis is 0.039. The corresponding mean values in unbounded domain were found to be very similar. Although viscous effects cause almost a 22% increase in rms value, the increase is only 7.2% of the maximum peak value. In the present case, where free-surface effects are considered, the corrected BEM predictions differ only at the third significant digit from the corresponding CFD. Moreover, the empirical corrections reduce the discrepancies at the peaks between BEM and RANSE calculations. Therefore, in the present case of Froude number F f oil = 1.35 and mean submergence d/c = 1.71, free-surface boundary does not affect significantly the calculations, and calibration using only one unbounded domain simulation is sufficient.

Effects of Waves, Free Surface, Viscosity, and Numerical Convergence
In this section, detailed investigation of viscosity, free-surface, and wave effects in a wide range of mean submergence of the flapping thruster is performed. Moreover, a detailed convergence study of the two numerical schemes is presented. In all cases, a NACA0012 heaving foil is considered in forward motion that corresponds to F f oil = 1.35 simultaneously pitching about a pivot axis located on the chord at distance x R /c = 0.33 from the leading edge, resulting in angle of attack a(t) and kinematics given in Figure 11.
The calculations in the unbounded domain are labeled as "inf" in Figure 18, while the calculations of the foil beneath the free surface without incident waves are denoted with "fs" and the calculations in irregular waves corresponding to H s /c = 0.86 and Strouhal number ω p,en H s /(2πU) = 0.067 are denoted with label "wave". The mean submergence parameter d/c varies from 1.14 to 18.29 and even in the "wave" calculations the same foil kinematics are used (ignoring the effect of the varying foil submergence to the ship-foil RAO). In all cases, BEM calculation with viscous corrections are also presented. The latter corrections are calibrated as described in Section 4.3 using data from unbounded domain BEM and RANSE simulations.
In Figure 18a-c, results concerning the lift, thrust, and moment coefficients are presented as functions of the mean submergence parameter d/c. Concerning viscosity effects, BEM and RANSE "fs" calculations present a similar trend with respect to mean submergence parameter d/c with approximately constant deviation −12%, −28%, and +18% for the mean values of C L , C T , and C M , respectively. It is observed that the corrected BEM calculations are in very good qualitative and quantitative agreement with the RANSE calculations. The latter discrepancy is not significant compared to the total magnitude of the time signals; for more details see the discussion of Figures 11, 14 and 17 in previous sections. Numerical results indicate thrust production from the waves of the order of 20-30% of the calm water resistance in the examined cases.
As concerns the free surface effects, we observe that for foil submergence greater than 1-2 chord lengths they are not significant. To be more specific, the maximum increase with respect to the unbounded domain calculations is observed about d/c = 1.71 and is below 2% for all the load coefficients. Moreover, the maximum decrease is observed about d/c = 4.57 and is below 1.1% for the cases studied. However, it is noted that for the exact location of the maximum and minimum values more simulations with finer resolution of the d/c variation should be performed, and this is left as a subject for future work. The appearance of local extremes can be explained considering that the wave-making effects increase energy transfer from the foil to wave flow, reducing the load coefficients as the foil comes closer to the free surface. On the other hand, the free-surface boundary condition shows reverse effect as far as wave breaking does not occur. In Figure 18a  In addition, in the same figure results concerning load coefficients in presence of waves are plotted. With bold solid lines, the BEM calculations are indicated, and with bold dashed lines the corrected BEM ones. The comparison with RANSE calculations is left to be examined in future work. It is worth mentioning that calculations are in good agreement with the predictions in unbounded domain as the foil submergence increases. In general, the effect of the incident wave field in the present cases (ignoring the effect to the ship-foil RAO) is to reduce the load coefficients with maximum reduction 10%, 24%, and 18% for the mean C L , C T , and C M , respectively. The difference is not significant compared to the total magnitude of the time signals; however, it should be considered in the modeling. This could be approximately calculated by considering wave-like background incident velocity field in the form of an unsteady gust.
Next, a convergence study has been performed for both numerical schemes in terms of space and time resolution for all load coefficients and the whole range of values of foil submergence. The error indicator used is (C F − C F,dense )/C F,inf , where the index F stands for lift, thrust, and moment, C F,dense denotes the value of the coefficient obtained using the finer grid, and C F,inf is the converged value in unbounded domain. All BEM calculations until now were implemented using the following numerical parameters, which were found enough for convergence of results: N B = 240, dt/T P = 0.58 %, N F = 814 that corresponds to dl/λ = 2 %, where dl is the panel size on the free surface and λ the wavelength of the disturbance wave as predicted by the linear theory. Additionally, all RANSE calculations are performed using grid 3 (see Table 1) and dt/T P = 0.036%.
The most demanding goal was to achieve convergence with relative error significantly less than 1% in order to resolve the free-surface effects especially for the demanding RANSE simulations. This was achieved using the GRNET ARIS HPC infrastructure. In cases where the above goal is not achieved, the systematic study for the entire range of d/c and the good agreement between the corrected BEM and the RANSE calculations ensure that the trend of free-surface effects with respect to the submergence is correctly resolved. In BEM calculations, the number of panels in the body and the timestep are connected in order to ensure reasonable panel size ratio in the foil and wake panels near the trailing edge. Therefore, BEM convergence is presented only with respect to N B . In RANSE simulations, convergence is examined by using three different timesteps and grid 3.
Turning into more details, we present in Figure 19 a convergence study for a flapping foil in unbounded domain for the case studied in Section 4.3. BEM convergence results with respect to the number of panels on the foil N B and the size of the time step are plotted in Figure 19a-c. The relevant errors in all cases are reasonably small and mean sectional lift coefficient converges faster while the moment coefficient is found to be the most demanding. Moreover, a grid independence study of RANSE calculations is presented in Figure 19d-f using different grids (see Table 1). In this case, the level of error remains below 0.8% for all cases for grid 3. It was decided to present the finer grid predictions in the preceded analysis for the reasons explained before and since the wake is better resolved.
Furthermore, in Figure 20 we present a convergence study for a flapping foil beneath the free surface for the case studied in Section 4.4. BEM convergence results are plotted in Figure 20a-c with respect to the number of panels on the foil N B and the size of the time step and with respect to the number of panels on the free-surface boundary N F in Figure 20d-f, respectively. From the first case in Figure 20a-c, we conclude that convergence characteristics are similar to the unbounded domain study as expected with slightly faster convergence. From the second case in Figure 20d-f, it is seen that for the specific submergence ratio d/c = 1.71 the relative error reduces below 0.5% even for the coarser grids. However, we choose N F = 814 to resolve adequately the most demanding case of d/c = 1.14. A grid independence study of RANSE calculations is presented in Figure 20g-i, using different grids according to Table 1. The relative errors are very small, but the convergence is found to be slightly more demanding in comparison with the unbounded domain configuration due to the additional requirements imposed by the free-surface.  Table 1). Figure 20 we present a convergence study for a flapping foil beneath the free surface for the case studied in Section 4.4. BEM convergence results are plotted in Figure 20a,b,c with respect to the number of panels on the foil B N and the size of the time step and with respect to the number of panels on the free-surface boundary F N in Figure 20d,e,f, respectively. From the first case in Figure 20a,b,c, we conclude that convergence characteristics are similar to the unbounded domain study as expected with slightly faster convergence. From the second case in Figure 20d,e,f, it is seen that for the specific submergence ratio / . A grid independence study of RANSE calculations is presented in Figure 20g,h,i, using different grids according to Table 1. The relative errors are very small, but the convergence is found to be slightly more demanding in comparison with the unbounded domain configuration due to the additional requirements imposed by the free-surface.

Furthermore, in
Finally, in Figure 21, a convergence study for a flapping foil in waves for the case studied in Section 4.2 is presented. BEM convergence with respect to the number of panels on the foil B N and the size of the time step are presented in Figure 21a,b,c and with respect to number of panels on the free-surface boundary F N in Figure 21d,e,f. The findings are similar with those from the preceded figures; however, as explained in the beginning of the present section, the effect of waves causes a variation of reduction by 10% , 24% , and 18% for the mean sectional values of L C , T C , and M C , respectively, and not below 2% as in the case of free-surface effects. Therefore, if the interest is in the design of the flapping foil device extracting energy from the waves in mean submergence larger than 1 chord length and not to the resolution of the details of the free-surface effects, a coarser spatial and temporal grid could be adequate. The enhancement of the present method with incorporation of nonlinear effects associated with the free surface dynamics and the leading edge flow separation, as well as full treatment of 3D effects, are important future contributions, supporting also the design of the considered systems; see, e.g., [20,43].  Table 1).  BEM results are presented using solid lines and BEM corrected with dashed lines, respectively. Grid independence study of RANSE (g,h,i) with different grids (see Table 1).

Figure 20.
Convergence study for a flapping foil thruster beneath the free surface for the case studied in Section 4.4. (a-c) BEM convergence with respect to the number of panels on the foil N B and the size of the time step, and (d-f) with respect to number of panels on the free-surface boundary N F . BEM results are presented using solid lines and BEM corrected with dashed lines, respectively. Grid independence study of RANSE (g-i) with different grids (see Table 1).
Finally, in Figure 21, a convergence study for a flapping foil in waves for the case studied in Section 4.2 is presented. BEM convergence with respect to the number of panels on the foil N B and the size of the time step are presented in Figure 21a-c and with respect to number of panels on the free-surface boundary N F in Figure 21d-f. The findings are similar with those from the preceded figures; however, as explained in the beginning of the present section, the effect of waves causes a variation of reduction by 10%, 24%, and 18% for the mean sectional values of C L , C T , and C M , respectively, and not below 2% as in the case of free-surface effects. Therefore, if the interest is in the design of the flapping foil device extracting energy from the waves in mean submergence larger than 1 chord length and not to the resolution of the details of the free-surface effects, a coarser spatial and temporal grid could be adequate. The enhancement of the present method with incorporation of nonlinear effects associated with the free surface dynamics and the leading edge flow separation, as well as full treatment of 3D effects, are important future contributions, supporting also the design of the considered systems; see, e.g., [20,43]. BEM results are presented using solid lines and BEM corrected with dashed lines, respectively. Grid independence study of RANSE (g,h,i) with different grids (see Table 1).  (a-c) BEM convergence with respect to the number of panels on the foil N B and the size of the time step, and (d-f) with respect to number of panels on the free-surface boundary N F . BEM results are presented using solid lines and BEM corrected with dashed lines, respectively.

Conclusions
The concept of biomimetic flapping-foil thrusters for energy extraction and augmentation of ship prolusion has been further explored, focusing on the effects of free surface, incident wave, and viscosity. To this end, in-house GPU-BEM and HPC-RANSE solvers have been developed and results for cases of flapping foil in unbounded domain, beneath the free surface and in waves have been presented and discussed. For the determination of foil kinematics, a mixed frequency-time domain method based on linearized ship dynamics coupled with unsteady lifting line theory is applied, obtaining the response of the ship-foil system in irregular waves that correspond to specific sea states. After we have obtained the kinematics of the wing, an unsteady time domain boundary element method has been applied for the hydrodynamic analysis of flapping foils beneath the free surface in the presence of irregular incident waves. Numerical results concerning the thrust production have been presented, indicating significant energy extraction from the waves for the augmentation of ship's overall propulsion and enhancement of dynamic stability.
Free-surface, incident wave, and viscosity effects have been thoroughly studied indicating that in cases with foil submergence greater than a few chord lengths, the free-surface effects are minor. Moreover, BEM and RANSE calculations present quite similar trends with respect to mean submergence depth, presenting 12%, 28%, and 18% differences concerning the mean values of lift, thrust, and moment coefficients, respectively. From the simulations in unbounded domain, viscous corrections of the potential flow model have been derived, rendering the enhanced BEM solver cost-effective and useful for the systematic investigation and the optimization of the system operating relatively near the free surface in waves.
Future work is focused on the enrichment of the present models by the incorporation of nonlinear phenomena associated with the free surface dynamics and flow separation. Furthermore, full treatment of 3D effects including swept wings of low aspect ratio and effects of directional seas, as well as experimental verification of the examined ship flapping foil thruster are important directions for future research.
Author Contributions: This work is supervised by K.A.B. The BEM numerical scheme is developed and implemented in GPGPU by E.S.F. and the RANSE numerical scheme by G.P.P. The draft of the text is prepared by all authors and the numerical simulations were handled by E.S.F. and G.P.P. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Value of the load coefficient at the most refined case C F,inf Converged value of the load coefficient in unbounded domain x wing Horizontal location of the foil µ Surface dipole intensity or potential jump