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

Abstract

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.

1. 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 can be found in the review [15] with emphasis on aquatic locomotion as well as the works [16,17] regarding aerodynamics and aeroelasticity. Further investigation on the parameters that result in performance enhancement of flexible systems remains essential for the development of successful engineering applications in the future.

Over the years many researchers have presented mathematical models and numerical schemes to tackle this complex fluid–structure interaction (FSI) problem. One approach is to assume that the deformations are known a priori, see, e.g., [18]. The coupled problem can also be treated with strongly coupled methods; i.e., a fluid flow solver used in conjunction with a structural solver, employing iterative processes to solve the coupled problem of passive deformations, see, e.g., [19,20].

Notable is the approach by [21,22] to investigate thrust production in flexible wings with a potential-based model. Later, [23] presented a potential-based 2D flexible body vortex sheet model to estimate propulsive forces and optimal efficiency, which was further validated experimentally in [24]. Fully coupled FSI simulations for flapping wings with chord-wise and span-wise flexibility were carried out in the work of [25], with extensions in [26] allowing even the study of skeleton-strengthened fins. By employing the vortex-based method for the fluid flow problem along with a finite difference method for the structural response problem as presented in [27] it is shown that chord-wise flexibility enhances up to 10% the propulsive performance of a flapping-foil thruster for AUV’s.

The effects of flexibility in propulsive performance enhancement are also illustrated in experimental works, see, e.g., [28,29,30]. The effects of flexibility on the time-averaged thrust can be beneficial for plunging foils even when the flexible region is confined to a small section near the trailing edge, see, e.g., [31], compared to rigid foil cases. Moreover, proper selection of chord-wise flexibility characteristics in 2D foils in flapping regimes leads to up to 36% increase in efficiency compared to rigid ones, see, e.g., [32]. During experiments on elastic flat plates it has been observed that the more flexible plates showed enhanced propulsion characteristics, see, e.g., [33]. Other works [34,35] both theoretical and experimental, investigated the structural response of flexible flat plates under heaving motions. The experiments revealed that the thrust displayed peaks in motion frequency values coinciding with the resonance frequencies of the system comprised of the foil and the surrounding fluid. Following a different approach, see, e.g., [36], the propulsive performance of a flapping system that allows different inclinations of the robotic fin tip has been studied in terms of propulsive efficiency. Despite the increased efforts and contributions in this field, investigation on the parameters that result in performance enhancement of flexible systems remains essential for the development of successful engineering applications in the future.

Especially for fluid flow problems at high Reynolds numbers, moderate angles of attack and Strouhal number, which is the operating regime for marine thrusters, potential-based methods have proven to be suitable for the study of elasticity effects on thrust production. In that sense, lower-fidelity and cost-effective inviscid fluid flow simulations are a useful tool in the preliminary design phase of biomimetic thrusters, where emphasis is given on parametric studies, whereas high-fidelity CFD simulations, see, e.g., [37], are more computationally intensive and resource demanding compared to potential-based solvers. A direct comparison between the two approaches can be found in [38].

The scope of the present work is to propose a non-linear fully coupled BEM–FEM scheme for the solution of the FSI problem of chord-wise flexible flapping foils operating as marine thrusters, including thickness and flexural rigidity profile variation. The major contribution of this work is that the developed method and software could serve, after enhancement and further validation, as a useful cost-effective tool for the preliminary design and optimum control of flexible biomimetic thrusters. Although the proposed cost-effective method is less accurate than high-fidelity CFD methods at specific regions of the operational parameters (large angles of attack and Strouhal number), it is suitable for the investigation of the effects of flexibility on the propulsive performance of such systems, the accurate prediction of hydrodynamic loads as well as the structural response of the foil.

In the present formulation, the flexible foil performs prescribed oscillatory motions and is free to deform under inertial and reactive forces caused by its motion and hydrodynamic pressure. The hydro-elastic analysis is based on a non-linear time domain boundary element method (BEM) for the unsteady lifting flow problem and a high-order finite element method (FEM) for the prediction of the deformed body geometry using the Kirchhoff–Love thin plate theory for cylindrical bending under plane strain conditions. The foil is considered clamped at a position along chord length, while its trailing edge (TE) and leading edge (LE) act as free ends. The structural and hydrodynamic aspects of the problem are coupled in an implicit manner, encapsulating forms of non-linearity.

The structure of the paper is as follows: In Section 2, we describe the physical problem of a flexible flapping-foil thruster. Section 3 presents the mathematical formulation of the hydromechanical problem and Section 4 the numerical methods as well as the proposed BEM–FEM coupling schemes. In Section 5 we present numerical results consisting of a convergence study and comparisons of the method with experimental data. The effects of elasticity over a range of design parameters, including Strouhal number, heaving and pitching amplitudes are presented and discussed in a series of parametric numerical studies, illustrating that chord-wise flexibility and flexural rigidity profile variations can significantly improve the propulsive efficiency of the biomimetic thruster. Conclusions and a short discussion are presented in Section 6.

2. 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 $\tau (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^{\prime },z^{\prime })$ positioned at the foil’s center of rotation with $x^{\prime }-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.

Additionally, the foil is subjected to a combination of harmonic heaving $h(t)$ and pitching $\theta (t)$ motions,

$$h(t)=h_o\sin (2\pi f(t-t_o)), \theta (t)={\theta }_o\sin (2\pi f(t-t_o)+\psi ),$$

(1a)

where $h_o,{\theta }_o$ denotes the motion amplitudes, $\psi$ the phase difference and $f$ the oscillation frequency. In that sense, the effective angle of attack $a_{eff}$ is

$${\alpha }_{eff}(t)=\arctan [\dot{h}/U]-\theta (t).$$

(1b)

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 $\theta (t)$ compared to the maximum angle of attack induced by the heave motion ${\alpha }_{eff}(t)$); (c) the phase angle between heave and pitch $\psi$; (d) the relative position of the pitching axis $x_R$ (center or rotation) and (e) Strouhal number as a measure of unsteadiness $St=fA/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/{\rho }_{s}gc$ is also introduced, with ${\rho }_s$ denoting the material density and $g$ the acceleration of gravity.

3. 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.

3.1. 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 $\Omega =[x{}_{LE},x{}_{TE}]$, 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)\cdot {\partial }_{tt}w(x;t)+{\partial }_{xx}\left(D(x){\partial }_{xx}w(x;t)\right)=q(x;t), \mathrm{x}\in \Omega ,t>0,$$

(2)

$$q(x;t)=0.5{\mathsf{\rho }}_f{\mathrm{U}}_{}^2\delta {\mathrm{C}}_p-m(x)\left(x\ddot{\theta }(t)+\ddot{h}(t)\cos \theta (t)\right),$$

(3)

where $m(x)={\rho }_s\tau (x)$ denotes the mass distribution, ${\rho }_f$ the fluid density, $D(x)=E{\tau }^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, $\delta {\mathrm{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.

$${w(x;t)|}_{x=x_R}={{\partial }_{x}w(x;t)|}_{x=x_R}=0.$$

(4)

Additionally, at the free edges, conditions of vanishing moment and shear force are applied, as follows

$$\begin{array}{l}{D(x){\partial }_{xx}w|}_{x=x_{LE}}={{\partial }_x\left(D(x){\partial }_{xx}w\right)|}_{x=x_{LE}}=0\\ {D(x){\partial }_{xx}w|}_{x=x_{TE}}={{\partial }_x\left(D(x){\partial }_{xx}w\right)|}_{x=x_{TE}}=0.\end{array}$$

(5)

supplemented by the following initial conditions,

$${w(x;t)|}_{t=0}={{\partial }_{t}w(x;t)|}_{t=0}=0.$$

(6)

The equivalent weak formulation of the IBVP can be derived by multiplying Equation (2) by the test functions $\phi (x;t)\in {\mathsf{{\rm H}}}^2\left(\mathsf{\Omega }\right)$ 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 $\forall \phi \in {\mathsf{{\rm H}}}^2\left(\mathsf{\Omega }\right)$ it holds,

$${\displaystyle {\int }_0^L{\partial }_{xx}\phi D(x){\partial }_{xx}wdx}+{\displaystyle {\int }_0^L\phi m(x){\partial }_{tt}w\phi }dx={\displaystyle {\int }_0^L\left(q_{hydro}+q_{inertial}\right)\phi dx},$$

(7)

and ${\langle w(x,0),\phi \rangle |}_{L^2(\Omega )}={\langle {\partial }_{t}w(x,0),\phi \rangle |}_{L^2(\Omega )}=0,$ where ${\langle {\phi }_1,{\phi }_2\rangle |}_{L^2(\Omega )}={\displaystyle {\int }_0^L{\phi }_1,{\phi }_2}dx$ denotes the $L^2(\Omega )$ inner product.

3.2. 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\subseteq I{R}^2$ is an open domain with time-dependent boundaries assumed to be smooth everywhere except the TE $\partial D(t)=\partial D_B(t)\cup \partial D_W(t).$ The first component $\partial D_B(t)$ refers to the foil’s deformable exterior surface and the latter $\partial 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 $r{}_o$, 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,

$$\Delta \Phi (x;t)=0,x=(x,z)\in D,$$

(8)

supplemented by the no-entrance boundary condition,

$${\partial }_n\Phi (x;t)=V_B(x;t)\cdot n(x;t), x\in \partial {\mathrm{D}}_B,$$

(9)

where ${\partial }_n\Phi (x;t)=\nabla \Phi (x;t)\cdot 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,

$${\partial }_n{\Phi }_W^{+}(x;t)={\partial }_n{\Phi }_W^{-}(x;t), x\in \partial D_W,$$

(10)

$$p_W^{+}(x;t)=p_W^{-}(x;t) , x\in \partial D_W,$$

(11)

with superscripts $\{+,-\}$ denoting the upper and lower side of the wake, respectively, stating that the pressure $p_W$ and normal velocity ${\partial }_n{\Phi }_W$ are continuous through the wake $\partial D_W$. Using Equations (10) and (11) in conjunction with Bernoulli’s theorem

$$\frac{p(x;t)}{\rho }+{\partial }_t\Phi (x;t)+\frac{1}{2}{\left[\nabla \Phi (x;t)\right]}^2=0, x\in \partial D,$$

(12)

we obtain

$$\frac{D{\mu }_W(x;t)}{Dt}=0, x\in \partial D_W,$$

(13)

where ${\mu }_W=〚{\Phi }_W^{}〛={\Phi }_W^{+}-{\Phi }_W^{-}$ is the potential jump on the wake and $D/Dt=\partial /\partial t +V_m\cdot \nabla$ the material derivative based on the mean velocity $V_m=0.5\left(\nabla {\Phi }_{}^{+}+\nabla {\Phi }_{}^{-}\right)$ on the trailing vortex sheet. Under this approach, $\partial 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].

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],

$${\partial }_t\left({\Phi }_{}^{+}-{\Phi }_{}^{-}\right)+0.5\left(\nabla {\Phi }_{TE}^{+}-\nabla {\Phi }_{TE}^{-}\right)\left(\nabla {\Phi }_{TE}^{+}+\nabla {\Phi }_{TE}^{-}\right)=0.$$

(14)

Applying the representation theorem to Equations (8)–(11), for every point $x_0\in \partial D_B$ the following boundary integral equation (BIE) is obtained,

$$\begin{array}{l}\frac{1}{2}{\Phi }_B(x_0;t)+{\displaystyle \underset{\partial D_B(t)}{\int }{\Phi }_B(x;t){\partial }_{n}G(x_0|x)ds(x)}=\\ ={\displaystyle \underset{\partial D_B(t)}{\int }\left(V_B\cdot n_B\right)G(x_0|x)}ds(x)-{\displaystyle \underset{\partial D_W(t)}{\int }〚{\Phi }_W〛(x;t){\partial }_{n}G(x_0|x)}ds(x),\end{array}$$

(15)

$$\mathrm{with} 〚{\Phi }_W〛={\Phi }_W^{+}-{\Phi }_W^{-}={\mu }_W,$$

(16)

where $b\left(x;t\right)=V_B\cdot n_B$, and $G(x_0|x)$ denotes the fundamental solution of the Laplace equation

$$G(x_0|x)=\frac{1}{2\pi }\ln r(x_0|x) , r=\left|x_0-x\right| .$$

(17)

Next, ${\partial }_{n}G(x_0|x)$ denotes the directional derivative and ${\mu }_W$ the potential jump or dipole intensity on the wake, i.e., a quantity that changes over time, thus representing the history of circulation.

3.3. Hydrodynamic Pressure and Force

From Equation (12) we can derive the non-dimensional pressure coefficient along the body boundary,

$$C_P=(p-p_{\infty })/(0.5\rho U_{}^2),$$

(18)

where $p_{\infty }$ 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

$$C_L=L(t)/\left(0.5\rho U_{}^{2}c\right)=-\frac{1}{c}{\displaystyle \underset{\partial D_B}{\int }(C_p\cdot n)\cdot \widehat{y}ds},$$

(19)

$$C_T=T(t)/\left(0.5\rho U_{}^{2}c\right)=\frac{1}{c}{\displaystyle \underset{\partial D_B}{\int }(C_p\cdot n)\cdot \widehat{x}ds},$$

(20)

$$C_M=M(t)/\left(0.5\rho U_{}^2{c}^2\right)=-\frac{1}{c^2}{\displaystyle \underset{\partial D_B}{\int }(C_p\cdot n)\cdot r(s|s^{\ast };t)ds},$$

(21)

where $r(s|s^{\ast };t)$ denotes the reference vector for moment calculation. In addition, the instantaneous power input coefficient is defined as

$$C_{P_{in}}(t)=P_{in}(t)/\left(0.5\rho U_{}^{3}c\right),$$

(22)

$$P_{in}(t)=L(t)\dot{h}+M(t)\dot{\theta } .$$

(23)

The Froude efficiency is calculated as follows, where the bar denotes the mean value in time

$$\eta =U\overline{T}/{\overline{P}}_{in}={\overline{C}}_T/{\overline{C}}_{Pin}.$$

(24)

4. Numerical Methods

4.1. Finite Element Method (FEM)

4.1.1. 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,

$$w^h={\displaystyle \sum _i^6{w}_i^h(t)H_i(x)},$$

(25)

where $w_i^h(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,

$$M_{glob}\ddot{U}+K_{glob}U=F,$$

(26)

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 ${\Omega }^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

$$M_{glob}\ddot{U}+C\dot{U}+K_{glob}U=F,$$

(27)

$$C=a_1{M}_{glob}+a_2{K}_{glob},$$

(28)

where ${\alpha }_1,{\alpha }_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.

4.1.2. Time Integration

We proceed with implementing order reduction in Equation (27) and deriving the following system of non-linear first order differential equations,

$$A\dot{Q}=BQ+F,Q=[U{\dot{U}]}^T,$$

(29)

where $A=[\begin{array}{c}0\\ -I\end{array} \begin{array}{c}{M}_{glob}\\ 0\end{array}],$$B=-[\begin{array}{c}{K}_{glob}\\ 0\end{array} \begin{array}{c}C\\ I\end{array}],$$F=[\begin{array}{c}{F}_{glob}\\ 0\end{array} ]$ and the identity matrix denoted as $I$. For numerical time integration of Equation (29) we use the Crank–Nicolson time integration scheme,

$$\left(A-\frac{1}{2}\Delta t\cdot B\right)Q(t_{n+1})-\left(A+\frac{1}{2}\Delta t\cdot B\right)Q(t_n)-\frac{1}{2}\Delta t[F(t_{n+1})+F(t_n)]=0.$$

(30)

4.2. Boundary Element Method (BEM)

4.2.1. Boundary Integral Equation (BIE) & Discretization

Following a low-order panel method, see, e.g., [46], the boundary $\partial 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,

$$\begin{array}{l}{\Phi }_{{\rm B}}(x;t)={\Phi }_{{\rm B}i}(t),i=1\dots ,N_B,\\ {\partial }_n{\Phi }_{{\rm B}}(x;t)={\partial }_n{\Phi }_{{\rm B}i}(t),i=1\dots ,N_B,\\ {\mu }_W(x;t)={\mu }_{Wk}(t),k=1\dots ,N_W(t).\end{array}$$

(31)

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,

$$A\cdot {\Phi }_{{\rm B}}=S\cdot b+W\cdot {\mu }_W+W_K{\mu }_{W1},$$

(32)

where

$$\begin{gathered} A=\left\{\frac{{\delta }_{ij}}{2}+B_{ij}^{}\right\},S=\left\{A_{ij}^{}\right\}, \mathrm{i}\in \{1,\dots {\mathrm{N}}_B\}, \mathrm{j}\in \{1,\dots {\mathrm{N}}_B\}, \\ W=\left\{-B_{ik}^{}\right\}, W_K=\left\{-B_{i1}^{}\right\}, i\in \left\{1,\dots ,N_B\right\}, k\in \left\{2,\dots ,N_W(t)\right\}. \end{gathered}$$

In the above equations, ${\delta }_{ij}$ is the Kronecker delta, ${\Phi }_{{\rm B}}=\left\{{\Phi }_{Bi}\right\},$$b=\left\{{\partial }_n{\Phi }_{Bi}\right\},$${\mu }_W=\left\{{\mu }_{Wk}\right\}.$ 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

$$A_{ij}^{}={\displaystyle \underset{panel j}{\int }{G}_s^{}(x_i|x)dS}(x),B_{ij}^{}={\displaystyle \underset{panel j}{\int }{\partial }_n{G}_s^{}(x_i|x)dS}(x).$$

(33)

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

$${\Phi }_{{\rm B}}=D\cdot b+P\cdot {\mu }_W+Z{\mu }_{W1},$$

(34)

$$D=A^{-1}S, P({\mu }_W)=A^{-1}\left(W\cdot {\mu }_W\right), Z=A^{-1}{W}_K.$$

(35)

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).

4.2.2. 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 Section 5.1 and Section 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 ${\mu }_{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 $\upsilon ={\mu }_{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 $\upsilon ={\mu }_{W1}.$ The latter is finally put in the following form,

$$\frac{d\upsilon }{dt}=\mathrm{f}(\upsilon ), \mathsf{\upsilon }(t_o)={\upsilon }_o,$$

(36)

$$\mathrm{f}(\upsilon )=\frac{-1}{〚{\rm Z}〛}\frac{d\left(〚{\rm Z}〛\right)}{dt}\cdot \upsilon + \frac{-1}{〚{\rm Z}〛}\frac{d\left(〚D〛\cdot b+〚P〛\right)}{dt} + \frac{N\left(\upsilon \right)+L\left(\upsilon \right)}{〚{\rm Z}〛},$$

(37)

where $L\left({\mu }_{W1}\right),N\left({\mu }_{W1}\right)$ are linear and non-linear terms coming from the discretized version of the pressure-type Kutta condition, $〚D〛={〚{\left\{D\right\}}_{ij}〛}_{i=1}^{i=N_B}$ is a vector and $〚P〛={〚P_i〛}_{i=1}^{i=N_B}$, $〚Z〛={〚Z_i〛}_{i=1}^{i=N_B}$ are scalars. Note that $〚\Psi 〛$ symbolizes the difference of a function $\Psi$ 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 $(\Delta t^5)$, where $\Delta 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 ${\Phi }_B$ and ${\mu }_{W1}$ as unknowns. A set of $N_B+1$ equations can be solved for the unknown values of ${\Phi }_{Bi}$ and ${\mu }_{w1}$ at each time step, which can be written in a compact form

$$f\left(x\right)=0,x={[{\Phi }_{{\rm B}1} \dots {\Phi }_{{\rm B}{{\rm N}}_{{\rm B}}} {\mu }_{w1}]}^T.$$

(38)

For the present problem, a Newton–Raphson (NR) method is implemented at each time step as

$$x_{n+1}=x_n-J{(x_n)}^{-1}f(x_n),$$

(39)

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 $\Phi$ in the domain.

4.3. 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

$$\begin{gathered} V_B=V_{rigid}+{\partial }_{t}w\cdot n_{rigid} \\ \mathrm{with} n_{rigid}={[-\sin \theta (t), \cos \theta (t)]}^T, \end{gathered}$$

(40)

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 \theta (t), \cos \theta (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,

$$G\left(Q_{n+1}\right)=0.$$

(41)

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

$$Q_{n+1}^{q+1}=Q_{n+1}^q-J^{-1}\left(Q_{n+1}^q\right)\cdot G\left(Q_{n+1}^q\right), q=1,2,\dots ,$$

(42)

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\left(Q\right), i=1,\dots ,2N_T$ of the function $G\left(Q\right),$ 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,

$$\frac{\partial G_i}{\partial Q_j}\approx \frac{G_i\left(Q_j+{\epsilon }_j\right)-G_i\left(Q_j-{\epsilon }_j\right)}{2{\epsilon }_j},$$

(43)

where ${\epsilon }_j$ is sufficiently small and in practice it is selected as a small percentage of $\left|Q_j\right|$. Therefore, by defining a perturbation vector ${\epsilon }_n{}^q={\left[{\epsilon }_1^q, {\epsilon }_2{}^q,\dots , {\epsilon }_{2N}^q\right]}^T$ we have

$$J\left(Q_{n+1}{}^q\right)\approx \left[G\left(Q_{n+1}^q+{\epsilon }_{n+1}^q\right)-G\left(Q_{n+1}^q-{\epsilon }_{n+1}^q\right)\right]/2{\epsilon }_{n+1}^q.$$

(44)

5. 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. Modal analysis for double-tapered cantilever beams.

	Mabie et al. [50]		Relative Error with FEM (% ×10−3)
frequency	α = 5.0	α = 10.0	α = 5.0	α = 10.0
Ω1	30.9820	72.0487	0.1451	0.0231
Ω2	91.9273	186.802	0.0410	0.2007
Ω3	199.1682	371.238	0.0197	0.1413
Ω4	356.2088	635.049	0.0198	0.1609
Ω5	564.1394	981.657	0.0835	0.5628

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 \mathrm{m}$ of variable thickness under tip load forcing $F=-100 \mathrm{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].

5.1. 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 Figure 6 and Figure 7 the relative error concerning the calculated thrust coefficient and the efficiency of a rigid flapping foil against the time discretization parameter $\Delta t/{\rm 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 $\lambda =U\Delta t/\Delta 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 $\lambda =U\Delta t/\Delta x=3.5$,$\Delta t/{\rm T}=0.35\%$, the error is of the order of 0.5%, and thus the latter values are selected for the simulations.

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 $\lambda =U\Delta t/\Delta x=2.8,$$\Delta t/{\rm 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 $\Delta t/{\rm 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.

5.2. 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 ${\mathrm{x}}_R=1/3c$ and the following kinematic parameters $Str=0.3, {\alpha }_{eff}={15}^{\circ }, U=0.4 \mathrm{m}/\mathrm{s}, h_o/c=0.75, \psi ={90}^{\circ }.$ 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(\cdot {10}^7) \mathrm{Pa},$ Poisson’s ratio $v=0.4$ and material density ${\rho }_s=1100 {\mathrm{kg}/\mathrm{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 Figure 9 and Figure 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.

5.3. 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 \mathrm{Nm}$, ${\rho }_s=1200 {\mathrm{kg}/\mathrm{m}}^3$ is actuated at the LE with heaving amplitude $h_o/c=0.033,$ oscillating frequency within the interval $\omega /{\omega }_o=[0.3,8]$ and $\mathrm{Re}=6000$. For a plate immersed in fluid, see [34], it is reported that the first resonance frequency is equal to ${\omega }_o=4.71 \mathrm{rad}/\mathrm{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 Figure 11, Figure 12, Figure 13 and Figure 14. The trailing edge/leading edge (TE/LE) amplitude response $(A_{TE}/A_{LE})$ as a function of the non-dimensional frequency $\omega /{\omega }_o$ 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.

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 $\omega /{\omega }_o>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 $\omega /{\omega }_o=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].

5.4. 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.2. To further investigate the effect of flexibility on Froude efficiency and thrust, additional simulations were performed, and results are shown and discussed in this section. For this study, a NACA 0012 hydrofoil with $c=0.1$ and center of rotation at 1/3c is considered. The kinematic parameters of the flapping foil are the following $Str=0.3, U=0.4 \mathrm{m}/\mathrm{s}, h_o/c=0.75, \psi ={90}^{\circ }.$ The average thrust coefficient and Froude efficiency are studied as functions of the non-dimensional flexural rigidity $E/({\rho }_{s}gc)$ for a selection of pitch amplitudes ${\theta }_o=\{0^{\circ }, {10}^{\circ }, {20}^{\circ }, {30}^{\circ }\}.$ Results are presented in 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 ${\theta }_o=30\mathrm{deg}$, 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 ${\theta }_e=\theta (t)={\tan }^{-1}\left[(z_{LE}-z_{TE})/(x_{TE}-z_{LE})\right]$.

For the most stiff foil with $E\sim {10}^8\mathrm{Pa}$, the estimated maximum effective angle of attacks, that correspond to ${\theta }_o=\{0^{\circ }, {10}^{\circ }, {20}^{\circ }, {30}^{\circ }\}$ are $a_m=\{{43}^{\circ }, {33}^{\circ }, {23}^{\circ }, {13}^{\circ }\}.$ 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}^{\circ }$ corresponding to zero pitching amplitude for the most flexible foil $E\approx {10}^5 \mathrm{Pa}$ explains the behavior of the results in Figure 15. A qualitative explanation is also provided in Figure 16, where typical configurations of the instantaneous curvature of the flexible foil compared to the camber line of the rigid equivalent are presented on the trajectory of the LE presented.

In Figure 17, numerical results are presented to illustrate the effect of chordwise flexibility on the thrust and Froude efficiency over a range of Strouhal numbers with heaving amplitude $h_o/c=\{0.25, 0.5, 0.75, 1.0\}$ as a parameter. The simulations are performed for a flapping NACA 0012 hydrofoil with $c=0.12 \mathrm{m}, x_R=1/3c$ and the rest of parameters: $U=0.3 \mathrm{m}/\mathrm{s}, {\theta }_o={10}^{\circ }, \psi ={90}^{\circ }.$ The chosen material corresponds to one of the most elastic ones examined before, characterized by $E=3.45\cdot {10}^5 \mathrm{Pa}.$ It is clearly observed in this figure that chordwise flexibility leads to decrease of thrust as Strouhal number increases, especially for higher amplitudes of heaving motion. On the contrary, flexibility enhances the propulsive efficiency with increasing amplitude of heaving motion and higher value of Strouhal number.

6. 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 Figure 6, Figure 7 and Figure 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 Figure 9 and Figure 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 Figure 11, Figure 12, Figure 13 and Figure 14. 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 Figure 15 and Figure 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 ${\theta }_o=30\mathrm{deg}$, 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].