2.1 Geometry definition

Kulfan [Reference Kulfan and Bussoletti30, Reference Kulfan31] proposed a unique and powerful geometry representation method called class/shape transformation (CST). The analytical CST geometry representation methodology provides a unified and systematic approach to present a wide variety of 2D and 3D geometries, encompassing a significantly large design space with relatively few parametres. Zhu and Qin [Reference Zhu and Qin32] extended the CST method and proposed a new aerofoil parameterisation method called intuitive class/shape transformation (iCST), which combines the flexibility and accuracy of the CST method with the intuitiveness of the PARSEC method. Christie et al. [Reference Christie, Heidebrecht and MacManus33] developed a system that analytically calculates the transformation matrix for a set of geometric constraints based on the iCST method and successfully applied it to construct a 3D axisymmetric nacelle. The iCST method has been successfully applied to parameterise various aerodynamic shape optimisation problems [Reference Goulos, Stankowski and Otter34].

In this study, the iCST method was applied to parameterise the surface of the fan cowl and intake. This method consists of shape and class functions. The class function is used to define the general classes of geometry, whereas the shape function is used to define specific shapes within the geometry class. The general mathematical expression of CST is of the following form:

(1) \begin{align}\zeta \!\left( \psi \right) = C_{{N_2}}^{{N_1}}\!\left( \psi \right)S\!\left( \psi \right) + \psi {\rm{\Delta }}{\zeta _{TE}};{\rm{\;\;\;\;}}\zeta = \frac{z}{c};{\rm{\;\;\;\;}}\psi = \frac{x}{c}\end{align}

where $C_{{N_2}}^{{N_1}}$ is the class function, $S\!\left( \psi \right)$ is the shape function, and $\psi {\rm{\Delta }}{\zeta _{TE}}$ controls the trailing-edge thickness. The mathematical expression of the class function is of the following form:

(2) \begin{align}C_{{N_2}}^{{N_1}}\!\left( \psi \right) = {\psi ^{{N_1}}}{(1 - \psi )^{{N_2}}}\end{align}

where ${N_1}$ and ${N_2}$ are tuned to obtain different shapes, as shown in Fig. 2. For a round-nose aerofoil, usually ${N_1} = 0.5$ and ${N_2} = 1.0$ .

Figure 2. Illustration of the effects of ${N_1}$ and ${N_2}$ on the class function.

The shape function was constructed using a Bernstein polynomial with different weight coefficients. The mathematical expression of the Bernstein polynomials is

(3) \begin{align}BP\!\left( \psi \right) = \mathop \sum \limits_{i = 0}^N \!\left[ {{K_{i,N}} \cdot \!\left( {{\psi ^i} \cdot {{(1 - \psi )}^{N - i}}} \right)} \right];\; {K_{i,N}} = {{N!} \over {i!\!\left( {N - i} \right)!}}\end{align}

In addition, the shape function used in this study can be expressed as

(4) \begin{align}S\!\left( \psi \right) = \mathop \sum \limits_{i = 0}^N \!\left[ {{A_i} \cdot {K_{i,N}} \cdot \!\left( {{\psi ^i} \cdot {{(1 - \psi )}^{N - i}}} \right)} \right]\end{align}

In this work, the open source common research model (CRM) was used and was modified based on the fan diameter to meet the short nacelle $L/D = 0.35$ requirement. The target geometry can be obtained using appropriate weight coefficients. For an axisymmetric nacelle geometry, the aero-lines of the fancowl and intake were first created using the iCST method and then revolved to construct the surface. The aero-lines of the fancowl and intake were defined by several variables, as shown in Fig. 3. To simplify the optimisation process and comply with the specific design requirements, six design variables ( ${r_{if}},{r_{{\rm{max}}}},{f_{{\rm{max}}}},{\beta _{{\rm{nac}}}},{r_i}$ , and ${f_i}$ ) were chosen to parameterise a single aero-line, whereas the remaining parametres were set as constants during the optimisation process. As shown in Fig. 3, ${r_{if}}$ denotes the initial forebody radius of curvature, ${r_{{\rm{max}}}}$ denotes the maximum radius of the fancowl, ${f_{{\rm{max}}}}$ denotes the location of the maximum radius, expressed as a fraction, ${\beta _{{\rm{nac}}}}$ represents the boat-tail angle, ${r_i}$ denotes the intake throat radius, and ${f_i}$ denotes the location of the intake throat, which is expressed as a fraction. The size of the fan and the specific engine limit some of the parametres to being constant. Specifically, when designing an engine nacelle, the ${r_{{\rm{fan}},}},{r_{{\rm{hi}}}},{r_{{\rm{te}},}}\;{l_{{\rm{nac}}}}$ and ${l_{{\rm{int}}}}$ should match the size of the engine because the aerodynamic shape optimisation of the nacelle is usually performed after the design of the engine. The non-asymmetry of the nacelle was achieved based on three different aero-lines at ${\rm{\Phi }} = {0^ \circ },{\rm{\Phi }} = {90^ \circ }$ , and ${\rm{\Phi }} = {180^ \circ }$ , as shown in Fig. 4, with ${\rm{\Phi }} = {0^ \circ }$ being the top aero-line of the nacelle. It is worth mentioning that the aero-lines at ${\rm{\Phi }} = {90^ \circ }$ and ${\rm{\Phi }} = {270^ \circ }$ were the same, owing to the left-right symmetry.

Figure 3. Design variables of aero-lines.

Figure 4. Schematic of the boundary conditions.

2.2 Computational method

In this study, the numerical simulations were performed using an in-house CFD solver. The finite-volume method was applied to discretise the governing equations. The implicit Gauss–Seidel method was employed for the time discretisation. Scalar dissipation coupled with the central difference was used to discretise the spatial derivative. The shear stress transport (SST) k-w turbulence model was adopted to close the Reynolds-averaged Navier–Stokes (RANS) equation. Ideal compressible air was assumed for all the simulations. The freestream condition was specified using a pressure farfield boundary condition by setting the static pressure ${P_{{\rm{st}}}}$ , static temperature ${T_{{\rm{st}}}}$ , and Mach number ${M_\infty }$ . The takeoff and cruise flight conditions were considered in this study. The exit-pressure boundary condition was applied at the duct ends, extended by 1.8 fan radii downstream of the fan location. It’s understood that fan has a significant effect on the flow around with the inlet, especially under the high angle-of-attack. Usually, the presence of fan will improve the performance of nacelle inlet under the large angle. Hence, the current optimisation framework of short nacelle will provide a conservative design. In the future, a fan model will be implemented to improve the optimisation framework for short nacelle. The inlet conditions for the bypass ducts were determined by using ${P_{\rm{t}}}$ and ${T_{\rm{t}}}$ .

To verify the accuracy and validity of the in-house solver, two cases were used as test cases: one was the ONERA M6 wing [Reference Diskin, Anderson and Pandya35] and the other was the NACA-1-Series inlets [Reference Re and Abeyounis36]. A description of the ONERA-M6 case is presented in Table 1. The surface pressure coefficients calculated using the in-house solver are shown in Fig. 5 and were compared with the data calculated using the FUN3D solver [Reference Anderson and Bonhaus37–Reference Nielsen39] and the experimental data. Here, $\eta $ refers to the span station (0 = root; 1 = tip), and the x-location was non-dimensionalised by the local chord. The skin friction coefficient calculated using the in-house solver is shown in Fig. 6 for comparison with the data calculated using the FUN3D solver. Evidently, in both figures, the results calculated using the in-house solver show good agreement with the FUN3D results. Moreover, in Fig. 5, the in-house solver performs even better than the FUN3D solver when $\eta = 0.44$ and $\eta = 0.65$ , and the in-house pressure coefficient distribution is smoother when $x/c = 0.56$ and $x/c = 0.48$ than the FUN3D pressure coefficient distribution.

Table 1. ONERA M6 case description

Figure 5. Comparisons of surface-pressure coefficients along the wing sections.

An axisymmetric NACA 1-85-43.9, with a contraction ratio of 1.250, was adopted as the second test case. Examples of the computed wall-pressure coefficient ${C_p}$ distribution along the chord normalised abscissa $x/c$ are shown in Fig. 7 for four Mach numbers and mass flow capture ratios (MFCRs) [Reference Magrini and Benini40]. Evidently, in all four different boundary conditions, the results calculated using the in-house solver show good agreement with the experimental results. The above results show that the in-house solver is suitable for the subsequent optimisation procedures.

To mesh the 3D non-axisymmetric nacelle, ANSYS ICEM was used to create a structured mesh for the computational domain. The first cell height was set to achieve $y_1^ + \approx 1$ with a stretching ratio of 1.2. The mesh details are shown in Fig. 8. As shown in Fig. 9, there are five main parametres for controlling grid variation. The first parametre $a$ indicates the number of nodes in the boundary layer normal to the wall. The second parametre, $b$ , indicates the average spacing on the fancowl surface in the axial direction, which also influences the inner surface grid distribution owing to the structured block. The number of nodes in the circumferential and radial directions was determined using the third parametre $c$ and fourth parametre $d$ , respectively. The number of nodes from the nacelle to the farfield was determined using the fifth parametre, $e$ . The other parametres remained unchanged. Table 2 summarises these parametres. Tables 3 and 4 summarises the boundary conditions used in this study. A mesh-independence study was performed for three different mesh sizes (3.9 M, 8.67 M, and 15.8 M), as shown in Fig. 10. It can be observed that the drag coefficient at cruise condition and pressure recovery at takeoff condition reached the grid convergence at a mesh size of 8.67 M. Therefore, this mesh size was adopted for the following high-fidelity CFD simulation, whereas the mesh size of 3.9 M was adopted for the low-fidelity CFD simulation.

Figure 6. Comparisons of skin-friction coefficients along the wing sections with results from FUN3D.

Figure 7. Comparisons of wall-pressure coefficients on the inlet forebody exterior and lip over a range of mass-flow ratios.

Figure 8. Mesh details of the nacelle surface and symmetry plane.

Figure 9. Illustration of grid independence varying parametres.

2.3 Multi-fidelity analysis

In engineering optimisation problems, high-fidelity samples are often difficult to obtain, whereas low-fidelity samples are easier to obtain. With the help of low-fidelity samples, the prediction accuracy of the surrogate model was improved. To make good use of the low-fidelity samples, multi-fidelity analysis was performed, and a multi-fidelity surrogate model was constructed. In this study, cokriging was adopted as the metamodel in the optimisation process. The mathematical formula for cokriging is detailed below.

Table 2. Parametres in grid independence

Table 3. Cruise boundary conditions

Table 4. Takeoff boundary conditions

Figure 10. Grid-independence study results showing the cruse and takeoff condition.

The samples used in cokriging are expressed using Equation (5), which consists of low- and high-fidelity samples.

(5) \begin{align}{\bf{X}} = \left( {\begin{array}{l}{{{\bf{X}}_{{\rm{lo}}}}}\\[5pt] {{{\bf{X}}_{{\rm{hi}}}}}\end{array}} \right) = {\!\left( {{\bf{x}}_{{\rm{lo}}}^{\!\left( 1 \right)} \cdots {\bf{x}}_{{\rm{lo}}}^{\!\left( {{n_{{\rm{lo}}}}} \right)}{\bf{x}}_{{\rm{hi}}}^{\!\left( 1 \right)} \cdots {\bf{x}}_{{\rm{hi}}}^{\!\left( {{n_{{\rm{hi}}}}} \right)}} \right)^{\rm{T}}}\end{align}

where ${{\bf{X}}_{{\rm{lo}}}}$ denotes the low-fidelity sample locations, ${{\bf{X}}_{{\rm{hi}}}}$ denotes the high-fidelity sample locations, and the high-fidelity sample locations coincided with a subset of the low-fidelity sample locations ${{\bf{X}}_{{\rm{hi}}}} \subset {{\bf{X}}_{{\rm{lo}}}}$ . In addition, ${n_{{\rm{lo}}}}$ represents the number of low-fidelity samples and ${n_{{\rm{hi}}}}$ represents the number of high-fidelity samples. The value at ${\bf{X}}$ was treated as if it were a realisation of a stochastic process.

(6) \begin{align}{\bf{Y}} = \left( {\begin{array}{l}{{{\bf{Y}}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}}} \right)}\\[5pt] {{{\bf{Y}}_{{\rm{hi}}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}}} \right)}\end{array}} \right) = {\!\left( {{{\bf{Y}}_{{\rm{lo}}}}\!\left( {{\bf{X}}_{{\rm{lo}}}^{\!\left( 1 \right)}} \right) \cdots {{\bf{Y}}_{{\rm{lo}}}}\!\left( {{\bf{X}}_{{\rm{lo}}}^{\!\left( {{n_{{\rm{lo}}}}} \right)}} \right){{\bf{Y}}_{{\rm{hi}}}}\!\left( {{\bf{X}}_{{\rm{hi}}}^{\!\left( 1 \right)}} \right) \cdots {{\bf{Y}}_{{\rm{hi}}}}\!\left( {{\bf{X}}_{{\rm{hi}}}^{\!\left( {{n_{{\rm{hi}}}}} \right)}} \right)} \right)^{\rm{T}}}\end{align}

The local features of the low-fidelity and high-fidelity simulations were represented as the Gaussian process ${{\rm{Z}}_{{\rm{lo}}}}\!\left( . \right)$ . and ${{\rm{Z}}_{{\rm{hi}}}}\!\left( . \right)$ . Based on the auto-regressive model, the high-fidelity simulation was approximated as a low-fidelity simulation multiplied by a scaling factor $\rho $ plus a Gaussian process ${{\rm{Z}}_{\rm{d}}}\!\left( . \right)$ , which represented the difference between $\rho {{\rm{Z}}_{{\rm{lo}}}}\!\left( . \right)$ and ${{\rm{Z}}_{{\rm{hi}}}}\!\left( . \right)$ . The approximation is expressed as

(7) \begin{align}{{\rm{Z}}_{{\rm{hi}}}}\!\left( {\bf{x}} \right) = \rho {{\rm{Z}}_{{\rm{lo}}}}\!\left( {\bf{x}} \right) + {{\rm{Z}}_{\rm{d}}}\!\left( {\bf{x}} \right)\end{align}

The covariance matrix is constructed as

(8) \begin{align}{\bf{C}} = \left( {\begin{array}{c@{\quad}c}{\sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)} & {\rho \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)}\\[5pt] {\rho \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{\rm{h}}},{{\bf{X}}_{{\rm{lo}}}}} \right)} & {\rho \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right) + \sigma _{\rm{d}}^2{{\boldsymbol{\Psi }}_{\rm{d}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)}\end{array}} \right)\end{align}

where ${{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)$ denotes a matrix of correlations of the form ${\psi _{{\rm{lo}}}}$ between data ${{\bf{X}}_{{\rm{lo}}}}$ and ${{\bf{X}}_{{\rm{hi}}}}$ . Similarly, ${{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)$ denotes a matrix of correlations of the form ${\psi _{{\rm{lo}}}}$ between data ${{\bf{X}}_{{\rm{lo}}}}$ and ${{\boldsymbol{\Psi }}_{{\rm{hi}}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)$ denotes a matrix of correlations of the form ${\psi _{{\rm{lo}}}}$ between data ${{\bf{X}}_{{\rm{hi}}}}$ . Furthermore, ${{\boldsymbol{\Psi }}_{\rm{d}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)$ denotes a matrix of correlations of the form ${\psi _{\rm{d}}}$ between data ${{\bf{X}}_{{\rm{hi}}}}$ . Here, the squared exponential kernel was selected to describe the correlation between random variables, where ${p_j}$ was fixed at ${p_{\!\left( {1,2, \ldots ,k} \right)}} = 2$ . The specific expression is as follows:

(9) \begin{align}cor\!\left[ {Y\!\left( {{{\bf{x}}^{\!\left( i \right)}}} \right),Y\!\left( {{{\bf{x}}^{\!\left( l \right)}}} \right)} \right] = {\rm{exp}}\!\left( { - \mathop \sum \limits_{j = 1}^k {\theta _j}{{\left| {x_j^{\!\left( i \right)} - x_j^{\!\left( l \right)}} \right|}^{{p_j}}}} \right)\end{align}

where $k$ is the number of design variables, ${\theta _j}$ represents the correlation coefficient of $j$ -th design variable, and $x_j^{\!\left( i \right)}$ is the $j$ -th design variable of $i$ -th sample.

The basis of cokriging is that our prediction of a new expensive sample is assumed to be consistent with the observed data and the maximum-likelihood-estimations for the model parametres. Therefore, we augmented the observed data with a predicted value and maximised the likelihood of this augmented data by varying our prediction while keeping the model parametres fixed. The augmented dataset is defined as follows:

(10) \begin{align}\tilde{\textbf{X}} = \left( {\begin{array}{c}{{{\bf{X}}_{{\rm{lo}}}}}\\[5pt] {{{\bf{X}}_{{\rm{hi}}}}}\\[5pt] {\bf{X}}\end{array}} \right) = {\!\left( {{\bf{x}}_{{\rm{lo}}}^{\!\left( 1 \right)} \cdots {\bf{x}}_{{\rm{lo}}}^{\!\left( {{n_{{\rm{lo}}}}} \right)}{\bf{x}}_{{\rm{hi}}}^{\!\left( 1 \right)} \cdots {\bf{x}}_{{\rm{hi}}}^{\!\left( {{n_{{\rm{hi}}}}} \right)}{\bf{x}}} \right)^{\rm{T}}}\end{align}

In addition, the values at locations with a predicted value are defined as

(11) \begin{align}\tilde{\textbf{y}} = {\!\left( {{\bf{y}}_{{\rm{lo}}}^{\rm{T}}{\bf{y}}_{{\rm{hi}}}^{\rm{T}}{{{\hat{\textbf{y}}}}_{{\rm{hi}}\left( {\bf{x}} \right)}}} \right)^{\rm{T}}}\end{align}

The new covariance matrix $\widetilde{\boldsymbol{C}}$ can be expressed as

(12) \begin{align}\tilde{\boldsymbol{C}} = \left( {\begin{array}{c@{\quad}c@{\quad}c}{\hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)} & {\rho \hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)} & {\rho \hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{\bf{x}}} \right)}\\[5pt] {\rho \hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{hh}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)} & {\rho _{{\rm{lo}}}^2\hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right) + \hat \sigma _{\rm{d}}^2{{\boldsymbol{\Psi }}_{\rm{d}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)} & {\!\left( {{\rho ^2}\hat \sigma _{{\rm{lo}}}^2 + \hat \sigma _{\rm{d}}^2} \right){{\boldsymbol{\Psi }}_{\rm{d}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{\bf{x}}} \right)}\\[5pt] {\rho \hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}{{\left( {{{\bf{X}}_{{\rm{lo}}}},{\bf{x}}} \right)}^{\rm{T}}}} & {\!\left( {{\rho ^2}\hat \sigma _{{\rm{lo}}}^2 + \hat \sigma _{\rm{d}}^2} \right){{\boldsymbol{\Psi }}_{\rm{d}}}{{\left( {{{\bf{X}}_{{\rm{hi}}}},{\bf{x}}} \right)}^{\rm{T}}}} & {{\rho ^2}\hat \sigma _{{\rm{lo}}}^2 + \hat \sigma _{\rm{d}}^2}\end{array}} \right)\end{align}

Based on the Gaussian process theory, we obtain

(13) \begin{align}{\hat{\textrm{y}}_{{\rm{hi}}}}\!\left( {\bf{x}} \right) = \hat \mu + {{\bf{c}}^{\rm{T}}}{{\bf{c}}^{ - 1}}\!\left( {{\bf{y}} - \textbf{1}\hat \mu } \right)\end{align}

where

(14) \begin{align}{\bf{c}} = \left( {\begin{array}{c}{\rho \hat \sigma _{{\rm{lo}}}^2{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{\bf{x}}} \right)}\\[5pt] {\!\left( {{\rho ^2}\hat \sigma _{{\rm{lo}}}^2 + \hat \sigma _{\rm{d}}^2} \right){{\boldsymbol{\Psi }}_{\rm{d}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}},{\bf{x}}} \right)}\end{array}} \right)\end{align}

(15) \begin{align}\hat \mu = \frac{{{\textbf{1}^{\rm{T}}}{{\rm{C}}^{ - 1}}{\bf{y}}}}{{{\textbf{1}^{\rm{T}}}{{\rm{C}}^{ - 1}}\textbf{1}}}\end{align}

Here $\textbf{1}$ is an $n \times 1$ column vectors of ones. To obtain ${\mu _{{\rm{lo}}}}$ , $\hat \sigma _{{\rm{lo}}}^2$ , and ${\theta _{{\rm{lo}}}}$ for the low-fidelity Gaussian process ${{\rm{Z}}_{{\rm{lo}}}}\!\left( {\bf{x}} \right)$ , we can maximise the ln-likelihood as

(16) \begin{align} - \frac{n}{2}{\rm{ln}}\!\left( {\sigma _{{\rm{lo}}}^2} \right) - \frac{1}{2}{\rm{ln}}{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right) - \frac{{{{\!\left( {{{\bf{y}}_{{\rm{lo}}}} - \textbf{1}{\mu _{{\rm{lo}}}}} \right)}^{\rm{T}}}{{\boldsymbol{\Psi }}_{{\rm{lo}}}}{{\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)}^{ - 1}}\!\left( {{{\bf{y}}_{{\rm{lo}}}} - \textbf{1}{\mu _{{\rm{lo}}}}} \right)}}{{2\sigma _{{\rm{lo}}}^2}}\end{align}

The MLEs of

(17) \begin{align}{\hat \mu _{{\rm{lo}}}} = \frac{{{\textbf{1}^{\rm{T}}}{{\boldsymbol{\Psi }}_{{\rm{lo}}}}{{\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)}^{ - 1}}{{\bf{y}}_{{\rm{lo}}}}}}{{{\textbf{1}^{\rm{T}}}{{\boldsymbol{\Psi }}_{{\rm{lo}}}}{{\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)}^{ - 1}}\textbf{1}}}\end{align}

(18) \begin{align}\hat \sigma _{{\rm{lo}}}^2 = \frac{{{{\!\left( {{{\bf{y}}_{{\rm{lo}}}} - \textbf{1}{{\widehat \mu }_{{\rm{lo}}}}} \right)}^{\rm{T}}}{{\boldsymbol{\Psi }}_{{\rm{lo}}}}{{\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)}^{ - 1}}\!\left( {{{\bf{y}}_{{\rm{lo}}}} - \textbf{1}{{\widehat \mu }_{{\rm{lo}}}}} \right)}}{{{n_{{\rm{lo}}}}}}\end{align}

and ${\widehat{\boldsymbol\theta}_{{\rm{lo}}}}$ could be obtained by maximising the following equation after substituting Equations (17) and (18) in Equation (16):

(19) \begin{align} - \frac{{{n_{{\rm{lo}}}}}}{2}{\rm{ln}}\!\left( {\hat \sigma _{{\rm{lo}}}^2} \right) - \frac{1}{2}{\rm{ln}}\!\left| {{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)} \right|\end{align}

To estimate ${\mu _{\rm{d}}},\sigma _{\rm{d}}^2,{\boldsymbol\theta _{\rm{d}}}$ , and $\rho $ , we defined the following:

(20) \begin{align}{\bf{d}} = {{\bf{y}}_{{\rm{hi}}}} - \rho {{\bf{y}}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{hi}}}}} \right)\end{align}

Similarly, ${\mu _{\rm{d}}}$ and $\sigma _{\rm{d}}^2$ were obtained by maximising the ln-likelihood of the expensive sample as

(21) \begin{align} - \frac{n}{2}{\rm{ln}}\!\left( {\sigma _{\rm{d}}^2} \right) - \frac{1}{2}{\rm{ln}}{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right) - \frac{{{{\!\left( {{\bf{d}} - \textbf{1}{\mu _{\rm{d}}}} \right)}^{\bf{T}}}{{\boldsymbol{\Psi }}_{\rm{d}}}{{\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)}^{ - 1}}\!\left( {{\bf{d}} - \textbf{1}{\mu _{\rm{d}}}} \right)}}{{2\sigma _{\rm{d}}^2}}\end{align}

Equation (21) yields the MLEs of the following as:

(22) \begin{align}{\hat \mu _{\rm{d}}} = \frac{{{\textbf{1}^{\rm{T}}}{{\boldsymbol{\Psi }}_{\rm{d}}}{{\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)}^{ - 1}}{\bf{d}}}}{{{\textbf{1}^{\rm{T}}}{{\boldsymbol{\Psi }}_{\rm{d}}}{{\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)}^{ - 1}}\textbf{1}}}\end{align}

(23) \begin{align}\hat \sigma _{\rm{d}}^2 = \frac{{{{\!\left( {{\bf{d}} - \textbf{1}{{\widehat \mu }_{\rm{d}}}} \right)}^{\rm{T}}}{{\boldsymbol{\Psi }}_{\rm{d}}}{{\left( {{{\bf{X}}_{{\rm{hi}}}},{{\bf{X}}_{{\rm{hi}}}}} \right)}^{ - 1}}\!\left( {{\bf{d}} - \textbf{1}{{\widehat \mu }_{\rm{d}}}} \right)}}{{{n_{{\rm{hi}}}}}}\end{align}

Furthermore, ${\widehat{\boldsymbol\theta}_{\rm{d}}}$ could be obtained by maximising the following equation after substituting Equations (22) and (23) in Equation (21)

(24) \begin{align} - \frac{{{n_{{\rm{hi}}}}}}{2}{\rm{ln}}\!\left( {\hat \sigma _{\rm{d}}^2} \right) - \frac{1}{2}{\rm{ln}}\!\left| {{{\boldsymbol{\Psi }}_{{\rm{lo}}}}\!\left( {{{\bf{X}}_{{\rm{lo}}}},{{\bf{X}}_{{\rm{lo}}}}} \right)} \right|\end{align}

With the above derivation, it is possible to use cokriging to make effective predictions. For more details about derivation and implementation, resources could be found in Refs [Reference Sobester, Forrester and Keane41, Reference Bouhlel, Hwang and Bartoli42].

2.4 Infill strategy

To improve the prediction accuracy of the surrogate models, the low- and high-fidelity samples must be updated after each iteration until the termination condition is achieved. Infilling low- and high-fidelity samples were generated from the current Pareto front to refine the surrogate models. Specifically, to update the high-fidelity samples, the current ${X_{{\rm{Parete}}}}$ was first clustered using the K-means clustering method. The Pareto fitness function (MG) was then calculated for each point in ${X_{{\rm{Parete}}}}$ [Reference Shi, Long and Baoyin27, Reference Kitayama, Srirat, Arakawa and Yamazaki43]:

(25) \begin{align}{G_i} = 1 - \mathop {{\rm{max}}}\limits_{j \ne i} \!\left[ {{\rm{min}}\!\left( {f_{s1}^i - f_{s1}^j, \ldots ,f_{sk}^i - f_{sk}^j} \right)} \right],i = 1,2, \cdots ,m\end{align}

where ${G_i}$ denotes the Pareto fitness at the $i$ -th sampling point, and $f_{sk}^i$ is the normalised $k$ -th objective function at the $i$ -th sampling point, given by

(26) \begin{align}f_{sk}^i = \frac{{{f_{k,i}} - {f_{k,{\rm{min}}}}}}{{{f_{k,{\rm{max}}}} - {f_{k,{\rm{min}}}}}}\end{align}

where ${f_{k,i}}$ is the $k$ -th objective function at the $i$ -th sampling point, and ${f_{k,{\rm{max}}}}$ and ${f_{k,{\rm{min}}}}$ denote the maximum and minimum values of the $k$ -th objective function among all sampling points, respectively. When the Pareto fitness ${G_i}$ is greater than one, the sampling point is a Pareto-optimal solution. A modified Pareto fitness function was calculated to construct a function that generates local minima around a set of Pareto-optimal solutions.

(27) \begin{align}M{G_i} = \left\{ {\begin{array}{c@{\quad}c}{\!\left( {1 - {G_i}} \right) + 2 \times {G_{{\rm{max}}}}} & {{G_i} \lt 1}\\[5pt] {{G_i}} & {{G_i} \geq 1}\end{array}\;\;\; i = 1,2, \cdots ,m} \right.\end{align}

A lower MG value indicates that the sample is closer to the Pareto front, and the sample point with the minimum MG value in each cluster was selected as the high-fidelity infill sample point, as follows:

(28) \begin{align}x_{{\rm{infill}}}^{\!\left( k \right)} = {\rm{arg\;min}}\;MG\!\left( x \right){\rm{with}}\;x \in {X_k},k = 1,2, \cdots ,{k_c}\end{align}

where ${X_k}$ is the sample set of the $k$ -th cluster and $\chi _{{\rm{infill}}}^{\!\left( k \right)}$ are the selected high-fidelity infill samples. Specifically, all Pareto fronts were selected to update the low-fidelity samples.

With the high-fidelity and low-fidelity infill strategies described above, both exploration and exploitation could be performed.

2.5 Genetic algorithm

The NSGA-II, originally proposed by Deb, was adopted in this study because of its proven capabilities for global optimisation in flow aerodynamic applications. NSGA-II generates offspring using specific types of crossover and mutation and then selects the next generation according to a non-dominated sorting and crowding distance comparison. The initial seed was set at 36, with 36 individuals in the following generation and a total of 5,000 generations. The mutation probability was set at 0.1 and the crossover rate was set at 0.9. The selection operator was set as a tournament, the mutation operator was set as a polynomial mutation, and the crossover operator was set as the simulated binary crossover.

2.6 Validation and verification

Two test functions were used to validate the effectiveness of the MO-ACK framework. To demonstrate the role of multi-fidelity, a multi-objective optimisation framework based on kriging (MO-AK) was also tested. Comprehensive evaluation indicators were introduced to evaluate the performance of the multi-objective optimisation process, such as the hypervolume (HV) indicator [Reference Zitzler and Thiele44] and inverted generational distance (IGD) indicator [Reference Czyz and Jaszkiewicz45]. In this study, the HV indicator was calculated through optimisation iterations to measure the convergence and diversity of the Pareto front solutions. A higher value of HV indicates better capability for exploring the Pareto front. As shown in Fig. 11, the two main factors in the HV indicator are the reference point and the non-dominated solution. When non-dominated solutions are obtained, a reference point must be validated. A schematic of a 2D HV indicator is shown in Fig. 11. The HV indicator used in higher dimensional problems follows this analogy.

2.6.1 FON test function optimisation

The definition of the FON two-objective test function, composed of high-fidelity and low-fidelity functions, is as follows [Reference Huband, Hingston, Barone and While46]:

(29) \begin{align}\eqalign{ {\rm{Minimise}}:\;\;{f_1}\!\left( {{x_1},{x_2}, \cdots ,{x_5}} \right) = 1 - {\rm{exp}}\!\left( { - \mathop \sum \limits_{i = 1}^5 {{\!\left( {{x_i} - {1 \over {\sqrt 5 }}} \right)}^2}} \right) \cr {\rm{Minimise}}:\;\;{f_2}\!\left( {{x_1},{x_2}, \cdots ,{x_5}} \right) = 1 - {\rm{exp}}\!\left( { - \mathop \sum \limits_{i = 1}^5 {{\!\left( {{x_i} + {1 \over {\sqrt 5 }}} \right)}^2}} \right) \cr {\rm{Minimise}}:\;\;{f_{11}}\!\left( {{x_1},{x_2}, \cdots ,{x_5}} \right) = 1 - {\rm{exp}}\!\left( { - \mathop \sum \limits_{i = 1}^5 {{\!\left( {0.5{x_i} - 0.1 - {1 \over {\sqrt 5 }}} \right)}^2}} \right) \cr {\rm{Minimise}}:\;\;{f_{12}}\!\left( {{x_1},{x_2}, \cdots ,{x_5}} \right) = 1 - {\rm{exp}}\!\left( { - \mathop \sum \limits_{i = 1}^5 {{\!\left( {0.6{x_i} + 0.1 + {1 \over {\sqrt 5 }}} \right)}^2}} \right) \cr} \end{align}

where ${f_1}$ and ${f_2}$ denote the high-fidelity functions, and ${f_{11}}$ and ${f_{12}}$ denote the low-fidelity functions. The design space of the test function is ${x_1},{x_2}, \cdots ,{x_5} \in \!\left[ { - 4,4} \right]$ .

Figure 11. Schematic of a 2D HV indicator.

The optimisation routine was initialised with LHS to create 20 initial high-fidelity samples and 40 initial low-fidelity samples. During each iteration, the high-fidelity samples were updated by adding five samples using the high-fidelity infill strategy. Meanwhile, the low-fidelity samples were updated by adding 10 samples, including the samples added to the high-fidelity samples using the low-fidelity infill strategy. The optimisation routine ended when the prescribed 24 iterations were met. This was consecutively tested in eight trials. The results of the eight trials are shown in Fig. 12. Moreover, 500,000 samples were generated by LHS, and the HV was calculated as the reference, as shown in Fig. 12.

As observed from Fig. 12, the MO-ACK model outperformed the MO-AK model in terms of achieving higher diversity because its HV was larger after five iterations in all trials. Eight trials were initialised with random LHS; therefore, the initial HV values were different. Both MO-ACK and MO-AK performed better than the reference. It can also be observed from the boxplots of the HV values from the eight trials that the MO-ACK has a larger mean HV value than the MO-AK. The NSGA-II was also included for comparison. The population size was set to 28 and the maximum generation was set to 20. As shown in Fig. 13, the NSGA-II exhibited the worst performance.

In comparison to the MO-AK model, we can conclude that the MO-ACK model can effectively explore the Pareto front of a relatively low-dimension optimisation problem.

2.6.2 DTLZ2 test function optimisation

The definition of the DTLZ2 three-objective test functions, composed of high-fidelity and low-fidelity functions, is as follows [Reference Huband, Hingston, Barone and While46]:

Figure 12. HV results of the FON test function.

Figure 13. Boxplots of the HV values of the FON test function from the eight trials.

(30) \begin{align}\eqalign{ {\rm{Minimise}}:\;\;{f_1}\!\left( {{x_1},{x_2}, \cdots ,{x_{18}}} \right) = \left( {1 + g} \right){\rm{cos}}\!\left( {0.5{x_1}\pi } \right){\rm{cos}}\!\left( {0.5{x_2}\pi } \right) \cr {\rm{Minimise}}:\;\;{f_2}\!\left( {{x_1},{x_2}, \cdots ,{x_{18}}} \right) = \left( {1 + g} \right){\rm{cos}}\!\left( {0.5{x_1}\pi } \right){\rm{sin}}\!\left( {0.5{x_2}\pi } \right) \cr {\rm{Minimise}}:\;\;{f_3}\!\left( {{x_1},{x_2}, \cdots ,{x_{18}}} \right) = \left( {1 + g} \right){\rm{sin}}\!\left( {0.5{x_1}\pi } \right) \cr g = \mathop \sum \limits_{i = 3}^{18} {\!\left( {{x_i} - 0.5} \right)^2} \cr {f_{11}}\!\left( {{x_1},{x_2}, \cdots ,{x_{18}}} \right) = \left( {0.5 + 1.5g} \right){\rm{cos}}\!\left( {0.6{x_1}\pi } \right){\rm{cos}}\!\left( {0.4{x_2}\pi } \right), \cr {f_{12}}\!\left( {{x_1},{x_2}, \cdots ,{x_{18}}} \right) = \left( {0.3 + 1.2g} \right){\rm{cos}}\!\left( {0.4{x_1}\pi } \right){\rm{sin}}\!\left( {0.6{x_2}\pi } \right), \cr {f_{13}}\!\left( {{x_1},{x_2}, \cdots ,{x_{18}}} \right) = \left( {0.4 + 1.3g} \right){\rm{sin}}\!\left( {0.5{x_1}\pi } \right) \cr {g_c} = \mathop \sum \limits_{i = 3}^{18} {\!\left( {0.8{x_i} - 0.3} \right)^2} \cr} \end{align}

where ${f_1}$ , ${f_2}$ , and ${f_3}$ denote the high-fidelity functions, and ${f_{11}}$ , ${f_{12}}$ , and ${f_{13}}$ denote the low-fidelity functions. Furthermore, the design space of the test function is ${x_1},{x_2}, \cdots ,{x_{18}} \in \!\left[ {0,1} \right]$ .

The optimisation routine was initialised with LHS to create 120 initial high-fidelity samples and 240 initial low-fidelity samples. During each iteration, the high-fidelity samples were updated by adding 18 samples using the high-fidelity infill strategy. Meanwhile, the low-fidelity samples were updated by adding 36 samples, including the samples added to the high-fidelity samples using the low-fidelity infill strategy. The optimisation routine ended when the four prescribed iterations were met. This was consecutively tested in eight trials. The results of the eight trials are shown in Fig. 14. Moreover, 500,000 samples were generated by LHS, and the HV was calculated as the reference, as shown in Fig. 14.

Figure 14. HV results of the DTLZ2 test function.

As observed from Fig. 14, the MO-ACK model outperformed the MO-AK model in terms of achieving higher diversity because its HV was larger after five iterations in all trials. Eight trials were initialised with random LHS; therefore, the initial HV values were different. All eight trials of the MO-ACK obtained higher HV values than the reference, whereas all eight trials of the MO-AK obtained lower HV values. It can also be observed from the boxplots of the HV values from the eight trials that the MO-ACK has a larger mean HV value than the MO-AK. The NSGA-II was also included for comparison. The population size was set to 32, and the maximum generation was set to 6. As shown in Fig. 15, the NSGA-II exhibited the worst performance.

Figure 15. Boxplots of the HV values of the DTLZ2 test function from the eight trials.

We may infer that the MO-ACK model can more effectively explore the Pareto front of a high-dimensional optimisation problem than the MO-AK model.

3.1 Results under cruise condition

Designs A1 and B1 exhibited significant reductions in terms of ${C_{D - {\rm{cruise}}}}$ compared to the configuration of the base sample. While design A1 exhibited a reduction of 23.6% while design B1 exhibited a reduction of 22.2%. As shown in Fig. 19, the large variation in the nacelle drag characteristics was caused by a fundamental change in the associated nacelle transonic flow aerodynamics. The pressure distributions of all designs are shown in Fig. 20. As shown in Fig. 21, the isentropic Mach number distributions of the base sample configuration on all four aero-lines were almost the same because of their axisymmetric characteristics. The location of the peak isentropic Mach number at the top aero-line ( $\psi = {0^ \circ }$ ) was significantly different for the three selected designs and varied by 2.474. Design A1, which had the highest $IP{R_{{\rm{takeoff}}}}$ found throughout the optimisation process, had an initial acceleration to terminate with a strong shockwave at $X = 0.14$ . Design A1 had a well-defined shockwave at $X = 0.14$ , whereas the configuration of the base sample had a similar shockwave at $X = 1.697$ . Compared to the other two configurations, design B1 had a shockwave in the shortest position at $X = - 0.777$ . Unlike the other two, which experienced accelerations at the start, design B1 experienced a stable process at the beginning. For the azimuthal section $\psi = {45^ \circ }$ , design A1 had a similar benign transonic flow structure compared to design B1. For the azimuthal section $\psi = {90^ \circ }$ , design A1 had a well-defined shock wave at $X = - 0.325$ . Design B1 and the configuration of the base sample had similar flow structures except for the intensity of the associated transonic flow aerodynamics. For instance, the configuration of the base sample had an isentropic Mach number of 1.207, whereas design B1 had a Mach number of 1.177. For the azimuthal section $\psi = {180^ \circ }$ , design B1 still had flow structures similar to those of the base sample. Design B1 had a double shockwave structure with the first shockwave located at $X = - 0.370$ and the second shockwave at $X = 1.510$ . The base sample had a double shockwave structure with the first shockwave located at $X = - 0.322$ and the second shockwave at $X = 1.870$ . Design A1 had an initial acceleration and terminated with a shockwave at $X = 0.245$ .

The design with the largest radius of curvature at the azimuthal section $\psi = {0^ \circ }$ and $\psi = {90^ \circ }$ has the highest peaks in isentropic Mach number, which is consistent in Figs 18 and 21. For the azimuthal section $\psi = {180^ \circ }$ , design A1 and design B1 have similar geometric configurations, but there is a significant difference at isentropic Mach number distributions, which indicates the sensitivity of the flow field to geometry and the strong nonlinearity of the transonic flow aerodynamics.

Figure 19. Comparison of the isentropic Mach number contours for the selected designs.

Figure 20. Comparison of the pressure contours for the selected designs.

Figure 21. Comparisons of the isentropic Mach number for the selected designs at different azimuthal sections.

3.2 Results under takeoff condition

Designs A1 and B1 had significant improvements in terms of $IP{R_{{\rm{takeoff}}}}$ compared to the configuration of the base sample. While design A1 had a reduction of 1.63%, design B1 had 1.59%. In Fig. 22, the total pressure distribution on the fanface under the takeoff condition is clearly illustrated. For the base sample, there was a distinct low-pressure area at the lower part of the fanface. This may occur because of the flow separation induced by the shockwave or adverse pressure gradient. This is an undesirable situation during takeoff. If internal separation occurs, the resulting fan entrance distortion could be critical to causing excessively high fan-blade stresses and may even cause core-compressor stalls. Designs A1 and B1 had significant improvements in $IP{R_{{\rm{takeoff}}}}$ compared to the base sample. As shown in Fig. 22, the low-pressure area at the lower part of the fanface of designs A1 and B1 became significantly smaller. This also indicates that the $IP{R_{{\rm{takeoff}}}}$ was higher. As shown in Fig. 18, design A1 and design B1 are similar in all azimuthal sections and have a smaller thickness compared to the base sample. Combined with the flow field analysis in Fig. 22, these designs allow for a lower total pressure loss when the airflow reaches the fan surface. Compared to design B1, the low-pressure area of design A1 was smaller and more concentrated in the lower part.

Figure 22. Contours of the total pressure at the fanface for the selected designs.

Figure 23. Parallel coordinates plot of optimised designs coloured by ${C_{D - {\rm{cruise}}}}$ .

Figure 24. Parallel coordinates plot of optimised designs coloured by $IP{R_{{\rm{takeoff}}}}$ .

3.3 Influence of the design variables on the nacelle performance

Parallel coordinate plots represent the regions of the design space where the non-dominated solutions are concentrated [Reference Tejero, Robinson, MacManus and Sheaf15]. Here the ten samples of MO-ACK closest to the origin are selected as in Fig. 1, and the distance criterion formula is defined as follows:

(33) \begin{align}distance = \sqrt {{{\!\left( {{C_{D - {\rm{cruise}}}} - 0.03} \right)}^2} + {{\left( {1 - IP{R_{{\rm{takeoff}}}}} \right)}^2}} \end{align}

Parallel coordinates plot of optimised designs coloured by ${C_{D - {\rm{cruise}}}}$ is shown in Fig. 23. Parallel coordinates plot of optimised designs coloured by $IP{R_{{\rm{takeoff}}}}$ is shown in Fig. 24. The figures show that the range of variation of design variables ${r_{if}}$ and ${r_i}$ is small in all azimuthal sections and converges to the design space boundary. This provides some guidelines for design space considerations for subsequent preliminary designs. The overall trend of the design parametres is consistent in the azimuthal section $\psi = {0^ \circ }$ and azimuthal section $\psi = {180^ \circ }$ , and differs significantly from azimuthal section $\psi = {90^ \circ }$ . In the azimuthal section $\psi = {0^ \circ }$ , design variable ${f_{{\rm{max}}}}$ has the largest range of variation, which is consistent in the azimuthal section $\psi = {180^ \circ }$ , but in the azimuthal section $\psi = {90^ \circ }$ , design variable ${r_{{\rm{max}}}}$ has the largest range of variation. Overall, the range of design variables changed the least in the azimuthal section $\psi = {90^ \circ }$ and the most in the azimuthal section $\psi = {180^ \circ }$ . This is likely because throughout the multi-objective optimisation process, the takeoff condition was taken into account but the crosswind condition was not. This also demonstrates the complexity of multi-objective optimisation of non-axisymmetric nacelles.