1 Introduction

Helicopters are very important aerial vehicles in civil services such as aerial photography, fire-fighting, and delivering medical aid because they can take off vertically and fly forward and backward. Helicopters are also used in military tasks, such as transporting supplies and escorting wounded individuals. During flight over desert surfaces, a helicopter will accumulate electrostatic charge on its body due to friction with suspended ionized particles. Charge transfer is enhanced when the atmosphere contains dust and sand particles due to triboelectric effects as the volume fraction of particles increases^[1]. These flight conditions are further complicated by the presence of particles because the helicopter interacts with the turbulent wake generated by the particles. Therefore, the hydrodynamics of gas-particle mixtures resulting from revolving rotors and the effects of hovering altitude on the performance of helicopters should be investigated.

The flow behavior between the helicopter and a building obstacle was investigated within the aerodynamic wake of a rotor^[2]. Experimental gas circulation measurements were gathered using a particle image velocimetry system, and the velocity field was obtained between the helicopter and the wall obstacle at the ship deck^[3]. Gas circulation between the helicopter and the vertical wall was measured using the particle image velocimetry system at the ship deck^[4]. The flow behavior of downwash and outwash gas from the rotor was measured as the helicopter hovered near a land-based hangar^[5]. The influence of the downwash airflow was simulated using the sliding mesh model in order to study variations in the forces at the airfoil surface^[6]. A comparison between the simulated results and the experimental data confirmed that the numerical computation model can be used to effectively analyze the aerodynamics and design of helicopter rotors.

Despite the number of numerical simulations and experimental measurements reviewed in the literature^[7], a systematic investigation of the dynamic phenomena, to the best of the authors' knowledge, has not been done. Present numerical simulations attempt to investigate the effect of the ground on the flow dynamics of gas and particles induced by the low-level rotary-winged aircraft hovering over a particle bed. In this study, the flow field of particles is simulated in the presence of a rotor, and collisional interactions at the particle-blade surface and between particles are modeled using the kinetic theory of granular flow (KTGF). The aerodynamics around the rotor is numerically modeled using the multiple reference frame (MRF) method. The effects of the hovering altitude on the velocity and volume fraction of gas and particles are analyzed along the particle bed. The simulated gas pressure coefficient is compared with that in the experimental measurements presented in the literature.

2 Two-fluid model of gas and particle flow induced by rotor

The hydrodynamics of gas and particles in a real helicopter is very complex as different body components, such as the fuselage, main rotor, and tail rotor, influence the aerodynamics. An understanding of the influence of these components on flow behavior can be studied separately, and useful information can be obtained for use in flight simulators. Thus, the present numerical simulations will only focus on the main rotor. A schematic diagram used in numerical simulations is shown in Fig. 1, which is based on experimental rotational wing facilities^[8]. The rotor diameter is 2.286 m. The chord is 0.190~5 m, and the aspect ratio is 6. The attack angle is set as 8°.

2.1 Governing equations of fluid-solid phase flow

The present study aims to investigate the flow phenomena of a gas and particles from a particle bed induced by rotation of the main rotor. Therefore, the governing equations used must be able to describe the motion of the gas and particles, as well as variations in the parameters of the gas and particles. Liu and Pan^[6] simulated the rotational motion of rotor with an NACA0012 airfoil profile using the k-ε turbulence model. As a result, the two equations in the k_g-ε_g turbulence model and KTGF are selected to calculate the gas turbulence and collisional interactions. The mass conservation equation of the gas and solid phases is expressed as follows:

(1)

where u is the velocity, and α and ρ are the volume fraction and the density, respectively. Conservation of momentum for gas and solid phases can be defined as follows:

(2)

(3)

where the subscripts g and s denote the gas and solid phases, respectively, and Ω and R are the angular spin rate and the rotation radius, respectively. The stress tensors for the gas and solid phases are

(4)

(5)

where μ_eg and μ_s are the effective gas viscosity and solid viscosity, respectively. λ_s is the bulk viscosity of the particles. The effective viscosity of the gas phase is the sum of the gas dynamic viscosity μ_g and the gas turbulent viscosity μ_gt which is expressed as μ_gt=C_μρ_gk_g²/ε_g. The gas turbulent kinetic energy k_g and its dissipation rate ε_g are modeled using the k_g-ε_g turbulence model^[9],

(6)

(7)

where G_kg is the turbulent kinetic energy generation due to the velocity gradients and can be expressed as follows:

(8)

and C₁, C₂, C_μ, σ_g, and σ_ε are empirical constants equal to 1.44, 1.92, 0.09, 1.0, and 1.3, respectively.

The solid shear viscosity includes kinetic and collisional contributions due to inelastic collisions, as well as a frictional contribution due to frictional interactions with the particles. This can be expressed as follows^[9]:

(9)

where d_s and e are the particle diameter and the coefficient of restitution of the particles, respectively. g₀ is the radial distribution function during contact. I_2D is the second invariant of the deviatoric stress tensor. The internal friction angle of the particles ϕ is 28.5°. The solid pressure is expressed as follows:

(10)

(11)

(12)

where the maximum solid volume fraction due to packing α_{s, max} is 0.635. θ is the granular temperature, defined as θ=C²/3, where C defines the velocity fluctuations of the particles^[10]. The granular temperature is modeled using the KTGF as follows:

(13)

where the collision energy dissipation γ_s and the diffusivity of granular temperature k_s are modeled using the KTGF^[10].

Gas-solid drag forces are modeled by considering the interface momentum transfer coefficient β and the slip velocity (u_g-u_s). In the present numerical simulations, the Huilin-Gidaspow model is used to calculate the interface momentum transfer coefficient in ANSYS^[11-12],

(14)

(15)

where C_D is the drag coefficient of the particles.

2.2 Boundary conditions and simulation strategy

The flow field of gas and particles induced by revolving rotors can be predicted using the MRF model, and denser unstructured grids are deployed near the blades. The computational domain has a length of 12 m, and width and height of 12 m and 5 m, respectively. This model is shown with the central point of the rotor hub in Fig. 2. Liu and Pan^[6] used a cylindrical calculation domain with a radius of 6.0 m to simulate the same type of blade. The simulated results show that the computational domain is sufficient to fully capture the influence range of the blades, and the predicted gas pressures agree with the experiments.

Table 1 shows the simulation parameters used in the three-dimensional rotational wing flow simulation. The particle bed is located at the bottom with a particle diameter and a density of 50 μm and 2 000 kg/m³, respectively. The particle bed is initially filled with particles. The initial bed height H₀ and the solid volume fraction are 0.5 m and 0.55, respectively. The initial velocity of the particles is zero. Pressure inlet boundary conditions are defined at the top inlet, and pressure outlet boundary conditions are defined at the domain sides. A wall boundary condition is defined at the bottom of the particle bed. Periodic boundary conditions are defined at the interface of the particle bed.

The hub has a solid boundary as it rotates. The gas phase has a no-slip boundary condition, and the velocity of the fluid along the boundary is equal to the velocity of the hub. The solid phase has a slip boundary condition at the surface of the hub^[13],

(16)

(17)

where n is the component normal to the wall, g_o is the radial distribution function, and u_sw and θ_w are the particle velocity in the direction parallel to the wall and the granular temperature at the wall, respectively. e_w is the wall-particle coefficient of restitution. The wall specularity coefficient ϕ is 0.1.

The computational domain is divided into an internal grid block with a rotational wing and an external grid block. The internal grid block revolves when the hub rotates along its axis. The MRF method is used to simulate the rotation of the hub. Node discontinuity arises between the internal and external grid blocks. Thus, the interface boundary condition is assigned at this location, and the physical quantities of the block are transferred using interpolation through the boundary with the other block. Thus, the interface between the two blocks prevents discontinuities in physical quantities across the bordering nodes. The numerical simulations are conducted using FLUENT^[11].

3 Simulations and discussion 3.1 Comparison between simulation and experimental results

The distribution of the pressure coefficients is calculated along the span of the blade from the suction side and the pressure side of the blade surface using the simulated gas pressure. Figure 3 shows the pressure coefficient distribution along the dimensionless distance at a span r/R of 0.5, where r is the radial distance from the center of the rotor, and R is the radius of the rotor blade. The dimensionless distance is defined as x/c, where x is the distance from the leading edge to the trailing edge, and c is the chord length of the blade. 128 556 coarse grids and 266 200 finer grids are tested in order to study the impact of the grid size on the results. The simulated pressure coefficients using coarse grids and finer grids exhibit similar trends when compared with the measured pressure coefficients along the chord length. The simulation results using the finer grid are consistent with the experimental measurements, indicating that a finer grid facilitates a more accurate prediction of the gas pressure along the curved structure of the surface of the blade. The simulated and experimental results show that the force on the pressure side is greater than that on the suction side, resulting in a fluid pressure difference between the two sides of the blade surface, which eventually provides upward lift.

3.2 Instantaneous velocity and volume fractions

Variations in the velocity and volume fraction of gas are predicted from the numerical simulations as the hovering altitude of the rotor changes. The streamlines of the gas axial velocity are shown in Fig. 4 at two hovering altitudes. The rotor is set to revolve at a constant rotation speed, and the gas pressure is received on both sides of the hub surface, thus providing lift. Gas flows from the inlet at the top, passes through the hubs, and moves toward the particle bed. When gas impinges on the surface of the particle bed, it moves horizontally along the particle bed surface as the particle bed obstructs the vertical flow. It is found that the axial velocities of the gas are low at the hub and high near the tip along the lateral direction Z. A symmetric gas vortex forms at the center of the rotor, and gas recirculates beneath the blades. This flow behavior is characterized as a symmetric vortex at the front of the rotor. The negative axial velocity indicates that gas flows from the top to the particle bed. The downward flow of gas produces an interaction force between the gas phase and the particles located at the surface of the particle bed. It is also found that the gas axial velocity is larger at a 1.5 m hovering altitude than that at a 0.75 m hovering altitude at the same rotation speed. Thus, the lifting force provided by the rotor increases as the hovering altitude decreases.

Variations in the volume fraction and velocity vector of the particles are shown in Fig. 5 at two different hovering altitudes. Particles are entrained by the gas flow along the horizontal direction. The shear stress accumulates in the gas as it makes contact with the particles at the surface of the particle bed, which changes their configuration. Particles obtain energy from the gas phase, which increases the drag force between the gas and particles. Finally, particles located at the top of the particle bed move horizontally along the surface of the particle bed. Meanwhile, particles within the particle bed maintain their positions due to large contact forces between adjacent particles. It is also found that the high velocity gas stream impacts on the particles located at the particle bed surface. The volume fraction of particles is lower at a greater distance from the surface of the particle bed. Thus, the blades suffer from collisions with particles, which influences the dynamics of the rotor.

Figures 4 and 5 show that the gas and particles have high velocities in the jet-wake regime at both hovering altitudes. The flow of the gas-particle mixture is similar to an impinging gas jet, which impinges on the surface of the particle bed and fans out radially along the top surface of the particle bed. Experimental measurements show that the particles acquire much higher velocity when the helicopter travels over the ground^[1]. The disturbance of the soil has a dramatic impact on the strength of particle entrainment because particles are removed from the top surface. The simulation reveals the aerodynamic behavior of the gas and particles at the top surface of the particle bed.

The distributions of surface gas pressure and velocity streamlines are shown in Fig. 6 at the upper and lower surfaces of the blade at two different hovering altitudes. Both cases show that the surface gas pressure is high at the root and low at the tip of the blade. Gas will flow from the root to the tip of the blade. The surface gas pressure at the suction side of the blade is large at the leading edge and low at the trailing edge. The results also show that the surface gas pressure increases from the trailing edge to the leading edge. The distribution of the surface gas pressure at the suction side and pressure side produces a pressure gradient when the attack angle is at 8° and the rotational velocity is 1 000 r/min. In contrast, the surface gas pressure undergoes a series of significant changes along the span of the blade, indicating that significant variations in the surface gas pressure and pressure gradient occur at the suction side of the rotor. When the surface gas pressures at different hovering altitudes are compared, the effect of the hovering altitude on the surface gas pressure is minor on both sides of the blade.

3.3 Distribution of velocity and volume fractions

The axial velocity distributions of the gas and particles are shown in Fig. 7 along the lateral position Z/R at different hovering altitudes, where Y is the distance from the center of the rotor. The gas and particle axial velocities are low at the hub, increase along the lateral distance, and eventually reach a maximum. They subsequently decrease at the tip of the rotor. The axial velocity of the gas is larger than the axial velocity of the particles. This indicates that the tip of the blade provides a large relative velocity between the gas and particles, and the drag force increases along the lateral position of the blade.

The gas turbulent kinetic energy induced by the rotation of the blades is predicted using the k_g-ε_g turbulence model. The turbulent kinetic energy and dissipation rate distributions of the gas are shown in Fig. 8 at two hovering altitudes. Roughly, both the gas kinetic energy and the dissipation rate increase toward the tip of the blade because the rotation of the rotor provides energy to the gas. Simulations showed that the blade tip vortex dominates the flow field as the helicopter hovers^[14]. They found that the kinetic energy of air is low in the middle portion of the rotor blade, which is similar to the simulated kinetic energy distribution. It is also found that the gas kinetic energy and dissipation rate decrease as the hovering altitude of the rotor increases. This indicates that the distance between the rotor and the surface of the particle bed affects turbulent kinetic energy transfer in the gas, thus changing the gas dynamics of the rotor.

The total force component acting on a blade wall zone is computed by summing the dot product of the pressure and viscous forces with the specified force vector. This means that the total force component consists of a pressure force component and a viscous force component in the blade wall zone. The total force component at the surface of the blade is shown in Fig. 9 as a function of hovering altitude. The total force component F_z acting on the blade surface along the z-direction is smallest, while the total force component F_y along the vertical direction is the largest at four different hovering altitudes. The differences in the total force component are obvious, and the lift force on the helicopter varies as the hovering altitude changes. Using a six-component balance method, the ground effect (characterized by the thrust coefficient) was measured with respect to the obstacle^[15]. The measured results show that the gas recirculates when the rotor is close to the obstacle. When the helicopter is close to the ground, the thrust coefficient of the blade increases, indicating that the loads on the rotor increase. This shows that the height from the ground affects the force on the blade, which is consistent with the simulation results.

4 Conclusions

The aerodynamics of gas and particles is predicted with the Euler-Euler two-fluid model in a revolving rotor over a particle bed. The collisional interactions at the particle-blade wall and between particles are modeled using the KTGF. The flow field of a revolving rotor is simulated using the MRF method. The simulated pressure coefficient is compared with that in the experimental measurements.

The velocity, volume fractions, and gas pressure distributions are predicted as the hovering altitude is varied. The simulated gas pressure at the pressure side of the blade is greater than that at the suction side, thus producing lift. The largest rotation speed occurs at the tip of the blades. The gas turbulent kinetic energy and dissipation rate are high at the blade tip. The rotor induces the horizontal flow of the gas and particles as the rotational speed of the rotor increases. The solid volume fraction decreases as the hovering altitude increases.

Although there are some discrepancies between the numerical results and the experimental results in terms of the pressure coefficient, the results show the same trends in curve development, which further confirms the reliability of the numerical computation when applied to a helicopter.