Conference Paper
2022

Development of a Finite Element Solver Including a Level-Set Method for
Modeling Hydrokinetic Turbines

Ahmed A. Hamada1, Mirjam Fürth1
1Department of Ocean Engineering, Texas A$\&$M University, College Station,
TX, 77843, USA

## Abstract

Hydrokinetic flapping foil turbines in swing-arm mode have gained considerable
interest in recent years because of their enhanced capability to extract
power, and improved efficiency compared to foils in simple mode. The
performance of foil turbines is closely linked to the development and
separation of the Leading-Edge Vortex (LEV). To accurately model the formation
and the separation of the LEV on flapping foils, a purpose-built 2D numerical
model was developed. The model is based on the weighted residual Finite
Element Method (FEM); this is combined with an interface capturing technique,
Level-Set Method (LSM), which was used to create a reliable and high-quality
numerical solver suitable for hydrodynamic investigations. The solver was
validated against well-known static and dynamic benchmark problems. The effect
of the mesh density was analyzed and discussed. This paper further covers an
initial investigation of the hydrodynamics of flapping foil in swing-arm mode,
by studying the structure of the vortex around a NACA0012 foil. The presented
method helps to provide a better understanding of the relation between the
Leading-Edge Vortex creation, growth, and separation over the flapping foil in
swing-arm mode and the extracted power from a hydrokinetic turbine.

Keywords: Finite Element Method; Level-Set Method; current/tidal energy
harvesting; hydrokinetic foil turbine; swing-arm mode.

## 1\. Introduction

The demand for clean electricity from alternative power sources is rising as
more consumers value sustainability [1]. Ocean currents have great potential
to generate reliable clean power [2], and a range of market-ready hydrokinetic
turbines have been developed over the last twenty years. For example, the
Oscillating Marine Current Energy Converter (“Stingray”) [3, 4, 5], and the
1kW Oscillating Hydrofoil Energy Harvesting System [6]. Furthermore, current
energy has less environmental impact than conventional hydroelectric
facilities, such as dams, diversions, and pumped storage [7].

Flapping foils allow for energy extraction in free surface waves, open
channels, and uniform currents [8]. Wu et al. [9] introduced the concept of
generating energy using an oscillating wing in free-surface waves in the
1970s. They showed that a flapping foil generates power by performing two main
motions; heave and pitch. Oscillating foil turbines generate a leading-edge
vortex (LEV), which leads to higher forces. As a result, they create higher
lift forces than conventional rotary turbines in which the flow remains
attached to the turbine blades [10]. Furthermore, by adopting a swing-arm
mode, in which the heave motion is performed by an arm rotation, the
efficiency can further be improved compared to simple mode, in which the heave
motion is performed simply up and down [11]. The swing-arm mode flapping foil
(a) is shown in Figure 1 in comparison with a simple mode foil (b). Karbasian
et al. [12] showed that the swing-arm mode changes the LEV’s

Figure 1: A flapping foil system with different motion patterns: (a) Swing-arm
power extraction mode, and (b) Simple power extraction mode. In each mode, the
associated angle of attack $\alpha$, and force directions during the flapping
cycle are shown

creation, growth, separation, and delay (in time). Furthermore, more power can
be extracted by selecting the optimum value of the swing factor. Numerical
simulations allow for deep investigations of the formation LEV around a
flapping foil, which may require very advanced motion capture equipment and
other experimental setups [13, 14]. Thus, numerical methods allow faster
evaluation of new designs compared to experiments.

In this paper, the water flow around the flapping foil is modeled with the
incompressible Navier–Stokes (NS) equations. To numerically solve the NS
equations, a discretization technique has to be used, such as Finite
Difference Method (FDM), Finite Element Method (FEM), or Finite Volume Method
(FVM). FEM works on the weak (integral) form of the partial differential
equations (PDE), with a lower derivative order compared to the strong
(differential) form. The order of the approximated solution in the weak form
is half of that in the strong form [15]. This makes FEM more powerful than
FDM, which discretizes the strong form directly. FEM has an advantage compared
to other numerical methods, such as FDM, and FVM, in the application of
boundary conditions [16]. Because FEM applies only physical boundary
conditions directly for the nodal values of the field variables, there is no
need for numerical boundary conditions. In addition, the natural and mixed
boundary conditions are automatically inherent in the weak form of the PDE
[17]. Further, FEM is better at handling unstructured meshes [18].

There are two main types of meshes, moving and fixed [19]. In moving meshes,
the coordinate nodes follow the moving boundaries; i.e. the mesh is
reprocessed at each time step which consumes computational power [20].
Further, the moving meshes are more likely to lose their quality due to the
element inversion after several time steps. This requires a re-meshing that
increases computational time [21]. The coordinate nodes do not change in fixed
meshes, making them preferred [22]. Instead, the boundaries are defined with
an interface capture technique, such as Volume of Fluid [23, 24], Marker and
Cell technique [25], Level-Set technique [26], and Diffuse interface technique
[27]. During the last two decades, Level-Set Method (LSM) has experienced
increased usage in hydrodynamics due to its simplicity in capturing the moving
boundaries [28, 29], and the accuracy of obtained results [30, 31]. Kandasamy
et al. [32] showed the best prediction of the wave-induced separation and the
wave profile on a surface piercing NACA0024 foil by using the LSM with the
Unsteady Reynolds Averaged Navier Stokes (URANS) and Detached Eddy Simulation
(DES). The reason is that the LSM is suitable to simulate the shear layer
instability problems, such as turbulent/splashing free-surface, unlike the
surface tracking method [33]. Further, Chung [34] Used the LSM to capture the
free surface near a flapping flat plate. Moreover, Bergmann and Iollo [35]
showed the applicability of using the LSM in simulating bio-inspired swimmers,
such as fish.

Building a purpose-built code is very beneficial because it offers the
controllability and flexibility to integrate the powerful solvers, FEM and
LSM. This allows for simulations of incompressible flow while capturing the
flapping motion of the foil. Thus, a better understanding of how LEV formation
and separation affect the performance of swing arm flapping foils and their
relation to the energy harvested can be provided by studying the hydrodynamics
and vortex structure around the flapping hydrokinetic energy converter.

## 2\. Computational Model

Foil in simple mode (Figure 1 (b)) moves only in the $y$ direction. However,
in swing-arm mode (Figure 1 (a)), the foil rotates around the base pivoting
point, $P_{B}$, and the radius of rotation is the length of the arm, $R$.The
swing-arm mode increases the movement in $x$ and $y$ directions. Up-stroke and
down-stroke define the motion in the $y$ direction, whereas forward and
backward motions are in the $x$ direction. The hydrofoil moves on an arc
during the cycle. The origin, $O$, is chosen to be at the middle of the line
joining the extreme positions. The maximum distance from $O$ vertically is the
heave amplitude, $h_{o}$, and the maximum horizontal distance is $R-Sh_{o}$,
where $S$ is the swing factor. In addition, the foil performs a second
rotational motion around its pivoting point, $P_{F}$, during its cycle. Hence,
the hydrodynamic and flywheel inertia forces enable the foil to rotate the
connected swing-arm, which rotates the generator shaft at the base pivoting
point, $P_{V}$, to generate power.

### 2.1 Governing Equations

The influence of using swing-arm mode on the LEV, created over an oscillating
foil was modeled using 2D Navier-Stokes equations for unsteady incompressible
laminar flow [36]. The dimensionless governing equations for the flow around
an oscillating foil may be expressed as:

$\nabla\cdot\textbf{V}=0$ (1)
$\frac{k}{2\pi}\frac{\partial\textbf{V}}{\partial
t^{*}}+\left(\textbf{V}\cdot\nabla\right)\textbf{V}=-\nabla
p+\frac{1}{Re}\nabla^{2}\textbf{V}$ (2)

where V is the dimensionless velocity vector of the flow in 2D Cartesian
coordinates with respect to the free-stream velocity, $U$, $t^{*}$ is the
dimensionless time with respect to the flapping period, $T$, $p$ is the
dimensionless pressure of the flow with respect to the dynamic pressure, $\rho
U^{2}$, ($\rho$ is the density of flow), $Re$ is the free-stream Reynolds
number, and $k$ is the non-dimensional reduced frequency parameter, which
results from the non-dimensionalization of time, and is defined as:

$k=\frac{\omega c}{U}$ (3)

where $c$ is the foil’s chord (the characteristic length), and $\omega=2\pi f$
is the angular frequency of the oscillating foil.

### 2.2 Incompressibility Constraint

The steady-state solution is obtained using an iterative process over the
whole domain to satisfy the incompressible continuity equation [37]. In this
study, the pressure stabilization technique is used to restore the coupling
between continuity and momentum equations [38, 39]. This technique modifies
the continuity equation by adding a Laplacian term as follows:

$\nabla\cdot\textbf{V}=\epsilon\nabla^{2}p$ (4)

where $\epsilon$ is the pressure dissipation parameter; its order of magnitude
is the same as the dimensionless time step, $\Delta t^{*}$.

### 2.3 Kinematics

As shown in Figure 2, the flapping foil in swing-arm mode with a harmonic
sinusoidal function experiences simultaneous pitching motion, $\theta(t)$, and
plunge motion in Cartesian coordinates, $x(t)$ and $y(t)$, which are
respectively defined as:

$\theta(t)=\theta_{o}\sin(\omega t)$ (5a) $x(t)=Sh_{o}\left|\sin\left(\omega
t\right)\right|$ (5b) $y(t)=-h_{o}\sin\left(\omega t+\psi\right)$ (5c)

where $\theta_{o}$ is the pitch amplitude, $\psi$ is the phase difference
angle between pitching and plunging motions, $h_{o}$ is the heave amplitude,
and $S$ is the swing factor whose value is in the range of $\left(0,1\right]$
for the swing-arm mode. When $S$ is zero, the motion reverts to a simple
flapping foil. The pitching velocity, $\Omega(t)$, and the upstream velocity
components, $V_{x}(t)$ and $V_{y}(t)$, of a flapping foil due to the swing-arm
mode are respectively expressed as:

$\Omega(t)=\frac{d\theta(t)}{dt}=\theta_{o}\omega\cos\left(\omega t\right)$
(6a) $\begin{split}V_{x}&(t)=\frac{dx(t)}{dt}\\\
&=\left\\{\begin{matrix}Sh_{o}\omega\cos\left(\omega
t\right)&0\leq\frac{t}{T}\leq\frac{1}{2}\\\ -Sh_{o}\omega\cos\left(\omega
t\right)&\frac{1}{2}<\frac{t}{T}\leq 1\end{matrix}\right\\}\end{split}$ (6b)
$V_{y}(t)=\frac{dy(t)}{dt}=-h_{o}\omega\cos\left(\omega t+\psi\right)$ (6c)

Figure 2: Kinematics (left) and foil motion (right) of flapping foil in Swing-
arm mode. The right figure shows the foil position and orientation during the
up-stroke and down-stroke. The arc path indicates the location of the foil’s
pivoting point, $P_{F}$. The extreme positions show the foil at maximum heave
amplitude.

The plunge motion of the hydrofoil introduces an induced angle of attack.
Thus, the effective angle of attack, $\alpha_{eff}$, and the effective
upstream velocity, $V_{eff}$, are:

$\alpha_{eff}(t)=\theta(t)+\arctan\left(\frac{V_{y}}{U+V_{x}}\right)$ (7)
$V_{eff}(t)=\sqrt{\left(U+V_{x}\right)^{2}+\left(V_{y}\right)^{2}}$ (8)

The maximum effective angle of attack, $\alpha_{max,eff}$, which is
approximately the modulus of the angle of attack at quarter-period, highly
affects the generation of peak forces over the hydrofoil and occurrence of
dynamic-stall [40]. The maximum effective upstream velocity, $V_{max,eff}$,
also happens at the same time as the maximum effective angle of attack,
$\alpha_{max,eff}$, and they are defined as:

$\begin{split}\alpha_{max,eff}=\left|\alpha_{eff}\left(\frac{T}{4}\right)\right|=\left|\theta_{o}-\arctan\left(\frac{h_{o}\omega}{U}\right)\right|\end{split}$
(9)
$\begin{split}V_{max,eff}=V_{eff}\left(\frac{T}{4}\right)=\sqrt{\left(U\right)^{2}+\left(h_{o}\omega\right)^{2}}\end{split}$
(10)

The feathering parameter, $\chi$, that indicates power extraction or
propulsion, is defined as:

$\chi=\frac{\theta_{o}}{\arctan\left(\frac{h_{o}\omega}{U}\right)}$ (11)

The necessary but not sufficient condition for power extraction mode is
$\chi>1$, i.e., $\alpha_{eff}\left(\frac{T}{4}\right)>0$. However, the
condition for propulsion mode is $\chi<1$, i.e.,
$\alpha_{eff}\left(\frac{T}{4}\right)<0$.

### 2.4 Extracted Power And Efficiency

For power extraction mode, $\chi>1$, the instantaneous power and thrust
coefficients, $C_{P}$ and $C_{T}$, depend on the contributions from both the
rotational plunge and pitching motions and are respectively defined as:

$C_{P}(t)=\frac{F_{y}V_{y}+F_{x}V_{x}+M\Omega}{0.5\rho cU^{3}}$ (12)
$C_{T}(t)=\frac{F_{x}V_{x}}{0.5\rho cU^{2}}$ (13)

where $F_{x}$, $F_{y}$, and $M$ are the magnitudes of the instantaneous forces
and pitch moment at the foil pivoting point, $P_{F}$, in the Cartesian
coordinates; as a result, the mean of either power or thrust coefficients,
$\overline{C}_{P}$ or $\overline{C}_{T}$, can be calculated over one cycle.
Finally, the efficiency, $\eta$, which indicates the capability of the
oscillating foil system to generate power from the flow, is defined as:

$\eta=\frac{\overline{P}_{extracted}}{P_{available}}=\frac{c}{2h_{o}}\overline{C}_{P}$
(14)

## 3\. Numerical solver

### 3.1 Finite Element Formulation

The numerical model is based on the weighted residual FEM. FEM redefines the
variables by separating their spatial and temporal dependence using the shape
functions, $N_{i}(.)$ [41]. The weighted residual approach minimizes the
weighted error of the solution for the weak form of the PDE [42]. Furthermore,
the residual approach varies by changing the value of the weighting parameter
[43]. Hamada et al. [44] showed that the up-winding Galerkin$/$Least-Square
(GLS), which was introduced by Hughes et al. [45], is the best Weighting
approach to solve second-order differential equations with convection terms,
such as the momentum Navier-Stokes equations. This is because of its stability
and accurate diffusivity compared to the Least-Square, Standard Galerkin,
Collocation, Galerkin$/$Least-Squares, Collocation$/$Galerkin,
Collocation$/$Least-Squares, and Collocation$/$Galerkin$/$Least-Squares. On
the other hand, there are not any convective terms in the modified continuity
equation and its type is elliptic. Therefore, the Standard Galerkin (SG)
formulation is enough to solve the modified continuity equation [46]. A
forward first-order explicit scheme [47] was chosen to conduct the time
marching term. The weight function, $w_{i}(.)$, of SG and GLS are:

$w_{i}(.)=\left\\{\begin{matrix}N_{i}(.)&SG\\\ N_{i}(.)+\tau
L(N_{i}(.))&GLS\end{matrix}\right\\}$ (15)

where $N_{i}(.)$ is the shape functions within the elements, $L(N_{i}(.))$ is
the Least-Square operator over the shape functions which is constructed from
the residual of the governing equation, and $\tau$ is the local mesh
stabilization parameter, which is defined as:

$\tau=\frac{1}{\sqrt{\left(\frac{k}{\pi\Delta
t^{*}}\right)^{2}+\left(\frac{2V_{e}}{L_{e}}\right)^{2}+\left(\frac{4}{L_{e}^{2}Re}\right)^{2}}}$
(16)

where $V_{e}$ and $L_{e}$ are the total velocity and characteristic length of
each element, respectively.

Integration by parts was performed for the Laplacian term to reduce the order
of partial differential equations which enables the usage of bi-linear shape
functions and allows the use of the Neumann boundary condition for the
pressure if needed [41]. Moreover, the mass lumping technique is used for the
unsteady term matrix to reduce computational time [48, 49]. The weak forms of
the continuity and momentum equations in indicial notation are:

$\frac{-1}{\epsilon}\iint N_{i}\nabla\cdot\textbf{V}dA=\iint\nabla
N_{i}\cdot\nabla pdA-\oint N_{i}\frac{\partial p}{\partial n}dS$ (17)
$\begin{split}\frac{k}{2\pi}&\iint N_{i}\frac{\Delta\textbf{V}}{\Delta
t^{*}}dA+\iint\widetilde{N}_{i}\left(\textbf{V}_{e}\cdot\nabla\right)\textbf{V}dA=-\iint
N_{i}\nabla pdA\\\ &-\frac{1}{Re}\left(\iint\left(\nabla
N_{i}\cdot\nabla\right)\textbf{V}dA+\oint
N_{i}\frac{\partial\textbf{V}}{\partial n}dS\right)\end{split}$ (18)

### 3.2 Level-Set Method

LSM is an interface capturing method which defines the nodes of the body in
the fixed meshes using an implicit function, $\phi$ [26]. The function, $\phi$
is constant along the contour levels around the foil’s surface, see Equation
(19). The value of the LSM function classifies the domain into three regions;
inside, on, and outside the foil surface, as seen in Figure 3. Thus, the
boundaries of the foil will be determined with non-positive contour levels.
Further, LSM prevents element inversion and maintains the quality of the mesh,
unlike the dynamic meshes which require re-generation with time [50]. The
boundary of the foil is located dynamically by a Signed-Distance Function
(SDF) at each time step [51], see Figure 4 right.

$\frac{D\phi}{\partial t}=\frac{\partial\phi}{\partial
t}+\bm{V_{LS}}\cdot\nabla\phi=0$ (19)

where $\bm{V_{LS}}$ is the convective velocity vector of the level set
function, $\phi$.

Figure 3: Contour levels of foil using LSM; the levels inside the foil are
negative, the levels outside the foil are positive, and the level on the foil
has a zero value.

### 3.3 Mesh Design

The mesh (Figure 4 left) is a structured square shape, does not include the
foil, and has a length of $20c$, where $c$ is the chord length. Mesh
independence was checked to avoid grid dependent results. Six meshes, with a
number of elements of $120\times 320$, $180\times 480$, $240\times 640$,
$300\times 800$, $360\times 960$, and $420\times 1120$ were selected to
perform the mesh independence, by simulating the oscillating NACA0012 foil in
swing-arm mode with $Re=1\times 10^{5}$, $k=0.08$, and $S=0.25$. Figure 5
shows the instantaneous lift, drag, and power coefficients for the considered
test during a cycle. The results of the meshes $300\times 800$, $360\times
960$, and $420\times 1120$ are very similar, which means that the model
becomes mesh independent, see Table 1. Hence, the change in results, with
finer mesh, can be neglected. Thus, a mesh of $360\times 960$ was selected.

Figure 4: Characteristics of the computational domain (top), showing the mesh, the clustered foil region, and the boundary condition. The flapping foil (bottom) at certain time steps using LSM, indicating the refinement of the mesh at the foil region and the quality of the LSM. Figure 5: Independent mesh test: the lift, drag, and power coefficients for oscillating NACA0012 foil in swing-arm mode using different grid sizes at heave amplitude of $1.25$, pitch amplitude of $70^{o}$, swing factor of $0.25$, non-dimensional reduced frequency of $0.08$, and Reynolds number of $1\times 10^{5}$. Table 1: Independent Mesh test: Absolute percentage error of the mean of power coefficient, negative mean of drag coefficient, and root mean square of lift coefficient during one flapping cycle with respect to the mesh of $420\times 1120$ elements for oscillating NACA0012 foil in swing-arm mode. Number of elements | $|Error|\%$
---|---
$\overline{C_{P}}$ | -$\overline{C_{D}}$ | $RMS\left(C_{L}\right)$
$120\times 320$ | $6.7284$ | $8.4181$ | $7.8392$
$240\times 640$ | $0.5818$ | $4.1344$ | $4.3438$
$360\times 960$ | $0.5466$ | $0.3605$ | $1.3376$

The area of the flapping foil region is $2c\times 6c$ and the number of
elements within that region is $240\times 900$. The distance from the first
node to the foil surface is around $6.67\times e^{-3}c$. The top and bottom
are located $7c$ from the flapping foil region, while the inlet boundary is
located $5c$ from the flapping foil region. All these three boundaries are
specified with a horizontal unity free stream and a zero pressure gradient.
The downstream boundary is $13c$ from the foil region, with a condition of
zero velocity gradient and single-point gauge pressure. The oscillating foil
points have a no-slip condition and zero pressure gradient.

## 4\. Validation Study

The numerical solver was validated using four well-known benchmark problems,
stationary circular cylinder, stationary foil, oscillating circular cylinder,
and flapping foil with simple sinusoidal motion.

### 4.1 Stationary Circular Cylinder

The stationary circular cylinder problem is well-known and often used for
studying the characteristics of the vortex shedding phenomenon, where the flow
behavior is controlled by the Reynolds number [52]. As the Reynolds number
increases, the flow changes from generating symmetric wakes to periodic vortex
shedding behind the circular cylinder. Two simulations of flow past a
stationary circular cylinder were performed to validate the FEM solver using a
body-fitted mesh; two reference cases were chosen for the validation. Huang et
al. [53] showed that the vortex length to cylinder length ratio is about
$0.991$, while the present solver obtained a vortex length to cylinder length
ratio of $0.941$ at $Re=20$. Williamson [54] experimentally obtained a
Stroughal number, $St=\frac{fd}{U}$, where $f$ is the shedding frequency, $d$
is the diameter of the circular cylinder, and $U$ is the free-stream fluid
velocity, of $0.202$ at a Reynolds number of $300$. The Stroughal number using
the present solver was $0.17$. The streamlines around the cylinder at the two
Reynolds numbers are shown in Figure 6, showing the validity of the purpose-
built FEM solver to simulate the problems of incompressible flows.

(a) at $Re$ = 20.

(b) at $Re$ = 300.

Figure 6: Streamlines of flow around a stationary circular cylinder using FEM
at different Reynolds numbers.

### 4.2 Stationary Foil

The developed model’s applicability to accurately determine the incompressible
low-Reynolds number flow over a stationary NACA foil with different angles of
attack is paramount, and therefore used as a validation case. The NACA0012
foil’s location was obtained by implementing the LSM, see Figure 7 (left), and
the vorticity contours of the flow around the foils are shown in the right of
Figure 7. Figure 8 shows the validation of the numerical model by comparing
the variation of lift coefficient and drag coefficients at different angles of
attack with the numerical work of Srinath and Mittal [55]. Furthermore, the
comparison shows the reliability and accuracy of the numerical capturing
technique, LSM, as the maximum absolute error percentage of lift and drag
coefficients is within an acceptable margin, $7\%$.

### 4.3 Oscillating Circular Cylinder

Synchronization or lock-in, between the vortex shedding and vibration
frequencies, is one of the main features of an oscillating circular cylinder,
and the lock-in range becomes larger

Figure 7: Stationary foils points (left) and vorticity contours (right) around
the foils at different angles of attacks, $AOA=0^{o}$, $8^{o}$, $16^{o}$, $\&$
$24^{o}$, and different Reynolds numbers, $Re=100$, $250$, $\&$ $500$.

(a) Drag coefficient.

(b) Lift coefficient.

Figure 8: Variation of hydrodynamic coefficients of a stationary NACA0012 foil
with the angle of attack in an incompressible low-Reynolds number flow.

as the amplitude of oscillation increases [56]. Furthermore, this phenomenon
appears in flapping foil turbines within a certain regime of the reduced
velocity and the pivot point location [57]. To compare with Dekhatawala and
Shah [58], the flow around a circulating cylinder was determined at $Re=185$
with a frequency ratio of 0.8 and amplitude ratio of 0.5. The location of the
oscillating cylinder was obtained using the LSM, see Figure 9 (right).
Dekhatawala and Shah [58] obtained a root mean square of the lift coefficient
of $0.608$ and a mean drag coefficient of $1.5557$, while the results of the
present investigation are $0.6482$ and $1.2758$, respectively. The reasons for
the deviation in the results could be the difference in the solver type used,
the element type in the mesh, and the mesh quality near the cylinder. The
vorticity contours for this case at the upper extreme position are shown in
Figure 9 (left).

Figure 9: Vorticity contours around an oscillating circular cylinder at a
Reynolds number of 185, the frequency ratio of 0.8, and amplitude ratio of 0.5
at the upper extreme position (left). Oscillating cylinder points at different
time steps (right) using LSM.

### 4.4 Flapping-Foil with Simple Sinusoidal Motion

Simple hydrokinetic turbines use sinusoidal flapping foils. LSM was used to
capture the locations of the flapping foil, see the middle of Figure 10, and
the vorticity contours of the flow around the flapping foil during the up- and
down-stroke are shown in the right and left of the same figure, respectively.
Figure 11 shows the power coefficient for a simple flapping NACA0012 foil at
$Re=1\times 10^{4}$, with a Stroughal number of $0.25$, heave amplitude of
$0.8$, and pitch amplitude of $54^{o}$. Lu et al. [59] and Karbasian et al.
[12] obtained an efficiency (Equation 14) of $0.147$ and $0.17$, respectively.
With the present solver, the efficiency was $0.216$. The deviation in results
might be due to the difference in the selected numerical solvers’ type between
the present work (FEM solver) and their work (FVM solvers).

## 5\. Result and Discussion

In this study, a FEM model for the advancement of hydrokinetic flapping foil
turbines was developed. To accurately predict the performance of flapping foil
in the swing-arm correct capturing of the structure of the vortex around the
foil is key. The flow around a NACA0012 foil during the down-stroke phase of
its cycle is shown in Figure 12. The relation between the vortex formation,
the lift, drag forces, and the extracted power are shown in Figure 13. Table 2
shows the used parameters used in this study.

Figure 10: flapping foil (middle) and the Vorticity contours of the flow (right and left) at certain time steps using FEM and LSM at heave amplitude of $0.8$, pitch amplitude of $54^{o}$, Stroughal number of $0.25$, and Reynolds number of $1\times 10^{4}$. Figure 11: Power coefficient as a function of time for a NACA0012 flapping foil with simple sinusoidal motion at heave amplitude of $0.8$, pitch amplitude of $54^{o}$, Stroughal number of $0.25$, and Reynolds number of $1\times 10^{4}$. Table 2: Parameters used in the numerical simulation. Foil’s chord | $c$ | 0.2 m
---|---|---
Heave amplitude | $h_{o}$ | 1.25 m
Pitch amplitude | $\theta_{o}$ | $70^{o}$
Phase difference angle | $\psi$ | $90^{o}$
Swing factor | $S$ | 0.25
Non-dimensional reduced frequency | $k$ | 0.08
Reynolds number | $Re$ | $1\times 10^{5}$

At the start of the down-stroke, the foil has a zero effective angle of
attack, Figure 12 (a), and the power coefficient is at its lowest value,
Figure 13 (a). Here, there is a thin clockwise vortex strip at the leading-
edge, which is created during the end of the up-stroke (the effective angle
goes from positive to zero), and a small trailing-edge vortex (TEV) is
developed, which reduces the lift force [60]. Subsequently, as the effective
angle of attack increases, an adverse pressure gradient appears. This leads to
the creation of the separated shear layer, Figure 12 (b), and a small increase
in the power coefficient, see Figure 13 (b). Figure 12 (c) shows that the flow
reversal layer grows with time, creating the Laminar Separation Bubble (LSB).
The LSB is located at the leading-edge and is thus a small LEV. The LSB is the
main reason for the increase in lift, drag forces and, as a result the
extracted power as shown in Figure 13 (c). As the LSB enlarges, the LEV forms.
Figure 12 (d) shows how the LEV travels towards the trailing edge during its
growth. The stall phenomenon begins when the LEV is well-developed. During
this phase, the flapping foil produces the highest possible lift and extracted
power, Figure 13 (d). The break of LEV into small vortices causes detachment
from the surface of the foil, and sequentially its strength decreases, as
shown in Figure 12 (e). This happens near the maximum value of the effective
angle of attack. Figure 13 (e) shows a slight decrease in the generated
energy. As the value of the effective angle of attack continues to decrease,
the TEV size quickly widens. This pushes the strip of LEV away from the foil’s
surface, as shown in Figures 12 (f) and (g). Here, the lift, drag, and power
coefficients significantly decrease. When the foil is close to finishing the
down-stroke, the TEV moves away from the foil’s surface and the LEV strip gets
close the foil again, as shown in Figures 12 (h) and (i). Here, a small bump
in the extracted power occurs before reaching its lowest value again at a zero
effective angle of attack (end of down-stroke), as shown in Figures 13 (h) and
(i).

Figure 12: Vorticity contours around a NACA0012 flapping foil in swing-arm
mode during the down-stroke phase at Reynolds number of $1\times 10^{5}$,
swing factor of $0.25$, and non-dimensional reduced frequency parameter of
$0.08$.

## 6\. Conclusion

The purpose of this study was to introduce a reliable numerical purpose-built
code for the investigation of the hydrodynamics of hydrokinetic energy
converters. The purpose-built solver is based on Finite Element code with the
Level-Set technique. LSM is used to specify the foil nodes in the fixed grid,
while FEM is applied to obtain the flow field around the flapping foil at each

(a) Lift coefficient.

(b) Drag coefficient.

(c) Power coefficient.

Figure 13: Lift, drag, and power coefficients of a NACA0012 flapping foil in
swing-arm mode during the down-stroke phase at a Reynolds number of $1\times
10^{5}$, swing factor of $0.25$, and non-dimensional reduced frequency
parameter of $0.08$.

time step. LSM reserves the quality of the grid while avoiding re-meshing that
can lead to element inversion. The numerical model was validated using
existing experimental and numerical results. Finally, a detailed vortex
structure study of the flapping foil turbine in swing-arm mode is performed to
show the relationship between the Leading-Edge Vortex and extracted power.

The main indication of the power extraction capability for the flapping foil
is the strength of LEV. This is in line with the findings of Karbasian et al.
[12] and Xu et al. [61]. Further, the gained power starts to fade out when the
stall phenomenon appears. In other words, the flapping foil provides energy as
long as the separation of LEV is delayed. Further numerical investigations,
using FEM and LSM, to study the effect of the shape of the foil or the active-
controlled high lift devices, such as flaps and slats, on power extraction
from flapping foil hydrokinetic turbines in swing-arm mode is recommended.

## References

  * Armaroli and Balzani [2011] Nicola Armaroli and Vincenzo Balzani.  Towards an electricity-powered world.  _Energy & environmental science_, 4(9):3193–3222, 2011.
  * Freris and Infield [2008] Leon Freris and David Infield.  _Renewable energy in power systems_.  John Wiley & Sons, 2008.
  * Limited [2002] The Engineering Business Limited.  Stingray tidal stream energy device - phase 1.  Technical Report ETSU-T/06/00211/00/REP, URN No 02/1400, 2002.
  * Limited [2003] The Engineering Business Limited.  Stingray tidal stream energy device - phase 2.  Technical Report ETSU-T/06/00218/00/REP, URN No 03/1433, 2003.
  * Limited [2005] The Engineering Business Limited.  Stingray tidal stream energy device - phase 3.  Technical Report ETSU-T/06/00230/00/REP, URN No 05/864, 2005.
  * Cardona et al. [2016] Jennifer L Cardona, Michael J Miller, Thomas Derecktor, Steven Winckler, Kenneth Volkmann, Aryssa Medina, Sarah Cowles, Rebecca Lorick, Kenneth S Breuer, and Shreyas Mandre.  Field testing of a 1 kw oscillating hydrofoil energy harvesting system.  In _Marine Energy Technology Symposium_ , 2016.
  * Lago et al. [2010] LI Lago, FL Ponta, and L Chen.  Advances and trends in hydrokinetic turbine systems.  _Energy for sustainable development_ , 14(4):287–296, 2010.
  * Fernandes et al. [2013] Antonio Carlos Fernandes, Mohammadmehdi Armandei, et al.  Marine current energy extraction using torsional galloping based turbine.  In _Offshore Technology Conference_. Offshore technology conference, 2013.
  * Wu et al. [1972] Th Wu et al.  Extraction of flow energy by a wing oscillating in waves.  _Journal of ship research_ , 16(01):66–78, 1972\.
  * Young et al. [2014] John Young, Joseph CS Lai, and Max F Platzer.  A review of progress and challenges in flapping foil power generation.  _Progress in aerospace sciences_ , 67:2–28, 2014.
  * Wu et al. [2015] J Wu, JP Zhan, X Wang, and N Zhao.  Power extraction efficiency improvement of a fully-activated flapping foil: With the help of an auxiliary rotating foil.  _Journal of fluids and structures_ , 57:219–228, 2015.
  * Karbasian et al. [2016] HR Karbasian, Javad Abolfazli Esfahani, and E Barati.  The power extraction by flapping foil hydrokinetic turbine in swing arm mode.  _Renewable energy_ , 88:130–142, 2016.
  * Boudreau et al. [2018] Matthieu Boudreau, Guy Dumas, Mostafa Rahimpour, and Peter Oshkai.  Experimental investigation of the energy extraction by a fully-passive flapping-foil hydrokinetic turbine prototype.  _Journal of Fluids and Structures_ , 82:446–472, 2018.
  * Duarte et al. [2019] Leandro Duarte, Nicolas Dellinger, Guilhem Dellinger, Abdellah Ghenaim, and Abdelali Terfous.  Experimental investigation of the dynamic behaviour of a fully passive flapping foil hydrokinetic turbine.  _Journal of Fluids and Structures_ , 88:1–12, 2019.
  * Martin and Carey [1973] Harold Clifford Martin and Graham F Carey.  _Introduction to finite element analysis: theory and application_.  McGraw-Hill College, 1973.
  * Lewis et al. [2004] Roland W Lewis, Perumal Nithiarasu, and Kankanhalli N Seetharamu.  _Fundamentals of the finite element method for heat and fluid flow_.  John Wiley & Sons, 2004.
  * Donea and Huerta [2003] Jean Donea and Antonio Huerta.  _Finite element methods for flow problems_.  John Wiley & Sons, 2003.
  * Argyris et al. [1979] JH Argyris, H Balmer, J St Doltsinis, PC Dunne, M Haase, M Kleiber, GA Malejannakis, H-P Mlejnek, M Müller, and DW Scharpf.  Finite element method—the natural approach.  _Computer Methods in Applied Mechanics and Engineering_ , 17:1–106, 1979.
  * Mastin et al. [1985] CW Mastin, JF Thompson, and ZUA Warsi.  Numerical grid generation: Foundations and applications, 1985.
  * Liseikin [1999] Vladimir D Liseikin.  _Grid generation methods_ , volume 1.  Springer, 1999.
  * Carey [1997] Graham F Carey.  _Computational grids: generations, adaptation & solution strategies_.  CRC Press, 1997.
  * Edelsbrunner [2001] Herbert Edelsbrunner.  _Geometry and topology for mesh generation_.  Cambridge University Press, 2001.
  * Hirt and Nichols [1981] Cyril W Hirt and Billy D Nichols.  Volume of fluid (vof) method for the dynamics of free boundaries.  _Journal of computational physics_ , 39(1):201–225, 1981.
  * Ashgriz and Poo [1991] Nasser Ashgriz and JY Poo.  Flair: Flux line-segment model for advection and interface reconstruction.  _Journal of computational physics_ , 93(2):449–468, 1991.
  * Harlow and Welch [1965] Francis H Harlow and J Eddie Welch.  Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface.  _The physics of fluids_ , 8(12):2182–2189, 1965\.
  * Osher and Fedkiw [2006] Stanley Osher and Ronald Fedkiw.  _Level set methods and dynamic implicit surfaces_ , volume 153.  Springer Science & Business Media, 2006.
  * Anderson et al. [1998] Daniel M Anderson, Geoffrey B McFadden, and Adam A Wheeler.  Diffuse-interface methods in fluid mechanics.  _Annual review of fluid mechanics_ , 30(1):139–165, 1998.
  * Di Mascio et al. [2007] A Di Mascio, R Broglia, and R Muscari.  On the application of the single-phase level set method to naval hydrodynamic flows.  _Computers & fluids_, 36(5):868–886, 2007.
  * Ali et al. [2021] Karim M Ali, Mohamed Madbouli, Hany M Hamouda, and Amr Guaily.  A stress mapping immersed boundary method for viscous flows.  _Journal of Advanced Research in Fluid Mechanics and Thermal Sciences_ , 87(3):1–20, 2021.
  * Sharma [2015] Atul Sharma.  Level set method for computational multi-fluid dynamics: A review on developments, applications and analysis.  _Sadhana_ , 40(3):627–652, 2015.
  * Elbadry et al. [2021] Yusuf T Elbadry, Amr Gamal Guaily, Mohammed A Boraey, and Mohamed M AbdelRahman.  Active morphing control of airfoil at low reynolds number using level-set method.  In _2021 3rd Novel Intelligent and Leading Emerging Sciences Conference (NILES)_ , pages 265–268. IEEE, 2021.
  * Kandasamy et al. [2005] Manivannan Kandasamy, Tao Xing, Robert Wilson, and Frederick Stern.  Vortical and turbulent structures and instabilities in unsteady free-surface wave induced separation.  In _5th Osaka Colloquium on Advanced Research on Ship Viscous Flow and Hull Form Design by EFD and CFD Approaches, Osaka Prefecture University, Osaka, Japan_ , 2005.
  * Brocchini [2002] M Brocchini.  Free surface boundary conditions at a bubbly/weakly splashing air–water interface.  _Physics of fluids_ , 14(6):1834–1840, 2002.
  * Chung [2016] Meng-Hsuan Chung.  Propulsive performance of a flapping plate near a free surface.  _Journal of Fluids and Structures_ , 65:411–432, 2016.
  * Bergmann and Iollo [2016] Michel Bergmann and Angelo Iollo.  Bioinspired swimming simulations.  _Journal of Computational Physics_ , 323:310–321, 2016\.
  * Wang et al. [2010] Shengyi Wang, Derek B Ingham, Lin Ma, Mohamed Pourkashanian, and Zhi Tao.  Numerical investigations on dynamic stall of low reynolds number flow around oscillating airfoils.  _Computers & fluids_, 39(9):1529–1541, 2010\.
  * Kwak et al. [2005] Dochan Kwak, Cetin Kiris, and Chang Sung Kim.  Computational challenges of viscous incompressible flows.  _Computers & fluids_, 34(3):283–299, 2005.
  * Langtangen et al. [2002] Hans Petter Langtangen, Kent-Andre Mardal, and Ragnar Winther.  Numerical methods for incompressible viscous flow.  _Advances in water resources_ , 25(8-12):1125–1146, 2002.
  * Elhanafy et al. [2017] Ahmed Elhanafy, Amr Guaily, and Ahmed Elsaid.  Pressure stabilized finite elements simulation for steady and unsteady newtonian fluids.  _Journal of applied mathematics and computational mechanics_ , 16(3):17–26, 2017.
  * Simpson et al. [2008] Bradley J Simpson, Franz S Hover, Michael S Triantafyllou, et al.  Experiments in direct energy extraction through flapping foils.  In _The Eighteenth International Offshore and Polar Engineering Conference_. International Society of Offshore and Polar Engineers, 2008.
  * Hughes [2012] Thomas JR Hughes.  _The finite element method: linear static and dynamic finite element analysis_.  Courier Corporation, 2012.
  * Connor and Brebbia [2013] Jerome J Connor and Carlos Alberto Brebbia.  _Finite element techniques for fluid flow_.  Newnes, 2013.
  * Segerlind [1976] Larry J Segerlind.  _Applied finite element analysis_ , volume 316.  Wiley New York, 1976.
  * Hamada et al. [2020] Ahmed A. Hamada, Mahmoud Ayyad, and Amr Guaily.  Variants of the finite element method for the parabolic heat equation: Comparative numerical study.  In Mohamed Hesham Farouk El-Sayed and Maha Amin Hassanein, editors, _Recent Advances in Engineering Mathematics and Physics_. Springer International Publishing, 6 2020.
  * Hughes et al. [1989] Thomas JR Hughes, Leopoldo P Franca, and Gregory M Hulbert.  A new finite element formulation for computational fluid dynamics: Viii. the galerkin/least-squares method for advective-diffusive equations.  _Computer methods in applied mechanics and engineering_ , 73(2):173–189, 1989.
  * Thomée [1984] Vidar Thomée.  _Galerkin finite element methods for parabolic problems_ , volume 1054\.  Springer, 1984.
  * Strikwerda [2004] John C Strikwerda.  _Finite difference schemes and partial differential equations_.  SIAM, 2004.
  * Bicanic and Máca [1991] N Bicanic and J Máca.  Mass lumping by spectral matching.  In _The finite element method in the 1990’s_ , pages 346–356. Springer, 1991.
  * Wendland and Schulz [2005] E Wendland and HE Schulz.  Numerical experiments on mass lumping for the advection-diffusion equation.  _Revista Minerva_ , 2(2):227–233, 2005.
  * Osher and Fedkiw [2001] Stanley Osher and Ronald P Fedkiw.  Level set methods: an overview and some recent results.  _Journal of Computational physics_ , 169(2):463–502, 2001.
  * Osher and Fedkiw [2003] Stanley Osher and Ronald Fedkiw.  Signed distance functions.  In _Level set methods and dynamic implicit surfaces_ , pages 17–22. Springer, 2003.
  * Fredsoe and Sumer [1997] J Fredsoe and MB Sumer.  Hydrodynamics around cylindrical structures, 1997.
  * Huang et al. [2003] Yu-qi Huang, Jian Deng, and An-lu Ren.  Research on lift and drag in unsteady viscous flow around circular cylinders.  _Journal-Zhejiang university engineering science_ , 37(5):596–601, 2003.
  * Williamson [1989] Charles HK Williamson.  Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low reynolds numbers.  _Journal of Fluid Mechanics_ , 206:579–627, 1989.
  * Srinath and Mittal [2009] DN Srinath and S Mittal.  Optimal airfoil shapes for low reynolds number flows.  _International journal for numerical methods in fluids_ , 61(4):355–381, 2009.
  * Koopmann [1967] GH Koopmann.  The vortex wakes of vibrating cylinders at low reynolds numbers.  _Journal of fluid mechanics_ , 28(3):501–512, 1967.
  * Wang et al. [2017] Zhuo Wang, Lin Du, Jisheng Zhao, and Xiaofeng Sun.  Structural response and energy extraction of a fully passive flapping foil.  _Journal of Fluids and Structures_ , 72:96–113, 2017.
  * Dekhatawala and Shah [2019] Ankit Dekhatawala and Rupesh Shah.  Numerical investigation of two-dimensional laminar flow past various oscillating cylinder.  In _Advances in Fluid and Thermal Engineering_ , pages 261–271. Springer, 2019.
  * Lu et al. [2014] Kun Lu, Yonghui Xie, and Di Zhang.  Nonsinusoidal motion effects on energy extraction performance of a flapping foil.  _Renewable energy_ , 64:283–293, 2014.
  * Gharali and Johnson [2013] Kobra Gharali and David A Johnson.  Dynamic stall simulation of a pitching airfoil under unsteady freestream velocity.  _Journal of fluids and structures_ , 42:228–244, 2013.
  * Xu et al. [2019] Jianan Xu, Haiyang Zhu, Daitao Guan, and Yong Zhan.  Numerical analysis of leading-edge vortex effect on tidal current energy extraction performance for chord-wise deformable oscillating hydrofoil.  _Journal of marine science and engineering_ , 7(11):398, 2019.