Pipeline systems in many industries and power plants are usually characterized by its layout complexity, which consists of many main long straight pipes and secondary pipes connected by sharp bends. In nuclear power plants, the flow downstream of a 90˚ bend is essential for the primary and secondary cooling systems, where many sharp bends are used to interconnect the components. However, fluid flows through a 90˚ sharp bends are a very complex phenomenon. As the fluid flows through a bend, the centrifugal force acting on the fluid develops a radial pressure gradient. Because of the pressure gradient in the fluid, the secondary flow is generated downstream of the bend. In addition, a significant pressure gradient produces swirling flow downstream the bend pipe  . Under the swirling flow condition, the significant pressure fluctuation and the high-velocity fluctuation occur in the bends. These are sources of flow-induced vibration (FIV). Consequently, FIV can cause the fatigue fracture in a piping system of the industries and power plants  . Therefore, the investigation of fluid flow field and velocity fluctuation is essential issues for the pipeline safety.
The structure of the secondary flow in the bends is dependent on the bend curvature radius (Rc) and Reynolds number (Re). Some researchers have investigated the flow characteristics in the 90˚ single bend pipe flow with curvature ratio of Rc/D > 1.3 (Rc: curvature radius, D: pipe diameter) both experimentally and numerically      . When the non-dimensional curvature ratio of the bend is greater than 1.5, the secondary flow which consisting of a pair of counter-rotating vortices (Dean vortices) are generated. Meanwhile, the velocity profile of the primary stream wise flow is distorted and shifted away from the center of the curvature of the elbow  . For the small curvature ratio, Kawamura et al.  studied the flow separation using the bends with curvature ratio of Rc/D = 1.1, 2, 4. The experiments were done by using Laser Doppler Velocimetry (LDV) at Reynolds numbers of 5 × 104 < Re < 1 × 105. They revealed that the flow separation occurred at Rc/D = 1 and did not occur at Rc/D = 2 and Rc/D = 4. Also, they found that the power spectrum of the turbulence intensity downstream near the elbow with any curvature ratio and Reynolds number had a distinct peak at the reduced frequency of about 0.5. Later on, Ono et al.  investigated the water experiments in two bends with different curvature ratio (Rc/D = 1 and Rc/D = 1.5) in order to investigate the interaction between flow separation and the secondary flow due to the elbow curvature. Particle Image Velocimetry (PIV) was used in their experiment, and they confirmed that the flow separation always occurred in the short-elbow (Rc/D = 1) while the flow separation occurred intermittently in the long-elbow case (Rc/D = 1.5).
The experiments on the above previous studies mainly were conducted under the inlet condition of fully developed flow and flat velocity profile. However, non-uniform velocity profile on the inlet condition might be appeared in a specific condition. Kubo et al.  confirmed that a swirling flow occurred in the dual elbow. The experiments were done by two-dimensional Particle Image Velocimetry (2D-PIV) at Reynolds number Re = 3 × 105 with three types of curvature ratio Rc/D = 1, 1.5, 2. They found that the swirl intensity of the swirling flow, which was generated in the dual elbow, became high and fluctuated largely as the curvature ratio was small. The influence of inlet swirling flow on 90˚ bend pipe was studied in order to understand the flow structure downstream the bend pipe flow  . Sudo et al.  measured the pressure and velocity distribution at Re = 6 × 104 with a curvature ratio of Rc/D = 2 under the weak inlet swirling flow condition. They found that inlet swirl flow affected the behavior of secondary flow generating at the downstream of the bend pipe flow. Later, Kalpakli and Orlu  experimentally studied the influence of swirling intensity at the inlet on the 90˚ bend pipe flow. The investigation was done using stereoscopic particle image velocimetry. When the swirl number was increased to S = 1.2, it was found that the Dean vortices at the downstream of bend pipe became a single large-scale motion located almost at the center of the pipe and flow became a dominated swirl flow. Swirling inlet flow condition is also encountered in the cold-leg piping system in a nuclear reactor. Therefore, Yamano et al.   investigated the effect of swirl flow at the inlet on the 90˚ double bend pipe both experimentally and numerically, to evaluate the flow-induced vibration for primary cooling pipes in the Japan sodium-cooled fast reactor (JSFR). The experiment was done by PIV method with curvature ratio Rc/D = 1. They found that the flow separation region was deflected at the downstream from the bend when the inlet condition on the first bend was swirling flow. Later, Mizutani et al.  investigated the influence of inlet condition upstream the triple elbow. Particle Image Velocimetry (PIV) was used for the velocity flow mapping, and the tangential injection method was used to generate swirling flow on the inlet. The triple elbow with curvature ratio Rc/D = 1 was used to be close to an actual condition and to accumulate knowledge towards optimization of a prospective piping layout in the conceptual design of Japan Sodium Fast Reactor (JSFR)  .
In the previous studies, the researchers mostly used the optical system to measure the velocity field and velocity fluctuation. However, the optical system has some challenges to apply in non-transparent wall channel or pipe. Also, it is difficult for the applications in actual plant process. Thus, another measurement technique should be developed to evaluate the velocity field and velocity fluctuation. In addition, the influence of the strong swirling inlet on the double bend with curvature ratio Rc/D < 1 has not been done yet. So in this study, the influence of the strong swirling inlet condition on the velocity fluctuation at the 90˚ double bend pipe with curvature ratio Rc/D = 0.5 is investigated with ultrasound technique. Ultrasound technique measurement is used because its advantages which can be applied in the non-transparent pipe and opaque liquid flow  .
Initially, conventional Ultrasonic Velocity Profiler (UVP) method only measures one-dimensional velocity profile in the measurement line. In the case of two-dimensional velocity vector measurement, Takeda and Kikura  investigated velocity field of the mercury flow using UVP system with multiple transducers. Nevertheless, the measurement system using multiple transducers is quite large as the number of transducers is increased. To minimize this problem, the authors had developed a phased array UVP system. A phased array sensor has multiple ultrasonic piezoelectric elements; an ultrasound beam can be steered to a specific angle by controlling time delay of ultrasound transmission from each piezoelectric element. Therefore, velocity profiles can be measured in multiple measurement lines. Also, Phased Array UVP can measure instantaneous two-dimensional velocity profile. The performance of Phased Array UVP for velocity flow mapping was confirmed by Fukumoto et al.  for a detecting a water leakage in the tank.
The main objectives of this study are to clarify the reattachment point of secondary swirling flow that occurs just downstream of the double bend pipe and to observe the influence of strong swirling inlet flow on the velocity fluctuation of the 90˚ double bend pipe flow using ultrasound technique. The bend curvature caused the secondary swirling flow, and the high-velocity fluctuation occurs near the reattachment point region. To achieve this purpose, Phased Array UVP system is utilized to measure the two-dimensional velocity and the axial and tangential velocity fluctuation at the secondary flow region.
2. Measurement Principle and Experimental Apparatus
2.1. Two-Dimensional Velocity by Phased Array UVP
The working principle of Phased Array UVP system based on Doppler shift frequency detection along ultrasound beam lines. Phased array sensor emits an ultrasonic pulse, and each piezoelectric element of sensor receives the echo reflected from the surface of a particle. The exciting element emits a spherical ultrasonic wave. When adjacent elements emit within a close second, interference of wave fronts occurs as shown in Figure 1(a). Figure 1(b) shows a schematic diagram of phased array sensor which a is element size, d is inter-element spacing, and c is element width. The pattern of the interference depends on the time delay, the steering angle of ultrasonic beam can be changed at the time delay Δt. The steering angle θs and the time delay Δt is related with the speed of sound in a medium c and inter-element spacing d as shown in Equation (1).
Figure 1. Beam steering principle of Phased Array sensor and its schematic diagram. (a) Beam steering; (b) The schematic diagram of phased array sensor.
Basic equation of Doppler shift is derived from Doppler equation as shown below:
where fd is the Doppler frequency, s is the speed at which object is approaching the transducer, f0 is the basic frequency of the transducer, and c is the speed of sound in the medium, i.e., water (c = 1480 m/s at 20˚C).
If the object is moving at an angle θ to the transducer, then s = V *cosθ. By substitution, we get the Doppler shift equation for a single transducer:
Equation (3) can be rewritten as:
where 2 *cosθ applies to a roundtrip Doppler shift using a single transducer.
If two transducers are used, one receiver and one transmitter, as shown in Figure 2. Then, 2 *cosθ becomes cosθ + cosγ, where θ is the angle between the transmitter and the vector V, and γ is the angle between the receiver and the vector V. Moreover, α = θ − γ, where α is the angle between the transmitter and receiver, then cosγ = cos(θ − α).
The Doppler equation becomes:
Using the trigonometric identity:
The Doppler Equation (5) becomes:
Figure 2. A two-transducer ultrasound system produces a Doppler signal from velocity vector V.
In this paper, phased array sensor is used as transmitter and receiver (transceiver). As shown in Figure 3, the number of piezoelectric elements is eight elements. These eight elements transceiver transmit ultrasound beam at the same time or delayed for certain time to get steering angle, then the beams from these transceivers will interference pattern forming (beamforming) into one beam.
For measuring two-dimensional velocity vector, the development system uses two piezoelectric elements as transceivers to calculate the actual velocity magnitude and angle from the returned signal at specific measuring volume, i.e., Channel 4 (Ch 4). The Doppler shift equation for the piezoelectric element number 8 is identical to Equation (7):
where fd8 is the Doppler frequency received by the piezoelectric element number 8 at angle α (measured clockwise from the axis of the transmitting beam), f0 is basic frequency, V is the magnitude of the velocity of the particle travelling at angle θ (measured clockwise from the axis of the transmitting beam), and c is the speed of sound in fluid i.e. water.
The Doppler shift equation for the piezoelectric element number 1 is:
where fd1 is the Doppler frequency received by the piezoelectric element number 1 at angle β (measured clockwise from the axis of the transmitting beam), and all other terms are identical to those of the piezoelectric element number 8.
The signals from both of these transceivers are demodulated with the transmitted frequency f0 to produce four quadrature signals. It should be noted
Figure 3. Velocity vector reconstruction.
that at a given depth (channel) from the surface of both piezoelectric elements f0, c, α, and β are all-constant due to the fixed geometry of the piezoelectric elements. The real velocity vector V is calculated by multiplying the four-quadrature signals from the two transceivers elements number 1 and 8. The multiplying produces two subcomponents, which are the sum of the frequencies and the difference between the frequencies. Using trigonometry:
From Equation (8) and Equation (9), we will get the equations for the sum and difference as follow:
If sin β = −sinα and cosβ = cosα, then Equation (11) becomes
If sinβ = −sinα and cosβ = cosα, then Equation (13) becomes
If Vx = V * sinθ and Vy = V * cosθ, then from Equations (12) and (14):
Since these Vx and Vy are orthogonal, the real magnitude can be determined by vector addition, and simple trigonometry can determine the angle:
The spatial resolution or channel distance is defined as:
where Δy is channel distance, Ncycle is a number of cycles per pulse, c is the speed of sound, and f0 is the basic frequency of the transducer. If Ncyle = 2, c = 1480 m/s, and f0 = 2 MHz, then channel distance (Δy) is 0.74 mm.
Phased Array UVP system is showed in Figure 4. National Instrument LabVIEW program is used to control Phased Array UVP system and reconstruct two-dimensional velocity vector. The measurement system consists of a 2 MHz phased array sensor with eight piezoelectric elements, 8-channel pulse receiver, analog to digital converter and personal computer to control the pulse receiver and analyze echo signal from the digitizer.
2.2. Near-Field Effect of Phased Array Sensor
In Phased Array UVP system, we have to consider the effect of the near field oscillation. The high oscillation burst occurs near the active elements surface. It influences the accuracy of the measurement close to the sensor (near-field region). The near-field boundary of phased array sensor has been numerically investigated as shown in Figure 5. The numerical calculation is based on the work by Ocheltree and Frizzel  and performed using MATLAB®. Near-field oscillation boundary depends on the width of piezoelectric element b and the wavelength of the ultrasound pulse λ. From Equation (20), we can estimate near-field oscillation length Noscillation of phased array sensor.
where, ka is correction factor, and b is piezoelectric total element width. The correction factor can be calculated from the rectangular element ratio (b/c) = 1
Figure 4. Phased array UVP measurement system.
as shown in Figure 6. From the calculation, the near-field oscillation boundary is 5.78 mm from the element surface. Hamdani et al.  confirmed the near- field effect on the velocity measurement. Velocity data on the near-field region showed a bad accuracy. Nevertheless, the velocity data beyond the near field showed a good agreement with conventional UVP, and they confirmed that the angle uncertainty for two-dimensional velocity was ±1˚.
2.3. Experimental Apparatus
The experiment was conducted in a horizontal water circulation system at atmospheric pressure, which consisting of the cooling system, electromagnetic flow meter, pump, ball valve, a bypass pipe, and flow conditioner as shown in Figure 7. Water is supplied to the test channel from the water storage tank by a 3 phase induction motor pump (65 × 50 FSED EBARA CORPORATION). The flow rate of the liquid phase is controlled using the control ball valve. The flow rate is measured by magnetic flow meter (MAGMAX MGC010 Tokyo Keiso Co., Ltd.) with accuracy ±0.5%. Fluid enters the straight pipe test section after an elbow. Therefore a flow conditioner is installed to eliminate the flow distortion
Figure 5. Sound pressure simulation of 2 MHz phased array sensor.
Figure 6. Phased array UVP measurement system. (a) Rectangular active element; (b) Correction factor graph (ka).
effects of the elbow. The length of straight pipe before entering the first bend is 42D (D = 50 mm) in order to achieve fully developed turbulent pipe flow condition at the inlet of the first bend. In this experiment, the double bend pipe is utilized to investigate the secondary swirling flow. The double bent pipe has a bent angle of 90 degrees and curvature ratio Rc/D = 0.5. The Reynolds number is Re = 1 × 104, based on the bulk velocity and the pipe diameter. The cooling system is used to control and maintain a constant water temperature. The temperature is recorded using thermocouple during the measurements, and it is confirmed that the water temperature fluctuation is ±1˚C.
2.4. Swirling Generator
For the generation of swirling flow, some different methods exist (e.g., pipe rotation, tangential injection, guide-vane, twisted tape, helical turbulators, and propeller-type), which have a different effect on the main flow  . In this present study, we used the rotary swirler   . The rotary swirler is used as swirl generator because it is easy to control the swirl intensity. The enlarged top view of the test section with therotary swirler is shown in Figure 8.
The swirl generator is installed at 12D upstream of the double bent pipe. The rotary swirler consists of a 150 mm long aluminum pipe with inner diameter 50
Figure 7. Experimental apparatus of water circulation system.
Figure 8. Enlarged top view of rotary swirler and test section.
mm. Small tubes of diameter 3 mm and 50 mm long are inserted into the aluminum pipe. The small tubes are packed as tightly as possible, and their number is approximately 95. The pipe can be rotated about its axis at speed varying from 15 to 1,100 r.p.m. by induction motor (5IK90SW-5 Oriental Motor Co., Ltd.) and a timing belt (K40L50BF) connects it. With the pipe rotating, the growth of the boundary layer on the pipe walls establishes an azimuthal velocity distribution corresponding to the solid-body rotation in the core. Whereas the rotary swirler allows the solid rotation of the small tubes in peripheral direction, while the axial velocity distribution is made uniform by the function of the small tubes structure  . Therefore, the swirl intensity can be defined as the ratio of the circumferential momentum to the axial momentum.
For axisymmetric flow, swirl intensity is usually defined by its swirl number (S). A parameter S that is used by several researchers   is defined as follow:
where r is the radial distance from a pipe axis, vx is the streamwise mean velocity and vy is the circumferential mean velocity.
An alternative to the parameter S on the rotating swirleris defined as follow  :
Equation (22) indicates that the swirl intensity can be evaluated directly from the angular velocity of the rotary pipe, the diameter D of the pipe and the bulk velocity Um of the flow through the pipe.
2.5. Experimental Procedures
A phased array sensor, which has basic frequency 2 MHz, is installed through the pipe wall. Thus, there is a direct contact between sensor and fluid to overcome the refraction in the pipe wall. The cross-sectional plane measurement is performed at 7D downstream of the swirling generator to observe developed swirling flow. The angular velocity of rotating pipe is measured by a digital optical tachometer (AD-5172 A & D Company) with accuracy ± 0.01% ± 1 digit (10 - 6000 r.p.m.). According to the literature  , S = 1 corresponds to the downstream development of initially strong swirl in pipe flow. Therefore, in this experiment, the swirl number is set to S = 1 in order to know the influence of strong swirl intensity at the inlet. For the cross-sectional plane measurement, to include the whole 360 degrees of the pipe, one measurement position will be taken 20 degrees, which adds up to 18 measurement lines.
Table 1 shows the experimental conditions of the measurements. In case of axial plane measurements, there are four measurement positions at x/D = 0.2,
Table 1. Experimental conditions and phased array UVP parameters.
x/D = 0.6, x/D = 1 and x/D = 1.5 downstream of the double bend. After the measurement is done at the first position, the sensor is moved to another position. The maximum steering angle of the phased array sensor is −10 degree to 0 degree. In each measurement position, we measured 3 measurement lines, and the interval between each measurement lines is 5 degree as shown in Figure 9.
3. Results and Discussion
The experiments were done in the case of with and without inlet swirling flow condition. In the case of without inlet swirling flow condition, Figure 10 and Figure 11 shows the experiment results of one-dimensional normalized axial and tangential velocity with their standard deviations in four positions downstream of the double bend pipe.
Ten thousand instantaneous velocities are averaged in each velocity profile. The vertical axis indicates the dimensionless distance of measurement line through the pipe. The horizontal axis is the dimensionless axial and tangential velocity normalized by the average axial velocity and the average tangential velocity respectively. According to Figure 10, the low velocity magnitude was observed at x/D = 0.2 to x/D = 1 in the extrados pipe wall region (y/D = 0.7 to y/D = 1). This region (y/D = 0.7 to y/D = 1) is the secondary flow region. As we can see the velocity profiles in secondary flow region are slightly different in every axial position. These differences are related to the flow structure of the fluid, which flow direction is reverse to the main flow. Nevertheless, the velocity magnitude starts to increase at x/D = 1. Therefore, we can assume there is a no reverse flow beyond x/D = 1. Figure 11 shows one-dimensional time-average tangential velocity at each measurement positions. The tangential velocity magnitude is relatively high on the secondary flow region, start from x/D = 0.6 to x/D = 1.5.
Figure 9. Schematic view of the measurement positions.
Figure 10. Time-average axial velocity without inlet swirling flow (S = 0).
Figure 11. Time-average tangential velocity without inlet swirling flow (S = 0).
At the center of the pipe (y/D = 0.5), the velocity magnitude is relatively lower compared to the end wall regions. The reason is that the bend pipe curvature generate the centrifugal force and it enhances tangential velocity magnitude near the pipe walls.
Figure 12 shows the one-dimensional time-average axial velocity with their standard deviations downstream of the double bend pipe with inlet swirling flow condition (S = 1). In the case without and with inlet swirling flow condition, at each axial measurement positions, there is no significant difference on the velocity profiles in the main flow region (0 < y/D < 0.7). However, in case of with swirling flow, the velocity magnitude in the secondary flow region (0.7 < y/D < 1) is stronger than without inlet swirling flow condition. These differences are due to the inlet swirling flow enhance the centrifugal force on the curvature bend. Figure 13 shows the one-dimensional time-average tangential velocity profile of inlet swirling flow with their standard deviations. The strong tangential velocity occurs at x/D = 0.6 and x/D = 1 because the swirling enhance the tangential velocity in the main flow region (0.1 < y/D < 0.4) and the secondary flow region (0.7 < y/D < 1). These results indicate that tangential velocity in swirling flow is more dominant on main flow regime than without swirling flow condition. According to the standard deviation, the tangential velocity fluctuation in the inlet swirling condition is higher than without inlet swirling flow. It means that the inlet swirling flow also increases the tangential velocity fluctuation.
Figure 12. Time-average axial velocity with inlet swirling flow (S = 1).
Figure 13. Time-average tangential velocity with inlet swirling flow (S = 1).
3.1. Two-Dimensional Velocity with and without Inlet Swirling Flow
Two-dimensional velocity measurement results are plotted to know the flow structure downstream of 90˚ double bend pipe in the condition of with and without inlet swirling flow. Figure 14 shows two-dimensional time average axial velocity without inlet swirling flow at four downstream locations of the second bend pipe. The velocity vector cannot be measured enough near the oscillation boundary of the phased array sensor from 0 < y/D < 0.08 for all positions due to the ultrasound field oscillation effects. From the measurement result at x/D = 0.2, the direction of velocity vectors separate from the edge of the bend. It indicates the flow separation starts from the edge of the bend at x/D = 0.2 (i.e., separation point). Then, the separated fluid attaches near the first line measurement of x/D = 1.5 (i.e., reattachment point). The secondary flow region is located between the separation point and the reattachment point. It means that there is a no reverse flow downstream the reattachment point and the flow acceleration occurs beyond the reattachment point.
Firstly, to confirm the inlet swirling condition, we visualized two-dimensional radial velocity as shown in Figure 15. The measurement is taken at 7D
Figure 14. Two-dimensional time average axial velocity without inlet swirling flow (S = 0).
Figure 15. Two-dimensional time average radial velocity with the inlet swirling flow (S = 1).
downstream from the swirling generator. It is confirmed that the generated swirling flow is symmetric, and the velocity distribution is homogenous. In the core region of the pipe, the velocity magnitude is lower than near wall region. The highest velocity magnitude is between the core and near wall region.
After confirming inlet swirling flow, we investigated the influence of swirling flow on the flow structure downstream of the double bend. Figure 16 shows the measurement results in case of swirling flow condition. The result shows the flow separation phenomenon is same as the condition of without inlet swirling flow. However the reattachment point shift to upstream around 25 mm when swirling flow is introduced as inlet condition. The reattachment point occurs near the first measurement line of x/D = 1. Therefore, the reverse flow region is narrow. At x/D = 1.5, the fluid becomes the accelerated swirling flow.
3.2. Comparison of Velocity Fluctuation between with and without Inlet Swirling Flow
For the investigation of velocity fluctuation, we calculated turbulent intensity for axial velocity and tangential velocity in the condition of without inlet swirling flow and with inlet swirling flow respectively. The x-axis represents the normalized distance, and the y-axis represents turbulent intensity. The turbulent intensity of axial velocity and tangential velocity fluctuation is very effective to analyze the velocity fluctuation. Figure 17 and Figure 18 show the comparison of turbulent intensity in axial and tangential velocity at the position of x/D = 1 downstream of the double bend pipe. According to axial turbulent intensity graph, the axial velocity fluctuation is relatively high between y/D = 0.6 and y/D = 0.8 (see Figure 17). It means that the high-velocity fluctuation occurs at the boundary of the main flow and secondary flow. Moreover, the influence of reattachment point exists in this region. From the comparison, the axial velocity fluctuation of without inlet swirling is higher than inlet swirling. The reason is that the inlet swirling flow suppress axial velocity fluctuation.
Figure 18 shows tangential turbulent intensity analysis. This graph indicates
Figure 16. Two-dimensional time average axial velocity with the inlet swirling flow (S = 1).
Figure 17. Comparison of the axial velocity fluctuation between with and without inlet swirling flow.
Figure 18. Comparison of the tangential velocity fluctuation between with and without inlet swirling flow.
the tangential velocity fluctuation is high near the extrados sidewall and at the center of the pipe. When the fluid touches the extrados sidewall, the flow direction changes. Some fluids flow as the main flow and some fluids change as the secondary flow. Therefore in the main flow region (near the extrados) and at the core of the secondary swirling flow (at the pipe center) occurs the high tangential velocity fluctuation. Here, the maximum turbulent intensity is around 40% with inlet swirling flow, and the maximum turbulent intensity is around 25% without swirling flow. Therefore, the velocity fluctuation of inlet swirling flow is stronger than without inlet swirling flow. It seems that swirling flow enhances tangential velocity fluctuation.
Phased Array UVP system was applied for two-dimensional velocity measurements in the condition of without inlet swirling flow and with inlet swirling flow on the double bend pipe flow. According to one-dimensional velocity profiles, the velocity magnitude of inlet swirling flow is much stronger than without inlet swirling flow in the secondary flow region. These differences are due to the influence of the inlet swirling flow. In two-dimensional velocity of without inlet swirling flow, the flow separation occurs around x/D = 0.1, and the reattachment point is located at x/D = 1.5. In case of inlet swirling flow, the flow separation phenomenon is same as the condition of without inlet swirling flow, but the reattachment point is located at x/D = 1. Therefore, the reverse flow region is narrow. At x/D = 1.5, the fluid becomes the accelerated swirling flow. In the condition of without inlet swirling flow, the axial velocity fluctuation is higher than with inlet swirling flow, but tangential velocity fluctuation is lower than inlet swirling flow.
We are grateful to Professor Nobuyuki Fujisawa in Niigata University, Japan, for his advice on the design of swirling generator.