# Numerical Investigation on Vortex Shedding from a Hydrofoil with a Beveled Trailing Edge

Authors: Seung-Jae Lee; Jun-Hyeok Lee; Jung-Chun Suh
Journal: Modelling and Simulation in Engineering
 (2015)

---

## Abstract

To better understand the vortex shedding mechanism and to assess the capability of our numerical methodology, we conducted numerical investigations of vortex shedding from truncated and oblique trailing edges of a modified NACA 0009 hydrofoil. The hybrid particle-mesh method and the vorticity-based subgrid scale model were employed to simulate these turbulent wake flows. The hybrid particle-mesh method combines the vortex-in-cell and the penalization methods. We have implemented numerical schemes to more efficiently use available computational resources. In this
study, we numerically investigated vortex shedding from various beveled trailing edges at a Reynolds number of 106. We then compared the numerical results with the experimental data, which show good agreement. We also conducted numerical simulations of wakes behind the hydrofoil at rest in periodically varying flows. Results reveal that vortex shedding is affected by the periodicity of a free-stream flow, as well as the trailing-edge shape.


---

## Body

## 1. Introduction

The vortex shedding phenomenon is encountered in many practical engineering applications and physical sciences, and it is an important characteristic of flows past a bluff body. A vortex sheet shed from a solid body consists of alternating vortices of strengthΓ that produce their own velocity field, superimposed on the free-stream velocity. These shedding vortices cause an oscillating force component perpendicular to the direction of flow. This force can induce a vibration in the body. Vortex shedding behind circular and square cylinders is a benchmark problem that has been well investigated and is addressed in a vast amount of literature. Several scholarly reviews [1–4] are entirely devoted to the state of the art of this problem. In contrast, vortex shedding from a streamlined body such as a hydrofoil has been studied to a much lesser extent despite its being of direct relevance to practical engineering problems in hydraulic turbines, pumps, and marine propellers. Propeller-singing is well known as a critical vibration phenomenon generated by the interaction between a Kármán vortex-shedding mechanism from the trailing-edge of a blade and its natural frequencies [5].Some laboratory experiments have been conducted on vortex shedding from the trailing edge of a hydrofoil. Bourgoyne et al. [6, 7] experimentally investigated the dominant features of flows over a two-dimensional (2D) hydrofoil with an antisinging trailing edge at Reynolds numbers ranging from 1
0
6 to 5
×
1
0
7. They concluded that turbulent fluctuations in the near wake are Reynolds-number dependent because of the varying strengths of the structured near-wake vortex shedding. Mulvany et al. [8] numerically simulated turbulent flows around a modified NACA 16 hydrofoil using four turbulent models available in the commercial computational fluid dynamics (CFD) code FLUENT. They compared their numerical results with those from the experimental data of Bourgoyne et al. [6]. The authors mainly discussed the capabilities of the four different models and did not report any vortex shedding phenomenon. Ausoni et al. [9] performed laboratory-scale experiments to investigate the von Kármán vortex street generated in the wake of a 2D hydrofoil with a truncated trailing edge. Their experiments were performed at a zero angle of attack for a Reynolds number ranging from 5.0
×
1
0
5 to 2.0
×
1
0
6, based on the hydrofoil chord length [9, 10]. Their experimental results showed that a lock-in phenomenon occurs, where the vortex shedding frequency is locked onto the natural frequency of a hydrofoil. Interestingly, it has been reported that cavitation has a minor effect on the wake dynamic. Recently, Zobeiri et al. [11] investigated two NACA 0009 hydrofoils with blunt and oblique trailing edges, respectively, for Reynolds numbers ranging from 5.0
×
1
0
5 to 2.9
×
1
0
6 and conducted high-speed visualization and flow-induced vibration measurements. They confirmed experimentally that flow-induced vibration was significantly reduced with an oblique trailing edge compared with a truncated edge. They also concluded that the collision between upper and lower vortices and the resulting vorticity redistribution were the main reasons for the obtained vibration reduction with the oblique trailing edge.This study numerically investigated vortex shedding from truncated and oblique trailing edges of a hydrofoil to better understand the vortex shedding mechanism and to assess the capability of our numerical methodology. We employed the hybrid particle-mesh method and the vorticity-based subgrid scale model to undertake a numerical flow simulation. The hybrid particle-mesh method is a combination of the vortex-in-cell (VIC) method (see [12–15] and references therein) and the penalization technique [16], which has been developed in recent years (see [17–24] for a review). The hybrid particle-mesh method enables the use of fast and efficient techniques for computing differential operators, thereby enabling large scale simulations. However, techniques to save computational memory and CPU time consumption are still required. Thus, we have introduced the implementations of numerical schemes to more efficiently use available computational resources [25, 26]. With these numerical schemes, we numerically investigate vortex shedding from various beveled trailing edges of a NACA 0009 hydrofoil at a Reynolds number of 1
0
6. We then compare our numerical results with experimental data [10, 11] to assess the capability of our numerical methodology. In addition, we investigate the influence of periodically varying free-stream flows on vortex shedding to better understand the vortex shedding mechanism.The organization of the remainder of this paper is as follows: we present a brief description of the numerical methods used for the vortex shedding simulation in Section2. We describe our implementation of numerical schemes to efficiently use available computational resources in Section 3 and the computational procedure followed in Section 4. We present out vortex shedding results from a variety of beveled trailing edges and discuss the influence of periodic free-stream flows in Section 5. In Section 6, we provide a summary of this study and some perspectives for future work.


## 2. Governing Equations and Numerical Methods

The vorticity-velocity formulation of the Navier-Stokes equations allows a purely kinematical problem to be decoupled from the pressure term, which is eliminated by applying the curl operator. Pressure which can be evaluated in an explicit manner using the identified vorticity and velocity fields [27] is not part of the solution algorithm. For a 2D flow parallel to the x
y-plane, the vorticity transport equation in terms of the vorticity can be expressed as(1)
∂
ω
z
∂
t
+
u
·
∇
ω
z
=
ν
∇
2
ω
z
,where ω
z is the scalar plane component of the vorticity vector; that is, ω
≡
ω
z in two dimensions. The evolution of a flow is considered in discrete time steps. The viscous splitting algorithm can be expressed in a Lagrangian framework as(2a)
D
x
D
t
=
u
(2b)
D
ω
z
D
t
=
ν
∇
2
ω
z
,where D
/
D
t
=
∂
/
∂
t
+
u
·
∇ is the material derivative. N
p discrete Lagrangian fluid particles with a finite core size ϵ are linearly superposed to approximate the vorticity field as(3)
ω
z
x
,
t
=
∑
p
=
1
N
p
Γ
p
ζ
ϵ
x
-
x
p
,where ζ
ϵ is a mollification of the Dirac-delta function [12]. Finite core size vortices are used instead of point vortices. Each particle is characterized by its position x
p and its strength Γ
p. For a given vortex particle, the circulation Γ
p is identical to the product of the vorticity and the area of the vortex particle ϵ
2 which also represents the contribution of the vortex particle to the vorticity field, Γ
p
=
∫
‍
ω
z
d
S
≈
ω
z
ϵ
2. The vorticity-carrying fluid particle is advanced with the velocity u and is gradually diffused because of viscous effects.


### 2.1. Vortex-in-Cell (VIC) Method for Convection

Since the continuity equation and the definition of vorticity lead to∇
·
u
=
0 and ∇
×
u
=
ω, the velocity u of (2a) can be expressed as the gradients of the stream-function ψ. The velocity vector is expressed as(4)
u
=
U
∞
+
u
ω
=
U
∞
+
∇
×
ψ
,where U
∞ is the free-stream velocity and u
ω is the rotational velocity. The Poisson equation for the stream-function is ∇
2
ψ
=
-
ω. In two dimensions, the stream function is scalar, ψ
z, and only the vorticity component perpendicular to the plane is nonzero ω
=
ω
z
e
k. The vector Poisson equation reduces to a scalar Poisson equation:(5)
∇
2
ψ
z
=
-
ω
z
.Rotational velocity is defined by u
x
=
∂
ψ
z
/
∂
y and u
y
=
-
∂
ψ
z
/
∂
x. In the VIC method, stream functions to evaluate velocities are computed using a uniform Cartesian mesh. The fast Fourier transform- (FFT-) based Poisson solver reduces the computational cost of obtaining the velocity field to O
(
N
l
o
g
N
) where N is the number of grid points [12]. To compute the stream functions using the FFT-based Poisson solver, the vorticity ω
z is interpolated onto equally spaced Cartesian grid using(6)
ω
z
x
g
=
∑
p
=
1
N
p
ω
z
x
p
W
x
p
-
x
g
h
,where h is the mesh spacing, using the following third order interpolation kernel [28, 29]: (7)
W
x
=
M
4
′
x
=
0
f
o
r
x
>
2
1
2
2
-
x
2
1
-
x
f
o
r
1
<
x
≤
2
1
-
5
2
x
2
+
3
2
x
3
f
o
r
x
≤
1in each coordinate direction. This kernel conserves the 0th, 1st, and 2nd order moments. The subscripts g and p denote the grid and particle quantities, respectively. A vortex particle contributes to the nearest 16 nodes through the M
4
′ scheme, and the total vorticity at each node is obtained by summing the vorticity contributions of all the vortex particles.Boundary conditions for the stream function are required to solve the Poisson equation. If the computational domain boundaries are far enough from the particles, homogeneous Dirichlet boundary conditions (ψ
=
0) may be used. However, using a larger domain is inefficient since it requires a regular grid with too many points. To compactly place a square FFT-domain Ω, we used nonhomogeneous Dirichlet boundary conditions in this study. Unknown stream functions at the domain boundaries (∂
Ω) can be directly obtained from Green’s function solution [27], whereby(8)
ψ
z
x
b
=
-
1
2
π
∑
p
=
1
N
p
Γ
p
ln
⁡
x
b
-
x
p
b
=
1
,
…
,
M
,where x
b
∈
∂
Ω, x
p
∈
Ω, and x
g
≠
x
p. M denotes the number of domain boundary points and the circulation strength of each vortex particle Γ
p
=
ω
z
(
x
p
)
ϵ
2. We note here that all particle vorticity values should be considered. Fast evaluation of the stream functions at the domain boundaries is further discussed in Section 3.3.


### 2.2.  Brinkman Penalization Method for Diffusion

The penalization method was initially designed to take into account solid obstacles in fluid flows [16]. The main point of the penalty term is to replace the usual vorticity creation algorithm for enforcing the no-slip condition for a solid body. By adding the penalization term to (1), the vorticity transport equation becomes(9)
D
ω
z
D
t
=
ν
∇
2
ω
z
+
∇
×
λ
χ
u
s
-
u
,where u
s is the velocity of the solid body and χ denotes a mask function that yields 0 in a fluid and 1 in a solid [16, 18]. λ indicates the penalization parameter with the dimension [
s
-
1
]. To avoid too small time step, the penalization term is evaluated by the implicit expression [16, 18, 30]:(10)
Δ
ω
z
Δ
t
=
∇
×
λ
χ
u
s
n
-
u
~
n
+
1
,where(11)
u
~
n
+
1
=
u
n
+
λ
Δ
t
χ
u
s
n
1
+
λ
Δ
t
χ
.This scheme is unconditionally stable [19]. Equation (9) is rewritten by replacing its terms with those in (10) and (11) as follows:(12)
ω
z
n
+
1
-
ω
z
n
Δ
t
=
ν
∇
2
ω
z
n
+
∇
×
λ
χ
1
+
λ
Δ
t
χ
u
s
n
-
u
n
.The penalization and diffusion terms of (12) are consecutively evaluated on a temporary grid. To reduce any error due to the penalization term, we used a mask function with second-order accuracy, for which the boundary is located at the midpoint of the grid, where the mask function jumps from 0 to 1 [18, 31].


### 2.3.  Smagorinsky Turbulence Model

The filtered vorticity transport equation for a Lagrangian, vorticity-based large-eddy simulation (LES) in two-dimensions can be expressed as [32](13)
D
ω
¯
z
D
t
=
ν
∇
2
ω
¯
z
+
∇
·
Φ
,where the bar denotes spatial filtering at length scale ϵ. The subgrid-scale vorticity stress term Φ, which accounts for the effect of unresolved velocity and vorticity fluctuations, can be defined as Φ
=
ν
T
∇
ω
¯
z. Smagorinsky eddy-viscosity model is one of most commonly used subgrid-scale (SGS) turbulence models. In spite of its simplicity, the Smagorinsky model is still extensively used in practice and has served as a basic tool for the developments of LES modeling [33]. In this model, eddy viscosity is defined by(14)
ν
T
=
C
S
Δ
2
S
¯
,where C
S is a nondimensional Smagorinsky constant and Δ is a subgrid characteristic length scale. S
i
j
=
(
∂
u
i
/
∂
x
j
+
∂
u
j
/
∂
x
i
)
/
2 is the resolved strain-rate tensor and S
¯
=
(
2
S
¯
i
j
S
¯
j
i
)
1
/
2 is its magnitude. Here S
¯
2
=
ω
¯
z
2 in incompressible flows (∇
·
u
=
0), using the Caswell formula [34]. Finally, the vorticity transport equation for the LES in two dimensions is rewritten as (15)
D
ω
¯
z
D
t
=
ν
+
C
S
Δ
2
ω
¯
z
∇
2
ω
¯
z
.In 2D LES models, the filtering process is applied in only two directions, so Δ
=
(
h
x
h
y
)
1
/
2 is defined by the grid spacing. The Smagorinsky constant C
S may vary with the type of flow and the Reynolds number. Values ranging from 0.1 to 0.24 have been suggested in the literatures [33, 35, 36]. The constant C
S is prescribed and its standard value is typically 0.15 [37, 38].


## 2.1. Vortex-in-Cell (VIC) Method for Convection

Since the continuity equation and the definition of vorticity lead to∇
·
u
=
0 and ∇
×
u
=
ω, the velocity u of (2a) can be expressed as the gradients of the stream-function ψ. The velocity vector is expressed as(4)
u
=
U
∞
+
u
ω
=
U
∞
+
∇
×
ψ
,where U
∞ is the free-stream velocity and u
ω is the rotational velocity. The Poisson equation for the stream-function is ∇
2
ψ
=
-
ω. In two dimensions, the stream function is scalar, ψ
z, and only the vorticity component perpendicular to the plane is nonzero ω
=
ω
z
e
k. The vector Poisson equation reduces to a scalar Poisson equation:(5)
∇
2
ψ
z
=
-
ω
z
.Rotational velocity is defined by u
x
=
∂
ψ
z
/
∂
y and u
y
=
-
∂
ψ
z
/
∂
x. In the VIC method, stream functions to evaluate velocities are computed using a uniform Cartesian mesh. The fast Fourier transform- (FFT-) based Poisson solver reduces the computational cost of obtaining the velocity field to O
(
N
l
o
g
N
) where N is the number of grid points [12]. To compute the stream functions using the FFT-based Poisson solver, the vorticity ω
z is interpolated onto equally spaced Cartesian grid using(6)
ω
z
x
g
=
∑
p
=
1
N
p
ω
z
x
p
W
x
p
-
x
g
h
,where h is the mesh spacing, using the following third order interpolation kernel [28, 29]: (7)
W
x
=
M
4
′
x
=
0
f
o
r
x
>
2
1
2
2
-
x
2
1
-
x
f
o
r
1
<
x
≤
2
1
-
5
2
x
2
+
3
2
x
3
f
o
r
x
≤
1in each coordinate direction. This kernel conserves the 0th, 1st, and 2nd order moments. The subscripts g and p denote the grid and particle quantities, respectively. A vortex particle contributes to the nearest 16 nodes through the M
4
′ scheme, and the total vorticity at each node is obtained by summing the vorticity contributions of all the vortex particles.Boundary conditions for the stream function are required to solve the Poisson equation. If the computational domain boundaries are far enough from the particles, homogeneous Dirichlet boundary conditions (ψ
=
0) may be used. However, using a larger domain is inefficient since it requires a regular grid with too many points. To compactly place a square FFT-domain Ω, we used nonhomogeneous Dirichlet boundary conditions in this study. Unknown stream functions at the domain boundaries (∂
Ω) can be directly obtained from Green’s function solution [27], whereby(8)
ψ
z
x
b
=
-
1
2
π
∑
p
=
1
N
p
Γ
p
ln
⁡
x
b
-
x
p
b
=
1
,
…
,
M
,where x
b
∈
∂
Ω, x
p
∈
Ω, and x
g
≠
x
p. M denotes the number of domain boundary points and the circulation strength of each vortex particle Γ
p
=
ω
z
(
x
p
)
ϵ
2. We note here that all particle vorticity values should be considered. Fast evaluation of the stream functions at the domain boundaries is further discussed in Section 3.3.


## 2.2.  Brinkman Penalization Method for Diffusion

The penalization method was initially designed to take into account solid obstacles in fluid flows [16]. The main point of the penalty term is to replace the usual vorticity creation algorithm for enforcing the no-slip condition for a solid body. By adding the penalization term to (1), the vorticity transport equation becomes(9)
D
ω
z
D
t
=
ν
∇
2
ω
z
+
∇
×
λ
χ
u
s
-
u
,where u
s is the velocity of the solid body and χ denotes a mask function that yields 0 in a fluid and 1 in a solid [16, 18]. λ indicates the penalization parameter with the dimension [
s
-
1
]. To avoid too small time step, the penalization term is evaluated by the implicit expression [16, 18, 30]:(10)
Δ
ω
z
Δ
t
=
∇
×
λ
χ
u
s
n
-
u
~
n
+
1
,where(11)
u
~
n
+
1
=
u
n
+
λ
Δ
t
χ
u
s
n
1
+
λ
Δ
t
χ
.This scheme is unconditionally stable [19]. Equation (9) is rewritten by replacing its terms with those in (10) and (11) as follows:(12)
ω
z
n
+
1
-
ω
z
n
Δ
t
=
ν
∇
2
ω
z
n
+
∇
×
λ
χ
1
+
λ
Δ
t
χ
u
s
n
-
u
n
.The penalization and diffusion terms of (12) are consecutively evaluated on a temporary grid. To reduce any error due to the penalization term, we used a mask function with second-order accuracy, for which the boundary is located at the midpoint of the grid, where the mask function jumps from 0 to 1 [18, 31].


## 2.3.  Smagorinsky Turbulence Model

The filtered vorticity transport equation for a Lagrangian, vorticity-based large-eddy simulation (LES) in two-dimensions can be expressed as [32](13)
D
ω
¯
z
D
t
=
ν
∇
2
ω
¯
z
+
∇
·
Φ
,where the bar denotes spatial filtering at length scale ϵ. The subgrid-scale vorticity stress term Φ, which accounts for the effect of unresolved velocity and vorticity fluctuations, can be defined as Φ
=
ν
T
∇
ω
¯
z. Smagorinsky eddy-viscosity model is one of most commonly used subgrid-scale (SGS) turbulence models. In spite of its simplicity, the Smagorinsky model is still extensively used in practice and has served as a basic tool for the developments of LES modeling [33]. In this model, eddy viscosity is defined by(14)
ν
T
=
C
S
Δ
2
S
¯
,where C
S is a nondimensional Smagorinsky constant and Δ is a subgrid characteristic length scale. S
i
j
=
(
∂
u
i
/
∂
x
j
+
∂
u
j
/
∂
x
i
)
/
2 is the resolved strain-rate tensor and S
¯
=
(
2
S
¯
i
j
S
¯
j
i
)
1
/
2 is its magnitude. Here S
¯
2
=
ω
¯
z
2 in incompressible flows (∇
·
u
=
0), using the Caswell formula [34]. Finally, the vorticity transport equation for the LES in two dimensions is rewritten as (15)
D
ω
¯
z
D
t
=
ν
+
C
S
Δ
2
ω
¯
z
∇
2
ω
¯
z
.In 2D LES models, the filtering process is applied in only two directions, so Δ
=
(
h
x
h
y
)
1
/
2 is defined by the grid spacing. The Smagorinsky constant C
S may vary with the type of flow and the Reynolds number. Values ranging from 0.1 to 0.24 have been suggested in the literatures [33, 35, 36]. The constant C
S is prescribed and its standard value is typically 0.15 [37, 38].


## 3. Efficient Implementation of Numerical Methods

### 3.1. Particle-Based Domain Decomposition

Parallel machines and algorithms enable significant reduction in computing time and greater amount of available memory. For high-performance computing, we used the Message Passing Interface (MPI) for our distributed systems, in which each processor has its own local memory. The appropriate decomposition of particles and/or grids should be employed. In domain or spatial decomposition, a given domainΩ is geometrically decomposed by splitting it into several subdomains. Typically, the subdomains are uniform and particles are not evenly distributed among the subdomains. Differences in the number of particles assigned to each processor cause computational load imbalances. As some processors will sit idly while waiting for others, poor performance can result.We introduce a simple idea to achieve a better load balance during the domain decomposition. In vortex particle methods, fluid particles must periodically be redistributed since their accumulation leads to inaccuracies in numerical solutions. This process is called particle remeshing or redistribution. The remeshing step creates a new set of particles on a uniform Cartesian grid, and then the randomly spaced old particles are removed. All the new particles are reindexed byi
=
1,2
,
…
,
N
p from upstream to downstream; that is, higher-indexed particles are located further downstream. Due to the position-dependent index of the particles, the entire spatial domain can be easily split to ensure an equal amount of particles in each subdomain (each processor). As a result, the horizontal extent of the domain is unequal, but the particles are evenly distributed in subdomains of different sizes. In every remeshing step, the entire domain is dynamically partitioned. In this way, each subdomain contains an equal number of particles and better load balance is maintained during simulation. This approach is very simple and effective. No additional algorithms are required to evenly distribute the particles or to identify the position of the particles. As shown in Figure 1, we obtained an approximately 30% reduction in computation time simply by load balancing during parallel computation. This technique is also valuable for the efficient use of available computational memory.Figure 1
Computation times for equal- and unequal-sized subdomains.


### 3.2.  Multidomain Approach

To avoid an excessive domain size, we considered multiple domains. In this study, the entire domainΩ is defined as the union of all the domains, covering all the fluid particles; that is, x
p
∈
Ω where Ω
=
Ω
1
∪
Ω
2
∪
⋯
∪
Ω
N
d. The domain number N
d is not constant and is dependent on a spatially evolving flow. When the domain size exceeds a certain limit, a new domain is created, situated downstream from the domain. As a result, the first domain Ω
1 includes a solid body and the others are sequentially located downstream from the body. For the total circulation conservation, the grid spacing of the ith domain can be defined as h
i
=
ϵ (single-resolution) or h
i
=
2
h
i
-
1 (multi-resolution) where h
1
=
ϵ. Here ϵ denotes the particle size. To interpolate the vorticity values of the particles into the nodes on a grid with different grid spacing, (6) is rewritten as(16)
ω
z
x
g
=
ϵ
2
h
i
2
∑
p
=
1
N
ω
z
x
p
W
x
p
-
x
g
h
i
,where x
p
∈
Ω
i and x
g
∈
Ω
i. The nodes are declared to be a 2D array so that the position of each node can be found straightforwardly from its index. The velocity and vorticity of each domain are sequentially and independently evaluated. Hence, relatively small amounts of working memory are required. This multidomain approach also allows more freedom in the choice of variable arrangements in the domains. The versatility of this approach was demonstrated in our previous study [26]. It can be achieved at the expense of a significant increase in computation time. This issue is further discussed in the next section.


### 3.3.  Approximation of Far-Field Conditions

A substantial part of the overall computation time is spent on the calculation of the boundary stream functions to solve Poisson’s equation for the stream function. This requires the order ofO
(
M
N
p
/
P
) operations to obtain boundary values of domain Ω
i. Here, M denotes the number of nodes at the sides of domain Ω
i and N
p
/
P is the number of particles per processor. To achieve a nearly perfect load balance, particles are evenly distributed among the processors as discussed in Section 3.1 and each processor has all the boundary nodes M. It must also be very useful to avoid all unnecessary communications. N
p
/
P is reduced by increasing the number of processors involved in the computation whereas M is not. This is why the number of boundary nodes M is the most critical for determining the computation time.In this study, we attempted to reduce the number of boundary nodesM using the cubic spline approximation. This is made possible by the smooth variation of the stream functions along the boundary. Only one-fourth of the boundary nodes were chosen as equidistant points for the direct calculations. That is, the stream functions at the selected boundary nodes were directly computed using (8), and the values at boundary nodes not selected were evaluated by interpolation, using cubic splines [39]. In this study, the error between the interpolated value and the true value is typically less than 0.01
%. We have observed that the reduction in computational time is inversely proportional to the number of nodes selected for direct computation. We believe that this approach can be implemented with much lower computational effort because a complex data structure such as a quad-tree used in the fast multipole method [40] is not required. A rigorous comparison of the accuracy and computational costs is planned in the near future.


## 3.1. Particle-Based Domain Decomposition

Parallel machines and algorithms enable significant reduction in computing time and greater amount of available memory. For high-performance computing, we used the Message Passing Interface (MPI) for our distributed systems, in which each processor has its own local memory. The appropriate decomposition of particles and/or grids should be employed. In domain or spatial decomposition, a given domainΩ is geometrically decomposed by splitting it into several subdomains. Typically, the subdomains are uniform and particles are not evenly distributed among the subdomains. Differences in the number of particles assigned to each processor cause computational load imbalances. As some processors will sit idly while waiting for others, poor performance can result.We introduce a simple idea to achieve a better load balance during the domain decomposition. In vortex particle methods, fluid particles must periodically be redistributed since their accumulation leads to inaccuracies in numerical solutions. This process is called particle remeshing or redistribution. The remeshing step creates a new set of particles on a uniform Cartesian grid, and then the randomly spaced old particles are removed. All the new particles are reindexed byi
=
1,2
,
…
,
N
p from upstream to downstream; that is, higher-indexed particles are located further downstream. Due to the position-dependent index of the particles, the entire spatial domain can be easily split to ensure an equal amount of particles in each subdomain (each processor). As a result, the horizontal extent of the domain is unequal, but the particles are evenly distributed in subdomains of different sizes. In every remeshing step, the entire domain is dynamically partitioned. In this way, each subdomain contains an equal number of particles and better load balance is maintained during simulation. This approach is very simple and effective. No additional algorithms are required to evenly distribute the particles or to identify the position of the particles. As shown in Figure 1, we obtained an approximately 30% reduction in computation time simply by load balancing during parallel computation. This technique is also valuable for the efficient use of available computational memory.Figure 1
Computation times for equal- and unequal-sized subdomains.


## 3.2.  Multidomain Approach

To avoid an excessive domain size, we considered multiple domains. In this study, the entire domainΩ is defined as the union of all the domains, covering all the fluid particles; that is, x
p
∈
Ω where Ω
=
Ω
1
∪
Ω
2
∪
⋯
∪
Ω
N
d. The domain number N
d is not constant and is dependent on a spatially evolving flow. When the domain size exceeds a certain limit, a new domain is created, situated downstream from the domain. As a result, the first domain Ω
1 includes a solid body and the others are sequentially located downstream from the body. For the total circulation conservation, the grid spacing of the ith domain can be defined as h
i
=
ϵ (single-resolution) or h
i
=
2
h
i
-
1 (multi-resolution) where h
1
=
ϵ. Here ϵ denotes the particle size. To interpolate the vorticity values of the particles into the nodes on a grid with different grid spacing, (6) is rewritten as(16)
ω
z
x
g
=
ϵ
2
h
i
2
∑
p
=
1
N
ω
z
x
p
W
x
p
-
x
g
h
i
,where x
p
∈
Ω
i and x
g
∈
Ω
i. The nodes are declared to be a 2D array so that the position of each node can be found straightforwardly from its index. The velocity and vorticity of each domain are sequentially and independently evaluated. Hence, relatively small amounts of working memory are required. This multidomain approach also allows more freedom in the choice of variable arrangements in the domains. The versatility of this approach was demonstrated in our previous study [26]. It can be achieved at the expense of a significant increase in computation time. This issue is further discussed in the next section.


## 3.3.  Approximation of Far-Field Conditions

A substantial part of the overall computation time is spent on the calculation of the boundary stream functions to solve Poisson’s equation for the stream function. This requires the order ofO
(
M
N
p
/
P
) operations to obtain boundary values of domain Ω
i. Here, M denotes the number of nodes at the sides of domain Ω
i and N
p
/
P is the number of particles per processor. To achieve a nearly perfect load balance, particles are evenly distributed among the processors as discussed in Section 3.1 and each processor has all the boundary nodes M. It must also be very useful to avoid all unnecessary communications. N
p
/
P is reduced by increasing the number of processors involved in the computation whereas M is not. This is why the number of boundary nodes M is the most critical for determining the computation time.In this study, we attempted to reduce the number of boundary nodesM using the cubic spline approximation. This is made possible by the smooth variation of the stream functions along the boundary. Only one-fourth of the boundary nodes were chosen as equidistant points for the direct calculations. That is, the stream functions at the selected boundary nodes were directly computed using (8), and the values at boundary nodes not selected were evaluated by interpolation, using cubic splines [39]. In this study, the error between the interpolated value and the true value is typically less than 0.01
%. We have observed that the reduction in computational time is inversely proportional to the number of nodes selected for direct computation. We believe that this approach can be implemented with much lower computational effort because a complex data structure such as a quad-tree used in the fast multipole method [40] is not required. A rigorous comparison of the accuracy and computational costs is planned in the near future.


## 4. Computational Procedure

First, the entire simulation domainΩ, covering all the particles, is divided into small domains Ω
i (i
=
1
,
…
,
N
d) as discussed in Section 3.2. Each domain Ω
i is consecutively computed according to the following numerical procedure.(1)
Particle-based domain decomposition: as discussed in Section3.1, domain Ω
i is decomposed into subdomains to ensure that the particles are evenly distributed among the processors. If redistribution has not been done prior to the time step, particles that have been removed from a subdomain are assigned to an adjacent processor without performing domain decomposition.
(2)
Particle-to-grid interpolation: each processor interpolates the vorticityω
z of its own particles into the grid nodes through the M
4
′ interpolation kernel.
(3)
Calculation of the boundary condition: each processor computes the boundary stream functions by direct calculation and interpolation using cubic splines, as discussed in Section3.3.
(4)
Evaluation of the velocity field: the stream function on the grid is computed with the FFT-based Poisson solver, from an open-source library called FFTW (fastest Fourier transform in the west) [41, 42]. Thereafter, the rotational velocities on the grid nodes are computed from the definitions u
x
=
∂
ψ
/
∂
y and u
y
=
-
∂
ψ
/
∂
x using the fourth order central finite difference scheme.
(5)
Evaluation of the vorticity field: each processor evaluates the evolution of the vorticity field in time as follows:
(a)
The penalization term, that is, the second term on the right hand side of (12), is computed with the velocities on the grid. The spatial derivatives of the curl operators are evaluated using a centered difference approximation with a fourth-order error.
(b)
The diffusion term, that is, the first term on the right hand side of (12), is evaluated on the grid using a classical 9-point finite difference method.
(c)
The turbulent viscosity(
C
S
Δ
)
2
ω
¯
z in (15) is approximated by the vorticity strength of the grid. Finally, the vorticity field ω
z
n
+
1 is updated.
(6)
Grid-to-particle interpolation: each processor interpolates the velocity and vorticity on the grid back to its own particle positions through theM
4
′ interpolation kernel.Vortex particles are advanced in time using a midpoint predictor-corrector method as follows: particle positions are predicted by their velocities,x
∗
=
x
n
+
Δ
t
u
(
x
n
,
ω
z
n
), then the particle positions are corrected by using combinations of the predicted and previous positions and velocity values, x
n
+
1
=
x
n
+
0.5
Δ
t
u
(
x
n
,
ω
z
n
)
+
u
(
x
∗
,
ω
z
n
+
1
). Particle redistribution is conducted every three time steps to ensure particle overlap. If a new particle has a small vorticity magnitude, that is, Γ
i
<
0.0001
Γ
m
a
x, it is deleted to avoid having too many particles. Loss of circulation is equally distributed among the remaining particles to ensure the conservation of the total vorticity.


## 5. Numerical Results and Discussions

### 5.1. Model Description

We selected a NACA 0009 hydrofoil with a truncated trailing edge for the numerical simulations. This hydrofoil has the same cross-section as the experimental model used by Ausoni et al. [9], Ausoni [10], and Zobeiri et al. [11]. Ten percent of the original chord c
o was removed from the trailing-edge region of the NACA 0009 hydrofoil. The hydrofoil thickness distribution is written as (17)
0.0
≤
x
c
o
≤
0.5
y
c
o
=
0.1737
x
c
o
1
/
2
-
0.2422
x
c
o
+
0.3046
x
c
o
2
-
0.2657
x
c
o
3
0.5
<
x
c
o
≤
1.0
y
c
o
=
0.0004
+
0.1737
1
-
x
c
o
-
0.1898
1
-
x
c
o
2
+
0.0387
1
-
x
c
o
3
.The hydrofoil geometry is further detailed in Ausoni [10]. The maximum thickness, as normalized by the chord length c, is t
m
a
x
/
c
=
0.1 and the thickness at the trailing edge t
T
E
/
c is 0.0322. The bevel angle β is defined in Figure 2. The trailing edge with β
=
0
∘ denotes that the trailing edge is truncated in the vertical direction.Figure 2
Trailing-edge shape. Note that the points A, B, and C indicate the knuckles of the chamfer and trailing-edge wedge, and the bevel angle is defined asβ.


### 5.2.  Numerical Setup

Numerical parameters can be determined using the stability conditionΔ
t
=
ϵ
2
/
4
ν [18]. The grid convergence was achieved for ϵ
=
0.00030, and Δ
t
=
0.00015 was chosen in order to guarantee numerical accuracy for a Reynolds number of 1
0
6. The grid spacing of the ith domain Ω
i is defined as h
i
=
ϵ
2
i
-
1 to reduce computation time. The penalization parameter λ in (12) is set to be 1
0
8. The hydrofoil at a zero angle of attack is immersed in a Cartesian grid (the first domain) that does not conform to the surface of the body. The first domain, with 4096
×
1024 lattice nodes, is fixed and the other domains are dynamically determined depending on the spatially evolving flow. In this study, numerical simulations were conducted until T
=
10. A new domain is created downstream of a body once the number of horizontal grid points exceeds 2048. Although there is no limit on the number of vertical grid points, it was typically less than or equal to 1024. During simulation, less than 0.8 million particles and four computational domains were created. The entire computation in each case, including start-up and printing times, took approximately 170 hours on 16 Intel Xeon64 3.3 GHz CPUs with 6 GB memory per processor.


### 5.3.  Truncated Trailing-Edge

For an impulsively started flow past a hydrofoil at a chord-based Reynolds number of1
0
6, the wake pattern can be classified into three regimes, based on the time evolutions shown in Figure 3. In the very early stages, the wake remained very symmetrical and consisted primarily of two large vortices in the near wake as shown in Figure 3(a). At that time, the drag force gradually decreased. As time went on, the wake became asymmetrical, as shown in Figure 3(b), and then shed a vortex, as shown in Figure 3(c). After a number of shedding periods, vortices were shed regularly from the symmetric trailing edge in the form of two trains of oppositezsign but equal-strength vortices as shown in Figure 3(d). Such regular vortex shedding induces periodic loading on the structure as shown in Figure 4. It is well known that drag-force oscillation during vortex shedding is much smaller than the lift force. Numerical results show that oscillations in the drag force occurred at twice the vortex-shedding frequency, due to the fact that two vortices were shed from alternate sides during one full period of wake oscillation.Figure 3
Wake patterns of the NACA 0009 hydrofoil with a truncated trailing-edge (β
=
0
∘) for R
e
=
1
0
6.
(a)
Symmetric wake atT
=
0.396
(b)
Symmetry-breaking transition atT
=
0.900
(c)
Nonperiodic vortex shedding atT
=
1.098
(d)
Fully developed periodic shedding atT
=
3.500Figure 4
Time evolution of the lift and drag coefficients for the impulsively started NACA 0009 hydrofoil with a truncated trailing edge (β
=
0
∘) at R
e
=
1
0
6. Regime I denotes the developing regime of a symmetric wake, II a symmetry-breaking transition, and III periodic vortex shedding.From our simulation, the Strouhal number based on the trailing-edge thickness,S
t
t
=
f
s
t
T
E
/
U
∞, is about 0.22 for R
e
=
1
0
6. Ausoni [10] (also refer to Zobeiri et al. [11]) had a similar finding and experimentally gave S
t
t
=
0.25, which may be approximately read from their figure. They used a hydrofoil with a 100 mm chord length and 150 mm span length. The aspect ratio, defined as the ratio of the span length to the chord length, was 1.5. A lower aspect ratio leads to a three-dimensional (3D) effect, even though the hydrofoil is 2D. Ausoni [10] also observed experimentally that the shedding process was not induced in phase along the span. It was therefore speculated in [25] that the slight underestimation of the Strouhal number is due to 3D effects which were absent in our simulation.


### 5.4.  Oblique Trailing-Edges

We conducted numerical simulations to observe vortex shedding from trailing edges with six different bevel angles ofβ
=
0, 10, 20, 40, 60, and 70 degrees. Figure 5 shows the instantaneous vorticity distributions at T
=
10 for R
e
=
1
0
6. At β
=
0
∘, the vortex shedding is regular and periodic. As the bevel angle β increases, vortex shedding becomes increasingly disorganized. Such irregular shedding of vortices demonstrates the periodicity of oscillatory forces. Figure 6 shows the power spectral density function (PSD) of the drag and lift coefficients for different bevel angles. Compared with β
=
0
∘, the amplitude of the lift force oscillation for asymmetric trailing edges is significantly reduced, compared with the symmetric trailing edge at β
=
0
∘. The lift force amplitude is gradually decreased as the bevel angle increases, while the drag force amplitude increases and then decreases to between β
=
1
0
∘ and 6
0
∘. The periodicity of both forces is almost lost at β
=
6
0
∘, and there is no definite periodicity of the induced force components at β
=
7
0
∘.Figure 5
Instantaneous vorticity distributions atT
=
10 for R
e
=
1
0
6.
(a)
Bevel angleβ
=
0
∘
(b)
Bevel angleβ
=
1
0
∘
(c)
Bevel angleβ
=
2
0
∘
(d)
Bevel angleβ
=
4
0
∘
(e)
Bevel angleβ
=
6
0
∘
(f)
Bevel angleβ
=
7
0
∘Figure 6
Power spectral density (PSD) of the lift and drag coefficients.f
∗ denotes the frequency normalized by the chord length and the undisturbed flow velocity.
(a)
Truncated trailing edge,β
=
0
∘
(b)
Beveled trailing edge,β
>
0
∘At a bevel angle of1
0
∘ as shown in Figure 6(b), it is interesting that the dominant peak of the drag force is identical to that of the lift force. From a vibration standpoint, the forcing frequency is important for the prevention of resonance. The distance between the two vortices in Figure 5(b) does not stay the same, but the difference between the shedding vortices between β
=
0
∘ and 1
0
∘ (Figures 5(a) and 5(b)) is not statistically significant. A possible reason for this is the phase difference between the vortices shed from the upper and lower sides. In the wake formation region, the two vortices shed from the symmetric trailing edge at β
=
0
∘ are almost exactly out of phase (180 degrees), whereas at β
=
1
0
∘ there is a phase difference of approximately 30 degrees. This means that before one vortex is fully shed from the body another vortex is created. At higher bevel angles, β
>
2
0
∘, the negative vortices from the upper side of the hydrofoil become stronger than the positive ones from the lower side. The oscillation force is induced mostly by the vortices shed from the acute (upper) side. As a result, the dominant frequency of the drag force oscillation appears to be the same as the vortex shedding frequency. We numerically observed that differences in strength and/or phase between two adjacent vortices cause nonperiodic and irregular vortex shedding.Figure7 shows the shedding frequency of vortices from different beveled trailing edges, compared with experimental results by Zobeiri et al. [11]. The numerically simulated shedding frequencies were approximately 6.8 and 6.4 for β
=
0
∘ and 6
0
∘, respectively. Numerical results show good qualitative agreements with the experimental results, but the numerical simulations underestimated the vortex-shedding frequency by about 10
%. Figure 8 shows the mean value and standard deviation of the velocity components. We compared our numerical results with the experimental results measured by Zobeiri et al. [11]. Thetruncated andoblique trailing edges in Figure 8 indicate bevel angles of β
=
0
∘ and 6
0
∘, respectively. For the truncated trailing edge, the mean velocity profiles are more symmetric than those in the experimental data. The results from numerical simulation are in good agreement with the experimental results. We conclude that the vortex shedding phenomenon from the trailing edge of a 2D hydrofoil can be reasonably simulated using 2D LES which considerably reduces the computational cost.Figure 7
Shedding frequency of vortices from the different beveled trailing edges forR
e
=
1
0
6. The experimental data by Zobeiri et al. [11] were derived from LDV measurements, and the approximate shedding frequency, as read from their figure, is normalized by the chord length and the undisturbed flow velocity.Figure 8
Mean and standard deviation of velocities atx
/
c
=
0.516 for R
e
=
1
0
6. The experimental data by Zobeiri et al. [11] were measured at x
/
c
=
0.52 and 0.528 for truncated (β
=
0
∘) and oblique (β
=
6
0
∘) trailing edges, respectively, and the reported Reynolds number of the flow is about 2
×
1
0
6.
(a)
Normalized mean of streamwise velocity
(b)
Normalized mean of transverse velocity
(c)
Normalized standard deviation of streamwise velocity
(d)
Normalized standard deviation of transverse velocity


### 5.5.  Effect of a Periodically Varying Flow on Vortex Shedding

The propeller is generally located behind the ship’s hull and is subjected to a nonaxisymmetric wake field. Considering a 2D flow, the propeller section effectively experiences periodic variation in flow incidence. We conducted numerical simulations to investigate vortex shedding with respect to the sinusoidal motions of the free-stream flow. The hydrofoil was at rest and the flow became the oscillator under these conditions. In the current study, the angle of attack varied sinusoidally with timet; α
=
α
o
s
i
n
(
2
π
f
∞
t
) where f
∞ was the frequency of the free-stream flow oscillation and its oscillation amplitude was restricted at α
o
=
2
∘. The magnitude of the free-stream velocity was kept constant. The tested model is the NACA 0009 hydrofoil with a beveled trailing edge of β
=
6
0
∘.Figure9 shows the instantaneous vorticity distributions at T
=
10 and the PSD of the drag coefficients at four different frequencies of the free-stream flow, that is, f
∞
=
15, 20, 25, and 30. There are two peaks in the PSD functions; the first peak appears around f
∗
≈
6 and the second one is identical to the frequency of the free-stream flow oscillation. The obvious feature recognized in all the PSD functions, as shown in Figure 9, is that the first peak frequency is nearly identical to the vortex-shedding frequency from the hydrofoil in the uniform free-stream flow. Vortex formation in the near wake is thought to be less affected by small inflow oscillations. Most interestingly, at f
∞
=
20 vortices are regularly shed from the trailing edge, as shown in Figure 9(c). In the PSD functions, drag oscillation induced by vortex shedding is clearly observed. This means that a particular oscillation frequency in the free-stream velocity can cause regular and periodic vortex shedding. This is significant with respect to changing a marine propeller’s angle of attack over time.Figure 9
Effect of the incoming flow frequency on vortex shedding from a beveled trailing edge ofβ
=
6
0
∘ at f
∞
=
15, 20, 25, and 30. The contour indicates the instantaneous vorticity distribution at T
=
10, and the power spectral density (PSD) is calculated from the drag coefficients in the range T
=
[
3,10
].
(a)
Steady flow,f
∞
=
0
(b)
Incoming flow frequencyf
∞
=
15
(c)
Incoming flow frequencyf
∞
=
20
(d)
Incoming flow frequencyf
∞
=
25
(e)
Incoming flow frequencyf
∞
=
30


## 5.1. Model Description

We selected a NACA 0009 hydrofoil with a truncated trailing edge for the numerical simulations. This hydrofoil has the same cross-section as the experimental model used by Ausoni et al. [9], Ausoni [10], and Zobeiri et al. [11]. Ten percent of the original chord c
o was removed from the trailing-edge region of the NACA 0009 hydrofoil. The hydrofoil thickness distribution is written as (17)
0.0
≤
x
c
o
≤
0.5
y
c
o
=
0.1737
x
c
o
1
/
2
-
0.2422
x
c
o
+
0.3046
x
c
o
2
-
0.2657
x
c
o
3
0.5
<
x
c
o
≤
1.0
y
c
o
=
0.0004
+
0.1737
1
-
x
c
o
-
0.1898
1
-
x
c
o
2
+
0.0387
1
-
x
c
o
3
.The hydrofoil geometry is further detailed in Ausoni [10]. The maximum thickness, as normalized by the chord length c, is t
m
a
x
/
c
=
0.1 and the thickness at the trailing edge t
T
E
/
c is 0.0322. The bevel angle β is defined in Figure 2. The trailing edge with β
=
0
∘ denotes that the trailing edge is truncated in the vertical direction.Figure 2
Trailing-edge shape. Note that the points A, B, and C indicate the knuckles of the chamfer and trailing-edge wedge, and the bevel angle is defined asβ.


## 5.2.  Numerical Setup

Numerical parameters can be determined using the stability conditionΔ
t
=
ϵ
2
/
4
ν [18]. The grid convergence was achieved for ϵ
=
0.00030, and Δ
t
=
0.00015 was chosen in order to guarantee numerical accuracy for a Reynolds number of 1
0
6. The grid spacing of the ith domain Ω
i is defined as h
i
=
ϵ
2
i
-
1 to reduce computation time. The penalization parameter λ in (12) is set to be 1
0
8. The hydrofoil at a zero angle of attack is immersed in a Cartesian grid (the first domain) that does not conform to the surface of the body. The first domain, with 4096
×
1024 lattice nodes, is fixed and the other domains are dynamically determined depending on the spatially evolving flow. In this study, numerical simulations were conducted until T
=
10. A new domain is created downstream of a body once the number of horizontal grid points exceeds 2048. Although there is no limit on the number of vertical grid points, it was typically less than or equal to 1024. During simulation, less than 0.8 million particles and four computational domains were created. The entire computation in each case, including start-up and printing times, took approximately 170 hours on 16 Intel Xeon64 3.3 GHz CPUs with 6 GB memory per processor.


## 5.3.  Truncated Trailing-Edge

For an impulsively started flow past a hydrofoil at a chord-based Reynolds number of1
0
6, the wake pattern can be classified into three regimes, based on the time evolutions shown in Figure 3. In the very early stages, the wake remained very symmetrical and consisted primarily of two large vortices in the near wake as shown in Figure 3(a). At that time, the drag force gradually decreased. As time went on, the wake became asymmetrical, as shown in Figure 3(b), and then shed a vortex, as shown in Figure 3(c). After a number of shedding periods, vortices were shed regularly from the symmetric trailing edge in the form of two trains of oppositezsign but equal-strength vortices as shown in Figure 3(d). Such regular vortex shedding induces periodic loading on the structure as shown in Figure 4. It is well known that drag-force oscillation during vortex shedding is much smaller than the lift force. Numerical results show that oscillations in the drag force occurred at twice the vortex-shedding frequency, due to the fact that two vortices were shed from alternate sides during one full period of wake oscillation.Figure 3
Wake patterns of the NACA 0009 hydrofoil with a truncated trailing-edge (β
=
0
∘) for R
e
=
1
0
6.
(a)
Symmetric wake atT
=
0.396
(b)
Symmetry-breaking transition atT
=
0.900
(c)
Nonperiodic vortex shedding atT
=
1.098
(d)
Fully developed periodic shedding atT
=
3.500Figure 4
Time evolution of the lift and drag coefficients for the impulsively started NACA 0009 hydrofoil with a truncated trailing edge (β
=
0
∘) at R
e
=
1
0
6. Regime I denotes the developing regime of a symmetric wake, II a symmetry-breaking transition, and III periodic vortex shedding.From our simulation, the Strouhal number based on the trailing-edge thickness,S
t
t
=
f
s
t
T
E
/
U
∞, is about 0.22 for R
e
=
1
0
6. Ausoni [10] (also refer to Zobeiri et al. [11]) had a similar finding and experimentally gave S
t
t
=
0.25, which may be approximately read from their figure. They used a hydrofoil with a 100 mm chord length and 150 mm span length. The aspect ratio, defined as the ratio of the span length to the chord length, was 1.5. A lower aspect ratio leads to a three-dimensional (3D) effect, even though the hydrofoil is 2D. Ausoni [10] also observed experimentally that the shedding process was not induced in phase along the span. It was therefore speculated in [25] that the slight underestimation of the Strouhal number is due to 3D effects which were absent in our simulation.


## 5.4.  Oblique Trailing-Edges

We conducted numerical simulations to observe vortex shedding from trailing edges with six different bevel angles ofβ
=
0, 10, 20, 40, 60, and 70 degrees. Figure 5 shows the instantaneous vorticity distributions at T
=
10 for R
e
=
1
0
6. At β
=
0
∘, the vortex shedding is regular and periodic. As the bevel angle β increases, vortex shedding becomes increasingly disorganized. Such irregular shedding of vortices demonstrates the periodicity of oscillatory forces. Figure 6 shows the power spectral density function (PSD) of the drag and lift coefficients for different bevel angles. Compared with β
=
0
∘, the amplitude of the lift force oscillation for asymmetric trailing edges is significantly reduced, compared with the symmetric trailing edge at β
=
0
∘. The lift force amplitude is gradually decreased as the bevel angle increases, while the drag force amplitude increases and then decreases to between β
=
1
0
∘ and 6
0
∘. The periodicity of both forces is almost lost at β
=
6
0
∘, and there is no definite periodicity of the induced force components at β
=
7
0
∘.Figure 5
Instantaneous vorticity distributions atT
=
10 for R
e
=
1
0
6.
(a)
Bevel angleβ
=
0
∘
(b)
Bevel angleβ
=
1
0
∘
(c)
Bevel angleβ
=
2
0
∘
(d)
Bevel angleβ
=
4
0
∘
(e)
Bevel angleβ
=
6
0
∘
(f)
Bevel angleβ
=
7
0
∘Figure 6
Power spectral density (PSD) of the lift and drag coefficients.f
∗ denotes the frequency normalized by the chord length and the undisturbed flow velocity.
(a)
Truncated trailing edge,β
=
0
∘
(b)
Beveled trailing edge,β
>
0
∘At a bevel angle of1
0
∘ as shown in Figure 6(b), it is interesting that the dominant peak of the drag force is identical to that of the lift force. From a vibration standpoint, the forcing frequency is important for the prevention of resonance. The distance between the two vortices in Figure 5(b) does not stay the same, but the difference between the shedding vortices between β
=
0
∘ and 1
0
∘ (Figures 5(a) and 5(b)) is not statistically significant. A possible reason for this is the phase difference between the vortices shed from the upper and lower sides. In the wake formation region, the two vortices shed from the symmetric trailing edge at β
=
0
∘ are almost exactly out of phase (180 degrees), whereas at β
=
1
0
∘ there is a phase difference of approximately 30 degrees. This means that before one vortex is fully shed from the body another vortex is created. At higher bevel angles, β
>
2
0
∘, the negative vortices from the upper side of the hydrofoil become stronger than the positive ones from the lower side. The oscillation force is induced mostly by the vortices shed from the acute (upper) side. As a result, the dominant frequency of the drag force oscillation appears to be the same as the vortex shedding frequency. We numerically observed that differences in strength and/or phase between two adjacent vortices cause nonperiodic and irregular vortex shedding.Figure7 shows the shedding frequency of vortices from different beveled trailing edges, compared with experimental results by Zobeiri et al. [11]. The numerically simulated shedding frequencies were approximately 6.8 and 6.4 for β
=
0
∘ and 6
0
∘, respectively. Numerical results show good qualitative agreements with the experimental results, but the numerical simulations underestimated the vortex-shedding frequency by about 10
%. Figure 8 shows the mean value and standard deviation of the velocity components. We compared our numerical results with the experimental results measured by Zobeiri et al. [11]. Thetruncated andoblique trailing edges in Figure 8 indicate bevel angles of β
=
0
∘ and 6
0
∘, respectively. For the truncated trailing edge, the mean velocity profiles are more symmetric than those in the experimental data. The results from numerical simulation are in good agreement with the experimental results. We conclude that the vortex shedding phenomenon from the trailing edge of a 2D hydrofoil can be reasonably simulated using 2D LES which considerably reduces the computational cost.Figure 7
Shedding frequency of vortices from the different beveled trailing edges forR
e
=
1
0
6. The experimental data by Zobeiri et al. [11] were derived from LDV measurements, and the approximate shedding frequency, as read from their figure, is normalized by the chord length and the undisturbed flow velocity.Figure 8
Mean and standard deviation of velocities atx
/
c
=
0.516 for R
e
=
1
0
6. The experimental data by Zobeiri et al. [11] were measured at x
/
c
=
0.52 and 0.528 for truncated (β
=
0
∘) and oblique (β
=
6
0
∘) trailing edges, respectively, and the reported Reynolds number of the flow is about 2
×
1
0
6.
(a)
Normalized mean of streamwise velocity
(b)
Normalized mean of transverse velocity
(c)
Normalized standard deviation of streamwise velocity
(d)
Normalized standard deviation of transverse velocity


## 5.5.  Effect of a Periodically Varying Flow on Vortex Shedding

The propeller is generally located behind the ship’s hull and is subjected to a nonaxisymmetric wake field. Considering a 2D flow, the propeller section effectively experiences periodic variation in flow incidence. We conducted numerical simulations to investigate vortex shedding with respect to the sinusoidal motions of the free-stream flow. The hydrofoil was at rest and the flow became the oscillator under these conditions. In the current study, the angle of attack varied sinusoidally with timet; α
=
α
o
s
i
n
(
2
π
f
∞
t
) where f
∞ was the frequency of the free-stream flow oscillation and its oscillation amplitude was restricted at α
o
=
2
∘. The magnitude of the free-stream velocity was kept constant. The tested model is the NACA 0009 hydrofoil with a beveled trailing edge of β
=
6
0
∘.Figure9 shows the instantaneous vorticity distributions at T
=
10 and the PSD of the drag coefficients at four different frequencies of the free-stream flow, that is, f
∞
=
15, 20, 25, and 30. There are two peaks in the PSD functions; the first peak appears around f
∗
≈
6 and the second one is identical to the frequency of the free-stream flow oscillation. The obvious feature recognized in all the PSD functions, as shown in Figure 9, is that the first peak frequency is nearly identical to the vortex-shedding frequency from the hydrofoil in the uniform free-stream flow. Vortex formation in the near wake is thought to be less affected by small inflow oscillations. Most interestingly, at f
∞
=
20 vortices are regularly shed from the trailing edge, as shown in Figure 9(c). In the PSD functions, drag oscillation induced by vortex shedding is clearly observed. This means that a particular oscillation frequency in the free-stream velocity can cause regular and periodic vortex shedding. This is significant with respect to changing a marine propeller’s angle of attack over time.Figure 9
Effect of the incoming flow frequency on vortex shedding from a beveled trailing edge ofβ
=
6
0
∘ at f
∞
=
15, 20, 25, and 30. The contour indicates the instantaneous vorticity distribution at T
=
10, and the power spectral density (PSD) is calculated from the drag coefficients in the range T
=
[
3,10
].
(a)
Steady flow,f
∞
=
0
(b)
Incoming flow frequencyf
∞
=
15
(c)
Incoming flow frequencyf
∞
=
20
(d)
Incoming flow frequencyf
∞
=
25
(e)
Incoming flow frequencyf
∞
=
30


## 6. Concluding Remarks

We simulated impulsively started flows past a modified NACA 0009 hydrofoil at a Reynolds number of1
0
6, using the hybrid particle-mesh method in conjunction with the vorticity-based subgrid scale model. We have achieved significant savings in both memory usage and computation time by making improvements in the numerical schemes: the particle-based domain decomposition, the use of multiple domains, and the approximation of domain boundary conditions. These techniques are simple and easy to implement but also effective in reducing the computational costs. The particle-based domain decomposition scheme is suitable for a distributed system in which each processor has its own private memory (distributed memory).In this study, the numerical simulation results nearly match the experimental results. This demonstrates that vortex shedding from the trailing edge of a 2D hydrofoil can be reasonably simulated with a 2D LES. It is due to the fact that the flow is characterized by quasi-two dimensionality with vortex shedding from the trailing edge of the hydrofoil. We have considered that 2D vortex shedding simulations should not be dismissed so easily. With two-dimensionality, it is possible to solve the problem much less expensively, since 3D problems require significant computational resources. However, the LES should be 3D since turbulence is inherently 3D. While we accept that more realistic 3D simulation is required physically, a 2D LES can be used as a basis for comparison prior to 3D implementation. Once the capability of a numerical algorithm is evaluated in two dimensions, it can be easily extended into three dimensions. This paper is intended as a guide for researchers and designers who are working on a similar problem: vortex shedding from a hydrofoil at high Reynolds numbers.We found in this study that the vortex formation characteristic in periodically varying inflows remains unchanged. A certain frequency of free-stream flow oscillation tends to make vortex shedding more regular again. This means that the periodicity of flow incidence can stimulate an inherent frequency of vortex shedding from a hydrofoil. However, the effect of inflow oscillation on vortex shedding is not yet fully understood. Future work will include experimental investigations, as well as a 3D flow simulation, to better understand the relationship between inflow oscillation and vortex shedding.



---

*Source: 565417-2015-09-01.xml*