*2.2. Optimization Problem*

The surrogate-based framework is applied to the I-31T aircraft air-intake duct multipoint optimization problem. The I-31T flight demonstrator is a retrofitted make of the I-23 "Manager" airplane. The plane was designed in Łukasiewicz Research Network—Institute of Aviation [67] and modified under the EU R&D program ESPOSA (Efficient Systems and Propulsion for Small Aircraft) [68]. The redesign had in scope a replacement of the piston motor with the PBS TP100 turboprop engine [69–72]. Integration of the nominally pusherpurpose engine with the tractor-type airplane required the introduction of a tailored air delivery U-shaped duct. The duct was designed by Stalewski and Z˙ ółtak [73] using a parametric design strategy [74,75] relying on the NURBS-based in-house code PARADES [76]. This study takes the outcome of the abovementioned work as an initial point and aims to further improve the duct's performance characteristics in various flight conditions.

Figure 3 presents the initial air-intake duct geometry. The upstream portion of the channel is slightly S-shaped and transitions into a U-shaped duct with the flow direction. Further downstream, the U-duct interfaces with the engine compressor front face at the location marked as the Aerodynamic Interface Plane (AIP). The intake channel has a slightly converging character with a hydraulic diameter equal to 0.180 m at the duct's entrance and 0.147 m at the AIP. The dummy extension, of a length of two diameters, is an artificial portion introduced to secure the flow solver's stability. The duct's inlet and AIP shapes and positions are constrained, i.e., not subject to shape deformation.

The surface of the duct in Figure 3 shows the initial computational mesh. The nine points surrounding the geometry indicate the baseline location of the mesh morpher control points subject to displacement during optimization. They are located in regions where presumed deformation should have a significant impact on the duct performance. The horizontal and vertical translation vector components form a set of design variables. Three pairs of points are bounded so that each duplet shares the exact translation. This strategy guarantees the maintenance of duct symmetry. Such a formulation results in twelve design variables, shown in Figure 3, whose ranges are defined in Table 1. The ranges corresponding to control points close to the AIP are narrower to avoid excessive

deformations in the vicinity of the constrained section. Such distorted geometries would have a low probability of generating superior solutions.

**Figure 3.** Initial mesh with an indication of flow stations and mesh morpher control points' locations.

**Table 1.** Ranges of design variables.


The amount of pressure loss along the air intake and the level of circumferential pressure distortion at the compressor entrance substantially influence the ultimate performance and operability of the airplane engine. These two objectives are subject to minimization in the present study:

$$\text{pressure loss coefficient}: f\_1 \equiv dP = \frac{P\_t^{IN} - P\_t^{AIP}}{P\_t^{IN}} \tag{25}$$

$$\text{Listortion coefficient}: f\_2 \equiv D \mathcal{C}\_{60} = \frac{P\_t^{AIP} - \min\_{\theta \in \left[0, 2\pi\right]} P\_t^{AIP} \left(\theta - \frac{\pi}{6}, \theta + \frac{\pi}{6}\right)}{q^{AIP}} \tag{26}$$

In the above equations, *PIN <sup>t</sup>* and *PAIP <sup>t</sup>* denote total pressure values mass-averaged over the inlet and AIP cross-sections, respectively, and *qAIP* indicates dynamic pressure mass-averaged across the AIP. The term min *θ*∈[0,2*π*] *PAIP t <sup>θ</sup>* <sup>−</sup> *<sup>π</sup>* <sup>6</sup> , *<sup>θ</sup>* <sup>+</sup> *<sup>π</sup>* 6 refers to the lowest average pressure in the 60◦ sector, which is proven to have the most detrimental impact on the compressor performance [77]. This evaluation is realized by clocking a virtual sector around the AIP with Δ*θ* = 10◦ angular intervals—as given schematically in Figure 4.

**Figure 4.** Illustration of min *θ*∈[0,2*π*] *PAIP t <sup>θ</sup>* <sup>−</sup> *<sup>π</sup>* <sup>6</sup> , *<sup>θ</sup>* <sup>+</sup> *<sup>π</sup>* 6 evaluation process at the AIP. The blue segment represents the pressure averaging region.

The multi-point optimization problem considers three I-31T airplane flight conditions: DP1—nominal cruise, DP2—low-altitude climb, and DP3—high-altitude cruise. Details of the design points and related environmental properties are gathered in Table 2. The corresponding targets and *nadir* points are given in Table 3. The goal for the nominal cruise has a priority over the off-design conditions to reflect the contribution of this design point to the aggregated efficiency over the entire flight mission. Both performance objectives are considered to be of equal importance and share similar targets for particular design points. On the grounds of the designer's experience, the possible falloff of 200% for each objective gives a basis for the function's upper bound estimation.

**Table 2.** Details of operating conditions for considered design points.


Ambient conditions were evaluated with use of the U.S. Standard Atmosphere model [78].


**Table 3.** Target and nadir points for objectives under selected flight conditions.

*2.3. Evaluation of Objectives*

2.3.1. Flow Solver Governing Equations

The values of the objective functions are derived from the flow field solution of the air-intake duct CFD simulations. For this purpose, we use a Reynolds-averaged Navier– Stokes (RANS) solver capable of modeling three-dimensional turbulent viscous flow to accurately predict flow features characterized by a secondary flow motion, strong pressure gradients, and occurrence of separation. Such flow features are highly expected in channels of high curvature [79]. The solver is implemented using a commercial CFD code: ANSYS CFX Release 18.0.

In this study, the code solves the following Favre-averaged Navier–Stokes conservation equations:

• The continuity equation:

$$\frac{\partial \overline{\rho}}{\partial t} + \frac{\partial}{\partial x\_i} (\overline{\rho} \overline{u}\_i) = 0 \tag{27}$$

• The momentum conservation equation:

$$\frac{\partial}{\partial t}(\overline{\rho}\overline{u}\_i) + \frac{\partial}{\partial \mathbf{x}\_j} \left(\overline{\rho}\overline{u}\_i\overline{u}\_j\right) = \frac{\partial}{\partial \mathbf{x}\_j} \left(-p\delta\_{ij} + \overline{\mathbf{r}\_{ij}} + \overline{\rho}\mathbf{r}\_{ij}\right) \tag{28}$$

In the above formulations, *ρ*, *ui*, and *p* denote fluid density, velocity components, and static pressure, respectively. The tensor *tij* indicates the viscous shear stress component. The operators (·) and ( !·) stand for the Reynolds-averaged and Favre-averaged mean variables, respectively [80]. The presence of the turbulent (Reynolds) stress tensor *ρτij* = −*ρu <sup>i</sup> u <sup>j</sup>* requires additional closure for the system of RANS equations. By employing the Boussinesq hypothesis, which is the foundation of eddy viscosity-based turbulence models, the Reynolds stress tensor is approximated as:

$$\overline{\rho}\tau\_{\overline{i}\overline{j}} = -\overline{\rho u\_i'' u\_j''} \cong \mu\_l \left( \frac{\partial \overline{u}\_i}{\partial x\_{\overline{j}}} + \frac{\partial \overline{u}\_{\overline{j}}}{\partial x\_{\overline{i}}} - \frac{2}{3} \frac{\partial \overline{u}\_k}{\partial x\_k} \delta\_{\overline{i}\overline{j}} \right) - \frac{2}{3} \overline{\rho}k \delta\_{\overline{i}\overline{j}} \tag{29}$$

where *μ<sup>t</sup>* is the turbulent viscosity determined using the selected turbulence model. Here, we use the two-equation SST k-ω turbulence model by Menter [81,82]. This model computes *μ<sup>t</sup>* from a solution of additional transport equations for two variables: the turbulent kinetic energy *k* and the specific dissipation rate *ω* (Equations (30) and (31)).

$$\frac{\partial \rho}{\partial t}(\overline{\rho}k) + \frac{\partial}{\partial x\_j}(\overline{\rho}\overline{u}\_j k) = \overline{\rho}\tau\_{i\bar{j}}\frac{\partial \overline{u}\_i}{\partial x\_j} - \beta^\*\overline{\rho}k\omega + \frac{\partial}{\partial x\_j}\left[ (\mu + \sigma\_k \mu\_t)\frac{\partial k}{\partial x\_j} \right] \tag{30}$$

$$\frac{\partial}{\partial t}(\overline{\rho}\omega) + \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}}(\overline{\rho}\overline{u}\_{\dot{j}}\omega) = \frac{\mathbf{a}}{\nu\_{l}}\overline{\rho}\mathbf{r}\_{\overline{\eta}}\frac{\partial \overline{u}\_{\dot{i}}}{\partial \mathbf{x}\_{\dot{j}}} - \beta\overline{\rho}\omega^{2} + \frac{\partial}{\partial \mathbf{x}\_{\dot{j}}}\bigg[ (\mu + \sigma\_{\omega}\mu\_{l})\frac{\partial \omega}{\partial \mathbf{x}\_{\dot{j}}} \bigg] + 2(1 - F\_{l})\frac{\overline{\rho}\sigma\_{\omega 2}}{\omega} \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_{\dot{j}}} \frac{\partial \omega}{\partial \mathbf{x}\_{\dot{j}}} \tag{31}$$

Details of the SST k-ω model formulation used in this paper and a specification of the closure coefficients are available in [83].

The system of transport equations is closed with the ideal gas law equation, solved for the density (Equation (32)).

$$
\overline{\varphi} = \overline{\varphi}RT\tag{32}
$$

Here, *R* is the specific gas constant for air, and *T* is the fluid temperature. This study employs an isothermal model, for which the temperature is selected and kept constant according to the specifics of particular design points. Such an assumption is well grounded for a low Ma number flow regime and the absence of heat sources, which is accurate for the intake duct conditions.

#### 2.3.2. Validation of Turbulence Modeling Technique

The selection of the SST k-ω turbulence model was validated prior to the optimization study. The model's predictive quality was compared with the eddy viscosity-based k-ε turbulence model and the non-linear k-ε Explicit Algebraic Reynolds Stress Model (k-ε EARSM). All three models were considered in two options: with and without the curvature correction (CC) term. The flow field predictions for the six abovementioned models were compared with the experimental data for a benchmark U-duct provided by Azzola et al. [84,85]. Details of the validation study are available in [37].

Figure 5 groups the study results for the longitudinal (U) and circumferential (V) velocity components normalized by the bulk velocity (Ub) at the entrance to the duct's curved portion. The k-ε turbulence model showed the most drastic underprediction in the core flow (r/D = 0) deceleration at the duct bend (see Figure 5a). Quite oppositely, the k-ε EARSM model failed to predict the circumferential flow motion at the station downstream of the curved portion. Overall, the study showed the best agreement between the flow simulations and the experimental reference data for the SST k-ω turbulence model. Moreover, no significant beneficial influence of the curvature correction term was detected. Following this reasoning, the SST k-ω model with deactivated CC term is employed in this work.

**Figure 5.** Prediction of normalized longitudinal (U) and circumferential (V) velocity components for various turbulence models in a 180◦ curved duct with a circular cross-section. Station: (**a**,**b**) at 90◦ bend; (**c**,**d**) two diameters downstream of the curved section. Adapted from [37].

### 2.3.3. Details of the Flow Field Modeling Approach

Steady Navier–Stokes equations are discretized by the element-based finite volume method, with a second-order upwind approximation of convective fluxes and a second-order central representation of diffusive fluxes. An additional limiter using the Barth and Jespersen principles [86] is activated to maintain the solution's boundedness. The solution to the equations is obtained using the coupled pressure-based algorithm iterated until the normalized residuals of the momentum and transport equations drop below 10<sup>−</sup>5.

The boundary conditions mimic the mission characteristics given in Table 2. The duct's inlet flow conditions are defined, with the total pressure value calculated using Bernoulli's equation for specific flight speed and altitude. Furthermore, the inlet turbulence intensity is set to *Tu* = 10%, and the eddy-to-molecular viscosity ratio equals *νt*/*ν* = 100. The amount of airflow accepted by the engine at specific design points is reflected through the mass flow rate condition at the domain outlet. The duct's walls obey the non-slip boundary condition.

The flow domain corresponding to the duct's interior is discretized using a structured hexahedral mesh. To secure a high-level resolution in the near-wall region accompanied by quality elements in the flow core, the mesh is built in the hybrid O-H grid topology (see Figure 6). On average, the non-dimensional wall distance is *<sup>y</sup>*<sup>+</sup> <sup>∼</sup> 1.7 with extreme variation up to *<sup>y</sup>*<sup>+</sup> <sup>≤</sup> 3.5. These values hold throughout the optimization process.

**Figure 6.** O-H grid at the duct's AIP.

The execution of the comparative study secures the optimization results' independence from the mesh size. Four different grid sizes are compared against the duct-wise pressure distribution. The mass-averaged total pressure values are normalized using the duct inlet pressure. The results are compared in Figure 7. The two finest grids (1,451,000 and 8,007,000 cells) show no difference in predictions for the vast majority of the duct length (~95%). A minor deviation is observable close to the duct's outlet. Hence, the mesh with 1,451,000 elements is used as it has sufficient prediction quality and protects the computational economy of the simulations.

**Figure 7.** Results of the duct grid independence study. L = duct centerline length. Numerical values indicate the number of mesh elements.

#### **3. Results and Discussion**

#### *3.1. Metamodel Fit Validation*

Before the optimization began, the Kriging surrogate was assessed for its ability to generalize given data, i.e., to predict objective values at unobserved locations. For this purpose, we used the leave-one-out cross-validation (LOO-CV) procedure, described in detail in Appendix A.

The validation process was executed for the nominal ASF formulation and its mathematical transformations given in Table 4. This action anticipates potential deficiencies in the surrogate's prediction accuracy. In such a scenario, a deterministic mathematical function can be applied to the data set to deliberately impact its distribution and potentially improve the metamodel's fit quality.



The outcome of the LOO-CV is presented graphically using diagnostic plots in Figures 8–10. In Figure 8, values on the horizontal axis denote exact evaluations of each observation, while the vertical axis shows the surrogate prediction when the corresponding sample is omitted. The markers would be distributed perfectly on a diagonal in a utopian metamodel with zero-error predictions. A realizable high-performing surrogate would be characterized by samples located in close proximity to the imaginary line. Figure 8a shows the results for the nominal ASF. The samples approximately follow the diagonal and are reasonably distributed along the reference line. Only a few observations slightly depart from the ideal arrangement. The results for transformed data (see Figure 8b–f) present a similar pattern, with a slight tendency in some functions to cluster around the higher values (e.g., Figure 8d,e).

**Figure 8.** LOO−CV plots for various objective function transformations. Blue lines indicate locations of perfect predictions.

**Figure 9.** Plots of standardized cross-validated residuals for various objective function transformations. Red dashed lines indicate interval of +/− 3 standard deviations.

Figure 9 shows the distribution of the cross-validated residuals, which—as discussed in Appendix A—should be bounded in a ±3 interval. Inspection of the residuals for nominal ASF (see Figure 9a) reveals one observation that slightly exceeds the desired limits with a residual value of ~3.5. Except for that fault, the residuals are randomly distributed along the sample values with no apparent pattern. The abovementioned deficiency does not impose a rejection of the model, although it might adversely affect its prediction accuracy. Most of the considered data transformations (see Figure 9b–e) do not remove the outlier, although the negative inverse function reduces it to the threshold value. Quite oppositely, the cube root transformation, given in Figure 9f, reduces the error variance to a healthy range without affecting randomized distribution.

The quantile-quantile plot (Q-Q plot) in Figure 10 helps to examine the standardized residuals' normality and homoscedasticity. The assumption of normal distribution should be satisfied for the legitimate use of MLE in the surrogate's hyperparameter estimation. Moreover, following the Kriging principles, the residuals should be homoscedastic, i.e., the variance ought to be homogeneous across the design space [87]. The Q-Q plot draws the residuals' quantiles sorted in ascending order (vertical axis) against the quantiles from the theoretical Gaussian distribution (horizontal axis). For normally distributed residuals, markers on the plot would follow the diagonal. A random spread of points along this line manifests residuals' homoscedasticity.

**Figure 10.** Q−Q plot for various objective function transformations. Blue lines indicate target normal distribution.

The Q-Q plot for the original ASF (see Figure 10a) reveals that the extreme positive residuals tend to deviate counter-clockwise relative to the reference line. Such a pattern indicates that more data are located on the range's right bound, i.e., the distribution is 'heavy-tailed'. Moreover, the lack of a similar pattern on the left tail suggests asymmetry in the distribution with a tendency to right-skewness. Inspection of the remaining Q-Q plots reveals that only two transformations, negative inverse square root (Figure 10d) and cube root (Figure 10f), normalize the distribution of the standardized residuals.

As the cube root transformation is superior in the normalization of residuals and maintains the error variance in a healthy span, we use it to transform the ASF in further studies. Moreover, this transformation is insensitive to the sign of input, which supports the robustness of the process. Although this is an optimistic scenario, the ASF formulation potentially produces negative values for solutions superior to the presumed target. Although less popular than other considered transformations, the cube root has been recognized in applications for long-tailed data [88].

#### *3.2. Sensitivity Analysis*

A sensitivity analysis was executed to assess how strong is the dependence of the surrogate's output (i.e., the objective function prediction) on variance in particular design variables. As a result, less impactful variables could be identified and removed from further consideration, thus reducing the computational cost of the optimization.

We used Functional Analysis of Variance (FANOVA), described in detail in Appendix B. This method assesses the contribution of particular inputs and interactions between them using the notion of Sobol' indices [89,90]. FANOVA is particularly beneficial

in use with surrogates that might be non-linear and non-monotonic but are inexpensive in the output estimation.

Figure 11 displays graphically the outcome of the FANOVA analysis applied to the Kriging surrogate constructed on the DoE observations. The height of the bars represents the value of the total effect Sobol' indices, while black and gray bars' fractions denote the main effect and aggregated higher-order interactions, respectively.

**Figure 11.** The effect of FANOVA analysis on model considering 12 variables. Bar height indicates the total effect Sobol' index value; black color represents the main effect index, and gray color denotes all higher-order interactions.

Undoubtedly, the design variable *x*<sup>11</sup> has the most substantial influence on variance in the objective function. The second most impactful is the variable *x*12. Both factors correspond to the spatial motion of the control point located in the duct's inner wall area (see Figure 3). Except for the high total index value, these variables are characterized by a strong direct effect, suggesting their crucial influence on the final optimization results.

Quite oppositely, the impact of variables {*x*1, *x*5, *x*7, *x*10} is less significant. Their already low total effect value is related mainly to the higher-order interactions rather than to the direct impact. For this reason, these variables are removed from the design space, notably reducing its dimensions and increasing the process's cost-efficiency. Variable *x*<sup>6</sup> presents a similar level of total impact to *x*7, although its main influence is slightly more prominent. Contrarily, the direct effect of *x*<sup>9</sup> is minor; however, strong higher-order interactions result in a meaningful value of the total effect index. It depends on the designer's intuition whether such variables should be further considered; here, we preserve them in the study.
