4.1. Cruise propulsion performance

This study focuses on the effects of the dimensionless submergence depth $H/L$ and dimensionless wavelength $\lambda /L$ on hydrofoil propulsion performance. A total of eight dimensionless submergence depths and nine dimensionless wavelengths are considered, with the parameter settings summarised in table 1. In the simulations, the undulation frequency is set to $f = 2\,\textrm{Hz}$ , and the trailing-edge amplitude is set to $0.1L$ . This combination corresponds to a Strouhal number of $St=0.4$ , which lies within the range $0.2 \leq St \leq 0.4$ commonly associated with high propulsive efficiency in aquatic locomotion (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998). The Reynolds number is fixed at $\textit{Re}=4000$ , which lies within the representative range of $10^3$ – $10^4$ for the present study. Previous studies have shown that, within this Reynolds number range, the large-scale wake topology remains essentially unchanged, while the hydrodynamic forces exhibit similar trends in response to variations in kinematic parameters (Anderson et al. Reference Anderson, Streitlien, Barrett and Triantafyllou1998; Zhang et al. Reference Zhang, Pan, Chao and Zhang2018). Therefore, $\textit{Re}=4000$ is sufficient to capture the dominant flow physics and investigate the underlying propulsion mechanisms.

Table 1.

The specific values of the governing parameters applied in current simulations.

Table 1. Long description

Panel A: Table. The table lists parameters and their specifications. Parameters include Reynolds number, dimensionless submergence depth, dimensionless wavelength, undulation frequency, trailing-edge amplitude. Specifications are provided for each parameter, such as Reynolds number set to 4000, dimensionless submergence depths ranging from 0.4 to 2.5, and dimensionless wavelengths ranging from 0.5 to infinity. The undulation frequency is set to 2 Hz, and the trailing-edge amplitude is set to 0.1.

The dimensionless lift coefficient acting on the hydrofoil is defined as follows:

(4.1) \begin{align} C_L=\frac {F_L}{0.5\rho {U_{ref}}^2L}. \end{align}

The variables $F_T$ and $F_L$ denote the resultant forces acting on the hydrofoil in the forward swimming and lateral directions, respectively. Following Liu et al. (Reference Liu, Yu and Tong2012), the swimming efficiency is quantified using the energy-utilisation ratio $\eta$ . This metric is adopted because the present study focuses on self-propelled cruising. Under these conditions, the cycle-averaged net thrust approaches zero because thrust and drag are balanced, which makes the conventional Froude efficiency unsuitable. The physical meaning of the energy-utilisation ratio is defined as the ratio of the swimming distance $l$ travelled by the hydrofoil over one tail-beat cycle to the total energy consumption per unit mass $W_t/m$ , which can be expressed as follows:

(4.2) \begin{align} \eta =\frac {l}{W_t/m}=\frac {U_c}{P_t/m}. \end{align}

The dimensionless energy-utilisation ratio is defined as $\tilde {\eta } = L\eta /T^2$ . The total power consumption $P_t$ is defined as follows:

(4.3) \begin{align} P_t=\frac {W_t}{T}=\frac {1}{T}\int _{t_0}^{t_0+T}\oint _S-(\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{u_b}\,\textrm{d}S\,\textrm{d}t. \end{align}

Here, $S$ denotes the surface of the hydrofoil, $\boldsymbol{\sigma }$ represents the surface stress tensor acting on the hydrofoil, $\boldsymbol{n}$ denotes the outward unit normal vector on the body surface and $\boldsymbol{u_b}$ denotes the surface velocity of the hydrofoil. The input power coefficient $C_P$ is defined as follows:

(4.4) \begin{align} C_P=\frac {P_t}{0.5\rho U_{ref}^3L}. \end{align}

All cycle-averaged quantities are evaluated after the initial acceleration stage, once the self-propelled motion has reached a statistically steady state, using data from four subsequent undulation cycles. Figure 8(a) shows the variation of the cycle-averaged cruising speed $\overline {U_c}$ with submergence depth $H/L$ for four selected dimensionless wavelengths $\lambda /L$ . Overall, the cruising speed $\overline {U_c}$ increases monotonically with increasing submergence depth $H/L$ for all wavelengths considered. For $\lambda /L = 0.5$ , 1 and 2, the increase in $\overline {U_c}$ is more pronounced within $0.4 \lt H/L \lt 0.6$ than over the broader range $0.6 \lt H/L \lt 2.5$ . When $\lambda /L = \infty$ , the hydrofoil deformation can be regarded as pure pitching motion, and within $0.4 \lt H/L \lt 1$ , the increase in $\overline {U_c}$ with $H/L$ becomes even more pronounced. For $\lambda /L = 0.5$ , 1 and 2, at a given submergence depth $H/L$ , a larger wavelength consistently results in a higher cruising speed $\overline {U_c}$ . Interestingly, an exception is observed for $\lambda /L = \infty$ . Specifically, at $H/L = 0.6$ , the cruising speed $\overline {U_c}$ for $\lambda /L = \infty$ lies between those for $\lambda /L = 1$ and $\lambda /L = 2$ , whereas at $H/L = 0.4$ , it falls between $\lambda /L = 0.5$ and $\lambda /L = 1$ . These observations indicate that the presence of the free surface is generally detrimental to hydrofoil propulsion speed, and that its weakening effect becomes more pronounced with increasing body-wave wavelength. As the submergence depth increases, the inhibiting influence of the free surface on cruising speed is gradually reduced. This trend is qualitatively consistent with the experimental observations of Huera-Huarte (Reference Huera-Huarte2023), who reported that the thrust generated by a pitching foil decreased as it approached the free surface.

To further investigate the influence of body-wave wavelength on the cruising speed $\overline {U_c}$ , four representative dimensionless submergence depths $H/L$ are selected, and the variation of $\overline {U_c}$ with $\lambda /L$ is examined. The results are presented in figure 8(b), where open symbols denote $\lambda /L = \infty$ . It can be observed that, for a given wavelength $\lambda /L$ , a larger submergence depth $H/L$ consistently corresponds to a higher cruising speed $\overline {U_c}$ . At $H/L = 2.5$ , $\overline {U_c}$ increases monotonically with $\lambda /L$ . When $\lambda /L \lt 2$ , $\overline {U_c}$ is relatively sensitive to changes in $\lambda /L$ , whereas for $\lambda /L \gt 2$ , the influence of further increases becomes weaker. The cruising speed $\overline {U_c}$ for $\lambda /L = \infty$ is slightly higher than that for $\lambda /L = 8$ , but the two values are very close, indicating convergence of $\overline {U_c}$ with increasing $\lambda /L$ . These trends are consistent with observations reported by Thekkethil et al. (Reference Thekkethil, Sharma and Agrawal2018) under unbounded flow conditions. When the hydrofoil operates closer to the free surface, the trend changes markedly. For $H/L = 0.8$ , $\overline {U_c}$ increases with $\lambda /L$ over $0.5 \lt \lambda /L \lt 5$ , but decreases as $\lambda /L$ increases further, and the cruising speed for $\lambda /L = \infty$ becomes lower than that for $\lambda /L = 8$ . Similar non-monotonic trends are observed for $H/L = 0.6$ and $0.4$ , where $\overline {U_c}$ first increases and then decreases with $\lambda /L$ . For $H/L = 0.6$ , values of $\lambda /L$ between 3 and 5 yield the highest cruising speeds, whereas for $H/L = 0.4$ , optimal performance is achieved at $\lambda /L = 2$ . These results indicate that, under near-free-surface conditions, an optimal wavelength $\lambda ^*/L$ exists at each submergence depth $H/L$ to ensure high-speed propulsion. Moreover, as the submergence depth $H/L$ decreases, the optimal wavelength corresponding to the maximum cruising speed exhibits a decreasing trend. In addition, for $\lambda /L \lt 2$ , the rate at which $\overline {U_c}$ increases with $\lambda /L$ gradually diminishes as $H/L$ decreases, suggesting that the free surface suppresses the speed gains associated with increasing wavelength.

Figure 8.

(a) Variation of the time-averaged cruising speed $\overline {U_c}$ with $H/L$ at selected values of $\lambda /L$ ; (b) variation of the time-averaged cruising speed $\overline {U_c}$ with $\lambda /L$ at selected values of $H/L$ (open symbols denote $\lambda = \infty$ ).

Figure 8. Long description

Panel A: A line graph shows the variation of the time-averaged cruising speed with H/L at selected values of λ/L. The x-axis represents H/L, and the y-axis represents the time-averaged cruising speed normalized by a reference speed. Different symbols and colors represent different values of λ/L. Panel B: Another line graph shows the variation of the time-averaged cruising speed with λ/L at selected values of H/L. The x-axis represents λ/L, and the y-axis represents the time-averaged cruising speed normalized by a reference speed. Different symbols and colors represent different values of H/L. Open symbols denote specific conditions.

Figure 9(a) presents the variation of the cycle-averaged input power coefficient $\overline {C_P}$ with submergence depth $H/L$ for four dimensionless wavelengths $\lambda /L$ . Overall, at a given submergence depth, a larger wavelength consistently results in a higher $\overline {C_P}$ , while $\overline {C_P}$ decreases monotonically as the submergence depth decreases. This trend suggests that swimming near the free surface has the potential to improve propulsion efficiency, consistent with previous results obtained using tethered models (Zheng et al. Reference Zheng, He, Zhao, Shen and Zhu2024; Gao et al. Reference Gao, Sun, Zhang and Liu2025). In viscous fluid dynamics, added mass arises from the additional kinetic energy imparted to the surrounding fluid by a moving solid body. For a sphere moving near a free surface, the added mass can be expressed as $M_H=(2/3)K_H\pi r^3\rho _w$ , where the coefficient $K_H$ depends on the dimensionless submergence depth $H/(2r)$ . As $H/(2r)$ decreases, $K_H$ decreases, indicating that the added mass is reduced when the body approaches the free surface (Brennen Reference Brennen1982). Although the present problem involves prescribed body deformation, unsteady vortex shedding and nonlinear free-surface motion, for which the classical added-mass theory is not expected to provide a quantitative description, the above analogy offers a useful qualitative interpretation of the observed trend. This mechanism can partially explain the reduction in $\overline {C_P}$ observed for hydrofoils undulating near the free surface. Because the pressure at the free surface is zero and the surface is deformable, fluid motion normal to the surface becomes easier as the hydrofoil approaches it, which is qualitatively consistent with a reduction in the corresponding added mass. The classical added-mass analogy suggests that a smaller added mass reduces the inertial work required to accelerate the surrounding fluid during undulatory motion, thereby providing a possible physical interpretation for the observed reduction in the input power.

To further investigate the influence of $\lambda /L$ on $\overline {C_P}$ , the variation of $\overline {C_P}$ with $\lambda /L$ at four submergence depths $H/L$ is shown in figure 9(b), where open symbols denote $\lambda /L = \infty$ . Overall, for all values of $\lambda /L$ , a deeper submergence depth $H/L$ consistently leads to a higher cycle-averaged input power coefficient $\overline {C_P}$ . In addition, at a given submergence depth, $\overline {C_P}$ increases monotonically with $\lambda /L$ . For all submergence depths considered, $\overline {C_P}$ for $\lambda /L = \infty$ is significantly higher than that for $\lambda /L = 8$ . Under near-free-surface conditions, the effect of $\lambda /L$ on $\overline {C_P}$ is consistent with that observed in unbounded flows (Thekkethil et al. Reference Thekkethil, Sharma and Agrawal2018). As $\lambda /L$ increases, the active undulatory motion progressively approaches rigid-body pitching. This transition results in a larger effective contact area and stronger interaction between the hydrofoil and the surrounding fluid during flapping, thereby increasing the required input power.

Figure 9.

(a) Variation of the time-averaged input power coefficient $\overline {C_P}$ with $H/L$ at selected values of $\lambda /L$ ; (b) variation of the time-averaged input power coefficient $\overline {C_P}$ with $\lambda /L$ at selected values of $H/L$ (open symbols denote $\lambda = \infty$ ).

Figure 10(a) shows the variation of the efficiency $\tilde {\eta }$ with submergence depth $H/L$ for four wavelengths $\lambda /L$ . Overall, at a given submergence depth $H/L$ , higher efficiency $\tilde {\eta }$ consistently corresponds to a smaller wavelength $\lambda /L$ . Notably, for all wavelengths considered, the efficiency $\tilde {\eta }$ exhibits a non-monotonic dependence on submergence depth, first increasing and then decreasing, indicating the existence of an optimal submergence depth $H^*/L$ that maximises propulsion efficiency. As $\lambda /L$ increases, the optimal submergence depth $H^*/L$ increases. Specifically, $H^*/L \approx 0.6$ for $\lambda /L \approx 0.5$ , $H^*/L \approx 0.8$ for $\lambda /L \approx 1$ , and $H^*/L \approx 1$ for both $\lambda /L = 2$ and $\lambda /L = \infty$ . These results indicate that an undulating hydrofoil can achieve efficiency gains when cruising at an appropriate distance below the free surface, compared with unbounded-flow conditions. Since the submergence depths are sampled discretely, with a spacing of $0.2L$ in the shallow-submergence range and coarser spacing at larger depths, the exact optimal submergence depth cannot be directly obtained from the computed data. Therefore, for each wavelength, the propulsive-efficiency curve is fitted with a local quadratic polynomial, from which the corresponding optimal submergence depth is estimated. In addition to the four wavelength cases discussed above, five additional cases ( $\lambda /L=0.75$ , 1.5, 3, 5 and 8) are analysed in the same manner. The resulting nine estimated optimal submergence depths are then fitted using an empirical exponential function with saturation characteristics. The functional form is given as follows:

(4.5) \begin{align} H^*/L = 1.09 - (1 - a) \textrm{e}^{-b \lambda /L }. \end{align}

Here, $a$ and $b$ are fitting parameters. The fitting procedure yields $a = 0.310$ and $b = 0.718$ . Figure 11 illustrates the variation of the optimal submergence depth $H^*/L$ with increasing $\lambda /L$ . It can be observed that when $\lambda /L \lt 2$ , $H^*/L$ is relatively sensitive to variations in $\lambda /L$ , whereas for $\lambda /L \gt 2$ , $H^*/L$ gradually approaches $1.1$ . The observed dependence on submergence depth also has clear biological relevance. Recent animal-tracking studies have shown that non-foraging marine reptiles, birds and mammals tend to travel at fixed subsurface depths that reduce wave drag and the associated cost of transport (Stokes et al. Reference Stokes, Esteban, Casale, Chiaradia, Kaska, Kato, Luschi, Ropert-Coudert, Stokes and Hays2024). This optimum occurs neither immediately beneath the free surface, where wave-induced energy losses are pronounced, nor under deep-water conditions. This behaviour is qualitatively consistent with the present results, which reveal an optimal cruising depth at which propulsive efficiency is maximised. The present findings therefore provide a hydrodynamic explanation for why marine animals may select an intermediate submergence depth during sustained travel.

Figure 10(b) shows the variation of the efficiency $\tilde {\eta }$ with wavelength $\lambda /L$ at five submergence depths $H/L$ . Overall, variations in $\lambda /L$ have a more pronounced effect on the efficiency $\tilde {\eta }$ than changes in $H/L$ . For a fixed $\lambda /L$ , comparison between submergence depths indicates that when $\lambda /L = 0.5$ , the case with $H/L = 0.6$ yields the highest efficiency $\tilde {\eta }$ , whereas the efficiency at $H/L = 2.5$ is significantly lower than that at the other four depths. However, when $\lambda /L$ increases to 8 and beyond, the efficiency $\tilde {\eta }$ at $H/L = 2.5$ consistently exceeds that at the other four submergence depths. This result further indicates that the optimal submergence depth for maximum cruising efficiency depends strongly on the body-wave wavelength. Hydrofoils with shorter wavelengths tend to achieve their maximum efficiency at shallower submergence depths, whereas longer wavelengths generally require larger submergence depths to attain optimal efficiency. At all submergence depths considered, $\tilde {\eta }$ decreases monotonically with increasing $\lambda /L$ . It is worth noting that when $\lambda /L \lt 1$ , decreasing $\lambda /L$ at $H/L = 2.5$ results in a slower increase in $\tilde {\eta }$ compared with the other submergence depths. Although decreasing the wavelength generally enhances cruising efficiency, beyond a certain submergence depth and sufficiently small wavelength, further reductions yield diminishing efficiency gains. This behaviour is consistent with observations reported for unbounded flows (Thekkethil et al. Reference Thekkethil, Sharma and Agrawal2018). In practice, the efficiency gain associated with decreasing wavelength is bounded; once the wavelength is reduced beyond a certain threshold, the efficiency reaches a maximum. The additional result at $H/L=1.5$ further indicates that the diminishing efficiency gain associated with decreasing wavelength becomes evident at approximately this submergence depth. As highlighted in the magnified view, although reducing $\lambda /L$ continues to improve the propulsion efficiency, the corresponding efficiency gain becomes noticeably smaller than that at shallower submergence depths.

Figure 10.

(a) Variation of the efficiency $\tilde {\eta }$ with $H/L$ for selected values of $\lambda /L$ ; vertical dotted lines indicate the optimal submergence depths. (b) Variation of $\tilde {\eta }$ with $\lambda /L$ for selected values of $H/L$ ; open symbols denote the case with $\lambda /L=\infty$ .

Figure 10. Long description

Panel A: A line graph shows the variation of efficiency with submergence depths for selected values of a dimensionless parameter. The x-axis represents the ratio of submergence depth to length (H/L), and the y-axis represents the efficiency (η). Vertical dotted lines indicate optimal submergence depths. Different markers and line styles represent different values of the dimensionless parameter (λ/L). Panel B: Another line graph shows the variation of a parameter with the dimensionless parameter for selected values of the submergence depth ratio. The x-axis represents the dimensionless parameter (λ/L), and the y-axis represents the parameter (η). Open symbols denote the case with infinite submergence depth ratio (H/L = ∞). Different markers and line styles represent different values of the submergence depth ratio (H/L).

Figure 11.

Evolution of the optimal submergence depth $H^*/L$ with the wavelength $\lambda /L$ .

4.2. Cruise stability analysis

Figure 12(a) shows the variation of the cycle-averaged lift coefficient $\overline {C_L}$ with submergence depth $H/L$ for four wavelengths $\lambda /L$ . It is observed that when the hydrofoil cruises close to the free surface, it is consistently subjected to an attractive force directed towards the surface. For all wavelengths considered, this attractive force decreases with increasing submergence depth. In essence, the presence of a free surface induces asymmetry in the computational domain, which gradually diminishes as the submergence depth increases, causing the attractive force to approach zero. To further examine the influence of wavelength, figure 12(b) presents the variation of the cycle-averaged lift coefficient $\overline {C_L}$ with wavelength $\lambda /L$ at four submergence depths $H/L$ . With a denser sampling of wavelength values, it is found that for all four submergence depths, $\overline {C_L}$ increases monotonically with $\lambda /L$ . At $H/L = 0.4$ , where the free-surface-induced flow asymmetry is relatively strong, $\overline {C_L}$ for $\lambda /L = \infty$ is significantly larger than that for $\lambda /L = 8$ . In contrast, for $H/L = 0.6$ , 0.8 and 2.5, $\overline {C_L}$ for $\lambda /L = \infty$ and $\lambda /L = 8$ are nearly identical. This indicates that, at each submergence depth, $\overline {C_L}$ ultimately converges to an approximately constant value with increasing $\lambda /L$ . During cruising, maintaining the hydrofoil at a fixed depth requires a vertical force balance to counteract this attractive force. In general, two strategies can be employed: adjusting the angle of attack or modifying the displaced volume. As the submergence depth $H/L$ decreases and the wavelength $\lambda /L$ increases, a larger negative angle of attack or greater fluid intake into the vessel’s compartment is necessary to maintain cruising stability.

Figure 12.

(a) Variation of the cycle-averaged lift coefficient $\overline {C_L}$ with $H/L$ at selected values of $\lambda /L$ ; (b) variation of the cycle-averaged lift coefficient $\overline {C_L}$ with $\lambda /L$ at selected values of $H/L$ (open symbols denote $\lambda = \infty$ ).

Figure 12. Long description

Panel A: A line graph shows the variation of the cycle-averaged lift coefficient with H/L at selected values of λ/L. The x-axis represents H/L ranging from 0 to 2.5, and the y-axis represents the cycle-averaged lift coefficient ranging from 0 to 0.6. Different symbols and line styles represent different values of λ/L (0.5, 1, 2, and ∞). The graph shows a general trend of decreasing lift coefficient with increasing H/L. Panel B: Another line graph shows the variation of the cycle-averaged lift coefficient with λ/L at selected values of H/L. The x-axis represents λ/L ranging from 0 to 9, and the y-axis represents the cycle-averaged lift coefficient ranging from 0 to 0.7. Different symbols and line styles represent different values of H/L (0.4, 0.6, 0.8, and 2.5). The graph shows a general trend of increasing lift coefficient with increasing λ/L.

Figure 13 shows the temporal evolution of the lift coefficient $C_L$ over one undulation cycle at $\lambda /L = 1$ for four submergence depths, namely $H/L = 0.4$ , 0.6, 0.8 and 2.5. Within a single cycle, $C_L$ exhibits one peak and one trough, with the absolute magnitude of the peak consistently exceeding that of the trough. The undulatory motion of the hydrofoil during the first half of the cycle is also illustrated, revealing that the trough of $C_L$ occurs when the trailing edge flaps towards the free surface. Because the free surface can deform in response to the hydrofoil motion, a pressure-release-like effect is induced, weakening the interaction between the hydrofoil and the surrounding fluid. For different submergence depths $H/L$ , the peak lift coefficient remains nearly unchanged, whereas the primary differences arise from the trough values. As $H/L$ decreases, the trough of $C_L$ increases gradually. This trend is qualitatively consistent with the reduced added mass near the free surface, although unsteady vortex dynamics and free-surface deformation may also contribute. It is this increase in the trough of $C_L$ during the first half of the cycle, as $H/L$ decreases, that causes the hydrofoil to experience an attractive force towards the free surface, which increases monotonically with decreasing $H/L$ .

Figure 13.

Time histories of the lift coefficient $C_L$ over one undulation cycle for four different submergence depths $H/L$ at $\lambda /L = 1$ .

During self-propelled undulatory motion, the net streamwise force acting on the hydrofoil is nearly zero, and the input power is strongly correlated with the lift force. Figure 14 shows the temporal evolution of the input power coefficient $C_P$ over one undulation cycle for $\lambda /L = 1$ at four submergence depths, namely $H/L = 0.4$ , 0.6, 0.8 and 2.5. Overall, consistent with the behaviour of $C_L$ , the input power coefficient $C_P$ exhibits two distinct peaks within a single cycle. During the first half of the cycle, as illustrated in figure 14, the trailing edge of the hydrofoil flaps upward towards the free surface. For all values of $H/L$ , the first peak of $C_P$ is consistently smaller than the second peak, because the free surface deforms in response to the hydrofoil motion and releases pressure in the fluid above the hydrofoil. Among different submergence depths, the most pronounced differences are observed in the first peak of $C_P$ . As $H/L$ decreases, the first peak of $C_P$ decreases monotonically, which is qualitatively consistent with the reduction in the effective added mass near the free surface.

Figure 14.

Time histories of the input power coefficient $C_P$ over one undulation cycle for four different submergence depths $H/L$ at $\lambda /L = 1$ .

4.3. Vortical structures and stream velocities

The wake vortex structures generated by the hydrofoil are closely linked to its self-propelled performance. Figure 15(a–p) presents the vorticity ( $\omega _Z$ ) fields of cruising wakes for a total of 16 cases with different submergence depths $H/L$ and dimensionless wavelengths $\lambda /L$ . Specifically, $\lambda /L$ takes values of 0.5, 1, 2 and $\infty$ , whereas $H/L$ takes values of 0.4, 0.6, 0.8 and 2.5, with all snapshots taken at the same phase instant $t/T = 0.75$ . For the instantaneous flow field results, $t/T=0$ is defined as the beginning of the selected cycle during statistically steady cruising. Figure 15(a–p) reveals a phase difference in the trailing-edge undulation of the hydrofoil between different wavelength cases. Moreover, as shown in figure 15(a–c), at finite body-wave wavelengths, the free-surface deformation closely follows the phase of the hydrofoil undulation, such that the free-surface profile above the hydrofoil broadly conforms to its instantaneous shape. By contrast, this correspondence is no longer apparent in figure 15(d), where the pure-pitching motion at $\lambda /L=\infty$ intensifies the hydrofoil–fluid interaction and induces pronounced wave overturning and breaking, thereby obscuring the geometric similarity between the free-surface profile and the hydrofoil shape.

From figures 15(a,e,i,m), 15(b, f, j,n) and 15(c,g,k,o), it can be observed that for $\lambda /L = 0.5$ , 1 and 2, decreasing the submergence depth $H/L$ causes the hydrofoil wake to transition from a relatively symmetric ‘2S’-type reverse Kármán vortex street to an asymmetric ‘P+S’-type vortex street, consistent with previous tethered-model simulations (Zheng et al. Reference Zheng, He, Zhao, Shen and Zhu2024; Gao et al. Reference Gao, Sun, Zhang and Liu2025). In a ‘2S’-type reverse Kármán vortex street, anticlockwise positive vortices form above the wake centreline and clockwise negative vortices below it, inducing a backward jet that provides thrust for propulsion. In contrast, as illustrated in figure 16, with the formation of a ‘P+S’-type vortex street, the positive vortices in the original ‘2S’ pattern undergo splitting and form a row of weaker positive vortices on the lower side of the wake. These smaller positive vortices interact with nearby negative vortices, inducing a forward jet (blue) that weakens thrust generation and hinders high-speed propulsion. From figure 15(d,h,l,p), it can be observed that for $\lambda /L = \infty$ , as $H/L$ decreases from 2.5 to 0.8, the wake transitions from a ‘2S’ to a ‘P+S’-type vortex street. As shown in figure 15(d,h), with further reductions in $H/L$ , positive vortices undergo more pronounced splitting, generating a larger number of smaller and stronger vortices. These vortices interact with negative vortices and further weaken thrust production. A theoretical analysis of the thrust-reducing mechanism associated with the ‘P+S’-type vortex street is provided in our previous work (Gao et al. Reference Gao, Sun, Zhang and Liu2025). Overall, as the undulating hydrofoil approaches the free surface from deeper submergence, the evolution can be divided into two stages. In the first stage, the hydrofoil remains sufficiently far from the free surface, resulting in only weak surface disturbances. In the second stage, pronounced free-surface deformation and even wave breaking emerge, accompanied by intensified wave activity.

Figure 15.

Vorticity field snapshots corresponding to different submergence depths $H/L$ and wavelengths $\lambda /L$ : (a–d) $H/L = 0.4$ , (e–h) $H/L = 0.6$ , (i–l) $H/L = 0.8$ and (m–p) $H/L = 2.5$ ; (a,e,i,m) $\lambda /L = 0.5$ , (b, f,j,n) $\lambda /L = 1$ , (c,g,k,o) $\lambda /L = 2$ and (d,h,l,p) $\lambda /L = \infty$ . All snapshots are taken at the same phase, $t/T=0.75$ , during statistically steady cruising.

Figure 15. Long description

Panel A: A heat map showing vorticity fields for a submergence depth of 0.4 and wavelength to length ratio of 0.5. The color scale ranges from -8 to 8, with blue indicating negative vorticity and red indicating positive vorticity. Panel B: A heat map showing vorticity fields for a submergence depth of 0.4 and wavelength to length ratio of 1. Panel C: A heat map showing vorticity fields for a submergence depth of 0.4 and wavelength to length ratio of 2. Panel D: A heat map showing vorticity fields for a submergence depth of 0.4 and wavelength to length ratio of infinity. Panel E: A heat map showing vorticity fields for a submergence depth of 0.6 and wavelength to length ratio of 0.5. Panel F: A heat map showing vorticity fields for a submergence depth of 0.6 and wavelength to length ratio of 1. Panel G: A heat map showing vorticity fields for a submergence depth of 0.6 and wavelength to length ratio of 2. Panel H: A heat map showing vorticity fields for a submergence depth of 0.6 and wavelength to length ratio of infinity. Panel I: A heat map showing vorticity fields for a submergence depth of 0.8 and wavelength to length ratio of 0.5. Panel J: A heat map showing vorticity fields for a submergence depth of 0.8 and wavelength to length ratio of 1. Panel K: A heat map showing vorticity fields for a submergence depth of 0.8 and wavelength to length ratio of 2. Panel L: A heat map showing vorticity fields for a submergence depth of 0.8 and wavelength to length ratio of infinity. Panel M: A heat map showing vorticity fields for a submergence depth of 2.5 and wavelength to length ratio of 0.5. Panel N: A heat map showing vorticity fields for a submergence depth of 2.5 and wavelength to length ratio of 1. Panel O: A heat map showing vorticity fields for a submergence depth of 2.5 and wavelength to length ratio of 2. Panel P: A heat map showing vorticity fields for a submergence depth of 2.5 and wavelength to length ratio of infinity.

As shown in figure 15(a–p), with $\lambda /L$ held fixed, the amplitude of free-surface wave generation increases as the submergence depth $H/L$ decreases. As observed in figure 15(a,e,i), for $\lambda /L = 0.5$ , no vorticity is generated near the free surface throughout the simulations. This indicates that the free surface undergoes only mild oscillations without overturning or breaking. As shown in figure 15(b), for $\lambda /L = 1$ , free-surface breaking occurs at $H/L = 0.4$ , accompanied by the generation of positive vorticity (shown in red) near the free surface. As shown in figure 15(g,h), when $\lambda /L = 2$ and $\lambda /L = \infty$ , wave breaking initiates at $H/L = 0.6$ and becomes particularly intense; high-vorticity regions and vortex pairs form near the free surface and are subsequently convected into the wake. As shown in figure 8(a), where the cruising speed $\overline {U_c}$ decreases monotonically with decreasing submergence depth, the combined effects of wake asymmetry and intensified free-surface activity are associated with reduced near-surface cruising speed. In particular, for $\lambda /L = \infty$ , both wake asymmetry and free-surface breaking are significantly intensified as the hydrofoil approaches the free surface. Consequently, as observed in figure 8(a), the reduction in $\overline {U_c}$ with decreasing $H/L$ is markedly more rapid over the range $0.4 \lt H/L \lt 1$ for $\lambda /L = \infty$ than for the other wavelength cases.

With the submergence depth $H/L$ fixed, the vortex strength in the hydrofoil wake increases markedly as $\lambda /L$ increases, indicating a stronger interaction between the hydrofoil and the surrounding fluid. In general, the thrust generated by the hydrofoil is positively correlated with the wake vortex strength; a detailed theoretical derivation is provided in our previous work (Yu & Huang Reference Yu and Huang2021). This relationship explains why, over most of the parameter range shown in figure 8(b), the cruising speed $\overline {U_c}$ increases as $\lambda /L$ increases. In unbounded flow conditions, a larger $\lambda /L$ generally leads to stronger foil–fluid interactions and, consequently, a higher cruising speed $\overline {U_c}$ . However, in near-free-surface environments, excessively strong hydrofoil–fluid interactions induce large-amplitude free-surface waves, which exacerbate flow field boundary asymmetry and generate stronger free-surface responses. Both effects are detrimental to achieving high-speed and high-efficiency propulsion. For a fixed $H/L$ , free-surface oscillations intensify with increasing $\lambda /L$ , eventually leading to wave overturning and breaking. Consistent with figure 8(b), when the submergence depth is small (i.e. $H/L = 0.4$ , 0.6 and 0.8), increasing $\lambda /L$ beyond a certain threshold leads to a reduction rather than an increase in $\overline {U_c}$ . One contributing factor is that pronounced free-surface motions indicate enhanced energy transfer from the hydrofoil to surface waves, which may reduce the fraction of input energy available for effective propulsion. In addition, increasing $\lambda /L$ enhances flow field boundary asymmetry, leading to a significantly stronger attractive force acting on the hydrofoil towards the free surface, as shown in figure 12(b).

Furthermore, figure 15 shows that no wave breaking occurs at $H/L = 0.8$ for any considered $\lambda /L$ ; at $H/L = 0.6$ , free-surface wave breaking emerges when $\lambda /L = 2$ ; and at $H/L = 0.4$ , high-vorticity regions are generated near the free surface when $\lambda /L = 1$ . These observations indicate that decreasing the submergence depth $H/L$ increases the sensitivity of the free surface to variations in the wavelength $\lambda /L$ .

Figure 16.

Schematic of a forward jet (blue) induced by a ‘P+S’-type vortex street.

Figure 17 presents snapshots of the dimensionless streamwise velocity field, $u/U_{ref}$ , in the hydrofoil wake for 16 combinations of submergence depth $H/L$ and dimensionless wavelength $\lambda /L$ . With $\lambda /L$ fixed, a pronounced high-speed region emerges above the undulating hydrofoil as the submergence depth $H/L$ decreases, confirming that the presence of the free surface induces flow asymmetry. The formation of this high-speed region is attributed to the so-called channelling effect: as the hydrofoil upper surface approaches the free surface, the flow passage between them becomes progressively narrower, resulting in local flow acceleration. Similar channelling-induced acceleration has been widely reported in previous studies of hydrofoil propulsion near solid boundaries (Deng et al. Reference Deng, Wang, Kandel and Teng2022; Liu et al. Reference Liu, Zhang and Wang2025). In contrast, the velocity field beneath the hydrofoil remains largely unchanged across different submergence depths. With $H/L$ fixed, increasing $\lambda /L$ leads to a larger effective contact area and stronger interactions between the hydrofoil surface and the surrounding fluid. Meanwhile, the free surface exhibits stronger oscillations and improved phase-following behaviour relative to the hydrofoil upper surface, resulting in a narrower flow channel. Consequently, the intensity of the high-speed region above the hydrofoil increases markedly as $\lambda /L$ increases. For cases with small $H/L$ and large $\lambda /L$ , free-surface wave overturning and breaking occur above the rear portion of the hydrofoil, generating a low-speed region, as indicated by the black dashed boxes in figures 17(b–d) and 17(g–h). In addition, the velocity circulation associated with the trailing vortices increases significantly with increasing $\lambda /L$ , indicating stronger hydrofoil–fluid interactions.

Figure 17.

Snapshots of the dimensionless horizontal (streamwise) velocity field $u/U_{ref}$ for different submergence depths $H/L$ and wavelengths $\lambda /L$ : (a–d) $H/L=0.4$ , (e–h) $H/L=0.6$ , (i–l) $H/L=0.8$ and (m–p) $H/L=2.5$ ; (a,e,i,m) $\lambda /L=0.5$ , (b, f,j,n) $\lambda /L=1$ , (c,g,k,o) $\lambda /L=2$ and (d,h,l,p) $\lambda /L=\infty$ . All snapshots are taken at the same phase, $t/T=0.75$ , during statistically steady cruising.

To quantitatively investigate momentum transport in the wake, figure 18 presents the time-averaged streamwise velocity profiles $\overline {u}/U_{ref}$ at $X/L = 0.7$ for different submergence depths. In figure 18, the horizontal dashed line denotes the extension of the hydrofoil midline in its initial stationary configuration, while the vertical dashed line represents $\overline {u}/U_{ref} = 0$ . The quantity $\overline {u}/U_{ref}$ is directly related to momentum transport in the flow: $\overline {u}/U_{ref} \lt 0$ indicates a momentum deficit, whereas $\overline {u}/U_{ref} \gt 0$ corresponds to jet-like behaviour associated with strong momentum transport. According to previous studies, a momentum deficit leads to drag generation, whereas a momentum surplus contributes to thrust production (Jones et al. Reference Jones, Dohring and Platzer1998).

As shown in figure 18, a jet peak is observed near the wake centreline downstream of the hydrofoil for all submergence depths considered. As $H/L$ decreases, the peak magnitude of the jet gradually diminishes, which is associated with the formation of a ‘P+S’-type vortex street. In addition, a pronounced momentum deficit is observed in the fluid above the hydrofoil for $H/L = 0.4$ , followed by $H/L = 0.6$ , while the smallest momentum deficit occurs at $H/L = 0.8$ . Taken together, these results indicate that proximity to the free surface disrupts wake symmetry and is detrimental to net thrust generation by the hydrofoil.

Figure 18.

Time-averaged streamwise velocity profiles $\overline {u}/U_{ref}$ at $X/L = 0.7$ for different submergence depths $H/L$ at $\lambda /L = 1$ .

4.4. Pressure distributions

Figure 19(a–l) presents snapshots of the pressure field $C_p$ at different phase instants for three representative cases. The pressure coefficient is defined as

(4.6) \begin{align} C_p=\frac {p}{0.5\rho U_{ref}^2}, \end{align}

where $p$ denotes the gauge pressure, $\rho$ the reference fluid density and $U_{ref}$ the reference velocity. By comparing figures 19(a,d,g,j) and 19(b,e,h,k) corresponding to different values of $\lambda /L$ , it is evident that variations in $\lambda /L$ affect the phase of the hydrofoil trailing-edge motion. By comparing figures 19(a,d,g,j) and 19(b,e,h,k), it can be observed that the pressure contours for the cases $\lambda /L = 1$ and $\lambda /L = 2$ exhibit an approximate half-period phase difference. A comparison between the cases $\lambda /L = 1, H/L = 0.4$ and $\lambda /L = 2,{} H/L = 0.4$ shows that, at corresponding phase instants, the pressure magnitude in the former case is consistently lower than that in the latter, while the overall evolution of the pressure contours remains qualitatively similar. Comparing the cases $\lambda /L = 2, H/L = 0.4$ and $\lambda /L = 2, H/L = 0.6$ (figures 19 b,e,h,k and 19 c, f,i,l), it is evident that the pressure distribution and intensity on the lower surface of the hydrofoil are nearly identical at all phase instants, with the primary differences appearing near the upper surface. From figure 19(b,c,e, f), it can be observed that at $t/T = 0.25$ and $t/T = 0.5$ , when the trailing edge approaches the free surface, the pressure on the hydrofoil upper surface for $H/L = 0.4$ is significantly lower than that for $H/L = 0.6$ . In contrast, the pressure on the lower surface shows little difference. From figure 19(h,i,k,l), it can be observed that at $t/T = 0.75$ and $t/T = 1$ , the pressure distributions on both the upper and lower surfaces of the hydrofoil are nearly identical for the two submergence depths. The only notable difference is that, for $H/L = 0.4$ , the pressure in the depressed free-surface region highlighted by the red dashed box is slightly lower than that for $H/L = 0.6$ . These observations indicate that the free-surface-induced attractive force acting towards the free surface is primarily generated during the phase when the hydrofoil trailing edge approaches the free surface.

Figure 19.

Pressure field ( $C_p$ ) snapshots at different phase instants for three representative cases: (a,d,g,j) $\lambda /L = 1$ , $H/L = 0.4$ ; (b,e,h,k) $\lambda /L = 2$ , $H/L = 0.4$ ; (c, f,i,l) $\lambda /L = 2$ , $H/L = 0.6$ ; (a–c) $t/T = 0.25$ , (d– f) $t/T = 0.5$ , (g–i) $t/T = 0.75$ and ( j–l) $t/T = 1$ .

Figure 19. Long description

Panel A: A heat map showing pressure field at t/T = 0.25 for λ/L = 1, H/L = 0.4. The color scale ranges from -0.50 to 0.50, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel B: A heat map showing pressure field at t/T = 0.25 for λ/L = 2, H/L = 0.4. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel C: A heat map showing pressure field at t/T = 0.25 for λ/L = 2, H/L = 0.6. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel D: A heat map showing pressure field at t/T = 0.5 for λ/L = 1, H/L = 0.4. The color scale ranges from -0.50 to 0.50, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel E: A heat map showing pressure field at t/T = 0.5 for λ/L = 2, H/L = 0.4. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel F: A heat map showing pressure field at t/T = 0.5 for λ/L = 2, H/L = 0.6. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel G: A heat map showing pressure field at t/T = 0.75 for λ/L = 1, H/L = 0.4. The color scale ranges from -0.50 to 0.50, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel H: A heat map showing pressure field at t/T = 0.75 for λ/L = 2, H/L = 0.4. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse, with red dashed circles indicating regions of interest. Panel I: A heat map showing pressure field at t/T = 0.75 for λ/L = 2, H/L = 0.6. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse, with red dashed circles indicating regions of interest. Panel J: A heat map showing pressure field at t/T = 1 for λ/L = 1, H/L = 0.4. The color scale ranges from -0.50 to 0.50, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse. Panel K: A heat map showing pressure field at t/T = 1 for λ/L = 2, H/L = 0.4. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse, with red dashed circles indicating regions of interest. Panel L: A heat map showing pressure field at t/T = 1 for λ/L = 2, H/L = 0.6. The color scale ranges from -1.0 to 1.0, with blue indicating lower values and red indicating higher values. The pressure field shows a central concentration around the black ellipse, with red dashed circles indicating regions of interest.

Figure 20(a–p) presents the dimensionless pressure fields ( $C_p$ ) of the hydrofoil cruising wakes for 16 cases with different submergence depths $H/L$ and dimensionless wavelengths $\lambda /L$ . First, the effect of the submergence depth $H/L$ on the pressure contours is examined while the wavelength $\lambda /L$ is held fixed. As shown in figures 20(a,e,i,m) and 20(b, f, j,n), the pressure-field variations with submergence depth are concentrated primarily over the upper surface of the hydrofoil and in the region between the hydrofoil and the free surface, consistent with our previous tethered-model results (Gao et al. Reference Gao, Sun, Zhang and Liu2025). For $\lambda /L=0.5$ and 1, the low-pressure region over the upper surface progressively extends towards the free surface as $H/L$ decreases, forming a more pronounced low-pressure passage between the two boundaries. Figure 20(c,g,k,o) shows that the instantaneous hydrofoil configuration at $\lambda /L=2$ differs markedly from those in figures 20(a,e,i,m) and 20(b, f, j,n) because of the wavelength-dependent phase distribution along the hydrofoil. At $H/L=0.8$ , a weak low-pressure region forms immediately beneath the free surface and above the hydrofoil, as indicated by the dashed ellipse in figure 20(k). The resulting pressure deficit draws the free surface downward, thereby strengthening its phase-following response to the instantaneous hydrofoil deformation. As $H/L$ is further reduced to 0.6 and 0.4, the hydrofoil–free-surface interaction intensifies markedly, accompanied by expansion of the low-pressure region, stronger free-surface deformation, and wave overturning. The development of this near-surface low-pressure region therefore explains the increasing upward attraction on the hydrofoil as $H/L$ decreases and reflects the strengthened phase-following response of the free surface. This phase-following response reduces the inertial load associated with accelerating the surrounding fluid, thereby lowering the input power required for cruising. For the pure-pitching case ( $\lambda /L=\infty$ ), figure 20(d,h,l,p) show pressure-field characteristics that differ markedly from those observed in the finite-wavelength cases. At $H/L=0.4$ , the low-pressure region over the hydrofoil and adjacent to the free surface is slightly intensified, whereas the high-pressure region beneath the hydrofoil expands substantially. This behaviour is attributed to the strong unsteadiness induced by large-amplitude free-surface deformation and the associated disruption of the wake topology.

At a fixed $H/L$ , the pressure fluctuations over the hydrofoil surface intensify and the wake low-pressure region strengthens as $\lambda /L$ increases. This trend is associated with stronger fluid acceleration, trailing-edge vortex shedding and hydrofoil–free-surface interactions at larger wavelengths, thereby helping to explain the monotonic increase in the cycle-averaged input power $\overline {C_P}$ with $\lambda /L$ at a fixed submergence depth.

Figure 20.

Snapshots of the pressure field ( $C_p$ ) for different submergence depths $H/L$ and dimensionless wavelengths $\lambda /L$ : (a–d) $H/L = 0.4$ , (e–h) $H/L = 0.6$ , (i–l) $H/L = 0.8$ and (m–p) $H/L = 2.5$ ; (a,e,i,m) $\lambda /L = 0.5$ , (b, f,j,n) $\lambda /L = 1$ , (c,g,k,o) $\lambda /L = 2$ and (d,h,l,p) $\lambda /L = \infty$ . All snapshots are taken at the same phase, $t/T=0.75$ , during statistically steady cruising.

4.5. Free surface effect

Figure 21(a–l) illustrates the evolution of the free surface over one undulation cycle for three different wavelengths $\lambda /L$ , with the submergence depth fixed at $H/L = 0.4$ . The horizontal dashed line marks the free-surface level, while the instantaneous hydrofoil profile is superimposed to show its shape and position relative to the free surface. As observed in figure 21(a–l), overall, the free-surface oscillations become increasingly intense with increasing $\lambda /L$ , which is associated with stronger interactions between the hydrofoil and the free surface at larger wavelengths. A pronounced wave trough forms on the free surface above the hydrofoil; as $\lambda /L$ increases, the trough deepens and its location shifts slightly downstream. This local depression of the free surface effectively narrows the flow passage, thereby accelerating the fluid in the vicinity of the free surface. In contrast, the wave crest formed downstream of the trough widens the flow passage and suppresses the local flow velocity. Consequently, negative pressure is observed beneath the free-surface troughs, whereas positive pressure appears beneath the crests in figure 20. In addition, for the case $\lambda /L = 3$ , cavity-like region is observed in the fluid beneath the free surface. When the submergence depth $H/L$ is small, pronounced wave breaking and associated free-surface disturbances are observed, which coincide with a substantial reduction in the cruising efficiency of the hydrofoil.

Figure 21.

Evolution of the free surface over one oscillation cycle for three representative wavelengths $\lambda /L$ : (a,d,g,j) $\lambda /L = 1$ , (b,e,h,k) $\lambda /L = 2$ and (c, f,i,l) $\lambda /L = 3$ ; snapshots are shown at (a–c) $t/T = 0$ , (d–f) $t/T = 0.25$ , (g–i) $t/T = 0.5$ and ( j–l) $t/T = 0.75$ .

To assess whether free-surface wave breaking necessarily leads to a lower cruising efficiency for a near-free-surface hydrofoil compared with that at a larger submergence depth ( $H/L = 2.5$ ), a normalised efficiency metric $\eta /\eta _0$ is introduced, defined as $\eta /\eta _0 = \eta (\lambda , H)/\eta (\lambda , 2.5)$ . Here, $\eta (\lambda , H)$ denotes the cruising efficiency at a given submergence depth $H/L$ and wavelength $\lambda /L$ , while $\eta _0 = \eta (\lambda , 2.5)$ represents the corresponding reference efficiency in deep water. When $\eta /\eta _0 \gt 1$ , cruising near the free surface yields a higher efficiency than cruising at $H/L = 2.5$ , whereas $\eta /\eta _0 \lt 1$ indicates a lower efficiency relative to the deep-water case. Figure 22 illustrates the variation of $\eta /\eta _0$ in the $\lambda /L$ – $H/L$ parameter space. Green circles indicate cases without free-surface wave breaking, whereas red crosses denote cases with wave overturning and breaking. For $\lambda /L = 0.5$ , no free-surface wave breaking occurs over the entire range considered, and the cruising efficiency remains higher than that in deep water even when $H/L$ decreases to 0.4. In contrast, for $\lambda /L \gt 1$ , the cruising efficiency at $H/L = 0.4$ is consistently lower than the corresponding deep-water value. Overall, for most non-breaking cases, the cruising efficiency exceeds the deep-water reference efficiency $\eta _0$ . Notably, the boundary between wave-breaking and non-breaking regimes does not coincide with the isocontours of $\eta /\eta _0$ . For example, in the cases $\lambda /L = 2, H/L = 0.6$ and $\lambda /L = 3, H/L = 0.6$ , free-surface wave breaking occurs, yet $\eta /\eta _0$ remains greater than unity. It is also important to note that the absence of wave breaking does not necessarily lead to improved propulsion efficiency. For example, the case of $\lambda /L=8$ and $H/L=0.8$ remains within the non-breaking regime, yet exhibits lower efficiency than the corresponding deep-water case due to enhanced wave radiation and stronger hydrofoil–free-surface interactions. These observations indicate that wave breaking alone is insufficient to explain efficiency variations near the free surface. Instead, the efficiency change is governed by the combined effects of free-surface wave radiation and nonlinear wave deformation.

Figure 22.

Evolution of the normalised efficiency $\eta /\eta _0$ in the $\lambda /L$ – $H/L$ parameter space, where green circles indicate non-breaking wave conditions and red crosses denote cases with wave overturning and breaking.

In near-free-surface cruising of a hydrofoil, a clear competitive mechanism is present. As the hydrofoil approaches the free surface, the reduction in the cycle-averaged input power $\overline {C_P}$ creates the potential for improved propulsion efficiency. However, when the hydrofoil operates extremely close to the free surface, two adverse effects emerge. On the one hand, the thrust generated by the hydrofoil decreases, which is qualitatively consistent with the reduced added-mass effect; on the other hand, intensified free-surface wave activity, including wave overturning and breaking, is associated with stronger free-surface responses and reduced propulsion efficiency. From a conceptual standpoint, an optimal submergence depth $H^*/L$ must therefore exist for each wavelength $\lambda /L$ to maximise cruising efficiency. This relationship has been obtained through empirical fitting, as shown in figure 11. For $\lambda /L = 0.5$ , 1, 2 and $\infty$ , the corresponding optimal submergence depths $H^*/L$ are 0.6, 0.8, 1.0 and 1.0, respectively. To characterise the free-surface features associated with maximum cruising efficiency, free-surface profiles at $H^*/L$ and at two neighbouring submergence depths are extracted and compared for these four cases. The corresponding results are presented in figure 23. It is observed that the optimal submergence depth $H^*/L$ consistently corresponds to conditions in which free-surface wave activity converges as $H/L$ increases. Further reductions in submergence depth beyond this point result in a pronounced increase in free-surface wave generation, and the free surface exhibits stronger phase-following behaviour relative to the undulating hydrofoil, consistent with previous findings (Zheng et al. Reference Zheng, He, Zhao, Shen and Zhu2024). For near-free-surface hydrofoil cruising, the convergence of free-surface wave activity with increasing submergence depth represents a balance point of this competitive mechanism, at which power savings due to the enhanced phase-following response of the free surface can be achieved without incurring excessive energy losses from free-surface wave generation. Therefore, from an engineering perspective, suppressing free-surface wave generation is essential for achieving and maintaining efficient hydrofoil cruising.

Figure 23.

Free-surface profiles corresponding to different submergence depths $H/L$ near the optimal submergence depth $H^*/L$ for various wavelengths $\lambda /L$ .

To further quantify the strength of the free-surface response, the excess potential energy associated with free-surface elevation is evaluated. In water-wave theory, the gravitational potential energy induced by free-surface disturbances represents the energy stored in surface deformation and provides a quantitative measure of wave response and free-surface evolution (Stoker Reference Stoker1957; Henry Reference Henry2021). Therefore, the excess potential energy is used here as an indicator of free-surface deformation intensity.

The instantaneous excess potential energy of the free surface is defined as follows:

(4.7) \begin{align} E_s(t) = \frac {1}{2}\rho g \int _{\varOmega _s} \zeta ^2(x,t)\textrm{d}x , \end{align}

where $\zeta$ denotes the instantaneous free-surface elevation relative to the still-water level, and $\varOmega _s$ represents the selected free-surface region. To eliminate the phase-dependent fluctuations induced by the periodic hydrofoil motion, the cycle-averaged excess potential energy during the statistically steady stage is calculated as follows:

(4.8) \begin{align} \overline {E}_s = \frac {1}{T} \int _{t_0}^{t_0+T} E_s(t)\textrm{d}t . \end{align}

The corresponding dimensionless form is defined as follows:

(4.9) \begin{align} \widetilde {E}_s = \frac {\overline {E}_s}{0.5\rho {U_{ref}}^2L^2}. \end{align}

Figure 24 shows the variation of $\widetilde {E}_s$ with submergence depth $H/L$ for two representative wavelengths, $\lambda /L=1$ and $2$ . For both wavelengths, $\widetilde {E}_s$ decreases monotonically with increasing $H/L$ , indicating that the hydrofoil–free-surface interaction weakens as the hydrofoil moves farther away from the free surface. The rapid increase of $\widetilde {E}_s$ at shallow submergence depths indicates that stronger free-surface deformation and wave responses are generated when the hydrofoil operates closer to the interface. Moreover, the larger values of $\widetilde {E}_s$ observed for $\lambda /L=2$ indicate that longer body-wave wavelengths induce stronger free-surface responses, consistent with the enhanced wave activity observed in the instantaneous flow fields.

Figure 24.

Variation of the cycle-averaged dimensionless excess potential energy $\widetilde {E}_s$ , which characterises the free-surface response, with submergence depth $H/L$ for $\lambda /L=1$ and $2$ .

To further examine the dependence of the free-surface response on the hydrofoil–free-surface interaction strength, figure 25 shows the relationship between $\widetilde {E}_s$ and $\exp (-H/\lambda )$ . The exponential form is motivated by the attenuation of disturbances generated by submerged oscillating bodies, where the disturbance amplitude decreases exponentially with increasing distance from the free surface in the weakly nonlinear regime. The physical basis of this attenuation parameter and its incorporation into the reduced-order scaling relation are discussed further in § 4.6. For cases without wave breaking, the data exhibit an approximately collapsed trend, indicating that the free-surface response is primarily governed by the attenuation of hydrofoil-induced disturbances with increasing submergence depth. In contrast, cases involving wave breaking deviate from this trend and exhibit larger values of $\widetilde {E}_s$ , suggesting that additional nonlinear free-surface effects become significant beyond the weakly nonlinear attenuation mechanism.

Figure 25.

Relationship between the cycle-averaged dimensionless excess potential energy $\widetilde {E}_s$ and the exponential attenuation parameter $\exp (-H/\lambda )$ . Cases without and with wave breaking are distinguished to illustrate deviations induced by nonlinear free-surface effects.

These results provide quantitative evidence that the free-surface response intensifies progressively as the hydrofoil approaches the interface. For non-breaking cases, the free-surface response follows an exponential attenuation with increasing distance from the free surface, consistent with the characteristic depth dependence introduced later in the reduced-order scaling relation. However, when wave breaking occurs, additional nonlinear effects cause deviations from this simplified attenuation behaviour.

4.6. Reduced-order scaling relation for cruising speed

Building on the cruising-speed trends presented in § 4.1 and the hydrodynamic mechanisms identified above, a physics-informed reduced-order relation is developed to describe the coupled dependence of $U_c$ on the wavelength $\lambda$ and submergence depth $H$ . First, the prescribed body wave generates an intrinsic wavelength-dependent propulsive response that increases as the streamwise phase variation decreases and asymptotically approaches a finite limit at long wavelengths. Second, the free-surface pressure condition modifies this intrinsic response, whereas the hydrodynamic coupling between the hydrofoil and the free surface weakens with increasing relative submergence $H/\lambda$ . Third, at shallow submergence, enhanced wave generation, large-amplitude free-surface deformation and wake reorganisation result in additional losses of useful streamwise momentum. Together, these mechanisms motivate the development of the reduced-order scaling relation presented below.

The wavelength dependence of the propulsive response can be interpreted in terms of the reactive-force mechanism described by Lighthill’s elongated-body theory (Lighthill Reference Lighthill1971). To leading order, the cycle-averaged reactive thrust generated by an undulating body can be expressed schematically as follows:

(4.10) \begin{align} \overline {T}_r \sim \frac {m_a(L)}{2} \left [ \overline {h_t^2} - \overline {U}_c^{\,2}\overline {h_x^2} \right ]_{x=L}, \end{align}

where $m_a(L)$ denotes the added mass per unit length at the trailing edge, and $h(x,t)$ represents the prescribed lateral displacement. For the travelling-wave kinematics, $h(x,t)=A(x)\cos (kx-\omega t)$ with $k=2\pi /\lambda$ , the cycle-averaged temporal and spatial derivatives evaluated at the trailing edge are $\overline {h_t^2}=\omega ^2A_T^2/2$ and $\overline {h_x^2}=[(A'_T)^2+k^2A_T^2]/2$ , respectively, where $A_T=A(L)$ and $A'_T=A'(L)$ . The body-wave wavelength therefore enters the reactive-force estimate through the streamwise phase-gradient term $k^2A_T^2$ . As $\lambda$ increases, this phase-gradient term decreases, and the prescribed deformation approaches the phase-aligned pitching-like limit. Consequently, the wavelength-dependent propulsive response increases with $\lambda$ but remains finite in the long-wavelength limit as $k\rightarrow 0$ . Lighthill’s relation is invoked here to establish the expected monotonic increase and asymptotic saturation of the wavelength-dependent response, rather than to quantitatively predict the thrust generated by the present two-dimensional hydrofoil. Based on these limiting behaviours, the intrinsic wavelength-dependent response is represented by the bounded saturation function $1-\exp (-\lambda /\lambda _0)$ , where $\lambda _0$ denotes the characteristic wavelength controlling the approach to the long-wavelength limit. Nonlinear least-squares fitting of the data at $H/L=2.5$ yields $\lambda _0/L=2.0$ . The data at $H/L=2.5$ are taken to represent the weak free-surface-interaction limit within the parameter range considered, because free-surface effects are relatively small at this submergence depth. The resulting value of $\lambda _0/L$ is then held fixed in the global regression.

The intrinsic wavelength-dependent response is further modified by free-surface effects. According to linear potential-flow theory, the velocity-potential and pressure amplitudes induced by a harmonic subsurface disturbance with horizontal wavenumber $k_s$ decay exponentially with submergence as $\exp (-k_sH)$ (Wehausen & Laitone Reference Wehausen, Laitone and Flügge1960). The wavenumber of the radiated free-surface waves is determined by the actuation frequency, cruising speed and the gravity-wave dispersion relation. In the present reduced-order model, the imposed body-wave wavelength is adopted as the characteristic streamwise length scale, yielding $k_s=O(1/\lambda )$ . The hydrofoil–free-surface coupling is therefore represented by $\exp (-bH/\lambda )$ , where the dimensionless coefficient $b$ accounts for differences between the imposed body-wave and radiated-wave scales, as well as finite-geometry and dispersion effects. Accordingly, the fraction of the intrinsic wavelength-dependent response retained as the free-surface influence diminishes is represented by the complementary factor $1-\exp (-bH/\lambda )$ , which approaches unity with increasing $H/\lambda$ .

The exponential factor captures the depth dependence of the pressure coupling between the hydrofoil and the free surface but does not account for nonlinear effects arising at shallow submergence. As $H/L$ decreases and $\lambda /L$ increases, wave radiation and free-surface deformation become increasingly pronounced. Meanwhile, wave overturning, incipient breaking, vortex splitting and wake asymmetry increasingly redistribute streamwise momentum. In particular, the emergence of an asymmetric ‘P+S’ wake and a forward-directed jet indicates that an increasing fraction of the generated momentum is diverted away from useful propulsion. The relative contribution of these additional wave- and wake-related losses diminishes in the weak-interaction limit, $\lambda /H\rightarrow 0$ . Their leading-order contribution is therefore assumed to scale with $c\lambda /H$ , where $c$ is a dimensionless loss coefficient. By interpreting $c\lambda /H$ as the ratio of the additional loss to the useful cruising-speed response, the retained fraction is obtained as $1/(1+c\lambda /H)$ . This rational form remains positive and bounded over the parameter range considered and admits the asymptotic expansion $1-c\lambda /H+O[(\lambda /H)^2]$ as $\lambda /H\rightarrow 0$ .

The exponential and rational factors capture distinct aspects of the free-surface influence. The exponential factor describes the recovery of the intrinsic response as the hydrofoil–free-surface coupling weakens with increasing submergence. The rational factor accounts for additional losses of useful streamwise momentum associated with wave generation, large-amplitude free-surface deformation, and asymmetric wake dynamics. Combining these free-surface mechanisms with the intrinsic wavelength-dependent response gives

(4.11) \begin{align} \frac {\overline {U}_c(\lambda ,H)}{U_{ref}} = a \left [ 1-\exp \left (-\frac {\lambda }{\lambda _0}\right ) \right ] \frac { 1-\exp \left (-bH/\lambda \right ) }{ 1+c\lambda /H } +d . \end{align}

The multiplicative structure reflects the sequential effects of the mechanisms incorporated into the reduced-order model. The prescribed body wave first generates an intrinsic wavelength-dependent response in cruising speed. This response is subsequently modified by the depth-dependent hydrofoil–free-surface coupling, after which wave- and wake-related losses reduce the fraction available for effective forward propulsion. The coefficient $a$ determines the overall magnitude of the wavelength-dependent response, whereas the constant $d$ accounts for the component of cruising speed not captured by the leading-order wavelength–depth dependence. This offset is specific to the fixed Reynolds number and prescribed kinematics considered here and is therefore not given an independent asymptotic interpretation.

Figure 26 illustrates the scaling relationship for the cruising speed $\overline {U}_c$ as a function of the submergence depth $H$ and wavelength $\lambda$ over the ranges $0.5\leq \lambda /L\leq 8$ and $0.4\leq H/L\leq 2.5$ . With $\lambda _0/L$ fixed at $2.0$ , a global nonlinear least-squares optimisation over all finite-wavelength cases yields $a=1.50$ , $b=1.5$ , $c=0.28$ and $d=0.35$ . Here

(4.12) \begin{align} \frac {\overline {U}_c(\lambda ,H)}{U_{ref}} = 1.50 \left [ 1-\exp \left (-\frac {\lambda }{2L}\right ) \right ] \frac { 1-\exp \left (-1.5H/\lambda \right ) }{ 1+0.28\lambda /H } +0.35 . \end{align}

Figure 26.

Overall scaling collapse of the cruising-speed data based on the physics-informed reduced-order relation. The symbols denote numerical results for different wavelengths, while the dashed line represents the linear relationship predicted by the reduced-order scaling relation.

Although the reduced-order relation captures the overall scaling behaviour of the cruising speed, the deviations between the predictions and numerical results are not uniformly distributed across the parameter space. To identify the physical origins of these deviations, the reduced-order relation is reconstructed in the original parameter space of wavelength and submergence depth, as shown in figure 27.

Figure 27.

Comparison between numerical results (symbols) and reduced-order predictions (lines) of cruising speed as functions of wavelength and submergence depth. (a) Variation of $\overline {U}_c/U_{ref}$ with submergence depth $H/L$ for different body-wave wavelengths $\lambda /L$ ; (b) variation of $\overline {U}_c/U_{ref}$ with body-wave wavelength $\lambda /L$ for different submergence depths $H/L$ .

Figure 27. Long description

Panel A: A line graph shows the variation of time-averaged streamwise velocity profiles with submergence depth for different body-wave wavelengths. The x-axis represents the submergence depth (H/L) and the y-axis represents the normalized streamwise velocity (Uc/Uref). Different symbols and lines represent various body-wave wavelengths (λ/L = 0.5, 0.75, 1, 1.5, 2). The graph includes a horizontal dashed line denoting the extension of the hydrofoil midline and a vertical dashed line representing a specific value. The R-squared values for the fits are also shown. Panel B: Another line graph shows the variation of time-averaged streamwise velocity profiles with body-wave wavelength for different submergence depths. The x-axis represents the body-wave wavelength (λ/L) and the y-axis represents the normalized streamwise velocity (Uc/Uref). Different symbols and lines represent various submergence depths (H/L = 0.6, 0.8, 1.5, 2.5). The R-squared values for the fits are also shown.

As shown in figure 27, the deviations are mainly concentrated in two regimes. The first occurs at the shortest wavelengths ( $\lambda /L=0.5$ and $0.75$ ), where the prescribed deformation introduces stronger geometric nonlinearities. As the wavelength decreases, the streamwise curvature increases and the variation in the actual hydrofoil centreline arclength becomes more significant. Since the wavelength-dependent response in the present scaling relation is derived from Lighthill’s elongated-body framework using the undeformed streamwise coordinate as the characteristic length scale, the reduced-order description becomes progressively less accurate as geometric deformation increases. The second regime corresponds to shallow submergence depths ( $H/L=0.4$ ), where pronounced free-surface deformation and nonlinear wave dynamics become significant. The depth-dependent correction terms in the proposed scaling relation are motivated by the attenuation and loss mechanisms associated with hydrofoil–free-surface interactions. Specifically, the factor $1-\exp (-bH/\lambda )$ is motivated by the exponential attenuation of subsurface disturbances in linear free-surface theory and describes the variation of hydrofoil–free-surface coupling strength with submergence depth. Meanwhile, the rational correction factor $1/(1+c\lambda /H)$ is introduced as a low-order representation of additional wave-radiation and wave-breaking-related losses. However, when the hydrofoil operates close to the free surface, both approximations become less accurate because large-amplitude free-surface deformation, wave overturning and nonlinear wave–wake interactions introduce complex mechanisms that cannot be fully captured by the present reduced-order formulation.

To further examine the validity range of the reduced-order relation, the shallowest submergence cases ( $H/L=0.4$ ) are excluded from the regression, resulting in a more appropriate validity range and an improved $R^2$ value. Within the regime of $0.5\leq \lambda /L\leq 8$ and $0.6\leq H/L\leq 2.5$ , the free-surface response remains primarily governed by the mechanisms incorporated into the present model, including the depth-dependent attenuation of hydrofoil–free-surface coupling and the associated modification of wake dynamics. Therefore, this range is identified as the applicable regime for the proposed scaling relation. The corresponding scaling collapse shown in figure 28 is further improved within this regime, with the coefficient of determination increasing from $R^2=0.86$ to $R^2=0.93$ .

Figure 28.

Scaling collapse of the cruising-speed data based on the physics-informed reduced-order relation for $0.6\leq H/L\leq 2.5$ . Symbols denote numerical results for different wavelengths, whereas the dashed line represents the corresponding linear relationship predicted by the reduced-order scaling relation.