Numerical Simulation and PIV Experimental Investigation on Underwater Autorotating Rotor

Abstract

In this work, the flow field of an autorotating rotor in a water tunnel with various pitches and shaft backward angles was investigated via particle image velocimetry (PIV). The experiments were carried out on a free-rotating two-bladed single rotor. Computational Fluid Dynamics (CFD) based on moving overset grids were developed to study the hydrodynamic characteristics of an underwater autorotating rotor. The simulation results are in good agreement with the test results. The thrust and thrust coefficient of the underwater autorotating rotor were calculated by CFD simulation under different situations. The research demonstrates that rotational speed and thrust have a significant positive correlation with water velocity, pitch, and shaft back angle. In particular, the thrust coefficient scarcely varies with the shaft backward angle. An underwater autorotation rotor with a thin airfoil, negative torque, and a suitable number of blades can increase the thrust and thrust coefficient. The investigation is of significance in enriching the autorotation theory of rotors and helping to develop underwater autorotating rotors.

1. Introduction

Underwater vehicles are playing an important role in marine resource exploration and marine environmental investigation [1,2]. At present, the majority of underwater vehicles require power to provide propulsion for navigation [3]. It is important to develop an underwater vehicle that achieves long-time unpowered navigation with thrust provided by an autorotating rotor.

Autogiro is a mature aircraft. It has wide prospective applications in the future [4,5]. The autogiro was developed by Juan de la Cierva [6,7] in the 1920s, who introduced the mature age of rotorcrafts. Leishman [8,9] introduced the development history of rotorcraft in detail. Most of the research on autorotation was conducted in the early 20th century. Glauert [10] proposed a theory of autorotorcraft based on simplified assumptions, which established the basic theory for investigating the aerodynamic characteristics of autorotorcraft. Wheatley [11,12,13,14,15,16,17] carried out a lot of research on autorotorcraft. A scaled-down model of the VPM M16 rotorcraft was tested in a wind tunnel at Glasgow University, UK. Houston [18,19,20,21,22,23] investigated the advantages and disadvantages of this rotorcraft design scheme, the system identification technique, the maneuvering stability characteristics, the flight dynamics characteristics, and the flight quality. Bagiev [24,25,26] validated the mathematical simulation model of the rotorcraft in terms of leveling and maneuvering response by comparing numerical calculations with experimental results.

However, most of the research is mainly focused on autogiro in air, and less research has been done on autorotating rotors. With the development of technology, the performances of autorotation rotors can be better demonstrated. The autorotation state of the rotor can also be applied underwater, making so-called autonomous underwater vehicles [27,28,29], which is another meaning for this water tunnel experiment.

Particle image velocimetry (PIV) is a very mature technique developed in the 1970s, specializing in investigating flow features quantitatively [30,31]. For the helicopter rotor field, Yu et al. investigated the induced velocity and flow feature of the coaxial rotor in hovering and forward flight conditions in a water tunnel via PIV [32,33,34]. They compared the flow feature of coaxial rotors with single rotors in detail. With three-component PIV development, the flow features of full-equipped helicopters can be comprehensively investigated with multiple planes of observation windows [35,36,37].

Computational Fluid Dynamics (CFD) is a very mature technique developed in the 1960s, specializing in the flow field. Its advantages are low cost, short cycle time, and high accuracy. Many advanced CFD techniques have been used in the numerical simulation of helicopter rotors [38,39,40]. In particular, the overset grid technique [41,42,43,44,45,46,47,48] has been widely used in helicopter rotor simulation, and the method has been employed for various operating conditions such as hovering, forward flight, vertical flight, and near-ground flight of helicopters. PIV tests and CFD calculations for underwater rotorcraft are scarce [49]. Most of the studies have been conducted on rotorcraft flying in the air. Kim et al. [50] used computational fluid dynamics (CFD) for the numerical simulation of a six-degree-of-freedom (6DOF) free-running maneuver of an underwater vehicle. Dantas et al. [51] analyzed the hydrodynamic efforts generated on an AUV considering the combined effects of the control surface deflection and the angle of attack using CFD software based on the Reynolds-averaged Navier–Stokes formulations. Kilavuz et al. [52] conducted a numerical and experimental study of the flow characteristics of an unmanned underwater vehicle (UUV) with the commonly used Myring profile using computational fluid dynamics and particle image velocimetry techniques. Gao et al. [53] estimated the drag on an autonomous underwater vehicle (AUV) and optimized the shape using CFD. Hong et al. [54] conducted a numerical study of the hydrodynamic performance of a portable autonomous underwater vehicle (AUV). Current research on underwater vehicles is mainly focused on powered UAVs. There is a shortage of research on UAVs that rely on rotor autorotation to provide thrust and achieve unpowered navigation. In the future, unpowered underwater autorotates can be used for undersea mapping and military reconnaissance. Consequently, there is an enormous requirement to conduct experimental research and simulation analysis on underwater autorotating rotors.

The purpose of the paper is to investigate the hydrodynamic characteristics of the underwater autorotating rotor. Both an experimental and a numerical study have been performed. The rotational speed and flow field of the autorotating rotor was investigated in a water tunnel via the particle image velocimetry (PIV) technique. A three-dimensional overset grid numerical simulation was conducted on an underwater autorotating rotor. The test results were compared with the simulation data to validate the CFD simulations. The effect of water velocity, shaft backward angle, pitch, twist angle, airfoil, and blade number on the thrust and thrust coefficient was presented in detail.

2. Experiment Setup

2.1. Water Tunnel and Rotor System

The experiments were carried out in the low-speed water tunnel of the Key Laboratory of the Ministry of Education at Beihang University. As shown in Figure 1, the water tunnel was in a closed loop form, which had a main test section of $1.2\times 1.0\times 16{ \mathrm{m}}^3$ with a turbulence intensity of $0.3\%$. The walls of the test section were made of glass with good light transmission. The maximum velocity of the test section was approximately $0.366 \mathrm{m}/\mathrm{s}$ and the minimum velocity was $0.1 \mathrm{m}/\mathrm{s}.$ The testing water temperature was $20 ℃$.

The main part of the rotor system was a two-blade single rotor. The rotor was able to rotate freely in the clockwise direction (seen from above), and its rotational speed was recorded by a flange plate fixed on the rotation shaft that shaded a photoelectric switch fixed to the base. The photoelectric switch recorded the accurate time when the blade rotated to a specific phase position, which was an important reference of wake age. Each blade was untwisted and rigid and was connected rigidly to the hub. The rotor had a radius $\mathrm{R}=0.25 \mathrm{m}$, with a blade root cutout of $0.053 \mathrm{m}$ ($21\%$ of rotor radius). The blade had a rectangular planform and a chord of $0.025 \mathrm{m}$. The airfoil was NACA 0015. The surface of the blades was painted black for reducing reflection. In the rotor system, there were two variable parameters. One was the pitch, the other was the shaft backward angle. The shaft backward angle could be changed from $-60°$ to $60°$. The detailed rotor properties are shown in

Table 1. Properties of the test rotor.

Rotor Type	Hingeless
Number of blades of each rotor	2
Radius	0.25 m
Solidity	0.0512
Chord	0.025 m
Twist	0°
Root cutout	0.053 m
Airfoil	NACA0015

.

2.2. PIV Instrumentations

Figure 2 shows the entire PIV system and the rotor setup configuration. The complete system consisted of a pair of frequency-doubled Nd: YAG lasers, an optical arm, light-sheet optics, a digital CMOS camera of $2048\times 2048$ pixels resolution, a high-speed frame grabber, image acquisition and analysis software, and a computer. The fiber optic light arm transmitted the laser light from the dual Nd: YAG lasers to the rotor flow field. The lasers lighted an observation plane in the wake, by emitting a 1 mm thickness light sheet.

As shown in Figure 3, the longitudinal observation planes were determined by the shaft axes and the flow direction. The longitudinal plane contained the longitudinal window (LW). The observation planes obtained a spatial view of the flow field, which got a clear comprehension of the autorotation mechanism. The CMOS camera captured images of the flow field from a direction perpendicular to the observation plane lit by the laser sheet.

In Figure 3, ${\delta }_s$ means the shaft backward angle. It is a very important parameter of the autorotation that could influence the rotational speed greatly. The blade’s pitch angle $\mathsf{\theta }$ could also be changed. A positive $\theta$ means that the blade’s leading edge rotated downward direction to the bottom of the water tunnel, while a negative $\theta$ means that the blade leading edge is rotated in an upward direction toward the surface of the water tunnel. This is to say, the device was placed upside down in the water tunnel, which has the advantage of avoiding the flow feature being disturbed by the rotating shaft.

3. CFD Simulation Method

3.1. Governing Equation and Solution Method

The conserved Reynolds averaged Navier–Stokes equation [55] of three-dimensional, viscous, unsteady, and compressible flow in the inertial coordinate system can be obtained as:

$$\frac{\partial }{\partial \mathrm{t}}{\int \int \int }_V\overrightarrow{W}dV+{\int \int }_S\left({\overrightarrow{F}}_c-{\overrightarrow{F}}_v\right)·\overrightarrow{n}dS=0$$

(1)

where $\overrightarrow{W}={\left[\rho ,\rho u,\rho v,\rho w,\rho e\right]}^T$ is a conserved variable, ${\overrightarrow{F}}_c=\left(f,g,h\right)$ is the convective flux, and ${\overrightarrow{F}}_v=\left(a,b,c\right)$ is the viscous flux. The specific expression is:

$$f=\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ \rho uv\\ \rho uw\\ \rho ue+up\end{array}\right], g=\left[\begin{array}{c}\rho v\\ \rho uv\\ \rho v^2+p\\ \rho vw\\ \rho ve+vp\end{array}\right], h=\left[\begin{array}{c}\rho w\\ \rho uw\\ \rho vw\\ \rho w^2+p\\ \rho we+wp\end{array}\right]$$

(2)

$$\begin{gathered} a=\left[\begin{array}{c}0\\ {\tau }_{xx}\\ {\tau }_{yx}\\ {\tau }_{zx}\\ u{\tau }_{xx}+v{\tau }_{yx}+w{\tau }_{zx}-\dot{q_x}\end{array}\right] \\ b=\left[\begin{array}{c}0\\ {\tau }_{xy}\\ {\tau }_{yy}\\ {\tau }_{zy}\\ u{\tau }_{xy}+v{\tau }_{yy}+w{\tau }_{zy}-\dot{q_y}\end{array}\right] \\ c=\left[\begin{array}{c}0\\ {\tau }_{xz}\\ {\tau }_{yz}\\ {\tau }_{zz}\\ u{\tau }_{xz}+v{\tau }_{yz}+w{\tau }_{zz}-\dot{q_z}\end{array}\right] \end{gathered}$$

(3)

In the viscous flux, the stress can be expressed as follows:

$$\{\begin{array}{}{\tau }_{xx}=2\left({\mu }_l+{\mu }_t\right)\left[u-\frac{1}{3}\left(u+v+w\right)\right]\\ {\tau }_{yy}=2\left({\mu }_l+{\mu }_t\right)\left[v-\frac{1}{3}\left(u+v+w\right)\right]\\ {\tau }_{zz}=2\left({\mu }_l+{\mu }_t\right)\left[w-\frac{1}{3}\left(u+v+w\right)\right]\\ {\tau }_{xy}={\tau }_{yx}=\left({\mu }_l+{\mu }_t\right)\left(u+v\right)\\ {\tau }_{yz}={\tau }_{zy}=\left({\mu }_l+{\mu }_t\right)\left(v+w\right)\\ {\tau }_{zx}={\tau }_{xz}=\left({\mu }_l+{\mu }_t\right)\left(w+u\right)\end{array}$$

(4)

where, ${\mu }_l$ is laminar viscosity coefficient, ${\mu }_t$ is turbulent viscosity coefficient, $u$, $v, w$ is velocity in $x$, $y$, $z$ direction.

The heat flow in $x$, $y$, $z$ direction can be expressed as follows:

$$\{\begin{array}{c}{\dot{q}}_x=\frac{\gamma R}{\gamma -1}\left(\frac{{\mu }_l}{P_{rl}}+\frac{{\mu }_t}{P_{rt}}\right)\frac{\partial T}{\partial x}\\ {\dot{q}}_y=\frac{\gamma R}{\gamma -1}\left(\frac{{\mu }_l}{P_{rl}}+\frac{{\mu }_t}{P_{rt}}\right)\frac{\partial T}{\partial y}\\ {\dot{q}}_z=\frac{\gamma R}{\gamma -1}\left(\frac{{\mu }_l}{P_{rl}}+\frac{{\mu }_t}{P_{rt}}\right)\frac{\partial T}{\partial z}\end{array}$$

(5)

The relationship between temperature and pressure can be expressed as follows:

$$p=\rho \left(\gamma -1\right)\left[E-\left(u^2+v^2+w^2\right)/2\right]$$

(6)

$$T=p/\rho R$$

(7)

where, $V_r=V-V_t$, $V_t$ are, respectively, the normal and tangential velocities on the boundary, $n_x$, $n_y$,$n_z$ are the unit external normal vectors, $\rho$ is the fluid density, $p$ is the total pressure of the fluid, $E$ is the sum of the internal energy of the fluid per unit mass, $H$ is the total enthalpy per unit mass, T is the temperature, $\gamma$ is Specific heat capacity, $P_{rl}$ is the laminar Prandtl number, and $P_{rt}$ is the turbulent Prandtl number.

3.2. Turbulence Model

The Baldwin–Lomax (B-L) turbulence model [56] was adopted in this paper. In this model, the turbulent flow was divided according to the normal distance from the wall. The Prandtl–Van Driest model was used to calculate the inner region, and the Clauser model was used to calculate the outer region. The viscosity coefficient is as follows:

$${\mu }_t=\{\begin{array}{c}{\left({\mu }_t\right)}_{inner},y\le y_{crossover}\\ {\left({\mu }_t\right)}_{outer},y>y_{crossover}\end{array}$$

(8)

where y is the normal distance from the wall.

The judgment basis $y_{crossover}$ is expressed as:

$$y_{crossover}=\min \left\{y|{\left({\mu }_t\right)}_{inner}={\left({\mu }_t\right)}_{outer}\right\}$$

(9)

(1) Prandtl–Van Driest model

$${\left({\mu }_t\right)}_{inner}=\rho l^2\left|w\right|$$

(10)

where, $\left|w\right|$ is the magnitude of the vorticity.

The length $l$ can be expressed as:

$$l=ky\left[1-exp\left(-y^{+}/A^{+}\right)\right]$$

(11)

The dimensionless distance parameter $y^{+}$ is expressed as:

$$y^{+}=\frac{y\rho u_t}{\mu }=\frac{y}{v}\sqrt{\frac{{\tau }_w}{\rho }}$$

(12)

where, $u_t$ is the wall friction velocity, ${\tau }_w$ is the wall shear stress, $y$ is the wall distance, $k$ and $A^{+}$ are constant.

(2) Clauser model

$${\left({\mu }_t\right)}_{outer}=K{C}_{CP}\rho F_{WAKE}{F}_{KLEB}\left(y\right)$$

(13)

where K is the Clauser constant, $C_{CP}$ is an additional constant, and

$$F_{WAKE}=min\{y_{MAX}{F}_{MAX},C_{WK}{y}_{MAX}{u}_{DIF}{}^2/F_{MAX}\}$$

(14)

The quantities $y_{MAX}$ and $F_{MAX}$ are determined from the function

$$F\left(y\right)=y\left|w\right|\left[1-exp\left(-y^{+}/A^{+}\right)\right]$$

(15)

The quantity $F_{MAX}$ is the maximum value of $F\left(y\right)$ that occurs in a profile and $y_{MAX}$ is the value of $y$ at which it occurs.

The function $F_{KLEB}\left(y\right)$ is the Klebanoff intermittency factor given by

$$F_{KLEB}\left(y\right)={\left[1+5.5{\left(\frac{C_{KLED}}{y_{MAX}}\right)}^6\right]}^{-1}$$

(16)

The quantity $u_{DIF}$ is the difference between the maximum and minimum total velocity in the profile. $C_{WK}$ and $C_{KLED}$ are constants.

In this paper, the two-time advance method with higher computational efficiency is selected for time discretization, and the Runge–Kutta explicit solution method is used in the sub-iteration process.

Figure 4 shows the flow diagram of the dual time advance format using the motion nested grid. This time advance method is widely used in the unsteady CFD simulation of the rotor.

3.3. Boundary Condition

The boundary in the whole flow field can be divided into the far-field boundary and the surface boundary. In this paper, the flow field of the rotor was calculated using the moving overset grid method, choosing the velocity inlet as the inlet and pressure outlet as the outlet, with a non-slip grid used for the other boundaries. The non-slip grid was used on the rotor surface.

The thermodynamic boundary conditions are as follows:

$$\frac{\partial T}{\partial n}=0$$

(17)

The dynamic boundary conditions are as follows:

$$\frac{\partial p}{\partial n}=0$$

(18)

where $n$ is the normal gradient. In the numerical calculation, the pressure value on the blade surface is replaced by the pressure value at the center of the grid adjacent to the blade surface. The physical parameters: the medium is water, the temperature is 20 °C, density is 998.2 kg/m³, and kinematic viscosity is 1.0087 mPa∙s.

3.4. Dynamic Equation

In this paper, each blade with uniform density is evenly distributed around the rotor shaft, so the center of mass of the whole rotor is located at the center of the rotor shaft, which is the origin of the rotor disk coordinate system. According to Newton’s law of motion, the six-degree of freedom equation of the rotor in space is expressed as follows:

$$\begin{array}{}m\left(\frac{d{V}_x}{dt}+V_z{\omega }_y-V_y{\omega }_z\right)=F_x-mg\sin \theta \\ m\left(\frac{d{V}_y}{dt}+V_x{\omega }_z-V_z{\omega }_x\right)=F_y-mg\cos \theta \cos \gamma \\ m\left(\frac{d{V}_z}{dt}+V_y{\omega }_x-V_x{\omega }_y\right)=F_z-mg\cos \theta \sin \gamma \\ I_x\frac{d{\omega }_x}{dt}+\left(I_z-I_y\right){\omega }_z{\omega }_y+I_{xy}\left({\omega }_z{\omega }_x-\frac{d{\omega }_y}{dt}\right)=M_x\\ I_y\frac{d{\omega }_y}{dt}+\left(I_x-I_z\right){\omega }_x{\omega }_z+I_{xy}\left({\omega }_z{\omega }_y-\frac{d{\omega }_x}{dt}\right)=M_y\\ I_z\frac{d{\omega }_z}{dt}+\left(I_y-I_x\right){\omega }_y{\omega }_x+I_{xy}\left({\omega }_y{}^2-{\omega }_x{}^2\right)=M_z\end{array}$$

(19)

where $m$ is the blade mass; $\theta$ is the pitch angle; $\gamma$ is the roll angle; $I_x$, $I_y$, $I_z$ are the inertia moments; $I_{xy}$ is the inertia product; ${\omega }_x$, ${\omega }_y,$ ${\omega }_z$ are the three components of angular velocity; $V_x$, $V_y$, $V_z$ are the three components of water velocity; $F_x$, $F_y$, $F_z$ are the three components of the total aerodynamic force; and $M_x$, $M_y$, $M_z$ are the three components of the external moment. The convergence condition of the autorotation rotor is $M_z$= 0.

4. Grid System

In this paper, a structural overset grid system was developed, which can take into account the motion of the blade. The grid system is characterized by robustness and efficiency. The grid system consists of two kinds of grids: one is the blade body-fitted structural grid, and the second is the background structure grid, which belongs to the Cartesian grid and overlaps with the blade body-fitted grid to form a nested relationship.

The method of “plane to solid” was used to divide the blade’s grid. First, the blade was divided into grids on each blade element plane, and then the grids on blade element planes were combined to obtain a complete blade three-dimensional body-fitted grid. Finally, the algebraic method was used to generate the grids at the root and tip of the blade. Compared with the traditional elliptical grid generation method, the time required for mesh generation is reduced, although manual operations are increased. Furthermore, it can be extended to form high-quality meshes for blades with more complex shapes, which can demonstrate the tiny fluid flows on the blade surface.

As shown in Figure 5, the C-type mesh was selected as the blade body-fitted structure mesh. The two-dimensional cross section of the blade is distributed in an H-shape along the radial direction, resulting in a C-H blade body-fitted structure grid as shown in Figure 6. Due to the irregular shape of the blade tip, the flow field of autorotating rotors is inclined to have vortex dissipation. The advantage of using the H-shaped distribution is that the blade body-fitted structure grid has a high grid quality and is less capable of producing mesh distortion. In addition, this type of grid has high computational efficiency. In order to simulate blade tip vortex shedding, the mesh of the blade tip was encrypted in the radial direction.

There are two methods for the numerical CFD simulation of the rotor’s rotating flow field, a time multi-reference frame (MRF) model and a moving overset grid method.

In this paper, the moving overset grid method [57] was used to simulate the autorotating rotor. The flow field is divided into two parts: the motion region and the static region. The motion region is the blade grid. The 3D body-fitted grid of the blade is the moving grid. The static area is the background grid, as shown in Figure 7. In the process of grid nesting, it is necessary to dig holes and interpolate n the overlapping regions to ensure that the blade 3D body-fitted grid can be well-matched with the background grid near the blade.

5. Results and Discussion

In this section, the hydrodynamic characteristics of the underwater autorotating rotor are investigated. In the first three subsections of this section, the rotational speed, induced velocity, and tip vortex trajectory obtained from the test are compared with the CFD results to verify the accuracy of the CFD method. The fourth subsection calculates the thrust and thrust coefficients by the CFD method with different water velocities, shaft backward angles, blade pitches, blade twist angles, blade airfoils, and the number of blades.

5.1. Rotational Speed

The pitch $\theta$, water velocity $v$, and the shaft backward angle ${\delta }_s$ were changed in the experiments to investigate the rotational speed of the underwater autorotating rotor.

The variation of rotational speed with ${\delta }_s$ was shown in Figure 8. The rotational speed was linearly related to ${\delta }_s$.

The variation of rotational speed with $\theta$ is shown in Figure 9. A simple relationship between rotational speed and pitch is shown in Figure 8, first increasing and then decreasing, which was explained by the variable lift drag ratio with the angle of attack.

The variation of rotational speed with $v$ is shown in Figure 10. There is an apparent positive linear correlation between rotational speed and water velocity.

It can be seen from Figure 8, Figure 9 and Figure 10 that the rotational speeds obtained by CFD simulation are generally in agreement with the test data. As can be seen from

Table 2. Rotational speed from tests and CFD.

Condition No.	Pitch (°)	Shaft Backward Angle (°)	Test Speed (rpm)	CFD Speed (rpm)	Error
1	−3	40	104.2	100.5	−3.55%
2	0	34	98	93.2	−4.90%
3	0	40	109	113	3.67%
4	0	43	120	117	−2.50%

, the error is within 5%. As a result, the CFD simulations can be implemented for the calculation of the rotational speed.

5.2. Induced Velocity

In this paper, when v = 0.306 m/s, δ_s = 40°, and θ = −3°, the induced velocity of the LW plane was investigated.

The LW plane recorded the longitudinal view of the axial-induced velocity v_iz and radial-induced velocity v_ir. Figure 11 shows the axial induced velocity at ψ = 180°. The positive axial induced speed is the decreasing z direction.

Figure 12 shows the radial induced velocity at ψ = 180°. The positive velocity direction is toward the blade tip. ∆ψ indicates the rotation angle of the blade after the blade position coincides with ψ = 180°.

As can be seen from Figure 11, the axial induced velocity is around 0.05 m/s. It is noteworthy that the induced velocity varies considerably at the blade root and tip, which was influenced by the hub, shaft, and tip vortices. This is a remarkable feature that the induced velocity after the blade tip drops rapidly to zero, which was influenced by tip vortices.

As shown in Figure 11 and Figure 12, the induced velocities obtained from CFD simulations correspond to the test data. Consequently, the CFD simulation method can be performed to calculate the induced velocity.

There are some reasons for the variation between CFD results and experimental results:

(1) The stable rotational speed calculated by CFD has some errors with the experimental results;

(2) The CFD calculation has not considered the interference of the hub and shaft;

(3) The experimental results have very few bad points.

5.3. Tip Vortex Trajectory

In this paper, when v = 0.306 m/s, δ_s = 45°, and θ = −3°, the tip vortex trajectory was analyzed at ψ = 180°. The experimental results showed that the tip vortex trajectory at ψ = 0° was not significant, so it was not covered in this paper. Figure 13 illustrates the data provided by the experiments on tip vortex shedding trajectory at ψ = 180°.

Figure 14 and Figure 15 show the radial and axial displacements of the tip vortex, respectively. Both show a clear linear relationship with time. The results show a high degree of confidence that the CFD simulation data correspond to the test data. This demonstrates that the CFD simulations can be extended to the calculation of the tip vortex trajectory of the underwater autorotating rotor.

5.4. Thrust and Thrust Coefficient

5.4.1. Water Velocity

When δ_s = 40°, and θ = −3°, the thrust T and the thrust coefficient C_T was calculated by CFD simulation with water velocity v. As mentioned in the previous chapter, the rotational speed increases with the water velocity increase. Figure 16 shows the thrust and thrust coefficient with various water velocities. The thrust has a linear positive correlation with water velocity, while the thrust coefficient has a linear negative correlation with water velocity.

As the water speed increases, assuming that the rotational speed and induced velocity remain constant, the angle of attack of the blade will increase in order to achieve zero torque equilibrium conditions, which results in the lift component in the rotor disc exceeding the drag component. Thus, this results in an increase in the rotor driving region, which leads to increasing rotational speed and thrust.

The pressure distribution on the blade upper surface at various water velocities is illustrated in Figure 17. As the water velocity increases, the negative pressure area on the blade’s upper surface expands, while the positive pressure at the blade’s leading edge also appears to increase. The maximum and minimum negative pressure peaks on the blade’s upper surface occur at ψ = 135° and 315°, respectively. The smaller the peak negative pressure on the blade’s upper surface, the more favorable it is to delay the occurrence of air separation. As a result, when the water velocity exceeds a critical value, airflow separation will occur on the blade’s upper surface, and thus a stall zone will appear at the blade tip, and the thrust coefficient will drop rapidly.

5.4.2. Shaft Backward Angle

When v = 0.366 m/s and θ = 0°, the thrust T and thrust coefficient C_T was calculated by CFD simulation with the shaft backward angle δ_s. Figure 18 show the thrust and thrust coefficient with various shaft backward angles. As can be clearly seen in Figure 18 the thrust force is linearly and positively related to the shaft backward angle. However, the thrust coefficient scarcely varies with the shaft backward angle.

Figure 19 shows the force analysis of the blade profile. The increase in the shaft backward angle causes the vertical velocity U_p to rise. The tangential component of water velocity in the rotor disc is negligible relative to the blade’s rotational linear velocity, so that U_T is essentially close to the blade’s rotational linear velocity. Thus, at a fixed rotational speed, the angle of attack α increases as the shaft backward angle increases, resulting in an increase in the lift dL and drag dD. The increase in angle of attack α also leads to an increase in the lift-to-drag ratio, so that the increment of the lift component dL is greater than the increment of the drag component dD in the horizontal direction, which in turn generates a driving force that increases the rotor speed. An increase in rotor speed leads to a decrease in the angle of attack α. So, the increase in the vertical component of the total force dR causes an increase in the thrust. According to CFD simulation data, the rotational speed increases as the shaft backward angle increases, but the ratio U_p/ΩR remains essentially unchanged. This indicates that the angle of attack α remains essentially the same and so the thrust coefficient does not respond essentially to variation in the shaft backward angle.

5.4.3. Blade Pitch

When v = 0.366 m/s and δ_s = 40°, the thrust T and thrust coefficient CT was calculated by CFD simulation with θ. Figure 20 shows the thrust and thrust coefficient with various pitches. From Figure 20, we can see clearly that the thrust and thrust coefficient have a linear positive correlation with the pitch.

With the increase in the pitch, the inflow angle ϕ increases, which leads to the increase in the vertical component of the resultant force dR and the increase in the thrust. Because the ratio of U_p/ΩR increases, the angle of attack α of increases leading to an increase in the thrust coefficient.

As the pitch increases, the inflow angle ϕ increases, which leads to an increase in the vertical component of the resultant force dR. Ultimately, the thrust force increases. As the ratio U_p/ΩR increases, an increase in the angle of attack α induces an increase in the thrust coefficient.

5.4.4. Blade Twist Angle

The blade torsion angle is one of the most important geometrical features of a rotor and it has an influential effect on the hydrodynamic characteristics of underwater autorotating rotors. In this paper, the hydrodynamic characteristics with blade torsion angles of −8°, −4°, 0°, 4°, and 8° are investigated.

The variation of rotational speed with blade twist angles is shown in Figure 21. Within a certain range, the rotational speed under negative torsion is less than that without torsion.

When v = 0.366 m/s, θ = 0°, and δ_s = 40°, the thrust T and thrust coefficient C_T was calculated by CFD simulation with the blade twist angle. Figure 22 shows the thrust and thrust coefficient with various blade twist angle. The thrust has a positive correlation with the blade twist angle, while the thrust coefficient has a negative correlation with the blade twist angle.

Within a certain range, the angle of attack at the blade root decreases due to positive torsion, reducing the stall region. The blade tip lift-to-drag ratio increases, the driving region expands outward and the driven region decreases, thus increasing the rotational speed. However, when the speed exceeds a specific value, the blade tip stalls, resulting in a reduction in rotational speed. Due to the fact that the driving region makes up the majority of the rotor disc, the angle of attack in this region decreases as the blade torsion angle increases, so the thrust coefficient decreases gradually. However, excessive positive twisting of the blade can lead to a decrease in thrust and an increase in thrust coefficient.

5.4.5. Blade Airfoil

When v = 0.366 m/s, θ = 0°, and δ_s = 40°, the thrust T and thrust coefficient C_T with different blade airfoils was calculated by CFD simulation with various blade airfoils.

Table 3. The thrust and thrust coefficient with various blade airfoils.

Condition No.	Blade Airfoil	Rotational Speed (rpm)	Thrust (N)	Thrust Coefficient (×10−3)
1	NACA0009	131.1	8.35	7.23
2	NACA0012	125.6	7.97	7.53
3	NACA0015	113	7.05	8.22

shows the thrust and thrust coefficient with various blade airfoils.

This paper investigates the effect of a symmetrical airfoil on the hydrodynamic characteristics of underwater autorotating rotors. As the relative thickness of the blade airfoil increases, the speed and thrust decrease but the thrust coefficient increases. This is because, as the relative thickness of the blade airfoil increases, the lift-to-drag ratio decreases at low Reynolds numbers. This situation causes the zero-torque position of the rotor disc to move outward, resulting in an enlarged area of the driving region and an enhanced rotational speed.

The pressure distribution on the upper surface of the advancing blade at 0.9R is shown in Figure 23. The area of negative pressure distribution on both the upper and lower surfaces of the blade increases with increasing blade relative thickness, but the negative pressure on the lower surface increases more noticeably. At the trailing edge of the blade, the greatest positive pressure is found on the NACA0012 airfoil.

5.4.6. Number of Blades

When v = 0.366 m/s, θ = 0°, and δ_s = 40°, the thrust T and thrust coefficient C_T was calculated by CFD simulation with various numbers of blades.

Table 4. The thrust and thrust coefficient with various numbers of blades.

Condition No.	Number of Blades	Rotational Speed (rpm)	Thrust (N)	Thrust Coefficient (×10−3)
1	2	113	7.05	8.22
2	3	97.4	8.24	12.93
3	4	90.2	9.47	17.32

shows the thrust and thrust coefficient with various numbers of blades.

We can see from Table 4 that an increase in the number of blades results in a decrease in rotational speed but an increase in thrust and thrust coefficient. However, the aerodynamic interference between the blades increases with the number of blades. When the interference between blade aerodynamics is too large, the increase in the number of blades fails to enlarge the thrust and thrust coefficient.

The autorotation of an underwater rotor is the conversion of kinetic energy between the water and the rotor. In this process, the loss of kinetic energy of the water is limited. Therefore, as the number of blades increases, the amount of energy gained by each blade decreases, causing a reduction in the rotational speed. The increase in the number of blades means that the rotor solidity increases, which makes the thrust coefficient increase. The kinetic energy of the water will not be able to drive the rotor rotation when the rotor solidity reaches a critical value.

6. Conclusions

In this paper, PIV tests and CFD simulations were carried out to investigate the hydrodynamic characteristics of underwater autorotating rotors. The underwater autorotating rotor was tested and investigated in terms of three aspects: rotational speed, induced velocity, and tip vortex trajectory, and the test results were used to validate the accuracy of the CFD simulations. The CFD simulations were used to calculate the thrust and thrust coefficient of the underwater autorotating rotor in terms of six aspects: water velocity, shaft back angle, pitch, twist angle, and blade numbers. The following conclusions were drawn from this work.

• The experimental results show that the rotational speed has a significant positive correlation with the water velocity and shaft backward angle, but the relationship with the pitch is first increasing and then decreasing. The axial-induced velocity is around 0.05 m/s. The induced velocity fluctuates sharply around the blade tip and rapidly drops to zero outside the blade tip. The dramatic variations are due to the effect of the water velocity and the blade tip vortices. The radial and axial displacements of blade tip vortex trajectories are clearly linear with respect to time.
• A computational fluid dynamic (CFD) based on moving overset grids was developed to study the hydrodynamic characteristics of the underwater autorotating rotor. In this simulation, the Navier–Stokes equations were solved using the overset grids technique to calculate the flow field of the underwater autorotating rotor under various states. The experimental results demonstrate that CFD simulation method is suitable for investigating the hydrodynamic characteristics of underwater autorotating rotors. Induced velocity verifies the accuracy of the simulated velocity flow field. The position of the blade tip vortex trajectory confirms the accuracy of the simulated blade tip complex flow field.
• The thrust has a linear positive correlation with water velocity, but the thrust coefficient has a linear negative correlation with water velocity. When the water speed exceeds a certain value, waterflow separation occurs on the blade upper surface resulting in the emergence of a stall region at blade tip, which causes the thrust coefficient to drop rapidly.
• The thrust is linearly related to the shaft backward angle, but the thrust coefficient is almost a fixed value with δ_s. However, an excessively large shaft backward angle increases drag in the forward direction, so it is essential to select a properly angled shaft back. The thrust and thrust coefficient are linearly and positively correlated with the blade pitch as the blade pitch increases from a negative value to zero.
• As the blade twist angle changes from negative to positive, the thrust first increases and then tends to stable, in contrast to the decrease in the thrust coefficient. Therefore, a suitable negative blade torsion is more advantageous for underwater autorotating rotors. However, the excessive positive twist of the blade can lead to a decrease in thrust and an increase in thrust coefficient. For symmetrical airfoils, within a certain range, as the airfoil thickness increases, the rotational speed, and thrust decrease, but the thrust coefficient increases. Appropriately thin airfoils are more beneficial for underwater autorotating rotors. An increase in the number of blades increases thrust and thrust coefficient but reduces rotational speed. The rotor will not rotate when the number of blades reaches a specific value. When the rotor can be rotated, the number of blades can be increased appropriately to improve the performance of underwater autorotating rotors.
• All our preliminary results throw light on the fundamental characteristics of underwater autorotating rotors. A limitation of this study is that the blade deformation is not considered in the CFD simulations. Further research should be carried out on multiple blades to enrich the theory of underwater autorotating rotors.