Computational Fluid Dynamic Simulation of Vertical Axis Hydrokinetic Turbines

Hydrokinetic turbines are one of the technological alternatives to generate and supply electricity for rural communities isolated from the national electrical grid with almost zero emission. These technologies may appear suitable to convert kinetic energy of canal, river, tidal, or ocean water currents into electricity. Nevertheless, they are in an early stage of development; therefore, studying the hydrokinetic system is an active topic of academic research. In order to improve their efficiencies and understand their performance, several works focusing on both experimental and numerical studies have been reported. For the particular case of flow behavior simulation of hydrokinetic turbines with complex geometries, the use of computational fluids dynamics (CFD) nowadays is still suffering from a high computational cost and time; thus, in the first instance, the analysis of the problem is required for defining the computational domain, the mesh characteristics, and the model of turbulence to be used. In this chapter, CFD analysis of a H-Darrieus vertical axis hydrokinetic turbines is carried out for a rated power output of 0.5 kW at a designed water speed of 1.5 m = s, a tip speed ratio of 1.75, a chord length of 0.33 m, a swept area of 0.636 m 2 , 3 blades, and NACA 0025 hydrofoil profile.


Introduction
Nowadays, both the developed and developing countries have been actively promoting renewable energy for their environmental and economic benefits such as the reduction of greenhouse gas emissions, decarbonization, and the diversification of the energy supply in order to reduce the dependence on fossil fuels, especially for the production of electricity, creating economic development and jobs in manufacturing and installation of the systems [1][2][3][4]. Several countries have developed policies and programs to promote the use of renewable energy sources. In general, the policies and programs of developing countries are limited in scope, focusing mainly on regulatory measures, with little attention to other relevant aspects such as research and development, market liberalization, information campaigns, and training [3][4][5][6].
The most common renewable power technologies include solar (photovoltaic and solar thermal systems), wind, biogas, geothermal, biomass, low-impact hydroelectricity, and emerging technologies such as wave and hydrokinetic systems. These technologies have unquestionably started to transform the global energy landscape in a promising way [7,8].
The hydrokinetic systems convert kinetic energy of canal, river, and tidal or ocean water currents into mechanical energy and, consequently, in electrical energy without the use of large civil structures, being the energy provided by the hydrokinetic system more valuable and predictable than that supplied by wind and solar devices [9].
Even though the hydrokinetic systems are still in early stages of development, the referred technology is attractive as renewable energy sources, and it could be a response to the decentralization of energy markets due to the low operating costs, good predict ability, and minimal environmental impact associated [10]. Although the hydrokinetic turbines have relatively small-scale power production, they can be installed as multiunit arrays like wind farms to increase energy extraction. Therefore, these turbines can be used for rural electrification of isolated and river-side communities [10].
The two most common types of hydrokinetic turbines extract energy utilizing horizontal or vertical axis rotor with blade moving through the water [11]. The horizontal and vertical axis turbines have axes parallel and perpendicular to the fluid flow, respectively. Horizontal axis turbines are widely used in tidal energy converters [11,12]. In contrast, the hydrokinetic turbines with vertical axis are generally used for small-scale power generation because they are less expensive and require lower maintenance than horizontal axis turbines [12,13]. On the other hand, the rotor of vertical axis hydrokinetic turbines can rotate regardless of the flow direction, which constitutes an advantage. In turn, horizontal axis hydrokinetic turbines typically reach higher tip speeds, making them more prone to cavitation, which reduce the efficiency and create surface damage [11][12][13]. In spite of vertical axis turbine being not as efficient as the horizontal axis turbines because they exhibit a very low starting torque as well as dynamic stability problems, the interest in implementing vertical axis turbines is increased which drives further research into the development of improved turbine designs [11][12][13]. Gney and Kaygusu reported a detailed comparison of various types of hydrokinetic turbines and concluded that vertical axis turbines are more suitable for the cases where water flow rate is relatively limited [14].
There are two different types of vertical axis turbines: (i) those ones based on the drag force, which are included in the former group and (ii) those ones mainly based on the lift force, corresponding to the second group. Savonius turbines have a drag-type rotor and helical turbine (Gorlov turbine). In turn, Darrieus and Hshaped Darrieus turbines have a lift-type rotor since lift force is the main driving force for these kinds of machines. Generally, drag-based turbines are considered to be less efficient than their lift-based counterparts. The power coefficient of Darrieus hydrokinetic turbines is comparatively higher than that of the Savonius hydrokinetic turbines. However, Darrieus hydrokinetic turbines have poor starting characteristics [15].
Due to the high cost of the harvesting energy from water current associated with the use of these turbines, choosing a turbine with an optimum performance at the selected site is utmost importance. The hydrodynamics and performance of these turbines are governed by several parameters such as (i) the tip-speed ratio (TSR or λ ¼ Rψ=V, which is defined as the ratio between the blade tip speed and the fluid speed (V), being R the turbine radius and ψ the rpm); (ii) the aspect ratio AR ð Þ, defined as the ratio between the hydrokinetic turbine height H ð Þ and its diameter D ð Þ; and (iii) the solidity σ ¼ Bc=2R ð Þ , which refers to the ratio of the product of the blade chord length c ð Þ and the number of blades B ð Þ, the blade profiles, and the chord Reynolds number Re ¼ ρVc=μ ð Þ , where ρ and μ are the water density and viscosity, respectively. In the literature, there are several studies including experimental, numerical (computational fluid dynamics) and theoretical studies, focused on optimizing the hydrokinetic turbine geometric parameters, mainly the blade profiles, using different techniques [8,[16][17][18][19][20][21][22][23]. Figure 1 shows the typical power curves for the most common types of turbine. In Figure 1, it can be observed that horizontal axis hydrokinetic turbines are typically more hydrodynamically efficient and operate at much higher TSR values than vertical axis hydrokinetic turbines. Figure 1 is adapted from the behavior of wind turbines; from the best of the authors' knowledge, there is no overall graph regarding the behavior of such types of turbines considered as hydrokinetic turbines.
The objective of this chapter is to present a methodology for the efficient design and numerical simulation of a H-Darrieus vertical axis hydrokinetic turbine using CFD simulation and the commercial software ANSYS Fluent. The designed hydrokinetic turbine is expected to generate 500 W output power at 1.5 m=s water velocity. The hydrodynamic performance (power coefficient) was analyzed for several TSR values.

Vertical axis hydrokinetic turbine hydrodynamic models
For the design of vertical axis hydrokinetic turbine, several numerical models have been used, each of them with their own strengths and weaknesses, to accurately predict the performance of a hydrokinetic turbine depending on their configuration. The most common models are (i) blade element momentum (BEM) models, (ii) vortex models, and (iii) computational fluid dynamics models [24][25][26][27][28][29][30][31][32]. BEM are analytical models that combine the blade element and the momentum theory in order to study the behavior of the water flow on the blades and related forces. It is a technique that was pioneered by Glauert [24], Strickland [25], and Templin [26]. BEM models can be classified into three submodels: (a) single stream tube model, (b) multiple stream tube model, and (c) double multiple stream tube (DMS) model [28,30], respectively, according to the increasing order of complexity. The hydrokinetic turbine is placed inside a single streamtube, and the blade revolution is translated in an actuator disk when the single stream tube model is used. The water speed in the upstream and downstream sides of the turbine is considered constant. Additionally, the effects outside the streamtube are assumed negligible. It is highlighted that this model does not deliver good accuracy due to the assumptions required to be made. The obtained results in a number of cases provide higher estimate values. In contrast, it has fast processing time compared to other models.
A variation of the single streamtube modeling is the multiple streamtube modeling, which divides the single streamtube in several parallel adjacent streamtubes that are independent from each other and have their own undisrupted, wake, and induced velocities. Although this model is not accurate, the predicted performance is close to field test values, tending to give values a little higher for high TSR with a fast processing time [25,28].
On the other hand, a variation of the multiple streamtube modeling is the double multiple streamtube modeling. This model divides the actuator disk into two half cycles in tandem, representing the upstream and downstream sides of the rotor, respectively. This model has suffered several improvements throughout time and offers a good performance prediction and, however, presents convergence problems for high solidity turbines, giving high power prediction for high TSR and having a high processing time [31][32][33][34][35].
A number of studies have been carried out to simplify numerical models [36][37][38][39][40]. Nevertheless, the use of numerical methods associated with different algorithms of CFD to analyze complex fluid flow applications has grown in recent years due to the rapid development of computer technology. Therefore, the analysis of vertical axis hydrokinetic turbines is also being conducted using CFD simulation [18][19][20][21][22][23].
CFD is found in some commercial software tools such as ANSYS Fluent, FLACS, Phoenix, and COMSOL, among others. The CFD tools transform the governing equations of the fundamental physical principles of the fluid flow in discretized algebraic forms, which are solved to find the flow field values in time and space. The simulation results can be very sensitive to the wide range of computational parameters that must be set by the user; for a typical simulation, the user needs to select the variables of interest, turbulence models, computational domain, computational mesh, boundary conditions, methods of discretization, and the convergence criteria, among other simulation setup parameters [20][21][22][23]. In general, the CFD tools compute the flow field values by the equations of the fundamental physical principles of a fluid flow. The physical aspects of any fluid flow are governed by three principles: mass is conserved, Newton's second law (momentum equation) is fulfilled, and energy is conserved. These principles are expressed in integral equations or partial differential equation (continuity, momentum and energy equation), being the most common form of momentum equation, the Navier-Stokes equations for viscous flows and the Euler equations for inviscid flows [20][21][22][23].
The use of CFD simulation of the turbulent flow surrounding the hydrokinetic turbine is either done by statistical modeling of turbulence based on the Reynoldsaveraged Navier-Stokes (RANS) equations, large eddy simulation (LES), or direct numerical simulation (DNS) [41]. RANS are time-averaged equations of motion for fluid flow; therefore, when the flow is not a stationary one, RANS equation is not the same, due to the simulation must be a time-dependent one. This significantly increases the computational time over conventional RANS. The use of LES and DNS methods for this case requires a careful application and significant computational resources [42,43]. Additionally, since LES and DNS methods require very fine grids for wall-bounded flows, they are still in development stages for high Reynolds number flow; therefore, they are not practical for modeling the flow around the hydrokinetic turbine. An alternative for CFD simulation is the use of unsteady Reynolds-averaged Navier-Stokes (URANS) turbulence models, which are computationally economical and of great use among both industry and academia in comparison with LES models [41,42,[44][45][46][47][48][49][50].
Recently, hybrid unsteady RANS/LES method has been increasingly used for certain classes of simulations, including separated flows; nevertheless, the techniques that combine the near-wall RANS region with the outer, large-eddy simulation region need further development [51].
The specific case of the URANS equations that govern incompressible and isothermal turbulent flow around the turbine blade is given by Eqs. (1) and (2) [52]: where i ¼ 1, 2, 3; j ¼ 1, 2, 3; u and p represent the time-averaged velocity and pressure, respectively; μ is the dynamic viscosity of water; ρ represents the density of water; f i expresses the body force per unit of volume, which may represent the Coriolis and centrifugal contributions; x j denotes the spatial coordinate component; t is the time; and τ ij refers to the Reynolds stress, which is required to be modeled for approximating the above equations. Nowadays, turbulent flows may be divided into four groups ( Figure 2). The first one relies on the Boussinesq assumption according to which Reynolds stress is proportional to the mean strain rate. These models are often called eddy viscosity models. The second group includes the socalled nonlinear eddy viscosity models, which allow for the modeling of the turbulent stress as a nonlinear function of mean velocity gradients. Within this group, turbulent scales are determined by solving transport equations (usually, kþ one other physical quantity). The model is set to mimic response of turbulence to certain important types of strain. The third group involves modeling Reynolds stress transport equations known as a second-moment closure method. The last group concerns the DNS, LES, and the detached eddy simulation (DES). These models rely on a direct resolution of the Navier-Stokes equation for some class of the flows (DNS) or direct resolution of the Navier-Stokes equation for large scales and modeling sub-grid scale (LES, DES). DES is a hybrid model of URANS and LES models, being the former one for modeling the wall region and the latter one for the other regions. The motivation of researchers for utilizing DES is the high cost of LES in the boundary layer region; the strategy of using DES aims at accurately modeling the large vortices in the separated shed wake while keeping the simulation of the near-wall region at a limited cost [51].
Eddy viscosity models use the Boussinesq approximation, which is a constitutive relation to compute the turbulence stress tensor, τ ij , defined in Eqs. (3) and (4): In Eq. (5), k is the turbulent kinetic energy, μ t refers to the dynamic eddy viscosity, ρ is the water density, τ ij denotes the shear stress, δ ij is the Kronecker delta, and u expresses the velocity. Based on dimensional analysis, μ t can be determined from a turbulence time scale (or velocity scale) and a length scale: i ð Þ turbulent kinetic energy k ð Þ, represented by Eq. (6); (ii) turbulence dissipation rate ϵ ð Þ, described by Eq. (7); and (iii) specific dissipation rate ω ð Þ, expressed in Eq. (8): Each turbulence model calculates μ t differently. Within the group of eddy viscosity models, it is possible to identify the following one-and two-equation turbulence models: Spalart-Allmaras, standard k À ϵ model, renormalization group (RNG) k À ϵ model, realizable k À ϵ model, standard k À ω model, and shear stress transport (SST) k À ω model. The choice of turbulence model will depend on several considerations such as the physics encompassed in the flow, the established practice for a specific class of problem, the level of accuracy required, the available computational resources, and the amount of time available for the simulation. To select the most appropriate model for each application, the capabilities and limitations of such as options are needed to be understood [20,43,[52][53][54][55][56][57][58][59][60].
The turbulence models most usually utilized for the CFD analysis of a vertical axis hydrokinetic turbine are RANS approach because LES models require a highly refined mesh and are, therefore, extremely expensive from a computational point of view [53]. For vertical axis hydrokinetic design, the most commonly used turbulence models include the referred models mentioned previously, e.g., k À ϵ model, realizable k À ϵ model, standard k À ω model, and shear stress transport (SST) k À ω model. The SST model was developed to provide an answer to the necessity of models which can handle aeronautical flows with strong adverse pressure gradients and separation of the boundary layers. For this purpose, SST model combines the k À ω model in regions near to solid walls, with k À ϵ model for the free stream flow out of the boundary layers. SST model has been found to give very good results for CFD analysis of horizontal and vertical axis hydrokinetic turbines [20,[54][55][56][57][58][59][60].
For this reason, an SST turbulence model is adopted in the design of this vertical axis hydrokinetic turbine proposed here in the same way used by other authors [22,[61][62][63][64][65].
In general, BEM and vortex models previously described must adopt several assumptions and corrections to account for the full three-dimensional, turbulent flow dynamics around the turbine blades. The use of CFD simulation eliminates the need for many, but not all, of these assumptions. It also has the advantage of resolving the full time-dependent flow field allowing for a better understanding of the flow at the blade wall as well as in the wake of the turbine. Although a much greater computation cost is associated with the use of CFD simulation, it has demonstrated an ability to generate results that can be favorably compared with experimental data [19,58,[66][67][68][69].

Geometrical modeling of a Darrieus hydrokinetic turbine
In the vertical axis category, H-Darrieus turbines are the most common options because their installation is simple, they are less expensive, and they require less maintenance than other types of vertical axis hydrokinetic turbines [17,20]. During the design of Darrieus turbines, a high rotor efficiency is desirable for increased water energy extraction and should be maximized within the limits of affordable production. H-Darrieus turbine consists of two or more hydrofoil blades around a shaft perpendicular to the direction of the water flow. The available hydrokinetic power P ð Þ of this turbine can be computed from the water flow and the turbine dimension, as represented by Eq. (9) [11,12]: where A is the swept area and C p refers to the power coefficient of the turbine. It is widely known that a physical limit exists to quantify the amount of energy that can be extracted, which is independent on the turbine design. The energy extraction is maintained in the flow process through the reduction of kinetic energy and subsequent velocity of the water. The magnitude of energy harnessed is a function of the reduction water speed over the turbine. Total extraction, i.e., 100% extraction, would imply zero final velocity and, therefore, zero flow. The zero-flow scenario cannot be achieved; hence all the water kinetic energy may not be utilized. This principle is widely accepted and indicates that hydrokinetic turbine efficiency cannot exceed 59.3%. This parameter is commonly known as C p , being the maximum C p equal to 0.593 is referred to as the Betz limit. It is important to note that the Betz theory assumes constant linear velocity, inviscid, and without swirl flow. Therefore, any rotational forces such as wake rotation and turbulence caused by drag or vortex shedding (tip losses) will further reduce the maximum efficiency. Efficiency losses are generally reduced by avoiding low tip-speed ratios which increase wake rotation, selecting hydrofoils which have a high lift-to-drag ratio and specialized tip geometries. Making an analysis of data from literature, the value of maximum C p has been found to be usually ranging between 0.15 and 0.586 [22]. This coefficient only considers the mechanical energy converted directly from water energy, i.e., it does not consider the mechanical conversion into electrical energy [70,71].
The power output of the hydrokinetic turbine given by Eq. (9) is also limited by mechanical losses in transmissions and electrical losses [72,73]. Because of these losses and inefficiencies, one more variable is added to Eq. (9). The variable is η, which is a measure of the efficiency of the gearbox, the electrical inverter, and the generator. It takes into account all the friction, slippage, and heat losses associated with the interior mechanical and electrical components. Values for η can greatly differ among the turbine models. Some smaller turbines do not have transmission; therefore, the electric generator is directly moved by the turbine shaft, while larger turbines have transmission [74]. The range of values of η is presented in the literature, highlighting the study conducted by Hagerman et al., which states a range of efficiencies between 95 and 98%. However, for the design of the blade, a reasonable and conservative value of η around 70% was used in this work [75].
Adding η into Eq. (9), Eq. (10) is obtained: Eq. (10) is considered the power equation for the hydrokinetic turbine. P is the net power derived from the water after accounting for losses and inefficiencies. The hydrokinetic turbine parameters considered in the design process mentioned previously are (i) A ð Þ, (ii) P ð Þ and C p , (iii) TSR or λ, (iv) blade profile and c, (v) B, and (vi) σ and AR.
The swept area is the section of water that encloses the turbine in its movement. Its shape depends on the rotor configuration; thus, the swept area of horizontal axis hydrokinetic turbine is a circular shaped, while for a straight-bladed vertical axis hydrokinetic turbine, the swept area has a rectangular shape. Therefore, for a H-Darrieus vertical axis hydrokinetic turbine, the swept area depends on both D ð Þ and H ð Þ [76,77]. Therefore, A is given by Eq. (11): A limits the volume of water passing throughout the turbine. The rotor converts the energy contained in the water in a rotational movement, so as the area is larger, the power output under the same water operating conditions is also larger. Furthermore, C p is strongly dependent on λ, defined in Eq. (12). This parameter is the ratio between the tangential speed at the blade tip and the water speed: Each rotor design has an optimal TSR at which the maximum power extraction is achieved. On the other hand, the number of blades B ð Þ has a direct effect in the smoothness of the rotor operation as they can compensate cycled hydrodynamic loads [76][77][78]. For easiness of the turbine manufacturing, three blades have been considered in the rotor design. In general, three blades are used for the turbine system to keep the dynamic balance and minimize the fatigue effect [75][76][77].
In turn, σ is defined as the ratio between the total blade area and the projected turbine area [76,77]. It is a crucial non-dimensional parameter affecting selfstarting capabilities. For a straight-bladed hydrokinetic turbine, σ is calculated from Eq. (13): The factor AR is defined as the ratio between H and R [76,77]. It can be defined as represented by Eq. (14): The optimum range for the rotor AR for a vertical axis wind turbine was obtained using the double multiple streamtube analytical approach by Ahmadi-Baloutaki et al. [77], and it was found to be 1 < H=R < 4. The hydrokinetic turbine of the mentioned work was designed for a rotor with an AR of 1.5.
In relation to the blade profile, a great number of hydrofoil families and thicknesses have been informed to be suitable for Darrieus turbines. Since it is impossible to analyze all of them, the choices must be narrowed in some way. Symmetrical hydrofoils have been traditionally selected because the energy capturing is approximately symmetrical throughout the turbine axis.
In this sense, the majority of the previously conducted research works on vertical axis hydrokinetic turbines has been focused on straight-bladed vertical hydrokinetic turbines equipped with symmetric hydrofoils (such as NACA fourdigit series of 0012, 0015, 0018, and 0025). In the current work, NACA 0025 hydrofoil profile was selected [66,76,79]. The maximum C p of a turbine with NACA 0025 hydrofoil profile was numerical and experimentally analyzed by Dai and Lam [101]. The results showed that a maximum average C p equal to 25.1% was obtained for a TSR of 1745. Consequently, this C p value was used for the design of the blade.
Nevertheless, before using Eq. (10), V must be determined or assumed. In this research, a V value of 1.5 m=s was assumed since it constitutes an average speed of the great rivers in Colombia. In turn, P, η, and ρ equal to 500 W, 70%, and 997 kg=m 3 at 25°C ð Þwere calculated. R and H were equal to 1.13 m and 0.75 m, respectively. It is important to note that long blades can cause many natural frequency of vibration, which must be avoided during the operation [78,80,81].
Subsequently, given the rotor design parameters (e.g., R, H, hydrofoil profile and water current velocity, among others), the main task during the blade design procedure is to determine c, which can be obtained by equaling the thrust on the rotor, determined from the actuator disk theory, with the thrust obtained from the BEM theory, considering the drag different from zero [72,81]. From the actuator disk theory, the axial thrust (I) on the disk can be written as represented by Eq. (15): where a is the axial induction factor, which is defined as the fractional water velocity decrease between the free stream and the rotor plane. On the other hand, because the Darrieus turbine is primarily a lift force drive turbine, the resultant tangential components of the lift and the drag forces on the hydrofoil provide the driving force to the turbine. The value of a can be determined from Eq. (16) that relates C p with a [82]: Therefore, for C p equal to 0.25 and η of 70%, three values for axial induction can be obtained (1.1916, 0.76, 0.0483). The thrust coefficient for an ideal turbine is equal to zero when a is equal to 1. On the other hand, the momentum theory is no longer valid for axial induction factors greater than 0.5 because V in the far wake would be negative. In practice, as the axial induction factor increases above 0.5, the flow patterns through the hydrokinetic turbine become much more complex than those predicted by the momentum theory. However, when axial induction is equal to a small value, the chord length values can be also smaller, generating a decrease in the surface where the water acts and subsequently complicating the structural integrity of the blade; therefore, a value of axial induction equal to 0.76 was used for the design of the vertical axis hydrokinetic turbine proposed here. Figure 3 shows a basic schematic representation of the hydrodynamic design of a H-Darrieus turbine. A feature of H-Darrieus turbine is that the angle of attack (α) on a blade is a function of the blade azimuth angle (θ). In Figure 3, the tip velocity vector and the lift and drag vectors generated by the rotation of the turbine blade are illustrated. As it can be seen from Figure 3, the relative velocity (w) can be obtained from the tangential and normal velocity components given by Eq. (17): where v i is the induced velocity through the rotor. w can be written in a nondimensional form using free stream velocity. Therefore, Eq. (17) can be expressed as Eq. (18): The v i can be written in terms of a, as described by Eq. (19): Using Eqs. (12) and (19), Eq. (18) can be rewritten as Eq. (20): From the geometry of Figure 3, W can also be represented by Eq. (21): Using Eq. (19), Eq. (21) can be rewritten as Eq. (22): From the geometry of Figure 3, the local angle of attack can be expressed as Eq. (23): Eq. (23) can be put in a non-dimensional form using Eq. (19); therefore, Eq. (24) is obtained: As shown in Eq. (24), α, is affected by λ. In Figure 4, the variation of α as the blade rotates according to λ different values is illustrated. It is noteworthy that α is the angle between the vector of the relative velocity and the direction of c. As the vector of the relative velocity changes, a positive and a negative values in the upstream and downstream region of the rotor, respectively, are observed. The larger λ, the smaller α throughout the blade rotation. The larger the blade tip velocity vector is, the larger the tip speed ratio is.
The normal F N ð Þ and tangential F T ð Þ force for a single blade at a single azimuthal location can be expressed as Eqs. (25) and (26), respectively: where C L is the lift coefficient and C D is the drag coefficient for α. The instantaneous thrust force, which is the force of the water on the turbine experienced by one-blade element in the direction of the water flow can be calculate using Eq. (27): Substituting Eqs. (25) and (26) in Eq. (27), Eq. (28) is obtained, which represents the instantaneous thrust force: By equating the thrust value from Eq. (15) and Eq. (28) and substituting Eq. (22), c can be found using Eq. (29): C L and C D depend on the shape of the blade, α, and Reynolds number under a given operating condition. C L and C D values of 0.5097 and 0.0092 were used. These values were found in the literature for the hydrofoil NACA 0025 for α equal to 5° [79].
Although α is a variable in Eq. (29) and depends on θ and α, in order to calculate c, α was considered as a constant and was equaled to angle that maximizes C L and C D . In this case, α was equal to 5°in order to find the blade c for any θ. Therefore, with the function defined in Eq. (29) and once C L , C D , a, and R were found, the c response was possible to be evaluated at different θ values ( Figure 5). It is observed that c changes with azimuthal position, being a cyclic distribution every 180°with asymptotic values at 0°and 180°. For θ values between 0°and 60°and between 120°and 180°, c is very large, even much larger than R. Therefore, by averaging θ values between 60°and 120°, an optimal c value equal to 0.33 m is achieved.
Using Eq. (13), σ can be found, obtaining a value of 0.66. A lower σ requires less material and subsequently lower effective manufacturing cost. On the other hand, the blade aspect ratio (AB) is defined as the ratio between the blade length and c [76], as expressed by Eq. (30): For blades with smaller AB, the lift-to-drag ratio can be reduced. Any decrease in the lift-to-drag ratio will reduce tangential forces. This would decrease the turbine overall torque and, consequently, the turbine output power. Some researchers have demonstrated an increase in the maximum C p with the blade AB increasing a value of about 15. In spite of this fact, utilizing blades with relatively high AB has some drawbacks. Relatively long blades add a great amount of weight to the turbine, which increases the manufacturing and maintenance costs and creates a need for a more complex bearing. Moreover, such blades are exposed to larger bending moments. In this work, the blade AB parameter was fixed at 3.42.
After finding out the value of c for every section, the next step was to multiply this value by the non-dimensional coordinates of NACA 0025 hydrofoil profile. The values of x and y profile coordinates for each section were exported to parametric 3D design software. From the cross sections of the blade and using the Loft command, a 3D model of the whole blade was produced. The resulting image is shown in Figure 6. In general, the geometric specifications of the turbine are given in   Figure 7 can be followed.

Computational fluid dynamics of a Darrieus hydrokinetic turbine
A computational domain is a region in the space in which the numerical equations of the fluid flow are solved by CFD. In order to solve the physics of the flow field around the blade, dividing the flow domain in a set of small subdomains (rotating and stationary domains) is required, which implies the generation of a mesh of cells, also defined as control volumes. Three fundamental aspects to evaluate the accuracy and the resolution time of a simulation are the geometry and size of the cells coupled to the numerical method used to solve the governing equations. The mesh size must be sufficiently small to provide an accurate numerical approximation, but it cannot be so small that the solution is impractical to be obtained with the available computational resource. Thus, the mesh is usually refined in the regions of interest around the main obstacles affecting the flow [83]. In the particular case of the vertical axis hydrokinetic turbines, these obstacles are the blades.
The mesh can be structured or unstructured. In the first case, the mesh consists of quadrilateral cells in 2D or hexahedral cells in 3D. In turn, the unstructured mesh usually consists of triangles in 2D and tetrahedral in 3D, but cells of the mesh can adopt any form. Structured mesh usually implies shorter time resolution; nevertheless, the unstructured meshes may better represent the geometry. Yasushi [84] presents a discussion about the development of efficient computational analysis using unstructured mesh. On the other hand, Luo and Spiegel [85] propose a method to generate a hybrid mesh (coupling structured and unstructured mesh). In some occasions, the hybrid meshes consist of quadrilateral or hexahedral elements generated in layers near the wall surfaces to capture the boundary layer and of triangular or tetrahedral elements that fill the remainder of the domain. Several concepts related to mesh generation are found in [86], and a detailed discussion about the influence of the mesh in CFD applications is presented by Thompson et al. [83].
The CFD package, ANSYS Fluent version 18.0, was used for all the simulations performed in this study. 2DCFD model with less computational cost than 3D model was utilized to represent the vertical hydrokinetic turbine and the water domain. Based on the review of relevant works [78,[87][88][89], the use of a 2D model has been reported to be enough for revealing the factors that influence the performance of the turbine and the majority of factors affecting the flow physics surrounding the vertical axis hydrokinetic turbine, such as the hydrofoil profile, B, σ, and D. In the 2D numerical study of this work, the effects from supporting arms were not taken into consideration. URANS equations were solved using the SIMPLE algorithm for pressure velocity coupling [34,37,90,91].
The vertical hydrokinetic turbine studied was a three-bladed Darrieus rotor with a NACA 0025 blade profile. c was set at 0.33 m with R equal to 0.75 m. And central axis with a D of 0.025 m was placed in the rotation axis. As previously discussed, the mesh is a critical part of a CFD simulation for engineering purposes. It has to be coarse enough so that the calculation is affordable but also fine enough so that each important physical phenomenon is captured and simulated. In this sense, the domain mesh was created around each hydrofoil, and the surrounding water channel geometry was defined based on studies of the boundary extents. There is an inner circular rotating domain connected to a stationary rectangular domain via a sliding interface boundary condition that conserves both mass and momentum. The domain extents were also selected from a series of sensitivity tests to determine the appropriate distance of the walls, inlet and outlet boundaries from the R rotor. The domain extends three upstream and 8R downstream of the center of the turbine and two 5R laterally to either side of the turbine. The circular rotating domain has an overall diameter of two 8R. Because the computational domain is 2D, the turbine blades were implicitly assumed to be infinitely long.
Unstructured meshes were applied to both the rotor away from the near surface region and the outer grids. Finer meshes were used around the blades and regions in the wake of the blades. Particularly, regions at the leading and trailing edge and in the middle of blade were finely meshed in order to capture the flow field more accurately. The outer mesh was coarsened in regions expanding away from the rotor in order to minimize the central processing unit (CPU) time. The different mesh zones used for the present simulations are illustrated in Figure 8, while various mesh details are shown in Figure 9a and b. A number of simulations were carried out in order to determine how the mesh quality affected the CFD results. In this sense, the torque was calculated for each grid using the Fluent solver. The objective was to select the most appropriate mesh that can guarantee low computational costs and good result accuracy. It is widely known that the region near the blade plays an important role on the hydrokinetic operation, since it has the highest gradient of static pressure and velocity. Additionally, the near wake flow, which can extent up to downstream of the blade, has a great effect on the power [56,92,93]. It is well known that a refinement in the boundary layer and a sensitivity study of y þ are very important, since both of them have an effect on the turbine hydrodynamics. Therefore, a common parameter to identify the subparts of the boundary layer is the dimensionless distance from the wall y þ [94], defined by y þ ¼ Δyu þ =v ð Þ , where Δy is the distance of the first node from the wall, u þ is the wall shear velocity, and v is the kinematic viscosity. In this regard, the near-wall region can be roughly subdivided into three layers: the viscous layer y þ < 5 ð Þ, the buffer layer 5 < y þ < 30 ð Þ , and the fully turbulent layer y þ >30 ð Þ [95]. Therefore, the quality of the mesh was also checked, as well as the y þ values around the blades, which is important for the turbulence modeling.
Symmetrical boundary conditions were defined at the top plane and in the bottom plane. Additionally, a uniform pressure on the outlet boundary was set, and a uniform velocity on the inlet boundary with a magnitude of 1.5 m=s was used for the TSR and the corresponding ψ, as shown in Table 2. k À ω shear stress transport (SST) turbulence model was employed for turbulence modeling since it showed better performance for complex flows including adverse pressure gradients and flow separations like occurs in vertical axis hydrokinetic turbines. The no-slip boundary condition was applied on the turbine wall blades. To simulate the rotation of the rotor, the circular turbine mesh with embedded blades allowed the relative movement to the outer inertial fixed domain. An interface wall was introduced between the fixed and rotating domains. The origin of the reference frame is the center of the rotor. The simulation methods used in this study are similar to the methods used in other numerical studies. Here, the rotational speed of the turbine axis is specified by user input.
To solve the viscous sublayer of the turbulence model used, the values of y þ generally must be less than 1. Subsequently, several meshes were made by increasing the refinement of the computational domain near the blade, achieving a good resolution of the boundary layer. The mesh 4 reached a convergence of the results  with 298,878 nodes and 1.58% error on C p ( Table 3). In this regard, further refinement might not improve the numerical results. The numerical results obtained from mesh 4 were compared with several results available in the literature. The modeling of vertical axis hydrokinetic turbines can be done in steady state or transient modes depending on the objectives of the research and the available computational resources. If computational resources are scarce, relatively simple, steady flow models can be used to model the turbine blades in different azimuthal positions [96]. A more common approach is the transient modeling of the moving blades through the use of URANS approach [19,[97][98][99].
In this work, transient analyses were carried out to characterize the performance of the investigated NACA 0025 hydrofoil profile. Performances are described in terms of P, T, and C p , which are calculated according to Eq. (31). Because the flow over Darrieus turbine is a periodic one, sufficient temporal resolution is necessary to ensure proper unsteady simulation of the vertical axis hydrokinetic turbine and to obtain an independent solution on the time step. Different time step sizes Δt that are equivalent to specific rotational displacements along the azimuth were tested. The chosen time step size was Δt ¼ 0:1°, which properly captures the vortex shedding. Time step convergence was monitored for all conserved variables, and it was observed that acceptable levels of residuals (less than 1 Â 10 À6 ) were attained after six rotations of the hydrokinetic turbine. This means that periodic convergence was also achieved: C p variation over a range of TSR values was studied for the selection of the best TSR, leading to the optimal performance of NACA 0025 hydrofoil profile, which was used in the blade geometrical design. The power versus azimuthal position of different TSRs during rotation is shown in Figure 10. The instantaneous power generated by the turbine is equal to the product of the turbine ψ and T acting on it. A nearly sinusoidal curve was obtained with three positive maxima at each turn and three positive minima (when TSR is equal to 1.75) or negative minima (when TSR is equal to 0.50, 1.00, 2.00); meaning that during a revolution, there are periods of time where the turbine produces torque on the fluid. It is observed that C p increases and decreases approximately until 50 and 120°azimuthal position, respectively, i.e., the main power production occurs between 0 and 120°azimuthal positions for the first blade when TSRs are 1.75 and 2.00. In the figure, it can be seen that the maximum torque for the first blade is achieved at the azimuth angle around 50°. After the peak, the drag begins to increase as the blade enters into a dynamic stall, and the drag starts to be dominant up to a θ of 120°. Then, the second blade repeats the motion of the first blade, and the power production is completed with the same motion of the third blade for one rotation of the turbine. Plot of C p versus TSR  Table 3.
Mesh convergence study.
shows positive values of TSR close to 1.75, meaning that fluid is providing torque to the turbine. Beyond 1.75 TSR, C p is negative, indicating that the turbine, rotating at a constant ψ, exerts torque on the fluid. This can be explained because a high TSR implies a high turbine angular speed and, in such as case, kinetic energy contained in flow is not enough to deliver torque to the turbine and make it to rotate with the same ψ.
In the literature, optimum values of TSR for a Darrieus turbine were reported. For example, Kiho et al. [100], using a Darrieus turbine diameter of 1.6 m, found that the highest efficiency was achieved at 0.56 at 1.1 m=s V and two TSRs. In turn, Torri et al. [101] also found that the peak C p of a three-blade straight vertical axis hydrokinetic turbine was about 0.35 at a TSR of around 2. However, Dai et al. [22] performed numerical and experimental analyses with four sets of a three-bladed rotor where blades were designed with NACA 0025 with four different c values. These four sets of rotors were tested in a tank at different V and R. The highest C p was achieved by the rotor with a c of 162.88 mm, R of 450 mm, and TSR of 1.745. Furthermore, authors investigated the same rotor at different V and found 1.2 m=s as the most effective velocity. They concluded that larger rotors are more efficient, while flow velocity has little effect on the turbine efficiency. Figure 10 shows that there is a high fluctuation with each revolution. The fluctuation decreases when λ is equal to 1.75.
The same tendencies of C p are found in the case of TSR of 1.75 and 2.00, as it can be seen in Figure 10. The variation amplitude of C p is lower when TSR is equal to 1.75. This result confirms that the blade with TSR of 1.75 has better performance than any other TSRs. When TSR is equal to 1.75, the maximum C p is 0.62, and the variation amplitude of C p is near 0.37. On the other hand, when TSR is 2.00, the turbine achieved a maximum C p near 0.61 and an amplitude of only 0.83. Therefore, a turbine with a TSR of 1.75 has advantages over a turbine with a TSR of 2.00 in terms of the reduced fluctuation of its torque curve. C p decreases for lower values of TSR because at low values of TSR, the flow surrounding the blade is separated, implying low lift and high drag. As a result, transferred torque from fluid to turbine decreases. In general, the fluctuations in C p and, therefore, in T can produce a high amount of vibration in the turbine. The effect of these vibrations directly influences the fatigue life of the turbine blades and that of the generated power.  Figure 11 shows pressure contours at θ are equal to 0°and 30°. The contour plot shows that there is an increase of pressure from the upstream side to the downstream side across the blade. Additionally, the pressure in the outer area of the profile is observed to be larger (intrados) than that of the inner zone (extrados). The difference of pressure and velocity causes the overall lift for the turbine. The pressures are negative near the ends of the blade but positive at the exit of the computational domain.
In order to perceive the effect of TSR on the turbine performance, the average C p achieved for different TSRs can be found. C p is not constant because the torque and velocity are not constant in a Darrieus turbine. Hence, the average C p per cycle is calculated as the product of the average values of these terms per cycle. TSR and C p are in a direct relationship when the TSR is between 0.5 and 1.75. However, TSR and C p have an opposite tendency for TSR greater than 1.75. The maximum value of C p was achieved when TSR was 1.75; the value of the average C p was 44.33% per cycle. Similar results were obtained by Dai et al. [22] for NACA 0025 blade profile. Additionally, Lain and Osorio [102] used experimental data of Dai and Lams work [22] and developed numerical models. They performed analysis on CFX solver, DMS model, and Fluent solver and achieved efficiencies of 58.6, 46.3, and 52.8%, respectively, when TSR was 1.745. Results observed by using Fluent solver were quite accurate.

Conclusions
The design and simulation of a vertical axis hydrokinetic turbine were presented in this work. The most common modeling method for vertical axis hydrokinetic turbine is CFD. The numerical simulations allow analyzing many different turbine design parameters, providing an optimal configuration for a given set of design parameters. Before the simulations can be performed, the determination of the optimal grid and time step sizes required must be conducted. Finer grids and smaller time steps might give a more accurate solution, but they increase the computational cost. Thus, finding optimal values is also required. For this purpose, the turbulence model generally used is k-ω shear stress transport (SST) turbulence model.
The TSR is a significant parameter that affects the performance of hydrokinetic turbines. Consequently, the performance of the turbine was investigated with a simplified 2D numerical model. From the 2D model, C p was computed for various TSRs. During a turbine revolution, the blade of the turbine may experience large, as Figure 11. The contours of pressure magnitude for the rotor (a) θ ¼ 0°and (b) θ ¼ 30°.
well as rapid variation, in C p . A C p maximum of 62% was achieved when TSR was equal to 1.75. The value of the average C p was 44.33% per cycle. Though Darrieus turbine is very simple to be constructed, it has some disadvantages when compared to axial turbines since they exhibit a lower power coefficient and a variation in the torque produced within the cycle, leading to periodic loading on the components of the turbine.