Shape Sensitivity Analysis in Flow Models Using a Finite-Difference Approach

Abstract

Reduced-order models have a number of practical engineering applications for unsteady flows that require either low-dimensional approximations for analysis and control or repeated simulation over a range of parameter values. The standard method for building reduced-order models uses the proper orthogonal decomposition (POD) and Galerkin projection. However, this standard method may be inaccurate when used “off-design” (at parameter values not used to generate the POD). This phenomena is exaggerated when parameter values describe the shape of the flow domain since slight changes in shape can have a significant influence on the flow field. In this paper, we investigate the use of POD sensitivity vectors to improve the accuracy and dynamical system properties of the reduced-order models to problems with shape parameters. To carry out this study, we consider flows past an elliptic cylinder with varying thickness ratios. Shape sensitivities (derivatives of flow variables with respect to thickness ratio) computed by finite-difference approximations are used to compute the POD sensitivity vectors. Numerical studies test the accuracy of the new bases to represent flow solutions over a range of parameter values.

1. Introduction

Reduced-order modeling of incompressible flows plays an important role in academic and industrial research. In order to reduce the complexity of the governing equation, reduced-order models are often developed to represent the dynamical system with a few degrees of freedom. These models can provide analytical insight into the physical phenomena and enable application of dynamical systems theory and control methods.

Since Roshko [1] measured the vortex shedding period behind a bluff body, many researchers have investigated this phenomenon experimentally and numerically for a wide range of Reynolds numbers. The most frequently investigated bluff geometry is the circular cylinder. The flow behind a circular cylinder has become the canonical problem for studying such external separated flows [2–7]. Bishop and Hassan [2] were among the first to suggest a phenomenological model reduction by proposing Rayleigh/van der Pol self-excited oscillators to represent the forces over a cylinder due to vortex shedding. Several other analytical models have also been proposed along similar lines for stationary or moving cylinders [8–13].

On the other hand, engineering applications often involve flows over complex bodies, such as wings, submarines, missiles, and rotor blades, which can hardly be modeled as a circular cylinder. It would however be interesting to broaden the application of reduced-order models to other complex geometries. A systematic approach would be to start by studying a configuration that is more general than a circular cylinder and can characterize typical engineering flow configurations. Elliptic cylinders seem to provide such a configuration and changes in the eccentricity allow for shapes ranging from a circular cylinder to a flat plate. By tuning the damping coefficients, self-excited oscillator models [14] can be extended to elliptic cylinders.

In some cases, the reduced-order models represent the space variation with a limited number of basis functions, while still capturing the physics and dominant features of the flow. Most widely used reduced-order model techniques in fluid dynamics are derived from the proper orthogonal decomposition- (POD-) Galerkin projection approach [15, 16]. The POD provides a tool to formulate an optimal basis (or modes) to represent given trajectories of a dynamical system. The most common approach is to compute a “representative” collection of states of the system (commonly referred to as snapshots) through numerical simulations or experiments, then use the POD to find a low-dimensional set of basis functions [16, 17]. The POD processes time snapshot data to find the best basis functions on which to represent this data. Later, a low-dimensional dynamical system is typically obtained by Galerkin projection. These models give insight to the flow physics, often reproduce the data, and may be used for control purposes. POD-based reduced-order models have been successfully applied to model vortex shedding past a cylinder [18–23]. In fact, the flow past a cylinder has become a benchmark problem for reduced-order modeling.

Despite the accurate reproduction of the data from which it is originated, the POD-based reduced-order model lacks robustness away from the reference simulations. This is a serious limitation of the POD-based reduced-order models since most applications use reduced-order models in a predictive setting. In [24], the authors noted that the accuracy of the model predictions rapidly deteriorates as we move away from the decomposition value. Thus, the POD basis generated from an experiment or a numerical simulation is only useful within a narrow bandwidth of parameters close to that reference. In addition to the accuracy of the model away from the reference, numerical stability of these models for long-time integration is also an issue [23, 25–30].

Various methods have been proposed to extend the applicability and increase the accuracy of POD-based reduced-order models away from the reference data. Ma and Karniadakis [19] investigated the stability and dynamics of three-dimensional limit-cycle states near the Mode A [6] instability for flow past a cylinder using a low-dimensional modeling. They employed two procedures to construct this hybrid system, either by concatenating the two sets of snapshots at Re=182 and Re=185 or by extracting the POD modes at each state and orthonormalizing the entire set. They concluded that reference data must contain sufficient information about the flow field to predict any qualitative change within that flow regime. Noack et al. [20] suggested a generalization for POD-based Galerkin models to include the transient behavior. In their Galerkin approximation, they included an additional vector and termed it a shift-mode (also known as a centering trajectory of a mean-field correction). They obtained the steady solution with a Newton iteration, employing the discretized steady Navier-Stokes equations. The shift mode represents the shift of the short term averaged flow away from the POD space and is effectively a normalized mean-field correction. Within the control setting of POD-based reduced-order models, a similar strategy is required in cylinder wake problems to make the control effective.

In the application of reduced-order models to optimal control, one desires a model that is accurate over a wide range of actuation data. In a recent study, Graham et al. [21] used reduced-order models to develop optimal control strategies for flow past a circular cylinder at Re=100 through cylinder rotation. When the POD modes were generated from the snapshots at a given oscillation frequency, the active control model only worked near that particular frequency. Therefore, they incorporated multiple frequency information in their snapshot data set using a varying frequency sinusoid, or chirp. With this approach, one may require a large parametric data set as the reference data to broaden the applicability of POD-based reduced-order models over a wide range of parameter values.

In our previous work [31], we proposed the inclusion of parametric derivatives of the POD basis functions (computed using sensitivity analysis). The study considered two approaches for using the POD derivatives to define a reduced basis. One approach used linear extrapolation of the standard POD basis to the parameter value of interest (extrapolated basis). A second approach combined with the POD and POD derivatives to construct the basis (expanded basis). Both approaches provided promising results in the study of flow past a square cylinder where the viscosity (equivalently, the Reynolds number) was varied from the reference value.

A follow-up study [32] considered the development of reduced-order models for flows where the parameter described the orientation of the square cylinder relative to the free-stream flow direction. We demonstrated our shape sensitivity analysis approach by rotating the square cylinder. This problem also featured a “hidden” parametric dependence on the Reynolds number as the characteristic length scale changes with the orientation of the square cylinder to the mean-stream flow.

In this study, we use a similar basis improvement methodology. However, we consider modeling the flow as the obstacle is changed from a circular cylinder to an elliptic cylinder. The shape parameter is chosen as the thickness ratio (𝜏) between the minor axis to the major axis. It is important to note that the Reynolds number is defined using the projected area of the cylinder, as seen by the flow, and is kept constant in our analysis. Moreover, we use finite differences to compute the shape sensitivity as the ellipse parameter 𝜏 is varied. We employ an O-grid for the cylinders such that the outer domain is circular in all cases. Since the outer domains of both the circular and elliptic cylinders are the same, the shape sensitivity approaches to zero as it moves away from the body. In this regard, the finite-difference-based methodology provides a quick and easier way to compute the shape sensitivity and extend the application of POD basis to actual shape deformation. This study can find its application in the aerospace industry where the wing cross-section changes its actual shape, such as a morphing wing. This study also provides a validation to the basis selection approach used by Ma and Karniadakis [19]. In their study, they built a POD basis using two simulations with similar parameters. Their approach produces a similar basis as we would obtain with our expanded approach.

The paper is organized as follows. Section 2 explains the numerical methodology to solve the Navier-Stokes equations and simulate the flow past circular and elliptic cylinders. In Section 3, we present the POD-based reduced-order model procedure after collecting the snapshots of the flow field. In Section 4 we present the sensitivity analysis procedure and numerical results are discussed in Section 5.

2. Numerical Methodology

The Navier-Stokes and continuity equations are the governing equations for the present problem. For incompressible flow, they can be represented as follows:

𝜕𝑢_𝑗𝜕𝑥_𝑗𝜕𝑢_𝑖=0,(2.1)𝜕𝑡+𝜕𝜕𝑥_𝑗𝑢_𝑗𝑢_𝑖=−1𝜌𝜕𝜕𝑥𝑝_𝑖+𝜇𝜕²𝑢_𝑖𝜕𝑥_𝑗𝜕𝑥_𝑗,(2.2) where 𝑖,𝑗=1,2 for two-dimensional flows; the 𝑢_𝑖 represent the Cartesian velocity components (𝑢,𝑣); 𝑝 is the pressure; 𝜌 is the fluid density; 𝜇 is the fluid viscosity.

We employ curvilinear coordinates in an Eulerian reference frame; a planar view is shown in Figure 1(a). Here, 𝑋 and 𝑌 indicate the Cartesian coordinates, 𝜉 and 𝜂 represent a curvilinear coordinate system, and 𝐿_𝑥 and 𝐿_𝑦 are the major and minor axes, respectively. The thickness ratio of the ellipse is given by 𝜏=𝐿_𝑥/𝐿_𝑦 and is the parameter used in our parametric modeling study.

(a)

(b)

Figure 1 

(a) Geometry of an elliptic cylinder and (b) “O” grid distributed among 8 processors in the $\eta$ direction for $\tau =2.0$.

The Navier-Stokes equations are written in curvilinear coordinates and strong conservation form as follows:

𝜕𝑈_𝑚𝜕𝜉_𝑚𝜕𝐽=0,(2.3)₋₁𝑢_𝑖𝜕𝑡+𝜕𝜕𝐹𝜉^𝑖𝑚_𝑚=0,(2.4) where the flux is defined as

𝐹_𝑖𝑚=𝑈_𝑚𝑢_𝑖+𝐽⁻¹𝜕𝜉_𝑚𝜕𝑥_𝑖𝑝−1Re𝐺^𝑚𝑛𝜕𝑢_𝑖𝜕𝜉_𝑛,(2.5) where 𝐽⁻¹=det(𝜕𝑥_𝑖/𝜕𝜉_𝑗) is the inverse of the Jacobian or the volume of the cell, 𝑈_𝑚=𝐽⁻¹(𝜕𝜉_𝑚/𝜕𝑥_𝑗)𝑢_𝑗 is the volume flux (contravariant velocity multiplied by 𝐽⁻¹) normal to the surface of constant 𝜉_𝑚, and 𝐺^𝑚𝑛=𝐽⁻¹(𝜕𝜉_𝑚/𝜕𝑥_𝑗)(𝜕𝜉_𝑛/𝜕𝑥_𝑗) is the “mesh skewness tensor.” Equations (2.3) and (2.4) are nondimensionalized using the projected area 𝐿_𝑦 of the ellipse as the length scale and the freestream velocity (𝑈_∞) as the velocity scale. Thus, the Reynolds number is given by Re=𝐿_𝑦𝑈_∞/𝜈, where 𝜈=𝜇/𝜌. All of these parameters remain fixed in our computations.

In this study, a body conformal O-type grid is employed to simulate the flow over a body as shown in Figure 1(b). The grid is decomposed into different processors to achieve parallel performance. A non-staggered-grid layout is employed to solve the transformed Navier-Stokes equations. The Cartesian velocity components (𝑢,𝑣) and pressure (𝑝) are defined at the center of the control volume in the computational space, and the volume fluxes (𝑈,𝑉) are defined at the mid points of control volume. All of the spatial derivatives are approximated with second-order accurate central differences except for the convective terms. Using the same central differencing for the convection terms may lead to spurious oscillations in the coarser regions of the grid, thereby leading to erroneous results. In the present formulation, we discretize the convective terms using a variation of QUICK [33]; that is, we calculate the face values of the velocity variables (𝑢_𝑖) from the nodal values using quadratic upwinding interpolation. The upwinding of QUICK is carried out by computing the positive and negative volume fluxes ((𝑈_𝑚+|𝑈_𝑚|)/2) and ((𝑈_𝑚−|𝑈_𝑚|)/2), respectively, and using the generic stencil.

A semi-implicit scheme is employed to advance the solution in time. The diagonal viscous terms are advanced implicitly using the second-order accurate Crank-Nicolson method; whereas all of the other terms are advanced using the second-order accurate Adams-Bashforth method. In the present formulation, we apply a fractional-step method to advance the solution in time. The fractional-step method splits the momentum equation into

(a)an advection-diffusion equation—the momentum equation solved without the pressure term,(b)a pressure Poisson equation—constructed by implicit coupling between the continuity equation and the pressure in the momentum equation, thus satisfying the constraint of mass conservation.

The governing equations are solved using a methodology similar to that employed by Zang et al. [34]. However, the algorithm is extended to parallel computing platforms and the 1D domain decomposition technique (for two-dimensional flows) is employed to distribute the problem among different processors. Details of discretization, parallel implementation, validation, and verification can be found in [35].

The CFD results for the numerical methodology have been validated for three-dimensional [35] and two-dimensional [23] flows by comparing them with published experimental and numerical studies. We have also performed grid and domain dependence studies to verify the numerical results for circular and elliptic cylinders [14].

3. POD-Galerkin Reduced-Order Model

3.1. The Proper Orthogonal Decomposition

Mathematically, we compute Φ for which the following quantity is maximum: ||(𝐮,Φ)||²‖Φ‖²,(3.1) where ⟨⋅⟩ denotes either an ensemble or a time average. Applying variational calculus, one can show that (3.1) is equivalent to a Fredholm integral eigenvalue problem represented as _Ω𝑅_𝑖𝑗𝐱,𝐱^Φ^𝑗𝐱^d𝐱^=𝜆Φ^𝑗(𝐱),(3.2) where 𝑖,𝑗 are the number of velocity components and 𝑅_𝑖𝑗(𝐱,𝐱^) is the two-point spatial autocorrelation tensor.

In the classical POD or direct method, originally introduced by Bakewell and Lumley [36], 𝑅_𝑖𝑗 is a two-point spatial-correlation tensor and the eigenfunctions are the POD modes. In this approach, the averaging operator is performed in time. On the other hand, if the averaging operator is evaluated as a spatial average over the domain of interest, the method is known as the method of snapshots [17]. In this approach, we formulate a temporal-correlation function from the snapshots and transform it into an eigenvalue problem as follows:

𝐶_𝑖𝑗=𝐮^𝑖,𝐮^𝑗,(3.3) where (𝑎,𝑏)=∫_Ω𝑎⋅𝑏dΩ represents the inner product between 𝑎 and 𝑏. The POD modes are then computed by solving the eigenvalue problem

𝐂𝐐=𝐐𝜆,(3.4) where 𝐐 and 𝜆 are the eigenvectors and eigenvalues, respectively. Since 𝐂 is nonnegative Hermitian, 𝐐 is orthogonal by definition. The POD modes are computed as follows:

Φ_𝑖=1√𝜆_𝑖𝒲𝑄_𝑖.(3.5) An important characteristic of these modes is orthogonality; that is, Φ_𝑖⋅Φ_𝑗=𝛿_𝑖𝑗, where 𝛿_𝑖𝑗 is the kronecker delta. The optimality of the POD modes lies in capturing the greatest possible fraction of the total kinetic energy for a projection onto the given set of modes.

The flow data or snapshots of the steady-state velocity field are sampled with a constant time interval (Δ𝑇_𝑠). The velocity field data (𝑢,𝑣) are assembled in a matrix 𝒲_2𝑁×𝑆 as follows:

𝒲=⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣𝑢¹⁽¹⁾𝑢¹⁽²⁾⋯𝑢^1(𝑆)𝑢⋮⋮⋮_𝑁(1)𝑢^𝑁(2)⋯𝑢^𝑁(𝑆)𝑣¹⁽¹⁾𝑣¹⁽²⁾⋯𝑣^1(𝑆)𝑣⋮⋮⋮_𝑁(1)𝑣^𝑁(2)⋯𝑣^𝑁(𝑆)⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦.(3.6) Each column represents one time instant or a snapshot and 𝑆 is the total number of snapshots for N grid points in the domain. The vorticity field can also be used for POD; however, in the case of the velocity field, the eigenvalues of 𝒲 are a direct measure of the kinetic energy in each mode. Deane et al. [24] observed that 20 snapshots are sufficient for the construction of the first eight eigenfunctions at Re=100-200. In general, numerical studies [37] suggest that the first M POD modes, where 𝑀 is even, resolve the first 𝑀/2 temporal harmonics and require 2𝑀 number of snapshots for convergence.

We write the velocity field as the sum of the mean flow (𝐮) and the velocity fluctuations (𝐮^). The mean flow 𝐮=⟨𝐮⟩, where ⟨⋅⟩ is the time average of the assembled data, is subtracted from 𝒲. Then, the fluctuations are expanded in terms of the POD eigenfunctions (Φ_𝑖) as follows:

𝐮(𝐱,𝑡)≈𝐮(𝐱)+^𝑀_𝑖=1𝑞_𝑖(𝑡)Φ_𝑖(𝐱),(3.7) where 𝑀 is the number of POD modes used in the projection. The singular value decomposition (SVD) of this matrix provides the divergence-free velocity POD modes (Φ_𝑖).

3.2. Galerkin Projection

Substituting (3.7) in (2.4) and projecting this equation along the elements of the 𝑀-dimensional POD basis give

𝑞_𝑘̇(𝑡)=𝒜_𝑘+^𝑀_𝑚₌₁ℬ_𝑘𝑚𝑞_𝑚(𝑡)+^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝒞_𝑘𝑚𝑛𝑞_𝑛(𝑡)𝑞_𝑚(𝑡),(3.8) where

𝒜_𝑘=1Re_𝐿𝑦Φ_𝑘ℬ_𝑘𝑚,∇²𝐮−Φ_𝑘,𝐮⋅∇𝐮,=−Φ_𝑘,𝐮⋅∇Φ_𝑚−Φ_𝑘,Φ_𝑚⋅∇𝐮+1Re_𝐿𝑦Φ_𝑘𝒞_𝑘𝑚𝑛,∇²Φ_𝑚,=−Φ_𝑘,Φ_𝑚⋅∇Φ_𝑛.(3.9) In (3.8), 𝒜 is an 𝑀×1 vector resulting from the average flow field, ℬ is an 𝑀×𝑀 matrix corresponding to the linear part, and 𝒞 is a tensor that represents the quadratic nonlinearity of the dynamical system. From linear stability analysis, we can compute the eigenvalues of ℬ and hence ascertain the stability of the trivial solution. Details of the projection for each term are provided in the appendix.

3.3. Numerical Simulations and Reduced-Order Model

We perform numerical simulations of the flows past elliptic cylinders with thickness ratios ranging from 𝜏=0.5 to 𝜏=1.5. The flow is simulated over a 192×256 (𝜉×𝜂) grid with an outer domain size of 30𝐿_𝑦. In the present case, the grid is divided into eight domains/processors, such that each processor has a load of 192×32 grid points, as shown in Figure 1(b). All simulations are performed at Re=200. It is important to note that the Reynolds number is defined with the characteristic length of 𝐿_𝑦=1 and freestream velocity 𝑈_∞=1. The projected area 𝐿_𝑦 of the ellipse, as seen by the flow, is fixed for all values of 𝜏. Thus, the nondimensionalization is fixed for all shape variations.

We develop the reduced-order model for the flow past a circular cylinder (as an example) and use a similar procedure for the model reduction of the flow past elliptic cylinders with different thickness ratios. In the current study, we took 64 snapshots of the flow field over one vortex shedding cycle of the circular cylinder. The flow data or snapshots of the velocity and pressure fields are sampled in the steady-state of the dynamical system with a constant time interval (Δ𝑇_𝑠). The velocity field (𝑢,𝑣) is stored in 𝒲 after subtracting the mean flow. We compute the POD modes from respective snapshot data. In Figure 2, we plot the first six POD modes of u-velocity for the circular cylinder (𝜏=1). Projecting the Navier-Stokes equations onto these modes, we develop a reduced-order model and compute 𝒜,ℬ, and 𝒞 as given in (3.8). We integrate (3.8) to compute the 𝑞_𝑖(𝑡) and plot them in Figure 4. The results we obtained from the reduced-order model match well with the CFD numerical results found in [23].

(a) Mode 1

(b) Mode 2

(c) Mode 3

(d) Mode 4

(e) Mode 5

(f) Mode 6

Figure 2 

Streamwise POD modes for $\tau =1$.

4. Sensitivity Analysis

In this section, we derive the first-order total derivatives of the POD modes with respect to a generic shape parameter 𝛼. We refer to these as the Lagrangian sensitivities of the POD vectors. In our previous study [31], we presented POD sensitivities with respect to a flow parameter, such as the fluid viscosity and the geometry of the structure was fixed. In other words, the sensitivities were computed using the same domain and no grid deformation was required. In this study, the flow domain changes and this leads to a corresponding change in the grid with the elliptic parameter 𝜏. It is important to note that the outer domain of for all the cases is fixed with radius 30𝐿_𝑦.

The sensitivities with respect to the shape parameter 𝛼 (for our case) are computed by a second-order centered finite-difference (FD) approximation

𝐷Φ_𝑖𝐷𝛼𝐱𝛼⁰;𝛼₀||||_FD=Φ^𝑖𝐱𝛼₀+Δ𝛼;𝛼₀+Δ𝛼−Φ_𝑖𝐱𝛼₀−Δ𝛼;𝛼₀−Δ𝛼2Δ𝛼,(4.1) where 𝛼₀ is the parameter value at which the sensitivities are computed and Δ𝛼 is the step size in the finite-difference scheme. The notation 𝐷⋅/𝐷𝛼 represents the total derivative with respect to 𝛼. The parameter increment Δ𝛼 is chosen sufficiently small for the FD computation to be accurate and sufficiently large for the difference between the two nearby POD vectors to be at least one order of magnitude larger than the discretization error. It is also important to note that since the outer domain is the same for all elliptic cylinders and the flow is nearly freestream, these sensitivities become smaller as we move away from the cylinder and approach zero close to the outer domain of the cylinder.

The traditional approach in reduced-order modeling is to build the POD basis for one particular value of the parameter of the system. This will be referred to as the baseline approach and denoted by 𝛼₀, which refers to the baseline state, solution, and POD basis. We aim at producing reduced-order solutions at perturbed states for 𝛼=𝛼₀+Δ𝛼. As in our previous study [31], we define various basis elements used in our analysis.

(i)Baseline POD Basis (BL). This is the classical approach where the POD basis, built from the data at the baseline value 𝛼₀, is used in (3.7) to subsequently produce a reduced-order model at the perturbed value 𝛼. These spatial modes are only available on the baseline geometry but they can easily be mapped onto the perturbed geometry.(ii)Perturbed POD Basis (PR). The reduced-order model is constructed using the POD modes extracted from the solution data obtained by a full-order simulation at the perturbed value 𝛼. This is a costly approach since each new reduced-order simulation requires a full-order data at the perturbed parameter value. Thus, it has little interest in practice but will be used in the remainder of this study as the benchmark low-dimensional solution.

Following the previous studies for flow parameters in [31] and shape parameters in [32], we examine two different ideas for constructing improved reduced-order bases using the POD modes as well as the Lagrangian sensitivity of the POD modes at the baseline value 𝛼₀.

(i)Extrapolated Basis (ET). We treat each POD mode as a function of both space and parameter 𝛼: Φ_𝑖=Φ_𝑖(𝐱;𝛼). A change Δ𝛼 in the parameter from its baseline value 𝛼₀ is reflected in the modes through a first-order expansion in the parameter space: Φ_𝑖(𝐱;𝛼)=Φ_𝑖𝐱𝛼₀;𝛼₀+Δ𝛼𝐷Φ^𝑖𝐷𝛼𝐱𝛼⁰;𝛼₀+𝑂Δ𝛼².(4.2) The effectiveness of this approach clearly depends on whether or not the POD modes exhibit a (nearly) linear dependence with respect to the parameter 𝛼. However, the dimension of the reduced basis is preserved and the reduced approximation of the flow variables still uses (3.7). Once again, the spatial functions Φ and 𝐷Φ_𝑖/𝐷𝛼 are only available on the baseline geometry but they can be mapped on the perturbed geometry in a straight forward manner.(ii)Expanded Basis (EP). The sensitivities of the modes can be shown to span a different subspace than the POD modes (see, e.g., [38]). Thus, it is natural to expect that if the approximated solution is selected in the union of the two subspaces generated by the POD modes and their sensitivities, a broader class of solutions can be represented. The expanded basis consists of the 𝑀 first eigenfunctions with their 𝑀 sensitivities: [𝜙₁;…;𝜙_𝑀;𝐷Φ₁/𝐷𝛼;…;𝐷Φ_𝑀/𝐷𝛼]. The underlying assumption behind this approach is that the subspace spanned by the mode sensitivities is well suited to address the change in the solutions induced by a change in the parameter. However, the dimension of the reduced basis has doubled and the reduced approximation of the flow variables is now expressed as 𝐮(𝐱,𝑡)≈𝐮(𝐱)+𝐷𝐮(𝐱)𝐷𝛼+^𝑀_𝑖=1𝑞_𝑖(𝑡)Φ_𝑖(𝐱)+^2𝑀_𝑖=_𝑀+1𝑞_𝑖(𝑡)𝐷Φ^𝑖𝐷𝛼(𝐱).(4.3)

5. Numerical Results and Discussion

In the current study, we numerically simulate the flow past elliptic cylinders with varying thickness ratios ranging from 𝜏=0.5 to 1.5. In terms of the flow field, the projected area of the cylinder, as seen by the flow, is kept constant (i.e., 𝐿_𝑦=1). In other words, the width of the wake is relatively constant for all the elliptic cylinders. In Figure 3, we plot instantaneous spanwise vorticity contours for various elliptic cylinders. It is interesting to note that as the eccentricity decreases (𝜏<1), the vortex shedding pattern is no longer a von Karman vortex street as observed behind a circular cylinder. We observe different phenomena as the vortex shedding transitions to steady and unsteady secondary shedding. This pattern is the manifestation of bifurcations in the flow as geometry of the cylinder is varied. However, as 𝜏 increases (𝜏>1), the elliptic cylinder is transformed towards an elliptic foil. We still observe von Karman vortex street in the wake; however, the shedding period changes with 𝜏 [14]. The periodic nature of these von Karman vortex flows leads to a better approximation by reduced-order models (as the basis only needs to accurately predict the first period). Based on this observation, we can predict that the sensitivity analysis would tend to show better results for increasing 𝜏 rather than decreasing 𝜏.

(a) τ=0.5

(b) τ=0.6

(c) τ=0.7

(d) τ=0.8

(e) τ=0.9

(f) τ=1.0

Figure 3 

Instantaneous spanwise vorticity contours.

(a) q1  (solid)  and  q2  (dashed)

(b) q3  (solid)  and  q4  (dashed)

Figure 4 

The velocity coefficients $q_i\mathrm{\hspace{0.17em}\hspace{0.17em}}(i=\mathrm{1,2},\mathrm{3,4})$ at $\mathrm{Re}=200$.

From a reduced-order modeling point of view, variation in any geometric parameter would require a new set of flow data or snapshots, computation of the POD modes, and the projection onto these modes to develop a new model. This procedure corresponds to using the perturbed POD basis and is computationally expensive. This undermines the motivation of building a reduced-order model. However, it provides a mean for error calculations and will be used to determine how well the other approaches work. A simple solution could be the baseline approach, where we use the POD basis computed at the baseline, that is, for a circular cylinder (𝜏=1.0), map it onto the elliptic cylinder, and compute the modified flow field using the transformed baseline modes. Clearly, this classical approach where the solution at any parameter value is sought in the subspace spanned by the POD vectors generated at another parameter value performs poorly. This is because the baseline POD vectors are best suited to represent flow solutions at the parameter value for which they have been built.

Using the sensitivity analysis, we can modify the POD basis to include the effect of parameter variation in the flow field. To do so, we consider the thickness ratio of the ellipse (𝛼≡𝜏=𝐿_𝑥/𝐿_𝑦) defining the cylinder shape as a parameter. Initially, we vary 𝜏 by 1% (i.e., 𝜏=0.99) and compute the sensitivity from the finite-difference approach. The accuracy of the method is limited by the step size Δ_𝜏. Using the sensitivity data, we compute the extrapolated and expanded bases and project the snapshot data onto these bases to obtain the time-coefficients (𝑞_𝑖). We plot the two-dimensional projection of the phase portrait for one period in the (𝑞₁,𝑞₂), (𝑞₁,𝑞₃), and (𝑞₁,𝑞₄) planes and compare the results obtained with different basis functions in Figure 5. We observe that the data projected onto all of the bases compare well to the projection onto the perturbed POD basis. Using the coefficient data for each case, we compare the solutions to the full-order solution obtained by DNS at the perturbed state. The relative error in the velocity field 𝑢_mod on the modified geometry Ω_𝛼 for any basis is measured as

Figure 5 

Two-dimensional projections of the phase portraits for different bases.

error=^𝑇0_Ω𝛼𝑢_mod−𝑢_DNS²𝑑Ω_𝛼𝑑𝑡.(5.1) In Figure 6, we plot the relative error obtained for each set of basis vectors. Obviously, the perturbed approach has the least error while the baseline approach corresponds to the maximum error. On the other hand, we observe that the extrapolated and expanded bases perform better than the baseline approach. Comparing the two new bases, the expanded bases seem to work better than the extrapolated bases. The primary reason could be the additional subspace provided by the sensitivity vectors; however, it also doubles the degrees of freedom in this case. Thus, modification of the POD basis from geometric sensitivity broadens the general applicability of the reduced-order models.

Figure 6 

Relative error in reduced-order models at $\tau =0.99$.

We vary 𝜏 in both directions (i.e., 𝜏=0.5-1.5) and compute the corresponding shape sensitivities. We then use the projections and the reduced-order model to quantify the error with respect to perturbed basis (DNS data). We use 6, 12, and 24 POD modes for all cases. We compute the error for the baseline, extrapolated, and expanded basis in terms of velocity data obtained from projections (continuous curves) and reduced-order model (discrete points) as shown in Figure 7. Our first observation is that the error is small in the vicinity of 𝜏≈1. In fact, for a 10% perturbation in the geometric parameter, we gain nearly an order of magnitude improvement in the error over the baseline model. However, as we move away from unity, the error in all models increases and is weakly dependent on the number of modes used in the expansion. It is interesting to note that the error increases gradually as 𝜏 is increased above the value 1 while the error deteriorates quickly for 𝜏 less than 1. This fact can be explained from the flow field observed in the wake of elliptic cylinders. Lower values of 𝜏, while keeping 𝐿_𝑦=1, means more bluffness in the structure and the flow tends to separate due to high inverse pressure gradient. Earlier separation of the shear layer leads to transition from von Karman vortex shedding and the shedding pattern is altered in the wake. This corresponds to a bifurcation in the flow field and mere sensitivity analysis cannot capture the effect. On the other hand, increasing 𝜏, while keeping 𝐿_𝑦=1, smooths the flow and the separation is delayed. Vortex shedding frequency is changed while the shedding pattern remains the same. Thus, the flow structures are similar when 𝜏>1 and can be represented by first-order extrapolations from 𝜏=1.

(a) M=6

(b) M=12

(c) M=24

Figure 7 

Relative error in reduced-order approximations.

In general, the expanded basis tends to perform better than extrapolated basis. Since we assume linear dependency of the POD modes on the geometric parameter, finite-difference approach makes the extrapolated approach close to the baseline methodology. In general, the model based on expanded basis shows relatively better results close to 𝜏.

6. Conclusion

We investigated the possibility of using the POD sensitivity vectors corresponding to a change in shape to improve the accuracy and dynamical system properties of the reduced-order models. As a part of the ongoing research in this area, we modified a circular cylinder to an elliptic cylinder by changing its thickness ratio and computed the sensitivity in the POD modes with respect to this thickness ratio. We defined different POD bases functions, with and without sensitivity, and used them to approximate the velocity field. We then compared the performance of these bases and found that the inclusion of shape sensitivity information in the POD bases performs better than the baseline approach. However, the Sensitivity Analysis for shape parameters is more challenging than for flow parameters since the extrapolation explicitly requires the mapping from one domain to another. The results from the shape sensitivity are encouraging and require further investigation in this field, especially when the parameter changes lead to bifurcations which would require higher order sensitivities (and not merely the first order derivatives).

Appendix

The POD eigenfunctions are used as a basis in a Galerkin projection of the incompressible Navier-Stokes equations. The projection is performed by the inner product of the POD modes with (2.2) as

Φ^𝑘,𝜕𝜕𝐮𝑡+(𝐮⋅∇)𝐮+∇𝑝−R1e∇²𝐮=0,(A.1) where (𝑎,𝑏)=∫_Ω𝑎⋅𝑏dΩ represents the inner product between a and b, and 𝑘=1,2,…,𝑀. We substitute (3.7) into (A.1) and perform the inner product of each term. This leads to a reduced-order model comprising a set of 𝑀 ordinary-differential equations. The inner product reduces the two (or three) momentum equations. To elaborate various terms in the reduced-order model, we expand each term individually.

Time-Derivative Term
From the Galerkin expansion in (3.7), we see that the first term is Φ^𝑘,𝜕𝜕𝐮𝑡=Φ^𝑘,𝜕𝐮(𝐱)𝜕𝑡+𝜕𝐮^(𝐱,𝑡)=Φ𝜕𝑡_𝑘,0+𝜕∑^𝑀𝑚=1𝑞_𝑚(𝑡)Φ_𝑚(𝐱)=Φ𝜕𝑡_𝑘,Φ_𝑚𝑀𝑚₌₁𝑑𝑞_𝑚𝑑𝑡.(A.2) From the definition of the POD modes, the eigenfunctions Φ_𝑖(𝜙^𝑢𝑖,𝜙^𝑣𝑖,𝜙^𝑤𝑖) are orthogonal by construction; that is, Φ^𝑘,^𝑀_𝑚₌₁Φ_𝑚=Φ^𝑘,Φ_𝑘=𝜙^𝑢𝑘,𝜙^𝑢𝑘+𝜙^𝑣𝑘,𝜙^𝑣𝑘+𝜙^𝑣𝑘,𝜙^𝑣𝑘=𝜎_𝑘.(A.3) Substituting (A.3) into (A.2), we obtain Φ^𝑘,𝜕𝜕𝐮𝑡=𝜎^𝑘𝑑𝑞_𝑘𝑑𝑡.(A.4)

Convection Term
We substitute (3.7) into the second term of (A.1) and separate the terms, Φ_𝑘,𝐮⋅∇𝐮==ΦΦ_𝑘,𝐮+𝐮^⋅∇𝐮+𝐮^_𝑘=Φ_𝑘,𝐮⋅∇𝐮+𝐮⋅∇𝐮^+𝐮^⋅∇𝐮+𝐮^⋅∇𝐮^,𝐮⋅∇𝐮+𝐮⋅∇^𝑀_𝑚₌₁𝑞_𝑚Φ_𝑚+^𝑀_𝑚₌₁𝑞_𝑚Φ_𝑚⋅∇𝐮+^𝑀_𝑚₌₁𝑞_𝑚Φ_𝑚⋅∇^𝑀_𝑛=1𝑞_𝑛Φ_𝑛=Φ_𝑘,𝐮⋅∇𝐮+𝐮⋅^𝑀_𝑚₌₁∇Φ_𝑚𝑞_𝑚+^𝑀_𝑚₌₁Φ_𝑚⋅∇𝐮𝑞_𝑚+^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1Φ_𝑚⋅∇Φ_𝑛𝑞_𝑚𝑞_𝑛(A.5) We expand the inner product of each term on the right-hand side of (A.5). The first term is Φ_𝑘,𝐮⋅∇𝐮==𝜙𝜙^𝑢𝑘_𝑢𝑘,𝐮⋅∇𝑢+𝜙^𝑣𝑘,𝐮⋅∇𝑣+𝜙^𝑤𝑘,𝐮⋅∇𝑤𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑦+𝑤𝜕𝑢+𝜙_𝑣𝑘𝜕𝑧𝑢𝜕𝑣𝜕𝑥+𝑣𝜕𝑣𝜕𝑦+𝑤𝜕𝑣+𝜙_𝑤𝑘𝜕𝑧𝑢𝜕𝑤𝜕𝑥+𝑣𝜕𝑤𝜕𝑦+𝑤𝜕𝑤𝜕𝑧(A.6) The second term is Φ^𝑘,𝐮⋅^𝑀_𝑚₌₁∇Φ_𝑚𝑞_𝑚=^𝑀_𝑚=1𝜙^𝑢𝑘𝑢𝜕𝜙^𝑢𝑚𝜕𝑥+𝑣𝜕𝜙^𝑢𝑚𝜕𝑦+𝑤𝜕𝜙^𝑢𝑚𝜕𝑧𝑞^𝑚+^𝑀_𝑚₌₁𝜙^𝑣𝑘𝑢𝜕𝜙^𝑣𝑚𝜕𝑥+𝑣𝜕𝜙^𝑣𝑚𝜕𝑦+𝑤𝜕𝜙^𝑣𝑚𝜕𝑧𝑞^𝑚+^𝑀_𝑚₌₁𝜙^𝑤𝑘𝑢𝜕𝜙^𝑤𝑚𝜕𝑥+𝑣𝜕𝜙^𝑤𝑚𝜕𝑦+𝑤𝜕𝜙^𝑤𝑚𝜕𝑧𝑞^𝑚=𝜙^𝑢𝑘𝑢𝜕𝜙^𝑢𝑚𝜕𝑥+𝑣𝜕𝜙^𝑢𝑚𝜕𝑦+𝑤𝜕𝜙^𝑢𝑚+𝜙_𝑣𝑘𝜕𝑧𝑢𝜕𝜙^𝑣𝑚𝜕𝑥+𝑣𝜕𝜙^𝑣𝑚𝜕𝑦+𝑤𝜕𝜙^𝑣𝑚+𝜙_𝑤𝑘𝜕𝑧𝑢𝜕𝜙^𝑤𝑚𝜕𝑥+𝑣𝜕𝜙^𝑤𝑚𝜕𝑦+𝑤𝜕𝜙^𝑤𝑚𝜕𝑧∗𝐪.(A.7) The third term is Φ^𝑘,^𝑀_𝑚₌₁𝑞_𝑚Φ_𝑚⋅∇𝐮=^𝑀_𝑚₌₁𝜙^𝑢𝑘𝜙^𝑢𝑚𝜕𝑢𝜕𝑥+𝜙^𝑣𝑚𝜕𝑢𝜕𝑦+𝜙^𝑤𝑚𝜕𝑢𝜕𝑧𝑞^𝑚+^𝑀_𝑚₌₁𝜙^𝑣𝑘𝜙^𝑢𝑚𝜕𝑣𝜕𝑥+𝜙^𝑣𝑚𝜕𝑣𝜕𝑦+𝜙^𝑤𝑚𝜕𝑣𝜕𝑧𝑞^𝑚+^𝑀_𝑚₌₁𝜙^𝑤𝑘𝜙^𝑢𝑚𝜕𝑤𝜕𝑥+𝜙^𝑣𝑚𝜕𝑤𝜕𝑦+𝜙^𝑤𝑚𝜕𝑤𝜕𝑧𝑞^𝑚=𝜙^𝑢𝑘𝜙^𝑢𝑚𝜕𝑢𝜕𝑥+𝜙^𝑣𝑚𝜕𝑢𝜕𝑦+𝜙^𝑤𝑚𝜕𝑢+𝜙_𝑣𝑘𝜕𝑧𝜙^𝑢𝑚𝜕𝑣𝜕𝑥+𝜙^𝑣𝑚𝜕𝑣𝜕𝑦+𝜙^𝑤𝑚𝜕𝑣+𝜙_𝑤𝑘𝜕𝑧𝜙^𝑢𝑚𝜕𝑤𝜕𝑥+𝜙^𝑣𝑚𝜕𝑤𝜕𝑦+𝜙^𝑤𝑚𝜕𝑤𝜕𝑧∗𝐪.(A.8) The fourth term is Φ^𝑘,^𝑀_𝑚₌₁𝑞_𝑚Φ_𝑚⋅∇^𝑀_𝑛=1𝑞_𝑛Φ_𝑛=𝜙^𝑢𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛Φ_𝑛⋅∇𝜙^𝑢𝑚+𝜙_𝑣𝑘,^𝑀_𝑀_𝑚=1_𝑛=1𝑞_𝑚𝑞_𝑛Φ_𝑛⋅∇𝜙^𝑣𝑚+𝜙_𝑤𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛Φ_𝑛⋅∇𝜙^𝑤𝑚.(A.9) Equation (A.9) is further split into subterms for convenience as follows: 𝜙^𝑢𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛Φ_𝑛⋅∇𝜙^𝑢𝑚=𝜙^𝑢𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛𝜙^𝑢𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑢𝑚=𝐪_𝜕𝑧∗𝜙^𝑢𝑘𝜙^𝑢𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑧∗𝐪.(A.10) Similarly, the other subterms can be expanded as 𝜙^𝑣𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛Φ_𝑛⋅∇𝜙^𝑣𝑚=𝜙^𝑢𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛𝜙^𝑢𝑛∗𝜕𝜙^𝑣𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑣𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑣𝑚=𝐪_𝜕𝑧∗𝜙^𝑣𝑘𝜙^𝑢𝑛∗𝜕𝜙^𝑣𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑣𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑣𝑚𝜙𝜕𝑧∗𝐪,_𝑤𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛Φ_𝑛⋅∇𝜙^𝑤𝑚=𝜙^𝑢𝑘,^𝑀_𝑚₌₁0𝑥0200𝑑^𝑀_𝑛=1𝑞_𝑚𝑞_𝑛𝜙^𝑢𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑤𝑚=𝐪_𝜕𝑧∗𝜙^𝑤𝑘𝜙^𝑢𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑧∗𝐪.(A.11) Combining (A.10) and (A.11), we obtain Φ^𝑘,^𝑀_𝑚₌₁𝑞_𝑚Φ_𝑚⋅∇^𝑀_𝑛=1𝑞_𝑛Φ_𝑛=𝐪^∗𝜙^𝑢𝑘𝜙^𝑢𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑢𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑢𝑚+𝜙_𝑣𝑘𝜕𝑧𝜙^𝑢𝑛∗𝜕𝜙^𝑣𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑣𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑣𝑚+𝜙_𝑤𝑘𝜕𝑧𝜙^𝑢𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑥+𝜙^𝑣𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑦+𝜙^𝑤𝑛∗𝜕𝜙^𝑤𝑚𝜕𝑧∗𝐪.(A.12) Equation (A.12) represents the quadratic term in the model.

Pressure Term
The pressure term in the model is also projected onto the POD modes as follows: Φ_𝑘,∇𝑝=Φ_𝑘⋅∇𝑝=−∇⋅Φ_𝑘𝑝+_Ω𝑠𝑝𝐧⋅Φ_𝑘.(A.13) We note that using Green's theorem and the divergence-free property, the pressure term drops out from (3.8) for the case of 𝑝=0 on the outerflow boundary Ω_so [19]. The POD eigenfunctions are identically zero on the inflow boundary because the average flow is subtracted from the total flow. However, in case of Neumann boundary conditions on Ω_so, the contribution of the pressure term is not exactly zero for the cylinder wake. The outer domain is intentionally kept at 25𝐷 from the cylinder to minimize the pressure effects. Hence, the pressure is neglected on the outflow boundary so the pressure term vanishes in the reduced-order model [37]. Thus, Φ_𝑘,∇𝑝=0.(A.14)

Diffusion Term
We substitute (3.7) into the diffusion term to obtain Φ^𝑘,1Re_𝐷∇²𝐮=−1Re_𝐷Φ_𝑘,∇²𝐮+𝐮^=−1Re_𝐷Φ^𝑘,∇²𝐮+^𝑀_𝑚₌₁𝑞_𝑚∇²Φ_𝑚.(A.15) We expand each term on the right-hand side of (A.15) into its components. The term containing the mean velocity term is Φ_𝑘,∇²𝐮=𝜙^𝑢𝑘𝜕²𝑢𝜕𝑥²+𝜕²𝑢𝜕𝑦²+𝜕²𝑢+𝜙_𝑣𝑘𝜕𝑧²𝜕²𝑣𝜕𝑥²+𝜕²𝑣𝜕𝑦²+𝜕²𝑣+𝜙_𝑤𝑘𝜕𝑧²𝜕²𝑤𝜕𝑥²+𝜕²𝑤𝜕𝑦²+𝜕²𝑤𝜕𝑧².(A.16) Likewise, the term containing the POD mode is expanded as follows: Φ^𝑘,^𝑀_𝑚₌₁𝑞_𝑚∇²Φ_𝑚=^𝑀_𝑚=1𝑞_𝑚𝜙^𝑢𝑘𝜕²𝜙^𝑢𝑚𝜕𝑥²+𝜕²𝜙^𝑢𝑚𝜕𝑦²+𝜕²𝜙^𝑢𝑚+_𝑀𝜕𝑧²_𝑚₌₁𝑞_𝑚𝜙^𝑣𝑘𝜕²𝜙^𝑣𝑚𝜕𝑥²+𝜕²𝜙^𝑣𝑚𝜕𝑦²+𝜕²𝜙^𝑣𝑚+_𝑀𝜕𝑧²_𝑚₌₁𝑞_𝑚𝜙^𝑤𝑘𝜕²𝜙^𝑤𝑚𝜕𝑥²+𝜕²𝜙^𝑤𝑚𝜕𝑦²+𝜕²𝜙^𝑤𝑚=𝜙_𝑢𝑘𝜕𝑧²𝜕²𝜙^𝑢𝑚𝜕𝑥²+𝜕²𝜙^𝑢𝑚𝜕𝑦²+𝜕²𝜙^𝑢𝑚+𝜙_𝑣𝑘𝜕𝑧²𝜕²𝜙^𝑣𝑚𝜕𝑥²+𝜕²𝜙^𝑣𝑚𝜕𝑦²+𝜕²𝜙^𝑣𝑚+𝜙_𝑤𝑘𝜕𝑧²𝜕²𝜙^𝑤𝑚𝜕𝑥²+𝜕²𝜙^𝑤𝑚𝜕𝑦²+𝜕²𝜙^𝑤𝑚𝜕𝑧²∗𝐪.(A.17)

Acknowledgments

This research was partially supported by the Air Force Office of Scientific Research under Contract FA9550-08-1-0136. Numerical simulations were performed on Virginia Tech Advanced Research Computing—System X. The allocation grant and support provided by the staff are also gratefully acknowledged. Imran Akhtar would like to thank the Government of Pakistan for support during his graduate studies.