Abstract
A computational method is developed in order to predict the unsteady aerodynamic characteristics of the tilt-rotor aircraft in conversion mode. In this approach, the rotor is modeled as an actuator disk so that the effect of individual blades can be ignored. A novel predictor-corrector-based dynamic mesh method is presented for dealing with extremely large mesh deformation during a conversion process. The dual time-stepping approach and the finite volume scheme are applied to solve the unsteady N-S equation. A parallel algorithm is utilized in this work to improve the computational efficiency. By using the present method, quantitative and qualitative comparisons are made between the aerodynamic coefficients obtained in the quasi-steady fixed conversion mode and the time-accurate continuous transition flight condition. Both two-dimensional (2D) and three-dimensional (3D) computations are carried out. The influence of the tilt modes and the tilt period time on the unsteady aerodynamic forces are also studied. Numerical results demonstrate that the developed method is effective in simulating the aerodynamic characteristics of the tilt-rotor aircraft in conversion mode.
1. Introduction
The tilt-rotor aircraft combines the high-speed efficiency of a turboprop aircraft with the vertical takeoff and landing capabilities of a helicopter. It has a bright future in military and civilian areas [1–4]. During the past several decades, extensive researches have been conducted on its unique aerodynamic characteristics and technical complexity. In contrast to time-consuming experimental investigations, the computational fluid dynamics (CFD) method has become an important tool to predict the aerodynamic force in industrial practices.
In the framework of tilt-rotor aircraft design, a major challenge is how to account for the flow unsteadiness due to the rotating blades. Although there exist some low-fidelity rotor modeling tools, such as the blade element method and the panel method, they are not suited for an accurate computation of the viscous interaction between the rotor and fuselage. Such analysis demands high-fidelity methods, i.e., methods based on a solution of the Euler or N-S equations that compute both vortices and the pressure field accurately. In general, most computational studies conducted on modeling the rotor can be broadly classified into two categories. The first approach is to simulate the unsteady flow field by taking into account each individual blade. This is usually accomplished by using the overset or Chimera grid technique. Some notable works include those of Yoon et al. [5], Potsdam and Strawn [6], and Schwarz [7]. Although this approach has high fidelity, it is time consuming. Drastically reduced geometrical complexity and computational time can be achieved with the second approach, i.e., modeling the flow field by employing an actuator disk. In this approach, the rotor is represented as an infinitely thin disk capable of sustaining a pressure discontinuity. The actuator disk is essentially a limiting case in which the number of blades goes to infinity. Since helicopters possess a finite number of blades, the actuator disk assumption provides a time-averaged representation of the flow. This reduces the cost of an unsteady simulation down to a stationary one. An early actuator disk implementation was performed by Whitfield and Jameson [8] for propeller simulations. Barwey et al. [9] successfully employed this approach to simulate helicopter rotors in forward flight and hover. Renaud et al. [10] presented a comparison of different actuator disk CFD codes on both structured and unstructured grids. O’Brien and Smith [11] studied the rotor-fuselage interactions using various rotor models, and the computational results seemed to agree well with wind tunnel data.
A unique feature of the tilt-rotor aircraft is the existence of conversion mode. The transition from helicopter mode to airplane mode involves a change in orientation of the rotor systems, mounted at the wing tips, from a horizontal plane to a vertical plane. The rotors which predominantly provide the lifting force in helicopter mode gradually transitioned to provide aircraft thrust in airplane mode. This highly coupled aerodynamic interaction plays a major role in determining the aerodynamic characteristics and performance. Research on the tilt-rotor aircraft in conversion mode is still a challenging issue, although extensive research works on the aircraft in helicopter mode or in airplane mode have been conducted. For example, Ye et al. [12] calculated the performance of the V-22 in helicopter mode by using the moving overset gird methodology. Both the isolated rotor and the full configuration aircraft are considered, and the computed results have good agreements with the experimental data. The flow field around a quad tilt rotor vehicle was performed by Gupta and Baeder [13] by using a compressible solver. They modeled the rotor as an actuator disk and studied the hover characteristics in ground effect and out of ground effect. Abras and Narducci [14] performed independent OVERFLOW and FUN3D CFD analyses of the MV-22 in airplane mode over a range of angles of attack and then compared the numerical results with the data from a high-angle-of-attack wind tunnel test. The flow over a wing-fuselage-nacelle configuration in forward flight was simulated by Tai [15] with a multizone, thin-layer N-S method. Major flow features such as the three-dimensional flow separation due to viscous-vortex interactions observed experimentally were captured.
During the conversion flight condition, it is difficult to calculate the flow field and aerodynamic force of the full tilt-rotor aircraft configuration due to complicated rotor-wing-fuselage interactions. Limited studies were conducted for modeling the conversion mode. At present, there are mainly two kinds of numerical methods. One is the quasi-steady approach. For example, Zhang et al. [16] presented a method combining a rotor actuator disk model and embedded grid technique to predict the aerodynamic characteristics of the tilt-rotor aircraft in conversion mode, but only four different tilt angles are calculated. Sheng and Narramore [17] calculated the aerodynamic force of a quad tilt-rotor in a fixed conversion mode at the tilt angle of 75°. The quasi-steady method cannot consider unsteady effects, since the rotor and nacelle are fixed at a certain rotation angle. The second approach is the moving overset mesh technique. Although it can deal with large deformation, the so-called “hole cutting” and “donor searching” need to be used during grid assembling, which may decrease the efficiency and accuracy [18]. For example, Li et al. [19] used a new multilayer moving-embedded grid technique to predict the aerodynamic force of a rotor in conversion mode, where the wing and fuselage are not modeled. However, it still requires large computational costs.
In terms of the grid type, multiblock structured grid is perfectly suitable for numerical simulations of viscous flows around complex tilt-rotor configurations. However, the rotor, typically located less than a wing chord above the wing, increases the difficulty of dynamic mesh generation in conversion mode. In the previous work, we found that most existing dynamic mesh methods, including inverse distance weighting (IDW) [20] and radial basis function (RBF) [21], may not efficiently handle this extreme situation. Therefore, it is urgent to develop an effective dynamic mesh generation method for extremely large deformation, that is, the premise of unsteady calculation of the tilt-rotor aircraft in conversion mode.
In this work, we perform numerical calculations on the unsteady aerodynamic force of the tilt-rotor aircraft in conversion mode by using a CFD approach. The structure of the paper is as follows: Section 2 presents the multiblock gird topology and static mesh generation. A novel predictor-corrector-based dynamic mesh method for extremely large deformation is also introduced in this section. Section 3 describes the basic CFD solver and the used actuator disk approach. Section 4 shows the numerical results for the tilt-rotor aircraft in conversion mode. Particularly, the influences of tilt modes and tilt period time on the unsteady aerodynamic forces are studied. Finally, conclusions and future work are given in Section 5.
2. Static and Dynamic Multiblock Grid Generation
A simplified tilt-rotor configuration as shown in Figure 1 is used for numerical simulation. The fuselage length is 16.35 m, the wing semispan is 7.38 m, and the rotor radius is 5.23 m. The nacelle and empennage are currently not modeled. The tilt angle has a value of 0° in helicopter mode and a value of 90° in airplane mode. In the present work, the rotor is not modeled as a set of blades rotating in the shaft plane, but instead modeled as an actuator disk surface.
(a) Helicopter mode
(b) Airplane mode
A reasonable multiblock grid topology is elaborately designed by making a comprehensive analysis on aerodynamic configuration and characteristics of unsteady motion, as shown in Figure 2. On the whole, there are three separate zones: (1) a single block with the size of , which directly extends from the wingtip to the far field boundary, and (2) an inner O-H type block with the size of . The inner boundary of this block is the fuselage and wing surface, and the outer boundary is the actuator disk surface. The O-type mesh around the wall section is generated at various spanwise locations and then joined to each other to get a three-dimensional mesh. Fine spacing (10-6 chord) in a direction normal to the body surface is used in order to better capture the boundary layer. The third separate zone is (3) an outer grid enveloping the actuator disk surface and extending to a circular outer boundary approximately 30 wing chords from the wing surface. The size of the grids is . The axial spacing on the two sides of the actuator disk is approximately 0.004 times the radius and increases as a geometric progression in the axial direction away from the disk. The total number of grid cells used for CFD computation is 6,256,000. Figure 3 shows the surface mesh of fuselage, wing, and symmetry plane. A perspective cutaway view of portions of the three-dimensional grids in helicopter mode is given in Figure 4.
In the transition corridor, the actuator disk rotates 90 degrees about the stationary wing. Available dynamic mesh technologies, such as IDW and RBF, are not capable of treating this extremely large rotation problem. For a specific rigid body motion, if dynamic meshes are known at multiple dynamic positions beforehand, the large deformation will be transformed into some medium and small deformations between neighboring positions. Inspired by this new idea, a novel predictor-corrector-based dynamic mesh method is developed in this paper, which includes the following three steps.
2.1. Generating Multiblock Grids at Multiple Dynamic Positions
According to the prescribed motion rule of the actuator disk, the instantaneous geometric configuration at every time step is exactly evaluated. Different from the traditional methods which employ only one static mesh, this study is intended to obtain the dynamic meshes from multiple sets of pregenerated grids. For this case, three typical dynamic positions are selected at the tilt angles of 0°, 45°, and 90°, where multiblock grids can be in advance generated using the same topology. The resulted grid coordinates are denoted by , , and , respectively. Figures 5 and 6 show the static meshes at and 90°, respectively.
In order to describe the pregenerated computational mesh more clearly, Figures 7–9 show the sectional grids near the wingtip at different tilt angles, where the red line segment represents the location of the actuator disk.
(a) Grid topology
(b) A close-up view of the computational mesh
(a) Grid topology
(b) A close-up view of the computational mesh
(a) Grid topology
(b) A close-up view of the computational mesh
2.2. Predicting Dynamic Mesh
In this prediction step, the Lagrange interpolation method is adopted to predict the dynamic mesh at every dynamic position, which is formulated as follows: where represent the predicted grid coordinates and are the Lagrange interpolation coefficients. For a certain tilt angle , these coefficients can be computed from
2.3. Correcting Dynamic Mesh
Because the interpolated actuator disk position does not conform to the actual instantaneous position, a correction step is needed. Here, the deformation from the predicted geometric configuration to the actual disk position is small. Therefore, a rapid elastic deforming technique can be used to correct dynamic mesh.
The correction values at the actuator disk surface are given as follows: where represent the grid coordinates at the actual instantaneous disk position, which are exactly calculated using the given rule of the considered motion.
Then, the correction formulas of the whole mesh are written as follows: where is defined as a function of the grid node indices , , and as follows: where is 1 in the actuator disk boundary and 0 in the corresponding far field boundary.
Table 1 lists the maximum tilt angles of IDW, RBF, and present dynamic mesh method. Results demonstrate that the predictor-corrector-based method is powerful enough to create valid dynamic mesh from 0° to 90°, while IDW and RBF become invalid at 48° and 81°, respectively. Therefore, by using the present method, the deformation ability is significantly improved, which laying the foundation for unsteady aerodynamic force calculation in conversion mode. Figure 10 shows the dynamic meshes at 30° and 60° by the present method.
(a)
(b)
3. Computational Method
3.1. Dual Time-Stepping Approach
In the Cartesian coordinate system, the three-dimensional unsteady N-S equations in a conservative differential form can be written as follows [22]: where is the vector of flow variables; , , and are the convective flux vectors; and , , and are the viscous flux vectors. An optional source vector, , can be added to the governing equations to model other forces acting on the fluid volume. After being discretized in space with a finite volume scheme, the time dependent N-S equations can be written as the following semidiscrete form: where represents the volume, the residual, and the index denotes the particular control volume.
The dual time-stepping approach is employed to solve the unsteady N-S equation. The time derivative in Eq. (7) is approximated by an implicit backward difference formula of second-order accuracy in the following form: where is the global physical time step and superscript denotes the time level. A time-stepping methodology is used to solve Eq. (8), which can be written as follows: where is the approximation to , denotes a pseudotime variable, and is the unsteady residual. Theoretical deduction indicates that the steady solution of Eq. (9) corresponds to the flow variables at the new time level, i.e., .
In case of dynamic grids, besides the conservation of mass, momentum, and energy, the geometric conservation law (GCL) must be satisfied in order to avoid errors induced by the deformation of control volumes. The integral form of the GCL reads where and are the grid velocity vector and outward facing unit normal vector of the surface , respectively.
3.2. Actuator Disk Modeling
The actuator disk approach ignores the effect of individual blades, but captures the overall effect of the rotor flow field. In the literature, actuator disk implementations fall into two general categories: a boundary condition approach and a source term approach. A good survey of these two approaches is presented by Le Chuiton [23]. The primary difference between the two approaches is how the rotor is treated in the fluid control volume. In the source term approach, the fluxes are computed just as they were for a typical internal cell. To make the presence of the actuator disk known to the flow solver, an extra source term is added to the momentum and energy equations.
The source term reads where and are the Cartesian velocity and body force components in the directions , respectively.
In the present work, a uniform actuator disk is used to perform the different test cases. The uniform pressure jump, , can be calculated as follows: where , , and denote the total thrust, angular velocity, and radius of the rotor, respectively. is the thrust coefficient, and is the density.
3.3. Turbulence Modeling
For closure of the N-S equations, it is necessary to calculate the turbulent viscosity in addition to the conservation variables. The Spalart-Allmaras (S-A) one-equation turbulence model solves a single partial differential equation for a working variable related to the turbulence viscosity . The S-A model takes the following form [24]: where stands for the local strain rate and is the distance to the closest wall.
3.4. Parallelization
The N-S solver is parallelized for more efficient calculations. Based on the basic block of a multiblock system, the flow in each block is calculated on a single processor. The data exchange between adjacent blocks is done on the block interfaces using functions of MPICH2. The multiblock decomposition of the structured grid is used for the load balancing. Meanwhile, the decomposition should make the data exchanges as few as possible [25]. For each block, the interface mapping relationships with adjacent blocks and the orders of data sending and data receiving should be established prior to the iterative calculations, and then, they can directly be used to accomplish the data exchange during the iterative calculations.
4. Results and Discussions
4.1. Validation
The flow field around the tilt-rotor aircraft is highly complicated, and before attempting to simulate the transition aerodynamic force, it is necessary to apply the current methodology to simpler problems to gain confidence in the solution algorithm. The first validation case is the A821201 airfoil used in the V-22 aircraft, which has a thickness ratio of 23%. Comparison of the computed chordwise pressure distribution with the experimental measurements [26] is given in Figure 11. Very good agreement between the calculated results and the measurements is observed.
The second application is ROBIN airframe with rotor configuration. Figure 12 shows the fuselage and actuator disk plane mesh. The fuselage has the length of , and the rotor radius is , where represents the characteristic length. The advance ratio for this case is , which corresponds to a free stream Mach number of 0.087. The fuselage angle of attack is 2.86°. The rotor thrust coefficient is 0.005. The sectional pressure distributions at -locations are shown in Figure 13. Overall, the present numerical results have a good agreement with the experimental data [27].
(a)
(b)
(c)
4.2. Tilt-Rotor Aircraft in a Conversion Mode
The corridor is the most critical phase of the tilt-rotor flight, and it is very important to understand the aerodynamic interactions between the fuselage, wing, and rotor in order to improve the aircraft performance and stability. The basic research strategies of this section are from 2D to 3D, from a quasi-steady approach to time-accurate unsteady simulation. Meanwhile, the influence of the different tilt modes and the period time on the unsteady aerodynamic force are also discussed.
In this work, the rotor blades were represented by means of a uniform actuator disk. The rotor thrust coefficient is made as an input, and the distribution is shown in Figure 14. The thrust coefficient decreases gradually and smoothly with the increase of the tilt angle. The operating condition in the section is as follows: the free stream Mach number is 0.2, the Reynolds number per chord is , the fuselage angle of attack with respect to the free stream is 0°, and the rotor tip Mach number is 0.72.
4.2.1. Results in Two Dimensions
Firstly, quasi-steady computations have been carried out in helicopter mode. The wingtip section is selected as the 2D case, and the thrust coefficient is 0.03. For comparison purpose, Figure 15 shows the pressure contours and streamline distributions with and without actuator disk. As can be seen, there is an obviously pressure jump between the upper and lower actuator surfaces. The source terms located on the disc bottom side impart impulse and energy to the fluid, and the flow is deflected and accelerated through the disk. As shown in Figure 16, a clear difference in is noticeable for computations with and without actuator disk. Downwash induced by rotor decreases the actual angle of airfoil in helicopter mode.
(a) Without actuator disk
(b) With actuator disk,
In airplane mode, if there is no source term effect, the flow characteristics are the same as Figure 15(a). However, when the source terms act on the actuator disk, the flow presents a contraction due to the flow acceleration, as shown in Figure 17. The contraction ratio is determined by . Figure 18 shows a comparison of the difference in in the presence of the actuator disk.
Figure 19 shows the lift and drag coefficients in a fixed conversion mode. These results are from steady-state computations of the flow. Due to the reduction of the rotor thrust and the change of the source term action direction, with the increase of tilt angle, the lift coefficient first increases, then decreases, and increases again. However, the drag coefficient increases constantly. The critical value of lift coefficient is located at 29°, which draws the same conclusion as in Ref. [28]. At , the difference of with and without actuator disk is greater than 40%. Figure 20 shows the computed Mach contours at 29° and 40° tilt angles; as can be seen, due to the flow acceleration at the lower surface, lift coefficient curve shows a drop.
(a)
(b)
A time-accurate calculation on the unsteady aerodynamic force of the continuous conversion process is discussed below. Two tilt motion rules are employed in this work; there are constant angular velocity rotation and sinusoidal angular velocity rotation, respectively. Here, we assume the entire tilt period is 12 seconds. Figure 21 shows the tilt angle varying with time of these two conversion modes.
The computed unsteady lift and drag coefficients are compared with the quasi-steady data in Figure 22. Overall, the unsteady lift coefficient curves have the same variation trend with the quasi-steady data. However, the unsteady drag curves occur with difference. For different conversion modes, negligible differences are found before the critical tilt angle (29°). After that, the unsteady lift coefficient is greater than the quasi-steady value, while the unsteady drag coefficient is less than quasi-steady data. At the end of the conversion process, there is a step between constant angular velocity rotation and quasi-steady computed values. This is because the angular rate is not zero. Table 2 gives the quantitative comparisons of the aerodynamic coefficients at the tilt angle of 45°. The maximum lift difference is up to 7.36%.
The influence of the tilt period time on the unsteady aerodynamic force is also discussed. For this case, we adopt three different period times for transition flight, namely, 6 seconds, 12 seconds, and 24 seconds. The tilt motion rule is assumed to be constant angular velocity. Figure 23 shows the computed unsteady aerodynamic coefficients. As can be seen, before the critical tilt angle, differences among the conversion period time are very small. After the tilt angle of 30°, for , the lift coefficient is the largest and the drag coefficient is the smallest. Due to the unsteady effects, the step at the end of transition is also more obvious.
4.2.2. Results in Three Dimensions
Figure 24 shows the computed quasi-steady lift and drag coefficients at fixed tilt angles. Compared to 2D case, the 3D lift coefficient curve has the same variation tendency. Both of them first increase, then decrease, and increase again. The critical tilt angle is also almost the same, but the maximum lift coefficient of 3D is at the tilt angle of 90° which is different from 2D case. The lift coefficient value with actuator disk is greater than 12.2% of the value without actuator disk at 90°. The drag coefficient curve is totally different from 2D case due to the three-dimensional effects. It first decreases, then increases, and finally decreases quickly with the increase of the tilt angle. In other words, the interaction between rotor-wing-fuselage is more complicated than 2D condition, which relates to many factors.
Three spanwise locations are selected for analyzing the chordwise pressure variation. These locations are situated at 0.3R, 0.6R, and 0.9R from the rotor center and noted as Position 1, Position 2, and Position 3, respectively, as shown in Figure 25.
As we know, downwash induced by a rotor decreases the actual angle of the wing sections. Pressure distributions on the wing locations at different tilt angles are presented in Figure 26; the symbol on the horizontal axis represents the nondimensional airfoil chord length. It can be seen from the figure that at Position 1, the lowest upper surface pressure is at 30°, while at Position 2 and Position 3, the lowest upper surface pressure is at the tilt angle of 90°. At all positions, the upper surface pressure in helicopter mode is the largest. Meanwhile, at different tilt angles, the more close to the inside of the wing, the more lift will be generated at the spanwise stations.
(a) Position 1, 0.3R
(b) Position 2, 0.6R
(c) Position 3, 0.9R
Figure 27 shows the 3D time-accurate lift and drag coefficients for different tilt modes. It can be seen that the tilt modes have little influence on the unsteady aerodynamic coefficients in this 3D case. Table 3 gives the quantitative comparisons of the aerodynamic coefficients at the tilt angle of 45°. The maximum lift difference is 1.68%, and the unsteady drag coefficient has decreased by 1.19%. Figure 28 gives the 3D time-accurate aerodynamic coefficients for different conversion period times. It demonstrates that the unsteady effect is larger when less time is used for transition.
5. Conclusions
The developed CFD solver is able to capture essential flow characteristics and predict the unsteady aerodynamic force of the tilt-rotor aircraft in conversion mode. The effects of different tilt modes and the tilt period time are also conducted. These investigation results could provide a good foundation for tilt-rotor aircraft design in the future. The conclusions are as follows: (1)A novel predictor-corrector-based dynamic mesh method for the multiblock structured grid with extremely large deformation is proposed. Compared with the IDW and RBF methods, the present approach shows stronger deformation ability(2)Both 2D and 3D, quasi-steady and time-accurate are applied to the aerodynamic force computations of the tilt-rotor aircraft in conversion mode. With the increase of tilt angle, the lift curves first increase, then decrease, and increase again. However, the drag coefficient curves occur with difference between 2D and 3D(3)For constant angular velocity rotation and sinusoidal angular velocity rotation, the difference between 3D quasi-steady and time-accurate computations is comparatively small. Finally, the less time for transition, the more unsteady effect is observed
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Priority Academic Program Development of Jiangsu Higher Education Institutions.