Parametric Investigation of a Distributed Propulsion System on a Regional Aircraft

: With the mandatory requirement for more efﬁcient aircraft due to both economic and envi-ronmental purposes, academy and industry are exploring new aircraft design opportunities including the concepts of hybrid-electric and fully electric vehicles. Within this framework, distributed electric propulsion is a key technology for future aviation, as it allows the installation of a theoretically indeﬁnite number of small motors. The blowing effect induced by the propeller can be used to improve aerodynamic performance, hence, thanks to their reduced size, these small motors could be installed along the whole span covering the whole wing. This paper presents a study devoted to the investigation of the aerodynamic effects of distributed electric propulsion installation on a regional aircraft, computing the aerodynamic coefﬁcients using high-ﬁdelity CFD simulations via the RANS approach. Different propeller diameters and trust levels were analysed in climb and landing conditions, applying periodic boundary conditions on a ﬁnite span section of the wing, simulating an inﬁnite rectangular wing. The goal of the current study is to quantify the increase of aerodynamic coefﬁcients with reference to power off condition and report data as a function of the propeller’s characteristics. The objective is to identify and propose a simpliﬁed analytical formulation to be used in the phase of preliminary design. The implementation of such a formula in lower-ﬁdelity tools will allow fast and reliable procedures for preliminary conceptual design.


Introduction
Aviation has grown significantly in the last decades, and a rapid increase in total pollution attributable to air transport is evident, despite the development of more fuelefficient turbofan or turboprop engines. The reduction of greenhouse gas emissions is a key point in Flightpath 2050 vision [1] and represents a major goal in the framework of climate preservation. To this aim, regulations on the reduction of pollutants such as CO 2 and NOx and noise emissions are becoming more and more rigid. In light of this new development, air transport has undertaken a redesign of existing aircraft and the proposal of innovative architectures and technologies. One of these new technologies, which is studied in this paper, is Distributed Electric Propulsion (DEP). The concept behind DEP technology is to place several small electric motors in many different locations on the aircraft, thanks to the smaller size and weight compared to conventional engines. Conceptual designs have been proposed and studied, such as the Blended Wing Body (BWB) vehicle N3-X [2,3] or the newly proposed Airbus "ZEROe" [4]. Studies on NASA's N3-X regarding energy consumption [5], control system for propulsion [6], and noise and emissions [7] have also been published. Furthermore, ONERA started a project dedicated to DEP investigation [8], while the benefits of hybridization within the regional aircraft scale are reported in [9], as well as drag reduction achieved using DEP in [10]. A review of distributed electric concepts is also described in [11]. Two possible configurations are shown in Figure 1. The NASA SCEPTOR is depicted on the left, with propellers mounted in front of the wing leading edge, while the ONERA AMPERE on the right side shows a cluster of ducted fans mounted on the top of the wing. For the SCEPTOR program, the usage of distributed electric propulsion has been investigated by Dubois et al. [12] and Deere et al. [13], while Patterson et al. [14,15] formulated an improved method for the estimation of the lift augmentation induced by DEP. The leading edge propellers are viewed as high-lift devices rather than propulsive devices as they are operated only at low-speed flight conditions to increase the speed of the airflow over the wing. The advantage is that the wing chord can be reduced, increasing the aspect ratio with consequent benefits on the efficiency of the aircraft. Contributions to DEP are not limited to aerodynamics but range from structural aspects [16] to the proper design of batteries [17], as well as involving issues in control [18].
Preliminary design methods rely on empirical formulations that have been extensively validated and verified by experimental and numerical data [19][20][21]. The availability of such tools allows a fast and reliable evaluation of the aerodynamic performance of an aircraft, nevertheless, most of these tools cannot be applied for configurations with distributed propulsion. The codes that allow to simulate the presence of the propeller are based on the induced increment of velocity. The scope of this work is to correlate the effect on the aerodynamic coefficients to the propeller characteristics, i.e., thrust and advance ratio, identifying the key parameters that can affect the performances of a configuration in presence of distributed propulsion and, eventually, provide a draft formulation. In order to achieve this goal, a study of the effects of DEP on lift and drag coefficients is carried out on a finite wing span section. The idea is to simulate an ideally infinite span distribution of propellers with different diameters, applying a periodic boundary condition. The aim is to identify tendency lines using the results of RANS simulations, varying the diameter of the propeller and the thrust, in order to provide a simplified formula to be used by preliminary design tools. Particularly, in Section 2 the tools used to perform the computational analysis, namely grid generator and CFD flow solver, are described. In addition, a matrix providing the complete set of test performed is given. The main results regarding the two flap configurations studied in this paper are provided in Section 3. Finally, the most important conclusions extracted from this study are pointed out in Section 4. This article represents initial research in this area, aiming at the development of design tools capable of predicting aircraft performance in presence of DEP.

Computational Analysis
In this section, the computational tools used for the grid generation and the CFD evaluation are explained. In order to reduce the effort required for the CAD and meshing, an automatic procedure to generate the geometry and the computational mesh was established. The first step is performing a linear extrusion of the airfoil, without introducing twist or sweep angles, generating a finite span wing. In the second step, the propeller, modelled as an actuator disk, is positioned in the centre of the wing span with a prescribed diameter and offset with reference to the leading edge. The final step is the grid generation using an automatic ad-hoc procedure. The presence of infinite propellers is simulated applying periodic boundary conditions at both sides of the domain. The complete test matrix studied is also reported in this section. Moreover, it must be mentioned that two experimental campaigns are already scheduled to assess the aerodynamic and aeroacoustic effects inside CICLOP (Characterization of the Interaction between wing and Closely Operating Propeller, H2020 CS2-CFP11-2020-01 call, Grant agreement No: 101007567) and VENUS (inVestigation of distributEd propulsion Noise and its mitigation through wind tUnnel experiments and numerical Simulations, H2020 CS2-CFP10-2019-01 call, Grant agreement No 886019) European projects, respectively. Thus, these experimental data will be, then, used to validate this numerical study. In particular, the objective of CICLOP, led by Technische Universitat Braunschweig, is to close a gap in current predictions of aerodynamic effects of Distributed Electrical Propulsion and close wing coupling at high lift based on aerodynamic models and simulation, by providing high-fidelity experimental data set that will allow validation of theoretical tools and understanding of DEP aerodynamics. Additionally, VENUS, led by Università degli Studi Roma Tre, aims to understand the physics behind the aeroacoustics of DEP through deep theoretical, experimental, and numerical study. Appropriate numerical procedures for DEP noise assessment will be set up and an experimental dataset obtained in dedicated wind tunnel tests will be used both as experimental DEP noise validation reference and for providing support to the identification of the main parameters affecting DEP noise.

Grid Generation
To perform the CFD analysis, a parametric multiblock structured body-conformal hexahedral grid was generated using PARDOMO (parametric domain modeller) and the in-house-developed grid generator ENGRID. The domain is made of 69 blocks for three grid levels considered, and there are, approximately, 8 × The use of a parametric grid generation allows considerable time savings on the whole activity as, once the procedure is written, the user has to modify few parameters in the configuration file in order to obtain a new computational domain with different actuator disc size and position, as well as different flap settings. In Figure 2, the grid blocks near the airfoil, together with the airfoil surface, and the actuator disk are depicted.

CFD Evaluation
The numerical steady-state simulations were carried out using the in-house-developed code ZEN (Zonal Euler Navier-Stokes flow solver). ZEN is a multiblock structured flow solver for steady and unsteady RANS equations, which has been developed at CIRA for more than two decades [22][23][24]. It is based upon cell-centred, finite-volumes formulation, with central schemes. Convergence toward steady state is achieved by explicit multistage Runge-Kutta schemes, with acceleration techniques such as local time stepping, residual averaging, and multigrid [25]. Several turbulence models are available; all solutions in the present work were computed using the k − ω TNT two-equation model [26].
Moreover, actuator disk boundary conditions are available, both for steady and unsteady computations [27]. The term actuator disc refers to a simplified theory that allows simulating the propeller effects by applying an appropriate impulse to the fluid crossing the surface of a disc defined by the path of a propeller blade tip. This is equivalent to applying a source term to the transport equations of the flow variables [28]. Applying the actuator disc theory, it is possible to avoid the complexity of modelling the complete rotary wing. In case of propeller aircraft, the solution of Euler or RANS equations coupled with a technique based on the actuator disc theory allows to simulate reasonably well the effect of the propeller induced field on the wing. The propeller impulse, in general cases, depends upon the time and spatial position on the disc surface. In present simulations, it is considered that impulses change only with radial distance from the propeller axis. Specifically, one actuator disk per DEP was analysed. By applying periodic boundary conditions at the wing end planes, an infinite number of corotating actuator disks is simulated. Whereas, when symmetric boundary conditions are applied at the wind end planes, alternate counter-rotating propellers are simulated. In this study, periodic boundary conditions were considered.
Finally, to monitor and determine if a CFD solution has converged, a drop of three orders of magnitude of the flow solution residuals must be reached. Additionally, considering an interval of 10,000 iterations, it was required that the standard deviation of the drag coefficient be less than 0.001. Figure 3 shows an example of convergence history for the continuity equation (left) and aerodynamic coefficients (right).

Boundary Conditions
The main body and flap surfaces are treated as adiabatic no-slip wall, while the actuator disk boundary condition is used to simulate the propeller, as described in the previous section. On the farfield in the xand z-directions (see Figure 2), corresponding to a squared box of side approximately equal to 80 chords, the boundary condition based on the Riemann invariants is imposed. In the y-direction, the periodic boundary condition is applied on both sides of the domain.

Test Matrix
The goal of this work is to perform a parametric study to evaluate the effectiveness and performance of Distributed Electric Propulsion. However, to reduce the number of analyses to perform, it was decided to start from a realistic configuration. In particular, the reference configuration considered, in terms of free-stream values and wing sectional airfoil, is based on a typical regional aircraft with 40 passengers. The aerodynamic database was built considering four DEP configurations by varying the number of propellers along the wing.
The wing span (b) is set in order to include the propeller's diameter plus an additional 10% of the diameter, i.e., (b = d + 0.1d). As the centre of the propeller corresponds to the centre of the domain, the previous assumption means that the distance between the tip of two adjacent propellers is equal to 10% of the diameter itself. The diameter-to-chord ratio (d/c) is computed considering that the wing chord (c) is 2.57 m long. A summary of the DEP configurations is provided in Table 1.   A matrix for the different test cases is provided in Table 2. A total of 27 test cases were studied. Preliminary analyses (see Figure 5) showed that CTR, defined as CTR = T/(d 2 V 2 ρ) = C T /J 2 whereas C T = T ρn 2 D 4 , the Renard thrust coefficient, is a crucial parameter. As a consequence, it was decided to consider CTR as a variable for carrying out the parametric study. The second variable considered is the diameter-to-chord ratio d/c, as its effect on lift augmentation is already presented in [15].

Results
In this section, the main results obtained by CFD evaluations are presented divided in Flap A and Flap B configurations. The polar curves of the test cases with the highest d/c ratio are provided for both configurations. In addition, for each set of test (d/c = 0.802, 0.533, 0.401, 0.319), the increment in percentage of the aerodynamic coefficients is reported for three angles of attack 8 • , 10 • and 12 • . The increment is calculated with respect to the corresponding power-off condition. The increments are given with reference to the previously defined thrust coefficient (CTR). The motivation behind the definition of this parameter is shown in Figure 5. In particular, the parameter CTR is introduced to identify a tendency line for the increment of the aerodynamic coefficients, even if different values of J are used. Let us vary the advance ratio by halving the free-stream velocity, and keeping constant the number of propellers and the rotations per second. Thus, the working conditions of the propeller are changed, but the propeller characteristics are maintained. The reference case is the take-off condition with d/c = 0.802 and n equal to 38.79 rps at nominal free-stream velocity (V ∞ = 59.84 m/s). The thrust levels analysed are those in Table 2. The angle of attack is α = 8 • and the simulations are performed at standard sea level conditions. The same conditions are analysed with a different free-stream velocity (V ∞ = 29.92 m/s). An additional thrust level (T = 19,000 N) was also analysed for the case with free-stream velocity equal to 59.84 m/s. This thrust level was analysed in order to have a CTR point in common (CTR ≃ 1.0) for the test cases at both velocities. The current set of input values returns the same C T for both conditions, but different J and, as a consequence, different CTR. Figure 5 clearly shows the reason behind the choice of representing the increment of aerodynamic coefficients vs. CTR. In fact, the results tend to lay on a continuous curve. In this way, tendency lines are built based on the obtained results for lift and drag coefficients. For the lift coefficient, the authors propose a power trend, considering that the lift can increase with thrust up to a limit that is given by the performance of the airfoil and for CTR = 0, the condition of zero increment in lift is automatically satisfied. Differently, for the drag coefficient, a polynomial function of second order and zero intercept is selected since an increase in drag is always expected when the thrust rises, and at zero thrust, power-off conditions are found. It must be remarked that these tendency lines were selected taking into account only numerical data, hence, they have to be considered as a starting point to be improved with additional data and assessed when experiments become available.

Flap A Working Conditions
Firstly, the polar curves obtained for the configuration with d/c = 0.802 are provided in Figure 6. In particular, the polar curves for the test cases A1-G1 (see Table 2) are compared with the power-off conditions. Note that, data at angles of attack after stall conditions are not provided. In addition, only steady simulations were completed. Therefore, if at some angle of attack C l and C d are not given, this is due to the unsteady behaviour found in the CFD evaluation. A more in-depth analysis was carried out for some cases where unsteadiness was observed. At angles of attack prior to the stall condition, a timeaccurate convergence led to a steady-state solution. In the post-stall region, periodic and nonperiodic solutions were achieved. As an exhaustive analysis on the performance of the wing-propeller setup was beyond the scope of the current research, these cases were bypassed.
As expected, there is an increase in lift and drag coefficients at any power-on condition with respect to power-off, being at maximum at the highest thrust studied. Moreover, the lift polar curves of power-on conditions present a slope C l α higher than the power-off condition.  Regarding the lift coefficient (plots on the left in Figures 7-9), it is observed that the increment provided by the power-on conditions decreases with a decreasing d/c ratio. This result is in agreement with the one presented by Gallani et al. [29]. Hence, the configuration with d/c = 0.802 is the most promising configuration. In addition, the choice of a power trend to represent the increment on lift is suitable for the two higher diameter-chord ratios (d/c = 0.802 and 0.533) at all studied angles of attack. However, when the ratio further decreases (d/c = 0.401 and 0.319), some discrepancies appear for α = 10 • and 12 • .
Considering, now, the drag coefficient (plots on the right in Figures 7-9), a decrease in ∆C d /C d PO f f is observed when d/c decreases. Nevertheless, this decrease on the aerodynamic coefficient is not as obvious as in the case of the C l . Particularly, at α = 8 • , similar increments of the drag coefficient are obtained for configurations d/c = 0.802 and 0.533. This behaviour is also found between configurations d/c = 0.401 and 0.319 at an angle of attack equal to 12 • . Moreover, it must be mentioned that a proper fitting of the increment of drag coefficient by using a quadratic polynomial is confirmed. Only at α = 10 • for d/c = 0.319 some discrepancies appeared. From the previous considerations, it can be concluded that the configuration with the highest d/c ratio provides the best performance, greatest aerodynamic efficiency, at all angles of attack. In contrast, the configuration with the lowest d/c ratio is the worst, not only regarding performance, but also in terms of fitting. The selected curves, generally, do not fit properly the studied data.
Moreover, also interesting is the effect of the angle of attack on the increment of the aerodynamic coefficients at a fixed configuration. Specifically, this was studied at the most promising configuration (d/c = 0.802), and it is shown in Figure 10.
In terms of lift and drag coefficients, the angle of attack has barely any effect on the provided increment by the power-on conditions, as shown in Figure 10. For both aerodynamic coefficients, the maximum difference between increments at different angles of attack is found at the highest values of CTR, giving the highest increment the smallest angle of attack (α = 8 • ).    On the base of the actual results, for the Flap A condition, the following formula are proposed: The constants in Equations (1) and (2) are reported in Table 3.  The formula is tuned to better approximate the solutions at high d/c, losing accuracy towards low values of the diameter-chord ratio. It was preferred to match the results obtained at the highest values of d/c ratio as they provide the best increments of performance. A 3 − D plot of the formula in Equations (1) and (2) is shown in Figure 11. Figure 11. Surface plot of estimated increment of aerodynamic coefficients.

Flap B Working Conditions
For the second flap condition, only the configuration d/c = 0.802 is considered. This decision was made after studying Flap A conditions, since the increment of the aerodynamic coefficients of configurations with lower d/c ratios are less than those obtained for d/c = 0.802. In addition, some of the thrust was also discarded considering that the objective is to verify if the proposed trends at Flap A conditions are also anticipated for Flap B conditions.
In Figure 12, a comparison of the polar curves between power-on and power-off conditions is shown. As it was observed for the Flap A condition, both lift and drag coefficients increase at power-on conditions. Contrary to what was observed for the Flap A condition, an increase in C l α is not anticipated. Furthermore, the lift curves for the two highest thrust conditions are in a near-stall region. This behaviour is correlated to the different settings both for the larger flap deflection and to the increased gap and overlap distance. Furthermore, in Figure 13, the increment on the aerodynamic coefficients provided by the power-on conditions is going to be compared at three different angles of attack α = 8 • ,  In addition, their corresponding tendency lines are dashed.

Conclusions
This study aimed to identify and propose a simplified analytical formulation to be used during the preliminary design phase that characterises the aerodynamic performances of a wing in the presence of Distributed Electric Propulsion (DEP). Firstly, the parameters that most influence the maximum lift coefficient were identified. In order to allow a generalisation, these parameters were collected in a nondimensional representation, i.e., d/c and CTR. Then, the parametric study was carried out to find the semi-empirical formulae that characterise the effects of using DEP. Four different sets of propellers were chosen and, for each set, several thrust levels were selected. A finite wing span section with one propeller, simulated as an actuator disk, was analysed using periodic boundary conditions. With this hypothesis, an infinite wing with an infinite number of propellers was simulated.
For each d/c ratio, aerodynamic coefficients were reported as increment with reference to the corresponding power-off condition. The results showed that performance in terms of lift coefficient improves with increasing d/c ratio. A modified formulation of the Renard thrust coefficient was introduced. The increments of lift and drag coefficients plotted as function of this coefficient tend to overlap with each other at fixed d/c. With this representation, it is possible to build a fitting function.
Three angles of attack for the Flap A condition, namely, α = 8 • , 10 • , 12 • , were selected as reference and reported results showed that, for the lift coefficient, a power trend fitting curve approximates well the numerical data. A similar conclusion is achieved for drag coefficients, where a second-order polynomial function and zero intercept are proposed. A fitting surface as a function of d/c and CTR is also proposed based on available data. Flap B configuration was considered only for the highest d/c ratio. A comprehensive analysis of the results to be compared with Flap A configuration was not possible for the selected angles. In fact, for power-on condition, the data refer to the stall region of lift curve. Nevertheless, Flap B configuration presents a similar trend for the increase of the aerodynamic coefficients with reference to Flap A.
Experimental data are needed to assess this preliminary numerical activity and two experimental test campaigns are currently scheduled. Wind tunnel test campaigns will be carried out in the framework of European funded programme and, in particular, in the CICLOP project for the aerodynamic characterisation and VENUS project for the aeroacoustic one.