2.1 Computational grid and boundary conditions

The Rolls-Royce code PADRAM^{(Reference Milli and Shahpar19)} was used for both generating the geometry and for creating its corresponding computational grid, as shown in Fig. 3. The code is an algebraic-elliptic, multi-block grid generator based on a C-O-H topology. For all three fan blades, the mesh consisted of downstream and upstream H blocks for the stationary zones, upper and lower H blocks for the periodic boundaries, an O-mesh around the blade, and a C-mesh for the downstream splitter. The meshes were created with eleven blocks each; details of the meshes for the three blades is provided in Table 1.

Figure 3.

Representative mesh topology used for all the blades (shown here for Blade C).

A grid independence study was undertaken on all three meshes. For completeness, we provide details of the mesh convergence study undertaken on Blade C at cruise. The number of cells in the spanwise, circumferential and axial direction were varied for convergence in the pressure ratio and efficiency, see Fig. 4.

Figure 4.

Mesh convergence plots for Blade C at cruise. Shown are the objectives of (a) non-dimensional efficiency and (b) non-dimensional pressure ratio. The mesh with the red circular maker was used in our computations.

These are defined by

(1) \begin{equation}PR=\left(\frac{\dot{m}_{\text{byp}} P_{\text{byp}}+ \dot{m}_{\text{core}}P_{\text{core}}}{ \dot{m}_{\text{in}}P_{\text{in}}}\right)\end{equation}

and

(2) \begin{equation}\eta=\frac{\left(PR^{\left(\gamma-1\right)/\gamma}-1\right)}{TR-1}\times100\end{equation}

where the temperature ratio (TR) is given by

(3) \begin{equation}TR=\left(\frac{\dot{m}_{\text{byp}}T_{\text{byp}}+ \dot{m}_{\text{core}}T_{\text{core}}}{\dot{m}_{\text{in}}T_{\text{in}}}\right)\end{equation}

In Equations (1), (2) and (3): P indicates the stagnation pressure; T the stagnation temperature; $\dot{m}$ the mass-flow rate, and the subscripts ‘byp’, ‘core’ and ‘in’ indicate flow quantities at the bypass outlet, core outlet and inlet respectively. We use mass-averaged values (radially and circumferentially) of the stagnation temperature and pressure obtained at the three planes indicated above. In terms of boundary conditions, at the inlet we enforce a fixed stagnation pressure, and at the exit to both the core and the bypass we utilise a fixed exit capacity $\Gamma$ boundary condition:

(4) \begin{equation}\Gamma_{\text{core}}=\dot{m}_{\text{core}}\frac{\sqrt{T_{\text{core}}}}{P_{\text{core}}}\end{equation}

(5) \begin{equation}\Gamma_{\text{byp}}=\dot{m}_{\text{byp}}\frac{\sqrt{T_{\text{byp}}}}{P_{\text{byp}}}.\end{equation}

The endwalls (hub, blade, casing and splitter) are characterised by viscous wall boundary conditions with both the hub and the blade having an angular velocity based on the shaft speed.

The mesh convergence study involved 8 different meshes, each with a different number of grid points. The number of points in the grid were varied by increasing (i) axial nodes; (ii) circumferential nodes; (iii) spanwise nodes (including the tip gap) and their distribution. Careful attention was also paid to the nodes where the O-mesh of the blade mated with the upper and lower H-meshes. The results in Fig. 4(a) plot the $\Delta$ Efficiency on the vertical axis (logarithm-base 10 scale), defined as

(6) \begin{equation}\Delta \; \text{Efficiency} = \left| \text{Efficiency}_{\text{finest grid} } - \text{Efficiency}_{\text{current grid} } \right|.\end{equation}

Similarly, in (b), $\Delta$ Pressure ratio is given by

(7) \begin{equation}\Delta \; \text{Pressure ratio} = \left| \text{Pressure ratio}_{\text{finest grid} } - \text{Pressure ratio}_{\text{current grid} } \right| .\end{equation}

The red circles in these plots indicate the mesh used in all our subsequent computations. For this mesh, we have a confidence in the first two decimal places in efficiency and the first four decimal places in pressure ratio. Note that the total number of nodes along the horizontal axis in Fig. 4 varied parameters (i)–(iii) simultaneously.

2.2 Computational flow solver and constitutive equations

We solve the 3D Reynolds-averaged Navier-Stokes (RANS) equations on the aforementioned PADRAM meshes, using the HYDRA flow solver^{(Reference Crumpton, Moinier and Giles6,Reference Lapworth15)} . To close the set of equations when applying a Reynolds-based averaging to the Navier-Stokes, we utilise the one equation Spalart-Allmaras turbulence model^{(Reference Spalart and Allmaras28)}. Further details on the calculation of the inviscid fluxes, the viscous fluxes and convergence can be found in pages 10–21 of Ref. (Reference Östling and Wiberger20). We omit a detailed flow validation here and direct interested readers to Refs (Reference Seshadri, Parks and Shahpar24,Reference Zamboni and Xu29), where experiment vs CFD validations for the PADRAM-HYDRA suite of tools are provided for representative transonic fan blades. Boundary conditions and the inlet Reynolds number for the four design of experiments carried out in this study are given in Table 2.

Table 2

Reynolds number and inlet boundary conditions

Figure 5.

Projecting multi-dimensional samples of the function $\text{exp} \left(x_1 + x_2 + x_3 + x_4 \right)$ over a subspace: (a) ${\textbf{\textit{w}}} = [1,1,1,1]^{T}$ (b) $\boldsymbol{v} = [0.4, -0.3, 0.5, 0.4]^{T}$ .

2.3 Computing active subspaces via a quadratic model

Before we delve into active subspaces, it will be useful to introduce the concept of a sufficient summary plot (following the terminology in Ref. (Reference Cook and Weisberg5)). Consider the 4D function

(8) \begin{equation}f\!\left( {\textbf{\textit{x}}} \right) = \text{exp} \left(x_1 + x_2 + x_3 + x_4 \right),\end{equation}

defined over the hypercube $\mathcal{X} \subset [-1, 1]^{4}$ and where ${\textbf{\textit{x}}} = [x_1, x_2, x_3, x_4]^{T} \in \mathcal{X}$ . We can express this function as

(9) \begin{equation}f\!\left( {\textbf{\textit{x}}} \right) = \text{exp} \left( {\textbf{\textit{w}}}^{T} {\textbf{\textit{x}}} \right), \; \; \; \text{where} \; \; \; {\textbf{\textit{w}}} = [1, 1, 1, 1]^{T}.\end{equation}

The transformation induced by the vector ${\textbf{\textit{w}}}$ projects the 4D values in ${\textbf{\textit{x}}}$ onto the subspace—i.e. linear combination of the parameters—given by ${\textbf{\textit{w}}}^{T} {\textbf{\textit{x}}}$ . Thus, our 4D function is actually 1D over this subspace! Imagine we generate N random samples from $\mathcal{X}$ and project each ${\textbf{\textit{x}}}_{i}$ for $i=1, \ldots, N$ onto this subspace. Then, we evaluate (8) for each ${\textbf{\textit{x}}}_{i}$ and plot the resulting data $f_i = f\!\left( {\textbf{\textit{x}}}_{i} \right)$ . We would arrive at Fig. 5(a), where the horizontal axis of plots the subspace ${\textbf{\textit{w}}}^{T} {\textbf{\textit{x}}}_{i}$ and the vertical axis plots the function values $f_i$ .

Now, typically, we may not know ${\textbf{\textit{w}}}$ or even have access to the functional form of a model. Thus, we have to use techniques based on input-output pairs of the model to approximate this subspace. But what exactly are we looking for? Contrast Fig. 5(a) with (b) where the subspace used for projection is $\boldsymbol{v} = [0.4, -0.3, 0.5, 0.4]^{T}$ . It is readily apparent that although this subspace captures the monotonic trend of the function, any further inference is far more challenging compared to (a), owing to the scatter in the plot.

The above discussion raises the question of how one can discover such dimension reducing subspaces such that if the function admits low dimensional structure of the form of (8), then we can construct its sufficient summary plots. In what follows, we will assume that our functional output comes from CFD and that it does admit such low dimensional structure.

We now consider ${\textbf{\textit{x}}}$ as a set of design parameters and assume ${\textbf{\textit{x}}} \in \mathbb{R}^{25}$ . Without loss of generality, we will assume that ${\textbf{\textit{x}}}$ has been appropriately non-dimensionalised such that $\boldsymbol{-1} \leq {\textbf{\textit{x}}} \leq \boldsymbol{1}$ . Our goal in this subsection is to detail a strategy for reducing the dimensionality of our problem, via

(10) \begin{equation}f({\textbf{\textit{x}}}) \approx g( {\textbf{\textit{W}}}_{1}^T {\textbf{\textit{x}}})\end{equation}

where ${\textbf{\textit{W}}}_{1} \in \mathbb{R}^{d \times m}$ where d is the number of design parameters (here 25) and $m << d$ is the reduced dimension. In Equation (10), g is a suitably chosen m-dimensional function. Ideally, one would like m to be either 1 or 2, facilitating easy visualization. However, this may not necessarily be possible and depends on the relationship between ${\textbf{\textit{x}}}$ and f. So, how exactly do we find ${\textbf{\textit{W}}}_{1}$ ? Constantine^{(Reference Constantine3)} lays the theoretical foundations for approximating ${\textbf{\textit{W}}}_{1}$ and offers numerous practical computational heuristics. In essence, one seeks to estimate the covariance of the gradient, given by

(11) \begin{equation}{\textbf{\textit{K}}}=\int \nabla_{{\textbf{\textit{x}}}}\,f({\textbf{\textit{x}}})\nabla_{{\textbf{\textit{x}}}}\,f({\textbf{\textit{x}}})^{T} \rho ({\textbf{\textit{x}}}) d {\textbf{\textit{x}}}\end{equation}

where ${\textbf{\textit{K}}} \in \mathbb{R}^{d \times d}$ . Here the measure $\rho({\textbf{\textit{x}}})$ is the density of the inputs ${\textbf{\textit{x}}}$ over their parameter space $\mathcal{X}$ . As we can control our sampling strategy, we set $\rho$ to be the uniform measure and generate samples ${\textbf{\textit{x}}}_{i}$ randomly with respect to the uniform distribution over $\mathcal{X} \in [\boldsymbol{-1}, \boldsymbol{1}]$ . Then, one can approximate ${\textbf{\textit{W}}}_{1}$ using the first few eigenvectors associated with ${\textbf{\textit{K}}}$

(12) \begin{equation}{\textbf{\textit{K}}}={\textbf{\textit{W}}} \boldsymbol{ \Lambda} {\textbf{\textit{W}}}^{T}=\left[\begin{array}{c@{\quad}c}{\textbf{\textit{W}}}_{1} & {\textbf{\textit{W}}}_{2}\end{array}\right]\left[\begin{array}{c@{\quad}c}\boldsymbol{ \Lambda_{1} }\\ & \boldsymbol{ \Lambda_{2}}\end{array}\right]\left[\begin{array}{c@{\quad}c}{\textbf{\textit{W}}}_{1} & {\textbf{\textit{W}}}_{2}\end{array}\right]^{T}\end{equation}

Intuitively, one can interpret ${\textbf{\textit{W}}}_{1}$ as the directions along which $\nabla_{{\textbf{\textit{x}}}}\,f$ varies the most, and consequently ${\textbf{\textit{W}}}_{2}$ as the directions along which $\nabla_{{\textbf{\textit{x}}}}\,f$ is (close to) zero.

In the absence of gradients or adjoints (as in Ref. (Reference Seshadri, Shahpar, Constantine, Parks and Adams25)) one can approximate ${\textbf{\textit{W}}}_{1}$ by constructing a global quadratic model in $\mathbb{R}^{m}$

(13) \begin{equation}f\!\left({\textbf{\textit{x}}}_{i}\right) \approx \frac{1}{2} {\textbf{\textit{x}}}^{T} {\textbf{\textit{A}}} {\textbf{\textit{x}}} + {\textbf{\textit{c}}}^{T} {\textbf{\textit{x}}}+e\end{equation}

where the coefficients ${\textbf{\textit{A}}} \in \mathbb{R}^{d \times d}$ , ${\textbf{\textit{c}}} \in \mathbb{R}^{d}$ and $e \in \mathbb{R}$ are estimated via least squares. The number of coefficients M associated with this quadratic model is given by

(14) \begin{equation}N = \left(\begin{array}{c}d+k\\k\end{array}\right),\end{equation}

where k is the highest degree of the polynomial and d once again corresponds to the dimension of the space where the polynomial has to be constructed. As we seek to approximate rather than interpolate, the number of CFD evaluations N required to fit our model should be

(15) \begin{equation}N = \alpha M\end{equation}

where $\alpha \geq 1.5$ is selected to avoid over-fitting. Note that (14) corresponds to a total order index set, where the sum of the composite univariate polynomial orders in the multivariate expansion is less than or equal to k. Such a basis was used in Ref. (Reference Seshadri, Shahpar, Constantine, Parks and Adams25), resulting in a total of 351 coefficients (setting $d=25$ and $k=2$ in (14)) that needed to be estimated. Once obtained, computing gradients of this global quadratic is trivial, resulting in the following estimate of the covariance matrix

(16) \begin{equation}{\textbf{\textit{K}}} \approx \frac{4}{3}{\textbf{\textit{A}}}^{2}+ {\textbf{\textit{c}}} {\textbf{\textit{c}}}^{T}\end{equation}

from which ${\textbf{\textit{W}}}_{1}$ can be readily obtained following the eigendecomposition of this symmetric matrix. Note that in (16), we have assumed a uniform measure

(17) \begin{equation}\rho \left( {\textbf{\textit{x}}} \right) = \frac{1}{2^{d}}.\end{equation}

2.4. Results on the three blades

Design of experiments (DoE) were carried out on the three blades, with the total number of CFD evaluations shown in Table 3. Given the potential for noise in fitting a global quadratic to the data, some oversampling is necessary. Each design of experiment involved generating N uniformly distributed random samples, i.e. Monte Carlo samples, within the hypercube

(18) \begin{equation}\mathcal{X} \subset [-1, 1]^{25} \; \; \; \; \text{where} \; \; \; \; {\textbf{\textit{x}}}_{i} \in \mathcal{X} \; \; \; \text{for} \; i=1, \ldots, N.\end{equation}

Each sample ${\textbf{\textit{x}}}_{i} \in \mathcal{X}$ was then scaled based on the minimum and maximum values of the design parameters within the CFD workflow. Data-sets for these blades can be accessed at https://github.com/psesh/turbodata and were used to determine their active subspace via the global quadratic model.

Table 3

Number of CFD evaluations in DoE processes

The eigenvalues for the covariance matrix ${\textbf{\textit{K}}}$ for these two objectives across for Blade A is shown in Fig. 6. These plots are useful to ascertain the dimensionality of the active subspace; Constantine^{(Reference Constantine3)} suggests that we fix m based on the appearance of large gaps (see page 37 in^{(Reference Constantine3)}) in the eigenvalues. Typically, this is done in conjunction with a parametric bootstrap with 500–1000 replicates. Now, as the number of CFD evaluations in our three cases is far too low for a parametric bootstrap, we adopt a different approach to probe the accuracy of the active subspace. We compare the eigenvalues obtained using all 548 samples with those obtained by selecting 400 random samples without duplicates. In total a 1000 subsamples are selected; the range of obtained eigenvalues is reported in Fig. 6(a) and (b) for the two objectives. In Fig. 6(c) and (d) we increase this to 500 subsamples.

Figure 6.

Eigenvalues of the covariance matrix ${\textbf{\textit{K}}}$ for Blade A for: (a) efficiency with 400 samples; (b) pressure ratio with 400 samples; (c) efficiency with 500 samples; (d) pressure ratio with 500 samples. Shaded regions indicate minimum and maximum eigenvalues from 1000 random repetitions.

A similar study was undertaken for the other two blades; for brevity we only plot the eigenvalues obtained from using all the samples in Fig. 7. Based on these results, for both objectives, for Blades A and C, we can set $m=2$ and for Blade B, $m=1$ . Sufficient summary plots for the three studies are shown in Figs 8 and 9 for the efficiency and pressure ratio objectives respectively. For plotting convenience, we have set $m=2$ for Blade B too.

Figure 7.

Eigenvalues of the covariance matrix ${\textbf{\textit{K}}}$ for Blade B for: (a) efficiency (b) pressure ratio, and for Blade C for (c) efficiency (d) pressure ratio.

Figure 8.

Sufficient summary plots for efficiency values for blades (a) A, (b) B and (c) C, plotted on their active subspaces for efficiency.

Figure 9.

Sufficient summary plots for pressure ratio values for blades (a) A, (b) B and (c) C, plotted on their active subspaces for pressure ratio.

In these figures setting $m=2$ yields ${\textbf{\textit{W}}}_{1} \in \mathbb{R}^{25 \times 2}$ . The horizontal axes in each of the six sub-figures here denote the first and second eigenvectors in ${\textbf{\textit{W}}}_{1}$ , i.e.

(19) \begin{equation}{\textbf{\textit{W}}}_{1}=\left[\begin{array}{cc}{\textbf{\textit{w}}}_{1} & {\textbf{\textit{w}}}_{2}\end{array}\right]\end{equation}

The superscripts in the sub-figures denote the blade—i.e. (a), (b) or (c)—and the subscripts denote the objective—i.e. E for efficiency and P for pressure ratio. The markers in Figs 8 and 9 are then obtained by projecting each ${\textbf{\textit{x}}}_{i}$ onto these coordinates, with the vertical axis showing the values of the efficiency and pressure ratio.

It is apparent that on their respective subspaces, Blades A and C admit a quadratic variation in efficiency and a largely linear variation in pressure ratio. What is interesting is that for Blade B, which is a low-speed blade, the trend in efficiency is relatively linear. Further examination of these results reveals that the active subspace imparts a change in blade camber to the geometries and, depending on the amount of camber incorporated, the impact of the shock can either be advantageous or debilitating. For Blade B, given its relatively lower operating Mach number, no quadratic trend in efficiency is observed.

The results above attest to the notion that the active subspaces are driven by underlying physical phenomena and can be found across various blade geometries—note that the nominal geometries are different—at different operating conditions. But can we re-use the active subspaces obtained for one blade on another? In Fig. 10, we plot the efficiency and pressure ratio values of Blades B and C on the active subspace associated with Blade A. This sufficient summary plots clearly captures the trends we observed above and demonstrates that a new blade geometry with the same parameterisation can leverage the subspace discovered from another blade, even though each blade’s active subspace is different (see Fig. 11). This implies that for a new blade, we may not require as many CFD evaluations as for the three blades above as we can utilise the subspace previously computed.

Figure 10.

Sufficient summary plots for (a) efficiency and (b) pressure ratio values, for blades A, B and C plotted on the active subspace associated with Blade A.

Figure 11.

Components of the (a) first and (b) second eigenvectors of the covariance matrix for the efficiency active subspace for all three blades.

We end this section with a remark on the computational cost of approximating ${\textbf{\textit{W}}}$ and, by association, its active subspace. The cost of a global quadratic model scales poorly with dimension, motivating a different strategy to effectively compute the active subspace: for a 50-dimensional problem the number of coefficients that need to be estimated rises to 1326; for a 100-dimensional problem to 5151, etc.

3.1 Theory and computational approaches

As a necessary pre-processing step, the inputs to all the aforementioned methods need to be standardised, aimed at the reduction of multicollinearity (see Section 5.3 in Ref. (Reference Faraway7)). Analogous to ideas for computing the active subspace, SDR methods require the construction of a covariance matrix, the dominant eigenvectors of which define directions within the input parameter space that are important. While for active subspaces our metric of importance was related to $\nabla_{{\textbf{\textit{x}}}}\,f$ , SDR techniques adopt different metrics. In SIR, for instance, the covariance matrix is given by

(20) \begin{equation}\tilde{{\textbf{\textit{K}}}}=\frac{1}{N}\sum_{s=1}^{S}N_{s} \boldsymbol{\mu}_{s} \boldsymbol{\mu}_{s}^{T}\end{equation}

where

(21) \begin{equation}\boldsymbol{\mu}_{s}=\frac{1}{N_{s}}\sum_{j\in\mathcal{J}_{s}} {\textbf{\textit{x}}}_{j}\end{equation}

Here $s=1, \ldots, S$ represents the slices into which the output values $f\!\left( {\textbf{\textit{x}}}_{i} \right)$ for $i=1, \ldots, N$ must be uniformlyFootnote 2 placed, such that each slice contains $N_s$ samples, i.e. $N_{s}=\left|\mathcal{J}_{s}\right|$ , where $\mathcal{J}_{s} \in [ 1, \ldots, N ]$ . The term $\boldsymbol{\mu}_{s}$ is the sample mean of the input parameters in the s-th slice, making SIR a method based on the first moment only. Additionally, due to its dependence on the first moment, SIR fails to handle symmetric functions.

Once the covariance matrix in Equation (20) has been constructed, its first few eigenvectors can inform important directions. One apparent advantage of SIR is that it does not require gradient information. Further, while the number of slices S may alter the asymptotic variance of the output, the effect may not be significant (see Remark 4.3 in Ref. (Reference Li17)).

A similar slicing strategy is adopted in SAVE, a technique based on both the mean and the variance. Here the covariance matrix has the form

(22) \begin{equation}\tilde{{\textbf{\textit{K}}}}=\frac{1}{N}\sum_{s=1}^{S}N_{s}\left( {\textbf{\textit{I}}}- \boldsymbol{\Omega}_{s} \right)^{2}\end{equation}

where sample covariances are given by

(23) \begin{equation}\boldsymbol{\Omega}_{s}=\frac{1}{N_{s}-1}\sum_{j\in\mathcal{J}_{s}}\left( {\textbf{\textit{x}}}_{i}- \boldsymbol{\mu}_{s}\right)\left( {\textbf{\textit{x}}}_{i}- \boldsymbol{\mu}_{s}\right)^{T}\end{equation}

and where ${\textbf{\textit{I}}}$ is the identity matrix.

The next SDR approach we introduce is pHd, a technique that does not require slices, differentiating it from SIR and SAVE. As Li^{(Reference Li18)} states, “[pHd] begins with the observation that the eigenvectors for the Hessian matrices of the regression function are helpful in the study of the shape of the regression surface.” The method tries to locate the principal directions along which the regression surface, on average, exhibits the greatest curvature variation. Thus, its covariance matrix seeks to approximate the average Hessian of the function, i.e. $\textbf{K} = \mathbb{E}[\nabla^2 f({\textbf{\textit{x}}})]$ . To achieve this, Li^{(Reference Li18)} proposes the following approximation

(24) \begin{equation}\tilde{\textbf{K}}=\frac{1}{N}\sum_{i=1}^{N}\left(f\!\left(\textbf{\textit{x}}_{i}\right)-\bar{f}\left( {\textbf{\textit{x}}}\right) \right)\left({\textbf{\textit{x}}}_{i}-\bar{{\textbf{\textit{x}}}}\right)\left({\textbf{\textit{x}}}_{i}-\bar{{\textbf{\textit{x}}}}\right)^{T},\end{equation}

where

(25) \begin{equation}\bar{\textbf{\textit{x}}}=\frac{1}{N}\sum_{i=1}^{N}{\textbf{\textit{x}}}_i\end{equation}

is the sample mean for the inputs and

(26) \begin{equation}\bar{f} \left( {\textbf{\textit{x}}}\right) =\frac{1}{N}\sum_{i=1}^{N} f\!\left( {\textbf{\textit{x}}}_{i}\right)\end{equation}

is the sample mean for the outputs. It should be noted that, in general, pHd is less stable than SIR, due to its dependence on second-order moments.

The last SDR method we shall review here is CR^{(Reference Li, Zha and Chiaromonte16)}, where empirical directions are defined as the vectors

(27) \begin{equation}{\textbf{\textit{x}}}_i-{\textbf{\textit{x}}}_j \; \; \forall \; \; 1 \leqslant i\neq j \leqslant N\end{equation}

The method extracts a subset of empirical directions with small variations in $f\!\left({\textbf{\textit{x}}} \right)$ . Then the directions with largest variations in $f\!\left({\textbf{\textit{x}}} \right)$ become estimates for the dimension reducing subspace. The covariance matrix for CR is given by

(28) \begin{equation}\hat{ {\textbf{\textit{K}}}}(c)=\frac{2}{N(N-1)}\sum_{1\leqslant i\neq j \leqslant N}({\textbf{\textit{x}}}_j-{\textbf{\textit{x}}}_i)({\textbf{\textit{x}}}_j-{\textbf{\textit{x}}}_i)^{\rm T} \mathcal{I} \left( \lvert f\!\left({\textbf{\textit{x}}}_j \right) - f\!\left({\textbf{\textit{x}}}_i \right) \rvert \leqslant c \right)\end{equation}

where

(29) \begin{equation}\mathcal{I} \left( \lvert f\!\left({\textbf{\textit{x}}}_j \right) - f\!\left({\textbf{\textit{x}}}_i \right) \rvert \leqslant c \right)=\begin{cases}1 \; \; \;\text{for } \lvert f\!\left({\textbf{\textit{x}}}_j \right) - f\!\left({\textbf{\textit{x}}}_i \right) \rvert \leqslant c\\[3pt] 0 \; \; \;\text{for } \lvert f\!\left({\textbf{\textit{x}}}_j \right) - f\!\left({\textbf{\textit{x}}}_i \right) \rvert > c\\\end{cases}\end{equation}

and c is a user-defined tolerance. The dimension reducing subspaces are then found by selecting eigenvectors corresponding to the smallest eigenvalues. One drawback of CR is that the tolerance c needs to be chosen.

3.2 Results on Blade A

In what follows we apply the four SDR methods discussed above on the data-set associated with Blade A, focusing only on the efficiency objective. Recall that our overarching aim in this section is to ascertain whether SDR rooted methods can identify the dimension reducing subspaces—similar to those shown in Section 2—on a more parsimonious computational budget. To this end, we randomly select 200 input-output ( ${\textbf{\textit{x}}}_{i}, f\!\left( {\textbf{\textit{x}}}_{i} \right) $ ) pairs from the computed 548 pairs. These 200 pairs are then given as inputs to SIR, SAVE, pHd and CR. The computed subspaces obtained using these techniques are shown in Fig. 12 in the form of sufficient summary plots. Once the subspaces have been computed, we project the remaining 348 samples as well.

Figure 12.

Sufficient summary plots for the Blade A data-set for: (a) SIR; (b) SAVE; (c) pHd; (d) CR.

To understand the effect of the number of samples (and our random subsampling), we introduce a convergence metric

(30) \begin{equation}\phi=\arcsin \left( \left\Vert \tilde{{\textbf{\textit{W}}}}_{1}\tilde{{\textbf{\textit{W}}}}_{1}^{T}-{\textbf{\textit{W}}}_{1} {\textbf{\textit{W}}}_{1}^{T} \right\Vert _{2} \right)\end{equation}

defined as the angle between the subspace approximated by the various SDR methods, $\tilde{{\textbf{\textit{W}}}}_{1}$ , and the subspace obtained through the quadratic model via the eigendecomposition of Equation (16) yielding ${\textbf{\textit{W}}}_{1}$ . We report the change in $\phi$ as a function of the number of subsamples for the four different SDR strategies in Fig. 13 with 20 bootstrap replicates for a given number of samples. Our study shows that SAVE, CR and pHd fail to identify the underlying structure in the data, while SIR obtains an approximation that is relatively closer to that of the quadratic model.

Figure 13.

Subspace distance angles between the SDR and active subspace. The grey markers show the outcomes of 20 bootstrap replicates, while the red line shows their mean values.

We remark here that although we have used active subspaces as a basisFootnote 3 for comparison, we cannot state that the active subspace is optimal. What we can state is that the active subspaces offers a subspace where statements on the gradient of the residual $\left\Vert f\!\left( {\textbf{\textit{x}}}\right)-g\left({\textbf{\textit{W}}}_{1}^{T} {\textbf{\textit{x}}}\right)\right\Vert _{2}^{2}$ can be made. It can serve as an initial solution, if not local minima, to the problem of approximating a high dimensional function over an $m-$ dimensional subspace, but multiple other subspaces can exist. Thus, care should be taken when interpreting Fig. 13 in the absence of Fig. 12 as the latter is just as important.

While numerous examples of successful deployment of the statistical sufficient dimension reduction methods considered here exist in literature, its worth inquiring why in particular CR and pHd they have fallen short of discovering a visually identifiable dimension reducing subspace. One obvious reason for this is the rather limited number of samples that have been used here to find the dimension reducing subspace; in particular for pHD as we are estimating the Hessian, more samples are required.

4.1 Theory and computational approaches

Let us start by redefining Equation (10) as

(31) \begin{align}\begin{split}f\!\left( {\textbf{\textit{x}}} \right) & \approx g \left( {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}} \right) \\[3pt]& = \sum_{i=1}^{m}\alpha_{i} \boldsymbol{\psi}_{i} \!\left( {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}\right)\end{split}\end{align}

where g is a polynomial in $n<<m$ dimensions with a maximum order of p with cardinality $\rho$ , and where ${\textbf{\textit{U}}} \in \mathbb{R}^{m \times n}$ is an orthogonal matrix with $m=25$ for our studies. The matrix ${\textbf{\textit{U}}}$ belongs to the Grassman manifold of n-dimensional subspaces in $\mathbb{R}^{m}$ . We express g as a linear combination of coefficients $\alpha_{i}$ for $i=1, \ldots, n$ and multivariate $l_2$ orthogonal polynomial basis terms $\boldsymbol{\psi}_{i}$ . These multivariate polynomials are products of composite univariate polynomials $\psi$

(32) \begin{align}\begin{split} \boldsymbol{\psi}_{{\textbf{\textit{q}}}}\!\left( {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}\right) & =\prod_{j=1}^{n}\psi_{q_{j}} \left(\boldsymbol{u}_{j}^{T} {\textbf{\textit{x}}} \right) \\[3pt] \Rightarrow \boldsymbol{\psi}_{{\textbf{\textit{q}}}}\!\left({\kern1.2pt}\boldsymbol{y} \right) & =\prod_{j=1}^{n}\psi_{q_{j}} \!\left( y^{\left( j \right) } \right)\end{split}\end{align}

where the superscripts are used to designate the relevant dimension. Here the symbol ${\textbf{\textit{q}}} = (q_1, \ldots, q_n) \in \mathbb{N}^{n}$ is a finite set of indices that belongs to a multi-index set $\mathcal{I}$ . In this paper, we use both isotropic tensor product and isotropic total order bases, the difference between them being the number of basis terms, i.e. the cardinality of the index set. For a maximum degree of p, the former is given by the rule $\text{max} \left( q_j \right) \leq p$ , while for the total order basis $\sum_{i=1}^{n}q_{i}\leq p$ (see Section 1.2 in Ref. (Reference Seshadri, Narayan and Mahadevan22) for further details). The notation $\boldsymbol{u}_{j}$ denotes the j-th column of ${\textbf{\textit{U}}}$ and $y^{ \left( j\right) } = \boldsymbol{u}^{T}_{j} {\textbf{\textit{x}}}$ for $j=1, \ldots, n$ , represents the projected coordinates on the space spanned by the columns of ${\textbf{\textit{U}}}$ . Gradients of Equation (32) with respect to $\boldsymbol{y} = {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}$ can be easily computed via

(33) \begin{equation}\frac{\partial \boldsymbol{\psi}_{q}}{\partial y^{(l)}}=\frac{d\psi_{ql}}{dy^{(l)}}\prod_{j=1,j\neq l}\psi_{qj}\left(y^{\left(j\right)}\right).\end{equation}

These gradients will be required in subsequent computations.

Our main objective in this section is to solve the optimisation problem

(34) \begin{equation}\underset{ \stackrel{\alpha_{i}, \; i=1,\ldots,n}{\text{range}\; {\textbf{\textit{U}}} \in \mathbb{G}(n, \mathbb{R}^m) } }{ \text{minimise} } \; \; \left\Vert f\!\left({\textbf{\textit{x}}}\right)- \sum_{i=1}^{n}\alpha_{i} \boldsymbol{\psi}_{i} \!\left( {\textbf{\textit{U}}}^{T}x\right) \right\Vert^{2}_{2}\end{equation}

which can be expressed as

(35) \begin{equation}\underset{\boldsymbol{ \alpha}, {\textbf{\textit{U}}} }{\text{ minimise} }\left\Vert \boldsymbol{f}- {\textbf{\textit{P}}}\left({\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}\right) \boldsymbol{\alpha} \right\Vert _{2}^{2}\end{equation}

where $\boldsymbol{f}=\left[f\!\left({\textbf{\textit{x}}}_{1}\right) \dots f\!\left({\textbf{\textit{x}}}_{K}\right)\right]^{T}$ , $\boldsymbol{\alpha} =\left[\alpha_1 \dots \alpha_n \right]^{T}$ and

(36) \begin{equation}{\textbf{\textit{P}}}(i, j) = \boldsymbol{\psi}_{j} \left( {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}_{i} \right)\end{equation}

Using the notion of variable projection^{(Reference Golub and Pereyra11)} (also see^{(Reference Golub and Pereyra10)} and Section 8.4 in Ref. (Reference Hansen, Pereyra and Scherer12)), one can rewrite Equation (35) as an optimisation problem over ${\textbf{\textit{U}}}$ only

(37) \begin{align}\begin{split}\boldsymbol{h} &= \underset{ {\textbf{\textit{U}}} }{\text{ minimise} }\left\Vert \boldsymbol{f}- {\textbf{\textit{P}}}\left({\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}\right) {\textbf{\textit{P}}}^{\dagger}\left({\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}\right) \boldsymbol{f} \right\Vert _{2}^{2} \\[3pt] &= \underset{ {\textbf{\textit{U}}} }{\text{ minimise} }\left\Vert \left( {\textbf{\textit{I}}} - {\textbf{\textit{P}}}\left({\textbf{\textit{U}}}^{T} \right) {\textbf{\textit{P}}}^{\dagger}\left({\textbf{\textit{U}}}^{T} \right) \right) \boldsymbol{f} \right\Vert _{2}^{2} \\[3pt] &= \underset{ {\textbf{\textit{U}}} }{\text{ minimise} }\left\Vert {\textbf{\textit{r}}} \left( {\textbf{\textit{U}}} \right) \right\Vert _{2}^{2}\end{split}\end{align}

where ${\textbf{\textit{r}}}$ is known as the projector onto the orthogonal complement of the column space of the matrix ${\textbf{\textit{P}}} \left({\textbf{\textit{U}}}^{T} \right)$ and the matrix $ {\textbf{\textit{P}}}^{\dagger}$ is the generalised Moore-Penrose inverse of ${\textbf{\textit{P}}}$ . Now ${\textbf{\textit{r}}} \in \mathbb{R}^{M}$ where M represents the total number of samples. The first step in the variable projection approach is to formulate both the gradient and the Hessian (see Section 3 of Ref. (Reference Golub and Pereyra10))

(38) \begin{equation}\frac{1}{2}\nabla\left\Vert {\textbf{\textit{r}}}\left({\textbf{\textit{U}}}\right)\right\Vert _{2}^{2}= {\textbf{\textit{J}}}^{T}\left({\textbf{\textit{U}}}\right)r\left({\textbf{\textit{U}}}\right)\end{equation}

(39) \begin{equation}\frac{1}{2}\nabla^{2}\left\Vert {\textbf{\textit{r}}} \left( {\textbf{\textit{U}}}\right)\right\Vert _{2}^{2}= {\textbf{\textit{J}}}^{T}\left( {\textbf{\textit{U}}}\right){\textbf{\textit{J}}}( {\textbf{\textit{U}}})+\sum_{i=1}^{M}{\textbf{\textit{r}}}_{i}\left({\textbf{\textit{U}}}\right)\nabla^{2} {\textbf{\textit{r}}}_{i}\left({\textbf{\textit{U}}}\right)\end{equation}

The Jacobian ${\textbf{\textit{J}}} \in \mathbb{R}^{M \times m \times n}$ , a tensor, can be written as

(40) \begin{equation}{\textbf{\textit{J}}} \left(:, j, k \right) =-\left(\left({\textbf{\textit{I}}}- {\textbf{\textit{P}}} {\textbf{\textit{P}}}^{\dagger}\right)\frac{\partial {\textbf{\textit{P}}}}{\partial {\textbf{\textit{U}}}_{\left(j,k\right)}} {\textbf{\textit{P}}}^{-}\right)\boldsymbol{f}-\left(\left({\textbf{\textit{I}}}- {\textbf{\textit{P}}} {\textbf{\textit{P}}}^{\dagger}\right)\frac{\partial {\textbf{\textit{P}}}}{\partial {\textbf{\textit{U}}}_{\left(j,k\right)}} {\textbf{\textit{P}}}^{-}\right)^{T} \boldsymbol{f}\end{equation}

The symbol ${\textbf{\textit{P}}}^{-}$ denotes the generalised inverse; as stated in Ref. (Reference Archer, Saltelli and Sobol1), the full Moore-Penrose pseudoinverse is not required, and thus a symmetric generalised inverse, for which the two conditions

(41) \begin{equation}\left({\textbf{\textit{P}}} {\textbf{\textit{P}}}^{-} {\textbf{\textit{P}}}\right)= {\textbf{\textit{P}}} \; \; \; \; \text{and} \; \; \; \; \left({\textbf{\textit{P}}} {\textbf{\textit{P}}}^{-}\right)^{T}={\textbf{\textit{P}}} {\textbf{\textit{P}}}^{-}\end{equation}

hold, will suffice. We also remark here that the power of variable projection stems from the fact that the expression in Equation (40) does not involve gradients of the Moore-Penrose pseudoinverse. To compute the terms in Equation (40), one simply uses Equation (33) with Equations (35) and (36). Details on how these terms are computed, along with properties of the Jacobian are provided in Ref. (Reference Hokanson and Constantine14). Further, in Ref. (Reference Hokanson and Constantine14) the authors provide a step-by-step Gauss-Newton algorithm—where the Hessian ${\textbf{\textit{J}}}^{T}\left( {\textbf{\textit{U}}}\right){\textbf{\textit{J}}}( {\textbf{\textit{U}}})$ in Equation (39) is approximated—for solving the optimisation problem in Equation (37). Our implementation of their strategy is available in the Effective Quadratures open-source package^{(Reference Seshadri and Parks23)}.

Before moving on to the results of polynomial variable projection on the blade data sets, a few comments regarding computational cost are in order. The key advantage of polynomial variable projection is that the cost of finding a dimension reducing subspace scales as a function of the underlying data’s intrinsic dimension n and not its ambient dimension m. This growth in the number of data points is drastically different to the global polynomial approach in Ref. (Reference Seshadri, Shahpar, Constantine, Parks and Adams25), which scales with the ambient dimension m. It should also be noted that adjoint evaluations can be used to estimate the covariance matrix in (11), however, well-converged adjoint evaluations are often not readily available—as in our case—thus justifying the use of polynomial variable projection.

4.2 Results on Blade A

We apply variable projection on the data-set associated with Blade A in the same manner as we did with the SDR methods described in Section 3. Figure 14 shows the results for variable projection with different numbers of samples and their bootstrap replicates. For each experiment, the subspace angle obtained from variable projection is compared with the global quadratic active subspaces approach detailed before. Here we set $p=2$ with a total order basis. Compared to all of the aforementioned techniques, variable projection does yield the lowest angles, indicating that the approximated subspaces are close to the active subspace.

Figure 14.

Subspace distance angles between the variable projection and active subspace. The grey markers show the outcome of 20 bootstrap replicates, while the red line shows their mean value.

The results for a few of these experiments are shown in Fig. 15. Here sub-figures (a), (b) and (c) show the sufficient summary plots, where the 2D coordinates in the figure correspond to the columns of ${\textbf{\textit{U}}} = [\boldsymbol{u}_{1}, \boldsymbol{u}_{2}]$ when evaluating the algorithm with $M=150, 200$ and 250 samples, respectively. These results give us confidence in the utility of variable projection for finding dimension reducing subspaces for our turbomachinery problem. More importantly, the reduced sample count—150 to 200 samples—compared to the cost of a global quadratic model—above 351 samples—represents a significant computational saving.

Figure 15.

Sufficient summary plots of the variable projection (total order degree 2) approach with: (a) 150, (b) 200 and (c) 250 randomly subsampled input-output pairs for Blade A.

In the same vein as before, Fig. 16 shows sufficient summary plots for $p=3$ using a tensor order basis; the results are similar, with reduced scatter even when the number of samples used is 150. In Fig. 17(a) we show the polynomial fit for the result in Fig. 16(c) and compare the values of this approximation and the actual CFD yielded efficiencies in Fig. 17(b). The $R^2$ value of this fit is 0.9929, indicating that this response surface can be used for subsequent inference.

Figure 16.

Sufficient summary plots of the variable projection (tensor order degree 3) approach with: (a) 150, (b) 200 and (c) 250 randomly subsampled input-output pairs for Blade A.

Figure 17.

Polynomial (cubic) response surface on the subspace in Fig. 16(c) and a comparison between this polynomial and the CFD results in (b).

We repeated the aforementioned numerical results on the data-sets associated with Blades B and C, achieving similar results. In Section 5 we expand upon these findings.

5.1 Generating constrained designs

First, we define the matrix

(42) \begin{equation}{\textbf{\textit{W}}}=\left[\begin{array}{c@{\quad}c}{\textbf{\textit{U}}} & {\textbf{\textit{V}}}\end{array}\right]\end{equation}

where ${\textbf{\textit{V}}}=\textsf{null}\!\left({\textbf{\textit{U}}}^{T}\right)$ , in which $\textsf{null}$ represents the nullspace, and ${\textbf{\textit{U}}}$ is the dimension reducing subspace obtained from variable projection. As ${\textbf{\textit{W}}}$ is orthogonal, we can write

(43) \begin{align}\begin{split}& -\boldsymbol{1}\leq {\textbf{\textit{x}}} \leq \boldsymbol{1} \\\Rightarrow & -\boldsymbol{1}\leq {\textbf{\textit{W}}} {\textbf{\textit{W}}}^{T} {\textbf{\textit{x}}} \leq \boldsymbol{1} \\\Rightarrow & -\boldsymbol{1}\leq {\textbf{\textit{U}}} {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}} + {\textbf{\textit{V}}} {\textbf{\textit{V}}}^{T} {\textbf{\textit{x}}} \leq \boldsymbol{1} \\\Rightarrow & -\boldsymbol{1}\leq {\textbf{\textit{U}}} \boldsymbol{y} + {\textbf{\textit{V}}} \boldsymbol{z} \leq \boldsymbol{1}\end{split}\end{align}

where $\boldsymbol{y} = {\textbf{\textit{U}}}^{T} {\textbf{\textit{x}}}$ and $\boldsymbol{z} = {\textbf{\textit{V}}}^T {\textbf{\textit{x}}}$ . This implies that one can generate numerous ${\textbf{\textit{x}}}$ that have the same coordinates $\boldsymbol{y}$ but different coordinates $\boldsymbol{z}$ , under the constraint that ${\textbf{\textit{x}}}$ still lies within the original bounding box. This can be conveniently expressed as a linear programming problem where for the same coordinate $\boldsymbol{y}$ we wish to generate numerous $\boldsymbol{z}$ .

5.2 Multi-point results for Blade C

We apply such a linear programming strategy to generate samples that map to the different subspaces; our results are summarised in Fig. 19. In Fig. 19(a) we plot the contours of sea-level efficiency on the sea-level subspace (essentially a different perspective on Fig. 17(a), but for Blade C), with the initial geometry of Blade C denoted by the $\times$ marker and 7 sample 2D coordinates denoted by the white circular markers. For each marker we generate four 25D vectors that satisfy Equation (43), i.e. they have different values along $\boldsymbol{z}$ but the same value along $\boldsymbol{y}$ . These 28 new designs are then projected on the cruise efficiency subspace, shown in Fig. 19(b). We repeat this exercise for 10 high-efficiency designs in the cruise subspace, shown by the dark grey circular markers. It is clear that designs that attain higher efficiencies at cruise have a lower efficiency at sea-level and vice versa. Further, it is interesting to note that the original design reflects a healthy compromise between maximum efficiency at sea-level and at cruise. It is important to emphasise that we can infer the efficacy of various design configurations at these two operating points by visualizing the location of their $\boldsymbol{y}$ coordinates along the two subspaces.

In Fig. 19, we show the range (maximum - minimum) in efficiency for sea-level and cruise designs on both subspaces. On the sea-level subspaces, both the cruise and sea-level designs selected have an efficiency range of 0.6 $\%$ and 0.5 $\%$ respectively. These same design perturbations, when mapped on to the cruise subspace yield efficiency ranges of 0.1 $\%$ and 3.06 $\%$ respectively. There are two important remarks to make here. The first (and the obvious one) is that design perturbations aimed at improving the efficiency at sea-level are likely to have a detrimental effect at cruise. Second, higher efficiency designs at cruise (see dark grey markers in Fig. 19(b)) warrant further deviations from the datum design; inferred from their coordinates along $\boldsymbol{u}_{1}$ . Analogously, higher efficiency designs at sea level (see white markers in Fig. 19(a)) have comparatively reduced design excursion. However, despite this, these designs yield a greater efficiency penalty. This leads to the conclusion that greater design perturbations can potentially yield more favorable designs—a particularly important point when investigating the impact of manufacturing variations in design.

To gain a better understanding of the differences between the geometries, consider the parallel coordinates plot shown in Fig. 20. Here the vertical axis has been scaled by multiplying the dimensionalised design variables of the various blades with the absolute value of the cruise subspace $\boldsymbol{u}_{1}$ shown previously in Fig. 18(a). We remark here that this is purely to distinguish between the leading and trailing recambering dofs, which were found be to more important than the other three dofs. A distinct trend in parallel coordinates is readily apparent. Higher efficiency designs at cruise have both—on average—positive leading and trailing edge recambering. However, for higher efficiency designs at sea-level, the trend is reversed. The blade shapes for these cases are shown in Fig. 21. In addition to recambering, minor differences in dihedral are also evident.

Figure 20.

A parallel coordinates representation for higher efficiency designs at both cruise and sea-level conditions.

Figure 21.

Comparison of higher efficiency designs at sea-level in (a, b and c) and at cruise in (d, e and f). Higher efficiency designs are shown in grey with the baseline shown in magenta. Geometries are scaled for proprietary reasons.