Supporting Multi-point Fan Design with Dimension Reduction

10/20/2019
by   Pranay Seshadri, et al.
15

Motivated by the idea of turbomachinery active subspace performance maps, this paper studies dimension reduction in turbomachinery 3D CFD simulations. First, we show that these subspaces exist across different blades—under the same parametrization—largely independent of their Mach number or Reynolds number. This is demonstrated via a numerical study on three different blades. Then, in an attempt to reduce the computational cost of identifying a suitable dimension reducing subspace, we examine statistical sufficient dimension reduction methods, including sliced inverse regression, sliced average variance estimation, principal Hessian directions and contour regression. Unsatisfied by these results, we evaluate a new idea based on polynomial variable projection—a non-linear least squares problem. Our results using polynomial variable projection clearly demonstrate that one can accurately identify dimension reducing subspaces for turbomachinery functionals at a fraction of the cost associated with prior methods. We apply these subspaces to the problem of comparing design configurations across different flight points on a working line of a fan blade. We demonstrate how designs that offer a healthy compromise between performance at cruise and sea-level conditions can be easily found by visually inspecting their subspaces.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 3

page 4

page 5

page 25

page 26

02/01/2018

Dimension Reduction via Gaussian Ridge Functions

Ridge functions have recently emerged as a powerful set of ideas for sub...
02/12/2020

On sufficient dimension reduction via principal asymmetric least squares

In this paper, we introduce principal asymmetric least squares (PALS) as...
02/26/2021

Ensemble Conditional Variance Estimator for Sufficient Dimension Reduction

Ensemble Conditional Variance Estimation (ECVE) is a novel sufficient di...
04/14/2022

itdr: An R package of Integral Transformation Methods to Estimate the SDR Subspaces in Regression

Sufficient dimension reduction (SDR) is a successful tool in regression ...
08/10/2015

Model-based SIR for dimension reduction

A new dimension reduction method based on Gaussian finite mixtures is pr...
04/25/2022

A Time-Triggered Dimension Reduction Algorithm for the Task Assignment Problem

The task assignment problem is fundamental in combinatorial optimisation...
10/24/2017

Impossibility of dimension reduction in the nuclear norm

Let S_1 (the Schatten--von Neumann trace class) denote the Banach space ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In the modern turbomachinery design process, blade shapes are parametrized by a set of design variables useful from aerodynamic, structural and manufacturing perspectives. Typically, the number of such parameters can range from 20 to 250. This high-dimensional design space makes parametric studies—e.g., uncertainty quantification, optimization and sensitivity analysis—challenging and computationally prohibitive, motivating the use of dimension reduction strategies. Within the broad field of dimension reduction there exist two basic approaches: subspace- and subset-based dimension reduction. Subset-based dimension reduction aims to find an important222Importance here is usually characterized by the conditional variance on the output, although other metrics can also be adopted.

subset of original variables. Methods such as ridge regression, Lasso

[friedman2001elements] and ANOVA-based decompositions and the related Sobol’ indices [sobol2001global, archer1997sensitivity] can be categorized as methods for subset-based dimension reduction. Our interest in this paper, however, lies in ideas for subspace-based dimension reduction, where one seeks to identify specific linear combinations of all the parameters that are important with respect to the output.

This body of work by no means represents the first foray into turbomachinery-based dimension reduction. Prior papers that investigated the use of dimension reduction strategies in turbomachinery design include the works of Banamonde et al. [bahamonde2017active] and Qin et al. [qin2016flow]. In the former publication, the authors use active subspaces [constantine2015active]—a set of ideas for parameter space dimension reduction—to develop a reduced-order model for the design of a turbine power cycle. To compute the active subspaces, the authors used finite differences, a feasible approach given that their governing model was a mean-line code. In [qin2016flow]

the authors use a 2D CFD model of a NACA65-1210 airfoil and leverage a one-dimensional active subspaces heuristic (see page 5 in

[constantine2015active]) to estimate the impact of blade fouling and erosion.

More recently (and of greater relevance to this work), Seshadri et al. [seshadri2018turbomachinery] detail methods for generating turbomachinery active subspace performance maps. These 2D contour plots illustrate and quantify the variation of key flow performance metrics with different designs—even when the number of design variables is greater than two. The authors generate such maps for an isolated 3D fan blade (see Fig. 24 in [seshadri2018turbomachinery]), where the contours reflect the changes in pressure ratio, flow capacity and efficiency with a change in the design of the blade. Their work raises several critical questions—both from a computational perspective and from a turbomachinery design standpoint—-motivating our paper. As the study in [seshadri2018turbomachinery] was carried out on a single isolated blade, it is logical to ask whether those results can be generalized to other blades—even assuming the same design parametrization. We explore this line of questioning in Sect. 2 on three different modern fan blades. The computational workhorse for computing the dimension reducing subspaces in [seshadri2018turbomachinery] was an active subspace technique that leverages a global quadratic polynomial approximation. In Sect. 3 we explore a variety of alternative classical approaches. Motivated by the variable projection aided dimension reduction approach in [hokanson2017data], we devote Sect. 4 to its study, offering a few minor algorithmic variations. Then, we test this on one of our blades, being parsimonious in the number of samples (and hence CFD evaluations) used. As turbomachinery components need to operate at multiple points across a characteristic, in Sect. 5

we leverage variable projection and a pressure ratio interpolating strategy to obtain design maps at two points on the characteristic and discuss their utility within the design process.

2 Turbomachinery Active Subspace Performance Maps

In this study, three research fan blades are considered, shown in Fig. 1. Blades A and B are high-speed fan blades, with real engine representative pressure ratios and Mach numbers, while Blade C is a low-speed research fan blade.

Figure 1: The three fan blades considered in this study (images have been appropriately scaled for proprietary reasons).

Each blade uses the same 25D parametrization scheme in [seshadri2018turbomachinery]

, defined by 5 degrees of freedom (dof) specified at 5 spanwise locations. These dofs include leading and trailing edge recambering, skew, sweep and dihedral as shown in Fig. 

2.

Figure 2: The design space of Blades A, B and C. Shown are the five degrees of freedom, where the ranges have been exaggerated for clarity.

2.1 Computational grid and boundary conditions

The Rolls-Royce code PADRAM [PADRAM] 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).
Property Blade A Blade B Blade C
No. of cells in the tip gap 14 10 14
No. of cells in the spanwise direction 129 159 139
No. of cells in the circumferential direction 59 55 54
No. of cells in the axial direction 177 323 281
Total number of cells (millions) 1.7 3.345 3.153
Table 1: Mesh characteristics for the three blades.

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 Figure 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)

and

(2)

where the temperature ratio () is given by

(3)

In Eqs. (1), (2) and (3): indicates the stagnation pressure; the stagnation temperature; 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 utilize a fixed exit capacity boundary condition:

(4)
(5)

The endwalls (hub, blade, casing and splitter) are characterized 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 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 Figure 4(a) plot the Efficiency on the horizontal axis (logarithm-base 10 scale), defined as

(6)

Similarly, in (b), Pressure ratio is given by

(7)

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.

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 [crumpton1998unstructured, lapworth2004hydra]. To close the set of equations when applying a Reynolds-based averaging to the Navier-Stokes, we utilize the one equation Spalart-Allmaras turbulence model [spalart1992one]. Further details on the calculation of the inviscid fluxes, the viscous fluxes and convergence can be found in pages 10-21 of [ostling2011potential]. We omit a detailed flow validation here and direct interested readers to references [Giulio, Seshadri_LEAK], 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.

Property Blade A Blade B Blade C (cruise) Blade C (sea-level)
Table 2: Reynolds number and inlet boundary conditions.

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 [cook2009introduction]). Consider the 4D function

(8)

defined over the hypercube and where . We can express this function as

(9)

The transformation induced by the vector

projects the 4D values in onto the subspace—i.e., linear combination of the parameters—given by . Thus, our 4D function is actually 1D over this subspace! Imagine we generate random samples from and project each for onto this subspace. Then, we evaluate (8) for each and plot the resulting data . We would arrive at Fig. 5(a), where the horizontal axis of plots the subspace and the vertical axis plots the function values .

Now, typically, we may not know 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 . 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.

Figure 5: Projecting multi-dimensional samples of the function over a subspace: (a) (b)

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 as a set of design parameters and assume . Without loss of generality, we will assume that has been appropriately non-dimensionalized such that . Our goal in this subsection is to detail a strategy for reducing the dimensionality of our problem, via

(10)

where where is the number of design parameters (here 25) and is the reduced dimension. In Eq. (10), is a suitably chosen -dimensional function. Ideally, one would like to be either 1 or 2, facilitating easy visualization. However, this may not necessarily be possible and depends on the relationship between and . So, how exactly do we find ? Constantine [constantine2015active] lays the theoretical foundations for approximating and offers numerous practical computational heuristics. In essence, one seeks to estimate the covariance of the gradient, given by

(11)

where . Here the measure is the density of the inputs over their parameter space . As we can control our sampling strategy, we set to be the uniform measure and generate samples

randomly with respect to the uniform distribution over

. Then, one can approximate

using the first few eigenvectors associated with

(12)

Intuitively, one can interpret as the directions along which varies the most, and consequently as the directions along which is (close to) zero.

In the absence of gradients or adjoints (as in [seshadri2018turbomachinery]) one can approximate by constructing a global quadratic model in

(13)

where the coefficients , and are estimated via least squares. The number of coefficients associated with this quadratic model is given by

(14)

where is the highest degree of the polynomial and 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 required to fit our model should be

(15)

where 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 . Such a basis was used in [seshadri2018turbomachinery], resulting in a total of 351 coefficients (setting and 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)

from which can be readily obtained following the eigendecomposition of this symmetric matrix. Note that in (16), we have assumed a uniform measure

(17)

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 uniformly distributed random samples, i.e., Monte Carlo samples, within the hypercube

(18)

Each sample 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.

The eigenvalues for the covariance matrix

for these two objectives across for Blade A is shown in Figs. 6. These plots are useful to ascertain the dimensionality of the active subspace; Constantine [constantine2015active] suggests that we fix based on the appearance of large gaps (see page 37 in [constantine2015active]) 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 Figs. 6(a) and (b) for the two objectives. In Figs. 6(c) and (d) we increase this to 500 subsamples.

Figure 6: Eigenvalues of the covariance matrix 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 and for Blade B, . 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 for Blade B too.

Figure 7: Eigenvalues of the covariance matrix for Blade B for: (a) efficiency (b) pressure ratio, and for Blade C for (c) efficiency (d) pressure ratio.
Blade No. of CFD evaluations Description
A 548 High-speed fan
B 381 Low-speed fan
C 547 High-speed fan
Table 3: Number of CFD evaluations in DoE processes.
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 yields . The horizontal axes in each of the six sub-figures here denote the first and second eigenvectors in , i.e.

(19)

The superscripts in the sub-figures denote the blade—i.e. (a), (b) or (c)—and the subscripts denote the objective—i.e.  for efficiency and for pressure ratio. The markers in Figs. 8 and 9 are then obtained by projecting each 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 parameterization 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 utilize 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 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 Statistical Sufficient Dimension Reduction

In this section, we revisit well-known statistical strategies for sufficient dimension reduction (SDR), in the hope of addressing the computational cost of the quadratic approach detailed previously. It should be noted that these statistical techniques are rooted in classical regression, while ideas within the remit of active subspaces are rooted in the field of simulation-driven modeling. The first difference between these two schools of thought is that statistical methods often assume a zero-mean random error in the model. However, for computer experiments, which are usually deterministic, there is no such error. The second difference is that the joint probability density function of the input samples in computer experiments is known. We point the interested reader to Sect. 2.2 of

[seshadri2018dimension] for further details on these differences and to the work of Glaws et al. [glaws2018inverse, glaws2017inverse], who provide a theoretical and computational framework for utilizing two well-known SDR techniques with computer models.

With this preface, we discuss four classical dimension reduction approaches: sliced inverse regression (SIR) [li1991sliced], sliced average variance estimation (SAVE) [cook1991comment], principal Hessian directions (pHd) [li1992principal] and contour regression (CR) [li2005contour]. Our exposition of these approaches is paired with their application to the DoE results on Blade A, and closely follows the notation adopted in Glaws et al. [glaws2018inverse].

3.1 Theory and computational approaches

As a necessary pre-processing step, the inputs to all the aforementioned methods need to be standardized, aimed at the reduction of multicollinearity (see Sect. 5.3 in [faraway2016linear]). 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 , SDR techniques adopt different metrics. In SIR, for instance, the covariance matrix is given by

(20)

where

(21)

Here represents the slices into which the output values for must be uniformly333Li [li1991sliced] relaxes this requirement in Sect. 6 of his paper, i.e. each slice need not have the same number of samples. placed, such that each slice contains samples, i.e. , where . The term is the sample mean of the input parameters in the

-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 Eq. (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 may alter the asymptotic variance of the output, the effect may not be significant (see Remark 4.3 in [li1991sliced]).

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)

where sample covariances are given by

(23)

and where

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 [li1992principal] 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. . To achieve this, Li [li1992principal] proposes the following approximation

(24)

where

(25)

is the sample mean for the inputs and

(26)

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 [li2005contour], where empirical directions are defined as the vectors

(27)

The method extracts a subset of empirical directions with small variations in . Then the directions with largest variations in become estimates for the dimension reducing subspace. The covariance matrix for CR is given by

(28)

where

(29)

and 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 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 Sect. 2—on a more parsimonious computational budget. To this end, we randomly select 200 input-output () 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)

defined as the angle between the subspace approximated by the various SDR methods, , and the subspace obtained through the quadratic model via the eigendecomposition of Eq. (16) yielding . We report the change in 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 gray 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 basis444Literally. 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 can be made (see Definition 1 of [constantine2017near]). It can serve as an initial solution, if not local minima, to the problem of approximating a high dimensional function over an 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.

4 Polynomial Variable Projection for Dimension Reduction

The SDR methods applied to our problem have clearly been unable to identify similar structures to those we observed when using the active subspaces global quadratic model. Once again, motivated by the need to parsimoniously obtain dimension reducing subspaces, we study a new idea proposed by Hokanson and Constantine [hokanson2017data], which frames the subspace discovery problem as a non-linear least squares problem; our notation closely follows theirs.

4.1 Theory and computational approaches

Let us start by redefining Eq. (10) as

(31)

where is a polynomial in dimensions with a maximum order of with cardinality , and where

is an orthogonal matrix with

for our studies. The matrix belongs to the Grassman manifold of -dimensional subspaces in . We express as a linear combination of coefficients for and multivariate orthogonal polynomial basis terms . These multivariate polynomials are products of composite univariate polynomials

(32)

where the superscripts are used to designate the relevant dimension. Here the symbol is a finite set of indices that belongs to a multi-index set

. 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 , the former is given by the rule , while for the total order basis (see Sect. 1.2 in [seshadri2017effectively] for further details). The notation denotes the -th column of and for , represents the projected coordinates on the space spanned by the columns of . Gradients of Eq. (32) with respect to can be easily computed via

(33)

These gradients will be required in subsequent computations.

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

(34)

which can be expressed as

(35)

where , and

(36)

Using the notion of variable projection [golub1973differentiation] (also see [golub2003separable] and Sect. 8.4 in [hansen2013least]), one can rewrite Eq. (35) as an optimization problem over only

(37)

where is known as the projector onto the orthogonal complement of the column space of the matrix and the matrix is the generalized Moore-Penrose inverse of . Now where represents the total number of samples. The first step in the variable projection approach is to formulate both the gradient and the Hessian (see Sect. 3 of [golub2003separable])

(38)
(39)

The Jacobian , a tensor, can be written as

(40)

The symbol denotes the generalized inverse; as stated in [golub1973differentiation], the full Moore-Penrose pseudoinverse is not required, and thus a symmetric generalized inverse, for which the two conditions

(41)

hold, will suffice. We also remark here that the power of variable projection stems from the fact that the expression in Eq. (40) does not involve gradients of the Moore-Penrose pseudoinverse. To compute the terms in Eq. (40), one simply uses Eq. (33) with Eqs. (35) and (36). Details on how these terms are computed, along with properties of the Jacobian are provided in [hokanson2017data]. Further, in [hokanson2017data] the authors provide a step-by-step Gauss-Newton algorithm—where the Hessian in Eq. (39) is approximated—for solving the optimization problem in Eq. (37

). Our implementation of their strategy is available in the Effective Quadratures open-source package

[seshadri2017effective].

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 Sect. 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 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 gray 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 when evaluating the algorithm with and 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 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 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 Sec. 5 we expand upon these findings.

5 Supporting Multi-point Design via Dimension Reduction

Our prior analysis of Blade C was carried out on its sea-level condition, i.e. the inlet stagnation pressure and the exit flow capacity boundary conditions were based on operation at sea-level. We apply the same variable projection workflow along the cruise working line, requiring another 200 CFD evaluations. It is important to note that the baseline geometries for both cruise and sea-level conditions are different, and we match the hot running shape of the datum blade (no perturbations) at these conditions. Thus, in this section we are largely interested in a comparative study of the design space perturbations as opposed to absolute geometry descriptions.

A comparison of the components of the two subspaces is shown in Fig. 18, where the markers in the figures denote the 25 design variables. Although the subspaces are different (varying more so along than along ), we can compare designs that have high efficiencies at sea-level with designs that have high efficiencies at cruise and vice-versa.

Figure 18: Comparing the cruise and sea-level subspaces : components of are shown in (a) and components of in (b).

5.1 Generating constrained designs

First, we define the matrix

(42)

where , in which represents the nullspace, and is the dimension reducing subspace obtained from variable projection. As is orthogonal, we can write

(43)

where and . This implies that one can generate numerous that have the same coordinates but different coordinates , under the constraint that

still lies within the original bounding box. This can be conveniently expressed as a linear programming problem where for the same coordinate

we wish to generate numerous .

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 summarized 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 marker and 7 sample 2D coordinates denoted by the white circular markers. For each marker we generate four 25D vectors that satisfy Eq. (43), i.e. they have different values along but the same value along . 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 gray 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 emphasize that we can infer the efficacy of various design configurations at these two operating points by visualizing the location of their coordinates along the two subspaces.

Figure 19: Contours of efficiency as a function of its subspace for Blade C at two different conditions: (a) sea-level; (b) cruise. The difference between successive contour levels is 0.2

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 gray markers in Fig. 19(b)) warrant further deviations from the datum design; inferred from their coordinates along . 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 dimensionalized design variables of the various blades with the absolute value of the cruise subspace 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 gray with the baseline shown in magenta. Geometries are scaled for proprietary reasons.

6 Conclusions

Our first goal in this paper was to study numerous dimension reduction strategies that can identify dimension reducing subspaces for 3D turbomachinery simulations. We analyzed four sufficient dimension reduction methods—namely SIR, SAVE, pHd and CR—and polynomial variable projection. In our numerical experiments we found the latter to parsimoniously identify a dimension reducing subspace—by requiring approximately 200 fewer CFD evaluations compared to a global quadratic model based active subspace technique. Finally, we demonstrated how designs that offer a healthy compromise between cruise and sea-level conditions can be easily found by visually inspecting their subspaces. This can empower further design trade-off studies and serve as an important tool for the blade designer.

7 Acknowledgments

The authors are grateful to Rolls-Royce plc for permission to publish this work. The opinions and views expressed in this paper are those of the authors and not necessarily those of Rolls-Royce.

References