1. Introduction
Ocean waves contain massive amounts of energy and a suitably designed submerged flapping-foil system can be used to directly convert the energy of the waves into propulsion for a surface vehicle. Many different wave-assisted propulsion (WAP) system designs have been explored (Bøckmann & Steen Reference Bøckmann and Steen2014; Yang, Shi & Wang Reference Yang, Shi and Wang2019; Qi et al. Reference Qi, Jiang, Jia, Zou and Zhai2020; Xing & Yang Reference Xing and Yang2023) and here we focus on a simple WAP system which includes a single submerged flapping hydrofoil which is attached at a torsional spring loaded rotational hinge by a inflexible strut to a surface vehicle. As the surface vehicle moves up and down with the waves, the submerged foil attached via the strut to the surface vehicle also heaves up and down. A pitching motion of the foil is then generated ‘passively’ under the combined effect of hydrodynamic forcing and the restoring torque of the torsional spring. This combined pitching and heaving of the foil generates thrust which can be used to propel the surface vehicle or to augment the available propulsion system on the vehicle. The performance of the WAP system is therefore determined by the combined action of fluid–structure interaction and the hydrodynamics associated with pitching–heaving foils.
Combined pitching and heaving of control surfaces is commonly encountered in biological and bioinspired propulsion systems. Indeed, swimming using pectoral (Mittal et al. Reference Mittal, Dong, Bozkurttas, Lauder and Madden2006; Bozkurttas et al. Reference Bozkurttas, Mittal, Dong, Lauder and Madden2009; Dong et al. Reference Dong, Bozkurttas, Mittal, Madden and Lauder2010), caudal (Xiong & Liu Reference Xiong and Liu2019; Seo & Mittal Reference Seo and Mittal2022) and other fins essentially involves a combined pitching and heaving of the fin in a periodic manner. Similarly, the flapping flight of insects (Sane & Dickinson Reference Sane and Dickinson2002; Birch & Dickinson Reference Birch and Dickinson2003) and birds (Muijres et al. Reference Muijres, Johansson, Barfield, Wolf, Spedding and Hedenstrom2008; Song, Luo & Hedrick Reference Song, Luo and Hedrick2014; Jurado et al. Reference Jurado, Arranz, Flores and García-Villalba2022) is typically accomplished with combined pitching and heaving oscillations of the wings. Thus, the flow physics and dynamics of WAP systems are expected to have similarities to these biological propulsion systems. There are, however, two important differences. The first is the centrality of fluid–structure interaction (FSI) in WAP systems that drives the pitching motion. While FSI may be important for biological propulsion as well, pitching of wings and fins during biological propulsion are usually driven actively. Indeed, pitch in flying animals is essential for flight and manoeuvring (Haque et al. Reference Haque, Cheng, Tobalske and Luo2023; Liu, Wang & Liu Reference Liu, Wang and Liu2024). The second difference, which is a consequence of the first, is that, while in biological propulsion the phase difference between pitching and heaving is usually prescribed (and quite often, is 
$90^{\circ}$), for WAP systems, this phase difference is determined by the intrinsic properties of the system and can be quite different from 
$90^\circ$. The third possible difference is due to the fact that ocean waves can be highly stochastic and, since the entire system is driven by ocean waves, the dynamics of WAP systems could be highly stochastic as well. In contrast, flapping flight and animal swimming generally (but not always) involve control surface motions that are quite regular and often periodic.
The dynamics of WAP systems also has similarities to aeroelastic (or hydroelastic) flutter, which also manifests as combined pitching and heaving oscillations of control surfaces (Goyaniuk et al. Reference Goyaniuk, Poirel, Benaissa and Amari2023; Turner, Seo & Mittal Reference Turner, Seo and Mittal2023). The key difference between this phenomenon and WAP systems is that, while in elastic wing flutter both heave and pitch degrees-of-freedom (DoFs) are determined by FSI, in some WAP systems (including the one studied here), the heave is considered to be prescribed, and only the pitch DoF is determined by FSI. We note that there are other WAP system designs with elastically mounted struts/tethers, where the heave DoF is also subject to FSI effects. Another fundamental difference between wing/foil flutter and WAP systems is that, while in the former the objective of analysis is to determine causal mechanisms and reduce flutter, the objective in the study of WAP system is usually to find ways of increasing thrust generation.
Wave-assisted propulsion systems can also be contrasted to flapping-foil power generation systems (‘flapping foil turbines’) (Young, Lai & Platzer Reference Young, Lai and Platzer2014). In these systems, a foil capable of heaving and pitching is placed in a current and the heaving motion of the foil induced by the lift forces generated by the flow is used to extract energy from the flow. However, in these systems, pitch is controlled/driven via some mechanism and heave results from flow-induced forces, and this is the opposite of WAP systems. Recently, flapping-foil turbines with prescribed heave motion and passive-pitch motion have also been introduced (Boudreau, Gunther & Dumas Reference Boudreau, Gunther and Dumas2019). However, these flapping-foil turbines operate at very large pitch angles, whereas, as will be shown in this paper, WAP systems operate best at low to moderate pitch angles. Thus, the hydrodynamics of these two systems are quite different.
Several studies have been conducted on passive-pitch WAP systems similar to the one under consideration here (i.e. a system where heave is driven directly by wave motion and passive pitching is determined by torsional spring). Bøckmann & Steen (Reference Bøckmann and Steen2014) performed experiments and pointed out that a passive-pitch WAP system can obtain a higher efficiency as compared with a pitch-controlled flapping hydrofoil. Utilizing a boundary element method, Thaweewat, Phoemsapthawee & Juntasaro (Reference Thaweewat, Phoemsapthawee and Juntasaro2018) investigated the propulsion performance of a passive-pitch WAP system, and suggested that such a system has a function similar to an adjustable pitch propeller. Using computational fluid dynamics, Qi et al. (Reference Qi, Zou, Lu, Shi, Li, Qin and Zhai2019, Reference Qi, Jiang, Jia, Zou and Zhai2020) studied the effects of spring stiffness, reduced frequency and mass on a passive-pitch WAP configuration. Yang et al. (Reference Yang, Shi, Zhou, Guo and Wang2018, Reference Yang, Shi and Wang2019) investigated the hydrodynamic response of a flow-induced pitching foil, indicating that the torsional spring plays a key role in propulsion performance. Zhang et al. (Reference Zhang, Feng, Chen and Gao2022) studied the effect of pivot location on the semi-active flapping foil, demonstrating that the pivot location at 33 % chord might fall into a pitch unstable state at low oscillating frequency, leading to low efficiency.
Despite the growing interest in these systems, the physics behind thrust generation from WAP systems, including the effect of torsional spring stiffness and the pivot location on system behaviour as manifested through pitch amplitude and the phase difference between pitch and heave, which determines the thrust performance, is not fully understood. Such an analysis is complicated by the fact that FSI simulations are relatively expensive and introduce additional parameters. Adding to this is the complexity of the flow that is generated by these foils, which makes it difficult to pinpoint the hydrodynamic mechanisms that are key to the performance of the foil.
As pointed out earlier, ocean waves can be irregular and stochastic to a varying degree. It is expected that the addition of the stochastic component to the heave motion will modify the vortex dynamics (leading-edge vortex and wake vortices) of the flapping-foil WAP system and thereby change the hydrodynamics as well as thrust performance of the foils. The current study with sinusoidal heave motion is to be considered a starting point towards this more realistic situation. Sinusoidal heaving motion has been employed in several earlier studies related to flapping-foil propulsion (Thaweewat et al. Reference Thaweewat, Phoemsapthawee and Juntasaro2018; Yang et al. Reference Yang, Shi, Zhou, Guo and Wang2018; Qi et al. Reference Qi, Zou, Lu, Shi, Li, Qin and Zhai2019; Yang et al. Reference Yang, Shi and Wang2019; Qi et al. Reference Qi, Jiang, Jia, Zou and Zhai2020) and this also provides a useful context for the current results on WAP systems.
In the present study, we combine high-fidelity computational fluid dynamics with three tools to examine the hydrodynamics and thrust performance of these systems. The first is ‘energy maps’ that employ prescribed kinematics simulations to understand the energy exchange between the foil and the fluid. The second is the force-partitioning method (FPM), which is a data-enabled method for quantifying the contributions of various flow features to the pressure-induced forces (Menon & Mittal Reference Menon and Mittal2021a,Reference Menon and Mittalb). The third tool is a new phenomenological model that recognizes the dominant contribution of the leading-edge vortex to the pressure forces on the foil and provides a simple means of predicting the performance of these flapping foils. We use these tools to examine the effect of three key system parameters on the performance of the WAP system – the pivot axis location, the torsional spring stiffness and the flapping frequency.
2. Methodology
2.1. Problem set-up
The two-dimensional computational model in the current study employs a rigid NACA0015 hydrofoil with a slightly rounded trailing edge, immersed in an incompressible flow with free-stream velocity 
$U_{\infty }$. The slight rounding of the trailing edge ensures that the flow is well resolved around it, and is has been verified that this has no significant effect on the hydrodynamic characteristics of the hydrofoil (Menon & Mittal Reference Menon and Mittal2019). The governing incompressible Navier–Stokes equations expressed in the dimensionless form are as follows:
$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} ={-}\boldsymbol{\nabla p} + \frac{1}{Re} \nabla^2 \boldsymbol{u} \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}$$
where 
$Re=\rho U_{\infty }C/\mu$ is the chord-based Reynolds number, C is the chord length of the foil, 
$\rho$ is the density of the fluid, 
$\mu$ is the viscosity, u is the velocity and p is the pressure. The hydrofoil is subjected to heaving with a prescribed amplitude and pitching is generated via flow-induced action on the torsional spring mounted hinge as shown in figure 1. The leading edge of the foil is located at a distance of 
$4.5C$ from the inlet in the streamwise direction and the mean heaving location of the foil is kept at the centre of the domain in the cross-stream direction.
(a) Schematic of the hydroelastic system used in this study; (b) computational domain and close-up of the Cartesian computational grid.
The linear torsional spring has a spring constant 
$k$ and the rotational hinge (or the pivot axis) is at a non-dimensional chordwise location of 
$X_e^{*} = X_e/C$. Here, 
$X_e$ is measured from leading edge of the foil along the chord line towards the trailing edge. Thus, a negative value of 
$X_e$ corresponds to pivot locations ahead of the leading edge. The equilibrium angular position of the spring is denoted by 
$\theta _i$, which has been kept zero across all the simulations in this work. The hydrofoil is prescribed to heave in the 
$y$-direction given by equation 
$h(t)=h_0 \sin (2{\rm \pi} f_h t)$, where 
$h_0$ is the amplitude of heaving motion and 
$f_h$ is the frequency of heaving motion (non-dimensionally written as 
$H^*=H^*_0 \sin (2 {\rm \pi}\,St\, T^*)$, with 
$St=f_h C/U_\infty$, 
$T^*=tU_\infty /C$ and 
$H^*=h/C$). The flow-induced pitching motion is governed by the equation of a forced angular spring–mass oscillator. The equation is non-dimensionalized by the characteristic variables used for the flow (length: 
$C$; time: 
$C/U_{\infty }$), and the resulting non-dimensionalized equation governing the dynamics in pitch of the hydrofoil is given by
\begin{equation} I^* \ddot{\theta}+k^*(\theta - \theta_i) = C_M, \end{equation} where 
$C_M = M/(0.5 \rho U^2_{\infty } C^2 b)$ is the coefficient of pitching moment, and 
$I^*=2I/(\rho C^4 b)$, 
$k^*= 2k/(\rho U^2_{\infty } C^2 b)$ are the dimensionless moment of inertia and spring stiffness, respectively, 
$b$ is the spanwise width of the foil, M is the moment of inertia about the pivot location, k is the spring stiffness and 
$\theta$ is the pitch angle. Here, 
$C_M$ and 
$I^*$ are calculated with respect to 
$X_e^*$. The moment of inertia 
$I$ is calculated using the following volume integral 
$\rho \int _V r^2 \,{\rm d}V$, where 
$r$ is the radial distance from the pivot location 
$X_e$ and V represents volume. In this work, the torsional spring stiffness is expressed in terms of a non-dimensional frequency ratio 
$f_{\theta }/f_h$, where 
$f_{\theta }=(1/2{\rm \pi} )\sqrt {k/I}$ is the natural frequency of the spring–mass system and 
$f_h$ is the frequency of the heaving motion in the 
$y$-direction. The above equation can be written in terms of 
$f_{\theta }$ as
\begin{equation} \ddot{\theta}+\left(2{\rm \pi} \frac{f_\theta}{f_h} \,St\right)^2(\theta - \theta_i) = \frac{C_M}{I^*}. \end{equation} The above configuration has seven governing parameters (
$Re$, 
$I^*$, 
$X_e^*$, 
$f_{\theta }/f_h$, 
$\theta _i$, 
$H_0^*$, 
$St$) and this is indicative of the inherent complexity of this configuration. In the current study, we fix the values of 
$Re$, 
$I^*$, 
$\theta _i$ and 
$H_0^*$ to 10 000, 0.14, 0 and 0.5, respectively, and explore the effect of frequency ratio (
$\,f_{\theta }/f_h$), which is a surrogate for spring stiffness, the pitch axis location 
$X_e^*$ and Strouhal number 
$St$. The Reynolds number chosen is high enough so as to generate a robust vortex shedding phenomenon, which drives flutter, but is also low enough so as to allow resolved simulations at a reasonable computational expense. The choice of the moment of inertia about the elastic axis corresponds to a solid-to-fluid density ratio of 5. While the effect of varying 
$I^*$ in this system could be of significant engineering interest, we choose to use a fixed value throughout this study to reduce the number of independent parameters. However, it must be noted that varying 
$I^*$ is partially equivalent to changing the natural frequency as 
$f_{\theta } \sim 1/\sqrt {I^*}$, which is an important parameter in this study. For 
$St=0.2$, five different pitch-axis locations have been studied with 
$X_e^*=-0.1$, 0, 0.1, 0.25 and 0.5. For 
$X_e^*=0.1$, the effect of higher and lower Strouhal numbers (
$St = 0.3$ and 0.1) is also studied. The parameter 
$f_\theta /f_h$ represents the torsional spring stiffness (a smaller 
$f_{\theta }/f_h$ corresponds to a softer spring) and this parameter is varied from about 1 to 3 for all 
$X_e^*$ and 
$St$ considered.
2.2. Numerical method
The coupled fluid–structure system in this study has been simulated using the sharp-interface immersed boundary method solver ‘ViCar3D’ described in Mittal et al. (Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and Von Loebbecke2008) and Seo & Mittal (Reference Seo and Mittal2011). This method is particularly useful for flows around complex moving and deforming solid geometries as it allows for use of a simple body-non-conformal Cartesian grid to capture shapes and motions of the immersed body. This method leads to very accurate computations of surface quantities due to its ability of preserve a sharp interface around the immersed boundary (Mittal et al. Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and Von Loebbecke2008). The fluid equations are solved by using a second-order fractional-step method, and the pressure Poisson equation is solved using the bi-conjugate gradient method. Second-order finite differences are used for all spatial derivatives, and the time integration is performed using a second-order Adams–Bashforth and Crank–Nicolson methods. ViCar3D has been validated and verified extensively for a variety of stationary and moving boundary problems in earlier papers (Mittal et al. Reference Mittal, Dong, Bozkurttas, Najjar, Vargas and Von Loebbecke2008; Seo & Mittal Reference Seo and Mittal2011). More relevant to the current work is the computational study of flow-induced pitching and forced pitching oscillations of a foil in Menon & Mittal (Reference Menon and Mittal2019).
The FSI employs a loosely coupled approach wherein the flow and dynamical equations are solved sequentially. The aerodynamic forces/moments are calculated on the Lagrangian marker points on the surface of the solid body, and passed on to the solid dynamical equation (2.3). The angular velocity predicted from this equation is subsequently imposed on the marker points. The hydroelastic system is immersed in a large 
$10C\times 10C$ computational domain, where the hydrofoil is placed 4.5 chord lengths from the downstream boundary, and the isotropic grid resolution around the solid body corresponds to approximately 250 points along the chord. The grid is stretched in all directions away from the rectangular region that surrounds the foil and the near wake, resulting in a baseline grid of 
$640 \times 600$ points. A Dirichlet velocity boundary condition is used at the inlet boundary of the domain, and zero-gradient Neumann conditions are specified at all other boundaries. Grid refinement and domain dependence studies described in Appendix confirm that the results on this grid are well converged and domain-independent.
For all the flow-induced pitching motion cases simulated here, we initialize the foil to a zero angle-of-attack, 
$y=0$ (which represents the mean heaving location) and with zero initial angular and heaving velocity. For the prescribed pitching motion cases simulated here, the hydrofoil is initialized to 
$y=-h_0$ with zero heaving velocity. The initial angle and angular velocity are determined from the phase difference selected between the heaving and pitching motion 
$\psi$. The equilibrium angle of attack of the hydrofoil is set to zero. The constant free-stream flow is then imposed and the dynamics allowed to evolve until the system reaches a periodic state. For the cases simulated here, it takes between 10 and 15 oscillation cycles before this state is achieved for the flow-induced motion cases and approximately 5 cycles for the prescribed motion cases. Average quantities are computed for 6 cycles after the periodic state is reached.
2.3. Three-dimensional vs two-dimensional simulations
The foils employed in this study are assumed to be spanwise invariant and have a large (infinite) span but, despite the two-dimensional geometry, the flow over these foils at the chosen Reynolds number of 10 000 is expected to be intrinsically three-dimensional (3-D). However, the research approach to be employed here relies on using data from 
$O(100)$ simulations of these WAP foils and conducting this number of 3-D simulations is impractical. We instead employ 2-D simulations for the bulk of this research but we first verify that the 2-D simulations provide a reasonably accurate model for the 3-D flow. For this verification, 3-D simulations are performed with the same computational domain as that for the 2-D simulations (figure 1b) in 
$x$ and 
$y$. For the 3-D simulations we choose a span length of 2
$C$ and employ a uniform grid with 100 grid points in the spanwise direction.
These comparisons are made for prescribed kinematics simulations with the assumption that the overall 3-D nature of the flow is determined by the foil kinematics and not by the presence or absence of FSI. We compare the 2-D and 3-D simulations for three representative cases with 
$X^*=0.1$: the first two cases are at a high pitch amplitude of 
$\theta _\circ =25^\circ$ but with two different pitch–heave phases of 
$\psi =-180^\circ$ and 
$0^\circ$. These high amplitude cases are expected to generate very robust 3-D features and serve as a strict verification test for 2-D models. The third case is a low amplitude case with 
$\theta _\circ =5^\circ$ and 
$\psi =-70^\circ$), which as later simulations will show, is one of the high thrust generation cases.
Figure 2 shows several snapshots of the vortex structures for the case with 
$\theta _\circ =25$ and 
$\psi =-180$ and the flow clearly exhibits a robust 3-D structure. We note, however, that, while the vortices closer to the trailing edge and in the wake are highly three-dimensional, the vortices in the front region of the foil, including the leading-edge vortex (LEV), are mostly two-dimensional.
Contours of vorticity in the 
$z$ direction on the isosurface of non-dimensional 
$Q$ value 
$=10$ for 
$\theta _\circ = 25^\circ$ and 
$\psi = -180^\circ$ at 
$t/T_p =$ (a) 0, (b) 0.25, (c) 0.5 and (d) 0.75.
Figure 3 shows the comparison of the coefficients of drag (
$F_x$), lift (
$F_y$) and pitching moment (
$M_z$) between the 2-D and 3-D simulations for all the three cases and we note the very close match between the two sets of simulations. As will be shown later in the paper, the LEV dominates the generation of pressure forces on the flapping foil and the observation that the LEV remains mostly two-dimensional in its initial evolution provides an explanation for this concordance between the 2-D and 3-D predicted forces on the foil. The conclusion from this comparison is that 2-D simulations provide an acceptable level of fidelity for modelling the essential flow features and hydrodynamic force generation for these flapping foils and we therefore proceed with using 2-D simulations for the rest of this study.
Comparison of forces and moments between two-dimensional and three-dimensional simulations for several cases; (a) 
$\theta _\circ = 25^\circ$, 
$\psi = -180^\circ$, (b) 
$\theta _\circ = 25^\circ$, 
$\psi = 0^\circ$, (c) 
$\theta _\circ = 5^\circ$, 
$\psi = -70^\circ$.
2.4. Energy maps
Energy maps are tools that can be used to examine FSI and aeroelastic flutter via forced kinematics simulations (Morse & Williamson Reference Morse and Williamson2009; Bhat & Govardhan Reference Bhat and Govardhan2013; Menon & Mittal Reference Menon and Mittal2019). In this approach, simulations are conducted with prescribed kinematic parameters and the net energy imparted by the fluid to the elastically supported foil/body is estimated. Positive (negative) values of this quantity denote regions of the parameter space where the flow-induced oscillations are expected to grow (decay) and zero values denote equilibrium conditions. For the current WAP foil, where heave is prescribed and the pitch oscillations are flow induced, the non-dimensional energy imparted to the foil (
$E_f^*$) is given by
\begin{equation} E_f^* = \int^{t_n+T_p^*}_{t_n} C_M \dot{\theta} \,{\rm d}t, \end{equation}
where 
$\dot {\theta }$ is the angular velocity and 
$T_p^*$ is the non-dimensional period of oscillation. Maps of this exchanged energy can be made over a range of pitch amplitudes and other chosen parameters of interest wherein the above integral is estimated from the computed pitching moment on the foil. The energy maps so constructed provide the ability to predict growth, decay and limit-cycle behaviour of the flow-induced system. The validity of these energy maps obtained from forced oscillations for flow-induced oscillations rests primarily on the assumption that the free pitching oscillations are strictly sinusoidal at the system frequency. We examine this assumption for select cases and employ energy maps to gain insights into the dynamics of this coupled system.
2.5. Force-partitioning method
In order to quantify vortex-induced pressure force on bodies within fluid flows, the FPM is utilized (Menon & Mittal Reference Menon and Mittal2021b,Reference Menon and Mittalc). The objective is to quantify the different mechanisms that contribute to the pressure-induced hydrodynamic forces on a body with surface 
$B$, which is immersed in a flow domain given by 
$V_f$ as shown in figure 4. The domain is bounded externally by the surface 
$\varSigma$, and the unit normal pointing outward from the fluid domain into the internal and external boundaries, 
$B$ and 
$\varSigma$, is given by 
$\hat {\boldsymbol{n}}$.
Schematic of the set-up for the FPM, along with relevant symbols.
In order to quantify the pressure-induced loads in the 
$i$ direction on the body, the starting point for FPM is the calculation of an ‘influence field’, 
$\phi _i$, which satisfies the following equation:
\begin{equation} \nabla^2\phi_i = 0, \quad \text{in } V_f \text{ with } \hat{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\nabla} \phi_i =\left\{\begin{array}{@{}ll@{}} n_i, & \text{on } B \\ 0, & \text{on } \varSigma \end{array}\right.\quad \text{for } i=1,2. \end{equation}
Here, 
$n_i$ is the component of the unit vector 
$\hat {\boldsymbol{n}}$ in the 
$i$ direction. In the current problem, 
$i=1$ corresponds to the direction of thrust on the wing and the influence field of interest is therefore 
$\phi _1$. Note that 
$\phi _i$ is a surface-dependent field that does not depend on the flow field and can be obtained using standard numerical techniques for any given surface.
By performing an inner product between the momentum equation and the gradient of the influence field, integrating over the volume of the fluid domain and applying the incompressibility constraint and the divergence theorem, we obtain (as shown in detail in Menon & Mittal Reference Menon and Mittal2021c the following decomposition of the pressure-induced force on the control surface:
\begin{align} F_i =\int_B p n_i \,{\rm d}S &={-}\overbrace{\rho \int_B \hat{\boldsymbol{n}}\boldsymbol{\cdot}\left(\frac{{\rm d}\boldsymbol{U}_B}{{\rm d}t} \phi_i\right){\rm d}S}^{F_{KINE}} -\overbrace{2 \rho \int_{V_f} Q\phi_i \,{\rm d}V}^{F_{VIF}} + \overbrace{\mu \int_{V_f} (\nabla^2 \boldsymbol{u})\boldsymbol{\cdot}\boldsymbol{\nabla} \phi_i \,{\rm d}V}^{F_{VIS}}\nonumber\\ &\quad - \overbrace{\rho \int_{\varSigma} \hat{\boldsymbol{n}}\boldsymbol{\cdot}\left(\frac{{\rm d}u}{{\rm d}t} \phi_i\right) {\rm d}S}^{F_{BOUND}} \quad \text{for } i = 1,2, \end{align}
where 
$\boldsymbol{U}_B$ is the local velocity of the immersed surface. The four terms on the right-hand side of (2.7) represent different components of the pressure-induced force on the surface of the body. The first term 
$F_{KINE}$ contains the linear acceleration reaction force (or added mass in potential flow) (Menon & Mittal Reference Menon and Mittal2021a) as well as the centripetal acceleration reaction (Zhang, Hedrick & Mittal Reference Zhang, Hedrick and Mittal2015). The second term 
$F_{VIF}$ is referred to as the ‘vortex-induced force’ and is associated with the force exerted by flow structures within the fluid like the vortex or region surrounding a vortex. This is found to be a dominant component in a variety of flows including those over flapping wings (Zhang et al. Reference Zhang, Hedrick and Mittal2015), bluff bodies (Menon & Mittal Reference Menon and Mittal2021a), pitching foils (Menon & Mittal Reference Menon and Mittal2021b) and delta wings (Li, Zhao & Graham Reference Li, Zhao and Graham2020). The third term 
$F_{VIS}$ is the pressure force associated with viscous diffusion in the flow and this term is usually significant only for low Reynolds number flows. The fourth term 
$F_{BOUND}$ is associated with flow acceleration at the outer boundary of the domain and is negligible for large domains (Zhang Reference Zhang2015) and steady free-stream flows. More details on these forces can be found in the work of Zhang et al. (Reference Zhang, Hedrick and Mittal2015), Zhang (Reference Zhang2015) and Menon & Mittal (Reference Menon and Mittal2021a,Reference Menon and Mittalb,Reference Menon and Mittalc). As the main focus of this work is on thrust generation, we implement the above method in the 
$x$-direction to identify the main thrust generating mechanism.
3. Results
3.1. Flow-induced pitching motion of hydrofoil
For the flow-induced motion simulations, a sinusoidal heaving motion is prescribed for the hydrofoil (and given as 
$h^*=0.5 \sin (0.4{\rm \pi} T^*)$ with 
$h^*_0=0.5$) and the hydrofoil is allowed to pitch passively based on the pitching moment imposed on it by the flow and the torsional spring. We begin by presenting in figure 5 the results for the case with the pitch axis location 
$X_e^*=0.1$, a Strouhal number of 0.2 and a frequency ratio 
$f_{\theta }/f_h$ of 3, which corresponds to a relatively stiff spring. Figure 5(a) shows the variation of pitch and heave for this case after the simulation has achieved a periodic state and it shows that the pitching exhibited by the foil is mostly sinusoidal for this case. Given that the maximum heave nearly coincides with a zero-pitch angle, the overall phase angle between pitch and heave (
$\psi$) is expected to be close to 
$-90^\circ$ and, indeed, the calculated phase angle is found to be 
$-101^\circ$. The time series for the coefficients of lift (
$C_L=F_y/\frac {1}{2}\rho U^2_\infty C$), thrust (
$C_T=-F_x/\frac {1}{2}\rho U^2_\infty C$) and moment (
$C_M=M_z/\frac {1}{2}\rho U^2_\infty C^2$) corresponding to the same oscillation cycle are shown in figure 5(b).
Data for a representative case of flow-induced pitching oscillations with 
$X_e^*=0.1$ and 
$f_{\theta }/f_h=3$. Time variation of (a) heaving and pitching motions and (b) coefficient of lift (
$C_L$), coefficient of thrust (
$C_T$) and coefficient of pitching moment (
$C_M$) during one cycle of oscillation (corresponding to the cycle shown in the snapshots); (1–4) snapshots of the flow field coloured by contours of 
$z$-vorticity.
As the hydrofoil starts to heave in the upward direction, it experiences an increasing clockwise (negative) moment resulting in the pitching up of the foil. However, this pitch-up and heave-up motion increases the effective angle-of-attack of the foil and leads to the formation of LEVs on the suction surface of the hydrofoil, as seen in figure 5(c). These LEVs generate a suction pressure on the adjacent surface and this generates a positive thrust and a negative lift component. Consequently, thrust increases and lift become increasingly negative as the foil moves upwards. During the downstroke, the LEV formation process switches to the top surface, resulting in the generation of a positive lift and a positive thrust (figure 5(4)).
Figure 6(a) shows the foil kinematics for the two other frequency ratios (
$\,f_{\theta }/f_h = 1,2$) for this particular configuration. The decreasing spring stiffness associated with a reducing 
$f_{\theta }/f_h$ results in an increase in the pitch amplitude of the hydrofoil (pitch amplitude values are 
$24.6^\circ$, 
$14.2^\circ$ and 
$7.5^\circ$ for 
$f_{\theta }/f_h=1$, 2 and 3, respectively). There is also a significant change in the pitch–heave phase difference (
$\psi$), which is 
$-139^\circ$, 
$-104^\circ$ and 
$-101^\circ$ for these cases, respectively. Thus, the spring stiffness has a significant impact on the amplitude as well as phase of pitching, which is expected to have a concomitantly significant effect on the hydrodynamic forces and moment. Indeed, the mean thrust coefficients for 
$f_{\theta }/f_h=1$, 2 and 3 are estimated to be 
$-$0.06, 0.16 and 0.18, indicating that thrust increases monotonically with the spring stiffness. This behaviour will be explored and explained further in § 3.2. The pitching movement of the foil is also highly sinusoidal for these two additional cases and this is confirmed quantitatively from the Fourier transform of the pitching motion in figure 6(b), which shows a dominant peak at the flapping frequency.
(a) Results for flow-induced pitching motion of hydrofoil for 
$St =0.2$, 
$X_e^*=0.1$ and 
$f_{\theta }/f_h$ equal to 1, 2 and 3. Panel (b) shows the Fourier transform of the pitching motion for 
$X_e^*=0.1$ and different 
$f_\theta /f_h$. Inset in panel (b) shows the pitch-axis location. Similarly, panels (c) and (d) show the flow-induced pitching motion for 
$X_e^*=0$ and 
$X_e^*=-0.1$, respectively, and different 
$f_\theta /f_h$.
The location of the pitch axis is expected to be an important design parameter for WAP foil simulations and, in addition to the nominal case with 
$X^*_e=0.1$, simulations are also conducted for 
$X^*_e=-0.1, 0, 0.25$ and 0.5 with varying 
$f_\theta /f_h$. Figure 6(c,d) shows the results for 
$X^*_e=0.0$ and 
$-0.1$ and, qualitatively, these two cases seem to show a behaviour very similar to that of the 
$X^*_e=0.1$ case. All these cases generate pitch oscillations of similar magnitude and the oscillations are also quite sinusoidal. In table 1 we have tabulated all the key metrics for all the FSI cases simulated here and this confirms that, even quantitatively, these three cases generate very similar motion (as show by 
$\theta _o$ and the hydrodynamic thrust).
Data on the key metrics for the flow-induced pitching cases. Here, ‘IND’ indicates that the value is not determinable due to a high degree of non-sinusoidal motion and ‘rms’ indicates the root mean square value.
The results for the remaining two cases where the pitch axis is located further down the chord exhibit a very different behaviour. For 
$X_e^*=0.25$, the WAP foil exhibits very high amplitude oscillations with 
$\theta _0=63.8^\circ$ for 
$f_\theta /f_h=1$. As 
$f_\theta /f_h$ is increased to 2 and 3 (figure 7a), the pitching motion becomes chaotic and a periodic state is not achieved, even after 30 cycles. This case also exhibits a very significant drop in thrust when compared with 
$X^*_e=0.1$ (see table 1). For a pitch axis located at mid-chord, i.e. 
$X_e^*=0.5$, 
$f_\theta /f_h=1$ again results in very high amplitude oscillations but the system does not reach a periodic steady state. For higher 
$f_\theta /f_h$ (i.e. 2 and 3 (figure 7b), the system is able to reach a periodic steady state but the pitching motion is now ahead of the heaving motion, in contrast to all other 
$X_e^*$, where the pitching motion lags behind the heaving motion. This leads to the generation of high drag with thrust coefficients of 
$-$1.85 and 
$-$1.56 for 
$f_\theta /f_h=2$ and 3, respectively.
Results from FSI simulations of hydrofoil for 
$St =0.2$, 
$f_\theta /f_h=3$ and (a) 
$X_e^*=0.25$, (b) 
$X_e^*=0.5$.
This complex behaviour for 
$X_e^*=0.25$ is associated with the fact that for these foils, as the pitch axis gets close to the centre of pressure (CoP) (the CoP is estimated from our simulations to be located at around 
$x/C = 0.26$ for these foils), the pitch moment can undergo large variations and even undergo rapid changes in direction with the pitching and heaving of the foil. The coupling between the moment and motion of the foil can generate a complex flow behaviour that feeds back into the foil motion, resulting in chaotic motion. The pitch axis of 
$X_e^*=0.50$ is on the other side of the CoP and that results in the change of the pitching phase relative to heaving. The above set of simulations with different 
$X_e^*$ lead us to the conclusion that the pitch axis should be placed upstream near the leading edge and that the performance of the foil is quite insensitive to the precise location of the pitch axis as long as the above condition is satisfies.
The final set of FSI simulations are to examine the effect of Strouhal number on the foil dynamics. We note that the Strouhal number is determined by the forward velocity as well as the wave encounter frequency, and a WAP system could encounter a range of Strouhal numbers during its operation. We explore the effect of the Strouhal number by first increasing it by 50 % from its nominal value, i.e. to 
$St=0.3$, while keeping 
$X_e^*$ at 0.1 (figure 8a) and we find that there is a slight decrease in pitch amplitude and slight reduction in phase difference 
$\psi$ as compared with 
$St =0.2$. However, the important difference is that the thrust increases by 2.94 times and 1.46 times for 
$f_\theta /f_h=2$ and 3, respectively, as compared with 
$St =0.2$. The case 
$f_\theta /f_h=1$ also generates thrust (
$C_T=0.13$) here in contrast to 
$St =0.2$, which generates drag. The highest value of the thrust coefficient for this case (
$=$0.47) is found at 
$f_\theta /f_h=2$. We also note that the case with the highest stiffness exhibits pitch oscillations that are not purely sinusoidal. As will be shown later, this has implications for the predictions from the energy maps. Finally, simulations were also conducted for a 50 % lower Strouhal number, i.e. for 
$St = 0.1$, and the results are shown in figure 8(b) and we note all these cases generate small magnitude of drag or lift (table 1).
Sinusoidal heaving and flow-induced pitching motion of hydrofoil for (a) 
$St =0.3$, 
$X_e^*=0.1$, (b) 
$St =0.1$, 
$X_e^*=0.1$ and different 
$f_{\theta }/f_h$. Panel (c) shows the Fourier transform of the pitching motion for 
$St = 0.3$ and different 
$f_\theta /f_h$.
3.2. Energy and thrust maps
In order to construct energy maps that can provide insights into the flow-induced pitching of WAP foils, simulations are performed with a fixed heaving oscillation condition and a prescribed set of pitching motions. The pitching motion is prescribed in its non-dimensional form as 
$\theta (T^*)$ = 
$\theta _i +\theta _0 \sin (2{\rm \pi} \,St\, T^* + \psi )$, with 
$\psi$ being the phase difference between heaving and pitching. This configuration also has seven governing parameters (
$Re$, 
$X^*_e$, 
$\theta _i$, 
$H_0^*$, 
$\theta _0$, 
$St$ and 
$\psi$). As before, the values of 
$Re$, 
$\theta _i$ and 
$H_0^*$ are fixed at 10 000, 0 and 0.5, respectively. Given the results of the FSI simulations presented in the previous section, the frequency of pitching is kept the same as the frequency of heaving for these simulations. For 
$St =0.2$, five values of 
$X^*_e$ (i.e. 
$X^*_e = -0.1$, 0, 0.1, 0.25 and 0.5), six values of 
$\theta _0$ (i.e. 
$\theta _0 = 2.5^\circ$, 
$5^\circ$, 
$10^\circ$, 
$15^\circ$, 
$20^\circ$ and 
$25^\circ$) and 11 values of 
$\psi$ (i.e. 
$\psi =0^\circ$, 
$-30^\circ$, 
$-50^\circ$, 
$-60^\circ$, 
$-70^\circ$, 
$-80^\circ$, 
$-90^\circ$, 
$-100^\circ$, 
$-120^\circ$, 
$-150^\circ$ and 
$-180^\circ$) are employed. In addition, to identify the effect of the change in the Strouhal number, 66 simulations are performed for 
$St =0.3$ and 0.1 with 
$X^*_e$ kept at 0.1 taking six different 
$\theta _0$ and 11 different 
$\psi$ with values the same as before. A total of 462 prescribed pitching motion simulations are therefore conducted in order to generate these energy maps. Each case is simulated for 10 cycles to ensure that a periodic state is reached and 
$E^*_f$ is obtained by integrating over the last 5 cycles.
The pitching energy maps shown in the left column of figure 9 are for 
$St = 0.2$ and 
$X_e^* = 0.1$, 0.0 and 
$-$1.0. The dashed line in each pitching energy map indicates the zero pitching energy curve (
$E^*_f=0$) and represents a condition of equilibrium. As noted by Menon & Mittal (Reference Menon and Mittal2019), regions on the equilibrium curve where 
${\rm d}E_f^*/{\rm d} \theta _o < 0$ represent stable equilibria and the stationary states of FSI cases are expected to lie on these regions of stable equilibria. Indeed, as can be seen from the three plot, all three of the simulated FSI cases fall on the stable sections of the equilibrium curve, confirming the validity of these maps for these cases. It is also clear that these three cases exhibit very similar behaviour over the parameter space covered with the maps. We find that, as the spring stiffness (i.e. 
$f_\theta /f_h$) is increased for these cases, the amplitude of the oscillation decreases from approximately 
$25^\circ$ for 
$f_\theta /f_h = 1$ to 
$7.5^\circ$ for 
$f_\theta /f_h=3$. The map also indicates that further increase in spring stiffness could reduce the amplitude to approximately 
$5^\circ$ but further reductions are not possible while maintaining sinusoidal motion.
Pitch energy maps and thrust maps generated for 
$St =0.2$ and (a,b) 
$X_e^*=0.1$, (c,d) 
$X_e^*=-0.1$, (e,f) 
$X_e^*=0$. The dashed line in the energy and thrust map indicates the equilibrium curve. In the thrust map, the purple circles denote the boundary between thrust and drag and the dashed-dot curve indicates the region within 5 % of the maximum thrust value.
The pitch map also shows that, in the regions covered by the map, the equilibrium solutions achieve a phase angle between pitching and heaving ranging from 
$-150^\circ$ to 
$-95^\circ$, which is significantly different from 
$-90^\circ$, especially for the softer spring. Indeed, as mentioned earlier, this is one feature of WAP foils that makes them distinct from biomimetic flapping foils/wings since most of those configurations correspond to ‘normal’ flapping where the phase angle is 
$-90^\circ$.
The right column of figure 9 shows thrust coefficient maps generated from same set of simulations as in left column. The dashed line indicates the zero thrust coefficient curve (
$C_T=0$), with red (blue) contours indicating positive (negative) thrust values. Peak values of thrust coefficients are 0.20, 0.21 and 0.22 for the 
$X_e^* = 0.1$, 0.0 and 
$-$0.1 cases, respectively. For all three pitch-axis cases, high values of thrust can be found in the region with 
$\theta _0$ values between 
$7.5^\circ$ and 
$17.5^\circ$ and 
$\psi$ values that are nominally centred around 
$-90^\circ$. The equilibrium curves superposed on these plots show us the manifold of the FSI cases on this thrust map. For the 
$X_e^*=0.1$ case, the equilibrium curve does pass through the region of high thrust and this high thrust would correspond to a 
$f_\theta /f_h$ value between 2 and 3. For 
$X_e^*=-0.1$, the equilibrium curve does not pass through the region of highest thrust, indicating that, with this location of the pitch axis, the WAP foil cannot access the highest possible thrust that is available for this pitch axis.
Overall, the similarity between these three cases is quite remarkable, indicating that a pitch axis located close to the leading edge results in robust performance, with regular (sinusoidal) flapping, and thrust generation as long as the spring is not too soft, and the pitch amplitude is limited to lower values. All the thrust maps have a similar topology, with high thrust centred at intermediate pitch angles and for a phase difference between pitch and heave of 
$-90^\circ$. Cases with a combination of large pitch amplitude and a pitch–heave phase differences close to 
$0^\circ$ (or 
$-180^\circ$) tend to generate drag. For low spring stiffness (
$\,f_\theta /f_h$ = 1), the phase difference between pitching and heaving deviates significantly from 
$-90^\circ$, and this correlates with a diminished thrust performance.
Figure 10 shows the energy and thrust maps for cases with 
$X_e^* = 0.25$ and 0.5, which, as seen from the FSI simulations, exhibit a dynamics that is quite distinct from the first three cases. The thrust maps for these cases are quite similar in overall topology to the other three cases, suggesting that changes in pitch-axis location do not significantly change the underlying fluid dynamics and thrust production mechanisms. It is, however, in the pitch energy maps where we note the significant difference. For the 
$X_e^* = 0.25$, the manifold is dominated by negative pitch energy. Furthermore, the equilibrium curve is oriented nearly vertically and exhibits no significant regions of stability. This suggests that sinusoidal pitch oscillations would be difficult to achieve and sustain for this configuration and this is borne out by the FSI results for this configuration. Only for 
$f_\theta /f_h=1$ can a stable region be found from FSI simulation and this lies significantly outside the region for which prescribed motion simulations have been performed.
Pitch energy and thrust maps generated for 
$St =0.2$ and (a,b) 
$X^*=0.25$, (c,d) 
$X^*=0.5$, with 
$\theta _0$ varied from 
$2.5^{\circ }$ to 
$25^{\circ }$ and 
$\psi$ varied from 
$-180^{\circ }$ to 
$0^{\circ }$. The dashed line in the energy and thrust map indicates the equilibrium curve. In the thrust map, the purple circles denote the boundary between thrust and drag and the dashed-dot curve indicates a region within 5 % of the maximum thrust value.
The energy maps for the 
$X_e^* = 0.5$ case also show a nearly vertical zero pitching energy curve. The FSI simulations show that, while 
$f_\theta /f_h=1$ does not reach a stationary state, 
$f_\theta /f_h=2$ and 3 are able to reach stationary states but with high amplitudes of oscillation. Thus, for both 
$X_e^*=0.25$ and 0.5, the nearly vertical zero pitching energy curve leads to either non-sinusoidal behaviour or high amplitude pitching oscillation behaviour, which correspond to drag or very low values of thrust. These high amplitude cases fall far outside the range of the maps in figure 10 and are therefore not shown on the maps. We also note that the pitch energy map is also inverted from the other cases in that regions of negative (positive) pitch energy are to the left (right) of the equilibrium curve for this case. This is connected with the fact that the pitch axis for this case (
$X_e^*=0.5$) is aft of the CoP and therefore an upward heave generates a downward pitch of the foil.
Figure 12 shows the maps of pitching energy and thrust coefficient for 
$St =0.1$ and 
$X_e^*=0.1$. As compared with 
$St = 0.2$, the maximum thrust location here shifts downward to 
$\theta _0=5^\circ$ and 
$\psi =-100^\circ$, requiring a stiffer spring (
$\,f_\theta /f_h>3$) to achieve this value. The maximum thrust coefficient value decreases by a factor of approximately five from 0.2 to 0.04. Figure 11 shows the maps of pitching energy and thrust coefficient for 
$St =0.3$ and 
$X_e^*=0.1$. For 
$St =0.3$ and 
$f_\theta /f_h=3$, higher divergence from the zero pitching energy curve can be seen due to the strong influence of the natural frequency of the spring 
$f_\theta$ on the flow-induced pitching motion which generates a strong superharmonic, as can be seen in figure 8(c). The maximum thrust coefficient for this case is equal to 0.53 (a 2.7
$\times$ increase over the 
$St=0.2$ case) and the location can be seen to be shifted upwards towards higher pitching amplitude of 
$\theta _0=20^\circ$ and 
$\psi =-60^\circ$.
(a) Pitch energy map and (b) thrust coefficient map generated for 
$St =0.3$ and 
$X^*_e=0.1$ with 
$\theta _0$ varied from 
$2.5^{\circ }$ to 
$25^{\circ }$ and 
$\psi$ varied from 
$-180^{\circ }$ to 
$0^{\circ }$. The dashed line in the energy and thrust map indicates the equilibrium curve. In the thrust map, the purple circles denote the boundary between thrust and drag and the dashed-dot curve indicates the region within 5 % of the maximum thrust value.
(a) Pitch energy map and (b) thrust coefficient map generated for 
$St =0.1$ and 
$X^*_e=0.1$ with 
$\theta _0$ varied from 
$2.5^{\circ }$ to 
$25^{\circ }$ and 
$\psi$ varied from 
$-180^{\circ }$ to 
$0^{\circ }$. The dashed line in the energy and thrust map indicates the equilibrium curve. In the thrust map, the purple circles denote the boundary between thrust and drag and the dashed-dot curve indicates the region within 5 % of the maximum thrust value.
Decomposition of the pressure-induced force in the 
$x$-direction acting on a hydrofoil into individual components using the FPM for several cases with 
$X_e^*=0.1$. Panels show (a) 
$\theta _0=10^\circ$, 
$\psi = -100^\circ$ (
$St=0.2$), (b) 
$\theta _0=10^\circ$, 
$\psi = -180^\circ$ (
$St=0.2$), (c) 
$\theta _0=5^\circ$, 
$\psi = -100^\circ$ (
$St=0.1$), (d) 
$\theta _0=25^\circ$, 
$\psi = -60^\circ$ (
$St=0.1$), (e) 
$\theta _0=15^\circ$, 
$\psi = -90^\circ$ (
$St=0.3$), ( f) 
$\theta _0=20^\circ$, 
$\psi = -60^\circ$ (
$St=0.3$).
3.3. A LEV-based model for thrust estimation of flapping foils
The thrust maps shown in figure 9 show a characteristic topology for all the cases simulated here: first, peak thrust is at a pitching–heaving phase angle in the vicinity of 
$-90^\circ$ and for intermediate pitch angles. Second, regions of negative thrust occur for pitch–heave phase angles in the top left and left quadrants of the thrust maps, i.e. where pitch angles are large and the pitch–heave phase angle significantly deviates from 
$-90^\circ$. Such a universal behaviour in thrust is suggestive of a simple underlying mechanism and scaling law, and we explore this in the current section.
To understand the thrust generation mechanism, we first apply the FPM to flapping-foil flows in the 
$x$-direction. Figure 13 shows FPM applied to some selected cases for each Strouhal number of 0.1, 0.2 and 0.3. The two dominant components of vortex-induced force and the force induced due to the kinematics of the foil have been plotted and also a comparison is made between the total force calculated from the FPM and the force in the 
$x$-direction of the foil indicating the validity of FPM. Figure 13(a) shows the case with 
$\theta _0=10$ and 
$\psi =-100$, which corresponds to the case with the highest 
$x$-direction force of 
$-$0.086 (note that a negative value corresponds to thrust). The average values of 
$F_{VIF}$, 
$F_{KINE}$, over one cycle are computed to be 
$-$0.089 and 0.002, respectively. The average values of 
$F_{VIS}$ and 
$F_{BOUND}$ have been found to be very small for all cases in figure 13 and hence are not shown. Locally, these values can be significant, thereby contributing to 
$F_{total,FPM}$. Thus, FPM shows that it is the vortex-induced force which dominates thrust generation for all the simulations in figure 13. Contour plots of the 
$Q$ value and the vortex-induced force (VIF) at the instance of highest thrust (i.e. 
$t/T_p=7.79$) are shown in figure 14 for the case in figure 13(a) and it is evident that the LEV is responsible for this peak in thrust. The FPM was used to reach a similar conclusion for LEV on the caudal fin of a carangiform swimmer (Seo & Mittal Reference Seo and Mittal2022).
Contours of non-dimensional 
$Q$-value and non-dimensional VIF density at 
$t/T_p=7.79$ for 
$X_e^*=0.1$, 
$\theta _0=10$ and 
$\psi =-100$.
Schematic depicting the key features of the LEV-based model.
Given the dominant role of the LEV, it would be useful to develop a simple model of thrust that could be used for the rapid analysis of such flapping foils. Significant work has been done on the mathematical modelling of leading-edge separation to predict the response of pitching foils (Eldredge & Jones Reference Eldredge and Jones2019; Seo, Raut & Mittal Reference Seo, Raut and Mittal2023). Researchers have also used inviscid flow theory to predict the influence of the LEV on the flow and forces (Brown & Michael Reference Brown and Michael1954; Edwards Reference Edwards1954; Mangler & Smith Reference Mangler and Smith1959; Polhamus Reference Polhamus1966), extending classical thin airfoil theory with some representation of the separated flow. From these models we understand that LEV will depend on effective angle-of-attack and the local velocity. We combine all of these insights with the results from our FPM to derive a new phenomenological model for LEV-based thrust generation of flapping foils.
According to the Kutta–Joukowski theorem, the force perpendicular to the foil 
$F_N$(see figure 15) is given by the following equation:
\begin{equation} F_N = \rho V \varGamma, \end{equation}
here, 
$\varGamma$ is the circulation of the hydrofoil and 
$V$ is the net relative velocity of the foil with respect to the flow given by 
$V=\sqrt {U^2_\infty +\dot {h}^2}$. We next assume, based on the FPM analysis, that the circulation 
$\varGamma$ is dominated by the LEV. The formation and strength of the LEV is proportional to the component of velocity 
$V$ of the flow at the leading edge that is normal to the surface. This component is associated with the effective angle-of-attack (
$\alpha _{eff}$) induced by the flow and the movement of the foil. Thus, 
$\varGamma$ can be expressed as
\begin{equation} \varGamma \propto V \sin \alpha_{eff}, \end{equation}where the effective angle-of-attack from the surface of the foil, is given by
\begin{equation} {\alpha_{eff}} (t) = {\tan ^{ - 1}}({V_{LE}(t)/U_\infty} )- \theta_{LE}(t), \end{equation}
where the first term represents the angle-of-attack induced at the leading edge by the heaving and pitching of the foil with the lateral velocity of the leading edge obtained as 
$V_{LE}(t) = -( \dot h- \dot {\theta }(t) X_e\cos {\theta (t)})$. The second term 
$\theta _{LE}$ represents the inclination of the surface of the airfoil over which the LEV develops and induces it force on and is expressed as 
$\theta _{LE}(t) =\theta (t)- \theta _s \, \textrm {sgn}(\dot {h}(t))$. In this expression, the variable 
$\theta _s$, is a fixed angle offset between this surface and the chord of the airfoil and is related to the thickness and shape of the foil.
Given that the LEV-induced force is determined not just by 
$U_\infty$ but by 
$V$, which accounts for the velocity experienced by the flow at the leading edge of the flapping foil, we hypothesize the following modified normalization of the normal force:
\begin{equation} C_N = \frac{F_N}{\dfrac{1}{2}\rho V_{max}^2 C}, \end{equation}
where we employ the maximum value of 
$V$ instead of 
$U_\infty$. Then, using (3.1), (3.2) and (3.4) we arrive at 
$C_N \propto \sin ( \alpha _{eff} )$. From this, it follows that the thrust coefficient satisfies the following proportionality:
\begin{equation} C_T \propto \sin {( \alpha_{eff} )} \sin{( \theta_{LE} )}. \end{equation}
Denoting the time-dependent function on the right-hand side, which is written solely in terms of the kinematics and geometry of the foil by 
$\kappa _{LEV}$, the model for LEV-induced thrust can be expressed as
\begin{equation} C_T (t) = A \kappa_{LEV}(t;\theta_s), \end{equation}
where 
$\kappa _{LEV}(t;\theta _s)$ is given and can be expanded (using (3.3)) as
\begin{align} \kappa_{LEV}(t;\theta_s) &=\sin {( \alpha_{eff} )} \sin{( \theta_{LE} )} \end{align}
\begin{align} &=\sin {( {\tan ^{ - 1}}({V_{LE}(t)/U_\infty} )- \theta_{LE}(t) )} \sin{( \theta_{LE} )}. \end{align}
For a given kinematics of heaving and pitching, (3.6) has two unknown parameters: 
$A$, the constant of proportionality, and 
$\theta _s$, which is connected with the shape of airfoil. In order to determine these two parameters we apply the above model to the time-averaged value of the thrust coefficient i.e. 
$\bar {C}_T = A \bar {\kappa }_{LEV}(\theta _s)$, and the true thrust coefficient is estimated from our simulation data.
The values of 
$A$ and 
$\theta _s$ for the LEV-based model (LEVBM) are determined by using a least-square error fit for the 462 prescribed motion simulations that have been conducted for various pitch angle amplitudes, pitch-axis locations and Strouhal numbers. The fit results in 
$A=3.41$ and 
$\theta _s = 3.62^\circ$ and figure 16 show the correlation between the mean thrust coefficient 
$\bar {C}_T$ and the mean 
$\bar {\kappa }_{LEV}$ using these values of 
$A$ and 
$\theta _s$. The fit has a 
$R^2$ value of 0.91 indicating a relatively high degree of correlation between the data and the model. We note from figure 16 that there is a larger scatter in the data at the two ends of the data set. The lower end corresponds to cases with extreme 
$\psi$ values (close to 0 or 
$-$180) where LEV formation is suppressed, thereby generating drag. Hence, the LEVBM has a more difficult time in predicting the force for these cases. The scatter on the upper end is primarily associated with the high Strouhal number case and this will be discussed in the next section.
Plot showing the thrust coefficient (
$\bar {C}_T$) and thrust factor (
$\bar {\kappa }_{LEV}$) from each of 462 simulations and a linear fit through the data points.
One can rewrite the expression of 
$\kappa _{{LEV}}$ from (3.8) as
\begin{equation} \kappa_{{LEV}} = \frac{{ V_{LE}\sin (\theta_{LE}) \cos (\theta_{LE}) - U_\infty{{\sin }^2}(\theta_{LE}) }}{{\sqrt {{{V_{LE}}^2} + U^2_\infty} }}, \end{equation}
and it is clear from (3.5), (3.6) and (3.9) that, to maximize 
$\kappa _{{LEV}}$ , 
$-\dot {h}$ has to be in phase with 
$\theta$. For the sinusoidal pitching and heaving, this can be achieved with the pitching–heaving phase difference of 
$\psi =-{\rm \pi} /2$ and therefore the peak in the LEVBM predicted thrust will always be at 
$\psi =-{\rm \pi} /2$. This is also clear from the plot of 
$\bar {\kappa }_{LEV}$ in figure 17 for different values of 
$\theta _0$ and 
$\psi$ showing maxima at 
$\psi =-{\rm \pi} /2$. The influence of this on model prediction is also discussed in the next section.
Plot of 
$\bar {\kappa }_{LEV}$ for different values of 
$\psi$ and 
$\theta _0$.
3.4. Model comparison and prediction
Figure 18 shows the thrust coefficient map obtained from the above described model which is based on the LEV. The first figure shows the model prediction for the intermediate Strouhal number case of 
$St = 0.2$ and we note that, qualitatively, the model predicts the topology of the thrust map in figure 9(b) with reasonable accuracy. The peak thrust for both the direct numerical simulation (DNS) and the LEVBM is located in the region centred around 
$\theta _o$ of 
$10^\circ$ and 
$12^\circ$ and very near 
$\psi = -90^\circ$. The margin between thrust and drag generation is also predicted reasonably well. For the LEV model, we have also extended the thrust map to the region of 
$\psi$ between 
$0^\circ$ and 
$180^\circ$ and the model predicts that large drag would be generated for 
$\psi = 90^\circ$ and large values of 
$\theta _o$.
Thrust coefficient maps generated from the model for 
$X_e^*=0.1$ and (a) 
$St = 0.2$, (b) 
$St = 0.1$ and (c) 
$St = 0.3$. The zero thrust curve is shown by the dashed line and the dotted lines indicate the error bounds in the zero thrust curve predicted by the LEVBM. The empty circles represent the zero thrust curve from the DNS.
Figure 18(b) shows the thrust map obtained for 
$St = 0.1$ using LEVBM. When compared with the thrust map in figure 12(b), the model prediction is qualitatively reasonable for this case as well. For instance, the peak thrust now occurs at a lower 
$\theta _o$ and this is captured by the model. However, the comparison is clearly not as good as that for 
$St = 0.2$. This can be explained by noting that the overall magnitude of the thrust for this case is reduced compared with 
$St = 0.2$ by a factor of nearly fivefold. The LEV for this case is therefore not as dominant as for the higher Strouhal number case and thus the error inherent in the regression is much more apparent for this case. The boundary between thrust and drag has a complicated shape in the DNS generated map and, while this characteristic is not captured in the LEVBM, the DNS predicted boundary does fall within the error bounds of the LEVBM predicted boundary.
In figure 18(c), we plot the data for the highest Strouhal number (
$St = 0.3$) case. As for the 
$St = 0.3$, the boundary between drag and thrust is predicted very well by the LEVBM for this case. However, there are several features of the DNS thrust map in figure 11(b) that are not captured as well by the LEVBM. First, there are two values of 
$\theta _o$ where large values of thrust are obtained; one around 
$\theta _o=10^\circ$ and the other around 
$\theta _o=25^\circ$. The LEVBM model does not predict these two values but does predict the upward shift in the 
$\theta _o$ for peak thrust to approximately 
$18^\circ$, which is roughly the mean between the two peak values of the DNS map. Second, the DNS map indicates that peak thrust is obtained for 
$\psi \approx -60^\circ$, whereas the peak for all the LEVBM generated maps is precisely at 
$\psi =-90^\circ$, as discussed earlier. The fact that the peak in the DNS data is not exactly at 
$-{\rm \pi} /2$ for any of the cases and is particularly distant from 
$-{\rm \pi} /2$ for the 
$St = 0.3$ case is due to the fact that, at this high Strouhal number, the effective angle-of-attack is relatively high and this results in the formation and shedding of multiple LEVs in every cycle. This nonlinear effect cannot be accounted for in the LEVBM. This also explains the relatively large scatter in the data in figure 16 as well as a noticeable deviation from the model at the high end of the 
$\kappa _{{LEV}}$, which is all due to the 
$St = 0.3$ case.
Given that the model is constructed so as to predict the maximum thrust at 
$\psi =-{\rm \pi} /2$, we compare the performance of the model against the DNS data along 
$\psi =-{\rm \pi} /2$ for all three Strouhal numbers in figure 19. We note that, along this manifold of the parameter range, the model provides a reasonably good prediction of the qualitative and quantitative features of thrust. In particular, the regions of thrust and drag generation match well with the DNS data and so does the trend in thrust increase with increase in Strouhal number. The pitch angle corresponding to maximum thrust for a given Strouhal number along 
$\psi =-{\rm \pi} /2$ is also predicted reasonably well.
Plot of 
$\bar {C}_T$ at 
$\psi =-90^\circ$ and different 
$\theta _0$ and 
$St$ from the (a) DNS simulations and the (b) LEVBM. The dotted line indicates the maximum value of 
$\bar {C}_T$ for that particular 
$St$.
While the previous comparison between the DNS and the LEVBM over the entire parameter space such as in figure 18 provides useful insights into the predictive capability of the model, solutions for the full FSI system exist mostly (exceptions are when the response is highly non-sinusoidal such as for the case in figure 7a) over the manifold corresponding to zero pitch energy in the energy map. Thus, the performance of the LEVBM model over this manifold is particularly relevant for practical considerations of WAP systems. Figure 20 shows the comparison of mean thrust coefficient along this manifold obtained from the DNS with the corresponding thrust predicted by the LEVBM. The DNS prediction for both 
$St =0.1$ and 0.2 falls within the uncertainty bounds of the regression fit for the model. For 
$St=0.3$, while the trend is predicted well, some of the DNS data fall marginally outside the uncertainty bounds of the regression fit. Overall, the comparison indicates that the model performs reasonably well in predicting the mean thrust along the zero-pitch energy manifold for all three cases. This particular comparison therefore suggests that the model can serve as a simple but effective tool for the design and analysis of WAP systems. We explore this idea further in the final section below
Comparison of 
$\bar {C}_T$ along the zero pitching energy curve between the DNS simulations and the LEVBM for 
$X_e^*=0.1$ and (a) 
$St =0.1$, (b) 
$St =0.2$ and (c) 
$St =0.3$.
3.5. Model prediction for WAP design
It is well established that the most important parameter for the scaling of the performance of flapping foils is the Strouhal number based on the wake width (Streitlien & Triantafyllou Reference Streitlien and Triantafyllou1998; Schouveiler, Hover & Triantafyllou Reference Schouveiler, Hover and Triantafyllou2005), which for the WAP foil can be estimated as 
$St_w = 2f h_o/U_\infty$. The wave encounter frequency 
$f$ is given by 
$\lambda /U_\infty$, where 
$\lambda$ is the wavelength of the waves, and the Strouhal number can therefore be expressed as 
$St_w = 2h_o/\lambda _w$, indicating that this important parameter is determined by the wave conditions. A WAP powered vehicle will likely encounter a range of wave conditions (and therefore 
$St_w$) and one objective of the design of a WAP system could be to maximize the thrust for the given wave conditions. Here, we examine how the LEVBM predictions could be used for this purpose.
Plot of 
$\theta ^{max}_{0}$ at different 
$St_w$ for 
$X_e^*=0.1$ using the model for thrust estimation (dotted line). Square symbols show 
$\theta ^{max}_{0}$ along the zero pitching energy curve and the cross symbols show the overall 
$\theta ^{max}_{0}$ from DNS simulations for different values of 
$St_w$.
From (3.6), it is apparent that one can maximize the thrust by maximizing 
$\kappa _{LEV}$. In terms of kinematic parameters, 
$\kappa _{LEV}$ can be expressed as
\begin{equation} \kappa_{LEV} = \frac{ {St({\dot{h}^*_0}+\dot{\theta}_0X^*_e\cos(\theta_0))\sin {(\theta_0+\theta_s)}\cos {(\theta _0+\theta_s)} - {{\sin }^2}{(\theta _0+\theta_s)}} } {{\sqrt {St^2({\dot{h}^*_0}+\dot{\theta_0}X^*_e\cos(\theta_0))^2 + 1} }}. \end{equation}
From the above equation, for the condition 
$\dot {h}_0^*\gg \dot {\theta }_0 X_e^*$, which is valid as long as pitching amplitudes are not very large and/or the pitch axis being close to the leading edge, the pitching amplitude (
$\theta _0$) which maximizes 
$\kappa _{LEV}$ can be found as
\begin{equation} \theta^{{max}}_{0} = {\tan ^{ - 1}}\left( {\frac{{\sqrt {{{({\rm \pi} \,St_w)}^2} + 1} - 1}}{{{\rm \pi} \,St_w}}} \right) - \theta_s. \end{equation}
Equation (3.11) indicates that, firstly, for any wave conditions, the pitch amplitude that generates the maximum thrust is determined primarily by 
$St_w$ and, second, the dependence of this maximal condition is provided by the LEVBM. Indeed, if we ignore 
$\theta _s$, (3.11) is a universal law for selecting a pitch amplitude for a pitching–heaving foil undergoing ‘normal’ flapping (i.e. where the maximum pitch angle occurs at mid-heave) that maximizes thrust for a given wake Strouhal number.
We examine this prediction of the LEVBM model in figure 21(a), which shows the curve predicted by the LEVBM (3.11) and the corresponding values obtained from the three different Strouhal number cases simulated in the DNS for 
$X^*_e=0.1$. For each Strouhal number, we include the DNS predicted values of 
$\theta _o$ corresponding to the global maximum of thrust, the maximum thrust along the zero-pitch energy manifold (along which equilibrium solutions of the FSI system exist) and the 
$\psi = -{\rm \pi} /2$ manifold (along which the model predicts the maximum value). The prediction of the LEVBM is overall a reasonably match for the DNS data, thereby providing support for the use of this scaling law for WAP design. For instance, the scaling law in (3.11) could be used in a feed-forward manner to maximize thrust as follows: the wave height and encounter frequency could be determined via inertial or other sensors and the WAP system could designed so as to dynamically adjust the maximum pitch angle 
$\theta _o$ to a value based on (3.11) that maximizes thrust. This feed-forward control could be achieved in several different ways such as through dynamic adjustment of the spring stiffness or through the use of an adjustable pitch-limiting mechanism.
4. Conclusions
Motivated by application to wave-assisted hydrofoil propulsion, a systematic study has been carried out for sinusoidal heaving and flow-induced pitching motion of hydrofoil at a chord-based Reynolds number of 10 000 using time-accurate simulations of the fluid dynamics coupled with a linear structural model for deflection in pitch. The objectives of this work were twofold – to examine key parameters affecting the dynamics of pitch oscillations such as pitch-axis location and torsional spring stiffness, and to demonstrate the use of the energy exchange between the fluid and the structure as a way to analyse the pitching behaviour. The simulations indicate that placing the pitch axis very close to the leading edge is conducive to thrust generation over a wide range of conditions. Pitch-axis locations further down the chord result either in chaotic oscillations or oscillations where the pitching leads heaving. Both of these conditions are found to diminish thrust.
Energy maps have been utilized to explain the complex pitching behaviour of the hydrofoil. Energy maps are generated by computing the rate of energy exchange between the hydrofoil and the flow, for foils undergoing prescribed heaving and pitching oscillations over a range of pitching amplitude and phase differences between heaving and pitching. For the FSI cases with sinusoidal (or nearly sinusoidal pitching), the energy maps accurately predict the amplitude and phase of the flow-induced pitching. Thus the energy maps can serve as an effective tool for examining the performance and dynamic stability of these WAP foils.
The application of the FPM indicates that the LEV is the dominant feature responsible for thrust generation. Based on this observation, we develop a phenomenological LEVBM for thrust prediction which is parameterized based on over 400 distinct simulations. The predictions of the LEVBM are found to give a reasonable match to the DNS results, and a case is made that the LEVBM could be helpful in guiding the design of WAP systems. Considering the fact that the high-fidelity simulations can only be employed for a small sampling of the very large parameter space, the LEVBM provides a rapid and effective way of exploring the design and performance of WAP systems over a much wider range of the parameter space.
Acknowledgements
This work is supported by ONR Grants N00014-22-1-2655 and N00014-22-1-2770 monitored by Dr R. Brizzolara. This work used the computational resources at the Advanced Research Computing at Hopkins (ARCH) core facility (rockfish.jhu.edu), which is supported by the AFOSR DURIP Grant FA9550-21-1-0303, and the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1548562, through allocation number TG-CTS100002.
Declaration of interests
The authors report no conflict of interest.
Appendix. Grid and domain independence
Figure 22 shows data on grid convergence and domain independence for one of the high amplitude cases of 
$\theta _\circ = 25^\circ$ and 
$\psi =0^\circ$, which provides a robust test of this convergence. For the grid convergence study, a grid with 2X refinement i.e. a 
$1250 \times 1200$ grid, is employed while keeping the domain size the same. For the domain independence study, the domain size has been increased by a factor of two in all directions, while maintaining the same grid resolution. Tables 2 and 3 shows the per cent error in mean drag 
$F_{x,mean}$, root-mean-square (r.m.s.) lift force 
$F_{y,rms}$ and r.m.s. pitch moment 
$M_{z,rms}$ for all these cases. The maximum variation in these quantities is limited to approximately 3 %, which indicates that the simulation results are grid and domain independent. We therefore employ 
$640 \times 600$ grid on a 
$10C \times 10C$ domain for all the simulations in the study.
(a) Grid independence and (b) domain independence study for 
$X^*=0.1$, 
$\theta _\circ = 25^\circ$ and 
$\psi = 0^\circ$ with grid size of 
$1250 \times 1200$ used in the grid independence study and domain size of 
$20C \times 20C$ in the domain independence study.
Details of the domain dependence test with a domain size twice that of the baseline and the same grid arrangement.
Details of the grid dependence test with the domain size as shown in figure 1(b).
Open access
