Aerodynamic Shape Optimization of an Aircraft Propulsor Air Intake with Boundary Layer Ingestion

: The growth of the airline industry has highlighted the need for more environmentally conscious aviation, leading to the conceptualization of more fuel-efﬁcient aircraft. One concept that has received signiﬁcant attention and has been associated with improved fuel efﬁciency is the boundary layer ingesting (BLI) propulsion system, which refers to the ingesting of the aircraft wake by the propulsors. Although BLI has theoretically been proven to reduce fuel burn, this can potentially be offset by the reduced efﬁciency and stability experienced by the propulsor in the presence of distorted inﬂow. Therefore, engine intakes must be optimized in order to mitigate the effects of BLI on the propulsion system. In this work, the shape optimization of a BLI intake is investigated using a free-form deformation technique in combination with a multi-objective genetic algorithm, in order to minimize pressure losses and distortion at the engine inlet. The optimization is performed on an S-duct intake at a cruise altitude of approximately 37,000 feet and a free stream Mach number of 0.7. An optimization strategy was developed for the task which was able to produce a Pareto optimal set of designs with improved pressure recovery and distortion. The general trend of the optimal designs shows that to reduce distortion the optimizer accelerates the ﬂow to reduce the size of the low total pressure region and increase the dynamic pressure at the engine inlet. In contrast, the pressure recovery was increased by reducing velocity as well as shifting the maximum velocity region to the outlet, which reduces the viscous dissipation losses within the intake. The ﬁnal result is a fully autonomous optimization strategy resulting in reduced pressure losses and reduced distortion leading to higher efﬁciency BLI S-duct intake designs.


Introduction
The pylon-mounted under-wing engine configuration has long remained the most popular arrangement in civil aviation. However, given the need for a more fuel-efficient generation of aircraft, the research community has recently investigated alternate configurations for the potential improvements of reducing drag and increasing propulsive efficiency. One concept at the forefront for reduced fuel burn and noise reduction is the BLI propulsion system. A BLI propulsion system is one in which the boundary layer developed over the fuselage is ingested into the propulsor and re-energized [1][2][3], as depicted in Figure 1. This concept has been around since the advent of the jet engine and can be found in ship propulsion and torpedoes [4]. However, it has yet to become a prominent fixture in commercial aviation since the integration of the propulsion system with the fuselage begets design complications.
The viscous losses within the flow field create a momentum deficit in the wake and this rate of change in momentum contributes to the aircraft drag which must be overcome by the engine thrust. A propulsion system ingesting the boundary layer re-energizes the wake flow by reducing the momentum deficit [5][6][7]. In addition, the lower momentum at the engine inlet requires less kinetic energy input to the stream [8]. This conceptualization is for the condition of zero net momentum and cruise flight. Thus, BLI effectively reduces the drag acting on the aircraft and the wasted kinetic energy for a given thrust requirement [6,8]. It is difficult to quantify the benefit of BLI given it is configuration-dependent; however, there are conceptual designs that have been studied that leverage the benefits of BLI. The D8 "double bubble" aircraft is a design developed by MIT, NASA, Aurora Flight Sciences, and Pratt & Whitney. Several studies have been conducted on the D8 "double bubble" [2,9,10], which have concluded an 8-9% reduction in required power due to BLI. Blumenthal et al. [11] performed CFD simulations on a baseline aircraft geometry with and without BLI and found that the BLI configuration reduced the engine power requirements by 8% at cruise flight.
The major challenges associated with BLI are the reduced total pressure and total pressure non-uniformity present at the engine inlet, which could lead to reduced nominal thrust output and increased fuel burn that can offset the benefits of BLI. Lucas [12] conducted experimental testing on the effects of boundary layer ingestion on a Pratt & Whitney JT15D-1 turbofan engine; replicating the total pressure distortion produced by NASA's inlet A and found that it led to a 15.5% reduction in stream thrust. Hardin et al. [13] used the numerical propulsion system simulation computational model to determine the effects of BLI and found a 0.35% loss in inlet pressure which led to a 1.1% increase in fuel burn. The results also showed there would be a 1.3% decrease in fan efficiency, also resulting in a 1.1% increase in fuel burn. Gunn and Hall [14] found that the presence of BLI causes distorted flows, creating large variations in radial and axial velocity profiles leading to a 1% to 2% reduction in stage efficiency.
BLI is generally synonymous with blended wing body (BWB) and hybrid wing body (HWB) aircraft as these conceptual designs typically embed the engines within the airframe, thus ingesting the boundary layer. The propulsion system used by these aircraft uses an S-duct to redirect the airflow. When the S-duct is not used, the engines may experience a reduced mass flow which leads to reduced nominal thrust output, causing the surge of the compression system or stall. However, the intake curvature can aggravate the pressure losses and distortion present on the engine fan face. Thus, the intakes must be optimized to minimize the pressure losses and total pressure non-uniformity at the engine inlet, resulting in increased engine efficiency and reduced fuel burn. Therefore, this paper presents a shape optimization study of an S-duct in order to minimize the pressure distortion and pressure losses at the engine inlet. The paper is organized into the following sections: Section 2 will detail the methodology used to carry out the work, followed by a discussion of the results, which is presented in Section 3, and finally, the concluding remarks are given in Section 4.

Geometry Parameterization
The S-duct geometry is defined using the concepts of B-spline surfaces to allow for shape deformation. In this paper, the S-duct is optimized for cruise conditions with a zero-degree angle of attack. Thus, geometrical and flow field symmetry are assumed in the vertical plane (in this case the Y-Z plane) during optimization. As illustrated in Figure 2, nine control points are used to define the control polygon which produces the 2nd-degree B-spline curve representing the cross-sectional shape. The control points are defined such that the curve has C1 continuity. In addition, the first and last knots have a multiplicity of p + 1, where p is the degree of the curve, to ensure the curve interpolates the two end control points. Finally, the initial shape of the duct is arbitrary since the starting design does not impact the final result; however, the outlet has a fixed diameter of two meters as the duct is designed based on Airbus A320neo engines. In addition to the cross-sectional shape, S-ducts are characterized by their centerline curvature, which is also defined using a B-spline curve as shown in Figure 3. The curvature is derived from Wellborn et al. [15], and the control points are placed to reproduce the curvature in [15] and then scaled to the required length and offset used during the optimization. An outlet diameter-to-offset ratio of 0.5 is used based on the Boeing 727 tri-jet, which incorporates a center engine that is fed through an S-duct [16]. Once the centerline curvature is established, the control points defining the crosssectional shape are placed equidistantly along the centerline, perpendicular to the local curvature which forms the control net shown in Figure 4. With the control net defined, the B-spline function generates the surface points, which are imported into ICEM CFD (ANSYS Inc., Canonsburg, PA, USA) to generate the geometry and mesh of the S-duct.
There are several other features created in ICEM CFD which are required to assess the performance of the ducts. The S-duct surface is extended upstream from the duct inlet and downstream from the duct outlet, defining the pre-and post-domain, respectively. The pre-and post-domains are simply extrusions of the inlet and outlet cross-sectional shape and are established to ensure convergence of the CFD simulations, especially as the duct is deformed through the optimization process. The deformation of the duct can introduce perturbations within the flow field which propagate upstream and/or downstream, and thus the two domains are necessary to resolve these disturbances. Figure 5 illustrates the final three-dimensional computational domain.

Geometry Deformation
The deformation process is accomplished using a free-form deformation (FFD) technique [17]. The B-spline control points defining the geometry are embedded into the FFD lattice, which is defined by the blue control points as shown in Figure 6. A local coordinate system is established such that all points within the lattice have coordinates, s, t, u, within the local coordinate system, which is normalized to vary as, 0 ≤ (s, t, u) ≤ 1. Then, as the lattice is deformed, the object embedded inside is also deformed; the physical coordinates of the embedded object can be mapped from the parametric space using Equation (1) [17]. Figure 6. B−spline control points embedded inside the FFD volume. The initial B-Spline geometry's control points are red, and the FFD control points are blue.
The X f f d vector contains the Cartesian coordinates of the deformed object, l + 1, m + 1, and k + 1 are the number of control points in the s, t, u directions, respectively, and vector P ijk contains the coordinates of the control points (i.e., weights). The lattice is a Bezier volume and thus, the perturbations of the lattice results in global shape changes to the embedded object.
The FFD volume is defined using six, three, and two control points along the z, y, and x-directions, respectively. During the shape optimization, the locations of control points are allowed to move within ±20% (±0.2 units) from their original location. This ensures that the shape optimization can be performed with sufficient deformation without extremely large deformations. In addition, the control points in the symmetry plane (X-Y plane) are only allowed to move in the symmetry plane and the control points in outlets are fixed. The other control points are limited to moving radially (i.e., normal to the duct surface), thus reducing the number of design variables to one per control point. Thus, there are five planes each consisting of five control points responsible for the deformation of the cross sections; the length of the duct is also a design variable, and, therefore, the dimensionality of the optimization problem is reduced to 26. The six planes defining the FFD volume are placed equidistantly once the length has been set.
In order to ensure deformation normal to the surface, the cross sections are deformed as a "straight-through" duct, as depicted in Figure 6. As the optimizer moves the FFD control points, the perturbations are propagated to the B-spline control points and thus, deform the geometry embedded inside the volume. The displaced B-spline control points are reoriented along the centerline curvature to form the S-duct. An overview of the optimization problem is represented in Figure 7 and Table 1.

Optimization Algorithm
In this work, MATLAB's multi-objective genetic algorithm, gamultiobj, is employed which is a gradient-free optimization method [18]. This algorithm is derived from the non-dominating sorting genetic algorithm (NSGA-II) and can provide a Pareto optimal set of solutions. In contrast to the NSGA-II algorithm, gamultiobj, adheres to a controlled elitist technique which helps maintain a higher level of diversity within the population ensuring a more adequate exploration of the design space [18].
The progression of the algorithm towards optimality starts with the generation of the initial population [18]. See algorithm flow chart in Figure 8. As opposed to gradient-based methods, the genetic algorithms use an initial population rather than an initial design point, this effectively limits the dependence of the optimal solution on the initial starting design. The initial population (generation zero) is randomly generated, within the linear constraints of the given problem, to ensure proper representation of the design space. Next, the performance of each individual within the population is evaluated and then ranked based on their performance with a certain fraction of the top-ranking individuals (i.e., the elite) passed on to subsequent generations. To maintain the population level, the remaining individuals of the subsequent population are generated according to the selection, crossover, and mutation processes. In the selection process, individuals undergo a tournament to determine which individual is the top performer who will be selected as a parent. The parent variables are combined to form a child, in the crossover process, in the hopes that the child will also be an elite individual. Additionally, the mutation process introduces random perturbations to the design variables to ensure that there is an adequate representation of the design space [18]. These processes are repeated for subsequent generations until the stopping criterion is reached, which is a maximum of 50 generations in this work. The population size is also set to 50; both values were found to be adequate to effectively explore the design space. A numerical optimization strategy was implemented to optimize BLI S-duct intake. The length and cross-sectional shape are allowed to vary during the optimization whereas the offset and exit diameter are fixed. MATLAB was used to build the optimization framework using the multi-objective genetic algorithm. MATLAB was also used to execute, process, and collect results from other software such as ICEM, CFD Pre, CFX, and CFD Post. Figure 9 shows the general implementation of the optimization process used in this paper which was also used in [19]. The optimization algorithm modifies the design variables, i.e., the location of the FFD control points, which is used to generate new geometry and progresses the design towards optimality. The deformed cross sections generated by the FFD process are then used by ICEM to generate the mesh of the fluid domain. Then CFD Pre is used to apply the boundary conditions for the symmetry plane, inlet, and outlet. Next, CFX is used to solve the flow, and CFD post is used to calculate the pressure recovery and distortion coefficient. These values are then used by MATLAB to search for the optimal S-duct shape.

Computational Methods
The performances of all design iterations are evaluated using the steady-state Reynoldsaveraged Navier-Stokes (RANS) equations and the shear stress transport (SST) turbulence model of the commercial CFD software, ANSYS CFX. The CFD simulations are performed on the computational mesh shown in Figure 10, which consists of a H-O topology and is generated in ICEM CFD using a robust blocking strategy. The mesh consists of around 370,000 elements with a y+ value of around 10-20. Table 2 shows the mesh sensitivity study for an undeformed duct of eight meters, the change in PR and distortion with an increasing number of elements is minimal and thus, this mesh is determined to be adequate to represent the flow field within the S-duct while keeping the computational time low. Since the mesh defines the computational space and ultimately the design space, the characteristics of the mesh are consistent across all designs.  The S-duct optimization is carried out for a free stream Mach number of 0.7 with a boundary layer thickness of 0.3 m (~one foot). The boundary layer thickness of 0.24 m was obtained by calculating the boundary layer thickness of a flow over a flat plate with a length equivalent to an Airbus A320neo aircraft. This formula assumes a smooth surface without friction; thus, the boundary layer thickness was increased to 0.3 m. The boundary condition for the S-duct inlet is obtained by performing a simulation of a flow over a two-dimensional flat plate. The flow variables are then interpolated and translated slightly (along the y-axis) to fit the pre-domain inlet. In addition to the total pressure profile, the associated turbulent kinetic energy and turbulent eddy frequency profiles are also applied to ensure proper representation of the boundary layer. The outlet static pressure is initially set to the atmospheric pressure at the flight altitude, 37,000 feet. However, it was found that the outlet static pressure was strongly correlated to the inlet area. Thus, a relationship between the outlet static pressure and inlet area was established, as shown in Figure 11, which was used to set the outlet static pressure after the geometry was deformed. Finally, the pre-and post-domain are set to a no-slip condition as well. The total temperature at the inlet is set to the temperature at the cruise altitude. A summary of the boundary conditions is provided in Table 3. Table 3. Boundary conditions used for the simulations.

Boundary Boundary Condition Description
Inlet (Pre-domain) Total Pressure Profile (See Figure 12 The CFD simulation was set to run for a maximum of 200 iterations. The convergence criteria were set such that when both the residuals and system imbalance are below 10 −4 and 1%, respectively, the simulations would stop. A 2nd-order upwind biased scheme is used for spatial discretization. These parameters ensure an adequate convergence for the steady-state CFD simulation. Additionally, the timescale factor is set to four to reduce computational time without sacrificing solution accuracy. As with all RANS turbulence models, the inaccuracy introduced by turbulence modeling remains. However, it is believed that while there might be inaccuracies in the absolute magnitudes, the results are trendwise accurate, leading the optimization algorithm in the right direction to arrive at the optimized design.

Objective Functions
The objective of the optimization is to minimize the pressure losses within the duct and the pressure non-uniformity at the engine inlet to ensure minimal impact on the engine performance. Therefore, the performance of each design is assessed with respect to the pressure recovery, PR, and total pressure distortion, DC(60), at the engine inlet as the optimization algorithm progresses the designs towards optimality. Equation (2) is used to calculate the PR of the duct.
MATLAB's optimization function attempts to minimize the objective function. Thus, to maximize the PR (Equation (2)), the reciprocal of PR was used as an objective function to be minimized. The pressure values used to calculate the PR at the inlet and outlet are area-averaged to account for the non-uniformity of the mesh.
The second objective is to minimize the variation of the pressure at the engine inlet. A common performance metric to evaluate the total pressure uniformity is the DC(60) distortion coefficient given in Equation (3) [20].
where q is the dynamic pressure and P t, 60 • , min is the 60 • sector with the worst averaged total pressure. Minimizing the DC leads to more uniform total pressure distribution at the engine intake. Figure 13 shows how the 60 • sector is rotated to calculate the distortion coefficient. Initially, the area, pressure, and angular position for each element at the outlet face are calculated. The area-averaged pressure is calculated for a 60 • sector. Then the 60 • sector is rotated by 5 • increments. Then 60 • sector with the worst average pressure is used for P t, 60 • , min to calculate DC(60).

Results and Discussion
The optimization took a total of four days on a workstation with 16 cores and 64 gigabytes (GB) of RAM. In the following sections, the total pressure distribution at the duct outlet (also referred to as the engine inlet; Figure 5) is presented as a dimensionless value (normalized total pressure). The dimensionless value is obtained by dividing the nodal pressure values by the area-averaged pressure at the engine inlet. An experimental study on the baseline geometry was carried out and showed good agreement with the numerical results [21].

S-Duct Optimization
The optimization of the S-duct is carried out with respect to the parameters given in Table 1, which produced the corresponding evaluation history and Pareto front, as shown in Figure 14, after more than 1800 simulations. A subset of the Pareto front is marked and the normalized pressure distribution at the engine inlet is visualized. In addition, the objective function values for an undeformed duct of eight meters are also plotted to provide a reference for the improvement in performance; the reference point is not a starting point since genetic algorithms operate on a population (i.e., multiple designs) rather than an initial design point. The designs on the Pareto front have varying lengths with the length increasing as distortion decreases. The main reason that the results shown in Figure 14 seem to appear in columns is due to the lack of precision of the pressure recovery values and the smaller range of PR magnified over a large plot size. The subset of the Pareto optimal solutions is summarized in Table 4.  As evident from Table 4, the Pareto optimal designs have relatively substantial improvements in distortion and PR relative to the reference geometry. The optimal PR duct provides a 0.61% increase in PR compared to an undeformed duct of eight meters. On the other hand, the optimized duct for distortion improved the distortion metric by 36.4%, significantly reducing the overall pressure variation of the low total pressure region. Although the difference in pressure loss is minimal compared to the distortion metric, these small PR losses equate to improved fuel consumption. Hardin et al. [13] stated that a 0.35% increase in pressure losses can lead to an approximately 1.1% increase in fuel burn. Figure 15 shows a comparison of Mach number contours and pressure contours at the symmetry planes and engine inlet, respectively, for the baseline and optimal distortion duct. As seen in Figure 15, the low-pressure region at the bottom of the duct is similar in shape, but the magnitude of the area-averaged total pressure is greater for the optimal distortion duct. The increase in magnitude for each sector is reflected as a lower calculated distortion coefficient. The length of the duct was increased by the optimizer to reduce the adverse pressure gradient in the flow. This leads to a reduction in the momentum loss within the boundary layer flow which reduces the distortion present at the engine inlet. Additionally, the region of maximum velocity was shifted to the engine inlet, which increases the dynamic pressure. The higher kinetic energy of the flow leads to a lower distortion thus, reducing the susceptibility to stalls of the compression system.
The flow field within the optimal PR duct is similar to that of the optimal distortion duct as well as the total pressure distribution at the engine inlet. However, the low total pressure region is slightly larger due to the slight separation occurring at the inflection point on the inner curvature. As opposed to the undeformed duct of eight meters, the optimal PR duct is extended to nine meters to balance the acceleration of the flow over the inner and outer bends and the increased viscous dissipation due to the lengthening of the duct. The increase in length reduces the acceleration of the flow in the bends, which reduces the wall shear stresses and the pressure losses. In addition, similar to the optimal distortion duct, the region of maximum velocity is shifted to the duct outlet, this reduces the flow acceleration within the duct, leading to a higher PR. To the same extent, the velocity through the duct is generally slower compared to the undeformed duct which also contributes to the reduced viscous dissipation. The cross-sectional shapes of the optimal distortion and PR ducts are shown in Figure 16. The cross-sectional deformations of both ducts are very similar, the intermediate cross sections have been widened, which serves two purposes. First, as a consequence of the fixed outlet, the widening of the cross sections accelerates the flow as it progresses towards the outlet, increasing the dynamic pressure at the duct outlet and thus, reducing the distortion. Secondly, the wider cross-sectional area reduces the flow velocity within the duct relative to the undeformed case, as depicted in Figures 15 and 17, and shifts the region of maximum velocity to the duct outlet. Thus, this is the reason for the increased PR in both optimized ducts. The slight difference occurs towards the inlet, where the bottom of the cross section has been "pulled down". This narrowing of the bottom of the inlet area is to reduce the amount of BLI. However, if the area is narrowed too aggressively, the sudden expansion in the following cross sections will retard the boundary layer momentum and thus, lead to flow separation. The two objective functions seem to be linked in terms of cross-sectional deformation, as both require expansion in the intermediate areas. The difference between the two ducts is the length, a longer length for the optimal distortion duct reduces the adverse pressure gradient, however, would increase pressure losses due to viscous dissipation.

Effects of the Centerline Curvature
The S-duct centerline is defined using five control points as shown in Figure 3. The centerline curvature of the S-duct is investigated here because it has a large effect on the flow field. The centerline is deformed by moving the control points shown in Figure 3, which defines a B-spline curve. The control points CP2 and CP4 are constrained to move only in the z-direction (one degree of freedom), while CP3 can be moved in both the y-and z-directions (two degrees of freedom). The control points were allowed to move ±50% (of the total length of the duct), from their original location because greater freedom resulted in poor mesh quality and non-convergence of the CFD simulations. The deformation of the centerline is applied to the optimal distortion duct obtained from the S-duct optimization as this duct provided a good balance between PR and the distortion present at the engine inlet. As such, the cross sections and length were held fixed while only allowing deformation of the centerline.
The evaluation history for the centerline curvature and a subset of the Pareto front's engine inlet pressure visualizations are presented in Figure 18.  Figure 18 shows that there are two distinct developments of the Pareto front. The first Pareto front on the right side has substantially better PR; whereas the other Pareto front favors distortion. The Pareto front thumbnails show two distinct pressure distributions. The Pareto front at the right, which has higher PR, shows the distortion is at the six o'clock position; whereas the second Pareto front shows the distortion is shifted towards the center of the outlet. The optimized designs showed significant improvements in both maximizing PR and in the reduction of DC. A subset of the Pareto front shown in Figure 18 is tabulated in Table 5. In Figure 19, the Mach contours (at the symmetry plane) and pressure contours (at the engine inlet) are shown for optimal PR, optimal DC, and baseline S-duct geometry. As shown, the two optimal centerline profiles shown in Figure 19 are very different. The optimal distortion duct's centerline curvature is extremely aggressive almost resulting in a 90 • bend. The aggressive nature of the curvature causes flow separation at the lower curvature (before the inflection point) of the duct and a pocket of supersonic flow over the upper curvature (near the inflection point). The inflection point of the centerline curvature is moved near the duct inlet leading to an almost straight duct section prior to the engine inlet. The flow separation with the changes in the duct curvature causes a pressure differential that produces counter-rotating vortices (also known as Dean vortices), which shifts the low total pressure region towards the center of the outlet. The shift leads to a more uniform circumferential total pressure profile which is distributed more evenly across all the sectors of the outlet, improving the DC metric at the engine inlet. However, the flow separation and the pocket of supersonic flow both cause pressure loss at the engine inlet, which reduces the PR metric. As shown in Figure 19, the optimal PR duct has an elongated arc; closer to a straight path from the duct inlet to the engine inlet rather than an S-shaped curve. Interestingly, Gan and Zhang [22] also found the optimal centerline curvature to be an elongated arc when performing a shape optimization study of a diffusing S-duct. The changes to the curvature from the initial S-shaped curvature to an elongated arc reduce the acceleration experienced by the flow. Thus, the flow has a relatively constant velocity profile and has reduced pressure loss. Due to the design on the design of the downstream post-domain, which is just an extruded engine outlet profile, there is a sharp discontinuity in the fluid domain mesh. This results in flow separation and supersonic flow conditions at the outlet (engine inlet) and downstream section (post-domain) of the duct. This has minimal impact on the performance of the duct, but flow separation and supersonic flow would negatively impact the performance of the propulsor.

Conclusions
A fully autonomous optimization strategy has been implemented using commercially available software and applied to BLI S-duct intakes. The S-duct optimization was able to produce designs with improved PR and distortion. The reduction in distortion was accomplished by accelerating the flow and shifting the region of maximum velocity towards the duct outlet, increasing the dynamic pressure at the engine inlet. This is beneficial for the propulsor as a higher inlet velocity would reduce the likelihood of flow separation over the blades and thus the susceptibility to stalls. Additionally, the duct was extended to the maximum allowable length in order to reduce the adverse pressure gradient thus, reducing the deceleration of the boundary layer flow and the susceptibility to separation. The PR was increased by shortening the duct as much as possible before causing the flow to separate and widening the cross sections to reduce the velocity within the duct as well as shifting the maximum velocity region to the outlet, which reduces the viscous losses due to flow acceleration.
The sensitivity of the optimal designs to the centerline curvature was also investigated, which indicated a significant coupling between the centerline and the duct performance, as expected. This study revealed that for increased PR, the centerline is no longer "S-shaped" and is representative of an elongated arc, almost producing a straight path from the inlet to the outlet. This was done in order to reduce the flow acceleration through the duct, which reduces the pressure losses due to friction. In contrast, the optimal distortion duct had an aggressive S-shape with the inflection point moved upstream to allow the counter-rotating vortices to shift the low total pressure region towards the center of the outlet, providing a more uniform total pressure distribution.