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

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 10. 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.


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][2][3][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 10 6 to 5 × 10 7 .They concluded that turbulent fluctuations in the near wake are Reynoldsnumber 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 × 10 5 to 2.0 × 10 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 × 10 5 to 2.9 × 10 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 vortexin-cell (VIC) method (see [12][13][14][15] and references therein) and the penalization technique [16], which has been developed in recent years (see [17][18][19][20][21][22][23][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 10 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 Section 2. 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 freestream flows in Section 5.In Section 6, we provide a summary of this study and some perspectives for future work.

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 -plane, the vorticity transport equation in terms of the vorticity can be expressed as where   is the scalar plane component of the vorticity vector; that is,  ≡   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 where / = / + u ⋅ ∇ is the material derivative.  discrete Lagrangian fluid particles with a finite core size  are linearly superposed to approximate the vorticity field as 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  and its strength Γ  .For a given vortex particle, the circulation Γ  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, Γ  = ∫    ≈    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 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,   , and only the vorticity component perpendicular to the plane is nonzero  =     .The vector Poisson equation reduces to a scalar Poisson equation: Rotational velocity is defined by   =   / and   = −  /.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(log) where  is the number of grid points [12].To compute the stream functions using the FFT-based Poisson solver, the vorticity   is interpolated onto equally spaced Cartesian grid using where ℎ is the mesh spacing, using the following third order interpolation kernel [28,29]: in each coordinate direction.This kernel conserves the 0th, 1st, and 2nd order moments.The subscripts  and  denote the grid and particle quantities, respectively.A vortex particle contributes to the nearest 16 nodes through the   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 where x  ∈ Ω, x  ∈ Ω, and x  ̸ = x  . denotes the number of domain boundary points and the circulation strength of each vortex particle Γ  =   (x  ) 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.

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 where u  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 [ −1 ].To avoid too small time step, the penalization term is evaluated by the implicit expression [16,18,30]: where This scheme is unconditionally stable [19].Equation ( 9) is rewritten by replacing its terms with those in ( 10) and ( 11) as follows: 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 secondorder accuracy, for which the boundary is located at the midpoint of the grid, where the mask function jumps from 0 to 1 [18,31].

Smagorinsky Turbulence Model.
The filtered vorticity transport equation for a Lagrangian, vorticity-based largeeddy simulation (LES) in two-dimensions can be expressed as [32] 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 Φ = ]  ∇  .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 where   is a nondimensional Smagorinsky constant and Δ is a subgrid characteristic length scale.  = (  /  +   /  )/2 is the resolved strain-rate tensor and || = (2    ) 1/2 is its magnitude.Here || 2 = |  | 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 In 2D LES models, the filtering process is applied in only two directions, so Δ = (ℎ  ℎ  ) 1/2 is defined by the grid spacing.
The Smagorinsky constant   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   is prescribed and its standard value is typically 0.15 [37,38].

Particle-Based Domain Decomposition.
Parallel machines and algorithms enable significant reduction in computing time and greater amount of available memory.For highperformance 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 by  = 1, 2, . . .,   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.

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  ∈ Ω where Ω = Ω 1 ∪ Ω 2 ∪ ⋅ ⋅ ⋅ ∪ Ω   .The domain number   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 th domain can be defined as ℎ  =  (single-resolution) or ℎ  = 2ℎ −1 (multiresolution) where ℎ 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 where x  ∈ Ω  and x  ∈ Ω  .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.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 .It must also be very useful to avoid all unnecessary communications.  / is reduced by increasing the number of processors involved in the computation whereas  is not.This is why the number of boundary nodes  is the most critical for determining the computation time.

Approximation of Far-Field
In this study, we attempted to reduce the number of boundary nodes  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.

Computational Procedure
First, the entire simulation domain Ω, covering all the particles, is divided into small domains Ω  ( = 1, . . .,   ) as discussed in Section 3.2.Each domain Ω  is consecutively computed according to the following numerical procedure.
(1) Particle-based domain decomposition: as discussed in Section 3.1, domain Ω  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   of its own particles into the grid nodes through the   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 Section 3.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   = / and   = −/ 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 (  Δ) 2 |  | in ( 15) is approximated by the vorticity strength of the grid.Finally, the vorticity field  +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 the   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  + Δu(x  ,    ), then the particle positions are corrected by using combinations of the predicted and previous positions and velocity values, Particle redistribution is conducted every three time steps to ensure particle overlap.If a new particle has a small vorticity magnitude, that is, |Γ  | < 0.0001|Γ| max , 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.

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   was removed from the trailing-edge region of the NACA 0009 hydrofoil.The hydrofoil thickness distribution is written as The hydrofoil geometry is further detailed in Ausoni [10].The maximum thickness, as normalized by the chord length , is  max / = 0.1 and the thickness at the trailing edge  TE / 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.Δ = 0.00015 was chosen in order to guarantee numerical accuracy for a Reynolds number of 10 6 .The grid spacing of the th domain Ω  is defined as ℎ  = 2 −1 to reduce computation time.The penalization parameter  in ( 12) is set to be 10 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  = 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.

Truncated Trailing-Edge.
For an impulsively started flow past a hydrofoil at a chord-based Reynolds number of 10 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 vortexshedding frequency, due to the fact that two vortices were shed from alternate sides during one full period of wake oscillation.
From our simulation, the Strouhal number based on the trailing-edge thickness, St  =    TE / ∞ , is about 0.22 for Re = 10 6 .Ausoni [10] (also refer to Zobeiri et al. [11]) had a similar finding and experimentally gave St  = 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.

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  = 10 for Re = 10 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  = 10 ∘ and 60 ∘ .The periodicity of both forces is almost lost at  = 60 ∘ , and there is no definite periodicity of the induced force components at  = 70 ∘ .
At a bevel angle of 10 ∘ 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 10 ∘ (Figures 5(a  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  = 10 ∘ 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,  > 20 ∘ , 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.
Figure 7 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 60 ∘ , respectively.Numerical results show good qualitative agreements with the experimental results, but the numerical simulations underestimated the vortexshedding 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].The truncated and oblique trailing edges in Figure 8 indicate bevel angles of  = 0 ∘ and 60 ∘ , 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.

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 Modelling and Simulation in Engineering  sinusoidally with time ;  =   sin(2 ∞ ) where  ∞ was the frequency of the free-stream flow oscillation and its oscillation amplitude was restricted at   = 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  = 60 ∘ .Figure 9 shows the instantaneous vorticity distributions at  = 10 and the PSD of the drag coefficients at four different frequencies of the free-stream flow, that is,  ∞ = 15, 20, 25, and 30.There are two peaks in the PSD functions; the first peak appears around  * ≈ 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  ∞ = 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.

Concluding Remarks
We simulated impulsively started flows past a modified NACA 0009 hydrofoil at a Reynolds number of 10 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.[11] were measured at / = 0.52 and 0.528 for truncated ( = 0 ∘ ) and oblique ( = 60 ∘ ) trailing edges, respectively, and the reported Reynolds number of the flow is about 2 × 10 6 .
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 of O(  /) operations to obtain boundary values of domain Ω  .Here,  denotes the number of nodes at the sides of domain Ω  and   / is the number of particles per processor.

Figure 2 :
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 .
Fully developed periodic shedding at  = 3.500

Figure 6 :
Figure 6: Power spectral density (PSD) of the lift and drag coefficients. * denotes the frequency normalized by the chord length and the undisturbed flow velocity.

Figure 7 :
Figure7: Shedding frequency of vortices from the different beveled trailing edges for Re = 10 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 9 :
Figure 9: Effect of the incoming flow frequency on vortex shedding from a beveled trailing edge of  = 60 ∘ at  ∞ = 15, 20, 25, and 30.The contour indicates the instantaneous vorticity distribution at  = 10, and the power spectral density (PSD) is calculated from the drag coefficients in the range  = [3, 10].