12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

HYDRODYNAMICS OF SLENDER UNDULATING FOIL

Li-Ming Chao
Center for Turbulence Control
Harbin Institute of Technology (Shenzhen)
University Town, Shenzhen, 518055, China
chaoliming92@gmail.com

Md. Mahbub Alam
Center for Turbulence Control
Harbin Institute of Technology (Shenzhen)
University Town, Shenzhen, 518055, China
alam@hit.edu.cn

ABSTRACTPROBLEM DESCRIPTION

The thrust generation and wake structures of a slender undulating foil have been studied, where the Strouhal number St,
Reynolds number Re and dimensionless wavelength λ influences are considered. The relationship of time-mean thrustwith St, Re and λ is presented, suggesting that the propulsive
force increases with enlarging St, Re and λ . As such, the dragthrust boundary advances with increasing them. Seven types ofwake structures produced by the foil are identified, discussed,and connected to thrust generation, showing how St, Re and λaffect fluid dynamics, wake transition and vortex strength.

A schematic of the flow configuration is presented in figure 1(a), where a two-dimensional NACA0012 foil with chordlength of L is placed in a uniform flow of velocity U in thex-direction. The leading edge of the foil is located at the originof the Cartesian coordinate system. The following equation is
employed to control the foils undulation,


h
t i
x


 y(x,t) = A(x) sin 2π( λ − T )Z Ls (t) q{1 + (∂ y/∂ x)2 } dx = L(1)0

INTRODUCTION
After a long history of evolution, aquatic animals havemastered an exquisite capacity to control their body and theflow around themselves to efficiently cruise in water (Fish& Lauder, 2006). Slender swimmers, such as anguilliform,
carangiform and subcarangiform swimmers, usually employ
the undulating body-caudal-fin (BCF) propulsive strategy topropel themselves and overcome the fluid resistance, where
the swimmers undulating BCF strategy can be considered as
a combination of two waves. One is a cyclic change of thecurved body producing a lateral wave propagation in the caudal direction, and the other is every single point on the bodyundergoing a sinusoidal track in a horizontal plane as the consequence of the wave propagation (Gray, 1933). For different
slender swimmers, they may use different swimming strategies
to bend their bodies. Even one specific swimmer has different
deforming techniques at different conditions, such as the fish
uses C-start deformation during the survival behavior (Gazzola
et al., 2012). However, to the authors knowledge, no systematic investigation has been made to understand the effect of
deforming characteristics on the propulsive performance of an
undulating swimmer. In this work, we, therefore, attempt tosystematically investigate the hydrodynamic performance ofan undulating foil with varying foil kinematics (Strouhal number), fluid properties (Reynolds number), and foil deforming
characteristics (wavelength), with the foil length keeping constant to replicate the native slender swimmer.

where A(x) is the swimming amplitude function, usually considered as a polynomial A(x) = a0 + a1 x + a2 x2 , x is measured
from the leading edge of the foil, y is the lateral displacementof the foil midline, λ is the wavelength defined in figure 1(b),t is the instantaneous time, T is the undulating period, andLs (t) denotes projected length of the foil on the x-axis (streamwise direction), giving Ls (t) ≤ L. Without a loss of generality,a0 = a2 = 0 and a1 = 0.05 are employed to determine the foil
undulation.
To follow the kinematics of a real swimmer, both the foil
midline length and the included area are kept constant, withthe foil midline computed from equation (1). The foil surface
is given by maintaining the normal distance (e.g. hs at an arbitrary point on the midline) from the midline constant for the
deformed and undeformed foils. For example, consider that hs
is the distance between point A on the foil midline and point
B (normal to the midline) on the foil surface when the foil isstraight, undeformed (figure 1c). At each time step, the line
AB normal to the deformed midline can be computed from thedifferential form of equation (1). The new location of point B
on the foil surface is then pinpointed with AB = hs . A typicalenvelope of the undulation is illustrated in figure 1(d), where
it is easy to understand that Ls (t) is a time-dependent variable,with Ls (t) = L corresponding to the straight midline (dashedblue line), not undulating.

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 1. (a) Computational domain and mesh system, (b) traveling wavy motion, (c) method for achieving foil surface and (d) foilparameter.

THRUST GENERATION

Strouhal number and Reynolds number are defined as:
St = 2AL /UT(2)

Re = UL/ν

(3)

Figures 2(a)-(e) illustrate the dependence of C̄T on St, Re
and λ . The solid black line represents the drag-thrust boundary ( C̄T = 0). For a given λ , the C̄T increases from negative (lower-left corner) to positive (upper right corner) when
St and Re are increased. The traveling wavy foil cannot produce thrust when Re and St are low but can when Re and St are
sufficiently high, depending on λ . Interestingly, C̄T is equallydependent on St and Re in the drag regime ( C̄T < 0) but moresensitive to St in the thrust regime ( C̄T > 0). The observation
indicates that viscous (Re) and inertia (St) effects both play important roles in fluid dynamics in the drag regime while the inertia effect dictates the fluid dynamics in the thrust regime. For
a given set of St and Re, an increase in λ leads to an enhancement in C̄T . As such, the drag-thrust boundary shifts to thesmaller St and Re when λ is increased. With the same St and
λ , the C̄T gets higher at a higher Re, i.e. having acceleration iseasier at a higher speed. The reason is that the viscous effectweakens and the added mass lessens when Re is increased. On
the other hand, a foil with a smaller λ hinders the streamwise
flow over the surface, resulting in a smaller thrust generation.
When the stationary wavy foil (St = 0) remains travelingwavy shape (as sketched in figure 2a-e), it would experiencea thrust (C̄T 0 ) that is negative (i.e. drag). Apparently, C̄T 0
would dependent on Re as shown in figure 2(f), where declines with Re as C̄T 0 = −6.13Re−0.6 . Das et al. (2016)
showed that a straight NACA0012 foil at zero attack angle incurs C̄T 0 = −5.96Re−0.56 that is smaller in magnitude than thepresent counterpart. It is plausible as an undulated foil would
experience a higher drag than a straight foil.
The linear least-squares regression is used to find the relationship between C̄T , St and λ , giving C̄T = c1 St 3 λ +c2 . Afterreprocessing the data obtained from the numerical simulations,the relationship of C̄T with Re, St and λ can be unified as

andwhere AL (= 0.05L) is the tail amplitude and ν is the kinematic
viscosity of the fluid. The time-mean thrust coefficient is defined as
C̄T =

1
T

Z t+T
t

CT dt =

1
T

Z t+T
t

−

Fx (t)
dt
0.5ρU 2 L

(4)

where CT is the instantaneous thrust coefficient, Fx (t) is thedrag (x-direction) force acting on the foil, and ρ is the density of the fluid. Therefore, C̄T > 0 and C̄T < 0 represent the
foil generating time-mean thrust and experiencing time-meandrag, respectively, while C̄T = 0 is the drag-thrust transition
boundary.

NUMERICAL METHOD
The unsteady flow field around the foil is simulated using the commercially available CFD package Fluent 14.0. The
flow domain consists of three grid zones, including a regularzone with a high resolution in unstructured grids around thefoil (zone 1, figure 1a), zone 2 is of high resolution where alarge velocity gradient is expected and zone 3 with a medium
resolution. The locomotion of foil in the computational domain is achieved using a dynamic grid scheme implemented
in FLUENT. At each updated time instant, the foil motion isdetermined using equation (1). The inlet velocity boundary is
at a distance of 20L from the leading edge of the foil, withu = (U, 0) and ∂ p/∂ x = 0 while the outflow boundary is 36L
downstream, with ∂ u/∂ x = (0, 0) and ∂ p/∂ x = 0 (no-stressoutflow boundary conditions). The slip wall condition is used
for the upper and lower boundaries located symmetrically 40L
apart (Bos et al., 2008).
To investigate the thrust generation and wake structuresof the slender undulating foil, the numerical simulations are
performed at St = 0.1 1.0 with ∆St = 0.1, at Re = 50, 100, 250,
500, 750, 1000, 1250, 1500, 1750, and 2000, and at λ = 0.50,
0.67, 1.0, 1.5 and 2.0.

C̄T = 0.36Re0.208 St 3 λ − 6.13Re−0.6(5)Equation (5) displays how St, Re, and λ affect the propulsive
force on the foil. It is easy to understand from this equationthat the thrust generation (C̄T > 0) is largely determined by
St and λ , especially by St. The C̄T enhances slowly withincreasing Re, which essentially reflects the viscous flow resistance weakening with increasing Re. Particularly, when
Re is large enough, the second term on the right-hand side
approaches zero, which gives C̄T = 0.36Re0.208 St 3 λ . When2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 2. (a-e) Dependence of time-mean thrust on Re, St and λ for the undulating foil. (f) Dependence of time-mean thrust force ofa stationary foil on Re.

Figure 3. (a) Relationship of C̄T with St, Re and λ for the undulating foil. (b) Dependence of drag-thrust boundary on St, Re and λfor the undulating foil. The dash lines denote the proportional function with a coefficient equal to one.

the foil is stationary (St = 0), equation (5) degenerates intoC̄T = −6.13Re−0.6 , i.e. the C̄T 0 scaling. To see how much
equation (5) is capable of collapsing all simulated data, C̄T ispresented against 0.36Re0.208 St 3 λ − 6.13Re−0.6 in figure 3(a).
The thrust generation is found to be dependent on λ . It
is worth finding the relationship between St, Re, and λ corresponding to the drag-thrust boundary (figures 2a-e). The dragthrust boundary can be obtained from equation (5) by plugging
C̄T = 0 , i.e.
St = 2.57Re−0.27 λ −1/3

(2020) theoretically showed that the relationship between C̄T
and St can be considered as C̄T ∝ St 3 .
Figure 4 displays CT histories in one oscillation period atfour selected points: A (St = 0.2, Re = 100) in drag regime, B(St = 0.4, Re = 250) close to the drag-thrust boundary at high
λ , C (St = 0.6, Re = 600) close to the drag-thrust boundary
at small λ , and D (St = 0.8, Re = 1000) in thrust regime, allmarked in the St − Re map (figure 2). Similar to the previous works on the pitching foil (Koochesfahani, 1989), two CT
peaks are observed in the upstroke and downstroke, respectively. The traveling wavy foil experiences drag during theentire period at St = 0.2, Re = 100 (figure 4a). For the given
St and Re, the CT peak becomes larger when λ increases from0.50 to 2.0. At intermediate St and Re values (points B and C),
the foil undergoes both drag and thrust in the period (figures
4b, c). On the other hand, at a higher St and Re, only thrustacts on the foil in the entire period (figure 4d). The CT amplitude (the gap between maximum and minimum) grows withincreasing St, Re and λ (figures 4a-d). A smaller λ , therefore,
engenders a more steady thrust, which explains the swimming
behavior of the anguilliform swimmer. On the other hand, the(6)Equation (6) reflects that the drag-thrust boundary in the
St − Re domain advances with increasing λ , as illustrated infigure 3(b). The present results and relationship conform withthe biological observations (Gazzola et al., 2014) while the effect of the deforming strategy (λ ) on the cruising locomotion
is firstly presented in this work. Based on the theoretical model
involving the Euler angular, Coriolis, and centrifugal accelerations of a cone of fluid around the foil, Alam & Muhammad

3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 4.

The instantaneous thrust coefficient for the traveling wavy foil in one oscillation period.

Figure 5. The wake structure map with the black solid line denoting the drag-thrust boundary, red dash line denoting the KV-RKV
wake boundary and blue dash-dotted line denoting the RKV-s-RKV wake boundary.

hind a quiescent foil (St = 0). With a little increase in St and/or
Re, the steady wake transmutes to q-KV wake (figure 6b) or qRKV wake (figure 6c) representing a transition from a steadywake to a classical KV or RKV wake. The q-KV and q-RKV
wakes form at low and high St, respectively, the former prevailing for smaller λ (≤ 1.0) only (figure 5). Here, Kármánor reverse Kármán vortices form immediately behind the foil
while vanishing rapidly to generate a steady wake downstream(figures 6a-c). The generation of steady wake followed byKrmn wake was also found in the wake of a fixed ellipticalcylinder at Re = 100 − 150 in Shi et al. (2020a) and Shi et al.
(2020b). Depending on λ , when Re is further increased, KV
wake forms, featuring negative and positive vortices above andbelow the foil symmetry line correspondingly (figure 6d). The
2S aligned wake refers to the case where both positive andmaximum instantaneous thrust generated for a higher λ explicates why the slender swimmer usually employs the C-start
strategy during the survival behavior (Gazzola et al., 2012).

FLOW MAP
Seven distinguished wake structures are identified in the
St − Re domain examined: steady wake, quasi Kármán vortex (q-KV) wake, quasi reverse Kármán vortex (q-RKV) wake,Kármán vortex (KV) wake, 2S aligned wake, reverse Kármánvortex (RKV) wake, and slanted reverse Kármán vortex (sRKV) wake. Their presence in the St − Re domain and representative flow structures are shown in figures 5 and 6. The
steady wake is observed at St ≤ 0.6 and Re ≤ 250 (figures 6a,5), where the wake is steady, similar to the flow structure be-4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 6. Typical instantaneous vorticity structures: (a) steady wake; (b) quasi Kármán vortex (q-KV) wake; (c) quasi reverse Kármánvortex (q-RKV) wake; (d) Kármán vortex (KV) wake; (e) 2S aligned wake; (f) reverse Kármán vortex (RKV) wake; (g) downward
slanted reverse Kármán vortex (ds-RKV) wake; and (h) upward slanted reverse Kármán vortex (us-RKV) wake.

Figure 7.

Contours of spanwise vorticity ω ∗ = ωL/U. The red and blue colors denote positive and negative ω ∗ , respectively.

negative vortices lie on the symmetric line (figure 6e). The
2S aligned wake is considered as the boundary between KV
and RKV wakes, extensively discussed in the previous works(Godoy-Diana et al., 2008, 2009; Marais et al., 2012; Denget al., 2015, 2016; Andersen et al., 2017; Chao et al., 2019,2021; Muhammad et al., 2020). The emergence of 2S aligned
wake at the KV-RKV boundary (red dash line) precedes thedrag-thrust boundary (black solid line) in the St - Re map (figure 5). The same is true for a pitching foil (Godoy-Diana et al.,2008). With a further increase in St and/or Re, RKV wakeemerges, where the negative and positive vortices nestle oppositely to those of KV wake (figure 6f). In s-RKV wake, thesymmetric line of two vortex streets does not coincide withthe foil symmetric line. The slant direction can be upward
or downward. The downward slanted RKV is referred to as
ds-RKV wake (figure 6g) while upward slanted RKV refers
to us-RKV wake (figure 6h). The slant direction of s-RKV
wake is also affected by λ , see figure 6(g, h). Generally, theslant direction strongly depends on the first vortex pair generated by the foil (Zheng & Wei, 2012). Here, both ds-RKV
and us-RKV wakes are considered as s-RKV wake because
the physical mechanisms of the ds-RKV and us-RKV wakes
are identical. The formation of the slanted wake is also common in flow-induced vibrations of cylinders (Bhatt & Alam,2018).

a traveling wavy foil, particularly the transition between two
different types of flow structures. For example, the KV-RKV
wake transition bordered by 2S aligned wake advances withincreasing λ . Similarly, RKV-s-RKV wake transition (bluedash-dotted line) also shifts to smaller St and/or Re when λ isincreased. In addition, there is no s-RKV wake observed at λ =0.50 (figure 5a). Physically, λ describes how the foil locomotion propagates from the leading to trailing edges. Therefore,
the foils horizontal acceleration at the trailing edge declines
with the decrease in λ . As a result, the convection velocityof the vortices with respect to U also decreases with declining
λ . In the previous work, Godoy-Diana et al. (2009) showedthat the generation of s-RKV wake is significantly dependent
on the value of this convection velocity, where a larger value of
this convection velocity may induce s-RKV wake. A smaller λthus postpones the transition between RKV and s-RKV wakes.

EFFECT ON VORTEX STRENGTH
To evaluate quantitatively the influence of St, Re and λon vortex strength at different regimes, figure 7 presents dimensionless spanwise vorticity ω ∗ (= ωL/U) structures at thefour selected points located in the St − Re plane (figure 5). Theω ∗ increases with St, Re and λ , which is consistent with thechange in C̄T in figure 2. A similar relationship of the spatial wavelength ∆d with St, Re and λ is noticeable, having the

Figure 5 also surfaces how λ affects the wakes behind512th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
in hovering insect flight. Journal of Fluid Mechanics 594,
341–368.
Chao, L.-M., Alam, M.M., Ji, C. & Wang, H. 2021 Flow mapof foil undergoing fast and slow pitching. Physics of Fluids
33, 101902.
Chao, L.-M., Pan, G., Zhang, D. & Yan, G.-X. 2019 Numerical investigations on the force generation and wake structures of a nonsinusoidal pitching foil. Journal of Fluids andStructures 85, 27–39.
Das, A., Shukla, R. K. & Govardhan, R. N. 2016 Existence ofa sharp transition in the peak propulsive efficiency of a lowre pitching foil. Journal of Fluid Mechanics 800, 307–326.
Deng, J., Sun, L. P. & Shao, X. M. 2015 Dynamical featuresof the wake behind a pitching foil. Physics Review E 92,
063013.
Deng, J., Sun, L. P., Teng, L. B., Pan, D. Y. & Shao, X. M. 2016
The correction between wake transition and propulsive efficiency of a flapping foil: A numerical study. Physics of
Fluids 92, 094101.
Fish, F. E. & Lauder, G. V. 2006 Passive and active flow control by swimming fishes and mammals. Annual Review of
Fluid Mechanics 38, 193–224.
Gazzola, M., Argentina, M. & Mahadevan, L. 2014 Scaling
macroscopic aquatic locomotion. Nature Physics 10, 758–
761.
Gazzola, M., Van Rees, M. & Koumoutsakos, P. 2012 C-start:
optimal start of larval fish. Journal of Fluid Mechanics 698,
5–18.
Godoy-Diana, R., Aider, J.-L. & Wesfreid, J. E. 2008 Transitions in the wake of a flapping foil. Physics Review E 77,
016308.
Godoy-Diana, R., Marais, C., Aider, J.-L. & Wesfreid, J. E.
2009 A model for the symmetry breaking of the reverse
benard-von karman vortex street produced by a flapping
foil. Journal of Fluid Mechanics 622, 23–32.
Gray, J. 1933 Studies in animal locomotion: I, the movementof fish with special reference to the eel. Journal of Experimental Biology 10, 88–104.
Koochesfahani, M.M. 1989 Vortical patterns in the wake of an
oscillating airfoil. AIAA Journal 27, 1200–1205.
Marais, C., Thiria, J. E., Wesfreid, J. E. & Godoy-Diana, R.
2012 Stabilizing effect of flexibility in the wake of a flapping foil. Journal of Fluid Mechanics 710, 659–669.
Muhammad, Z., Alam, M.M. & Noack, B. R. 2020 Efficientthrust enhancement by modified pitching motion. Journalof Fluid Mechanics 933, A13.
Shi, X. Y., Alam, M. M. & Bai, H. L. 2020a Wakes of elliptical cylinders at low reynolds number. Journal of Fluid
Mechanics 30, 525–532.
Shi, X. Y., Alam, M. M., Bai, H. L. & Wang, H. 2020b The
effect of reynolds number on the elliptical cylinder wake.
International Journal of Heat and Fluid Flow 82, 108553.
Zheng, Z. C. & Wei, Z. 2012 Study of mechanisms and factorsthat influence the formation of vortical wake of a heaving
airfoil. Physics of Fluids 24, 103601.

same correspondence with C̄T . The growth of ∆d with St andRe is ascribed to the increased velocity of vortices when St andRe are increased. On the other hand, it is easy to understandthat a smaller λ would lead to a smaller velocity of vortices,hence a smaller ∆d and C̄T .

CONCLUSION
A systematic numerical study is performed on hydrodynamic performance of a fish-like foil undergoing travelingwavy undulation for Strouhal number St = 0.11.0, Reynolds
number Re = 502000, and non-dimensional wave length λ =0.502.0 with the foil length keeping constant to replicate thenative slender swimmer.
The time-mean thrust C̄T is mapped on St − Re plane fordifferent λ values. A unifying mechanistic relationship ofC̄T with St, Re and λ is obtained as C̄T = 0.36Re0.208 St 3 λ −6.13Re−0.6 . The C̄T enhances with increasing St, Re and λ . A
predator-prey game at a sufficiently high Re relies on St and λonly as Re effect on C̄T is negligible at a high Re. The unifying relationship also provides the scaling law of the drag-thrustboundary ( C̄T = 0) that is marked as St = 2.57Re−0.27 λ −1/3 .
An increase in λ advances the drag-thrust boundary to smaller
Re and St values. The decreasing λ leads to a smaller fluctuation in thrust. A swimmer may have a more steady thrust when
λ is smaller, consistent with the swimming behavior of the anguilliform swimmer. On the contrary, a larger λ enhances the
maximum instantaneous thrust which explains why a slenderswimmer adopts the C-start strategy during the survival behavior.
Seven distinct wake structures are identified, namely
steady, q-KV, q-RKV, KV, 2S, RKV and s-RKV wakes.
Among them, steady, q-KV, KV and 2S wakes are always dragproducing while q-RKV and RKV wakes can produce bothdrag and thrust depending on λ . On the other hand, s-RKV
wake is always thrust producing. When λ is smaller, s-RKV
wake transmutes into RKV, hence s-RKV wake disappears at
the λ = 0.50 examined. A smaller λ enhances the stabilityof the wake not to be slanted, deferring RKV-s-RKV wakeboundary, KV-RKV wake boundary, and drag-thrust boundary.

REFERENCES
Alam, M.M. & Muhammad, Z. 2020 Dynamics of flow arounda pitching hydrofoil. Journal of Fluids and Structures 99,
103151.
Andersen, A., Bohr, T., Schnipper, T. & Walther, J. H. 2017
Wake structure and thrust generation of a flapping foil in
two-dimensional flow. Journal of Fluid Mechanics 812, R4.
Bhatt, R. & Alam, M.M. 2018 Vibration of a square cylinder submerged in a wake. Journal of Fluid Mechanics 853,
301–332.
Bos, F. M., Lentink, D., Van Oudheusden, B. & Bijl, H. 2008
Influence of wing kinematics on aerodynamic performance6