Numerical Simulation of Water Flow through a Nano-Hydraulic Turbine of Waterfall-Type by Particle Method

Abstract

This study simulates the flow through an impulse-type small-scale hydraulic turbine utilizing a waterfall of extra-low head. The two-dimensional Moving Particle Semi-implicit (MPS) method is employed for the simulation. The fluid is discretized by particles, and the flow is computed by the Lagrangian calculation for the particle motion. When the distance between the particles discretizing the waterfall of a width 𝐵, 𝑙₀, is set at 𝑙₀/𝐵≤0.084, the flow can be simulated with the sufficiently high spatial resolution, and the rotor performance can also be favorably predicted. The present simulation also successfully analyzes the effect of the rotational speed of rotor on the flow and the turbine performance.

1. Introduction

Hydropower is one of the promising renewable energy resources. It is converted to electric energy through hydraulic turbines. In Japan in 2010, approximately 10% of all supplying electric energy is provided by hydraulic power generations. The percentage is expected to increase steadily on the basis of a government policy promoting the development of renewable energy. As large-scale hydroelectric plants require huge dams and long conduits, the places for the construction are hardly remained. Thus, expectations for the development of a small-scale hydropower, of which output is less than 1000 kW, have been increasing. When natural disaster occurs, the large-scale centralized hydraulic power plants may lose the power supplying ability due to the collapse of the power grid. Since the small-scale hydropower, existing in small-scale rivers, irrigation canals, and industrial drainages, distributes widely in Japan, it realizes the small-scale distributed power generation. The small-scale hydropower can contribute the local production for local consumption of electric power, which is more resistant to disaster. Consequently, the development is also of great worth from the viewpoint of constructing a disaster-resistant society.

To exploit effectively the small-scale hydropower, various types of hydraulic turbine have been presented [1–5]. Ikeda et al. [6] developed an impulse-type hydraulic turbine utilizing waterfalls of extra-low head which is 2 m or less in small rivers and agricultural canals and so forth. The nano-hydraulic turbine, of which output power is less than several kW, can be easily carried to the places where they are necessary, and produces electric power easily without damaging the environment. Ikeda et al. [6] investigated the power characteristic of the turbine through a laboratory experiment and made clear the flow inside the rotor by the experimental visualization.

This study proposes the simulation method for the flow in the nano-hydraulic turbine of Ikeda et al. [6], which can be favorably employed for the prediction of the turbine performance. The turbine performance is generally affected by the geometric conditions, such as the blade shape, the blade installation angle, and the number of blades. It is also influenced by the position of turbine relative to the water flow. The numerical simulation can examine the effect of each condition individually and promises to yield design guidelines for high-performance turbine. The nano-hydraulic turbine is driven by a waterfall. The falling water collides with the rotating blades, and it is scattered around the rotor. Therefore, the flow through the turbine includes free surfaces, and it is very complicated. The simulation of this study is based on the Moving Particle Semi-implicit (MPS) method [7–9], which is one of the particle methods for free-surface flows. The fluid is discretized with particles, and the flow is simulated by the Lagrangian computation of the particle motion. There are few numerical simulations of the flow through small-scale hydraulic turbines utilizing small-scale hydropower, except for a MPS simulation of an impulsive turbine driven by a nozzle jet [5]. But the parameters for the MPS method, such as the distance between the particles, were not examined, and accordingly the sufficient knowledge on the simulation of small-scale hydraulic turbines was not obtained. First, this study elucidates the appropriate distance between the particles for the flow and performance simulations of the hydraulic turbine. Secondly, it is demonstrated that the present simulation can analyze the effect of the rotational speed of rotor on the turbine performance.

2. Basic Equations and Simulation Method

2.1. Conservation Equations for Flow

If the flow through the hydraulic turbine is incompressible, it is governed by the mass and momentum conservation equations: D𝜌D𝑡=0,D𝐮1D𝑡=−𝜌∇𝑝+𝜈∇²𝐮+𝐅, (1)(2) where 𝜌 is the density, 𝑡 is the time, 𝐮 is the velocity, 𝑝 is the pressure, 𝜈 is the kinematic viscosity, and 𝐅 is the external forces such as the gravitational force and the surface tension.

This study analyzes (1) and (2) by the Moving Particle Semi-implicit (MPS) method [7–9], which is one of the particle methods. In the MPS method, the fluid is discretized by particles, and the particle motion is computed by the Lagrangian method. Equations (1) and (2) are discretized through the interactions between the particles.

2.2. Particle Interaction Model

The interactions between the particles are modeled with a weight function 𝑤 defined by the following equation: 𝑤(𝑟)=𝑟_𝑒𝑟−1𝑟≤𝑟_𝑒,0𝑟>𝑟_𝑒, (3) where 𝑟 is the distance between two particles and 𝑟_𝑒 stands for a kernel size.

The particle number density at the position of the 𝑖th particle, ⟨𝑛⟩_𝑖, is defined as ⟨𝑛⟩_𝑖=_𝑗≠𝑖𝑤𝐫_𝑗−𝐫_𝑖, (4) where 𝐫_𝑖 is the position vector of the 𝑖th particle.

As the fluid density remains unaltered in the incompressible flow, the particle number density is required to be constant. The incompressible flow condition in the MPS method is satisfied by maintaining ⟨𝑛⟩_𝑖 at a constant value 𝑛⁰.

The Laplacian operator, expressing the viscous term on the right-hand side of (2), is modeled with the weight function. The Laplacian operator at the position of the 𝑖th particle is given as ∇²𝜙_𝑖=2𝑑𝑛⁰Λ_𝑗≠𝑖𝜙_𝑗−𝜙_𝑖𝑤𝐫_𝑗−𝐫_𝑖, (5) where 𝜙 is a physical quantity. A parameter Λ is introduced so that the variance increase is equal to the analytical solution: ∫Λ=^𝑉𝑤(𝑟)𝑟²d𝑣∫_𝑉𝑤(𝑟).d𝑣 (6)

The gradient operator, expressing the pressure gradient term on the right-hand side of (2), is also modeled. The gradient operator at the position of the 𝑖th particle is modeled by setting the interparticle force at the repulsion to ensure the numerical stability [7–9]. Thus, the interacting pressure forces between two particles are not antisymmetric, and the momentum is not always conserved. To resolve this problem, this simulation employs the following model presented by Khayyer and Gotoh [10]: ∇𝜙_𝑖=𝑑𝑛⁰_𝑗≠𝑖𝜙_𝑖+𝜙_𝑗−𝜙_𝑖+𝜙_𝑗𝐫_𝑗−𝐫_𝑖²𝐫_𝑗−𝐫_𝑖𝑤𝐫_𝑗−𝐫_𝑖, (7) where 𝑑 is the number of space dimensions and 𝜙_𝑖 is defined as 𝜙_𝑖=min_𝑗∈𝐽𝜙_𝑖,𝜙_𝑗,𝐽=𝑗∶𝑤𝐫_𝑗−𝐫_𝑖≠0. (8)

2.3. Simulation Method

Equations (1) and (2) are solved by a semi-implicit method, which is used in the SMAC method [11]. If the particle velocity 𝐮^𝑘_𝑖 and position 𝐫^𝑘_𝑖 at time 𝑡=𝑘Δ𝑡 are known, the flow at time 𝑡=(𝑘+1)Δ𝑡 is simulated by the following two steps.

In the first step, the temporal velocity and position for the particle, 𝐮^∗_𝑖 and 𝐫^∗_𝑖, respectively, are calculated from (2) without considering the pressure gradient term. Then, the temporal particle number density ⟨𝑛^∗⟩_𝑖 is computed by using (4).

In the second step, the following Poisson equation is solved for the pressure 𝑝^𝑘+1 so that the mass conservation is satisfied or ⟨𝑛^∗⟩_𝑖 is made to coincide with 𝑛⁰: ∇²𝑝^𝑘+1_𝑖𝜌=−Δ𝑡²𝑛^∗_𝑖−𝑛⁰𝑛⁰. (9)

Then, the temporal velocity and position for the particle are corrected by the obtained pressure gradient: 𝐮_𝑖^𝑘+1=𝐮^∗_𝑖+𝐮^′_𝑖𝐫_𝑖^𝑘+1=𝐫^∗_𝑖+𝐮^′_𝑖Δ𝑡, (10)(11) where 𝐮^′_𝑖=−Δ𝑡𝜌∇𝑝^𝑘+1_𝑖. (12)

The viscous term in (2) and the left-hand side of (9) are computed by the Laplacian operator (5). The right-hand side of (12) is calculated by the gradient operator (7).

3. Simulation Conditions

The flow through an impulse-type nano-hydraulic turbine developed by Ikeda et al. [6] is simulated. The power characteristic of the turbine was investigated by a laboratory experiment, and the flow field was made clear by the experimental visualization [6]. Figure 1 shows the rotor and blade specifications. The rotor diameter 𝐷 is 200 mm, the axis diameter is 15 mm, and the number of blades is 12. The blade has a circular-arc shape. The curvature radius is 23 mm, the thickness is 3 mm, and the chord length is 40 mm. To visualize the flow inside the rotor, the blades are sandwiched between two circular plates made of transparent vinyl chloride.

Figure 1 

Rotor and blade specifications.

In the experiment, the rotor is placed at the bottom of a waterfall as shown in Figure 2. The origin of coordinates is set at the initial point of the waterfall. The 𝑥-axis is horizontal, and the 𝑧-axis is vertical. The flow through the rotor is simulated by a two-dimensional MPS method. The water flow rate 𝑄 is 0.0035 m³/s, and the head of waterfall 𝐻_𝐹 is 570 mm. The nondimensional horizontal distance between the initial point of the waterfall and the colliding point of the water with the blade 𝐿_𝐹, 𝐿_𝐹/𝐷, is 1.33.

Figure 2 

Positional relation between waterfall and rotor.

The square region of 500 mm × 500 mm around the rotor axis is chosen as the computational domain as shown in Figure 3, where the coordinates of the axis are supposed to be (𝑥_𝑎,𝑧_𝑎). The simulation is performed at the conditions of 0.56≤𝜆≤0.7, where 𝜆(=𝑉_𝑡/𝑈_𝐹) is defined by the rotor tip speed 𝑉_𝑡 and the impact velocity of waterfall with the blade 𝑈_𝐹.

Figure 3 

Computational domain.

The waterfall flows into the computational domain from the location A on the upper boundary, as indicated in Figure 3. The particles, discretizing the falling water, are released from the location A into the domain. The releasing velocity and angle are determined from the condition of the free fall. It is required that the released particles are arranged at a uniform interval in the horizontal and vertical directions to ensure the computational accuracy. The particles are arranged at an interval 𝑙₀ in the horizontal direction and released into the domain at a constant time interval Δ𝑡_𝑝, so that the width of the waterfall 𝐵 at the location A corresponds to the measured value 8.4 mm. The Δ𝑡_𝑝 value is prescribed so that the vertical distance between the particles coincides with 𝑙₀. The distance 𝑙₀ corresponds to the grid width for grid-based simulations such as a finite difference method. Therefore, the simulated results are considered to depend on 𝑙₀. This study investigates the effect of 𝑙₀ on the flow and the turbine performance to search for the appropriate value of 𝑙₀.

The rotor axis and the blades, which are the solid walls, are discretized by the particles having the same angular velocity as the rotor. The distance between these particles is also set at 𝑙₀. The Neumann boundary condition for the pressure gradient is imposed on the particles contacting with the fluid.

The pressure is set at zero on the free surface. The position of the free surface is detected according to the value of the particle number density. When the particle number density ⟨𝑛^∗⟩_𝑖 obtained by the first-step calculation of each time step satisfies the following relation, the 𝑖th particle is decided to be on the free surface: 𝑛^∗_𝑖<𝛽𝑛⁰, (13) where 𝛽 is a parameter of 𝛽<1.

The time increment Δ𝑡 is determined from the maximum particle velocity at each computational time step. In this simulation, the initial value is set at 0.001 s, resulting in the maximum value of 0.001 s. The value of 𝑟_𝑒 in (3) is generally chosen at 2≤𝑟_𝑒/𝑙₀≤4 [7, 8]. The value is 2.1𝑙₀ for the particle number density and the gradient operator, while it is 4𝑙₀ for the Laplacian operator [8]. The parameter 𝛽 in (13) is set at 0.97. It is reported that the simulation of a fragmentation of fluid scarcely depends on the 𝛽 value in the case of 0.8≤𝛽≤0.99 [7].

The simulation conditions are listed in Table 1.

Table 1 

Simulation conditions.

4. Results and Discussion

4.1. Results at Tip Speed Ratio of 𝜆=0.56

Figure 4 shows the time variation of the particle distribution for the fully developed flow, where the tip speed ratio 𝜆 is 0.56. The distributions at five time points during a period for the interaction between the waterfall and a blade are presented. The distance between the particles discretizing waterfall 𝑙₀ is set at 𝑙₀/𝐵=0.084. When the time is as in Figure 4(a), the waterfall collides directly with the tip of the concave surface for the blade A. At the subsequent time points of Figures 4(b), 4(c) and 4(d), the collision point moves toward the center of the concave surface for the blade A as the rotor rotates. When the time is as in Figure 4(e), a part of the waterfall gets in contact with the tip of the subsequent blade B. But the collision between the waterfall and the concave surface for the blade A is still maintained. The collision point is closest to the rotor axis during the one period for the interaction between the waterfall and a blade. On the concave surface of the blade A, the number of particles or the mass of water increases with the passage of time during the one period. When the time is as in Figure 4(d), the water on the concave surface of the blade A flows toward the subsequent blade B, and it collides with the convex surface of the blade B at the time of Figure 4(e). It should be noted that the blade C corresponds to a blade, which has just finished colliding with the waterfall. The water on the concave surface flows toward the inside and outside of the rotor. The water flowing toward the inside collides with the concave surface of the subsequent blade A, and then it enters into the rotor. The water directing toward the outside disperses markedly in the radial direction as the rotor rotates. The flow rate toward the outside is much larger. The particles always exist on the convex surface of the blades, which are on the opposite side of the waterfall. This demonstrates the stagnation of water inside the rotor. The blades give the angular momentum to the stagnant water and flick the water away from the rotor, causing the deterioration of the rotor performance.

(a) 𝑡=𝑡1[𝑠]

(b) 𝑡=𝑡1+0.005

(c) 𝑡=𝑡1+0.01

(d) 𝑡=𝑡1+0.015

(e) 𝑡=𝑡1+0.02

Figure 4 

Time variation for particle distribution when 𝜆=0.56.

Figure 5 shows the superposition for the particle distributions at every time interval 3Δ𝑡 during 30Δ𝑡 time period, where 𝜆=0.56 and 𝑙₀/𝐵=0.084. One can grasp the water dispersion around the rotor and the water stagnation in the rotor. The water scatters mainly toward the lower left direction and the waterfall direction.

Figure 5 

Distribution of superimposed particles when 𝜆=0.56.

The time-averaged water velocity for 𝜆=0.56 is shown in Figure 6. The water dispersion toward the outside of the rotor and the flow into the rotor are reconfirmed.

Figure 6 

Time-averaged velocity distribution when 𝜆=0.56.

The flow pattern inside and around the rotor, experimentally visualized by Ikeda et al. [6] at 𝜆=0.56, is presented in Figure 7. It was acquired by using a CCD camera and a strobe light sheet shaped through a 2 mm wide slit. The image visualizes vividly the water flow along the concave surface of the blades, the water dispersion and flick toward the outside of the rotor, and the water stagnation inside the rotor. The simulated flow, shown in Figure 4, agrees well with the experimental visualization, demonstrating the validity of the present simulation.

Figure 7 

Experimentally visualized flow pattern when 𝜆=0.56.

The distance between the particles discretizing the waterfall 𝑙₀ or the initial distance between the particles corresponds to the grid width in grid-based simulation methods, such as a finite difference method. The simulation with smaller 𝑙₀ has superior space resolution. The time-averaged number of particles, 𝑁, for the fully developed flow is plotted against 𝑙₀ in Figure 8. For the abovementioned simulation of 𝑙₀/𝐵=0.084, 𝑁 is 14437. The 𝑁 value increases greatly with the decrement of 𝑙₀. The improvement of the spatial resolution increases the number of particles, and therefore it causes the increment of the computational time. In the simulation of 𝑙₀/𝐵=0.084, 348 and 299 particles are used to discretize the axis and the single blade, respectively. These particles are included in Figure 8. For the flow simulation during the one revolution of the rotor, about 6.4 hours are required on a workstation (Processor: Intel Xeon X5660 2.8 GHz × 6, Memory: 12 GB).

Figure 8 

Change in number of particles 𝑁 due to distance between particles 𝑙0 when 𝜆=0.56.

The superimposed particles at 𝑙₀/𝐵=0.069 and 0.12 distribute as shown in Figure 9, where 𝜆=0.56. In the case of 𝑙₀/𝐵=0.12, the number of particles is low, and the water dispersion outside the rotor and the stagnation inside the rotor are not fully resolved. The simulation of 𝑙₀/𝐵=0.069, composed of more particles, yields almost the same result at 𝑙₀/𝐵=0.084 shown in Figure 5. It is discovered that the simulation of 𝑙₀/𝐵≤0.084 can ensure a sufficiently high spatial resolution.

(a)  In case of 𝑙0/𝐵=0.12

(b) In case of 𝑙0/𝐵=0.069

Figure 9 

Effect of distance between particles 𝑙0 on particle distribution when 𝜆=0.56.

Calculating the rotor power 𝑃 by estimating the water kinematic energy at the rotor inlet and outlet, the nondimensional value 𝐶_𝑝 changes as the function of 𝑙₀/𝐵 as plotted in Figure 10. When the spatial resolution is high enough (𝑙₀/𝐵≤0.084), the 𝐶_𝑝 value remains almost unaltered. It is slightly larger than the experimental result of 𝐶_𝑝=0.66. But it is considered to be valid in due consideration of the two-dimensional simulation. Ikeda et al. [6] calculated the 𝐶_𝑝 value at 𝜆=0.56 by estimating the force on a blade from the visualized flow pattern. It is also plotted in Figure 10, being in good agreement with the present simulation. Consequently, it is found that the simulation of 𝑙₀/𝐵≤0.084 can accurately predict the 𝐶_𝑝 value.

Figure 10 

Change in power coefficient 𝐶𝑝 due to distance between particles 𝑙0 when 𝜆=0.56.

Ikeda et al. [6] made it clear by their experiment that the horizontal distance between the initial point of the waterfall and the blade, 𝐿_𝐹 affects the rotor performance 𝐶_𝑝. The distance 𝐿_𝐹 varies with the water flow rate 𝑄. For the practical use of the hydraulic turbine, it is desirable that 𝐿_𝐹 is always set at the optimal value irrespective of 𝑄. Ikeda et al. [6] proposed a method to control the 𝐿_𝐹 value by installing a flat plate along the waterfall. The current simulation is effectively employed to search for the applicability of the method.

4.2. Results at Tip Speed Ratios of 𝜆=0.6 and 0.7

The flows at 𝜆=0.6 and 0.7 are simulated for the condition of 𝑙₀/𝐵=0.084. Figure 11 depicts the relation between 𝜆 and 𝐶_𝑝. The 𝐶_𝑝 value decreases with the increment of 𝜆. This change agrees with the experimental result of Ikeda et al. [6], indicating that the effect of 𝜆 on 𝐶_𝑝 is successfully analyzed by the simulation. The simulated 𝐶_𝑝 value is slightly larger. This is because the current simulation does not sufficiently resolve the turbulent flow, and therefore it ignores the losses caused by the turbulence on the blade surface as well as inside the rotor. It may be also because the current simulation employs the two-dimensional MPS method. Using the three-dimensional MPS method, the flow and the 𝐶_𝑝 value would be simulated more accurately. But the number of particles increases, and a longer computational time is required.

Figure 11 

Relation between power coefficient 𝐶𝑝 and tip speed ratio 𝜆.

The flow fields at 𝜆=0.7 are shown in Figures 12 and 13. The waterfall disperses due to the collision with the rotor. When compared with the results at 𝜆=0.56 (Figures 5 and 6), the divergence angle of the dispersed water lessens. One can grasp the decrement of the rotor angular momentum obtained from the waterfall. The amount of water inside the rotor is less than that at 𝜆=0.56.

Figure 12 

Distribution of superimposed particles when 𝜆=0.7.

Figure 13 

Time-averaged velocity distribution when 𝜆=0.7.

5. Conclusions

The flow through an impulse-type small-scale hydraulic turbine utilizing a waterfall of extra-low head is simulated by a two-dimensional MPS method. The rotor performance is also analyzed by using the simulated flow field. The results are summarized as follows.(1)When the distance between the particles discretizing the waterfall of a width 𝐵, 𝑙₀, is set at 𝑙₀/𝐵≤0.084, the flow simulated at the tip speed ratio 𝜆=0.56 is confirmed to agree well with the experimentally visualized one. Thus, the simulation of 𝑙₀/𝐵≤0.084 has a sufficiently high spatial resolution.(2)For the simulation of 𝑙₀/𝐵≤0.084 at 𝜆=0.56, it is also confirmed that the simulated power coefficient 𝐶_𝑝 agrees nearly with the experiment. Thus, the simulation of 𝑙₀/𝐵≤0.084 can also favorably predict the rotor performance.(3)The 𝐶_𝑝 values simulated at 𝜆=0.6 and 0.7 agree almost with the experimental results. Therefore, the present simulation can successfully analyze the effect of 𝜆 on the flow and the rotor performance.

Nomenclatures

𝐵:	Thickness of waterfall
𝐶𝑝:	Power coefficient =𝑃/𝜌𝑔𝑄𝐻𝐹
𝑑:	Number of space dimensions
𝐷:	Diameter of rotor
𝐅:	External force
𝐻𝐹:	Head of waterfall
𝑙0:	Distance between particles discretizing waterfall
𝐿𝐹:	Horizontal distance between blade and initial point of waterfall
𝑛:	Particle number density
𝑁:	Time-averaged number of particles
𝑝:	Pressure
𝑃:	Power output from rotor
𝑄:	Water flow rate
𝐫:	Position vector of particle
𝑡:	Time
𝐮:	Velocity
𝑈𝐹:	Impactvelocity of waterfall with blade =(2𝑔𝐻𝐹)1/2
𝑉𝑡:	Rotor tip speed =𝜔𝐷/2
𝑤:	Weight function
𝑥,𝑦,𝑧:	Spatial coordinates
Δ𝑡:	Time increment
𝜆:	Tip speed ratio =𝑉𝑡/𝑈𝐹
𝜈:	Kinematic viscosity
𝜌:	Density
𝜔:	Angular velocity of rotor.