1. Introduction
Gaining key insights into fluid flow behaviour enhances our understanding of fluid dynamics and further develops fluid mechanics for the accurate prediction and active control of fluid flows. Data-driven approaches for capturing the essential features of fluid flows have been intensively developed owing to significant improvements in experimental techniques and numerical simulations (Rowley & Dawson Reference Rowley and Dawson2017; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). Proper orthogonal decomposition (POD) is one of the most common approaches for extracting flow structures that characterise the flow field of interest (Lumley Reference Lumley1970). POD analysis seeks the basis, referred to as POD modes, which optimally represents the fluctuation component of the time-series dataset of the flow field obtained through experiments or numerical simulations. Dynamic mode decomposition (DMD) is a common modal-decomposition method that estimates a linear operator approximating the dynamics of interest from time-series snapshot data (Schmid Reference Schmid2010). Each DMD mode has a single frequency in terms of time and growth rate. Both modal-analysis techniques have been extensively employed to understand, predict and control fluid flows. The spectral modal analysis explored in recent years, such as spectral POD (Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018), has successfully extracted coherent structures from turbulent flows.
Modal decomposition, particularly for POD, postulates that the time series of snapshot fluid flow data, represented as a state vector in a high-dimensional Euclidean space, lies in an inherently low-dimensional subspace. Hereafter, the discussion assumes the use of the POD for modal analysis. As POD analysis identifies an optimal subspace, POD modes are ideally suited for constructing a reduced-order model (ROM) of the fluid-flow dynamics. The Galerkin projection-based ROM, which consists of ordinary differential equations (ODEs) of the expansion coefficients with respect to time is a well-known technique for constructing the ROM (Noack, Morzyński & Tadmor Reference Noack, Morzyński and Tadmor2011). Because the ROM approximates the dynamics of the fluid flow in a low-dimensional subspace, the computational cost is significantly lower than that required when solving the full-order model, which is typically the Navier–Stokes equation. Consequently, the ROM has been employed to predict the future state of the flow field and actively control fluid flows (Taylor & Glauser Reference Taylor and Glauser2004).
However, a major challenge is that a typical ROM fails to provide an appropriate solution for flow parameters (e.g. Reynolds number, Mach number and object shapes, such as a cylinder and aerofoils) that is different from that of the parameters under which the POD analysis is performed. This is because POD analysis obtains the optimal set of bases to represent time-series snapshot data under the specific flow conditions analysed. The optimal bases can vary with the flow conditions (Sato, Sakamoto & Ohnishi Reference Sato, Sakamoto and Ohnishi2021). The parametric ROM, whose objective is to perform multiple simulations over a wide range of parameter settings with a low computational cost, has been examined in numerous studies on applications such as optimal design, uncertainty quantification and active control (Benner, Gugercin & Willcox Reference Benner, Gugercin and Willcox2015). A straightforward approach for constructing a parametric ROM is to extract optimal modes from a mixed database that contains datasets with different parameters (Ma & Karniadakis Reference Ma and Karniadakis2002). Galletti et al. (Reference Galletti, Bruneau, Zannetti and Iollo2004) reported that a flow field around a square cylinder was appropriately reproduced at Reynolds numbers not included in the dataset used for modal decomposition. Nakamura, Sato & Ohnishi (Reference Nakamura, Sato and Ohnishi2024) demonstrated a parametric ROM that predicts flow fields around various object geometries by performing POD in computational space instead of physical space. Although these global modes can be used to construct parametric ROM, the global modes are extracted for global optimisation using datasets for all parameters, but are not locally optimal. Generally, the dimensions of the subspace increase with the number of parameters to be included, resulting in an increase in the computational cost of solving ODEs for the ROM.
To construct a parametric ROM with a low computational cost, the dimensions of the subspace should be kept low. This motivated us to investigate and model how the locally optimal subspace, which is obtained by conventional POD analysis, varies with the parameters because a parametric ROM can be constructed by estimating the locally optimal low-dimensional subspace as a function of these flow parameters. This study focuses on a geometric methodology that describes the parameter dependence of POD modes and subspace spanned by the POD modes. We assume that the subspace exhibits small variations for small parameter changes; thus, the subspace is a continuous function of the parameters. To discuss this continuity, it is important to define and compute the distance between subspaces.
To date, the subspaces extracted through the POD analysis have typically been examined independently for each parameter. By contrast, this study aims to analyse the relationships between subspaces corresponding to different parameters by defining the distance between them. The set of all
$r$
-dimensional subspaces of an
$n$
-dimensional vector space admits a manifold structure (
$r$
and
$n$
are natural numbers such that
$r\lt n$
), which is referred to as the Grassmann manifold (Edelman, Arias & Smith Reference Edelman, Arias and Smith1998). Consequently, this is a natural method to consider the distance, defined on the Grassmann manifold, between two subspaces spanned by POD modes with different parameters.
The objective of this study is to develop a geometric methodology for analysing the relationships between the local characteristic structures of fluid flows (i.e. POD modes) and the subspaces spanned by them, which are extracted from the datasets obtained by numerical simulations of the Navier–Stokes equations for parametric variations in the boundary conditions. The distance between two points on the Grassmann manifold enables a quantitative discussion of the similarity between the subspaces corresponding to two different parameters. Moreover, the parameter dependence of the subspace can be represented by modelling the curves or curved surfaces on the Grassmann manifold, enabling the construction of a parametric ROM derived from this representation.
The first part of this study focuses on investigating how the POD modes and the subspace spanned by them vary with flow parameters. We discuss the parameter dependence of the subspace spanned by the POD modes in terms of the geometric characteristics of a curve or curved surface on the Grassmann manifold. The sensitivity of the subspace is defined as the ratio of the subspace displacement along the curve to the parameter change. We show that the subspace sensitivity is closely related to the change in the behaviour of fluid flows with respect to the variation in the parameters. In addition, a method for visualising the distribution of subspaces over a wide range of parameters on the Grassmann manifold is explored.
Subsequently, we discuss the parameter dependence of the POD modes rather than the subspace. Parametric analysis of the set of POD modes for each parameter is conducted on the Stiefel manifold, which is the set of all
$n$
-by-
$r$
orthonormal matrices (Edelman et al. Reference Edelman, Arias and Smith1998); that is, the set of POD modes for each parameter is an element of the Stiefel manifold. Manifolds that have a natural representation of elements in the form of a matrix are referred to as matrix manifolds (Absil, Mahony & Sepulchre Reference Absil, Mahony and Sepulchre2008). Matrix manifolds provide insights into the intrinsic geometric properties of subsets of matrices defined by specific constraints (e.g. orthogonality). We examine the sensitivity of the set of POD modes with respect to parameter changes from the perspective of the Stiefel manifold. Hay, Borggaard & Pelletier (Reference Hay, Borggaard and Pelletier2009) derived the sensitivity of the POD modes by differentiating the POD modes using the equation used in the method of snapshots developed by Sirovich (Reference Sirovich1987), and demonstrated that considering the subspace spanned by the POD modes and POD-mode sensitivity improves the performance of the parametric ROM. In this study, the sensitivity of the POD modes is defined as the tangent vector along a curve on the Stiefel manifold, and the variation in the POD modes with the parameter changes is investigated. We then introduce sensitivity modes based on the analysis of the POD modes on the Stiefel manifold to represent the sensitivity of the flow field with respect to the parameter changes by superimposing the sensitivity modes. The analysis of flow sensitivities enables the estimation of the flow field at nearby parameters of interest, such as the operating or design parameters of the fluid machinery, as well as the estimation of the response of the flow field to input-parameter uncertainty (Pelletier et al. Reference Pelletier, Hay, Etienne and Borggaard2008). The sensitivity modes of the flow field can be interpreted as modes that represent the changes in the flow field in response to small variations in the parameters. We demonstrate that the sensitivity of the flow field can be characterised by analysing the sensitivity of the POD modes.
The last part of this study examines the parametric ROM, which employs the subspace-interpolation technique on the Grassmann manifold. Interpolation of the subspaces spanned by the POD modes at different parameters was employed to construct a parametric ROM. Lieu & Farhat (Reference Lieu and Farhat2007) developed a parametric ROM framework by interpolating two subspaces based on the principal angles between them. Amsallem & Farhat (Reference Amsallem and Farhat2008) generalised this framework as the subspace interpolation method based on the Grassmann manifold and applied the developed framework to an aeroelastic ROM using sets of POD modes obtained using the Euler equation for a full-order model of the fluid flow around an aircraft. Recently, Pawar et al. (Reference Pawar, Ahmed, San and Rasheed2020) developed a parametric ROM framework combining the subspace interpolation on the Grassmann manifold and a long short-term memory neural-network architecture to predict unknown physics that was not included in the training dataset through demonstrations of the Burgers and vorticity-transport equations. Hess, Quaini & Rozza (Reference Hess, Quaini and Rozza2023) employed the manifold-interpolation technique to interpolate the linear operator estimated by the DMD analysis and demonstrated that the Reyleigh–Bénard cavity flows over a wide range of the Grashof numbers can be reconstructed using the developed parametric ROM framework.
This study explores the relationship between the geometric characteristics of the distribution of subspaces on the Grassmann manifold and the reconstruction errors obtained by the parametric ROM using subspace interpolation over a wide range of flow parameters. We demonstrate that the reconstruction error of the flow field obtained by the parametric ROM is closely related to the estimation error of the subspace, which is associated with the sensitivity of the subspace to parameter variations. This finding is critical for constructing a parametric ROM using fewer samples.
The remainder of this paper is organised as follows. Section 2 describes the computational methods for the mode sensitivity analysis and ROM on the matrix manifolds. In § 3, we discuss the sensitivity analysis of the subspace, POD modes and flow field based on the geometric structure of the matrix manifolds using the flow field data around a rotating cylinder as a demonstration. The parametric ROM is performed to discuss the relationship between the reconstruction error and geometric characteristics of a curve or curved surface on the Grassmann manifold in § 4. Finally, § 5 presents the conclusions of this study.
2. Methodology
2.1. Overview of mode sensitivity analysis on matrix manifolds
Matrix manifolds play an important role in various applications, including parametric-model reduction for nonlinear dynamics, which is not limited to fluid dynamics (Son Reference Son2013; Liu & Liu Reference Liu and Liu2022), matrix-completion problems (Boumal & Absil Reference Boumal and Absil2015) and computer vision (Lui Reference Lui2012). Here, we discuss the basic idea of the parametric modal analysis for investigating mode sensitivity and constructing parametric ROM on matrix manifolds before discussing specific manifolds. Figure 1 shows a schematic representation of the sets of POD modes extracted from the flow field dataset over a wide range of flow parameters in terms of a matrix manifold. We consider the flow field around a rotating cylinder, where the flow parameters are the Reynolds number and rotation rate of the cylinder. The POD modes are extracted from the time series of the flow field data obtained from an experiment or numerical simulation for each flow parameter. A single matrix represents a set of POD modes for each flow parameter. The set of these matrices forms a subset of a matrix manifold and is the subject of the analysis in this study. The relationship between the two matrices can be defined by considering a matrix manifold and introducing an appropriate Riemannian metric. This enables the discussion of the relationship between the flow fields at different parameters through the POD modes or the subspace spanned by them. We then consider that the set of POD modes at a specific parameter (or the subspace spanned by them) varies with the Reynolds number or rotation rate. The variation of the POD modes (or subspace) as a function of Reynolds number or rotation rate corresponds to tracing a curve on the matrix manifold. The main objective of this study is to investigate the parameter dependence of fluid flows by analysing and modelling the curve represented by a subset of sampled sets of POD modes or subspaces.
Schematic of the representation of sets of POD modes extracted from the flow field dataset in a wide range of flow parameters in terms of a matrix manifold
$\mathcal{M}$
. The variation of the characteristic structures of the fluid flow around a rotating cylinder with the Reynolds number
${\textit{Re}}$
and rotation rate
$\alpha$
are described as curves on
$\mathcal{M}$
. The relationship between the matrix manifold
$\mathcal{M}$
and a tangent vector space at
$p$
, represented as
$T_{p}\mathcal{M}$
, is also described. Tangent vectors in the tangent-vector space in the Reynolds-number and rotation-rate directions are represented as
$\varDelta _{\textit{Re}}$
and
$\varDelta _{\alpha }$
, respectively.
When analysing or modelling a subset of elements in a matrix manifold (e.g. the Grassmann manifold or Stiefel manifold), considering the tangent vector space at an element on the manifold and performing the analysis in the tangent vector space is more practical in terms of the numerical computations. The Grassmann manifold is not a vector space because it does not admit vector addition and scalar multiplication in a natural way. However, the tangent vector space at each point of the Grassmann manifold has a vector-space structure. This structure allows vector addition and scalar multiplication in the tangent vector space, enabling the consideration of the linear interpolation of the tangent vectors. By considering a mapping, referred to as an exponential map, each tangent vector in the tangent-vector space can be associated with an element on the matrix manifold. Conversely, a mapping that maps an element on the matrix manifold to a vector in the tangent-vector space is referred to as a logarithmic map. A logarithmic map is the inverse of an exponential map.
2.2. Mode sensitivity analysis on Grassmann manifold
First, a brief overview of the POD methods is provided with the purpose of defining some variables. Let
$\boldsymbol{q}_i=\boldsymbol{q}(t_i)\in \mathbb{R}^N$
(
$i=1,2,\ldots ,N_t$
) be snapshot data indicating an instantaneous fluctuating flow field at a given time
$t_i$
, where
$N$
is the dimension of each snapshot data (the number of grid points multiplied by the number of variables) and
$N_t$
is the number of snapshots. A snapshot matrix
$Q$
is defined by arranging the snapshot data as column vectors as follows:
\begin{equation} Q := \begin{bmatrix} \boldsymbol{q}_1 & \boldsymbol{q}_2 & \ldots & \boldsymbol{q}_{N_t} \end{bmatrix} \in \mathbb{R}^{N\times N_t}. \end{equation}
In this study, the following inner product is introduced for two vectors
$\boldsymbol{a},\boldsymbol{b}\in \mathbb{R}^N$
to perform POD analysis:
\begin{equation} \langle \boldsymbol{a},\boldsymbol{b}\rangle := \boldsymbol{b}^T W\!\boldsymbol{a}, \end{equation}
where
$W$
is a diagonal matrix to account for the spatial grid size since this study uses a non-uniform grid for numerical simulation, corresponding to the consideration of the weighted 2-norm induced by the above-mentioned inner product (Schmidt & Colonius Reference Schmidt and Colonius2020).
The POD modes
$\boldsymbol{\phi }_i\in \mathbb{R}^N$
are obtained by solving the following eigenvalue problem:
\begin{equation} QQ^TW\boldsymbol{\phi }_i = \lambda _i\boldsymbol{\phi }_i\quad (i=1,2,\ldots ,r), \end{equation}
where the superscript
$T$
denotes the transpose and
$r$
is the number of POD modes to be employed. In practice, the method of snapshots (Sirovich Reference Sirovich1987) is typically employed when the size of
$N$
is large. The POD-mode matrix
$\varPhi$
is defined by arranging the POD modes as follows:
\begin{equation} \varPhi := \begin{bmatrix} \boldsymbol{\phi }_1 & \boldsymbol{\phi }_2 & \ldots & \boldsymbol{\phi }_r \end{bmatrix} \in \mathbb{R}^{N\times r}. \end{equation}
To account for the inner product introduced in (2.2) within matrix operations, the matrix
$U$
is defined as
\begin{equation} U:=W^{1/2}\varPhi , \end{equation}
where
$W^{1/2}$
denotes the diagonal matrix whose diagonal entries are the square roots of the diagonal entries
$w_1,\ldots ,w_N$
of
$W$
(i.e.
$W^{1/2}:=\textit{diag}(w_1^{1/2},\ldots ,w_N^{1/2})$
). Hence, the orthonormality of the POD modes are expressed as follows:
\begin{equation} U^TU = \varPhi ^TW\varPhi =I_r, \end{equation}
where
$I_r\in \mathbb{R}^{r\times r}$
is the identity matrix.
We then consider a subspace
$\mathcal{S}$
spanned by the POD modes
$\boldsymbol{\phi }_1,\boldsymbol{\phi }_2,\ldots ,\boldsymbol{\phi }_r$
and represent it as
$\mathcal{S}=\mathrm{span}(U)$
. As discussed in previous studies (Hay et al. Reference Hay, Borggaard and Pelletier2009; Sato et al. Reference Sato, Sakamoto and Ohnishi2021), the extracted POD modes depend on the flow parameters of the dataset; therefore, the subspace spanned by the POD modes also depends on the flow parameters. To discuss the dependence of the subspace on the flow parameter from the perspective of geometry, this study considers a set whose elements are subspaces, i.e. the Grassmann manifold. The Grassmann manifold
${\textit{Gr}}(N,r)$
is defined as the set of all
$r$
-dimensional subspaces of
$\mathbb{R}^{N}$
(see, e.g. Absil et al. Reference Absil, Mahony and Sepulchre2008):
\begin{equation} {\textit{Gr}}(N,r) := \{ \mathcal{S}\subset \mathbb{R}^{N}\,|\,\mathcal{S} \mathrm{\,is\,a\,subspace}, \,{\dim}(\mathcal{S}) =r \}. \end{equation}
A subspace spanned by a set of POD modes for a given parameter is represented as a point on the Grassmann manifold. Therefore, a point on the Grassmann manifold represents a linear subspace that can be specified using an orthogonal projector
$P\in \mathbb{R}^{N\times N}$
onto
$\mathcal{S}$
. The orthogonal projector onto
$\mathcal{S}$
is represented as
\begin{equation} P = \textit{UU}^T. \end{equation}
By using a specific orthogonal projector
$P$
, a point on the Grassmann manifold is uniquely specified; however, the matrix size of the orthogonal projector is
$N$
-by-
$N$
. Computations using these orthogonal projectors require a large memory size when the flow field data are used as snapshot data. Therefore, an alternative approach, referred to as the orthonormal-basis perspective (Edelman et al. Reference Edelman, Arias and Smith1998), is considered in this study. A subspace
$\mathcal{S}$
is specified by the matrix
$U$
, whose columns form an orthonormal basis of
$\mathcal{S}$
. This perspective requires an
$N$
-by-
$r$
matrix for computations on the Grassmann manifold, whereas the choice of the matrix is non-unique for determining points on the Grassmann manifold. Two different matrices,
$U_1$
and
$U_2$
, which consist of orthonormal bases as column vectors, span the same subspace when there exists a matrix
$R\in O(r)$
, where
$O(r)$
is the orthogonal group in dimension
$r$
, such that
$R$
satisfies
$U_2=U_1R$
. We consider the following equivalence class:
\begin{equation} [U] = \{ \textit{UR}\,|\,R\in O(r)\}. \end{equation}
Consequently, an alternative representation of the Grassmann manifold is obtained as follows:
\begin{equation} {\textit{Gr}}(N,r) = \{ [U]\,|\,U^TU=I_r \}. \end{equation}
A point on the Grassmann manifold can be identified by specifying
$[U]$
instead of the corresponding orthogonal projector.
The tangent vector space of the Grassmann manifold at point
$[U]\in {\textit{Gr}}(N,r)$
is represented as
\begin{equation} T_{[U]}{\textit{Gr}}(N,r) = \{ \varDelta \in \mathbb{R}^{N\times r}\,|\,U^T\varDelta = 0\}, \end{equation}
where
$\varDelta$
denotes a tangent vector. This study considers the following Riemannian metric:
\begin{equation} \langle \varDelta _1,\varDelta _2 \rangle _{[U]}^{{\textit{Gr}}} = \mathrm{trace}(\varDelta _1^T\varDelta _2), \end{equation}
where
$\varDelta _1,\varDelta _2\in T_{[U]}{\textit{Gr}}(N,r)$
(i.e.
$U^T\varDelta _i=0, (i=1,2)$
).
The subspaces spanned by sets of the POD modes at two different parameters are represented by two different points on the Grassmann manifold. Let us consider the shortest path between two points on the Grassmann manifold, which is referred to as a geodesic. Let
$\gamma (\xi )$
be the geodesic in the Grassmann manifold, parametrised by
$\xi \in \mathbb{R}$
. There exists an open interval containing
$0$
and a unique geodesic that satisfies
$\gamma (0)=[U_p]\in {\textit{Gr}}(N,r)$
and
$ {\rm d}/{{\rm d}\xi }\gamma (0)=\varDelta \in T_{[U_p]}{\textit{Gr}}(N,r)$
(see, e.g. Lee Reference Lee2018). The exponential map on the Grassmann manifold
$\mathrm{Exp}_{[U_p]}^{{\textit{Gr}}}(\varDelta )$
maps a tangent vector
$\varDelta \in T_{[U_p]}{\textit{Gr}}(N,r)$
to a point on the Grassmann manifold
$\gamma (1)\in {\textit{Gr}}(N,r)$
. Hence,
\begin{equation} \mathrm{Exp}^{{\textit{Gr}}}_{[U_p]}: T_{[U_p]}{\textit{Gr}}(N,r) \rightarrow {\textit{Gr}}(N,r),\quad \varDelta \mapsto \gamma (1). \end{equation}
The exponential map on the Grassmann manifold can be written explicitly as follows (Edelman et al. Reference Edelman, Arias and Smith1998):
\begin{equation} \mathrm{Exp}^{{\textit{Gr}}}_{[U_p]}(\varDelta ) = [U_pV\mathrm{cos}(\varSigma )V^T+X\mathrm{sin}(\varSigma )V^T], \end{equation}
where
\begin{equation} \varDelta \stackrel {\textrm {SVD}}{=} X\varSigma V^T, \end{equation}
where
$\stackrel {\textrm {SVD}}{=}$
indicates performing a singular value decomposition on the matrix on the left-hand side. The tangent vector is represented by an
$N$
-by-
$r$
matrix. The tangent vector indicates that it is an element of the tangent vector space
$T_{[U]}{\textit{Gr}}(N,r)$
, which is a vector space.
Let
$[U_p],[U_q]\in {\textit{Gr}}(N,r)$
. We consider a neighbourhood of
$[U_p]$
where there is a unique tangent vector
$\varDelta '\in T_{[U_p]}\in {\textit{Gr}}(N,r)$
such that
$\mathrm{Exp}_{[U_p]}^{{\textit{Gr}}}(\varDelta ')=[U_q]$
. The mapping that determines this tangent vector
$\varDelta '$
is defined in the neighbourhood and is referred to as a logarithmic map:
\begin{equation} \mathrm{Log}_{[U_p]}^{{\textit{Gr}}}:{\textit{Gr}}(N,r)\rightarrow T_{[U_p]}{\textit{Gr}}(N,r),\quad [U_q] \mapsto \varDelta ', \end{equation}
(see Bendokat, Zimmermann & Absil (Reference Bendokat, Zimmermann and Absil2024) for details regarding the neighbourhood where the logarithmic map is defined). The logarithmic map on the Grassmann manifold is computed as follows (Amsallem & Farhat Reference Amsallem and Farhat2008):
\begin{equation} \varDelta ' = X_q\mathrm{tan^{-1}}(\varSigma _q)V_q^T, \end{equation}
where
\begin{equation} (I-U_pU_p^T)U_q(U_p^TU_q)^{-1}=U_q(U_p^TU_q)^{-1}-U_p \stackrel {\textrm {SVD}}{=} X_q\varSigma _q V_q^T. \end{equation}
The exponential and logarithmic maps are inverses of each other. Therefore,
\begin{equation} \mathrm{Exp}_{[U_p]}^{{\textit{Gr}}}\circ \mathrm{Log}_{[U_p]}^{{\textit{Gr}}}([U_q]) = [U_q]. \end{equation}
Equation (2.19) indicates that the output matrix spans the same subspace as the input matrix. However, this does not indicate that the output matrix is identical to the input matrix (Bendokat et al. Reference Bendokat, Zimmermann and Absil2024).
The distance between two points on the Grassmann manifold
$[U_p],\,[U_q]\in {\textit{Gr}}(N,r)$
is calculated as follows (Edelman et al. Reference Edelman, Arias and Smith1998):
\begin{gather} \mathrm{dist}([U_p],[U_q]) = \sqrt {\sum _{i=1}^r\left (\mathrm{cos}^{-1}(\sigma _i)\right )^2}, \end{gather}
where
$\sigma _i$
denotes the
$i$
th singular value of
$U_p^TU_q$
. Note that
$\mathrm{cos}^{-1}(\sigma _i)$
is the principal (or canonical) angle between the two subspaces
$[U_p]$
and
$[U_q]$
. For any two points on the Grassmann manifold
${\textit{Gr}}(N,r)$
, the distance is bounded as (Wong Reference Wong1967)
\begin{equation} \mathrm{dist}([U_p],[U_q]) \leqslant \frac {\pi \sqrt {r}}{2}. \end{equation}
In addition to (2.20), several definitions of the distance between subspaces exist (Edelman et al. Reference Edelman, Arias and Smith1998).
One of the objectives of this study is to investigate the parameter dependence of subspaces spanned by the POD modes of the flow field on the Grassmann manifold by using the above-mentioned definitions. Let
$\xi$
be the flow parameter. The variation in the subspace with respect to
$\xi$
can be described as the motion of a point along the curve
$c$
on the Grassmann manifold,
\begin{equation} c: \mathbb{R}\supset I_\xi \rightarrow {\textit{Gr}}(N,r), \end{equation}
where
$I_\xi$
denotes the interval in which
$\xi$
is defined. In this study, we denote the displacement along the curve
$c$
due to a small change in parameter
$\xi$
by
${\rm d}s/{\rm d}\xi$
, and define it as the sensitivity of the subspace with respect to
$\xi$
. In practice, the sensitivity of the subspace is computed using two points:
$[U(\xi )]$
and
$[U(\xi +\varDelta \xi )]$
:
\begin{equation} \frac {{\rm d}s}{{\rm d}\xi } \approx \frac {\mathrm{dist}\left ( \left [U\left (\xi \right ) \right ],\left [U\left (\xi +\varDelta \xi \right ) \right ] \right )}{\varDelta \xi }. \end{equation}
2.3. Mode sensitivity analysis on Stiefel manifold
An
$r$
-dimensional linear subspace of an
$N$
-dimensional Euclidean space is considered as a point on the Grassmann manifold
${\textit{Gr}}(N,r)$
. The parameter dependence of the subspace spanned by the POD modes is discussed based on the geometric properties of the Grassmann manifold. On the Grassmann manifold, the focus is solely on the subspace without considering the POD modes themselves. However, from the perspective that the contribution of each POD mode to the flow field can be quantified based on the corresponding eigenvalues, the ordering of the modes can be considered meaningful. However, when considering a matrix
$\hat {U}=\begin{bmatrix} \boldsymbol{\phi }_r,\boldsymbol{\phi }_2,\ldots ,\boldsymbol{\phi }_1 \end{bmatrix}$
, which is obtained by interchanging the first POD mode with the
$r$
th POD mode of matrix
$U$
defined in (2.4), matrices
$U$
and
$\hat {U}$
are regarded as the same point on the Grassmann manifold. This means that, when the analysis is based on the Grassmann manifold, it becomes difficult to determine whether the observed subspace variations are induced by changes in physically important modes or by less important ones. To complement the analysis on the Grassmann manifold, this study also conducts analysis on the Stiefel manifold, which enables evaluation of the sensitivity of the POD modes themselves. The Stiefel manifold
${St}(N,r)$
is a set of rectangular, column-orthonormal
$N$
-by-
$r$
matrices (Edelman et al. Reference Edelman, Arias and Smith1998):
\begin{equation} {St}(N,r) := \{ U\in \mathbb{R}^{N\times r}\,|\,U^TU=I_r \}. \end{equation}
The Grassmann manifold is regarded as a quotient manifold of the Stiefel manifold: two column-orthonormal matrices are equivalent if they are related by the multiplication of an orthogonal matrix
$R\in O(r)$
, as indicated in (2.9). Therefore,
\begin{equation} {\textit{Gr}}(N,r) = {St}(N,r)/O(r). \end{equation}
The tangent vector space of
${St}(N,r)$
at
$U$
is represented by
\begin{equation} T_U{St}(N,r) = \{\varDelta \in \mathbb{R}^{N\times r}\,|\,U^T\varDelta = -\varDelta ^TU\}. \end{equation}
That is,
$U^T\varDelta$
is a skew-symmetric matrix. This study considers the following canonical metric (Zimmermann Reference Zimmermann2017):
\begin{equation} \langle \varDelta _1,\varDelta _2 \rangle _{U}^{{St}} = \mathrm{trace}\left ( \varDelta _1^T\left (I_N-\frac {1}{2}\textit{UU}^T \right )\varDelta _2 \right ), \end{equation}
where
$\varDelta _1,\varDelta _2\in T_U{St}(N,r)$
(i.e.
$U^T\varDelta _i=-\varDelta _i^TU,(i=1,2)$
). Note that other choices for the Riemannian metrics on the Stiefel manifold exist (Hüper et al. Reference Hüper, Markina and Leite2021).
As in the case of the Grassmann manifold, a geodesic defines an exponential map on the Stiefel manifold
$\mathrm{Exp}_U^{{St}}(\varDelta )$
that maps a tangent vector
$\varDelta \in T_U{St}(N,r)$
to a point on the Stiefel manifold
$\gamma (1)\in {St}(N,r)$
, where
$\gamma$
is a geodesic:
\begin{equation} \mathrm{Exp}_U^{{St}}:T_U{St}(N,r)\rightarrow {St}(N,r),\,\,\,\varDelta \mapsto \gamma (1). \end{equation}
The exponential map on the Stiefel manifold can be computed using the following matrix operations (readers may refer to Zimmermann (Reference Zimmermann2017) for a more detailed discussion of exponential and logarithmic maps). Let
$U_p$
and
$\varDelta$
be the base point on the Stiefel manifold and tangent vector, respectively. First, we perform QR decomposition of the following matrix:
\begin{equation} \left ( I_N-U_pU_p^T \right )\varDelta = YZ, \end{equation}
where
$Y$
and
$Z$
are the orthogonal and upper triangular matrices, respectively. We then compute the following matrix exponential:
\begin{equation} \begin{bmatrix}M\\N \end{bmatrix} = \mathrm{exp}_m\left ( \begin{bmatrix} U_p^T\varDelta & -Z^T \\ Z & 0 \end{bmatrix} \right ) \begin{bmatrix} I_r \\ 0 \end{bmatrix}, \end{equation}
where
$\mathrm{exp}_m(A):=\sum _{i=1}^\infty A^j/(j!)$
. Finally, by defining
$U_q:=\mathrm{Exp}_{U_p}^{{St}}( \varDelta )\in {St}(N,r)$
, we obtain
\begin{equation} U_q = U_pM+YN. \end{equation}
Let
$U_p,U_q\in {St}(N,r)$
. There is a neighbourhood of
$U_p$
where there is a unique vector
$\varDelta '\in T_{U_p}{St}(N,r)$
such that
$\mathrm{Exp}_{U_p}^{{St}}(\varDelta ')=U_q$
. Consequently, the logarithmic map on the Stiefel manifold is defined in the neighbourhood of
$U_p$
as the inverse map of the exponential map:
\begin{equation} \mathrm{Log}_{U_p}^{{St}}:{St}(N,r)\rightarrow T_{U_p}{St}(N,r),\,\,\,U_q \mapsto \varDelta '. \end{equation}
Computation of the logarithmic map on the Stiefel manifold equipped with a canonical metric requires iterative methods. In this study, we adopt the algorithm proposed by Zimmermann (Reference Zimmermann2017).
We define the sensitivity of the POD modes with respect to
$\xi$
as follows:
\begin{equation} \begin{bmatrix} \dfrac {\partial \boldsymbol{\phi }_1}{\partial \xi } & \dfrac {\partial \boldsymbol{\phi }_2}{\partial \xi } & \ldots & \dfrac {\partial \boldsymbol{\phi }_r}{\partial \xi } \end{bmatrix}:=\frac {\partial U}{\partial \xi }\in T_U{St}(N,r). \end{equation}
In practice, the sensitivity of the POD modes is computed as a tangent vector on the tangent-vector space
$T_U{St}(N,r)$
using two points
$U(\xi )$
and
$U(\xi +\varDelta \xi )$
,
\begin{equation} \frac {\partial U}{\partial \xi } \approx \frac {1}{\varDelta \xi }\mathrm{Log}_{U(\xi )}^{{St}}\left (U\left ( \xi +\varDelta \xi \right ) \right ). \end{equation}
Note that the matrix of the POD modes
$U(\xi +\varDelta \xi )$
is not approximated by
$U(\xi )+(\partial U/\partial \xi )\varDelta \xi$
. The matrix can be approximated using an exponential map as follows:
\begin{equation} U(\xi +\varDelta \xi ) \approx \mathrm{Exp}_{U(\xi )}^{{St}}\left ( \frac {\partial U}{\partial \xi }\varDelta \xi \right ). \end{equation}
By visualising the spatial distribution of POD mode sensitivities, one can identify the flow regions that play a significant role in response to parameter variations, thereby establishing a connection between the subspace sensitivity on the Grassmann manifold and the underlying physical behaviour of fluid flows. The two-step interpretation, from flow field variation, followed by mode sensitivity on the Stiefel manifold and finally subspace sensitivity on the Grassmann manifold, provides insight into understand how physical phenomena influence the (mathematical abstract) geometry on the Grassmann manifold.
Equation (2.34) represents the sensitivity of the POD modes and describes how the POD modes vary with respect to changes in parameter
$\xi$
. It is important to note, however, that the sensitivity of the flow field is not described as a superposition of these POD mode sensitivities. To gain further insight, we consider analysing the sensitivity of the flow field with respect to parameter variations as follows. Here, we explicitly denote the dependence of the snapshot data on the parameter
$\xi$
and represent it as
$Q(\xi )=\begin{bmatrix}\boldsymbol{q}_1(\xi ) & \boldsymbol{q}_2(\xi ) & \ldots & \boldsymbol{q}_{N_t}(\xi ) \end{bmatrix}$
. Our goal is to represent
$\partial {\boldsymbol{q}_i(\xi )}/\partial {\xi }\,(i=1,\ldots ,{N_t})$
as a superposition of several vectors, analogous to representing the flow field as a superposition of POD modes. First, we perform a singular value decomposition of the matrix
$Q(\xi )$
,
\begin{equation} Q(\xi ) \stackrel {\textrm {SVD}}{=} U(\xi )\varSigma (\xi )V^T(\xi ), \end{equation}
where
$U(\xi )\in {St}(N,r)$
is the POD-mode matrix,
$\varSigma (\xi )\in \mathbb{R}^{r\times r}$
is a diagonal matrix with singular values on the diagonal,
$V(\xi )\in {St}(N_t,r)$
is a column-orthonormal
$N_t$
-by-
$r$
matrix of POD expansion coefficients normalised by the corresponding singular values. Note that
$A(\xi ):=\varSigma (\xi )V^T(\xi )$
is a matrix of the expansion coefficients of the POD modes. Consider the partial derivative of (2.36) with respect to
$\xi$
:
\begin{equation} \frac {\partial Q}{\partial \xi } \approx \frac {\partial U}{\partial \xi }\varSigma V^T+U\frac {\partial \varSigma }{\partial \xi }V^T+U\varSigma \frac {\partial V^T}{\partial \xi }. \end{equation}
Equation (2.37) suggests that, in addition to the sensitivity of the POD modes, the sensitivities of matrices
$\varSigma$
and
$V$
are also considered for the sensitivity of the flow field. We can discuss the contributions of the change in the POD modes, singular values and expansion coefficients to the flow field sensitivity with respect to the flow parameter by evaluating the contribution of each term on the right-hand side of (2.37). In this study, we assume the following relationship:
\begin{equation} V(\xi ) \approx V^{\textit{ref}}R(\xi ),\,\,\,R(\xi )\in O(r), \end{equation}
where
$V^{\textit{ref}}$
is defined as the matrix
$V$
when the parameter is a specified reference parameter
$\xi ^{\textit{ref}}$
, i.e.
$V^{\textit{ref}}:=V(\xi ^{\textit{ref}})$
. In other words, we assume that the subspace spanned by matrix
$V$
is identical, regardless of the parameters. The validity of this assumption is discussed in the Appendix A. Substituting (2.38) into (2.37) yields the following equation:
\begin{equation} \frac {\partial Q}{\partial \xi } \approx \tilde {U}A^{\textit{ref}}, \end{equation}
where
\begin{align}& A^{\textit{ref}}:=\varSigma ^{\textit{ref}}\left(V^{\textit{ref}}\right)^T, \end{align}
\begin{align} \tilde {U}:=\frac {\partial U}{\partial \xi }\varSigma R^{T}\left(\varSigma ^{\textit{ref}}\right)^{-1}&+U\frac {\partial \varSigma }{\partial \xi }R^T\left(\varSigma ^{\textit{ref}}\right)^{-1}+U\varSigma \frac {\partial R^T}{\partial \xi }\left(\varSigma ^{\textit{ref}}\right)^{-1}, \end{align}
where
$\varSigma ^{\textit{ref}}:=\varSigma (\xi ^{\textit{ref}})$
. We finally obtain the following representation of the sensitivity of the flow field by defining
$\tilde {U}=:[ \tilde {\boldsymbol{\phi }_1} \quad \tilde {\boldsymbol{\phi }_2} \quad \ldots \quad \tilde {\boldsymbol{\phi }_r} ]\in \mathbb{R}^{N\times r}$
and
$A^{\textit{ref}}=:[\boldsymbol{a}^{\textit{ref}}_1 \quad \boldsymbol{a}^{\textit{ref}}_2 \quad \ldots \quad \boldsymbol{a}^{\textit{ref}}_{N_t} ] \in \mathbb{R}^{r\times N_t}$
(where
$\boldsymbol{a}^{\textit{ref}}_i=[ a^{\textit{ref}}_1(t_i) \quad a^{\textit{ref}}_2(t_i) \quad \ldots \quad a^{\textit{ref}}_r(t_i) ]^T$
):
\begin{equation} \frac {\partial \boldsymbol{q}_i}{\partial \xi } \approx \sum _{j=1}^ra_{\!j}^{\textit{ref}}(t_i)\tilde {\boldsymbol{\phi }_{\!j}}, \,\,\,(i=1,\ldots ,N_t). \end{equation}
Equation (2.42) shows that the sensitivity of the flow field represented by the superposition of
$\tilde {\boldsymbol{\phi }_{\!j}}$
is similar to the representation of the flow field by the superposition of the POD modes. Owing to the assumption of (2.38), the spatial structure, which includes the contributions of the sensitivities of the POD modes, singular values and expansion coefficients, can be represented by the time-independent vector
$\tilde {\boldsymbol{\phi }_{\!j}}$
. Equation (2.42) suggests that the sensitivity of the flow field can be represented by the vectors
$\tilde {\boldsymbol{\phi }_{\!j}}$
with the same expansion coefficients used to represent the flow field by the superposition of the POD modes. Therefore, the spatial structure of
$\tilde {\boldsymbol{\phi }_{\!j}}$
is expected to correspond to a mode that captures important features of the sensitivity of the flow field. In this study, we refer to
$\tilde {\boldsymbol{\phi }_{\!j}}$
as a sensitivity mode. The expansion coefficients
$a_{\!j}^{\textit{ref}}$
are already given; therefore, only the calculation of the sensitivity modes are required to obtain the time evolution of the flow field sensitivity. It is important to note that establishing a relationship between the matrices
$V(\xi )$
and
$V^{\textit{ref}}$
plays a key role. Even when (2.38) is not valid, it is expected that a similar formulation to (2.42) can still be derived if an appropriate relationship can be given between
$V(\xi )$
and
$V^{\textit{ref}}$
, for example,
$V(\xi ) \approx V^{\textit{ref}}B(\xi )$
, where
$B(\xi )$
is not necessarily an orthogonal matrix. To demonstrate that the parameter dependence of the normalised expansion coefficients can be represented by phase shifts in the flow fields considered in this study, we assume a specific relationship in which
$V(\xi )$
is approximated by the product of
$V^{\textit{ref}}$
and an orthogonal matrix
$R(\xi )$
.
2.4. Framework of parametric reduced-order modelling
A ROM based on the Galerkin projection is derived by projecting the Navier–Stokes equations onto the subspace spanned by the POD modes. The set of POD modes is mutually orthonormal, which leads to a well-known system of ODEs for the expansion coefficients. As the elements of the Grassmann manifold are represented as matrices whose column vectors form an orthonormal basis (orthonormal-basis perspective) in this study, a similar ROM based on the Galerkin projection can be derived using the element of the Grassmann manifold. This study considers a parametric-ROM framework by projecting the Navier–Stokes equations onto a locally optimal subspace according to a given flow parameter. The locally optimal subspace is estimated using the subspace-interpolation technique on the Grassmann manifold (Amsallem & Farhat Reference Amsallem and Farhat2008). The strategy considered in this study for constructing a parametric ROM consists of two steps: (1) finding an appropriate subspace to describe the fluid flow to be reconstructed using subspace interpolation on the Grassmann manifold; and (2) constructing a Galerkin projection-based ROM using the estimated subspace.
Here, we estimate the subspace
$ [ U(\xi ) ]\in {\textit{Gr}}(N,r)$
using interpolation with two subspaces
$ [U(\xi _1) ]$
and
$ [U(\xi _2) ]$
, where
$\xi _1 \lt \xi \lt \xi _2$
. Interpolation is performed on the tangent space
$T_{[U(\xi _1)]}{\textit{Gr}}(N,r)$
. First, we obtain the tangent vector
$\varDelta (\xi _2)\in T_{[U(\xi _1)]}{\textit{Gr}}(N,r)$
by mapping the subspace
$ [ U(\xi _2) ]\in {\textit{Gr}}(N,r)$
onto the tangent space using a logarithmic map (2.16). Then, we consider the linear interpolation of the tangent vector:
\begin{equation} \varDelta ([U(\xi )]) = \frac {\xi -\xi _1}{\xi _2-\xi _1}\varDelta ([U(\xi _2)]). \end{equation}
The tangent vector
$\varDelta ([U(\xi _1)])$
corresponds to the origin of the tangent vector space
$T_{[U(\xi _1)]}{\textit{Gr}}(N,r)$
. The interpolation methods commonly employed in vector spaces can be applied to the tangent-vector space. For example, the interpolation of tangent vectors using the Lagrange interpolation has been reported (Pawar et al. Reference Pawar, Ahmed, San and Rasheed2020). This study considers the linear interpolation on the tangent-vector space because we are mainly focused on gaining insight into the relationship between the subspace distribution from the perspective of the geometry of the Grassmann manifold and the errors in the flow field reconstruction using the developed parametric ROM, rather than providing a sophisticated technique for tangent-vector interpolation.
The subspace for parameter
$\xi$
is obtained using an exponential map:
\begin{equation} [U(\xi )] = \mathrm{Exp}_{[U(\xi _1)]}^{{\textit{Gr}}}\left ( \varDelta \left ( \left [U\left ( \xi \right ) \right ] \right ) \right ). \end{equation}
As a result, we have an
$N$
-by-
$r$
matrix
$U'(\xi )=\begin{bmatrix} \boldsymbol{\phi }^{\prime}_1(\xi ) & \boldsymbol{\phi }^{\prime}_2(\xi ) & \ldots & \boldsymbol{\phi }^{\prime}_r(\xi ) \end{bmatrix}\in {St}(N,r)$
, where
$[U'(\xi )]=[U(\xi )]$
is obtained.
The column vectors of
$U'(\xi )$
, which are obtained by subspace interpolation on the Grassmann manifold, do not necessarily coincide with the POD modes (i.e. in general,
$U'(\xi ) \ne U(\xi )$
). The POD modes for a given flow parameter are uniquely determined as the eigenvectors of the covariance matrix constructed from the snapshot data (except for the eigenvector signs). By contrast, subspace interpolation on the Grassmann manifold estimates a matrix
$U'(\xi )$
such that
$U'(\xi )U'^T(\xi )\approx U(\xi )U^T(\xi )$
, where
$U(\xi )$
denotes the POD mode matrix corresponding to a given parameter. In this context, by considering the product of
$U'(\xi )$
and the orthogonal matrix
$R\in O(r)$
, it holds that
$U'(\xi )R(U'(\xi )R)^T=U'(\xi )U'^T(\xi )\approx U(\xi )U^T(\xi )$
. This means that even if
$U'(\xi ) \ne U(\xi )$
, the matrices
$U'(\xi )$
and
$U(\xi )$
span the same subspace. Conversely, even when the matrix obtained via subspace interpolation on the Grassmann manifold spans the same subspace as that of the POD modes corresponding to a given parameter, its column vectors do not necessarily coincide with the POD modes. Nevertheless, the column vectors of the matrix
$U'(\xi )$
are mutually orthonormal and the subspace spanned by these vectors coincides with the subspace spanned by the POD modes. Taking this into account, we refer to the column vectors of
$U'(\xi )$
estimated via subspace interpolation on the Grassmann manifold as pseudo-POD modes.
We estimate the flow fields for the parameter
$\xi$
using the Galerkin projection-based ROM with pseudo-POD modes. Because the pseudo-POD modes, analogous to the POD modes, form orthonormal bases, the following ordinary differential equations can be employed as the Galerkin projection-based ROM:
\begin{equation} \frac {{\rm d}a_i(t;\xi )}{{\rm d}t} = \sum _{j,k=0}^rF_{\textit{ijk}}(\xi )a_{\!j}(t;\xi )a_k(t;\xi )+\sum _{j=0}^rG_{\textit{ij}}(\xi )a_{\!j}(t;\xi )\,\,\,(i=1,\ldots ,r), \end{equation}
where
$a_i$
is the expansion coefficient of the
$i$
th pseudo-POD mode. It is worth noting here that the relative importance of each
$i$
th pseudo-POD mode can be evaluated. In the (space-only) POD analysis, the eigenvalues associated with each POD mode corresponds to the expected value of the squared temporal expansion coefficient. Therefore, by computing the expected value of
$a_i^2 (t;\xi )$
, one can obtain the relative contribution of each pseudo-POD mode. While eigenvalue magnitudes are not directly incorporated when analysing subspace variations of the Grassmann manifold, the physical significance of each pseudo-POD mode can still be evaluated by applying a Galerkin projection within the estimated subspace.
The coefficients
$F_{\textit{ijk}}(\xi )$
and
$G_{\textit{ij}}(\xi )$
can be obtained as follows:
\begin{align}& F_{\textit{ijk}}(\xi ) = -\langle \boldsymbol{\phi }^{\prime}_i(\xi ),(\boldsymbol{\phi }^{\prime}_{\!j}(\xi )\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\phi }^{\prime}_k(\xi ) \rangle \quad (j,k=0,1,\ldots ,r), \end{align}
\begin{align}&\qquad G_{\textit{ij}}(\xi ) = \frac {1}{\textit{Re}}\langle \boldsymbol{\phi }^{\prime}_i(\xi ),{\nabla} ^2\boldsymbol{\phi }^{\prime}_{\!j}(\xi ) \rangle \quad (j=0,1,\ldots ,r). \end{align}
The coefficients
$F_{\textit{ijk}}$
and
$G_{\textit{ij}}$
depend on the flow parameter
$\xi$
because these coefficients are obtained by projecting the Navier–Stokes equations onto the subspace that depends on the parameter. As a result, the parametric ROM estimates appropriate coefficients
$F_{\textit{ijk}}$
and
$G_{\textit{ij}}$
according to the given subspace for the parameter of interest. Note that in this study, we impose Dirichlet boundary conditions on the velocity field: the wall boundary is given by (3.1) and the far-field boundary is set to the free stream conditions. Also, each POD mode is divergence free. Under these conditions, the pressure term drops out from the Galerkin projection-based ROM (Noack et al. Reference Noack, Morzyński and Tadmor2011). Additionally, the expansion coefficients depend on the selection of the basis for the subspace. Therefore, the coefficients
$F_{\textit{ijk}}$
and
$G_{\textit{ij}}$
also depend on basis selection. In particular, the values of
$F_{\textit{ijk}}$
and
$G_{\textit{ij}}$
obtained using the pseudo-POD modes do not necessarily coincide with those derived using the POD modes. This arises from the difference in representing the snapshot data as
$\boldsymbol{q}\approx \sum _{i}^ra_i\boldsymbol{\phi }_i$
or
$\boldsymbol{q}\approx \sum _{i}^ra^{\prime}_i\boldsymbol{\phi }^{\prime}_i$
, which merely reflects the difference in the coordinate system of the subspace. Readers may refer to Noack et al. (Reference Noack, Afansasievm, Morzynski, Tadmor and Thieke2003, Reference Noack, Morzyński and Tadmor2011) for a more detailed discussion on Galerkin projection-based ROM.
The vector
$\boldsymbol{\phi }'(\xi )_0$
indicates the mean field; hence,
$a_0=1$
. The mean field for parameter
$\xi$
is estimated by linear interpolation using
$\boldsymbol{\phi }^{\prime}_0(\xi _1)$
and
$\boldsymbol{\phi }^{\prime}_0(\xi _2)$
:
\begin{equation} \boldsymbol{\phi }^{\prime}_0(\xi ) = \frac {\xi -\xi _1}{\xi _2-\xi _1}\left ( \boldsymbol{\phi }^{\prime}_0(\xi _2)-\boldsymbol{\phi }^{\prime}_0(\xi _1) \right )+\boldsymbol{\phi }^{\prime}_0(\xi _1). \end{equation}
We estimate the initial condition for (2.45) using
$\boldsymbol{q}_1(\xi _1)$
and
$\boldsymbol{q}_1(\xi _2)$
. The expansion coefficients are computed using
$a_i(0)=\langle \boldsymbol{q}_1(\xi _{\!j}),\boldsymbol{\phi }^{\prime}_i \rangle ,\,(j=1,2)$
. The obtained expansion coefficients are used to determine the initial condition by performing a linear interpolation with respect to the parameter
$\xi$
.
3. Example of cylinder flow with varying Reynolds number and rotation rate
In this section and the subsequent one, we consider the mode sensitivity analysis and reduced-order modelling of the flow field around a rotating cylinder to demonstrate the utility of the proposed method. The flow parameters considered in this study are the Reynolds number and rotation rate, which are determined as the boundary conditions of the Navier–Stokes equations. The parameter dependencies of the subspace and POD modes are described on the Grassmann manifold and Stiefel manifold.
3.1. Full-order modelling
Two-dimensional compressible Navier–Stokes equations are considered as the governing equations for the full-order modelling of fluid flow past a cylinder. The governing equations are solved numerically using the finite-difference method. The sixth-order compact-difference scheme (Lele Reference Lele1992) is used to evaluate the spatial derivatives. An eighth-order compact filter is also used to stabilise the numerical dispersion (Visbal & Gaitonde Reference Visbal and Gaitonde2002). The lower-upper symmetric Gauss–Seidel implicit method (Yoon & Jameson Reference Yoon and Jameson1988) is employed for time integration.
An O-type grid is used with 331 and 301 grid points in the radial and circumferential directions, respectively. The computational domain has a radius of
$200D$
, where
$D$
is the cylinder diameter. The minimum grid width in the radial direction is set to
$0.01D$
. The Mach number of the inflow is set to
$0.2$
to ignore compressibility. In this study, the compressible Navier–Stokes equations are used as the governing equations of the full-order model. However, we consider a free stream Mach number that is sufficiently low for the compressibility effects to be negligible. Consequently, we employ a Galerkin projection-based ROM based on the incompressible Navier–Stokes equations by assuming that the obtained flow field data can be regarded as an incompressible flow field.
The following velocity is imposed on the surface of the cylinder:
\begin{equation} \boldsymbol{u}_{{w}} = \boldsymbol{\varOmega }, \end{equation}
where
$\boldsymbol{u}_{{w}}$
and
$\boldsymbol{\varOmega }$
are the velocity on the cylinder surface and angular velocity of the rotating cylinder, respectively. The following rotation rate
$\alpha$
is used as a parameter in this study:
\begin{equation} \alpha = \frac {\varOmega D}{2U_\infty }, \end{equation}
where
$\varOmega$
and
$U_\infty$
are the amplitudes of the angular and inflow velocity, respectively. Along the far-field boundary, free stream conditions are imposed in this study.
A schematic of the numerical-simulation conditions is shown in figure 2(a).
Schematic of the flow field to be described on matrix manifolds in this study: (a) sketch of the flow field around a rotating cylinder; (b,c) instantaneous spatial distributions of
$x$
-velocity component at
${\textit{Re}}=100$
and
$160$
without rotation; (d,e) with the rotation rate of
$\alpha =1.0$
.
Figure 2(b,c) shows the instantaneous spatial distributions of the
$x$
-velocity component without cylinder rotation after the flow field reaches a quasi-steady state at
${\textit{Re}}=100$
and
$160$
, respectively. In both cases, a well-known Kàrmàn vortex street can be observed. The spacing between vortex structures in the Kàrmàn vortex street becomes shorter as the Reynolds number increases. The instantaneous spatial distributions of the
$x$
-velocity component with a rotation rate of
$\alpha =1.0$
for
${\textit{Re}}=100$
and
$160$
are shown in figure 2(d,e). The Kàrmàn vortex street is observed as in the case without rotation and the vortex shedding is deflected upward owing to the cylinder rotation.
Comparison of the Strouhal–Reynolds number between the results obtained by the numerical simulation (circle symbol) and an empirical theory (solid line).
Figure 3 shows the Strouhal number (
$St=fD/U_\infty$
, where
$f$
is the shedding frequency) obtained from the numerical simulation as a function of the Reynolds number, ranging from 60 to 160. The curve obtained from the relationship derived from the following empirical theory (Williamson & Brown Reference Williamson and Brown1998) is also shown:
\begin{equation} St \approx 0.2665-\frac {1.018}{\sqrt {\textit{Re}}}. \end{equation}
The Strouhal–Reynolds number relationship obtained by numerical simulation is in good agreement with the empirical theory.
3.2. Sensitivity analysis on Grassmann manifold
To extract an optimal subspace to represent the dynamics of the fluid flow in a low-dimensional space for each parameter, POD analysis was conducted. For the POD analysis, 1000 snapshots of the spatial distributions of the
$x$
- and
$y$
-velocity components obtained after the flow field reached a quasi-steady state were used. The sampling time corresponds to ten cycles of vortex shedding, that is, the sampling-time width is
$1/(100St)$
. The sampling-time width is adjusted according to the Strouhal number so that the matrix
$V$
satisfies (2.38).
Figure 4 shows the spatial distributions of the extracted first POD modes corresponding to the streamwise (
$x$
-) velocity component in the domain of
$(Re,\alpha )\in [100,160]\times [0.0,1.6]$
. In the following, discussions on the spatial structure of the extracted modes are based on the POD modes corresponding to the
$x$
-velocity component. It should be noted, however, that the POD analysis was performed using snapshot vectors constructed by stacking the spatial distributions of both the
$x$
- and
$y$
-velocity components into a single column vector (Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017). We have confirmed that the similar findings can be observed in the spatial structure of the POD modes associated with the
$y$
-velocity component. For the flow field without cylinder rotation (
$\alpha =0.0$
), an antisymmetric structure with respect to the
$x$
-axis (
$y=0$
) is observed downstream of the cylinder to reconstruct the Kàrmàn vortex street, as shown in figure 4(a–d). The structure of the first POD mode varies smoothly with the Reynolds number. At higher Reynolds numbers, the characteristic structure becomes finer. This corresponds to the gap between the shedding vortices becoming shorter at higher Reynolds numbers. When the cylinder is rotated at the rate of
$\alpha =0.8$
, a similar characteristic structure to that of the non-rotating case is obtained in the first POD mode, whereas the characteristic structure is deflected upward, as shown in figure 4(e–h). Focusing on the Reynolds number dependence of the first POD mode at a rotation rate of
$0.8$
, the characteristic structure becomes finer with an increasing Reynolds number, as in the non-rotating case. Figure 4(i–l) show the first POD modes when the rotation rate is
$1.6$
for different Reynolds numbers. The characteristic structure is further deflected upward. Moreover, the variation in deflection angle is larger when the rotation rate is increased from
$\alpha =0.8$
to
$1.6$
compared with the variation from
$0.0$
to
$0.8$
.
Spatial distributions of the first POD modes corresponding to the
$x$
-velocity component for (a,b,c,d)
$\alpha =0.0$
, (e, f,g,h)
$\alpha =0.8$
and (i,j,k,l)
$\alpha =1.6$
. Panels (a,e,i), (b, f,j), (c,g,k) and (d,h,l) show the POD modes for
${\textit{Re}}=100$
,
$120$
,
$140$
and
$160$
, respectively.
Figure 5 shows the subspace sensitivity with respect to the Reynolds number obtained by (2.23) as a function of the Reynolds number when the dimensions of the subspace are fixed at 12. The rotation rate is fixed at
$\alpha =0.0$
. Numerical simulations using the full-order model and POD analysis are performed in the range of
${\textit{Re}}=60$
to
$160$
with
$\Delta Re = 5$
to calculate the subspace sensitivities. The sensitivity of the subspace decreases smoothly and monotonically as the Reynolds number increases. This suggests that at lower Reynolds number, the variation of the subspace becomes larger per unit Reynolds number, whereas the subspace variation is less significant at higher Reynolds numbers.
Sensitivity of the subspace with respect to the Reynolds number variation as a function of the Reynolds number when the dimension of the subspace is 12. The rotation rate is fixed at
$\alpha =0.0$
.
Relationship between the inverse of the subspace sensitivity and Roshko number
$Ro(=Re\boldsymbol{\cdot }St)$
for different subspace dimensions: (a) inverse of the subspace sensitivity calculated with (2.23); (b) inverse of the subspace sensitivity using normalised line elements by subspace dimension
$\varDelta \hat {s}=\varDelta s/r$
instead of
$\varDelta s$
in (2.23). Fitting curves of the form
$\Delta Re/\varDelta s = a\boldsymbol{\cdot }Ro+b$
are also indicated.
To discuss this in more detail, the inverse of the subspace sensitivity as a function of the Roshko number,
$Ro=Re\boldsymbol{\cdot }St$
for different subspace dimensions is shown in figure 6(a). Fitting curves of the form
$\Delta Re/\Delta s = a\boldsymbol{\cdot }Ro+b$
are also plotted. Curve fitting is performed in the ranges of
$60\leqslant Re \leqslant 70$
and
$60\leqslant Re \leqslant 100$
for subspace dimensions of eight and ten, respectively. For subspace dimensions of 12 and 14, curve fitting is performed in the range of
$60\leqslant Re \leqslant 160$
. At a low Roshko number, the inverse of the subspace sensitivity is approximately proportional to the Roshko number, regardless of the subspace dimension. When the subspace dimension is equal to or greater than 12, the linear relationship shows good agreement with the obtained results in the range of Roshko numbers considered in this study.
Additionally, this result implies that the inverse of the subspace sensitivity becomes zero when
$Ro\approx 2.7$
, which corresponds to
${\textit{Re}}\approx 31$
according to (3.3), regardless of the subspace dimensions. If we consider moving along the curve on the Grassmann manifold in the direction of the decreasing Roshko number, figure 6(a) indicates that the point (subspace) corresponding to
$Ro\lt 2.7$
is inaccessible when moving along a curve from a point corresponding to
$Ro\gt 2.7$
. In other words, the subspace that characterises the flow field for
$Ro\lt 2.7$
cannot be obtained by the continuous deformation of the subspace that characterises the flow field considered in this study. In view of fluid dynamics, no vortex shedding can be excited for
${\textit{Re}}\lt 25$
, which is referred to as the diffusion dominated regime for fluid flow around a cylinder (Ahlborn, Seto & Noack Reference Ahlborn, Seto and Noack2002). Noack & Eckelmann (Reference Noack and Eckelmann1994) reported that no distinct complex conjugate eigenvalue pair defining a characteristic frequency of the flow field is obtained for
${\textit{Re}}\lt 30$
in their stability analysis. These results suggest that the considerable changes in the properties of the eigenvalues and eigenvectors are closely related to the bifurcation. Our results imply that the subspace sensitivity approaches infinity when
${\textit{Re}}\approx 31$
, which corresponds to a significant change in the properties of the POD modes for the flow field around a cylinder. This suggests the potential to detect the existence of bifurcation points by exploring the curve on the Grassmann manifold in the direction towards which the subspace sensitivity approaches infinity. Furthermore, the inverse of the subspace sensitivity collapses to a unified line when the line element is normalised by the subspace dimension (i.e.
$\varDelta \hat {s}:=\varDelta s/r$
), as shown in figure 6(b). These results imply that the geometric features of the curve on the Grassmann manifold are closely related to the features of the fluid flow or the subspace spanned by the POD modes. This motivates the modelling of the parameter dependence of the subspace for the dynamics of fluid flow based on the geometric features of the matrix manifolds.
Figure 7(a) shows the inverse of the subspace sensitivity with respect to the rotation rate as a function of the rotation rate for different subspace dimensions. The Reynolds number is fixed at 100. The subspace sensitivity to the rotation rate is obtained using the numerical simulation of the full-order model, ranging from
$\alpha =0.0$
to
$2.0$
with
$\varDelta \alpha = 0.2$
. The subspace sensitivity is approximately constant for
$\alpha \lt 1$
regardless of the subspace dimension. By contrast, for
$\alpha \gt 1$
, the inverse of the sensitivity rapidly decreases as the rotation rate increases (hence, the subspace sensitivity increases). This trend is consistent with the fact that the variation in the deflection angle of the first POD mode structure is larger when the rotation rate increases from
$\alpha =0.8$
to
$1.6$
compared with the variation from
$\alpha =0.0$
to
$0.8$
, as observed in figure 4.
Relationship between the inverse of subspace sensitivity and rotation rate
$\alpha$
for different subspace dimensions
$r$
: (a) inverse of the subspace sensitivity calculated with (2.23); (b) inverse of the subspace sensitivity using normalised line element. The Reynolds number is fixed at 100.
Figure 7(b) shows the inverse of subspace sensitivity, where the line element is normalised by the subspace dimension as a function of the rotation rate. The inverse of the normalised subspace sensitivities for different dimensions seems to converge to a unified curve at a higher rotation rate and reaches zero at approximately
$\alpha =2.2$
. The flow past a rotating cylinder exhibits a steady state when the rotation rate increases (Hopf bifurcation) and the critical rate is approximately
$\alpha =2.0$
for
${\textit{Re}}=100$
(e.g. Sierra et al. Reference Sierra, Fabre, Citro and Giannetti2020). These results indicate that the inverse of the subspace sensitivity decreases to zero as the parameters approach the vicinity of the Hopf bifurcation point in both the Reynolds number and rotation-rate directions. The points at which the inverse of the sensitivity appears to approach zero are slightly away from the Hopf-bifurcation point (on the side at which a steady solution is obtained). This is because even in a steady flow field, the modes characterising vortex shedding can be defined; however, these modes are stable and decay. This is also observed in the sensitivity of the subspaces with respect to variation in the Reynolds number (figure 6). The point at which the inverse sensitivity appears to approach zero corresponds to the critical Reynolds number for the diffusion-dominated regime. This is different from the bifurcation point (
${\textit{Re}}\approx 49$
), at which the vortex-shedding modes become unstable (Williamson Reference Williamson1996).
In addition to investigating the subspace sensitivity based on the variation in distance with the variation in the parameters, it is useful to visualise the distribution of the subspaces for different Reynolds numbers and rotation rates to understand the parameter dependence of the subspaces. However, visualising the subspace distribution on the Grassmann manifold is difficult when the dimension of the Grassmann manifold is high (
$\mathrm{dim}({\textit{Gr}}(N,r))=(N-r)r$
, Absil et al. Reference Absil, Mahony and Sepulchre2008). This study considers a two-dimensional visualisation of the relative positions of subspaces on the Grassmann manifold for different Reynolds numbers and rotation rates. This is achieved by using a two-dimensional polar-coordinate system defined by the norm of a tangent vector in the tangent-vector space at
$({\textit{Re}}_0,\alpha _0)$
and by the angle, which is determined by the inner products between the two tangent vectors. The norm of the tangent vector
$\varDelta \in T_{({\textit{Re}}_0,\alpha _0)}{\textit{Gr}}(N,r)$
is obtained from (2.12) as follows:
\begin{equation} \|\varDelta \| = \sqrt { \mathrm{trace}\left (\varDelta ^T\varDelta \right ) }. \end{equation}
The angle between tangent vectors
$\varDelta _1,\varDelta _2\in T_{({\textit{Re}}_0,\alpha _0)}{\textit{Gr}}(N,r)$
is defined as follows:
\begin{equation} \theta := \cos ^{-1}\left (\frac {\mathrm{trace}(\varDelta _1^T\varDelta _2)}{\|\varDelta _1\|\|\varDelta _2\|}\right )\quad (0\leqslant \theta \leqslant \pi ). \end{equation}
Subspace distribution in a domain of
$(Re,\alpha )\in [100,160]\times [0.0,1.6]$
based on the norm and angle of the tangent vector in the tangent-vector space at
$(Re,\alpha )=(100,0.0)$
. The curves along with the Reynolds number at
$\alpha =0.0$
(circle symbol) and rotation rate at
${\textit{Re}}=100$
(square symbol) are also shown.
Figure 8 shows the subspace distribution in the domain
$(Re,\alpha )\in [100,160]\times [0.0,1.6]$
(
$\Delta Re=10,\Delta \alpha =0.2$
) based on the norm and angle of the tangent-vector space at
$({\textit{Re}}_0,\alpha _0)=(100,0.0)$
. The norm of the tangent vector is normalised by the norm of the tangent vector for
$(Re,\alpha )=(110,0.0)$
. The angle used for visualisation is defined as the angle between the tangent vector for each parameter and the tangent vector for
$(Re,\alpha )=(110,0.0)$
. The subspace dimension is set to 12 (hereafter, unless otherwise noted, the subspace dimension is 12). Note that the visualisation of the subspace distribution shown in figure 8 depends on the tangent-vector space under consideration because the norm and angles are defined in the tangent-vector space. The corresponding subspace is far from the tangent point when the norm is large. The subspaces are aligned along the geodesic when the angles are constant.
First, we focus on the line describing the subspace variation in the Reynolds number direction when
$\alpha =0.0$
(circle symbol). The angle determined by (3.5) increases with the Reynolds number, suggesting that the variation in the subspace with respect to the Reynolds number is not along a geodesic. Instead, the subspace varies along a curved path with non-zero curvature. In contrast, the variation in the subspace with respect to the rotation rate is relatively along with a geodesic when the Reynolds number is fixed at 100 with a range of
$0.0\leqslant \alpha \leqslant 0.8$
(square symbol in figure 8). The region in which the subspaces are distributed along a straight line corresponds to the range of the rotation rate, where the subspace sensitivity in the
$\alpha$
-direction is constant (figure 7). The angle increases with increasing rotation rate for
$\alpha \geqslant 1.0$
, corresponding to the region in which the inverse of the subspace sensitivity starts to decrease. Moreover, the angle between the tangent vectors corresponding to
$(Re,\alpha )=(100,0.2)$
and
$(110,0.0)$
is almost orthogonal. This suggests that the subspace variation in the
$\alpha$
-direction is completely distinct from that observed in the
${\textit{Re}}$
direction. Consequently, visualisation of the subspace distribution based on the norm and angle of the tangent vectors provides insight into the dependency of the subspace on the parameters. This method is particularly useful for examining whether subspaces are distributed along a geodesic or a curve on the Grassmann manifold in a simple manner.
3.3. Sensitivity analysis on Stiefel manifold
The dependence of the subspace spanned by the POD modes on the flow parameters was previously discussed. In this subsection, we focus on the variation in the POD modes themselves with respect to the flow parameters based on the Stiefel manifold.
The spatial distributions of the sensitivity of the POD modes to the Reynolds number for
$\alpha =0.0$
are shown in figure 9. Figure 9(a–c) show the sensitivity of the first POD modes, which correspond to the first column of the matrix
$\partial U/\partial Re$
defined in (2.33), when
${\textit{Re}}=100$
,
$120$
and
$150$
, respectively. The sensitivities of the POD modes are evaluated based on (2.34), with
$\Delta Re=10$
. The spatial distributions of the first POD-mode sensitivity exhibit an antisymmetric structure with respect to the
$x$
-axis, as observed in the spatial distribution of the first POD modes themselves (figure 4). These characteristic distributions of the POD-mode sensitivity represent the variation in the Reynolds number of the first POD mode, which appears to be dilating in the
$x$
-direction, thereby changing its spatial wavelength while maintaining its asymmetric structure. The sensitivity of the first mode for
${\textit{Re}}=150$
indicates that the variation in the first mode is smaller than those for lower Reynolds numbers. This corresponds to the fact that the subspace sensitivity is lower for higher Reynolds numbers, as observed in the mode sensitivity analysis of the Grassmann manifold.
Spatial distributions of the sensitivity of the POD modes with respect to variation in the Reynolds number at
${\textit{Re}}=100$
,
$120$
and
$150$
: (a,b,c) first POD modes; (d,e, f) third POD modes.
The spatial distributions of the sensitivity of the third POD mode to the Reynolds number are shown in figure 9(d–f). The third-mode sensitivity shows a symmetric pattern with respect to the
$x$
-axis, representing the dilation of the third mode, which is also characterised by a symmetric structure, downstream as the Reynolds number increases. Similar antisymmetric and symmetric structures of the sensitivities of the first and third POD modes are observed in a previous study (Hay et al. Reference Hay, Borggaard and Pelletier2009). As observed in the first-mode sensitivity, the sensitivity of the third mode decreases slightly as the Reynolds number increases.
Figure 10 presents the spatial distributions of the sensitivities of the POD modes to the rotation rate at a fixed Reynolds number of 100. The sensitivities of the POD modes are evaluated using
$\varDelta \alpha =0.2$
in (2.34). The sensitivities of the first and third POD modes at
$\alpha =0.0$
are shown in figure 10(a,d), indicating approximately symmetric and antisymmetric structures with respect to the
$x$
-axis, respectively. These symmetric and antisymmetric structures can be interpreted as representing the shifts along the
$y$
-axis of the antisymmetric and symmetric structures with respect to the
$x$
-axis, respectively. The spatial distributions of the first- and third-mode sensitivities at
$\alpha =0.6$
(figure 10
b,e) show amplitudes of sensitivity similar to those at
$\alpha =0.0$
. In contrast, both the first and third modes at
$\alpha =1.4$
show higher sensitivities than those at
$\alpha =0.0$
and
$0.6$
. This indicates that the structures of these POD modes vary significantly with variations in rotation rate at
$\alpha =1.4$
. These results are consistent with the finding that the subspace sensitivity to the rotation rate is approximately constant for
$0.0\leqslant \alpha \leqslant 1.0$
, while increasing for
$\alpha \gt 1.0$
, as shown in figure 7.
Spatial distributions of the sensitivity of the POD modes with respect to variation in rotation rate at
$\alpha =0.0$
,
$0.6$
and
$1.4$
: (a,b,c) first POD modes; (d,e, f) third POD modes.
Spatial distributions of the sensitivity modes with respect to variation in the Reynolds number at
${\textit{Re}}=100$
,
$120$
and
$150$
: (a,b,c) first sensitivity modes
$\tilde {\boldsymbol{\phi }}^{\textit{Re}}_1$
; (d,e, f) third sensitivity modes
$\tilde {\boldsymbol{\phi }}^{\textit{Re}}_3$
.
Thus far, the sensitivities of the POD modes have been analysed. We now focus on analysing the flow field sensitivities using the POD mode sensitivities defined in (2.42). Figure 11(a–c) show the spatial distribution of the first sensitivity mode with respect to the Reynolds number
$\tilde {\boldsymbol{\phi }}^{\textit{Re}}_1$
at
${\textit{Re}}=100$
,
$120$
and
$150$
when the rotation rate is fixed at
$0.0$
(where superscript indicates the direction of the parameter change). The derivative terms with respect to the Reynolds number in (2.37) are evaluated using the dataset at
${\textit{Re}}_1=100,120,150$
and
${\textit{Re}}_2=110,130,160$
, respectively (i.e.
$\Delta Re=10$
). In this study, the reference Reynolds number is set to
${\textit{Re}}_1$
. As observed in the first POD mode sensitivity distributions, the spatial distributions of the first sensitivity modes exhibit an antisymmetric structure with respect to the
$x$
-axis regardless of the Reynolds number. However, the sensitivity modes exhibit different spatial patterns compared with the sensitivity of the POD modes. This indicates that the contribution of the second and third terms in (2.41) to the sensitivity modes is not negligible.
The contributions of the second and third terms in (2.41) to the first sensitivity mode with respect to the Reynolds number at
${\textit{Re}}=100$
are shown in figure 12(a,b), respectively. The contribution of the first term corresponds to the sensitivity of the POD mode, which is shown in figure 9(a). The second term in (2.41), which is associated with the sensitivity of singular values to variations in the Reynolds number, has smaller amplitudes than the first and third terms, resulting in a slight effect on the sensitivity mode. By contrast, the third term, which is related to the sensitivity of the semi-orthogonal matrix
$V$
, has an amplitude comparable to that of the first term. This suggests that the sensitivity of the POD modes and the sensitivity of matrix
$V$
are necessary to represent the sensitivity modes. The sensitivity of
$V$
is evaluated using the sensitivity of the matrix
$R(Re)$
, as indicated in (2.41). Matrix
$R(Re)$
represents the phase shift of the trajectory of the expansion coefficients (normalised by the singular values) owing to the variation in the Reynolds number (see the Appendix A for details). Therefore, the third term in (2.41) can be interpreted as representing the sensitivity of the flow field caused by the phase shift of the trajectory of the expansion coefficients owing to the Reynolds number variation.
Spatial distributions of the contributions of the (a) second term and (b) third term in (2.41) to the first sensitivity mode with respect to the variation in the Reynolds number at
${\textit{Re}}=100$
.
The first sensitivity mode clearly shows that flow sensitivity decreases as the Reynolds number increases. In particular, figure 11(c) indicates that at
${\textit{Re}}=150$
, the variation in the Reynolds number results in a smaller variation in the flow field compared with the sensitivity at
${\textit{Re}}=100$
. The third sensitivity mode also suggests that the magnitude of the flow field variations owing to the variation in the Reynolds number decreases as the Reynolds number increases, as shown in figure 11(d–f).
The spatial distributions of the sensitivity modes with respect to the rotation rate at
$\alpha =0.0$
,
$0.6$
and
$1.4$
when the Reynolds number is fixed at 100 are shown in figure 13. The derivative terms are evaluated using the datasets at
$\alpha _1=0.0,0.6,1.4$
and
$\alpha _2=0.2,0.8,1.6$
(i.e.
$\varDelta \alpha =0.2$
). The reference rotation rate is set to
$\alpha _1$
. The first and third sensitivity modes at
$\alpha =0.0$
present structures similar to the spatial distributions of the sensitivities of the first and third POD modes (figure 13
a,d). This indicates that the effect of the sensitivity of the POD modes is predominant compared with the second and third terms in (2.41). As the rotation rate increases, the structures of the first and third sensitivity modes evolve into spatial structures that differ from the spatial structures of the corresponding POD mode sensitivities. In particular, the sensitivity modes at
$\alpha =1.4$
suggest that the effects of the second and third terms in (2.41) become significant, as shown in figure 13(c, f).
Spatial distributions of the sensitivity modes with respect to the variation in the rotation rate at
$\alpha =0.0$
,
$0.6$
and
$1.4$
: (a,b,c) first sensitivity modes
$\tilde {\boldsymbol{\phi }}^{\alpha }_1$
; (d,e, f) third sensitivity modes
$\tilde {\boldsymbol{\phi }}^{\alpha }_3$
.
Spatial distributions of the contributions of (a,c) second term and (b,d) third term in (2.41) to the first sensitivity mode with respect to variation in the rotation rate: (a,b) at
$\alpha =0.0$
; (c,d) at
$\alpha =1.4$
.
Figure 14(a,b) shows the contribution of the second and third terms to the first sensitivity mode with respect to the rotation rate at
$\alpha =0.0$
, respectively. In contrast to the symmetric structure of the first term, an antisymmetric structure appears in both the second and third terms. However, the amplitudes of the second and third terms are significantly smaller than those of the first term, resulting in a slight contribution to the sensitivity mode. However, the contribution of the second and third terms are non-negligible when the rotation rate increases. In particular, the contribution of the third term is comparable to that of the first term, as shown in figure 14(c,d). This indicates that the sensitivity of the matrix
$R(\alpha )$
, which represents the phase shift of the expansion coefficient, plays an important role in the sensitivity of the flow field, in addition to the sensitivity of the POD modes when the rotation rate is high.
As discussed previously, parameters associated with high sensitivity in the POD modes also exhibit high subspace sensitivity. This correspondence can be attributed to the fact that the Grassmann and Stiefel manifolds share a closely related geometric structure. Moreover, as shown in figure 8, regions where the subspace distribution deviates from a geodesic path and instead follows a curved trajectory correspond to the parameters that exhibit high POD mode sensitivity. These results suggest that visualising the spatial distribution of the POD mode sensitivity allows us to identify which regions in the flow field contribute to high subspace sensitivity, as well as to the curvature of the trajectory on the Grassmann manifold. Therefore, the framework developed in this study provides a methodology for linking the geometric characteristics of matrix manifolds with the underlying fluid dynamics.
Following the previous discussion of the flow field sensitivity modes, we examine the sensitivity of the flow field by superimposing the sensitivity modes. Figure 15 shows the instantaneous spatial distributions of the
$x$
-velocity component and corresponding spatial distributions of its sensitivity to the Reynolds number at
$(Re,\alpha )=(100,0.0)$
. Note that, in addition to (2.42), the sensitivity distributions (figure 15
b,d, f,h) include the sensitivity of the mean flow field. The spatial distribution of the sensitivity at each time can be interpreted as representing how the flow field at the corresponding time varies with small variation in the Reynolds number. The region of high sensitivity corresponds to the areas in which the Kármán vortex street forms. As the vortex structures advect downstream, the regions of high sensitivity advect accordingly. Therefore, this sensitivity distribution represents the modification of the Kármán vortex street structure with the variation of the Reynolds number. Note that the variation of the Reynolds number affects not only the flow field structures immediately behind the cylinder, but also the structure of the flow field in the downstream region.
Spatial distributions of the flow field sensitivity with respect to variation in the Reynolds number at
$(Re,\alpha )=(100,0.0)$
: (a,c,e,g) instantaneous spatial distributions of the
$x$
-velocity component at
$t/T=0.00$
,
$0.25$
,
$0.50$
and
$0.75$
; (b,d, f,h) distributions of the
$x$
-velocity component sensitivity.
Spatial distributions of the flow field sensitivity with respect to variation in the rotation rate at
$(Re,\alpha )=(100,1.4)$
: (a,c,e,g) instantaneous spatial distributions of the
$x$
-velocity component at
$t/T=0.00$
,
$0.25$
,
$0.50$
and
$0.75$
; (b,d, f,h) distributions of the
$x$
-velocity component sensitivity.
Figure 16 shows the spatial distribution of the
$x$
-velocity component and corresponding spatial distribution of its sensitivity to the rotation rate at
$(Re,\alpha ) = (100,1.4)$
. Unlike the spatial distribution of the sensitivity to the Reynolds number, the sensitivity distribution indicates that the variation in the rotation rate slightly affects the structures of the Kármán vortex structure. In particular, varying the rotation rate has little effect on the flow field in the downstream region (
$x/D\gt 5$
). Instead, it indicates that the flow field has high-sensitivity regions around and immediately behind the cylinder.
As discussed previously, the mode sensitivity analysis on the Stiefel manifold allows the visualisation of the variation in the flow field with the flow parameters. A distinctive feature of this method is its ability to visualise the instantaneous sensitivity of the flow field and represent the temporal evolution of the spatial distribution of the flow field sensitivity based on a superposition of sensitivity modes. The visualisation of the temporal evolution of regions, indicating high sensitivity to variations in the flow parameters, provides meaningful insights for applications in the optimal design and active flow control of fluid machinery.
4. Parametric reduced-order modelling
This section evaluates the performance of the parametric ROM using subspace interpolation on the Grassmann manifold. First, we compare the reconstructed flow fields obtained by the parametric ROM, which employs POD modes estimated by direct interpolation and the parametric ROM, which employs pseudo-POD modes obtained by the method outlined in § 2.4. We then discuss the subspace-estimation and flow-field reconstruction errors over a wide range of flow parameters when using the parametric ROM developed in this study.
4.1. Comparison of parametric reduced-order models based on direct interpolation, global POD and subspace interpolation on Grassmann manifold
We first discuss the performance comparison of the parametric ROMs for estimating the subspace and flow field at a given target Reynolds number
${\textit{Re}}$
using sets of POD modes at two different Reynolds numbers
${\textit{Re}}_1$
and
${\textit{Re}}_2$
(
${\textit{Re}}_1\lt Re\lt {\textit{Re}}_2$
). In this subsection, the target Reynolds number is fixed at 90. The values of
${\textit{Re}}_1$
and
${\textit{Re}}_2$
are determined such that their average is
${\textit{Re}}=90$
. An intuitive direct interpolation of the POD modes to estimate the POD modes at the target Reynolds number is examined for comparison with the results obtained by subspace interpolation on the Grassmann manifold. The direct interpolation of the POD modes is defined as the linear interpolation of the POD matrix
$U({\textit{Re}}_1),U({\textit{Re}}_2)\in {St}(N,r)$
:
\begin{equation} U(Re) = \frac {Re-{\textit{Re}}_1}{{\textit{Re}}_2-{\textit{Re}}_1}\left ( U({\textit{Re}}_2)-U({\textit{Re}}_1)\right )+U({\textit{Re}}_1). \end{equation}
In general, the POD modes estimated by (4.1) do not satisfy orthonormality. We can easily confirm that the directly interpolated POD set lacks orthonormality as follows:
\begin{align} U^T(Re)U(Re)&=\left \{ (1-c)U({\textit{Re}}_1)+cU({\textit{Re}}_2)\right \}^T\left \{ (1-c)U({\textit{Re}}_1)+cU({\textit{Re}}_2) \right \} \nonumber \\ &=(c^2-2c+1)I_r-c(c-1)\left ( U^T\left ({\textit{Re}}_1\right )U\left ({\textit{Re}}_2\right )+U^T\left ({\textit{Re}}_2\right )U\left ({\textit{Re}}_1\right ) \right ) \nonumber \\ &\neq I_r, \end{align}
where
$c=(Re-{\textit{Re}}_1)/({\textit{Re}}_2-{\textit{Re}}_1)$
. This indicates that the set of POD matrices is not closed under the operation of addition, i.e. the summation of the column-orthonormal matrices does not necessarily become a column-orthonormal matrix. Another clear example is the interpolation of POD matrices
$U$
and
$-U$
. Because
$[U]=[-U]$
, the result obtained by an appropriate subspace interpolation should be
$[U]$
. However, direct interpolation yields a matrix in which all the POD modes are zero vectors. In this study, a ROM was constructed using POD modes obtained by interpolating via (4.1), followed by orthonormalisation through QR decomposition (hereafter referred to as ‘direct interpolation’). In contrast to the direct interpolation method, the orthonormality condition with respect to the estimated pseudo-POD modes is rigorously satisfied. The matrix obtained by subspace interpolation on the Grassmann manifold is also an element of the Stiefel manifold, whose elements
$U\in {St}(N,r)$
are defined as
$U^TU=I_r$
; that is, the column vectors are mutually orthonormal. This property is preferable for constructing a Galerkin projection-based ROM, where the orthonormality of the (pseudo-) POD modes plays an important role.
Furthermore, in addition to the direct interpolation method, a performance comparison is conducted between a parametric ROM based on global POD modes and that based on manifold interpolation. The global POD approach is a widely used technique for estimating POD modes in the construction of parametric ROMs (Benner et al. Reference Benner, Gugercin and Willcox2015). Let
$Q_{{\textit{Re}}_1}$
,
$Q_{{\textit{Re}}_2}\in \mathbb{R}^{N\times N_t}$
be the snapshot matrices corresponding to Reynolds numbers
${\textit{Re}}_1$
and
${\textit{Re}}_2$
, respectively. The global snapshot matrix
$Q_{global}$
is defined as follows:
\begin{equation} Q_{\textit{global}} := \begin{bmatrix} Q_{{\textit{Re}}_1} & Q_{{\textit{Re}}_2} \end{bmatrix}\in \mathbb{R}^{N\times 2N_t}. \end{equation}
The global POD modes
$\boldsymbol{\phi }_i\in \mathbb{R}^{N}$
are obtained by solving the following eigenvalue problem:
\begin{equation} Q_{\textit{global}}Q^T_{\textit{global}}W\boldsymbol{\phi }_i = \lambda _i\boldsymbol{\phi }_i. \end{equation}
A subspace spanned by the leading
$r$
global POD modes is then considered. A Galerkin projection-based ROM is constructed on this subspace using (2.45).
Note that for statistically weakly stationary flow, if the two Reynolds numbers
${\textit{Re}}_1$
and
${\textit{Re}}_2$
are identical, the global POD modes are expected to coincide with the POD modes obtained using only
$Q_{{\textit{Re}}_1}$
. This holds even if the snapshot matrices
$Q_{{\textit{Re}}_1}$
and
$Q_{{\textit{Re}}_2}$
are not exactly the same; for example, they may differ by a phase shift in time. To demonstrate this, we first observe that when
${\textit{Re}}_1={\textit{Re}}_2$
, each column of the snapshot matrices
$Q_{{\textit{Re}}_1}$
and
$Q_{{\textit{Re}}_2}$
can be represented as a linear combination of the same leading
$r$
POD modes (i.e. the leading
$r$
POD modes associated with
${\textit{Re}}_1$
). The phase difference between the two snapshot matrices is reflected in the phase shift of the corresponding expansion coefficients. Therefore, the following approximation holds when
${\textit{Re}}_1={\textit{Re}}_2$
:
\begin{equation} Q_{{\textit{Re}}_1}Q^T_{{\textit{Re}}_1} \approx \varPhi _{{\textit{Re}}_1}\varLambda _{{\textit{Re}}_1}\varPhi ^T_{{\textit{Re}}_1} \approx Q_{{\textit{Re}}_2}Q^T_{{\textit{Re}}_2}, \end{equation}
where
$\varPhi _{{\textit{Re}}_1}\in \mathbb{R}^{N\times r}$
and
$\varLambda _{{\textit{Re}}_1}$
are the POD-mode matrix and the diagonal matrix of eigenvalues corresponding to
${\textit{Re}}_1$
, respectively. Thus, we obtain the following:
\begin{equation} Q_{\textit{global}}Q^T_{\textit{global}} = \begin{bmatrix} Q_{{\textit{Re}}_1} & Q_{{\textit{Re}}_2} \end{bmatrix}\begin{bmatrix} Q^T_{{\textit{Re}}_1} \\[6pt] Q^T_{{\textit{Re}}_2} \end{bmatrix} = Q_{{\textit{Re}}_1}Q^T_{{\textit{Re}}_1}+Q_{{\textit{Re}}_2}Q^T_{{\textit{Re}}_2} \approx \varPhi _{{\textit{Re}}_1}\left (2\varLambda _{{\textit{Re}}_1}\right )\varPhi ^T_{{\textit{Re}}_1}. \end{equation}
Therefore, when both
$Q_{{\textit{Re}}_1}$
and
$Q_{{\textit{Re}}_2}$
yield the same POD modes, the global POD modes are expected to coincide with them as well.
The subspace estimated by each method is first examined. When the subspace corresponding to the target Reynolds number is accurately estimated, the POD modes at the target Reynolds number (hereafter, these will be referred to as the reference POD modes) can be represented as a linear combination of the estimated (pseudo-) POD modes (i.e.
$\boldsymbol{\phi }_i = \sum _{j=1}^ra_{\textit{ij}}\boldsymbol{\phi }^{\prime}_{\!j}$
, where
$a_{\textit{ij}}$
is a coefficient). Since the POD modes form an orthonormal basis, considering the squared norm of the
$i$
th reference POD mode yields
\begin{equation} \| \boldsymbol{\phi }_i \|^2 = \langle \boldsymbol{\phi }_i,\boldsymbol{\phi }_i \rangle = \sum _{j=1}^ra_{\textit{ij}}^2=1. \end{equation}
That is, the squared inner product between
$\boldsymbol{\phi }_i$
and
$\boldsymbol{\phi }^{\prime}_{\!j}$
can be interpreted as quantifying the contribution of the
$j$
th estimated (pseudo-) POD mode to the representation of the
$i$
th reference POD mode. Figure 17(a) shows the squared inner product between
$\boldsymbol{\phi }_i$
and
$\boldsymbol{\phi }^{\prime}_{\!j}$
, which is estimated using the direct interpolation method for the case of
$\Delta Re=10$
. There exist reference POD modes
$\boldsymbol{\phi }_i$
for which the inner products with all estimated POD modes
$\boldsymbol{\phi }^{\prime}_{\!j}$
are close to zero (e.g.
$i=2$
or
$4$
). This suggests the presence of reference POD modes that are inadequately approximated by linear combinations of the modes estimated via direct interpolation, as they are nearly orthogonal to the estimated subspace. As shown in figure 17(b), when the global POD method is used, the lower-order reference POD modes can be represented by the linear combination of the estimated global POD modes. In contrast, the inner products between the global POD modes and the higher-order reference POD modes are close to zero, suggesting that the estimated subspace fails to capture these higher-order POD modes. However, as shown in figure 17(c), the pseudo-POD modes estimated via manifold interpolation successfully capture the reference PDO modes across all orders. Moreover, the reference POD mode pairs
$\{\boldsymbol{\phi }_{2k-1},\boldsymbol{\phi }_{2k} \}$
lie in the span of the corresponding pseudo-POD mode pairs
$\{ \boldsymbol{\phi }^{\prime}_{2k-1},\boldsymbol{\phi }^{\prime}_{2k} \}$
for
$k=1,2,\ldots ,6$
, demonstrating that while the modes estimated via manifold interpolation do not directly approximate the reference POD modes themselves, they do span the same subspace. Figure 18(a) shows the sum of
$a_{\textit{ij}}^2$
over
$j$
(i.e.
$\sum _{\!j}a_{\textit{ij}}^2$
) for the case of
$\Delta Re=10$
. If the sum of
$a_{\textit{ij}}^2$
over
$j$
for a given
$i$
equals unity, the
$i$
th reference POD mode lies entirely within the estimated subspace, whereas a value of zero indicates that the
$i$
th reference POD mode lies entirely in the orthogonal complement of the estimated subspace. In the case of the direct interpolation method, the presence of modes with sums close to zero suggests that there exist reference POD modes that are not adequately captured by the span of the estimated modes. In the case of the global POD method, the sums remain close to unity up to the sixth mode, suggesting that these reference modes are effectively captured by the subspace spanned by the global POD modes. In contrast, higher-order modes are not adequately represented within this subspace. The manifold interpolation method, however, yields summations that remain close to unity across all orders, indicating that the estimated subspace agrees well with the subspace spanned by the reference POD modes. Figure 18(b) shows the subspace-estimation error as a function of
$\Delta Re$
for each method. We define the subspace-estimation error as the distance between the true subspace
$[U_{\textit{true}}]$
, which is obtained by the reference POD modes, and estimated subspace
$[U_{\textit{est}}]$
. Regardless of the value of
$\varDelta Re$
, the manifold interpolation method yields the smallest estimation error. Moreover, as
$\Delta Re$
decreases, the estimation error associated with the manifold interpolation method is significantly reduced, implying convergence of the estimated subspace towards
$[U_{\textit{true}}]$
.
Comparison of the squared inner product between the reference POD mode
$\boldsymbol{\phi }_i$
and the estimated mode
$\boldsymbol{\phi }^{\prime}_{\!j}$
: (a) direct interpolation; (b) global POD; (c) manifold interpolation.
Comparison of the subspaces estimated by the three methods (manifold interpolation, direct interpolation and global POD). (a) Evaluation of the representation accuracy of the reference POD modes based on (4.7) for the case of
$\Delta Re=10$
. A value of unity indicates that the
$i$
th reference POD mode lies entirely within the estimated subspace. (b) Subspace estimation error as a function of
$\Delta Re$
.
The errors in the flow field reconstruction obtained using the Galerkin projection-based ROM with estimated POD modes with each method as a function of
$\Delta Re$
are shown in figure 19(a). The error in the flow field reconstruction
$\varepsilon$
is defined as follows:
\begin{equation} \varepsilon = \overline {\int _{\mathcal{D}} \frac {\| \boldsymbol{u}_{\textit{true}}(\boldsymbol{x},t)-\boldsymbol{u}_{\textit{est}}(\boldsymbol{x},t) \|_2}{\|\boldsymbol{u}_{\textit{true}}(\boldsymbol{x},t) \|_2}\,{\rm d}\boldsymbol{x}}, \end{equation}
where
$\boldsymbol{u}_{\textit{true}}$
and
$\boldsymbol{u}_{\textit{est}}$
are the fluctuating component of the true and estimated velocity vectors, respectively, and
$\mathcal{D}$
denotes the simulation domain. The bar indicates time averaging. In this study, the ODEs for the Galerkin-projection ROM (2.45) are solved until the non-dimensional time is 30. While the reconstruction error of the direct interpolation method remains large regardless of
$\Delta Re$
, both the global POD and manifold interpolation methods show a clear reduction in error as
$\Delta Re$
decreases, indicating convergence of the estimated subspace to the reference subspace.
Comparison of reconstruction errors of the flow field estimation at the target Reynolds number (
${\textit{Re}}=90$
): (a) reconstruction errors as a function of
$\Delta Re$
for three methods: manifold interpolation, direct interpolation and global POD (with subspace dimension
$r=12$
); (b) reconstruction errors as a function of
$\Delta Re$
for different subspace dimensions using the manifold interpolation method; (c) reconstruction errors as a function of
$\Delta Re$
for different subspace dimensions using the global POD method.
Although the global POD method and the manifold interpolation method exhibit different trends in subspace estimation error (figure 18
b), the trends in reconstruction error by ROM are similar. Figure 19(b) shows the reconstruction error of the flow field as a function of
$\varDelta Re$
for ROMs based on the manifold interpolation method, constructed using different numbers of POD modes. When
$\Delta Re$
is greater than 30, the reconstruction error shows little sensitivity to the number of modes; however, for
$\Delta Re$
less than 30, the error associated with
$r=4$
becomes larger than those observed for
$r=6$
and
$r=12$
. This result indicates the following: when
$\Delta Re$
is large, the subspace estimation error becomes large, which limits the effectiveness of increasing the number of modes in reducing the reconstruction error. In contrast, when
$\Delta Re$
becomes small, the estimated subspace better approximates the one spanned by the reference POD modes. As a result, increasing the number of modes used in the ROM leads to a more accurate reconstruction. The similarity in reconstruction errors between the cases of
$r=6$
and
$r=12$
can be attributed to the fact that the cumulative contribution ratio at the target Reynolds number exceeds 99.9 % when
$r=6$
, indicating the contribution of POD modes higher than the sixth is negligible in this case.
It should be noted that, in subspace interpolation on the Grassmann manifold, higher-order modes may sometimes be difficult to converge and can contain numerical noise, which may affect the results to a similar extent as noise contaminating lower-order modes, despite their relatively small eigenvalues. While all POD modes are treated equally on the Grassmann manifold, figure 19(b) demonstrates that the reconstruction error tends to decrease with increasing subspace dimension, particularly for small values of
$\Delta Re$
, supporting the validity of the framework employed in this study.
Figure 19(c) shows the reconstruction error obtained by ROM using the global POD modes with different numbers of modes. A similar trend to that observed with manifold interpolation can be seen overall. However, in contrast to the manifold interpolation, the reconstruction errors exhibit similar trends for
$r=4$
and
$r=6$
, whereas for
$r=12$
, a smaller error is achieved when
$\Delta Re\lt 30$
. These results highlight a key difference between the manifold interpolation and global POD methods. In the case of manifold interpolation, figure 19(b) implies that six POD modes are sufficient to accurately reconstruct the fluid flow. However, for the global POD method, figure 19(c) suggests that six modes are not sufficient and higher-order global POD modes are required to achieve accurate flow field reconstruction. Indeed, figure 17(b) shows that the fifth and sixth reference POD modes are represented by linear combinations of global POD modes of order higher than six, specifically the seventh and eighth modes.
Consequently, both the global POD and manifold interpolation methods can accurately predict flow fields at unseen parameter values when a sufficient number of POD modes is used. However, to achieve low reconstruction errors, the global POD approach requires more modes than manifold interpolation. This contrast reflects the difference between the two methods: the global POD method represents the flow field at a target parameter within a single linear subspace that is optimal for capturing the flow field across all parameters in the training dataset, whereas the manifold interpolation estimates parameter-dependent, locally optimal subspace to achieve more efficient reconstruction of the flow field at the target parameter. This difference becomes significant in the construction of parametric ROM to predict flow fields over a wide range of parameters as will be discussed in the following subsection.
Exploration of the spatial distribution of the flow field reconstructed by the ROM provides insights to characterising the property of the parametric ROM developed in this study. Figure 20(a) shows the time-averaged spatial distribution of the square of the velocity fluctuation of the
$x$
-component obtained using the full-order model. Figure 20(b–d) present the spatial distributions estimated by ROMs based on the direct interpolation method, the global POD method and the manifold interpolation method, respectively, under the condition of
$\Delta Re=30$
and
$r=6$
. While the spatial distribution estimated by the direct interpolation method deviates significantly from that of the full-order model, both the global POD and manifold interpolation methods capture spatial structures that are in good qualitative agreement with the result of the full-order model. Figure 20(e–g) show the spatial distributions predicted by each method when
$\Delta Re=10$
and
$r=6$
. Even with a small
$\Delta Re$
, the direct interpolation method still fails to reproduce the spatial structure observed in the full-order model. In contrast, as
$\Delta Re$
decreases, both global POD and manifold interpolation methods obtain more accurate reconstructions of the time-averaged spatial distribution of the squared
$x$
-velocity fluctuation.
Comparison of the time-averaged spatial distribution of the variance of the
$x$
-velocity fluctuation: (a) full-order model; (b,e) direct interpolation; (c, f) global POD; (d,g) manifold interpolation. Panels (b,c,d) and (e, f,g) show the results for
$\Delta Re=30$
and
$\Delta Re=10$
, respectively.
Comparison of the profiles of (a) variance of the
$x$
-velocity fluctuation, (b) variance of the
$y$
-velocity fluctuation and (c) covariance between the
$x$
- and
$y$
-velocity fluctuations at
$x/D=2.4$
obtained from manifold interpolation, direct interpolation, global POD and full-order model for
$\Delta Re=10$
and
$r=6$
.
To enable a more detailed comparison, figure 21(a–c) respectively present the profiles of the variance of the
$x$
-velocity fluctuation, the variance of the
$y$
-velocity fluctuation, and the covariance between the
$x$
- and
$y$
-velocity fluctuations at
$x/D=2.4$
, obtained from each method (
$\Delta Re=10$
and
$r=6$
). For all the profiles, the results obtained using the direct interpolation method lack symmetry and deviate significantly from the full-order model. In contrast, both the global POD and manifold interpolation methods yield profiles that are in good agreement with those obtained by the full-order model, with the manifold interpolation method tending to yield slightly improved results. The above-mentioned results are based on a case where the flow field at an unseen parameter is estimated using only two different training parameters, under which both the global POD and manifold interpolation methods yield similar performance; specifically, reducing
$\Delta Re$
leads to a decrease in reconstruction errors.
4.2. Error evaluations of subspace estimation and flow field reconstruction: one-parameter variation case
Here, we evaluate the errors of the subspace estimation and flow field reconstruction over a wide range of Reynolds numbers using a parametric ROM based on subspace interpolation on the Grassmann manifold. Parametric ROMs constructed using subspaces sampled at different
$\Delta Re$
values are examined to discuss convergence of the estimation error with respect to the subspace-sampling interval
$\Delta Re$
. Table 1 summarises the conditions for
$\Delta Re$
considered in this study, the Reynolds numbers used to interpolate the subspaces (denoted as the training dataset), and the Reynolds numbers used to estimate the subspaces and reconstruct the flow field (denoted as the test dataset). The Reynolds numbers for subspace and flow field estimation are determined so that their average coincides with the Reynolds number of the training dataset, i.e. the Reynolds number for the estimation is
${\textit{Re}}=({\textit{Re}}_1+{\textit{Re}}_2)/2$
when using the subspaces at
${\textit{Re}}_1$
and
${\textit{Re}}_2(={\textit{Re}}_1+\Delta Re)$
.
Sampling points
${\textit{Re}}$
used as the training dataset for POD modes and Reynolds numbers for evaluation (denoted as the test dataset) for each
$\Delta Re$
.
Figure 22(a) shows the subspace-estimation errors as functions of the Reynolds number for different
$\Delta Re$
values. The subspace-estimation error decreases with decreasing
$\varDelta Re$
for all Reynolds numbers considered in this study. In addition, the error decreases with increasing the Reynolds number, regardless of
$\varDelta Re$
, because the subspace sensitivity decreases with increasing the Reynolds number (see figure 5). The error in the flow field reconstructed using the parametric ROM as a function of the Reynolds number for different
$\Delta Re$
values is shown in figure 22(b). The reconstruction error decreases with decreasing
$\Delta Re$
and increasing Reynolds number for a fixed
$\Delta Re$
. This highlights that the reconstruction error is closely related to the subspace-estimation error.
Error evaluation of the parametric ROM across a wide range of Reynolds numbers using different
$\Delta Re$
; (a) subspace error; (b) flow field reconstruction error. The subspace dimension is fixed to 12.
Comparison of reconstruction errors of the parametric ROM across a wide range of Reynolds numbers: (a) reconstruction errors as a function of Reynolds number using the manifold interpolation, global POD and local POD methods, with a fixed subspace dimension of 12; (b) reconstruction errors obtained using the manifold interpolation method as a function of Reynolds number for different subspace dimensions.
We also compare the reconstruction errors of the flow field estimated using the global POD modes and the manifold interpolation method. Figure 23(a) shows the reconstruction errors of the flow field obtained by each method as a function of Reynolds number for the case of
$\varDelta Re=10$
. The global POD modes were computed from the covariance matrix constructed from the combined snapshot matrices at
${\textit{Re}}=60,70,\ldots ,140$
. For comparison, we also present the reconstruction error obtained using a ROM constructed with POD modes computed from the snapshot data at
${\textit{Re}}=95$
(referred to as local POD). While the reconstruction error for the local POD remains sufficiently small at the Reynolds number used for the POD analysis (i.e.
${\textit{Re}}=95$
), it increases significantly as the Reynolds number deviates from this value. In contrast, the global POD method demonstrates lower reconstruction errors compared with those of the local POD method over a broad range of Reynolds numbers. Moreover, the manifold interpolation method yields lower reconstruction errors compared with the global POD method, suggesting that manifold interpolation provides a more robust framework for problems involving flow field prediction in a wide parameter range. Figure 23(b) presents the reconstruction errors of the flow field estimated using the subspace interpolation on the Grassmann manifold as a function of the Reynolds number for different subspace dimensions. For all Reynolds numbers considered, the reconstruction error decreases as the subspace dimension increases. This indicates that, within the range of Reynolds numbers and subspace dimensions addressed in this study, incorporating higher-order modes leads to reduced reconstruction errors, resulting in more reliable parametric ROMs.
4.3. Parametric reduced-order modelling of flow field around a rotating cylinder, and error evaluations of subspace estimation and flow field reconstruction
We then evaluate the performance of the parametric ROM, which estimates the locally optimal subspace and reconstructs the flow field at a given Reynolds number and cylinder-rotation rate. Figure 24 shows the flow parameters (a combination of the Reynolds number
${\textit{Re}}$
and rotation rate
$\alpha$
) employed to interpolate the subspaces (training dataset) and parameters used to estimate the subspaces and reconstruct the flow field (test dataset). Two types of training datasets, coarse- and fine-sampling data sets, are used in this study. The subspace estimation and flow field reconstruction errors are compared between the two datasets. For the coarse dataset, the subspaces are sampled at intervals of
$\Delta Re=20$
and
$\Delta \alpha =0.4$
. For the fine dataset, the subspaces are sampled at intervals of
$\Delta Re=10$
and
$\Delta \alpha =0.2$
.
Sampling points
$(Re,\alpha )$
used as the coarse (circle symbol) and fine (triangle symbol) training datasets, and test-data set (diamond symbol) for evaluation of the errors of parametric ROM in two-dimensional parameter space.
The parametric ROM for two-dimensional parameter space is performed in two steps, as in the one-dimensional case: the subspace-estimation step and Galerkin projection-based ROM step. First, four subspaces corresponding to the conditions closest to the target parameters are selected from the training dataset. These subspaces are mapped onto the tangent-vector space using a logarithmic map (2.16). The subspace for the target-flow condition is estimated using bilinear interpolation in the tangent vector space. The initial conditions and mean field required for the Galerkin projection-based ROM are also estimated using bilinear interpolation. In the second step, using the estimated subspace, initial condition and mean field, the ODEs (2.45) for the Galerkin projection-based ROM are used to reconstruct the flow field. The second step follows the same procedure as that in the case of a one-dimensional parameter space.
Error distributions of the parametric ROM in two-dimensional parameter space: (a,b) subspace-estimation errors when using coarse and fine training datasets, respectively; (c,d) flow field reconstruction errors.
Figure 25(a,b) shows the error distributions of the subspace estimation in
$\alpha$
-
${\textit{Re}}$
space for the coarse and fine datasets, respectively. The subspace-estimation error decreases as the Reynolds number increases for both the coarse and fine datasets. This is related to the dependence of the subspace sensitivity on the Reynolds number, as discussed in § 3.2. As shown in figure 5, the subspaces are not uniformly distributed on the Grassmann manifold when sampled at a constant interval of
$\Delta Re$
. The interval between two different subspaces on the Grassmann manifold
$\Delta s$
increases at lower Reynolds numbers. Therefore, the subspace-estimation error increases as the Reynolds number decreases when the subspaces are sampled at a constant
$\Delta Re$
. In addition, the error increases as the rotation rate increases, particularly in the region where
$\alpha \gt 1.1$
, which corresponds to the subspace sensitivity increasing with respect to the rotation rate for
$\alpha \gt 1.1$
(figure 7). Furthermore, the estimation error decreases when the subspace is estimated using the fine dataset compared with the coarse dataset. This indicates that reducing the sampling interval of the subspaces used for interpolation also leads to a decrease in the estimation error.
The error distributions of the flow field reconstruction by the parametric ROM using the coarse and fine datasets are shown in figure 25(c,d), respectively. The flow field reconstruction error follows a trend similar to that of the subspace-estimation error: the reconstruction error decreases as the Reynolds number increases and the rotation rate decreases. In addition, the use of a finer dataset reduces the reconstruction error because the use of a finer dataset reduces the subspace-estimation error. This result suggests that, when constructing a parametric ROM, the subspace-estimation error plays a dominant role in determining the flow field reconstruction error. This highlights that minimising the subspace-estimation error is important for reducing the error in the flow field estimated by the parametric ROM. A parametric ROM that reproduces flow fields with a small error over a wide range of parameters can be efficiently constructed by sampling the subspaces finely in regions with high subspace sensitivity and coarsely in regions with low sensitivity, because the subspace-estimation error is related to the subspace sensitivity.
5. Conclusions
In this study, the sensitivity analysis of POD modes and the subspace spanned by them with respect to changes in flow parameters were analysed from the perspective of matrix manifolds. This approach provides insight into how the locally optimal structures and fluid flow evolve with parameter variations, leading to the construction of a reliable parametric ROM. The sets of POD modes and subspaces spanned by them over a wide range of flow parameters were represented as subsets on the Stiefel manifold and Grassmann manifold, respectively. The relationship between the POD modes or subspaces at different parameters were analysed using the geometric features of the curves or curved surfaces on the matrix manifolds. The set of POD modes and subspaces were characterised by defining the Riemannian metric, and distance on the Stiefel manifold and Grassmann manifold, respectively. The results obtained from the geometric analysis on the matrix manifolds in this study provided the following insights.
First, the tangent vector along the curve on the Grassmann manifold can be interpreted as the sensitivity of the subspace. This is closely related to the variation in the dynamics of the fluid flow due to changes in the flow parameters. The sensitivity of the subspace, which is spanned by the POD modes for the flow field around a cylinder, increased with decreasing Reynolds number. Our results indicate that the inverse of the subspace sensitivity increased linearly with the Roshko number, especially for higher subspace dimensions. This relationship collapsed into a unified line when the line element of the curve on the Grassmann manifold was normalised by the subspace dimension. In addition, the results obtained in this study indicate that the Reynolds number, at which the inverse of the subspace sensitivity becomes zero, was in good agreement with the lower bound of Reynolds number, where the characteristic frequency of the Kármán vortex street exists (see figure 6). The inverse of the subspace sensitivity with respect to the rotation rate of the cylinder decreases as the rotation rate increases and is approaching zero as the rotation rate approaches the vicinity of the Hopf bifurcation point (see figure 7). These results imply that the flow parameter at which the subspace sensitivity approached infinity corresponded to the parameter at which the properties of the POD modes and fluid flow change significantly. Consequently, the geometric features of the curve on the Grassmann manifold provided insights into the parameter dependence of the fluid-flow dynamics. Additionally, the distribution of the subspaces as a function of the Reynolds number and rotation rate was visualised using the norm of the tangent vector and angle between tangent vectors for two different parameters. The distribution obtained in the tangent vector space whose base point corresponded to
$(Re,\alpha )=(100,0.0)$
indicated that the variation in the subspace to the Reynolds number is not along a geodesic. Instead, the subspace varies along a curved path with non-zero curvature. In contrast, the subspaces were almost aligned with a geodesic with a rotation rate in the range of
$0.0\leqslant \alpha \leqslant 0.8$
.
Second, the sensitivity of the POD modes was represented as a tangent vector on the Stiefel manifold. This enabled the analysis of the flow field sensitivity by superposing the sensitivity modes, which were defined using the tangent vector and sensitivity of the matrices associated with the expansion coefficients. The sensitivity of the POD modes with respect to the Reynolds number showed the displacement of the POD modes in the
$x$
-direction (main streamwise direction). The sensitivity with respect to the rotation rate represented a distribution indicating the displacement of the POD modes in the
$y$
-direction. The sensitivity mode of the flow field with respect to the Reynolds number exhibited a distribution similar to that of the POD mode sensitivity. We found that not only the POD mode sensitivity, but also the expansion coefficient sensitivity played an important role in representing the flow field sensitivity. In particular, the phase shift in the expansion coefficients due to the change in the Reynolds number had a crucial contribution to the flow field sensitivity. Regarding the sensitivity modes with respect to the rotation rate, the sensitivity of the POD modes was predominant when the rotation rate was low, indicating that the influence of sensitivity related to the expansion coefficient was negligible. In contrast, when the rotation rate was high, the phase shift of the expansion coefficients affected the sensitivity of the flow field. We also visualised the spatial distribution of the flow field sensitivity by superposing the obtained sensitivity modes. The spatial distribution of the flow field sensitivity with respect to the Reynolds number showed how the structure of the Kármán vortex street evolved as the Reynolds number varied, that is, the high-sensitivity region appeared along with the structure of the Kármán vortex street. In contrast, the spatial distribution of the flow field sensitivity with the rotation rate indicates that the downstream region exhibited low sensitivity. The flow field around and immediately behind the cylinder showed high sensitivity to changes in the rotation rate. The spatial distribution of the flow field sensitivity with respect to the parameter changes enabled the investigations of the variation in the flow field when deviating slightly from the design point of fluid machinery and the effect of parameter uncertainties on the uncertainties in the flow field at a low computational cost.
Third, the subspace estimation error, which was closely related to subspace sensitivity, influences dominantly the reconstruction error of the parametric ROM based on the subspace interpolation on the Grassmann manifold. This suggested that clarification of the parameter dependence of subspace sensitivity leads to the realisation of a parametric ROM with a small reconstruction error using a limited number of subspace samples. In this study, we compared the performance of parametric ROMs constructed using the direct interpolation method, the global POD method and the manifold interpolation method. Our results demonstrated that, for estimating the subspace at an unseen parameter, the direct interpolation and global POD methods failed to accurately capture the subspace spanned by the reference POD modes, whereas the manifold interpolation method successfully identified the locally optimal subspace corresponding to the target parameter. For the flow field reconstruction, the direct interpolation method failed to reproduce the physically consistent flow field, even when
$\Delta Re$
was small. In contrast, both the global POD and manifold interpolation methods accurately reproduced the flow field at the target parameter, with smaller
$\Delta Re$
resulting in lower reconstruction errors. The difference between the global POD and manifold interpolation methods lies in the number of POD modes required to achieve sufficiently small reconstruction errors. The manifold interpolation method resulted in smaller reconstruction errors with fewer modes, reflecting its ability to capture parameter-dependent, locally optimal subspaces, whereas the global POD method does not necessarily provide the locally optimal subspace for each parameter. We also confirmed that the effectiveness of the manifold interpolation method becomes more pronounced when the parametric ROM is constructed using training data that span a wider parameter range. Furthermore, the reconstruction error of the flow field estimated using the parametric ROM with pseudo-POD modes decreased by decreasing the interval for subspace sampling in both the Reynolds-number and rotation-rate directions. The distribution of the reconstruction error in the parameter space exhibited a trend similar to that of the subspace-estimation error, which was affected by the subspace sensitivity. This result indicates that obtaining an appropriate subspace for the target parameter was essential to improve the performance of the parametric ROM. Specifically, our results implied that finely sampling subspaces in regions with high subspace sensitivity and coarsely sampling in regions with low sensitivity was preferable for efficiently minimising the reconstruction error with a limited number of subspace samples. Also, when interpolating subspaces on the Grassmann manifold, the magnitude of the eigenvalues of each POD mode is disregarded. As a result, in the case where higher-order modes exhibit poor convergence or contain noise, often due to limiting sampling or more complex unsteady dynamics, the subspace interpolation may introduce inaccuracies. Within the range of flow fields and subspace dimensions considered in this study, the developed framework yielded physically consistent results. However, in situations where higher-order modes are not sufficiently converged, exploring extensions such as the order-preserving interpolation method recently proposed by Goutaudier, Nobile & Schiffmann (Reference Goutaudier, Nobile and Schiffmann2023) could provide a promising direction for future research.
These findings were obtained from a geometric analysis of the POD modes on the matrix manifolds with the distance that relates the subspaces for different flow parameters. When discussing the features of fluid-flow dynamics over a wide range of flow parameters (e.g. the initial and boundary conditions of the Navier–Stokes equations), instead of attempting to extract a linear subspace in which the flow fields for the entire parameter space are described, we showed that the subspaces indicating the local features of the fluid flow evolved continuously with changes in the parameters, thereby providing another intuition for understanding fluid dynamics from a global perspective. This perspective provided a framework for analysing experimental and numerical simulation data using the geometric properties of the Grassmann manifold and Stiefel manifold, which were interpreted abstractly. Mode sensitivity analysis and parametric ROM will contribute to the extraction of comprehensive insights into fluid dynamics over a wide range of flow parameters from data. It will also lead to the development of novel techniques for the optimal design of fluid machinery, sensitivity analysis of fluid flows with respect to changes in parameters and active flow control, in which the evaluation of the flow field in a wide range of parameter spaces with significantly low computational costs is crucial.
Funding
This work was supported by JST PRESTO (Grant Number JPMJPR21O4) and JSPS KAKENHI (Grant Number 24K17442).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validity of rotation approximation of semi-orthogonal matrices
Here, we discuss the validity of the approximation in (2.38), which is used to define the sensitivity modes. To evaluate this validity, we demonstrate that the semi-orthogonal matrix consisting of right singular vectors at a specified reference parameter
$V^{\textit{ref}}$
can be represented as the product of the semi-orthogonal matrix
$V(\xi )$
and
$R^{T}(\xi )$
. The orthogonal matrix
$R(\xi )$
is obtained by solving the orthogonal Procrustes problem (Gulub & Loan Reference Golub and van Loan2013) as follows:
\begin{equation} \mathrm{minimise}\,\|V(\xi )-V^{\textit{ref}}R(\xi ) \|_{{F}},\quad \mathrm{subject\,to}\,R(\xi )\in O(r), \end{equation}
where
$\|\boldsymbol{\cdot }\|_{{F}}$
indicates the Frobenius norm. The SVD of
$(V^{\textit{ref}})^TV(\xi )$
is performed to minimise
$R(\xi )$
(A1):
\begin{equation} \big(V^{\textit{ref}}\big)^TV(\xi ) \stackrel {\textrm {SVD}}{=} U_{{P}}\varSigma _{{P}}V_{{P}}^T. \end{equation}
The product of
$U_{{P}}$
and
$V_{{P}}^T$
is the
$R(\xi )$
that minimises (A1).
Figure 26(a–e) show the phase portraits of the trajectories of
$\boldsymbol{v}(t;\xi )\in \mathbb{R}^r$
for different Reynolds numbers ranging from
$100$
to
$160$
when the rotation rate was fixed at
$0.0$
, where
$V^T(\xi )=:\begin{bmatrix}\boldsymbol{v}_1(\xi ) & \ldots & \boldsymbol{v}_{N_t}(\xi ) \end{bmatrix}\in \mathbb{R}^{r\times N_t}$
. The product of
$\boldsymbol{v}(t;\xi )$
and the singular values correspond to the expansion coefficients
$\boldsymbol{a}(t;\xi )\in \mathbb{R}^r$
. The reference parameter is set to
$(Re,\alpha )=(100,0.0)$
. The trajectory varies with the Reynolds number, whereas all trajectories exhibit periodic behaviour (i.e. limit cycle). This indicates that the difference in these periodic trajectories is owing to differences in their phases. The phase portraits of the trajectories determined by the matrix
$V(Re)R^T(Re)$
instead of
$V(Re)$
are shown in figure 26(f–j). These trajectories coincide with that determined by
$V^{\textit{ref}}$
regardless of the Reynolds number, which suggests that the trajectories of
$\boldsymbol{v}$
can be represented as
$V(Re)\approx V^{\textit{ref}}R(Re)$
.
Phase portraits of the trajectories of expansion coefficients normalised by corresponding singular values for different Reynolds numbers at
$\alpha =0.0$
: (a,b,c,d,e) trajectories determined by
$V(Re)$
(the parameter corresponds to the Reynolds number); (f,g,h,i,j) trajectories determined by
$V(Re)R^T(Re)$
. Panels (a,f), (b,g), (c,h), (d,i) and (e,j) show the phase portraits of the normalised expansion coefficients of the 1st–3rd, 1st–5th, 1st–7th, 1st–9th and 1st–11th modes, respectively.
Phase portraits of the trajectories of the normalised expansion coefficients for different rotation rates at
${\textit{Re}}=100$
: (a,b,c,d,e) trajectories determined by
$V(\alpha )$
; (f,g,h,i,j) trajectories determined by
$V(\alpha )R^T(\alpha )$
. Panels (a,f), (b,g), (c,h), (d,i) and (e,j) show the phase portraits of the normalised expansion coefficients of the 1st–3rd, 1st–5th, 1st–7th, 1st–9th and 1st–11th modes, respectively.
Figure 27(a–e) show phase portraits of the trajectory for rotation rates ranging from
$0.0$
to
$1.6$
when the Reynolds number is fixed at
$100$
. The trajectory, which shows a limited cycle, varies with the rotation rate, as in the case of the Reynolds number variation. The trajectories determined by the matrix
$V(\alpha )R^T(\alpha )$
show good agreement with the trajectory of
$V^{\textit{ref}}$
, where the reference parameter is
$(Re,\alpha )=(100,0.0)$
, as shown in figure 27(f–j). Thus, we can conclude that matrix
$R(\xi )$
can be interpreted as representing the phase differences between the trajectories determined by
$V(\xi )$
and
$V^{\textit{ref}}$
, and the matrix
$V(\xi )$
can be approximated as the product of
$V^{\textit{ref}}$
and
$R(\xi )$
. This approximation is valid for the flow parameters considered in this study because the difference in trajectories is caused by the phase difference in these trajectories.
Open access
