A Non-Linear BEM–FEM Coupled Scheme for the Performance of Flexible Flapping-Foil Thrusters

Recent studies indicate that nature-inspired thrusters based on flexible oscillating foils show enhanced propulsive performance. However, understanding the underlying physics of the fluid–structure interaction (FSI) is essential to improve the efficiency of existing devices and pave the way for novel energy-efficient marine thrusters. In the present work, we investigate the effect of chord-wise flexibility on the propulsive performance of flapping-foil thrusters. For this purpose, a numerical method has been developed to simulate the time-dependent structural response of the flexible foil that undergoes prescribed large general motions. The fluid flow model is based on potential theory, whereas the elastic response of the foil is approximated by means of the classical Kirchhoff–Love theory for thin plates under cylindrical bending. The fully coupled FSI problem is treated numerically with a non-linear BEM–FEM scheme. The validity of the proposed scheme is established through comparisons against existing works. The performance of the flapping-foil thrusters over a range of design parameters, including flexural rigidity, Strouhal number, heaving and pitching amplitudes is also studied. The results show a propulsive efficiency enhancement of up to 6% for such systems with moderate loss in thrust, compared to rigid foils. Finally, the present model after enhancement could serve as a useful tool in the design, assessment and control of flexible biomimetic flapping-foil thrusters.


Introduction
The remarkable propulsive and manoeuvring mechanisms of aquatic swimmers, that have fascinated researchers since the 1970s, inspire the design of modern autonomous underwater vehicles (AUV) as well as autonomous underwater gliders (AUG) for marine environmental data acquisition, see, e.g., [1], biomimetic swimming robots and novel propulsion devices with enhanced efficiency, see, e.g., [2]. Selection of the swimming mode to serve as inspiration for the artificial devices closely depends on hydromechanical aspects of the application itself. The thunniform swimming mode for example, where the caudal fin of the fish performs a combination of pitching and heaving motions, identifies as the most efficient and therefore suitable for nature-inspired propulsion systems operating at high cruising speeds, see, e.g., [3,4]. Devices based on flapping foils have also been studied as auxiliary thrusters augmenting the overall ship propulsion in waves, see, e.g., the works [5][6][7][8][9] and project BIO-PROPSHIP [10]. Moreover, oscillating foils are studied for the development of hydrokinetic energy devices, see, e.g., [11][12][13] or hybrid devices with enhanced performance exploiting combined wave and tidal energy resources, see [14].
Living organisms through natural selection have been able to further enhance their locomotion capabilities through passive or active deformations of fins. In this direction, computational and experimental work on the principal mechanisms for thrust production in flexible oscillating bodies

Problem Description
In the present work we consider the unsteady motion of a large-aspect-ratio rectangular foil with chord length c and thickness profile τ(x), see Figure 1. In general, a foil made of flexible material bends and twists in all directions. However, for large-aspect-ratio foils under the assumption of cylindrical bending, spanwise deformations are neglected. In this work, our aim is to predict the inertia and fluid-driven chord-wise deformations of a foil that undergoes large general motions. The following Cartesian coordinate systems are introduced in this work: • the space-fixed frame (x, z), with respect to which the foil moves in the negative direction of x − axis with constant cruising speed U • the body-fixed (non-inertial) (x , z ) positioned at the foil's center of rotation with x − axis in the direction of the un-deformed chord line • the body-fixed (non-inertial) (x , z ) position at the leading edge (LE). This frame is exclusively used for the structural response problem.

Mathematical Formulation
To investigate the coupled FSI problem, we initially consider the hydrodynamic and the structural response problems independently. The former focuses on the transformation of the fluid flow around the foil and the other on the determination of the structural response of the body under excitation. Coupling is achieved through the fluid-induced structural deformations and subsequent nonlinear variation of the body boundary condition governing the hydrodynamics.

Structural Dynamics of the Foil
The foil is represented by a perfectly elastic, homogeneous and isotropic thin elastic plate. The dynamic structural response of the plate under cylindrical bending is modelled using the classical plate theory (CPT) based on the Kirchhoff-Love hypothesis [39]. For the formulation, we consider Additionally, the foil is subjected to a combination of harmonic heaving h(t) and pitching θ(t) motions, where h o , θ denotes the motion amplitudes, ψ the phase difference and f the oscillation frequency. In that sense, the effective angle of attack a e f f is For the hydrodynamic performance of flexible oscillating foils with large aspect ratio, the following modelling parameters identify as primary: (a) the non-dimensional heaving amplitude h o /c; (b) the feathering parameter (ratio of pitching angle θ(t) compared to the maximum angle of attack induced by the heave motion α e f f (t)); (c) the phase angle between heave and pitch ψ; (d) the relative position of the pitching axis x R (center or rotation) and (e) Strouhal number as a measure of unsteadiness St = f A/U, where f is the flapping frequency, A = 2h o is the nominal trailing edge amplitude, see, e.g., [15]. In this study, the characteristic flexural rigidity E/ρ s gc is also introduced, with ρ s denoting the material density and g the acceleration of gravity.

Mathematical Formulation
To investigate the coupled FSI problem, we initially consider the hydrodynamic and the structural response problems independently. The former focuses on the transformation of the fluid flow around the foil and the other on the determination of the structural response of the body under excitation. Coupling is achieved through the fluid-induced structural deformations and subsequent nonlinear variation of the body boundary condition governing the hydrodynamics.

Structural Dynamics of the Foil
The foil is represented by a perfectly elastic, homogeneous and isotropic thin elastic plate. The dynamic structural response of the plate under cylindrical bending is modelled using the classical plate theory (CPT) based on the Kirchhoff-Love hypothesis [39]. For the formulation, we consider the body-fixed coordinate system, in Figure 1b, positioned at the leading edge (LE), such that the x −plane coincides with the geometric mid-plane of the plate and z −axis is pointing upwards. The domain in this formulation is Ω = [xLE, xTE], and the plate's fabrication is assumed to be symmetric about the mid-surface. The governing equation of the initial boundary value problem (IBVP) with respect to the transverse displacement on the mid-plane is as follows m(x) · ∂ tt w(x; t) + ∂ xx (D(x)∂ xx w(x; t)) = q(x; t), x ∈ Ω, t > 0, (2) q(x; t) = 0.5ρ f U 2 δC p − m(x) x ..
where m(x) = ρ s τ(x) denotes the mass distribution, ρ f the fluid density, D(x) = Eτ 3 /12(1 − v 2 ) the flexural rigidity, v Poisson's ratio and E Young's modulus. The first term in Equation (3) consists of the fluid-driven forces and the second of the inertia driven ones. For the fluid-driven forces, δC p denotes the non-dimensional pressure difference between the upper and the lower sides of the foil supplemented by the unsteady hydrodynamic problem, see Section 3.3. The inertia-driven (or fictitious) forces are included in the modelling due to the non-inertial motions enforced at the body-fixed reference; see, e.g., [27]. The thickness profile of the foil has finite values at the leading (LE) and trailing edge (TE), which is assumed to be 0.02% of the foil's maximum thickness. Regarding the boundary conditions, the foil is assumed to be clamped at x R , with the leading and trailing edges remaining free from loading. The center of rotation is assumed to be fixed with zero deflection and slope.
Additionally, at the free edges, conditions of vanishing moment and shear force are applied, as follows supplemented by the following initial conditions, The equivalent weak formulation of the IBVP can be derived by multiplying Equation (2) by the test functions ϕ(x; t) ∈ H 2 (Ω) and performing integration by parts using the appropriate boundary conditions in Equations (4) and (5), see, e.g., [40]. The variational problem is formulated as follows. Find w so that ∀ϕ ∈ H 2 (Ω) it holds, and

Modelling the Fluid Flow around the Foil
The mathematical formulation of the hydrodynamic problem is based on the theory of incompressible, inviscid, potential flow under the assumption that the rotational part of the fluid flow is contained in the trailing vortex sheet. The flow region D ⊆ IR 2 is an open domain with time-dependent boundaries assumed to be smooth everywhere except the TE ∂D(t) = ∂D B (t) ∪ ∂D W (t). The first component ∂D B (t) refers to the foil's deformable exterior surface and the latter ∂D W (t) to the trailing vortex sheet with respect to the earth-fixed reference frame, see Figure 2. The body-fixed Cartesian coordinate system denoted by (x', z') fixed at the foil's centre of rotation ro, along chord length with no inclination, undergoes large general motions. In the present work the flexible large-aspect-ratio foil is fully submerged into the surrounding fluid, while its fabrication is symmetric to the camber line. It is important to note that the camber line is free to deform under the fluid-driven loads.
supplemented by the following initial conditions, The equivalent weak formulation of the IBVP can be derived by multiplying Equation (2) by the test functions and performing integration by parts using the appropriate boundary conditions in Equations (4) and (5), see, e.g., [40]. The variational problem is formulated as follows.
Find w so that and

Modelling the Fluid Flow around the Foil
The mathematical formulation of the hydrodynamic problem is based on the theory of incompressible, inviscid, potential flow under the assumption that the rotational part of the fluid flow is contained in the trailing vortex sheet. The flow region 2

D IR
⊆ is an open domain with time-dependent boundaries assumed to be smooth everywhere except the TE refers to the foil's deformable exterior surface and the latter ( ) W D t ∂ to the trailing vortex sheet with respect to the earth-fixed reference frame, see Figure 2. The body-fixed Cartesian coordinate system denoted by (x', z') fixed at the foil's centre of rotation O r , along chord length with no inclination, undergoes large general motions. In the present work the flexible large-aspect-ratio foil is fully submerged into the surrounding fluid, while its fabrication is symmetric to the camber line. It is important to note that the camber line is free to deform under the fluid-driven loads. The governing equation for the potential field is, supplemented by the no-entrance boundary condition, x n denotes the normal derivative, with n the unit normal vector on the body and B V the instantaneous velocity of the body due to oscillatory motions and elastic displacements. We treat the above as an initial value problem, while it is assumed that the The governing equation for the potential field is, supplemented by the no-entrance boundary condition, where ∂ n Φ(x; t) = ∇Φ(x; t) · n denotes the normal derivative, with n the unit normal vector on the body and V B the instantaneous velocity of the body due to oscillatory motions and elastic displacements. We treat the above as an initial value problem, while it is assumed that the disturbance potential and velocity vanish at a large distance from the body. On the trailing vortex sheet, the following kinematic and dynamic conditions must hold, with superscripts {+, −} denoting the upper and lower side of the wake, respectively, stating that the pressure p W and normal velocity ∂ n Φ W are continuous through the wake ∂D W . Using Equations (10) and (11) in conjunction with Bernoulli's theorem we obtain W is the potential jump on the wake and D/Dt = ∂/∂t + V m · ∇ the material derivative based on the mean velocity V m = 0.5(∇Φ + + ∇Φ − ) on the trailing vortex sheet. Under this approach, ∂D W evolves in time as a material curve, whose motion is part of the solution introducing an implicit non-linearity. In the present study, a time-stepping method (TSM), namely the free wake method, is employed for the trailing vortex sheet modelling. The generated vortex curve emanates parallel to the bisector of the TE, and the hydrodynamics of the freely moving-trailing vortex sheet is based on [41], where the position of the vortices evolves in time and using the unsteady wake rollup mollifier filtering technique [42]. In Figure 3, we present a comparison between the time-evolution of the free wake trailing vortex sheet and the simplified wake model, see [43]. The latter assumes that the vortices remain were shed. This linearization provides satisfactory predictions on integrated quantities such as the thrust, lift and moment coefficients in cases of low to moderate unsteadiness; see, e.g., [42,44].
we obtain is the potential jump on the wake and Under this approach, evolves in time as a material curve, whose motion is part of the solution introducing an implicit non-linearity. In the present study, a time-stepping method (TSM), namely the free wake method, is employed for the trailing vortex sheet modelling. The generated vortex curve emanates parallel to the bisector of the TE, and the hydrodynamics of the freely movingtrailing vortex sheet is based on [41], where the position of the vortices evolves in time and using the unsteady wake rollup mollifier filtering technique [42]. In Figure 3, we present a comparison between the time-evolution of the free wake trailing vortex sheet and the simplified wake model, see [43]. The latter assumes that the vortices remain were shed. This linearization provides satisfactory predictions on integrated quantities such as the thrust, lift and moment coefficients in cases of low to moderate unsteadiness; see, e.g., [42,44]. The study of lifting flows around hydrofoils in the context of potential theory requires another condition to be enforced on the trailing edge. In the present work we implement a nonlinear pressure-type Kutta condition. This condition requires the pressure difference at the trailing edge (TE) to be zero, see, e.g., [44], The study of lifting flows around hydrofoils in the context of potential theory requires another condition to be enforced on the trailing edge. In the present work we implement a nonlinear pressure-type Kutta condition. This condition requires the pressure difference at the trailing edge (TE) to be zero, see, e.g., [44], Applying the representation theorem to Equations (8)- (11), for every point x 0 ∈ ∂D B the following boundary integral equation (BIE) is obtained, where b(x; t) = V B · n B , and G(x 0 x) denotes the fundamental solution of the Laplace equation Next, ∂ n G(x 0 x) denotes the directional derivative and µ W the potential jump or dipole intensity on the wake, i.e., a quantity that changes over time, thus representing the history of circulation.

Hydrodynamic Pressure and Force
From Equation (12) we can derive the non-dimensional pressure coefficient along the body boundary, where p ∞ stands for the ambient pressure at infinity. The forces and moments excited on the foil are given below in the form of non-dimensional coefficients for the instantaneous lift, thrust and moment, respectively where r(s|s * ; t ) denotes the reference vector for moment calculation. In addition, the instantaneous power input coefficient is defined as The Froude efficiency is calculated as follows, where the bar denotes the mean value in time

Discretization Scheme for FEM
For the numerical solution of the variational problem defined by Equation (7), the domain of interest is discretised, while the unknown response is approximated by 5th order Hermite polynomials.
The employed Hermite element features three nodes and six degrees of freedom. Hence, approximate solutions are taken as, where w h i (t) denote the time-dependent nodal unknowns and H i (x) are the Hermite shape functions. A second order system of ODEs is derived when the approximate solution in Equation (25) is substituted into a discretized Equation (7) and the resulting formula is tested with the shape functions. Finally, the discretized system is written in matrix form as, where M glob , K glob are the global mass and stiffness matrices of dimension N T , respectively, F glob is the global load vector and U is the vector containing the nodal unknowns for the partitioned domain Ω h . Finally, N T refers also to the total degrees of freedom (DOFS). The global load vector in our study of chord-wise flexible foil introduces an implicit non-linearity to the problem through the fluid-driven term, see Section 3.1. The integrals appearing in the coefficients of Equation (26) are calculated by Gaussian quadrature. Details concerning the numerical implementation of the FEM scheme can be found in [39]. The addition of the proportional damping terms yields an extended global equation where α 1 , α 2 denote the proportional (or Rayleigh) damping coefficients. In the present work, these coefficients are approximated using the procedure described in [45] and adjusted based on comparisons with previous experimental work; see Section 5.3.

Time Integration
We proceed with implementing order reduction in Equation (27) and deriving the following system of non-linear first order differential equations, where F glob 0 and the identity matrix denoted as I. For numerical time integration of Equation (29) we use the Crank-Nicolson time integration scheme,

Boundary Integral Equation (BIE) & Discretization
Following a low-order panel method, see, e.g., [46], the boundary ∂D is decomposed into piecewise linear boundary elements. Concerning the representation of the potential, its normal derivative and the potential jump at each time step are assumed to be represented by piecewise constant distributions, as follows, Particularly, the boundary element on the wake that is closest to the TE is denoted as Kutta element. Finally, following a collocation scheme, the BIE in Equation (16) is satisfied in a finite number of points (collocation points), and in order to avoid singularities, the centroids of the elements have been chosen as collocation points. The discretized BIE is as follows, where In the above equations, δ ij is the Kronecker delta, In the following sections we will denote with bold the quantities containing the values of piecewise constant hydrodynamic functions on the panels at various parts of the boundary. For the induction factors it holds In the case of straight-line panels, the integrals in Equation (33) are calculated analytically, see, e.g., Kress [47], Katz-Plotkin [48]. Multiplying Equation (32) with A −1 we obtain Equation (34) denotes the Dirichlet-to-Neumann (DtN) operator that sets a mapping between the boundary values of the potential and its normal derivative. In the present work we propose two approaches for the solution to the hydrodynamic problem based on the use of BIE presented in Equations (34) and (35).

Solution Schemes
The unsteady hydrodynamics problem is the core of the coupled BEM-FEM numerical scheme, as presented in Section 4.3 that follows. The numerical treatment of the coupled FSI problem is quite computationally demanding, and for that reason, we propose the use of the two solution approaches developed for this problem, based on Strouhal number. For low-frequency simulations with Str < 0.25, we employ the least computationally demanding methodology, a BEM based on the Adams-Bashford-Moulton scheme (BEM-ABM). For problems of high unsteadiness, a more numerically stable and accurate methodology is used, a BEM based on a Newton-Raphson scheme (BEM-NR). Details regarding the behavior of these methodologies are presented below in Sections 5.1 and 5.3.
(i) BEM-ABM: The DtN operator, derived from the BIE, acts as a constraint to the dynamical system evolution equation that is constructed using the pressure-type Kutta condition. We consider µ W1 as the dynamic variable of the problem, and thus, the formulation allows for the treatment of an initial value problem (IVP). In order to express the pressure-type Kutta condition as a function of υ = µ W1 , we use the DtN map in Equation (34), in conjunction with the discretized form of Equation (14), to obtain a (spatially and temporarily) nonlocal differential equation, with explicit and implicit nonlinearities with respect to υ = µ W1 . The latter is finally put in the following form, where L(µ W1 ), N(µ W1 ) are linear and non-linear terms coming from the discretized version of the pressure-type Kutta condition, are scalars. Note that Ψ symbolizes the difference of a function Ψ at the trailing edge. For the numerical solution of the IVP Equations (36) and (37), we implement a higher-order Adams-Bashford-Moulton (ABM) scheme that provides accuracy, stability and efficiency. The following scheme requires the calculation of only two derivative quantities at each time step and has error of order (∆t 5 ), where ∆t is the time step. More details about the evolution of Equations (36) and (37) and the use of the Adams-Bashford numerical scheme can be found in [49].
(ii) BEM-NR: The BIE along with the discretized form of the pressure-type Kutta condition, detailed in Appendix A, is used to construct the complete system of equations, with the boundary fields Φ B and µ W1 as unknowns. A set of N B + 1 equations can be solved for the unknown values of Φ Bi and µ w1 at each time step, which can be written in a compact form For the present problem, a Newton-Raphson (NR) method is implemented at each time step as where J(x n ) −1 denotes the inverse of the system's Jacobian, which can be analytically calculated for the present formulation, see Appendix B. Finally, the BIE can be used for the calculation of Φ in the domain.

Non-Linear BEM-FEM
The fluid dynamic equations and the structural response are treated numerically with the same temporal discretization. The BEM approximates the solution to the unsteady hydrodynamic problem at each time step and provides predictions for the pressure and velocity fields. The pressure difference between the upper and lower sides of the foil acts as the fluid-driven load for the structural problem treated by FEM. By utilizing the proposed iterative scheme (coupling), the BEM solver also receives data from the FEM to re-evaluate (a) the foil-body geometry and (b) the velocities for the no-entry boundary condition at each iteration loop to finalize the solution approximation at each time step of the simulation.
The FSI is implicitly non-linear, which is inherent to the present passive deformation problem. The forcing term in the right-hand side of Equation (2) is dependent on the mid-plane deformation, which coincides with the foil's camber line and vice versa. Thus, the pressure forcing term should be formally written as q(x, t; w). In that sense, the coupling mechanism between the two solvers is introduced through the no-entrance boundary condition, see Equation (9). The instantaneous velocity for the latter boundary condition is the following where V rigid is the velocity due to the oscillatory motions of a rigid foil, and the second term is deformation velocity (as calculated in the body-fixed reference frame with the FEM solver projected to the previously un-deformed camber line). The unit normal vector n rigid = [− sin θ(t), cos θ(t)] T depends only on the pitching motion, see, e.g., [27] for a similar approach. For the structural response problem, Equation (30) is written in a more compact form as, For the solution of Equation (41), we employ a Newton-Raphson iterative scheme. Having determined an initial guess Q n+1 0 , the unknown vector is recursively approximated using where J is the Jacobian of the function G, the n-index refers to the time marching and q to the Newton-Raphson iterations. The initial guess can be obtained under various assumptions. In the present version of the computational code, the foil is assumed to be un-deformed at the beginning of the simulations, and the calculation of Q n 0 is obtained with the (a) geometry and (b) velocity field data from the previous time step. The calculation of the Jacobian matrix requires knowledge of the partial derivatives of the scalar components G i (Q), i = 1, . . . , 2N T of the function G(Q), where N T denotes the total degrees of freedom (DOFs) for the FEM and depends on chosen shape functions and the boundary conditions (BC). For example, for discretization with N elem = 5, and the BC presented in Section 3.1, the total DOFs are N T = 20. For the numerical approximation of the partial derivatives we implement a central difference scheme, where ε j is sufficiently small and in practice it is selected as a small percentage of Q j . Therefore, by

Results
Regarding the BEM solver, as presented in Section 4.2 for the hydrodynamic analysis of rigid flapping foils, extensive validation against experimental data as well as calculations concerning the solver's numerical performance over a range of motion parameters can be found in [49] and [37]. An indicative comparison of the 2D-BEM and experimental data in the case of a rigid flapping foil is shown below in below in Section 5.2, used as a reference for illustration of the effect of elastic deformation.
As concerns the accuracy of the present FEM solver, as a first example results concerning the free vibration analysis of tapered cantilever beams with taper ratio a, are considered. The thickness distribution along the beam is linear. The relative error for the first five eigenfrequencies, between the present FEM for a mesh of N elem = 15, and the analytical solution from [50], is listed in Table 1 for two values of the taper ratio. The present method results are in excellent agreement with reference values and are further enhanced with refined discretization. Another example concerns the static behavior of cantilever beams of length L = 10 m of variable thickness under tip load forcing F = −100 kN; see [51]. In Figure 4, we present a comparison between the FEM and data obtained from [51] regarding the transverse displacement. Finally, the FEM solver is validated against a dynamic test case of a cantilever beam of a constant thickness profile, under transverse dynamic tip loading F. The beam's response in terms of tip transverse displacement profile is compared in Figure 5 against the analytic solution presented in [52].

Convergence Characteristics of the Numerical Scheme
In this section, results concerning the convergence characteristics of the proposed numerical schemes are presented. To begin with the hydrodynamic problem, we present in Figures 6 and 7 the relative error concerning the calculated thrust coefficient and the efficiency of a rigid flapping foil against the time discretization parameter

Convergence Characteristics of the Numerical Scheme
In this section, results concerning the convergence characteristics of the proposed numerical schemes are presented. To begin with the hydrodynamic problem, we present in Figures 6 and 7 the relative error concerning the calculated thrust coefficient and the efficiency of a rigid flapping foil against the time discretization parameter ∆t/T and the number of elements N B subdividing the foil contour. The former represents the ratio of the time stepping over the motion's period and the latter the characteristic panel length on the hydrofoil, whereas the contours correspond to constant values of the ratio λ = U∆t/∆x. The results were obtained using both solver options, as presented in Section 4.2.2, including free wake modelling. Both numerical schemes converge for rigid flapping-foil simulations, i.e., the relative error is close to zero for the finest space and time discretization. In Figure 6, we can observe that a coarse discretization in time corresponds to greater values of relative error (up to 4%). For the propulsive efficiency, the relative error is significantly lower for discretization domain (up to 0.3%). It is important to note that the Adams-Bashford method is not as stable, since simulations that correspond to a coarse discretization in time and a finer mesh in space lead to numerical instabilities, which explains the non-symmetric mesh-grid. Particularly, for λ = U∆t/∆x = 3.5,∆t/T = 0.35%, the error is of the order of 0.5%, and thus the latter values are selected for the simulations.
A convergence study for the case of a chord-wise flexible flapping foil is presented in Figure 8. These simulations are obtained with the proposed coupled BEM-FEM scheme based on the Newton-Raphson iterative scheme for the fluid-flow problem. The region that corresponds to greater values of the relative error appears for simulation with a coarse discretization in time. The same behaviour is also observed for the propulsive efficiency. The iso-λ curves for the coupled BEM-FEM are not correlated to relative error minimization regions; therefore, we introduce / 0.45 t ∆ Τ < as a parameter constraint that is used for the numerical results that follow. In the sequel, numerical results obtained by the present model are compared against experimental measurements and data from other methods for validation.   On the other hand, the BEM based on the Newton-Raphson iterative scheme shows a difference behaviour. In Figure 7, the relative error has its maximum value of 1% and coincides with the region of coarse mesh both in space and time. This behaviour is a significant advantage that this numerical scheme shows quantitively along with the fact that it is more stable. Similarly, for λ = U∆t/∆x = 2.8,∆t/T = 0.85%, the error is of the order of 0.25%, and thus the latter values are selected for simulations.
A convergence study for the case of a chord-wise flexible flapping foil is presented in Figure 8. These simulations are obtained with the proposed coupled BEM-FEM scheme based on the Newton-Raphson iterative scheme for the fluid-flow problem. The region that corresponds to greater values of the relative error appears for simulation with a coarse discretization in time. The same behaviour is also observed for the propulsive efficiency. The iso-λ curves for the coupled BEM-FEM are not correlated to relative error minimization regions; therefore, we introduce ∆t/T < 0.45 as a parameter constraint that is used for the numerical results that follow. In the sequel, numerical results obtained by the present model are compared against experimental measurements and data from other methods for validation.

The Case of a Flexible Flapping Foil
The effects of chord-wise flexibility on the propulsive efficiency of two-dimensional flapping foils have been investigated experimentally in [32], showing that when properly selected, it leads to significant increase in the efficiency with small loss on thrust compared to the rigid foil. Simulations were performed for a NACA 0014 flapping foil with x R = 1/3c and the following kinematic parameters For the flexible foil, the thickness distribution along chord length coincides with the hydrodynamic shape of the foil, whereas the material properties correspond to relatively hard rubbers with Young's modulus E = 3.4617(·10 7 ) Pa, Poisson's ratio v = 0.4 and material density ρ s = 1100 kg/m 3 .
The time history of flapping motion (pitch and heave), lift and thrust forces for both the rigid and the flexible foil are shown in Figures 9 and 10, respectively. These results were obtained using space and time discretization N B = 160, N elem = 5, TSR = 0.4. For the case of the rigid foil in Figure 9, the present BEM provides a very accurate prediction of the hydrodynamic forces as was expected. The experimental data in Figure 10 are also in good agreement with the present numerical predictions in terms of the time averaged values and the overall periodic behavior of the instantaneous lift and thrust. In both cases, the proposed method predicts the maximum thrust with 1% accuracy. The differences in the peak thrust values in Figure 10 are attributed to the non-linear structural behavior and/or viscous effects, which are not modelled. It is interesting to note here that, in this case, the maximum displacement is small and occurs at the TE of the foil. Thus, incorporating effects of chord-wise flexibility, leads to significant enhancement in propulsive efficiency with a slight decrease of the thrust coefficient. The effects of chord-wise flexibility on the propulsive efficiency of two-dimensional flapping foils have been investigated experimentally in [32], showing that when properly selected, it leads to significant increase in the efficiency with small loss on thrust compared to the rigid foil.  The time history of flapping motion (pitch and heave), lift and thrust forces for both the rigid and the flexible foil are shown in Figures 9 and 10, respectively. These results were obtained using space and time discretization 160, For the case of the rigid foil in Figure 9, the present BEM provides a very accurate prediction of the hydrodynamic forces as was expected.
The experimental data in Figure 10 are also in good agreement with the present numerical predictions in terms of the time averaged values and the overall periodic behavior of the instantaneous lift and thrust. In both cases, the proposed method predicts the maximum thrust with 1% accuracy. The differences in the peak thrust values in Figure 10 are attributed to the non-linear structural behavior and/or viscous effects, which are not modelled. It is interesting to note here that, in this case, the maximum displacement is small and occurs at the TE of the foil. Thus, incorporating effects of chord-wise flexibility, leads to significant enhancement in propulsive efficiency with a slight decrease of the thrust coefficient.

The Case of a Flexible Heaving Foil
In order to further validate the present method and illustrate its ability to capture the main hydro-elastic effects of chord-wise flexible foils, we performed another series of simulations based on the experimental work of [34], for which case semi-analytical predictions are also available in [35]. In the latter work, the response of flexible plates performing heaving-only motions across a range of frequencies and heaving amplitudes was experimentally studied.
In this case, a flat plate with material properties For a plate immersed in fluid, see [34], it is reported that the first resonance frequency is equal to 4.71 rad/s ο ω = , whereas for the same structure in vacuo the corresponding frequency is estimated to be 14.96 rad/s. A comparison between the present method and the experimental data from [34] is shown in . The first coefficient (mass) affects the plate's response near the first resonance frequency, while the second coefficient (stiffness) has a more significant role in the higher frequency regime.
It is observed that the present method, based on either the BEM-NR or BEM-ABM, displays general agreement with the experimental results, especially around the first resonant frequency. Interestingly, the BEM-ABM provides very accurate prediction up to some values. A second resonance frequency is also evident in the range of frequency examined, where the elastic response is considerably smaller. In this case, the behavior of the BEM-ABM becomes less efficient. This is due to the enforcement of the pressure-type Kutta condition concerning the pressure difference at the

The Case of a Flexible Heaving Foil
In order to further validate the present method and illustrate its ability to capture the main hydro-elastic effects of chord-wise flexible foils, we performed another series of simulations based on the experimental work of [34], for which case semi-analytical predictions are also available in [35]. In the latter work, the response of flexible plates performing heaving-only motions across a range of frequencies and heaving amplitudes was experimentally studied.
In this case, a flat plate with material properties D = 0.018 Nm, ρ s = 1200 kg/m 3 is actuated at the LE with heaving amplitude h o /c = 0.033, oscillating frequency within the interval ω/ω = [0. 3,8] and Re = 6000. For a plate immersed in fluid, see [34], it is reported that the first resonance frequency is equal to ω = 4.71 rad/s, whereas for the same structure in vacuo the corresponding frequency is estimated to be 14.96 rad/s. A comparison between the present method and the experimental data from [34] is shown in Figures 11-14. The trailing edge/leading edge (TE/LE) amplitude response (A TE /A LE ) as a function of the non-dimensional frequency ω/ω is presented in Figure 11a. For comparison purposes we performed simulations with both solution approaches for the hydrodynamics problem. Particularly, for the simulations performed with BEM-NR, the following space and time discretization are used: N B = 140, N elem = 5, TSR = 0.35. Moreover, the damping terms (Section 4.1.1) were tuned to a = 2.5, b = 0.03. The first coefficient (mass) affects the plate's response near the first resonance frequency, while the second coefficient (stiffness) has a more significant role in the higher frequency regime. very good agreement, justifying the hybrid time-integration numerical scheme presented in Section 4.2.2.
Finally, in Figure 14, the deflection plots for time instances in a period of the heaving motion are plotted. Results are presented for two values of the non-dimensional frequency. We observe that, for / 3 ο ω ω = , the response of the plate displays a neck at around 2/3 of the chord due to the second plate mode excitation. These results agree well with the experimental data presented in [34,35].

Effects of Flexural Rigidity on Froude Efficiency
Carefully chosen flexibility characteristics have the potential to further enhance the propulsive performance of flapping-foil thrust, as has been reported in [32] and confirmed from the simulations presented in Section 5.  Figure 15, where it is observed that as Young's modulus is reduced, the propulsion efficiency rises. Indeed, an efficiency increase of 6% is observed for , as compared to the rigid case. This, however, is at the cost of thrust reduction.
Especially, in cases when the kinematic parameters are not optimized, i.e., purely heaving motion, It is observed that the present method, based on either the BEM-NR or BEM-ABM, displays general agreement with the experimental results, especially around the first resonant frequency. Interestingly, the BEM-ABM provides very accurate prediction up to some values. A second resonance frequency is also evident in the range of frequency examined, where the elastic response is considerably smaller. In this case, the behavior of the BEM-ABM becomes less efficient. This is due to the enforcement of the pressure-type Kutta condition concerning the pressure difference at the trailing edge. The solution scheme based on the BEM-NR satisfies exactly the pressure-type Kutta condition for all frequencies; however, the BEM-ABM despite the fine discretization used (N B = 150, N elem = 5, TSR = 0.05) leads to a finite pressure difference at the TE, as shown in Figure 11b. This is further illustrated in the comparison between the distribution of pressure coefficient in Figure 2, for the case of BEM-ABM and BEM-NR. The pressure distributions are very similar, leading to compatible predictions of integrated forces and moments on the foil. However, the calculated pressure difference at the trailing edge by the BEM-ABM, affects the value of the shed vorticity and produces error in the numerical solution as the frequency increases further, see Figure 11b.
In this example, for ω/ω > 5, results have been obtained only by the BEM-NR model, and the second resonant peak is underestimated, a finding that agrees with similar predictions from the work by [35], which is attributed to the inability of potential-based methods to account for viscous effects manifesting at higher frequencies. In the present work, a proportional damping is employed, and the use of a more complex damping model, see, e.g., [27,35], that could further improve the results is left to be examined in future work.
A successful comparison between the numerical model and the experimental results regarding the TE/LE phase lag, as observed in the earth-fixed reference frame, is presented in Figure 13a. For higher values of the frequency, the error in the phase lag prediction increases. This is also the case for the TE/LE amplitude response ratio in Figure 11a. Nevertheless, the present method predictions are still within acceptable limits. Moreover, in Figure 13b it is shown that for thrust predictions the two solution approaches, namely the BEM-ABM and the BEM-NR, for the hydrodynamic problem are in very good agreement, justifying the hybrid time-integration numerical scheme presented in Section 4.2.2.
Finally, in Figure 14, the deflection plots for time instances in a period of the heaving motion are plotted. Results are presented for two values of the non-dimensional frequency. We observe that, for ω/ω = 3, the response of the plate displays a neck at around 2/3 of the chord due to the second plate mode excitation. These results agree well with the experimental data presented in [34,35].

Effects of Flexural Rigidity on Froude Efficiency
Carefully chosen flexibility characteristics have the potential to further enhance the propulsive performance of flapping-foil thrust, as has been reported in [32] and confirmed from the simulations presented in Section 5.  Figure 15, where it is observed that as Young's modulus is reduced, the propulsion efficiency rises. Indeed, an efficiency increase of 6% is observed for θ = 30deg, as compared to the rigid case. This, however, is at the cost of thrust reduction. Especially, in cases when the kinematic parameters are not optimized, i.e., purely heaving motion, carefully choosing flexibility has the potential to enhance the propulsion efficiency considerably. Motivated by the work of [25], we estimated the maximum effective angle of attack a m for the flexible foil based on Equation (1b) with For the most stiff foil with E ∼ 10 8 Pa, the estimated maximum effective angle of attacks, that correspond to θ = {0 • , 10 • , 20 • , 30 • } are a m = {43 • , 33 • , 23 • , 13 • }. Typically a decreasing a m leads to a decrease in thrust and an increase in efficiency for rigid foils. In that sense, the value a m = 40 • corresponding to zero pitching amplitude for the most flexible foil E ≈ 10 5 Pa explains the behavior of the results in Figure 15. A qualitative explanation is also provided in Figure 16 The chosen material corresponds to one of the most elastic ones examined before, characterized by

Conclusions
Flapping foils with chord-wise flexibility were studied in this work as unsteady thrusters with enhanced propulsive performance. To investigate the hydroelasticity effects on the thrust and propulsive efficiency of such systems, a mathematical model is proposed for the FSI problem. The fluid flow modelling is based on potential theory, whereas the elastic response of the foil is based on the Kirchhoff-Love theory for thin plates under cylindrical bending. A non-linear fully coupled BEM-FEM numerical scheme is developed to simulate the time-dependent structural response of the flexible foil undergoing large prescribed general motions. The proposed iterative scheme ensures stability and convergence of the coupled numerical simulation, as proven by the convergence study shown in Figures 6-8. The present method is also extensively compared against experimental data for validation, demonstrating its ability to capture the main aspects of the FSI problem.
The proposed method has been successfully compared with experimental data found in [32] for the case of a chord-wise flexible foil performing combined heaving and pitching motions; see Figures 9 and 10. The results indicate that incorporating chord-wise flexibility in flapping-foil design could lead to 13% enhancement of propulsive efficiency, as compared to rigid foils.
The response of flexible plates performing heaving-only motions across a range of forcing frequencies and heaving amplitudes, found in the experimental work of [34], has also been studied for comparison purposes. These numerical results were in good agreement with the experimental measurements for various aspects of the non-linear dynamic system; see  The present method is shown to satisfactorily predict both the TE/LE amplitude response as a function of the oscillating frequency, see Figure 11a, and the phase lag between the LE and TE as reported during the experiments; see Figure 13a. The first and second resonant frequencies were quite accurately predicted; however, our model slightly underestimates the TE amplitude response near the second resonance. Furthermore, the envelopes of the foil's elastic deflection, see Figure 14, agree with the predictions presented in the work of [35]. In that sense, the non-linear BEM-FEM scheme is shown to

Conclusions
Flapping foils with chord-wise flexibility were studied in this work as unsteady thrusters with enhanced propulsive performance. To investigate the hydroelasticity effects on the thrust and propulsive efficiency of such systems, a mathematical model is proposed for the FSI problem. The fluid flow modelling is based on potential theory, whereas the elastic response of the foil is based on the Kirchhoff-Love theory for thin plates under cylindrical bending. A non-linear fully coupled BEM-FEM numerical scheme is developed to simulate the time-dependent structural response of the flexible foil undergoing large prescribed general motions. The proposed iterative scheme ensures stability and convergence of the coupled numerical simulation, as proven by the convergence study shown in Figures 6-8. The present method is also extensively compared against experimental data for validation, demonstrating its ability to capture the main aspects of the FSI problem.
The proposed method has been successfully compared with experimental data found in [32] for the case of a chord-wise flexible foil performing combined heaving and pitching motions; see Figures 9 and 10. The results indicate that incorporating chord-wise flexibility in flapping-foil design could lead to 13% enhancement of propulsive efficiency, as compared to rigid foils.
The response of flexible plates performing heaving-only motions across a range of forcing frequencies and heaving amplitudes, found in the experimental work of [34], has also been studied for comparison purposes. These numerical results were in good agreement with the experimental measurements for various aspects of the non-linear dynamic system; see  The present method is shown to satisfactorily predict both the TE/LE amplitude response as a function of the oscillating frequency, see Figure 11a, and the phase lag between the LE and TE as reported during the experiments; see Figure 13a. The first and second resonant frequencies were quite accurately predicted; however, our model slightly underestimates the TE amplitude response near the second resonance. Furthermore, the envelopes of the foil's elastic deflection, see Figure 14, agree with the predictions presented in the work of [35]. In that sense, the non-linear BEM-FEM scheme is shown to successfully predict the hydrodynamic loads as well as the fluid-driven deformation of flexible flapping foils with general thickness profile and flexural rigidity.
Motivated by the propulsive performance enhancement offered by flexible foil-thrusters, we performed parametric studies, see Figures 15 and 17, in order to further investigate the effects of elasticity over a range of design parameters, including Strouhal number, heaving and pitching amplitudes. The results illustrate that chord-wise flexibility and flexural rigidity profile variations can significantly improve the propulsive efficiency of the biomimetic thruster. Particularly, it is shown in Figure 15 that as flexural rigidity is reduced, the propulsion efficiency rises, leading to an efficiency increase as large as 6% observed for θ = 30deg, as compared to the rigid case. This, however, is obtained at the cost of thrust reduction. The results in Figure 17 illustrate that chord-wise flexibility leads to a more pronounced decrease in the thrust as Strouhal number increases, especially for higher amplitudes of heaving motion. On the contrary, flexibility enhances the propulsive efficiency for high amplitudes of heaving motion and Strouhal numbers.
To conclude, future work is planned towards the detailed investigation and systematic examination of the structural response of the flexible foil over a range of design and operation parameters, including flexural rigidity profiles inspired by nature. Additional comparisons and benchmark studies between the present non-viscous BEM-FEM scheme and high-fidelity viscous CFD solvers is also left for future work. Direct extensions include modelling various nonlinearities associated with large deflections and viscous effects [37,53]. Another aspect concerns the code optimization using GPGPU programming and message passing interface (MPI) techniques to significantly reduce computational time and cost, see, e.g., [13,38,44]. This step will allow three-dimensional modelling as well as shape and material optimization, supporting applications concerning realistic designs. Finally, the present method could also find useful application to calculate the flexibility effects on the performance of novel marine renewable energy devices based on oscillating foils; see [13].

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

Appendix A
The discretized pressure-type Kutta condition is given by the following set of equations: where τ is the unit tangent vector on the body contour defined in the clockwise direction, and d j is the curvilinear distance between the midpoints of the ( j, j + 1) panels.

Appendix B
A finite difference method (FDM) is used for the temporal and spatial discretization of the pressure-type Kutta condition in order to form a system of nonlinear equations along with the BIE relation with respect to the unknown boundary fields Φ Bi and µ W1 at the vicinity of the trailing edge. The resulting system of equations can be solved numerically after the appropriate discretization at each time step of the simulation. Particularly, a backward finite difference scheme in time combined with forward and backward differences in space has been used for the discretization of the pressure-type Kutta condition in the set of Equations (A7)-(A11) as follows, The linearized form of the above equation is where c k = τ x,Bk i + τ y,Bk j k = 1, N B . In the relations above, τ in the above relations refers to the unit tangent vector on the body contour defined in the clockwise direction, and d j is the curvilinear distance between the midpoints of the ( j, j + 1) panels. In addition, Returning now to the discretized form of the boundary integral Equation (38), we derive the following expression by re-arranging terms, so that for (x i , y i ) , i = 1, . . . , N B : In this form, all the quantities in the rhs are known from the prescribed kinematics of the foil and the history of circulation of the foil that has been evaluated at previous time steps. Equations (A7a) and (A8)-(A12) form a set of N B + 1 equations, which can be solved for the unknown values of Φ Bi and µ w1 at each time step. Equations (A7b) and (A8)-(A12) consist of a linear system of equations that can be solved explicitly for the unknown values, that is, the initial guess for the solution of the nonlinear system of Equations (A7a) and (A8)-(A12) using a general iterative method.