Abstract
This paper presents the development of a tool integrated in the UNS3D code, proprietary of Alenia Aermacchi, for the simulation of external aerodynamic flow in a rotating reference frame, with the main objective of predicting propeller-aircraft integration effects. The equations in a rotating frame of reference have been formulated in terms of the absolute velocity components; in this way, the artificial dissipation needed for convergence is lessened, as the Coriolis source term is only introduced in the momentum equation. An Explicit Algebraic Reynolds Stress turbulence model is used. The first assessment of effectiveness of this method is made computing stability derivatives of a NACA 0012 airfoil. Finally, steady Navier-Stokes and Euler simulations of a four-blade single-rotating propeller are presented, demonstrating the efficiency of the chosen approach in terms of computational cost.
1. Introduction
The main objective of the present work is to calculate the flow field around a propeller. Unsteady Reynolds averaged Navier-Stokes equations (RANS) represent the state of the art for the numerical prediction of the viscous flow around propellers [1, 2]. However, unsteady Navier-Stokes simulations for such complex flows are computationally expensive. Therefore, these simulations are generally carried out only using reduced order methodologies [3]. In our case, to simulate the unsteady viscous flow around a propeller, we use the Navier-Stokes, or Euler, equations in a noninertial (rotating) reference frame in which the propeller is at rest. In this way, the simulations can be performed in steady-state mode. The advantage of this method is the relatively lower computational cost needed to obtain a CFD solution in comparison to a fully unsteady simulation. Moreover, the present method permits the insertion of the propeller as a building block into the full computational model of an aircraft, to investigate propeller-aircraft integration issues from an aerodynamics point of view.
When the Navier-Stokes equations for a steady viscous compressible flow are written in a rotating frame of reference, there are two choices regarding the velocity vector components. Either they can be the components with respect to the absolute (inertial) frame, hereafter called the absolute velocity components [4], or they can be the components with respect to the rotating (noninertial) frame, hereafter called the relative velocity components [5]. In this paper, the first formulation has been chosen because it has the advantage of enabling a steady-state formulation as the flow field can be viewed as a steady state in the rotating frame. In order to fully use the UNS3D code which is written in the absolute velocity components, source terms are added to take the coordinate-system rotation into account in the governing equations.
To assess this method, two applications are presented: the first application is the determination of the stability derivatives for a NACA 0012 airfoil and the second is the simulation of the flow field around a rotating propeller. The first application allows us to assess the method by comparing the results to a reference test case, for which several authors have obtained a solution [6–9]. In this case, the flow field is computed using Euler equations. The second application is a more complex test case, which has the purpose to evaluate the accuracy, efficiency, and robustness of the current method to predict the complex flow field of a rotating propeller. For this test case, a comparison with experimental results for cruise conditions in terms of thrust coefficient is also made.
2. Governing Equations
The governing equations which are numerically resolved using the UNS3D code are summarized below. The Navier-Stokes equations in a noninertial frame of reference are expressed as where is defined as and are the respective flux vectors: and the source terms are contained within : is equal to the following relation:
2.1. Formulation in Rotating Frame for the Absolute Velocity
Normally the governing equations for rotating systems are solved for the relative velocity components in the relative frame of reference. However, in an external aerodynamic context such as the propeller-wing configuration, two main advantages can be derived through taking instead the absolute velocity components in the relative reference frame. The first, and most important one, is that the grid block containing the propeller can be interfaced more easily within the full aircraft model. The second is that a smaller numerical effect of source terms is present, as will be shown below. To express (1) in terms of a relative reference frame, the following relations for substantial and local derivatives are used: where the prime ′ denotes the operation with respect to the relative reference frame. By using relations (7), the right-hand side of (1) becomes where is defined as
With this formulation, the source term vector (see (8)) contains only the contribution of the Coriolis force and the contribution of the centrifugal force is omitted. In this way the magnitude of the source term is greatly reduced and a smaller amount of artificial dissipation is required to ensure convergence.
To take into account the rotation of the coordinate system, (3) is modified in the following way: where is defined as
It is noteworthy to mention that the above expression for the Navier-Stokes equations is valid only for a rotating reference frame. Linear acceleration components should be added in case one wants to express the equations of motion in a general noninertial frame.
In addition, it should be pointed out that it is not possible in the present context to impose any radial equilibrium condition to the pressure, as it is normally done in turbomachinery.
An important thing to be noted here is that, except for the source term , the functional form of the noninertial Navier-Stokes equations is the same as the functional form of the standard conservative equations defined for inertial reference frames and including the Algebraic Lagrangian Eulerian (ALE) approach for generalized motion of the grid. Therefore, it is possible to implement a conservative formulation in terms of the conservative variables defined in (2) and the introduction of the ALE approach permits a local application of the noninertial frame of reference as a building block in a more complex configuration framework, without any interface between the noninertial and inertial part of the same mesh, because this formulation guarantees the flux conservation. In applications, it is sufficient to specify as input the bounds of the rotational region and in this way a rotational speed is imposed at the nodes inside this region.
2.2. Enforcement of Compatibility Conditions
The numerical source error due to the noninertial reference frame can be examined analytically by imposing the conservation of the freestream. In this case, all the flow derivatives are zero and the velocity vector is where is equal to
The continuity, momentum, and energy equations (see (1)) can then be reduced to the following expressions:
For the first application case, where we have a steady rotation parallel to the -direction, (14) can be reduced to the following expression: that is identically zero for any nonzero angular velocity , whereas, for the second case, with a rotation parallel to -direction, (14) becomes
In both cases, for the numerical formulation, the right-hand side is not exactly zero, however, producing a freestream error.
Using the results of (17) and (18), and denoting the right-hand side as , a simple and straightforward source term correction can be applied in (1). In particular, an additional source term can be included to exactly cancel the freestream error: where for the NACA 0012 is whereas for the propeller (20) becomes
3. Numerical Method
The computations have been performed using the code UNS3D. The solution algorithm is based on a finite volume, node centred approach operating on a hybrid unstructured grid. The artificial dissipation model is derived from the nonlinear scheme of Jameson [10], with no eigenvalue blending. Scalar or matrix dissipation can be chosen.
The Navier-Stokes equations are integrated in time with a second order backward difference scheme and dual time stepping. A five-stage Runge-Kutta scheme is used to drive toward zero the residual at each time step. With the use of residual averaging, a local CFL number of 4.9 could be employed in the multistage subiteration process. The Algebraic Lagrangian Eulerian approach for generalized motion of the grid is included [11].
The Weiss and Smith version of low Mach number preconditioning is implemented in the code [12]. A sensor depending on cell Reynolds number was also introduced to avoid applying the preconditioning inside boundary layers. For the computation of the present test case, its application was found to be beneficial in order to reduce numerical dissipation and enhance convergence.
Matrix dissipation was also found to be beneficial, allowing a strong reduction of the dissipation associated with convective eigenvalues, hence enabling a better resolution of vortices.
3.1. RANS Turbulence Model
The k-ω turbulence model proposed by Hellsten [13] has been employed. The model constants have been calibrated requiring consistent behavior near boundaries between turbulent and laminar flow, inside shear flows and for zero pressure gradient wall flows. In particular, the calibration has been considered taking into account a variable , as it is the case if an algebraic stress model (EARSM) is included.
The Wallin-Johansson Explicit Algebraic Stress Model (WJ-EARSM) [14] is implemented using Hellsten’s k-ω as the basic RANS model. The model is an exact solution of the corresponding ARSM in two-dimensional mean flow. In three dimensions there is still a complete, while approximate, solution.
The full anisotropic version of the model is used; that is, the anisotropic part of the Reynolds stress tensor is directly introduced in the momentum equations, while the isotropic part is taken into account in the form of an effective variable .
4. Boundary Conditions
The boundary conditions along solid walls for Navier-Stokes (viscous) flows are different from those for Euler flows. In the case of viscous flows, the velocity of the flow must vanish at the walls, while, in the case of Euler flows, it is only required that the flow does not go through the wall.
As a consequence of the foregoing statement, at the airfoil, the condition of nonpenetration has been imposed, whereas, on the blade surface, no-slip and no-penetration conditions are used by setting the absolute velocity equal to the absolute local blade velocity and the adiabatic wall condition and zero-normal pressure gradient condition at the wall are imposed at the blade surface.
In general, the boundary conditions applied at the far-field boundary are the same for Navier-Stokes and Euler flows; therefore, the far-field boundaries are treated by using characteristic boundary conditions.
5. Numerical Results
5.1. Model Validation: Steady Rotary NACA 0012 Airfoil
To validate the numerical model, the stability derivatives for a NACA 0012 were computed using finite differences and compared with the results obtained by Limache and Cliff [6]. In the experiment, an airfoil is submitted to a steady rotation performed at constant incidence for a given pitch rate , generating a steady flow field in a reference frame attached to the airfoil. The radius of the loop is inversely proportional to . Thus, as reduces to zero, the radius approaches infinity and steady level flight is recovered.
The results presented below are all computed for an angle of attack equal to zero, so we use the wind-axis reference frame for the computation of the derivatives.
5.2. Numerical Results
In Figure 1 is shown the 2D unstructured grid. The outer boundary is at a distance 30 times the length of the airfoil’s chord with respect to the grid center, coincident with the leading edge of the airfoil. The grid is made by 12334 nodes and 12096 elements.
(a)
(b)
To verify the implementation of 3D Navier-Stokes equations in terms of absolute velocities, we compare results for the NACA 0012 airfoil rotating at a finite to those produced by Limache [7] simulating inviscid flow around a NACA 0012 airfoil at Mach equal to 0.2 for nondimensional pitch rate equal to 0, 0.01, 0.03, and 0.05. In fact, at the present test conditions (low Mach number and low incidence) we do not expect that the integral quantities computed using viscous and inviscid methods, respectively, differ significantly.
The distributions around the airfoil are shown in Figure 2 where the computed distributions and streamlines of relative velocity are compared with those presented in [7] for = 0.01. In Figure 3, the pressure coefficient contours and streamlines of relative velocity in the whole computational domain are shown.
(a)
(b)
(a)
(b)
Finally, in Table 1, we compare the values of and from our implementation to the references results. The two implementations match quite well over a range of values.
The stability derivatives are calculated using finite differences:
In Table 2, the stability derivatives are compared with the results obtained by Limache.
It is possible to conclude that the results obtained by UNS3D are in good agreement with the numerical results obtained by Limache and Cliff [6].
5.3. Geometrical Model
In Figure 4, it is possible to see the geometry of the experimental model used by Biermann and Haetman [15]. The experimental results were performed for four- and six-blade single-rotating and dual-rotating propellers with and without the symmetrical wing in place. The maximum propeller speed was 550 rpm. The results for four-blade single-rotating propeller were made up with two two-way hubs mounted in tandem and the spacing between front and rear blades is not equal and therefore the front blade led the rear by 85.4 deg. In this paper, only the four-blade single-rotating propellers with and without wing are considered for the comparison with the experimental results in terms of thrust coefficient. The propeller, namely, a Hamilton Standard 3155-6, consists of four blades installed on a single way hub in front of a streamline body, or nacelle, housing the engine needed to spin it. The four blades are streamlined using Clark Y profiles and the angle between two blades is 90 deg.
In the report of Biermann and Haetman [15], several blade pitch angles, defined as the angle between the rotation plane and the airfoil chord at 75% of the radius of the propeller, in the range of 20 deg to 65 deg were studied. In our case, it was decided to investigate a propeller with a blade pitch angle of 45 deg. The propeller diameter is 3.08 m. Starting from the geometrical details reported in the mentioned report, a mathematical model describing the propeller has been created with CATIA V5 (Figure 5). The wing, shaped using NACA 0012 airfoils, is located in a midposition of the spinner and set at an angle of attack of 0 deg. Wing chord is 1828.8 mm and wing span is 4241.8 mm.
(a)
(b)
5.4. Results for the Four-Blade Single-Rotating Propeller + Spinner (Viscous)
A 3D unstructured grid of the propeller + spinner has been generated with ICEM-CFD (Figure 6). This grid is strongly refined in the region around the blades and on the blade surfaces. It is made by(i)5672824 nodes,(ii)16265544 elements,(iii)21 prismatic layers on solid surfaces to correctly match boundary layer behavior.
In the code UNS3D it is possible to specify an arbitrary velocity for a specific group of nodes within the mesh (ALE formulation). The resulting fluxes are automatically interfaced in order to ensure conservation at the boundary between rotating and nonrotating zones. The solution is then computed specifying a rotational velocity, as described in (11), only for the nodes inside the rotating block (Figure 7) and taking into account the source terms (see (8)) and the correction terms as in (23). It is worth noting that the formulation described above is steady. The solution is indeed an instantaneous snapshot of the flow field, which can be seen as frozen in correspondence with a fixed phase angle.
(a)
(b)
Five different operating conditions, shown in Table 3, were investigated. The axial undisturbed velocity has been set equal to 49.1744 m/s, corresponding to the maximum wind tunnel test speed of 110 mph [15].
In Figure 8, it is possible to note a general increasing level of the relative Mach number from the nacelle surface toward the tip, which is the result of the increasing rotational speed with increasing radius.
(a)
(b)
As the propeller rotates, it induces swirls in the slipstream and the blade tip vortices pass by periodically. This phenomenon is more evident when the advance ratio decreases. In fact, for lower advance ratios, we see a strong vortex shedding, which starts from each blade and travels downstream with the perturbation velocity creating strong spiral type regions in the rear wake for each blade (Figure 9). Furthermore, strong hub and tip vortices (Figure 10) are continuously shed from the respective blade regions and “absorb” the weaker vorticity regions at inner blade radii producing also spiral type patterns.
(a)
(b)
In Figures 11 and 12, the azimuthal velocity profiles downstream of the blades are shown at different position of and it is possible to note how the swirl induced by the rotation of the propeller vanishes with increasing distance from the blades. In particular, the fast decay of the slipstream swirl, in the region far from the spinner, is due to the coarse mesh and consequently to the high levels of numerical dissipation in that region.
With the aim of comparing the results with those obtained by Biermann and Haetman [15], the thrust coefficient, defined as has been calculated. Following the experimental procedure adopted by Biermann and Haetman [15], the thrust force has been obtained by integrating the forces along -direction on all the blade surfaces and subtracting the drag force due to the blades alone in case of zero thrust coefficient.
In Figure 13, the obtained thrust coefficients for five different advance ratios are plotted and compared with those obtained by Biermann and Haetman [15].
The computed thrust coefficients are in good agreement with the experimental values as listed in Table 4.
As the experimental errors are unknown, it is not possible to determine whether the computed results are within the range of the experimental uncertainty or to give a precise assessment of the quality of the results.
5.5. Results for the Four-Blade Single-Rotating Propeller + Spinner + Wing (Inviscid)
Steady Euler results for the propeller + spinner + wing are presented in this section. The grid is generated with ICEM-CFD (Figure 14) and it is made by(i)2173935 nodes,(ii)10718702 elements.
(a)
(b)
A rotational velocity has been imposed for the nodes inside the block around the blades, as indicated in Figure 14. The investigated operating conditions are identical to the previous case.
The Mach distributions in Figure 15 clearly show the effect of the propeller slipstream that washes the wing. In particular, the impact of the swirl velocity component is very pronounced. It should be noted that the Mach number distributions, much like the pressure distributions, as presented in Figure 16, are affected by both the local propeller induced flow angles and the dynamic pressure increases in the slipstream.
(a)
(b)
(c)
Another phenomenon that is clearly visible, due to the interference between the propeller and the wing, is the rise of vortices around the juncture of spinner and wing (Figure 17). A shedding of these vortices can be evident, which is indicative of the high gradient of spanwise load on the wing.
In fact, the streamwise and spanwise locations of blade vortices are staggered on the upper and lower surfaces of wing and interacted vortices are induced near spinner.
In order to make a further verification of the method the numerical calculations have been compared to experimental data (Figure 18).
Again, the computational results are in good agreement with the experiment and the maximum error is around 5% (Table 5). Also in this case the experimental errors are unknown.
6. Conclusions
The Alenia Aermacchi UNS3D code was modified, introducing the capability of flow simulations in a noninertial frame of reference. The modified code was at first applied to the computation of damping derivatives of a rotating airfoil and then to the prediction of the performance of a propeller, following the experimental test case described by Biermann and Haetman [15], for different rotational speeds. In the first case, good agreement has been obtained with the numerical results of [7]. In the second, the results are in good agreement with the experimental data within the propeller operating range. The computational results clearly showed the effect of the swirl velocity and the increased total pressure on the spinner and the wing. Therefore, this approach facilitates the identification of typical flow phenomena, like the deformation of the slipstream when passing the wing, being able to model aerodynamic phenomena linked to the propeller-airframe integration.
Nomenclature
| : | Generic scalar |
| : | Generic vector |
| : | Airfoil chord |
| : | Generic aerodynamic coefficient |
| : | Lift coefficient |
| : | Pitching moment coefficient |
| : | Pressure coefficient |
| : | Generic stability derivative |
| : | Thrust coefficient |
| : | Propeller diameter, [m] |
| : | Total energy per unit of mass, [J/kg] |
| : | Vector of external forces, [N] |
| : | Inviscid flux vector |
| : | Viscous flux vector |
| : | Advance ratio |
| : | Mach number |
| : | Propeller rotational speed, [rps] |
| : | Pressure, [Pa] |
| : | Pitch rate |
| : | Normalized pitch rate () |
| : | Heat flux, [W/m2] |
| : | Source vector |
| : | Position vector relative to the rotation center, [m] |
| : | Time in physical space, [s] |
| : | Effective thrust, [N] |
| , , : | Freestream velocity components, [m/s] |
| , : | Velocity vector, [m/s] |
| : | Rotational speed vector of the coordinate system, [m/s] |
| : | Generic component of the velocity vector, [m/s] |
| : | Freestream velocity vector, [m/s] |
| , , : | Components of velocity vector, [m/s] |
| : | Work of the external forces, [J] |
| : | Vector of conservative variables |
| , , : | Position vector components relative to the rotation center, [m] |
| : | Angle of attack |
| : | Kronecker symbol |
| : | Density, [kg/m3] |
| : | Generic component of the shear tensor, [Pa] |
| : | Rotational speed vector, [rad/sec] |
| , , : | Components of rotational speed vector, [rad/s]. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The present work has been performed through a close cooperation with Alenia Aermacchi S.p.A. and through the use of Confidential Information and Data, property of Alenia Aermacchi S.p.A. which remains the sole owner of all such relevant IP rights. The results of the present work shall be, therefore, property of Alenia Aermacchi S.p.A.