Active subspaces of airfoil shape parameterizations

TL;DR

This paper investigates low-dimensional active subspaces in airfoil shape parameterizations to simplify the understanding and optimization of transonic lift and drag, providing insights for multi-objective design.

Contribution

It identifies and relates active subspaces in two airfoil parameterizations, enabling low-dimensional approximations of lift and drag for improved design insights.

Findings

Active subspaces reveal important directions affecting lift and drag.

Low-dimensional approximations relate to physical airfoil properties.

Insights inform multi-objective airfoil design.

Abstract

Design and optimization benefit from understanding the dependence of a quantity of interest (e.g., a design objective or constraint function) on the design variables. A low-dimensional active subspace, when present, identifies important directions in the space of design variables; perturbing a design along the active subspace associated with a particular quantity of interest changes that quantity more, on average, than perturbing the design orthogonally to the active subspace. This low-dimensional structure provides insights that characterize the dependence of quantities of interest on design variables. Airfoil design in a transonic flow field with a parameterized geometry is a popular test problem for design methodologies. We examine two particular airfoil shape parameterizations, PARSEC and CST, and study the active subspaces present in two common design quantities of interest, transonic lift and drag coefficients, under each shape parameterization. We mathematically relate the two parameterizations with a common polynomial series. The active subspaces enable low-dimensional approximations of lift and drag that relate to physical airfoil properties. In particular, we obtain and interpret a two-dimensional approximation of both transonic lift and drag, and we show how these approximation inform a multi-objective design problem.

Tables

Table 1. Table 1: PARSEC airfoils are characterized by eleven parameters, each affecting a physical characteristic of the airfoil.

Parameter	Description
$x_1$	upper surface domain coord.
$x_2$	lower surface domain coord.
$x_3$	upper surface shape coord.
$x_4$	lower surface shape coord.
$x_5$	trailing edge shape coord.
$x_6$	trailing edge half-thickness
$x_7$	wedge angle
$x_8$	wedge half-angle
$x_9$	upper surface inflection
$x_{10}$	lower surface inflection
$x_{11}$	leading edge radius

Table 2. Table 2: Bounds on PARSEC parameters (see Table 1 ) chosen by first finding parameters that produce the NACA 0012 and then setting ranges to be ± 20 plus-or-minus 20 \pm 20 % of the optimized parameters.

Parameter	Lower	Upper
$x_1$	0.242	0.363
$x_2$	0.242	0.363
$x_3$	0.048	0.072
$x_4$	-0.072	-0.048
$x_5$	-0.004	0.004
$x_6$	0.008	0.012
$x_7$	-2.223	-3.335
$x_8$	7.40	11.10
$x_9$	-0.400	-0.600
$x_{10}$	0.400	0.600
$x_{11}$	0.012	0.018

Table 3. Table 3: Bounds on CST parameters (see section 2.5 ) chosen by first finding parameters that produce the NACA 0012 and then setting ranges to be ± 20 plus-or-minus 20 \pm 20 % of the optimized parameters.

Parameter	Lower	Upper
$x_1$	0.12	0.18
$x_2$	0.8	1.2
$x_3$	0.8	1.2
$x_4$	0.8	1.2
$x_5$	0.8	1.2
$x_6$	-0.12	-0.18
$x_7$	0.8	1.2
$x_8$	0.8	1.2
$x_9$	0.8	1.2
$x_{10}$	0.8	1.2

Equations

$$\mathcal{A}\subset {\mathbb{R}}^2,$$

$$s(\ell ;\{a_j\})=\sum _{j=1}^k{a}_j{\varphi }_j(\ell )\text{where}{a}_j\in \mathbb{R},\ell \in [0,1].$$

$$s_U(\ell )=\tau {\mathbf{a}}^T\mathbf{\varphi }(\ell ),{\varphi }_j(\ell )=\{\begin{array}{ll}\sqrt{\ell },& j=1\\ {\ell }^{j-1},& j=2,\dots ,5\end{array},\tau >0.$$

$$s_U(\ell )\approx \tau a_1\sqrt{\ell }\ell \in [0,ϵ).$$

$${\varphi }_j(t)=\{\begin{array}{ll}t,& j=1\\ t^{2(j-1)},& j=2,\dots ,5\end{array}.$$

$$\frac{d}{d\ell }s(\ell )=\frac{dt}{d\ell }\frac{d}{dt}s(t)⟹\frac{d}{d\ell }s(\ell )=\frac{1}{2\sqrt{\ell }}\frac{d}{dt}s(t).$$

$$s_U(\ell )=(x_1\mathbf{a}{)}^T\mathbf{\varphi }(\ell ),s_L(\ell )=-(x_2\mathbf{a}{)}^T\mathbf{\varphi }(\ell ).$$

$$\mathbf{x}⟼s_U,s_L⟼f,$$

$$s(\ell )=s(\ell ;\mathbf{a})=\sum _{j=1}^6{a}_j{\ell }^{j-1/2}={\mathbf{a}}^T\mathbf{\varphi }(\ell ),\ell \in [0,1],$$

$$s(t;\mathbf{a})=\sum _{j=1}^6{a}_j{t}^{2j-1},t\in [0,1].$$

$$\left[\begin{array}{cc}{\mathbf{M}}_U(x_1)& 0\\ 0& {\mathbf{M}}_L(x_2)\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_U(\mathbf{x})\\ {\mathbf{a}}_L(\mathbf{x})\end{array}\right]=\left[\begin{array}{c}{\mathbf{p}}_U(\mathbf{x})\\ {\mathbf{p}}_L(\mathbf{x})\end{array}\right],$$

$$s(\ell ;r_1,r_2,\mathbf{x})=c_f(\ell ;r_1,r_2)s_f(\ell ;\mathbf{x}),$$

$$c_f(\ell ,\mathbf{r})={\ell }^{r_1}(1-\ell {)}^{r_2},$$

$$s_f(\ell ,\mathbf{x})=\sum _{j=0}^{m-1}{x}_j{\ell }^j,\mathbf{x}=[x_0,\dots ,x_{m-1}{]}^T.$$

$$s(t;\mathbf{x})$$

$$=x_{0}t-x_{m-1}{t}^{2m+1}+\sum _{j=1}^{m-1}(x_{j+1}-x_j)t^{2j+1}.$$

$$\int \mathbf{x}\rho (\mathbf{x})d\mathbf{x}=0,\int \mathbf{x}{\mathbf{x}}^T\rho (\mathbf{x})d\mathbf{x}=\mathbf{I}.$$

$$\mathbf{C}=\int \nabla f(\mathbf{x})\nabla f(\mathbf{x}{)}^T\rho (\mathbf{x})d\mathbf{x}=\mathbf{W}\Lambda {\mathbf{W}}^T\in {\mathbb{R}}^{m\times m},$$

$$\mathbf{x}={\mathbf{W}}_1\mathbf{y}+{\mathbf{W}}_2\mathbf{z}\mathbf{y}\in {\mathbb{R}}^n,\mathbf{z}\in {\mathbb{R}}^{m-n}.$$

$$f(\mathbf{x})\approx g({\mathbf{W}}_1^T\mathbf{x}),$$

$$f(\mathbf{x})\approx \frac{1}{2}{\mathbf{x}}^T\mathbf{H}\mathbf{x}+{\mathbf{x}}^T\mathbf{v}+c,$$

$$\nabla f(\mathbf{x})\approx \mathbf{H}\mathbf{x}+\mathbf{v}.$$

$$\begin{gathered} \mathbf{C}\approx \int \left(\mathbf{H}\mathbf{x}+\mathbf{v}\right){\left(\mathbf{H}\mathbf{x}+\mathbf{v}\right)}^T\rho (\mathbf{x})d\mathbf{x} \\ ={\mathbf{H}}^2+\mathbf{v}{\mathbf{v}}^T. \end{gathered}$$

$${\mathbf{H}}^2+\mathbf{v}{\mathbf{v}}^T=\hat{\mathbf{W}}\hat{\Lambda }{\hat{\mathbf{W}}}^T.$$

$${\hat{e}}_n^{(k)}={\Vert \left({\hat{\mathbf{W}}}_1^{(k)}\right){\left({\hat{\mathbf{W}}}_1^{(k)}\right)}^T-\left({\hat{\mathbf{W}}}_1\right){\left({\hat{\mathbf{W}}}_1\right)}^T\Vert }_2,$$

$$P_{\mathbf{y}}=\left\{\mathbf{y}\in {\mathbb{R}}^2\mid \mathbf{y}=\gamma \left[\begin{array}{c}{y}_{1,\text{min}}\\ 0\end{array}\right]+(1-\gamma )\left[\begin{array}{c}0\\ y_{2,\text{min}}\end{array}\right],\gamma \in [0,1]\right\},$$

$$y_{1,\text{min}}=\underset{\mathbf{x}\in [-1,1{]}^{10}}{\mathrm{minimum}}{\hat{\mathbf{w}}}_1^T\mathbf{x},y_{2,\text{min}}=\underset{\mathbf{x}\in [-1,1{]}^{10}}{\mathrm{minimum}}{\hat{\mathbf{w}}}_2^T\mathbf{x}.$$

$$s({\ell }_{int})=\mathbf{\varphi }({\ell }_{int}{)}^T\mathbf{a}$$

$$s(1)=\mathbf{\varphi }(1{)}^T\mathbf{a}⟹1^T\mathbf{a}$$

$$0=\frac{ds({\ell }_{int})}{d\ell }={\left(\frac{d}{d\ell }\left[\mathbf{\varphi }({\ell }_{int})\right]\right)}^T\mathbf{a}$$

$$\frac{ds(1)}{d\ell }={\left(\frac{d}{d\ell }\left[\mathbf{\varphi }(1)\right]\right)}^T\mathbf{a}$$

Full text

Nomenclature

[TABLE]

1 Introduction

The approximation of transonic flow solutions modeled by Euler’s equations has remained a relevant test problem in aerospace engineering—specifically, the design and evaluation of airfoils in a two-dimensional inviscid compressible flow field. The numerical representation of airfoils in transonic flow regimes involves parametrized geometries that inform optimization procedures, as described by Vanderplaats[40]. Several applications—e.g., Garzon[16], Dow and Wang[13], and Lukaczyk et al.[28]—approximate flow field solutions with parametric representations of airfoil shapes to produce estimates of quantities of interest (QoI). The performance of a design is measured by studying changes in QoI (e.g., dimensionless force quantities like lift and drag or measures of efficiency) as shape parameters vary. The QoI are computed from sufficiently accurate flow field approximations. The flow field approximations depend on the parameterized airfoil, which defines the boundaries in the governing equations. Samareh[33] surveys parameterization techniques for high fidelity shape optimization. Hájek[19] also offers a concise taxonomy of parameterizations and optimization routines.

We seek to further understand the QoI dependence on airfoil shape parameters through transformations and approximations that inform algorithms, insights, and pedigree—and ultimately improve airfoil designs. Toward this end, we search for active subspaces in the map from airfoil shape parameters to design quantities of interest, lift and drag, in a two-dimensional transonic airfoil model. An active subspace is defined by a set of important directions in the space of input parameters; perturbing parameters along these important directions (when present) changes the outputs more, on average, than perturbing the parameters along unimportant directions (i.e., directions orthogonal to the important directions). Such structure, when present, can be exploited to reduce the input space dimension and subsequently reduce the cost of algorithmic approaches to the design problem. Our work builds on recent work exploiting active subspaces in shape optimization[28], vehicle design[30], jet nozzle design[2], and dimensional analysis[11], as well as a gradient-based dimension reduction approach from Berguin and Mavris[3].

The present work extends previous work as follows. We study two common shape parameterizations—PARSEC[36] and CST[24]—and we show how design-enabling insights derived from active subspaces are consistent across the two parameterizations. Moreover, we show how the active subspace reduces the dimension and enables simple visualization tools that lead to deep insights into the lift/drag Pareto front in a multi-objective optimization for the airfoil design problem.

2 Airfoil Design

2.1 An interpretation of airfoils

An airfoil is a subset of a two-dimensional Euclidean space,

[TABLE]

with a boundary, $\partial\mathcal{A}$, representing the surface of the airfoil. The airfoil surface can be partitioned into two shape functions, $s_{U}$ and $s_{L}$, representing the upper and lower surfaces, respectively. Each shape function maps a normalized spatial coordinate $\ell\in[0,1]$, running from leading to trailing edge, to its respective surface. The following are necessary conditions on $s_{U}(\ell)$ and $s_{L}(\ell)$ such that they represent an airfoil for the present study:

Bounded: $0\leq s_{U}(\ell)\leq U<\infty$ for an upper bound $U$ and $-\infty<L\leq s_{L}(\ell)\leq 0$ for a lower bound $L$. 2. 2.

Smooth: continuous and twice differentiable for $\ell\in(0,1)$. 3. 3.

Fixed endpoints: $s_{U}(0)=s_{L}(0)=0$ and, optionally, $s_{U}(1)=s_{L}(1)=0$. 4. 4.

Feasible: $s_{U}(\ell)>s_{L}(\ell)$ for $\ell\in(0,1)$.

These conditions can be interpreted physically. Boundedness ensures that the resulting airfoil has a minimum and a maximum, i.e., the shape must remain finite to be manufacturable. Smoothness ensures that the interior of the shape is manufacturable to manufacturing precision. Traditionally, airfoils have a round leading edge, which produces an unbounded first derivative in $s_{U}$ and $s_{L}$ at $\ell=0$. The fixed endpoints ensure the two surfaces meet at $\ell=0$ and $\ell=1$. A sharp trailing edge is not manufacturable to finite precision, though it is used in most computational models. Feasibility ensures surfaces do not intersect. We ensure property four is satisfied by assuming positive shapes, $s_{U}(\ell)>0$ and $-s_{L}(\ell)>0$ for $\ell\in(0,1)$.

2.2 Airfoils for design

Airfoil shapes can be parameterized by (i) deformation of a nominal shape[42], (ii) smooth perturbations to a nominal shape[20], (iii) Karhunen-Loève expansions[16], (iv) Joukowski[19] or Theodorsen-Garrick conformal mappings[21], and (v) piece-wise spline interpolation with relevant bounding[37]. Alternatively, engineering intuition has motivated a variety of parameterized shapes, specifically those from the National Advisory Committee for Aeronautics (NACA)[22], Sobieczky[36], and Kulfan[24]. Considering the engineering insights, we study partitions of $\partial\mathcal{A}$ into shape functions that can be written as a linear combination of basis functions, $\phi_{j}:\mathbb{R}\rightarrow\mathbb{R}$, such that, for $s$ as one of $s_{U}$ or $s_{L}$,

[TABLE]

Varying the coefficients $a_{j}$ modifies the shape functions. Hence, the coefficients $\boldsymbol{a}=[a_{1},\dots,a_{k}]^{T}$ define the shapes with basis $\boldsymbol{\phi}(\ell)=[\phi_{1}(\ell),\dots,\phi_{k}(\ell)]^{T}$. Vector notation enables a concise expression for separate expansions in (1). For example, the top surface of a NACA airfoil can be expressed as[22]

[TABLE]

For NACA airfoils, the coefficients $\boldsymbol{a}$ sum to zero, which implies the fixed endpoint at the trailing edge. And $\boldsymbol{a}$ are constrained to ensure feasibility ($s_{U}>s_{L}$). The basis expansion is smooth on the open interval, bounded on the closed interval, and admits a leading edge ($\ell=0$) derivative singularity to model the roundness. Figure 1 shows the well-known NACA 0012 as two shape functions.

Consider $0<\epsilon\ll 1$ and $\ell\in[0,\epsilon)$ such that powers of $\ell$ are negligibly small. In this neighborhood of $\ell=0$, the upper surface can be approximated by

[TABLE]

That is, the symmetric NACA airfoil behaves like a multiple of the square root and admits the round-nose leading edge (consequently, a singularity at $\ell=0$). Figure 2 shows this behavior along with a $t=\sqrt{\ell}$ change-of-variables for perspective. Writing the NACA 0012 basis functions $\phi_{j}(\ell)$ from (2) with the change-of-variables,

[TABLE]

The transformation from $t$ to $\ell$ helps compare shape representations in the following sections.

Figure 2 demonstrates the utility of approximating shapes with the $\sqrt{\ell}$ change of variables. This transformation also appears in other shape function representations that have circular leading edge (i.e., round-nose) airfoils: Parametric Sections (PARSEC)[36] and Class-Shape Transformations (CST)[24], which we discuss in the following sections. These parametric representations admit polynomial series in $\sqrt{\ell}$ that achieve the round-nose shapes. Under the change of variables $t=\sqrt{\ell}$, the derivative singularities at $\ell=0$ for such series is verified by the chain rule,

[TABLE]

It is possible to eliminate the singularity in this shape partition through conformal mapping or additional partitioning of $\partial\mathcal{A}$. A classical representation involves the Joukowski conformal mapping of a circle but the parameterization admits a limited space of shapes based on two parameter values[1]. Designing inverse conformal mapping parameterizations can be difficult and existing conformal mappings are often associated with shapes admitting undesirable aerodynamic features[29]. However, the Theodorsen-Garrick conformal mapping offers a more generic shape representation from the inverse of near-circles[34, 1].

2.3 Parametric shape functions

Consider the following parameterization of the NACA airfoil, where two parameters, $x_{1}>0$ and $x_{2}>0$, independently scale the coefficients of the series (2),

[TABLE]

We think of (6) as a map from a two-dimensional parameter domain to airfoils, where the parameters $\boldsymbol{x}=[x_{1},x_{2}]^{T}$ control the thickness. Using this parameterization, we can model asymmetric shapes with different values for $\boldsymbol{x}$ (that is, $x_{1}\not=x_{2}$). Given a shape, we can compute a designer’s quantity of interest from a flow field approximated on the domain containing the airfoil; we denote the quantity of interest $f$. Consider the composition map,

[TABLE]

which suggests the following numerical procedure:

Choose parameters $\boldsymbol{x}$. 2. 2.

Compute parametric shape functions, $s_{U}(\ell)$ and $s_{L}(\ell)$. 3. 3.

Create a mesh on the set $\mathcal{A}^{c}$, and compute the CFD flow field. 4. 4.

Approximate the quantity of interest, $f$, from the CFD flow field solution.

Typically, asymmetric NACA airfoil parameterizations use an integer encoding for thickness and camber functions. The thickness and camber are combined using a polar coordinate transformation to add thickness orthogonal to a camber curve. The simpler NACA parameterization in (6) is meant to introduce the function, $f:\mathbb{R}^{2}\rightarrow\mathbb{R}$, that maps shape parameters $\boldsymbol{x}=[x_{1},x_{2}]$ to the quantity of interest.

Given (i) a shape parameterization with $m$ parameters $\boldsymbol{x}\in\mathbb{R}^{m}$ and (ii) a composition (7), a design optimization seeks $\boldsymbol{x}^{\ast}\in\mathbb{R}^{m}$ such that $f(\boldsymbol{x}^{\ast})\leq f(\boldsymbol{x})$ subject to constraints[40]. However, if the number $m$ of parameters is too large, then optimization algorithms may be computationally prohibitive[28]. If accurate approximations for all quantities of interest (i.e., objectives and constraints) can be obtained using a subspace of sufficiently small dimension (e.g., dimension less than or equal to two) in the shape parameter space, then visualization tools can aid the search for the optimum by providing insights into the relationship between inputs and outputs. In section 4, we discover two-dimensional active subspaces that enable these visualizations.

2.4 PARSEC airfoils

The PARSEC parameterization defines a set of shape functions; this parameterization was popularized by Sobieczky[27] for applications in transonic flow regimes. There have been several adaptations to PARSEC—e.g., Derksen and Rogalsky[12] adapt PARSEC for an optimization routine, while Wu et al.[42] focus on comparisons for a specific optimization problem. We focus on the original PARSEC parameterization for simplicity. The PARSEC expansion for a generic shape $s(\ell)$ is

[TABLE]

with coefficients $\boldsymbol{a}=[a_{1},\dots,a_{6}]^{T}$ and basis functions $\boldsymbol{\phi}(\ell)=[\ell^{1/2},\dots,\ell^{11/2}]^{T}$. Each surface $s_{U}$ and $s_{L}$ has its own independent coefficient vector, $\boldsymbol{a}_{U}$ and $\boldsymbol{a}_{L}$.

If we substitute the change-of-variables $t=\sqrt{\ell}$, the basis expansion becomes an odd-power polynomial in $t$,

[TABLE]

This change-of-variables enables precise comparison with the NACA 0012 parameterization (6) with basis (4). In particular, PARSEC can not exactly reproduce a NACA 0012 airfoil because of the difference in bases; PARSEC uses odd powers of $t$ while NACA 0012 includes even powers of $t$ with fixed, nonzero coefficients. Moreover, the change-of-variables gives insight into the ability of PARSEC to produce round-nose airfoils; see Figure 2.

The coefficients of both PARSEC surfaces, upper $\boldsymbol{a}_{U}$ and lower $\boldsymbol{a}_{L}$, solve a system of 12 parameterized linear equations. Each equation represents a physically interpretable characteristic of the shape, and the parameters affect the shape in physically meaningful ways. Table 1 gives a description of each of the 11 parameters, denoted $\boldsymbol{x}=[x_{1},\dots,x_{11}]^{T}$, that characterize the PARSEC shapes.

In PARSEC, bounds on the parameters $\boldsymbol{x}$ restrict size and shape variation for design and optimization problems. Certain parameters—namely $x_{5}$ (trailing edge offset), $x_{6}$ (trailing edge half-thickness), $x_{7}$ (wedge angle), $x_{8}$ (wedge half-angle), and $x_{11}$ (leading edge radius)—affect both upper and lower surfaces. Setting $x_{5}=x_{6}=0$ produces a sharp trailing edge.The parameterization contains two nonlinear relationships: (i) trailing edge slope and (ii) leading edge radius.

Given $\boldsymbol{x}$, the PARSEC coefficients $\boldsymbol{a}_{U}=\boldsymbol{a}_{U}(\boldsymbol{x})$ and $\boldsymbol{a}_{L}=\boldsymbol{a}_{L}(\boldsymbol{x})$ satisfy

[TABLE]

where each diagonal block of the matrix is 6-by-6 and depends on either $x_{1}$ or $x_{2}$ from Table 1. All parameters from Table 1, except $x_{1}$ and $x_{2}$, affect the right hand side $[\boldsymbol{p}_{U}(\boldsymbol{x}),\boldsymbol{p}_{L}(\boldsymbol{x})]^{T}$. For reasonable choices of $x_{1}$ and $x_{2}$, the matrix in (10) is invertible. Further details on the equality constraints are in Appendix A.

2.5 Class-shape transformation airfoils

A more recent airfoil parameterization is the Class-Shape Transformation (CST), which was popularized by Kulfan[24]. The shape function is decomposed into the product of a class function $c_{f}$ and a smooth shape functions $s_{f}$:

[TABLE]

where $r_{1}$ and $r_{2}$ are between 0 and 1, and the vector $\boldsymbol{x}\in\mathbb{R}^{m}$ contains shape parameters as polynomial coefficients. Each class function takes the form,

[TABLE]

where $r_{1}$ and $r_{2}$ control the shape at the leading and trailing edges, respectively. For simplicity, we fix $r_{1}=1/2$ and $r_{2}=1$, which produces round nose airfoils with sharp trailing edges. Alternative class function forms have been proposed and implemented in other applications[25, 31, 32]. The shape function is a polynomial of degree $m-1$ in $\ell$,

[TABLE]

With (i) our choices $r_{1}=1/2$ and $r_{2}=1$ and (ii) the change of variables $t=\sqrt{\ell}$, each CST surface can be written

[TABLE]

Compare (14) to the PARSEC series (8); for $m=5$ in (14), both are polynomials of degree 11 in $t$ containing only odd powers of $t$. However, the parameterizations differ dramatically. The PARSEC parameters correspond to physically meaningful characteristics of the airfoil; see Appendix A. In contrast, the parameters of CST directly impact the coefficients in the polynomial in $t$—as opposed to through a system of equations (10).

The domain for the composition (7) is the space of coefficients $\boldsymbol{x}=[\boldsymbol{x}_{U},\boldsymbol{x}_{L}]^{T}\in\mathbb{R}^{2m}$ that parameterize the upper and lower surfaces as in (13). We choose $m=5$ in the numerical experiments in section 4 to maintain the comparison between polynomial representations of CST and PARSEC; thus, the parameter space for CST is 10-dimensional. Given a point in the 10-dimensional parameter space, the CST shape is constructed via (13), (12), and (11). Then a CFD solver estimates the flow fields and produces the quantity of interest.

3 Active subspaces

We closely follow the development of active subspaces as in Constantine[5]. Active subspaces are part of an emerging set of tools for dimension reduction in the input space of a function of several variables. Given (i) a probability density function $\rho:\mathbb{R}^{m}\rightarrow\mathbb{R}_{+}$ that identifies important regions of the input space and (ii) a differentiable and square-integrable function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$. Assume, without loss of generality, that $\rho(\boldsymbol{x})$ is such that

[TABLE]

That is, we assume the space of $\boldsymbol{x}$ has been shifted and transformed such that $\boldsymbol{x}$ has mean zero and identity covariance. Consider the following symmetric and positive semidefinite matrix and its eigendecomposition,

[TABLE]

where $\Lambda=\text{diag}(\lambda_{1},\dots,\lambda_{m})$ and the eigenvalues are ordered from largest to smallest. Assume that for some $n<m$, $\lambda_{n}>\lambda_{n+1}$ (that is, the $n$th eigenvalue is strictly greater than the $n+1$th eigenvalue). Then we can partition the eigenvectors into the first $n$, denoted $\boldsymbol{W}_{1}$, and the remaining $m-n$, denoted $\boldsymbol{W}_{2}$. The active subspace is the column span of $\boldsymbol{W}_{1}$, and its orthogonal complement (the column span of $\boldsymbol{W}_{2}$) is the inactive subspace. Any $\boldsymbol{x}$ can be decomposed in terms of the active and inactive subspaces[5],

[TABLE]

We call $\boldsymbol{y}$ the active variables and $\boldsymbol{z}$ the inactive variables. Perturbing $\boldsymbol{y}$ changes $f$ more, on average, than perturbing $\boldsymbol{z}$; see Constantine et al.[6]. If the eigenvalues $\lambda_{n+1},\dots,\lambda_{m}$ are sufficiently small, then we can approximate

[TABLE]

for some $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$. In effect, this approximation enables design studies to proceed on the $n$ active variables $\boldsymbol{y}=\boldsymbol{W}_{1}^{T}\boldsymbol{x}$ instead of the $m$-dimensional design variables $\boldsymbol{x}$—i.e., dimension reduction.

With basic terms defined, we remark on the given density $\rho(\boldsymbol{x})$. In real engineering applications, this density function is rarely if ever known completely; thus, the engineer must choose $\rho$ such that it models the known or desired variability in the inputs $\boldsymbol{x}$. In such cases, $\rho$ is a modeling choice, and different $\rho$’s produce different results (e.g., different active subspaces for the same $f$). The degree to which one choice of $\rho$ is better than another depends only on the suitability and appropriateness of each $\rho$ to represent the parameter variability in a particular problem. If data is available on the parameter variability, then the chosen $\rho$ should be consistent with that data. In the absence of data, the subjective choice of $\rho$ is informed by engineering intuition and design experience. There is no generic prescription in the active subspace-based toolbox for choosing $\rho$; the only constraint is that $\rho$ must be a probability density. From there, the active subspace depends on the chosen $\rho$.

The active subspace also depends on the function $f$. In general, two different functions (e.g., two different quantities of interest) contain different active subspaces, which distinguishes active subspaces from covariance-based dimension reduction techniques such as principal component analysis[23] (PCA) that only depend on $\rho$. Principal components can define a shape parameterization given airfoil measurements[16], which yields the left side of the composite map (7). Thus, it is possible to use a PCA-based shape parameterization, combined with a differentiable scalar-valued output of interest, to produce well defined active subspaces. The scalar-valued distinction is important; other dimension reduction techniques, such as the proper orthogonal decomposition[41] (POD), seek to reduce the dimension of the state space of a large dynamical system or discretized PDE system. Although both active subspaces and POD use eigenpairs to define low-dimensional linear subspaces, these subspaces relate to different elements of a model-based design problem—active subspaces on the input parameters and POD on the state vector of the discretized flow field solutions.

3.1 Global quadratic model for estimating active subspaces

When the gradient $\nabla f(\boldsymbol{x})$ is available as a subroutine (e.g., through an adjoint solver or algorithmic differentiation[18]), one can use Monte Carlo or other numerical integration methods to estimate $\boldsymbol{C}$ from (16)[9]; the eigenpairs of the numerically estimated $\boldsymbol{C}$ estimate those of the true $\boldsymbol{C}$. However, many legacy solvers do not have such gradient routines. Finite differences build a local linear model at the point of the gradient evaluation and use the model’s slope as the approximate gradient. But finite differences are expensive when each evaluation of $f(\boldsymbol{x})$ is expensive and the dimension of $\boldsymbol{x}$ is large. We pursue a different model-based approach: (i) fit a global quadratic model with least-squares and (ii) estimate $\boldsymbol{C}$ with the quadratic model. This approach is motivated as follows. Suppose $f(\boldsymbol{x})$ is well approximated by a quadratic function,

[TABLE]

for some symmetric matrix $\boldsymbol{H}\in\mathbb{R}^{m\times m}$, $\boldsymbol{v}\in\mathbb{R}^{m}$, and $c\in\mathbb{R}$. The approximation is in the mean-squared sense with respect to the given density $\rho(\boldsymbol{x})$. Additionally, suppose the gradient $\nabla f(\boldsymbol{x})$ is well approximated as

[TABLE]

Again, the approximation is in the mean-squared sense. Substituting (20) into (16) and applying (15) results in the approximation

[TABLE]

The eigenpairs of $\boldsymbol{H}^{2}+\boldsymbol{v}\boldsymbol{v}^{T}$ estimate those of $\boldsymbol{C}$ in (16),

[TABLE]

The quality of the estimation depends on how well $f(\boldsymbol{x})$ and its gradient can be modeled by the global quadratic model (19). We emphasize that the quadratic model (19) is not used as a response surface for any design or optimization problem. Its sole use is to identify important directions in the parameter space and estimate the active subspace.

We fit the coefficients $\boldsymbol{H}$ and $\boldsymbol{v}$ of the quadratic model with a standard discrete least-squares approach as analyzed by Chkifa et al.[4]. That is, we (i) draw $N$ independent samples of $\boldsymbol{x}_{i}$ according to the given density $\rho$, (ii) evaluate $f_{i}=f(\boldsymbol{x}_{i})$ for each sample, and (iii) use least-squares with the pairs $(\boldsymbol{x}_{i},f_{i})$ to fit $\boldsymbol{H}$ and $\boldsymbol{v}$. There are ${m+2\choose 2}$ coefficients to fit, so we take $N$ to be 2-to-10 times ${m+2\choose 2}$ for the least-squares fitting, which is a common rule-of-thumb.

A related quadratic model-based approach appears in Tipireddy and Ghanem[39] when adapting a Hermite polynomial chaos basis for low-dimensional surrogate construction. In statistical regression, Li [26] proposes using a least-squares-fit regression model to estimate the principal Hessian directions, which are derived from eigenvectors of an averaged Hessian matrix of $f(\boldsymbol{x})$—in contrast to the expression in (21).

We treat the estimated eigenpairs $\hat{\boldsymbol{W}}$, $\hat{\Lambda}$ the same way we treat any other estimates of $\boldsymbol{C}$’s eigenpairs. That is, we examine the eigenvalues in $\hat{\Lambda}$. If a large gap exists between eigenvalues $\hat{\lambda}_{n}$ and $\hat{\lambda}_{n+1}$, then we choose $n$ to be the dimension of the active subspace. Then the first $n$ columns of $\hat{\boldsymbol{W}}$, denoted $\hat{\boldsymbol{W}}_{1}$, define the active subspace, which we can use to build a low-dimensional approximation to $f(\boldsymbol{x})$ as in (18).

We use a bootstrap-based[15] heuristic to assess the variability and uncertainty in the estimated eigenpairs from (22). We apply the bootstrap to the pairs $\{(\boldsymbol{x}_{i},f_{i})\}$ used to fit the quadratic model. For $k$ from 1 to $N_{\text{boot}}$,

let $\pi^{k}_{j}$ be an independent draw from the set $\{1,\dots,N\}$ for $j=1,\dots,N$, 2. 2.

use the pairs $\{(\boldsymbol{x}_{\pi_{j}^{k}},f_{\pi_{j}^{k}})\}$ to fit the coefficients $\boldsymbol{H}_{k}$ and $\boldsymbol{v}_{k}$ of a quadratic model with least squares, 3. 3.

estimate the eigenpairs $\hat{\boldsymbol{W}}^{(k)}$, $\hat{\Lambda}^{(k)}$ of the matrix $\boldsymbol{H}_{k}^{2}+\boldsymbol{v}_{k}\boldsymbol{v}_{k}^{T}$.

The eigenpairs $\hat{\boldsymbol{W}}^{(k)}k$, $\hat{\Lambda}^{(k)}$ constitute one bootstrap sample. We typically choose $N_{\text{boot}}$ to be between 100 and 1000, depending on the dimension of $\boldsymbol{x}$; if $\boldsymbol{x}$’s dimension is less than 100, then we can easily run thousands of eigenvalue decompositions from step 3 on a laptop. Note that the pairs in step 2 likely contain duplicates, which is essential in the sampling-with-replacement bootstrap method. Given the eigenvalue bootstrap samples $\{\hat{\Lambda}^{(k)}\}$, we use the range over bootstrap samples to estimate uncertainty in the eigenvalues $\hat{\Lambda}$. To estimate uncertainty in the $n$-dimensional active subspace, let $\hat{\boldsymbol{W}}_{1}^{(k)}$ contain the first $n$ columns of the bootstrap sample of eigenvectors $\hat{\boldsymbol{W}}^{(k)}$, and compute the subspace error[17],

[TABLE]

where the norm is the matrix 2-norm. The mean and standard error of the set $\{e_{n}^{(k)}\}$ estimates the error and uncertainty in the active subspace defined by $\hat{\boldsymbol{W}}_{1}$ from (22).

There are two important caveats to the bootstrap-based heuristic for assessing uncertainty in the computed quadratic model-based active subspace. First, the pairs $\{(\boldsymbol{x}_{i},f_{i})\}$ (i.e., the data) used to fit the quadratic model do not satisfy the assumptions needed to allow a proper statistical interpretation of the bootstrap-based error estimates. Proper interpretations require structured noise on the $f_{i}$’s; when $f_{i}$ is the output of a deterministic computer simulation, such noise models are not valid. Second, the bootstrap-based error estimates only quantify the error due to finite sampling when estimating the quadratic model. They do not include errors from the motivating quadratic model assumption (19) and (20). Despite these caveats, the bootstrap-based error estimates give some indication of variability in the eigenpairs estimated in (22).

In section 4, we show numerical evidence of convergence as $N$ increases for active subspaces derived from the composite map (7) with both PARSEC and CST as the underlying airfoil shape parameterizations.

4 Active subspaces in PARSEC and CST parameterizations

The goal of the following numerical experiments is to study how the active subspaces of lift and drag change as the shape parameterization changes. That is, if we change the first component of the composite map (7), $\boldsymbol{x}\rightarrow s_{U},s_{L}$, (i) how do the active subspaces change for the composite map and (ii) how do insights into the design problem change? In particular, we compare active subspaces derived using both PARSEC and CST parameterizations for the airfoil surfaces.

In each case, we must set the density function $\rho(\boldsymbol{x})$ needed to define the active subspaces. We want to ensure that all parameters $\boldsymbol{x}$ in the support of $\rho$ produce shapes that satisfy the necessary conditions from section 2.1; our approach is somewhat conservative. For each parameterization, we first find the parameters that produce a shape nearest the NACA 0012 using nonlinear least-squares, which defines the respective center of the parameter space. We then set bounds on the parameters to be $\pm 20$% of the center. The density $\rho$ is uniform on the hyperrectangle defined by these ranges.

Given a shape, we estimate the quantities of interest lift and drag at fixed $0^{\circ}$ angle of attack and NIST conditions (ambient temperature $0^{\circ}$ C, ambient pressure $101.325$ KPa, Mach $0.8$) using the Stanford University Unstructured (SU2) CFD solver[14], which solves Euler equations for compressible flow on an unstructured mesh. We executed each run on the Colorado School of Mines Mio cluster. The data supporting these results can be found on GitHub (\urlhttps://github.com/zgrey/ASAP.git).

4.1 PARSEC shape parameterization results

The bounds for the PARSEC parameters (see Table 1) are in Table 2. These bounds are $\pm 20$% of fitted values for a NACA 0012 airfoil. To generate the pairs needed for the quadratic model-based approximation of the active subspaces, we draw $N=6500$ independent samples from the uniform density on the hyperrectangle defined by the ranges in Table 2. Each sample produces a shape that leads to a CFD solve and lift and drag coefficients, $C_{l}$ and $C_{d}$, respectively. The pairs $\{(\boldsymbol{x}_{i},C_{l}(\boldsymbol{x}_{i}))\}$ are used to estimate the active subspaces for lift as a function of PARSEC parameters following the procedure outlined in section 3.1. And the pairs $\{(\boldsymbol{x}_{i},C_{d}(\boldsymbol{x}_{i}))\}$ are similarly used to estimate the active subspaces for drag.

Figure 3 shows the 11 eigenvalues $\hat{\Lambda}$ from (22) for lift and drag on a logarithmic scale along with the bootstrap ranges that estimate uncertainty (see section 3.1); we use $N_{\text{boot}}=500$ bootstrap samples. For lift, there is a large gap between the first and second eigenvalues, which is evidence of an exploitable one-dimensional active subspace. For drag, a relatively large gap occurs between the second and third eigenvalues, which suggests an exploitable two-dimensional active subspace. Figure 4 shows, for both lift and drag, the means and ranges for the bootstrap-based subspace error estimates (23) as a function of the active subspace dimension. Notice that the estimated error increases as the subspace dimension increases. This is related to the decreasing gaps between subsequent eigenvalues from Figure 3. As the gap between two eigenvalues decreases, estimating the eigenspace whose dimension is associated with the smaller eigenvalue becomes more difficult for numerical procedures[38]. In the extreme case, if $\lambda_{n}=\lambda_{n+1}$, then there is no eigenspace of dimension $n$.

Figure 5 shows components of the eigenvectors $\hat{\boldsymbol{W}}$ from (22). Figure 5a shows the first two eigenvectors associated with lift. And Figure 5b shows the first two eigenvectors associated with drag. These eigenvector components can be used to compute sensitivity metrics for the PARSEC parameters with respect to the outputs of interest[7]. The indices on the horizontal axis labels correspond to the labels in Table 1. A large eigenvector component implies that output is more sensitive to changes in the corresponding parameter. The eigenvector components suggest that lift is most influenced on average by $x_{3}$, $x_{4}$, $x_{5}$, and $x_{6}$, which affect the upper and lower surface coordinates and trailing edge deflections. Additionally, the eigenvector components for upper and lower surface coordinates are nearly equal and opposite ($\pm 0.5$). This suggests that increasing the upper surface coordinate should be combined with reducing the lower surface coordinate to change lift the most. These coordinate changes lead to a change in asymmetry of the shape (camber). The eigenvector components for drag (Figure 5b) are largest for $x_{3}$ and $x_{4}$, corresponding to the upper and lower surface coordinates, respectively. These large components have the same sign, which suggests that simultaneous changes in the upper and lower surface coordinate affect drag more, on average. Simultaneously perturbing $x_{3}$ and $x_{4}$ changes thickness. This agrees with the design intuition that changes in thickness should influence drag the most, on average. Perturbing $\boldsymbol{x}$ along the second eigenvector for drag changes camber.

Given (i) the eigenvectors $\hat{\boldsymbol{W}}$ from (22) and (ii) the lift data $\{(\boldsymbol{x}_{i},C_{l}(\boldsymbol{x}_{i}))\}$ and drag data $\{(\boldsymbol{x}_{i},C_{d}(\boldsymbol{x}_{i}))\}$, we can create useful visualizations called shadow plots that elucidate the relationship between the shape parameters $\boldsymbol{x}$ and the quantities of interest $C_{l}$ and $C_{d}$. These shadow plots are motivated by sufficient summary plots developed by Cook[10] for sufficient dimension reduction. Given a data set and basis vectors defining a subspace of the inputs, the shadow plot’s construction is identical to the sufficient summary plot—but their interpretations differ. Proper interpretation of the sufficient summary plot relies on the deep theory of sufficient dimension reduction in statistical regression that justifies the label sufficient for the plot. In contrast, the shadow plot is a qualitative tool with no notion of statistical sufficiency. Interpreting the shadow plot requires the modeler’s subjective judgment. See Chapter 1 in Constantine[5] for more information on these plots using data derived from computational science simulation models.

Let $\hat{\boldsymbol{w}}_{1}$ be the first column of $\hat{\boldsymbol{W}}$. The one-dimensional shadow plot plots the output (e.g., $C_{l}(\boldsymbol{x}_{i})$) on the vertical axis versus a linear combination of the associated input vector $\hat{\boldsymbol{w}}_{1}^{T}\boldsymbol{x}_{i}$ on the horizontal axis, where the weights of the linear combination are the components of the vector $\hat{\boldsymbol{w}}_{1}$. If a one-dimensional summary plot shows a univariate relationship between the linear combination of inputs $\hat{\boldsymbol{w}}_{1}^{T}\boldsymbol{x}$ and the output, then an approximation of the form (18) with one vector is appropriate. If a univariate relationship is not apparent, then one can make a two-dimensional shadow plot. Let $\hat{\boldsymbol{w}}_{2}$ be the second column of $\hat{\boldsymbol{W}}$. The two-dimensional shadow plot puts the outputs on the $z$-axis, the linear combination of inputs $\hat{\boldsymbol{w}}_{1}^{T}\boldsymbol{x}_{i}$ on the $x$-axis, and the linear combination of inputs $\hat{\boldsymbol{w}}_{2}^{T}\boldsymbol{x}_{i}$ on the $y$-axis. If the output appears to be a function of these two linear combinations, then we can approximate the output as in (18) with two vectors. If the two-dimensional summary plot does not reveal a useful bivariate relationship, then the approximation (18) may be appropriate with more than two vectors, but metrics other than the shadow plots are needed to assess the approximation quality (e.g., residuals).

Figure 6 shows the one-dimensional shadow plots for lift and drag using the PARSEC parameterization with all 6500 samples used to estimate the quadratic model-based active subspaces. Note the strong monotonic relationship in both cases. Such a relationship may be useful if one wishes to estimate the range of outputs over the range of inputs; see Constantine et al.[8] for an example of this in a hypersonic scramjet model. However, the relationship is far from univariate. Figure 7 shows the two-dimensional shadow plots for the lift and drag. The bivariate relationship for drag is much more apparent.

4.2 CST shape parameterization results

We follow the same procedure to set $\rho$ for the CST parameterization. First find values for the CST parameters such that the resulting shape matches a NACA 0012 airfoil. Then choose bounds on the parameters to be $\pm 20$% of the fitted values, and set $\rho$ to be a uniform density on the resulting hyper-rectangle. Table 3 shows the bounds for the CST parameters. We draw $N=6500$ samples of the $m=10$ CST parameters independently according to $\rho$. For each parameter sample, we construct the airfoil, solve for the flow fields with SU2[14], and compute the outputs of interest (lift and drag). The resulting pairs are used to estimate the quadratic model-based active subspaces for lift and drag as described in section 3.1.

Figure 8 shows the eigenvalues from (22) with bootstrap ranges for lift and drag. In contrast to the eigenvalues using the PARSEC parameterization from Figure 3, the eigenvalues using CST show very pronounced gaps (i) between the first and second eigenvalues for lift and (ii) between the second and third eigenvalues for drag. These gaps suggest that a one-dimensional active subspace is useful for approximating lift, and a two-dimensional active subspace is useful for approximating drag. These conclusions are consistent with the conclusions we drew from the PARSEC-parameterized results. Figure 9 shows the bootstrap-based estimated subspace errors due to finite sampling as the dimension of the subspace increases. Similar to Figure 4, the subspace errors are inversely proportional to the associated eigenvalue gaps in Figure 3, which is consistent with well-known theory for numerical approximation of eigenspaces[38]. Thus, the errors in the estimated one-dimensional active subspace for lift and the one- and two-dimensional active subspaces for drag are small enough to be useful for design optimization.

Figure 10 shows (i) the first eigenvector from (22) for lift and (ii) the first two eigenvectors for drag using the CST shape parameterization. The nonzero eigenvector components correspond to the leading coefficient in the CST series (14) for each of the upper and lower surfaces—index 1 for the upper surface and index 6 for the lower surface. The eigenvector components suggest that the active subspaces are defined purely by these leading coefficients—both for lift and drag. In other words, CST parameters $x_{2},\dots,x_{5},x_{7},\dots,x_{10}$ have no impact on either lift or drag (for the chosen $\rho(\boldsymbol{x})$). Additionally, the eigenvector for lift is identical to the second eigenvector for drag.

Figure 11 shows (i) the one-dimensional shadow plot for lift and (ii) the two-dimensional shadow plot for drag. Lift is very well approximated by a linear function of one linear combination of the 10 CST parameters, where the weights of the linear combination are the components of the eigenvector in Figure 9a. This suggests that lift can be approximated as in (18) using the first eigenvector from (22) and a linear function for $g$. Drag is very well approximated by a function of two linear combinations of the 10 CST inputs; the linear combination weights are the components of the eigenvectors in Figure 9b. In this case, drag can be approximated as in (18) using the first two eigenvectors from (22) and a quadratic function for $g$. The fact that the relationship between the active variables and drag appears quadratic is independent from the quadratic model used to estimate the active subspaces. We emphasize that the drag values in Figure 11 are computed with the SU2 flow solver—not the quadratic response surface used to estimate the eigenvectors from (22).

Since the eigenvectors’ nonzero components correspond to the leading coefficients in the CST series (14), the active subspaces are directly related to the circular leading edge radii; see Figure 2. According to the eigenvector defining lift’s one-dimensional active subspace (Figure 10), the subspace relates to shape perturbations from equal and opposite changes in the upper and lower surface radii. That is, positive perturbations to the active variable increase the upper surface radius and decrease the lower surface radius by the same magnitude. This perturbation agrees with the intuition that changes in camber modify lift the most on average. Additionally, the first eigenvector for drag indicates exactly equal perturbations to the leading edge radii (i.e., positive perturbations to the first active variable increase the upper and lower surface radii equally). The resulting perturbation agrees with the intuition that changes in thickness modify drag the most on average. However, the PARSEC parameterization subspaces both indicated non-zero weights corresponding to trailing edge perturbations of the airfoil shape which were omitted from the CST parameterization. Therefore, repeating the CST study to incorporate trailing edge perturbations may reveal subspaces different from those defined by the eigenvectors in Figure 10.

4.3 Convergence of active subspaces with number of samples

The conclusions about the input/output relationships drawn from the active subspace analysis are dependent upon the number $N$ of samples used to estimate the active subspaces (in this case, via the least-squares fit global quadratic model). We study the convergence of the subspace errors as $N$ increases. Figure 12 shows convergence for both parameterizations using doubled sample sizes from 100 up to 6400.

The intervals in Figure 12 are the maximum and minimum from the bootstrap routine in 3.1. The point values were computed as averages over the bootstrap samples. For increasing sample size, $N$, there is approximately $1/\sqrt{N}$ convergence using the global quadratic approximation of the subspaces. The subspace errors for all doubled sample sizes were consistently less with CST compared to PARSEC. Therefore, for the same number of samples, the CST parameterization resulted in more accurate subspace estimates, on average.

4.4 Multi-objective design interpretation

For the purpose of design, we want to exploit the clear univariate and bivariate relationships between the active variables (linear combinations of the parameters) and the outputs present using the CST parameterization (see Figure 11). Leveraging the dimension reduction, we can visualize solutions to a multi-objective optimization problem. Assume the goal is to minimize drag and maximize lift. Recall that the second eigenvector for drag was the same as the eigenvector for lift; see Figure 10. Therefore, a two-dimensional shadow plot of lift versus the linear combinations defined by the drag eigenvectors shows variation in lift along the second active variable only; see Figure 13.

Figure 13 shows the two-dimensional shadow plot of lift using the two eigenvectors that define drag’s two-dimensional active subspace. Since lift’s one-dimensional active subspace is also defined by drag’s second eigenvector, the shadow plot shows that lift varies linearly along the second drag active variable. The contours overlaying the lift shadow plot show drag as a function of its active variables. The perceived relationship between the lift shadow and the drag contours immediately reveals the Pareto front for the multi-objective design problem[40]. In particular, the lower left edge of the shadow plot contains the set of drag active variables $\boldsymbol{y}=(\hat{\boldsymbol{w}}_{1}^{T}\boldsymbol{x},\hat{\boldsymbol{w}}_{2}^{T}\boldsymbol{x})$, constrained by the parameter space from Table 3, that either (i) maximize lift for fixed drag or (ii) minimize drag for fixed lift. Let $P_{\boldsymbol{y}}$ denote this set of active variables,

[TABLE]

where, for drag eigenvectors $\hat{\boldsymbol{w}}_{1}$ and $\hat{\boldsymbol{w}}_{2}$,

[TABLE]

The hypercube parameter space $[-1,1]^{10}$ is a shifted and scaled version of the hyperrectangle defined by the ranges in Table 3.

We can visualize the lift/drag Pareto front using response surfaces for lift and drag as functions of the CST parameters. The shadow plots in Figure 11 suggest simple response surfaces of the form (18) for each: (i) lift is well approximated by a linear function of its first active variable and (ii) drag is well approximated by a quadratic function of its first two active variables. The blue line in Figure 14 shows lift versus drag—computed with response surfaces—for all active variables in $P_{\boldsymbol{y}}$ from (24).

Note that for any value of the active variables $\boldsymbol{y}$, there are infinitely many $\boldsymbol{x}$ in the parameter space—and consequently infinitely many CST shapes—that map to the same $\boldsymbol{y}$. Each $\boldsymbol{x}$ that maps to a fixed $\boldsymbol{y}$ has different values of the inactive variables $\boldsymbol{z}$; see (17). However, by construction, we expect that perturbing the inactive variables does not change the outputs (lift and drag) very much. We verify this insensitivity by choosing several values for $\boldsymbol{z}$ for each $\boldsymbol{y}$ such that the resulting $\boldsymbol{x}$ are in the parameter space defined by the ranges in Table 3. Each $\boldsymbol{y},\boldsymbol{z}$ combination produces a unique airfoil shape. In Figure 14, several SU2-computed lift/drag combinations are shown; the CST shape for each combination has a parameter vector $\boldsymbol{x}$ whose active subspace coordinates $\boldsymbol{y}=\hat{\boldsymbol{W}}_{1}^{T}\boldsymbol{x}\in\mathbb{R}^{2}$ are in $P_{\boldsymbol{y}}$ from (24). The circles correspond to a particular choice of $\boldsymbol{z}=0$, and the asterisks correspond to random $\boldsymbol{z}$ values such that $\boldsymbol{x}=\hat{\boldsymbol{W}}_{1}\boldsymbol{y}+\hat{\boldsymbol{W}}_{2}\boldsymbol{z}$ is in the (scaled) parameter space. Note that all SU2-computed lift/drag combinations are remarkably close to the response surface-generated Pareto front. This verifies that (i) the inactive variables for lift and drag are indeed relatively inactive and (ii) the Pareto front estimated with quadratic model-based active subspaces and active subspace-exploiting low-dimensional response surfaces for lift and drag is accurate. We emphasize that this conclusion depends on several assumptions including the density $\rho$ derived from Table 3 and the angle-of-attack and NIST conditions for the flow field simulations. If any of these assumptions change, then the numerical study should be repeated.

In general, two nearby points on the Pareto front do not have similar input parameters or airfoil shapes[40]. There is no notion of continuity among the shapes whose lift/drag combinations define the Pareto front. Standard surrogate-based approaches[35] offer no means to reveal continuous perturbations to airfoil shapes that remain on the Pareto front. With active subspaces, we can define a conditional Pareto front as a function of the active variables by averaging over the inactive variables $\boldsymbol{z}$ or choosing a fixed $\boldsymbol{z}$ of interest. Then airfoil shapes are continuous with changes in the active variables. This continuity in the shape is emphasized in Figure 15 as five flow visualizations corresponding to the $P_{i}$ points in Figure 14.

For example, taking $\boldsymbol{z}=0$, the values $\boldsymbol{y}$ in the Pareto subset $P_{\boldsymbol{y}}$ from (24)—i.e., active variables along the lower left edge of the shadow plot in Figure 13—produce continuous changes to the airfoil shape parameters and consequently continuous changes in the airfoil shape. We expect this continuity, when enabled by the active subspaces, would be very helpful for airfoil designers.

5 Summary and conclusions

We have examined two transonic airfoil shape parameterizations. PARSEC uses a set of parameterized equality constraints to define the coefficients of a series that defines the upper and lower airfoil surfaces; each parameter affects a physical characteristic of the surfaces. CST’s parameters define the series’ coefficients directly, constrained to satisfy leading and trailing edge conditions. For each parameterization, we estimate active subspaces for airfoil performance metrics lift and drag—computed with CFD—as functions of the airfoil parameters. To estimate the active subspaces, we employ a global least-squares-fit quadratic approximation. The active subspaces derived from this analysis reveal one-dimensional structure in lift and two-dimensional structure in drag for each shape parameterization—the 11-parameter PARSEC and the 10-parameter CST. The vectors defining the active subspaces give insight into the important shape parameters for each performance metric, and these insights agree with first principles aerodynamics in both parameterizations. Moreover, the low-dimensional active subspaces enable simple visualization tools called shadow plots that aid in response surface construction. Additionally, the visualizations and low-dimensional active subspaces reveal insights into the Pareto front of a multi-objective design problem using both lift and drag as objective functions. For example, we can describe continuous perturbations to the airfoil that move smoothly along the Pareto front.

Appendix A Appendix

The appendix details the PARSEC constraints for a single shape function and an alternative leading edge radius relationship proof.

Constraints:

For a single interior interpolation point $\ell_{int}\in(0,1)$:

[TABLE] 2. 2.

Trailing edge point interpolation $\ell=1$:

[TABLE] 3. 3.

Interior interpolation point slope of zero $\ell_{int}\in(0,1)$:

[TABLE] 4. 4.

Trailing edge slope $\ell=1$:

[TABLE] 5. 5.

Interior point inflection at $\ell_{int}\in(0,1)$:

[TABLE] 6. 6.

Intermediate scale circular approximation $0<\epsilon\ll 1$, $\ell\in(0,\epsilon)$:

[TABLE]

For the six constraints, we construct a linear system to represent the separate shape function equality constraints

[TABLE]

and invertible

[TABLE]

Proof of leading coefficient and leading edge radius dependency:

Proof.

Determining the appropriate coefficient for the sixth constraint requires us to consider intermediate scale $t$. Consider a change of variables, $t=\sqrt{\ell}$, and recall the odd-power polynomial series expansion as the basis for PARSEC with six terms and CST with arbitrary basis,

[TABLE]

Additionally, construct a leading edge circular approximation with radius $\epsilon$ at a distance $\epsilon$ from the origin such that

[TABLE]

Substituting the intermediate scale change-of-variables and solving for $s_{c}^{2}(t)$ results in

[TABLE]

Truncate the odd polynomial basis expansion for $k>2$ and square the expression such that

[TABLE]

For $0<\epsilon\ll 1$, assume $t^{2}\in(0,\epsilon]$ (i.e., $\ell\in(0,\epsilon]$). Truncating powers of $t$ greater than $2$ results in the small scale approximations of the circular function and expansion,

[TABLE]

Equating the two approximations and solving for the coefficient when $t\neq 0$ implies

[TABLE]

The positive sign corresponds to the upper surface coefficient while the negative sign corresponds to the lower surface coefficient. ∎