Collocation Methods for Exploring Perturbations in Linear Stability Analysis

TL;DR

This paper introduces a stochastic collocation approach to pseudo-spectral analysis for better understanding of instabilities in fluid dynamics, revealing features missed by traditional eigenvalue analysis.

Contribution

It presents a novel, cost-effective method using collocation techniques to analyze perturbations and instability in dynamical systems beyond eigenvalue analysis.

Findings

Stochastic collocation methods can predict transient growth of perturbations.

The approach provides insights into unsteady flow behaviors.

It offers a quantitative tool for stability analysis in fluid dynamics.

Abstract

Eigenvalue analysis is a well-established tool for stability analysis of dynamical systems. However, there are situations where eigenvalues miss some important features of physical models. For example, in models of incompressible fluid dynamics, there are examples where linear stability analysis predicts stability but transient simulations exhibit significant growth of infinitesimal perturbations. This behavior can be predicted by pseudo-spectral analysis. In this study, we show that an approach similar to pseudo-spectral analysis can be performed inexpensively using stochastic collocation methods and the results can be used to provide quantitative information about instability. In addition, we demonstrate that the results of the perturbation analysis provide insight into the behavior of unsteady flow simulations.

Equations

$$\frac{\partial u}{\partial t}=f(u,t),$$

$$\frac{\partial u^{(s)}}{\partial t}=f(u^{(s)},t)=0\text{mbox}forallt.$$

$$M\frac{\partial \mathbf{u}}{\partial t}=\mathbf{f}(\mathbf{u},t),$$

$$J\mathbf{v}=\lambda M\mathbf{v},$$

$$\hat{J}({\mathbf{u}}^{(s)},\mathbf{\delta }(\mathbf{\xi }))\mathbf{v}=\hat{\lambda }(\mathbf{\xi })M\mathbf{v},$$

$$g(\mathbf{\xi }):=\text{mbox}rightmosteigenvalueof(\text{ref}perturbed-eig-linear-stability),$$

$${\mathbf{g}}^{(I)}(\mathbf{\xi }):=\sum _{k=1}^{n_{\mathbf{\xi }}}g({\mathbf{\xi }}^{(k)}){\ell }_k(\mathbf{\xi }),$$

$${\ell }_k({\mathbf{\xi }}^{(\ell )})={\delta }_{k\ell },1\le k,\ell \le m.$$

$${\begin{array}[]{rcl}\!\vec{\,u}_{t}-\nu\nabla^{2}\!\vec{\,u}+\!\vec{\,u}\cdot\nabla\!\vec{\,u}+\nabla p&=\ \!\vec{\,0},\\ -\nabla\cdot\!\vec{\,u}&=\ 0,\end{array}}$$

$$\vec{u}=\vec{w}\text{mbox}on\partial {\mathcal{D}}_D,\nu \frac{\partial \vec{u}}{\partial n}-\vec{n}p=\vec{0}\text{mbox}on\partial {\mathcal{D}}_N,$$

$${\mathcal{H}}_E^1:=\{\vec{u}\in {\mathcal{H}}^1(\mathcal{D}{)}^d|\vec{u}=\vec{w}\mathrm{on}\partial {\mathcal{D}}_D\},{\mathcal{H}}_{E_0}^1:=\{\vec{v}\in {\mathcal{H}}^1(\mathcal{D}{)}^d|\vec{v}=\vec{0}\mathrm{on}\partial {\mathcal{D}}_D\}.$$

$${\begin{array}[]{l}\int_{\mathcal{D}}\!\vec{\,u}_{t}\cdot\!\vec{\,v}+\nu\int_{\mathcal{D}}\nabla\!\vec{\,u}:\nabla\!\vec{\,v}+\int_{\mathcal{D}}(\!\vec{\,u}\cdot\nabla\!\vec{\,u})\cdot\!\vec{\,v}-\int_{\mathcal{D}}p\>(\nabla\cdot\!\vec{\,v})=\int_{\mathcal{D}}{\!\vec{\,f}}\cdot\!\vec{\,v}\\ \hskip 216.81pt\mbox{for all }\!\vec{\,v}\in\mathcal{H}^{1}_{E_{0}}\\ \hskip 72.26999pt-\int_{\mathcal{D}}q\,(\nabla\cdot\!\vec{\,u})=0\quad\mbox{for all }q\in L_{2}(\mathcal{D}).\end{array}}$$

$$(\vec{u}+\vec{\gamma })\cdot \nabla (\vec{u}+\vec{\gamma })=\vec{u}\cdot \nabla \vec{\gamma }+\vec{\gamma }\cdot \nabla \vec{u}+\vec{\gamma }\cdot \nabla \vec{\gamma }\approx \vec{u}\cdot \nabla \vec{\gamma }+\vec{\gamma }\cdot \nabla \vec{u},$$

$$a(\vec{\gamma },\vec{v};\vec{u}):=\nu {\int }_{\mathcal{D}}\nabla \vec{\gamma }:\nabla \vec{v}+{\int }_{\mathcal{D}}(\vec{u}\cdot \nabla \vec{\gamma })\cdot \vec{v}+{\int }_{\mathcal{D}}(\vec{\gamma }\cdot \nabla \vec{u})\cdot \vec{v}.$$

$${\begin{array}[]{l}\int_{\mathcal{D}}[\!\vec{\,u}_{h}]_{t}\cdot\!\vec{\,v}_{h}+\nu\int_{\mathcal{D}}\nabla\!\vec{\,u}_{h}:\nabla\!\vec{\,v}_{h}+\int_{\mathcal{D}}(\!\vec{\,u}_{h}\cdot\nabla\!\vec{\,u}_{h})\cdot\!\vec{\,v}_{h}-\int_{\mathcal{D}}p_{h}\>(\nabla\cdot\!\vec{\,v}_{h})=\int_{\mathcal{D}}{\!\vec{\,f}}\cdot\!\vec{\,v}_{h}\\ \hskip 216.81pt\mbox{for all }\!\vec{\,v}_{h}\in X_{0}^{h},\\ \hskip 72.26999pt-\int_{\mathcal{D}}q_{h}\,(\nabla\cdot\!\vec{\,u}_{h})=0\quad\mbox{for all }q_{h}\in Y^{h}.\end{array}}$$

$${\begin{array}[]{rcll}a(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,u}^{\,(s)}_{h})-\int_{\mathcal{D}}p_{h}(\nabla\cdot\!\vec{\,v}_{h})&=&\lambda_{h}\int_{\mathcal{D}}\!\vec{\,u}_{h}\cdot\!\vec{\,v}_{h}&\hbox{for all}\;\!\vec{\,v}_{h}\in X_{0}^{h},\\ \int_{\mathcal{D}}q_{h}(\nabla\cdot\!\vec{\,u}_{h})&=&0&\hbox{for all}\;q_{h}\in Y^{h}.\end{array}}$$

$$\left[\begin{array}[]{cc}F&B^{T}\\ B&0\end{array}\right]\left[\begin{array}[]{c}{\bf u}\\ {\bf p}\end{array}\right]=\lambda\,\left[\begin{array}[]{cc}-Q&0\\ \ 0&0\end{array}\right]\left[\begin{array}[]{c}{\bf u}\\ {\bf p}\end{array}\right].$$

$${\int }_{\mathcal{D}}{q}_h(\nabla \cdot {\vec{\delta }}_h)=0for\; all{q}_h\in Y^h.$$

$${\vec{\delta }}_h^{(k)}=curl{\varphi }_h^{(k)}={\left[-\partial {\varphi }_h^{(k)}/\partial x_2,\partial {\varphi }_h^{(k)}/\partial x_1\right]}^T,$$

$${\int }_{\mathcal{D}}{q}_h(\nabla \cdot {\vec{\delta }}_h)=\sum _{k\in T_h}{\int }_k{q}_h^{(k)}(\underset{=0}{\underbrace{\nabla \cdot curl{\varphi }_h^{(k)}}})=0.\square$$

$$c_h^{(k)}(u_h^{(k)},v_h^{(k)};{\vec{\delta }}_h^{(k)}):={\int }_k({\vec{\delta }}_h^{(k)}\cdot \nabla u_h^{(k)})v_h^{(k)}-\frac{1}{2}{\int }_{\partial k}{u}_h^{(k)}{v}_h^{(k)}{\vec{\delta }}_h^{(k)}\cdot \vec{n},$$

$$c_h(u_h,v_h;{\vec{\delta }}_h):=\sum _{k\in T_h}{c}_h^{(k)}(u_h^{(k)},v_h^{(k)};{\vec{\delta }}_h^{(k)})$$

$${\begin{array}[]{lcl}\int_{k}(\!\vec{\,\delta}_{h}^{(k)}\!\cdot\!\nabla u_{h}^{(k)})\,v_{h}^{(k)}&=&\int_{k}(v_{h}^{(k)}\!\vec{\,\delta}_{h}^{(k)})\!\cdot\!\nabla u_{h}^{(k)}=-\int_{k}\nabla\cdot(v_{h}^{(k)}\!\vec{\,\delta}_{h}^{(k)}\!)u_{h}^{(k)}+\int_{\partial k}u_{h}^{(k)}v_{h}^{(k)}\;\!\vec{\,\delta}_{h}^{(k)}\!\!\cdot\!\vec{\,n}\\ &=&-\int_{k}(v_{h}^{(k)}\nabla\cdot\!\vec{\,\delta}_{h}^{(k)}+\!\vec{\,\delta}_{h}^{(k)}\!\cdot\!\nabla v_{h}^{(k)})\,u_{h}^{(k)}+\int_{\partial k}u_{h}^{(k)}v_{h}^{(k)}\;\!\vec{\,\delta}_{h}^{(k)}\!\!\cdot\!\vec{\,n}\\ &=&-\int_{k}(\!\vec{\,\delta}_{h}^{(k)}\cdot\nabla v_{h}^{(k)})u_{h}^{(k)}+\int_{\partial k}u_{h}^{(k)}v_{h}^{(k)}\;\!\vec{\,\delta}_{h}^{(k)}\!\!\cdot\!\vec{\,n},\end{array}}$$

$${\int }_k({\vec{\delta }}_h^{(k)}\cdot \nabla u_h^{(k)})v_h^{(k)}-\frac{1}{2}{\int }_{\partial k}{u}_h^{(k)}{v}_h^{(k)}{\vec{\delta }}_h^{(k)}\cdot \vec{n}=-{\int }_k({\vec{\delta }}_h^{(k)}\cdot \nabla v_h^{(k)})u_h^{(k)}+\frac{1}{2}{\int }_{\partial k}{u}_h^{(k)}{v}_h^{(k)}{\vec{\delta }}_h^{(k)}\cdot \vec{n},$$

$${\begin{array}[]{rcll}a(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,u}^{\,(s)}_{h})+c_{h}(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,\delta}_{h})-\int_{\mathcal{D}}p_{h}(\nabla\cdot\!\vec{\,v}_{h})&=&\hat{\lambda}_{h}\int_{\mathcal{D}}\!\vec{\,u}_{h}\cdot\!\vec{\,v}_{h}&\hbox{for all}\;\!\vec{\,v}_{h}\in X_{0}^{h},\\ \int_{\mathcal{D}}q_{h}(\nabla\cdot\!\vec{\,u}_{h})&=&0&\hbox{for all}\;q_{h}\in Y^{h}.\end{array}}$$

$$\left[\begin{array}[]{cc}F+N(\boldsymbol{\xi})&B^{T}\\ B&0\end{array}\right]\left[\begin{array}[]{c}{\bf u}\\ {\bf p}\end{array}\right]=\hat{\lambda}\,\left[\begin{array}[]{cc}-Q&\alpha B^{T}\\ \ \alpha B&0\end{array}\right]\left[\begin{array}[]{c}{\bf u}\\ {\bf p}\end{array}\right],$$

$$(\mathbf{u},N\mathbf{v})=c_h(u_h,v_h,{\vec{\delta }}_h)$$

$$\mathbf{\varphi }(\mathbf{\xi }):=\sigma V{\Theta }^{1/2}\mathbf{\xi }=\sigma \sum _{j=1}^n\sqrt{{\theta }_j}{\mathbf{v}}_j{\xi }_j,$$

$$\mathbf{\varphi }(\mathbf{\xi }):=\sigma \sum _{j=1}^m\sqrt{{\theta }_j}{\mathbf{v}}_j{\xi }_j.$$

$$\mathcal{C}(x^{(1)},x^{(2)}):={\sigma }^2\exp \left(-\left[{\left(\frac{x_1^{(1)}-x_1^{(2)}}{c_1}\right)}^2+{\left(\frac{x_2^{(1)}-x_2^{(2)}}{c_2}\right)}^2\right]\right),$$

Full text

1 Introduction

This study is concerned with a refined understanding of the classic problem of stability of dynamical systems. Let

[TABLE]

represent a dynamical system, where $u:\hbox{{\msbm\char 82\relax}}^{d}\times[0,T]\to\hbox{{\msbm\char 82\relax}}$, $f:\hbox{{\msbm\char 82\relax}}\times[0,T]\to\hbox{{\msbm\char 82\relax}}$, and let $u^{(s)}$ denote a steady solution to (1), i.e.,

[TABLE]

Let $\gamma=\gamma(x,0)$ represent a small perturbation of $u^{(s)}$. Suppose the perturbed quantity $\hat{u}(x,0):=u^{(s)}(x)+\gamma(x,0)$ is taken as an initial condition for (1), for which integration leads to a solution $\hat{u}(x,t)=u^{(s)}(x)+\gamma(x,t)$. If $\hat{u}(x,t)$ reverts to $u^{(s)}(x)$ ($\gamma(x,t)\to 0$) as $t$ increases, then the steady solution is said to be stable; otherwise it is unstable. In typical applications, $f(u,t)=f_{\alpha}(u,t)$ depends on a parameter $\alpha$, as does the resulting steady solution $u^{(s)}_{\alpha}$, and we are interested in the set of values of $\alpha$ for which $u^{(s)}_{\alpha}$ is stable.

Spatial discretization of (1) leads to a discrete version of it, which has the form

[TABLE]

where ${\bf u}$ and ${\bf f}({\bf u},t)$ are finite-dimensional vectors of size $n_{u}$, the size of the spatial discretization. For finite-element discretization, $M$ is a mass matrix. As above, we wish to know if a steady solution ${\bf u}^{(s)}$ to (2) is stable.

Linear stability analysis addresses this question by examining the eigenvalues of the algebraic system

[TABLE]

where $J=\frac{\partial{\bf f}}{\partial{\bf u}}({\bf u}^{(s)})$ is the Jacobian matrix of ${\bf f}$ with respect to ${\bf u}$, evaluated at ${\bf u}^{(s)}$; see, for example, [13, Ch. 1]. A necessary condition for stability of ${\bf u}^{(s)}$ is that all eigenvalues $\lambda$ of (3) have negative real part. If any eigenvalue has positive real part, then there exists an arbitrary small perturbation $\boldsymbol{\gamma}$ such that if ${\bf u}^{(s)}+\boldsymbol{\gamma}$ is used as an initial condition for (2), the integrated solution will not revert to ${\bf u}^{(s)}$.

A problematic aspect of linear stability analysis is that it fails to account for transient effects that may take a long time to resolve. In particular, it may happen that the solution of the system (1) with initial condition $u(x,0)=u^{(s)}(x)+\gamma(x)$, consisting of a small perturbation of a steady solution, exhibits large growth over a significant period of time even if $u^{(s)}$ is linearly stable. This is discussed for models of flow in [21, Sections 2.3,4.1], [23, Sections 20,22]. It can be explained using pseudospectra: the $\epsilon$-pseudospectrum of the Jacobian matrix, defined for $M=I$ in (3), is the set of eigenvalues of $J+E$ for $\|E\|\leq\epsilon$. (A generalization to forms of $M$ considered in the present study is discussed in [11].) Transient growth is exhibited when some elements of this set protrude into the right-half of the complex plane [23].

Our aim in this study is to develop and explore a simple procedure to study the sensitivity of the eigenvalues of (3) when the dynamical system comes from models of incompressible flow. As observed in [11], it is not practical to compute pseudospectra for the large-scale systems that arise in this setting. This difficulty is addressed in [11] by projecting such systems into invariant subspaces of (shifted versions of) $J^{-1}M$, which have smaller dimension and for which computation of pseudospectra is feasible. It is shown in [11] that these pseudospectra estimates provide interior bounds on pseudospectra of (3) as well as insight into transient growth of solutions.

In this work, we develop a complementary approach to study the sensitivity of the eigenvalues of (3) for models of incompressible flow. The methodology derives from a two-fold procedure:

Introduce a simple way to construct perturbed versions of the eigenvalue problem (3) using spatial perturbations that depend on a finite number of parameters. 2. 2.

Approximate the critical eigenvalues of the perturbed problem using a surrogate function defined by interpolation.

This requires the solution of a relatively small number of perturbed eigenvalue problems determined from a special set of parameter values, using sparse-grid methods [1, 22]. The surrogate function interpolates the critical eigenvalues obtained from these eigenvalue problems and provides a means of approximating the critical eigenvalues for an additional set of perturbed problems. The surrogate function is very inexpensive to evaluate. As a result, it is possible to generate many samples of (approximate) eigenvalues in order to gain an understanding of the effects of perturbation. We apply this technique to the eigenvalue problems arising from stability analysis of the incompressible Navier-Stokes equations.

An outline of the remainder of the paper is as follows. In Section 2, we describe the collocation strategy and show in detail how it is developed for the Navier-Stokes equations. In Section 3, we describe two benchmark problems we use to test the methodology and show how the perturbed eigenvalues behave with respect to Reynolds numbers and sizes of perturbation, and in Section 4, we demonstrate that the behavior of perturbed eigenvalues predicts the behavior of transient solutions obtained from perturbed flow conditions. Finally, in Section 5, we make some concluding remarks.

2 Approach

In this section, we describe the methodology we will use to explore the sensitivity to perturbation of the eigenvalue problem (3), which is based on sampling. We first outline the approach in general terms in Section 2.1, and then we continue in Section 2.2 with a more detailed statement of how the ideas are applied to a specific benchmark problem, the incompressible Navier-Stokes equations.

2.1 General approach

Let ${\bf u}^{(s)}$ be a steady solution to (2) and let $\boldsymbol{\delta}$ be a small perturbation of ${\bf u}^{(s)}$. We will specify $\boldsymbol{\delta}=\boldsymbol{\delta}(\boldsymbol{\xi})$ to depend on a vector of parameters $\boldsymbol{\xi}:=(\xi_{1},\xi_{2},\ldots,\xi_{m})^{T}$ with $\boldsymbol{\delta}({\bf 0})={\bf 0}$, and we will explore a perturbed eigenvalue problem

[TABLE]

with the aim of understanding the impact of the perturbation $\boldsymbol{\delta}$ on the eigenvalues $\{\hat{\lambda}\}$. One way to define $\hat{J}$ is to evaluate the Jacobian at the perturbed velocity, $\hat{J}({\bf u}^{(s)},\boldsymbol{\delta}):=J({\bf u}^{(s)}+\boldsymbol{\delta})$. In this study, which concerns the incompressible Navier-Stokes equations, we will insist that the perturbation is not dissipative. Details on the structure of the perturbation and its parameter dependence are given in Section 2.2.

Remark 2.1. We call attention here to an important aspect of the issue under study. Classic linear stability analysis concerns the sensitivity of the steady solution ${\bf u}^{(s)}$ to perturbation. Our (different) concern here, like that of [23], is the sensitivity of the eigenvalues $\lambda$ to perturbation, and in particular whether the conclusions reached from stability analysis predict behavior. To highlight this distinction, we use different symbols for perturbation depending on context: $\gamma$ is used for perturbations arising in linear stability analysis, and $\delta$ for perturbations of eigenvalue problems as in (4).

Given the eigenvalue problem (4), let

[TABLE]

where, if there is a complex conjugate pair of rightmost eigenvalues, $g(\boldsymbol{\xi})$ can be taken to be the eigenvalue with positive imaginary part. One way to explore the sensitivity of (3) is by sampling $\boldsymbol{\xi}$, that is, to evaluate $g(\boldsymbol{\xi})$ for a large set of sample values of $\boldsymbol{\xi}$. If this function is very sensitive, that is, if small changes in $\boldsymbol{\delta}(\boldsymbol{\xi})$ lead to large changes in $g(\boldsymbol{\xi})$, then linear stability analysis may not provide an accurate assessment of stability; conversely, if $g$ is not sensitive to perturbation, then linear stability analysis is likely to yield insight.

The point of view here is that the study of perturbation is done by sampling a large number of nearby problems. A potential downside is that this approach requires the solution of many eigenvalue problems (4), one for each choice of $\boldsymbol{\xi}$ and resulting $\boldsymbol{\delta}(\boldsymbol{\xi})$, which tends to incur a high computational cost. To reduce this expense, instead of evaluating the function of (5) (by solving an eigenvalue problem), we will replace $g(\boldsymbol{\xi})$ with an approximation, a surrogate function $g^{(I)}(\boldsymbol{\xi})$, which is inexpensive to compute and therefore can be evaluated cheaply for many samples of $\boldsymbol{\xi}$. For this, we will use the method of collocation designed to construct approximations to functions on high-dimensional spaces [1, 22]. This entails evaluation of $g(\boldsymbol{\xi})$ at a relatively small number of special points, $\{{\bf\boldsymbol{\xi}}^{(1)},{\bf\boldsymbol{\xi}}^{(2)},\ldots,{\bf\boldsymbol{\xi}}^{(n_{\boldsymbol{\xi}})}\}$. The surrogate function is then taken to be the polynomial interpolant of $g$,

[TABLE]

where $\{\ell_{k}(\boldsymbol{\xi})\}$ are multidimensional Lagrange interpolation polynomials,

[TABLE]

For the interpolation points, we use sparse grids derived from the extrema of one-dimensional Chebyshev polynomials [1].

Remark 2.2. It might happen that there are multiple eigenvalues of (3) with the same rightmost real part and different imaginary parts. In this case, $g(\boldsymbol{\xi})$ of (5) would also be multi-valued or nearly so, and the ideas presented here would need to be applied to each of the rightmost eigenvalues. As long as there are not too many such values, this would have minimal impact on costs.

2.2 Application to the Navier–Stokes equations

We will explore these ideas when the dynamical system (1) comes from the incompressible Navier–Stokes equations, and we now describe a way to specify a perturbation $\boldsymbol{\delta}(\boldsymbol{\xi})$ for this benchmark problem for use in (4). To this end, consider the Navier–Stokes equations

[TABLE]

posed on a domain $\mathcal{D}\subset\hbox{{\msbm\char 82\relax}}^{d}$, $d=2$ or $3$, with boundary conditions

[TABLE]

for $\partial\mathcal{D}=\partial\mathcal{D}_{D}\cup\partial\mathcal{D}_{N}$ consisting of the portions of the boundary of $\mathcal{D}$ on which Dirichlet or Neumann boundary conditions hold. In a typical scenario (see [10, p. 413]), $\!\vec{\,w}$ is a time-dependent inflow function that rapidly goes to a steady state, and the Neumann boundary condition is applied at an outflow boundary. Let $\mathcal{H}^{1}(\mathcal{D})$ be the Sobolev space of functions on $\mathcal{D}$ with first derivatives in $L_{2}(\mathcal{D})$, and let

[TABLE]

For fixed time $t\in(0,\infty)$, the weak formulation of (7) is to find ${\!\vec{\,u}(\cdot,t)}\in\mathcal{H}^{1}_{E}$, ${p(\cdot,t)}\in L_{2}({\mathcal{D}})$ such that

[TABLE]

Linear stability analysis uses a linearized form of the first (momentum) equation of (7)–(8). Given a steady velocity field $\!\vec{\,u}$ (i.e., $\!\vec{\,u}_{t}=0$), consider a perturbation $\!\vec{\,u}+\!\vec{\,\gamma}$. Substitution of this perturbed velocity into the quadratic term from (7) gives

[TABLE]

where the approximation on the right is made under the assumption that $\!\vec{\,\gamma}$ is small. Addition of the diffusion operator and specification of a perturbed weak formulation then leads to a trilinear form associated with the linearized convection-diffusion operator,

[TABLE]

Mixed finite-element discretization of (8) uses finite-dimensional subspaces $X_{0}^{h}$ $\subset$ ${\mathcal{H}}_{E_{0}}^{1}$ and $Y^{h}\subset L_{2}(\mathcal{D})$ together with $X_{E}^{h}\subset\mathcal{H}_{E}^{1}$ containing functions that interpolate the Dirichlet boundary data at element nodes lying in $\partial\mathcal{D}_{D}$. We will assume that this discretization is div-stable [15, Sect. 2.2]. The discrete weak formulation is to find $\!\vec{\,u}_{h}\in X_{E}^{h}$ and $p_{h}\in Y^{h}$ such that

[TABLE]

Let $\!\vec{\,u}^{\,(s)}_{h}$ be a discrete steady solution to (10), i.e., $[\!\vec{\,u}^{\,(s)}_{h}]_{t}=0$. The eigenvalue problem (3) is derived from a linearized discrete formulation associated with (10) where the aim is to find eigenvalues $\lambda_{h}$ and associated eigenfunctions satisfying

[TABLE]

Here, we have linearized around a steady flow velocity field $\!\vec{\,u}^{\,(s)}_{h}$ satisfying (10).

Remark 2.3. A complete discussion of the development of the trilinear form $a(\cdot,\cdot;\cdot)$ of (9) and the derivation of (11) is given in [10, Sections 8.2–8.3]. This form also arises from use of Newton’s method for solving the nonlinear system of equations arising from implicit time discretization of (8).

Let the dimensions of $X_{0}^{h}$ and $Y^{h}$ be $n_{u}$ and $n_{p}$, respectively. Let ${\bf u}^{(s)}$ be the vector of coefficients of the steady finite-element solution $\!\vec{\,u}_{h}^{(s)}$ appearing in (11). Then the eigenvalue problem (3) has the structure

[TABLE]

Here, $F=F({\bf u}^{(s)})$ is the matrix of order $n_{u}$ derived from the bilinear form $a(\cdot,\cdot;\!\vec{\,u}_{h}^{(s)})$, $B$ and $B^{T}$ are matrix representations of negative-divergence and gradient operators, respectively ($B$ is of size $n_{p}\times n_{u}$), and $Q$ is a velocity mass matrix, also of order $n_{u}$.

Remark 2.4. The matrix on the right side of (12) is singular, and the resulting infinite eigenvalue can lead to instability in eigenvalue computations [19]. This can be avoided by replacing the matrix by \left[\!\!{\begin{array}[]{cc}-Q&\alpha B^{T}\\ \ \alpha B&0\end{array}}\!\!\right], which leaves the finite eigenvalues intact and maps the infinite eigenvalue to $1/\alpha$, see [4].111We use this variant of the mass matrix with $\alpha=-1/10$ in all of our computations. The resulting mapped eigenvalue $\lambda=-10$ is far enough from the near-critical ones that it does not affect the results.

We want to explore the sensitivity of our modified eigenvalue problem to perturbation. For this, we add a small perturbation $c(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,\delta}_{h})$ to $a(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,\delta}_{h})$ in (11). The perturbation is defined in two steps. First, we specify a discretely divergence-free vector field $\!\vec{\,\delta}_{h}$, that is, one satisfying

[TABLE]

This ensures that the perturbed velocity $\!\vec{\,u}_{h}+\!\vec{\,\delta}_{h}$ would be appropriate as an initial condition for testing the stability of $\!\vec{\,u}_{h}$. Second, the field $\!\vec{\,\delta}_{h}$ is used to generate a nondissipative perturbation $c_{h}(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,\delta}_{h})$ of $a_{h}(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,\delta}_{h})$. This means that the perturbation does not introduce any damping effects associated with numerical diffusion. To illustrate the construction, we will suppose that $T_{h}$ denotes a subdivision of $\mathcal{D}\subset\hbox{{\msbm\char 82\relax}}^{2}$ into triangular or rectangular elements. The extension to three-dimensional problems is perfectly straightforward.

Claim 2.1. Suppose that $\phi_{h}\in\mathcal{H}^{1}(\mathcal{D})$ is a finite element function defined on $T_{h}$ and that $\!\vec{\,\delta}_{h}$ is defined locally on every element $k\in T_{h}$ , $\!\vec{\,\delta}_{h}|_{k}:=\;\!\vec{\,\delta}_{h}^{\,(k)}$, via

[TABLE]

so that $\!\vec{\,\delta}_{h}^{\,(k)}$ is divergence-free on each element. Then $\!\vec{\,\delta}_{h}:=\sum_{k\in T_{h}}\!\vec{\,\delta}_{h}^{\,(k)}$ satisfies (13).

Proof. For any function $q_{h}\in Y^{h}$,

[TABLE]

Note that the local construction (14) generates a discontinuous velocity field so that in general $\!\vec{\,\delta}_{h}\not\in{\mathcal{H}^{1}(\mathcal{D})}^{d}$.

Claim 2.2. Let the perturbation operator on element $k\in T_{h}$ be given by

[TABLE]

Then $c_{h}^{(k)}$ is skew-adjoint on $T_{h}$, and

[TABLE]

is skew-adjoint on $\mathcal{D}$.222In this discussion, $u_{h}$ and $v_{h}$ are discrete scalar functions. For vector-valued arguments, e.g., $\!\vec{\,u}_{h}=([u_{h}]_{1},[u_{h}]_{2})^{T}$, $c_{h}(\!\vec{\,u}_{h},\!\vec{\,v}_{h};\!\vec{\,\delta}_{h})=\sum_{i}c_{h}([u_{h}]_{i},[v_{h}]_{i};\!\vec{\,\delta}_{h})$ is the sum of contributions from individual scalar components.

Proof. Since $\!\vec{\,\delta}_{h}\not\in{\mathcal{H}^{1}(\mathcal{D})}^{d}$, we need to apply Green’s theorem on each element:

[TABLE]

where the last equality follows from the fact that $\!\vec{\,\delta}_{h}^{(k)}$ is divergence-free on each element. It follows that

[TABLE]

that is, $c_{h}^{(k)}$ is skew-adjoint. Summation over all the elements establishes the same property for $c_{h}$. $\qquad\square$

The perturbed variant of (11) is

[TABLE]

This leads to the perturbed matrix eigenvalue problem (4)

[TABLE]

where the perturbation matrix $N=N(\boldsymbol{\xi})$ is determined from

[TABLE]

so that in particular $N$ is a skew-symmetric matrix, $N^{T}=-N$, for all parameter values $\boldsymbol{\xi}$ independent of the boundary conditions of the flow problem.

It remains to specify the finite element function $\phi_{h}$ used in Claim 2.1 to define the vector field $\!\vec{\,\delta}_{h}$. Following [20], we take $\phi_{h}(x,\boldsymbol{\xi})\in\mathcal{H}^{1}(\mathcal{D})$ to be a parameter-dependent scalar potential specified using a covariance function $\mathcal{C}(x^{(1)},x^{(2)})$ for $x^{(i)}\in\mathcal{D}$. In particular, given $\mathcal{C}$, let $C:=\mathcal{C}({\bf x},{\bf x})$ be the covariance matrix of order $n$ consisting of the vertices in the subdivision associated with $X_{0}^{h}$, so that $C_{ij}=\mathcal{C}\left(x_{i},x_{j}\right)$. Now let $\boldsymbol{\phi}$ be an $n$-dimensional zero-mean stationary random process with covariance matrix $C$, i.e., $C=\hbox{{\msbm\char 69\relax}}(\boldsymbol{\phi}\boldsymbol{\phi}^{T})$, where “𝔼” refers to expected value. If $C=\sigma^{2}V\Theta V^{T}$ is an eigenvalue–eigenvector decomposition (scaled by the variance), then $\boldsymbol{\phi}$ can be defined using a discrete Karhunen–Loève (KL) expansion

[TABLE]

where the eigenvector ${\bf v}_{j}$ is the $j$th column of $V$ and $\{\xi_{j}\}_{j=1}^{n}$ are uncorrelated random variables with zero mean and unit variance [18, Section 5.4]. It is often the case that many of the eigenvalues are small and some of the terms in (18) can be removed without significant loss of accuracy. We will choose $m<n$ such that $\left(\sum_{j=m+1}^{n}\theta_{j}\right)\left/\left(\sum_{j=1}^{n}\theta_{j}\right)\right.\leq 5/100$, and, in the sequel, $\boldsymbol{\xi}:=(\xi_{1},\ldots,\xi_{m})^{T}$ will represent an $m$-dimensional vector of parameters and $\boldsymbol{\phi}(\boldsymbol{\xi})$ is defined using the truncated KL-expansion

[TABLE]

This $\boldsymbol{\xi}$-dependent coefficient vector (of length $n$) then characterizes a piecewise-defined linear or bilinear function $\phi_{h}(\boldsymbol{\xi})$. For the computational results described in Section 3, we take the smooth covariance function

[TABLE]

where $c_{1}$ and $c_{2}$ are correlation lengths. We will also assume that $\{\xi_{j}\}$ in (19) are mutually independent, with each satisfying a truncated Gaussian distribution with range $[-3,3]$, so that $\xi_{j}$ has the density function

[TABLE]

A Matlab implementation of this distribution is given in [3] and described in [2].

We note two differences between this formulation of perturbations and traditional approaches based on pseudospectra. First, the perturbed eigenvalue problem (4), as specified in (16), is restricted to have a structure determined from that of the original problem, so that the resulting perturbed eigenvalues are closer in form to structured pseudoeigenvalues [23, Ch. 50]. Moreover, the perturbation itself derives explicitly from the nonlinear term $\!\vec{\,u}\cdot\nabla\!\vec{\,u}$ in (7). Indeed, if the original problem (1) were linear so that the Jacobian $J$ did not depend on ${\bf u}$, then the eigenvalues of $J$ might still be sensitive to perturbation but the ideas discussed here would not give insight into this. Despite these limitations, as will be shown in Sections 3-4, the behavior of the perturbed eigenvalues gives insight into transient growth of solutions and other features of transient solvers.

We conclude this section with an analytic result bounding the size of the eigenvalue perturbation in proportion to the perturbation size.

Theorem 1**.**

If $\delta_{h}$ is the perturbation defined using (14) and $\lambda$ is the rightmost eigenvalue of (3), then the eigenvalue $\hat{\lambda}$ of the perturbed problem (4) closest to $\lambda$ satisfies $|\hat{\lambda}-\lambda|\leq(c/h)\,\|\delta_{h}\|_{\infty}$.

Proof.

Let $\mathcal{M}$ denote the matrix on the right side of (16), and let $\mathcal{F}$ and $\widehat{\mathcal{F}}=\mathcal{F}+\mathcal{N}$ denote the unperturbed and perturbed matrices on the left sides of (12) and (16), respectively; here \mathcal{N}=\mathcal{N}(\boldsymbol{\xi})=\left[\begin{array}[]{cc}N(\boldsymbol{\xi})&0\\ 0&0\end{array}\right]. We are interested in the eigenvalue problems $\mathcal{F}v=\lambda\mathcal{M}v$ and $\widehat{\mathcal{F}}\hat{v}=\hat{\lambda}\hat{\mathcal{M}}v$. $\mathcal{M}$ can be factored as

[TABLE]

The velocity mass matrix $Q$ and Schur complement $\alpha^{2}BQ^{-1}B^{T}$ each admit a Cholesky decomposition, $Q=LL^{T}$, $\alpha^{2}BQ^{-1}B^{T}=RR^{T}$, so that (21) can be refined to

[TABLE]

Thus, we can consider the unperturbed and perturbed standard eigenvalue problems

[TABLE]

Let $\mathcal{A}=\mathcal{D}^{-1}\mathcal{L}^{-1}\mathcal{F}\mathcal{L}^{-T}$ and $\mathcal{E}=\mathcal{D}^{-1}\mathcal{L}^{-1}\mathcal{N}\mathcal{L}^{-T}$. If $\mathcal{A}=\mathcal{V}\Lambda\mathcal{V}^{-1}$ is diagonalizable, then we can use the Bauer-Fike Theorem [12, Ch. 7] to explore the eigenvalue perturbation: given an eigenvalue $\lambda$ of $\mathcal{A}$,

[TABLE]

We seek a bound on $\|\mathcal{E}\|_{2}=\|\tilde{\mathcal{E}}\|_{2}$, where $\tilde{\mathcal{E}}=\mathcal{L}\mathcal{N}\mathcal{L}^{-T}$; the equality here follows from the fact that $\mathcal{D}$ is unitary. The structure of $\mathcal{L}$ leads to

[TABLE]

so that $\|\tilde{\mathcal{E}}\|_{2}\leq\|L^{-1}NL^{-T}\|_{2}+O(\alpha)$. This holds for all $\alpha\neq 0$, and so it also holds in the limit as $\alpha\to 0$, giving $\|\tilde{\mathcal{E}}\|_{2}\leq\|L^{-1}NL^{-T}\|_{2}$. Since both $N$ and $L^{-1}NL^{-T}$ are skew-symmetric, it follows that $\|\mathcal{E}\|_{2}\leq\rho(L^{-1}NL^{-T})$, the spectral radius.

Thus, we require a bound on the Rayleigh quotient $\frac{|(y,Ny)|}{(y,Qy)}$. For this, we use (17) together with a standard bound on $|c_{h}|$ (see [10, p. 243]) to get

[TABLE]

The assertion then follows from the inverse estimate $\|\nabla y_{h}\|_{L_{2}(\mathcal{D})}/\|y_{h}\|_{L_{2}(\mathcal{D})}\leq O(1/h)$.

∎

3 Benchmark problems and structure of eigenvalues

We will illustrate these ideas for two benchmark problems. In this section, we specify the problems and their features of interest, the eigenvalues associated with linear stability analysis and the effect of perturbation of these eigenvalues. Each of these is a model of flow through a channel for which there are inflow and outflow boundaries. The position of the outflow boundary is far enough downstream that the flow is fully developed. The spatial approximation is done using $Q_{2}$–$P_{-1}$ (biquadratic velocity; discontinuous linear pressure) mixed approximation [10, Section 3.3.1], implemented in the ifiss software package [8, 9]. Unless otherwise specified, the discretization is done on a uniform grid with element width $h=1/32$, which gives $64$ elements along the vertical interval $[-1,1]$. This corresponds to “grid level” $\ell=6$ in ifiss, with element width $h=2/2^{\ell}$. For both benchmark problems, we will explore the stability of solutions obtained for choices of the viscosity parameter near the critical values for linear stability. We determined these critical values experimentally, as described for example in [7].

3.1 Expansion flow around a symmetric step

The domain $\mathcal{D}$ is a rectangular duct with a symmetric expansion, with boundary conditions

• •

parabolic profile $\!\vec{\,u}(-1,y)=(1-4y^{2},0)$ at the inflow boundary $(-1,y),\,|y|\leq.5$

• •

natural conditions $\nu\frac{\partial u_{x}}{\partial x}=p$, $\frac{\partial u_{y}}{\partial x}=0$ at the outflow boundary $(20,y),\,|y|\leq 1$

• •

no-flow conditions $\!\vec{\,u}=0$ along fixed walls

$(x,\pm 1),0\leq x\leq 20$; $(x,\pm.5),-1\leq x\leq 0$; $(0,y)$, $.5\leq|y|\leq 1$.

Details of the domain and a sample solution are shown in Figure 1. The discretization is defined on a uniform grid of square elements. The key feature of this solution is that it is reflectionally symmetric with respect to the centerline $y=0$, i.e., the stream function $\psi$ satisfies $\psi(x,y)=-\psi(x,-y)$. It follows that for the velocity,

[TABLE]

This flow problem exhibits a pitchfork bifurcation [6]: as the viscosity decreases through a critical value (approximately $\nu=1/220.5$), the rightmost eigenvalue of (12), which is real, changes from negative (indicating linear stability) to positive (instability). Figure 2 shows the ten rightmost eigenvalues for three values of $\nu$ in this range and the rightmost eigenvalues (detail in the inset) for each choice, whose values are also identified on the right.

We explored the sensitivity of the rightmost negative eigenvalues using the perturbed eigenvalue problem (16). This was derived using correlation lengths $c_{1}=c_{2}=2$ in (20), which resulted in a finite expansion (19) with $m=19$ terms. The surrogate function $g^{(I)}$ of (6) used to estimate eigenvalues was constructed from a two-level sparse grid on the $m$-dimensional parameter space, which in turn resulted in $n_{\boldsymbol{\xi}}=761$ sparse grid nodes. Thus, it is necessary to solve $760$ eigenvalue problems, that is, find the rightmost eigenvalues of $760$ perturbed systems (16), one for each sparse-grid node. (One of the sparse-grid nodes is $\boldsymbol{\xi}={\bf 0}$, which corresponds to an unperturbed system.) Once these are available, the estimates of eigenvalues for other choices of $\boldsymbol{\xi}$ can be obtained by evaluating $g^{(I)}$. We implemented the sparse-grid interpolation using the matlab toolbox spinterp [16, 17].

The dependence of the estimated perturbed eigenvalues on the perturbation size is illustrated in Figure 3. We found that insight can be provided using small values of the standard deviation in (19), and in these tests we used $\sigma=\beta h^{2}$ for $\beta=.1$, $.2$ and $.3$. For each $\sigma$, the figure shows the distribution of one million eigenvalue estimates, computed

using the interpolant (6), with results for $\nu=1/210$ shown on the left and for $\nu=1/220$ on the right. For both values of $\nu$, the rightmost unperturbed eigenvalue (center of the set of perturbations) is negative, showing that the associated steady solution is stable, but for the smaller, closer-to-critical value $\nu=1/220$, some of the estimated perturbed eigenvalues are positive, whereas all the perturbations are negative for $\nu=1/210$. The two figures have the same horizontal scaling, indicating that the magnitude of the perturbations does not depend on $\nu$. The bounding dashed lines show that the magnitude of perturbations varies linearly with $\sigma$, as the bound of Theorem 1 suggests.

Finally, Figure 4 shows the behavior of the eigenvalue perturbations as well as the critical eigenvalues (highlighted in the middle of each set of perturbations) as the discretization mesh size varies. Results are shown for four mesh sizes, $h=1/8$, $1/16$, $1/32$ and $1/64$ (corresponding to grid levels $\ell=4$ through $7$). For each $\ell$, $d_{\ell}$ is the maximal difference $|\hat{\lambda}-\lambda|$ between the surrogate rightmost perturbed eigenvalue and the rightmost true eigenvalue among all surrogate values. The range of perturbations appears to be independent of the discretization, which suggests that the $O(1/h)$ dependence of the bound in Theorem 1 is pessimistic. We attribute this to limitations on what can be obtained using the Bauer-Fike Theorem for the transformed problems (22). Also note that the critical eigenvalues move to the left with mesh refinement, indicating that stability of the discrete systems is enhanced with mesh refinement, although it is also clear that a limiting value is approached with refinement.

3.2 Flow around a square obstacle

In this case, the domain is a rectangular duct containing a square obstacle, with boundary conditions

• •

parabolic profile $\!\vec{\,u}(-1,y)=(1-4y^{2},0)$ at the inflow boundary $(0,y),\,|y|\leq 1$

• •

natural conditions $\nu\frac{\partial u_{x}}{\partial x}=p$, $\frac{\partial u_{y}}{\partial x}=0$ at the outflow boundary $(8,y),\,|y|\leq 1$

• •

no flow conditions $\!\vec{\,u}=0$ along the top and bottom walls $(x,\pm 1),0\leq x\leq 8$, and

on the obstacle, a square centered at $(2,0)$ with sides of length $0.5$.

For this example, we used a level-6 stretched grid with local refinement near the obstacle. A representative steady solution that retains the reflectional symmetry is shown in Figure 5. In this case, there is a symmetry-breaking Hopf bifurcation for $\nu\approx 1/186$; that is, for $\nu$ in this range there is a complex conjugate pair of rightmost eigenvalues whose real parts change from negative to positive as $\nu$ is reduced. Figure 6 shows the $100$ smallest eigenvalues for three values of $\nu$, two near critical ($\nu=1/175$ and $1/185.6$ and one super-critical ($\nu=1/200$), as well as a detail of the rightmost eigenvalues.

The behavior of the perturbed (estimated) eigenvalues is illustrated in Figure 7. Once again, we computed one million eigenvalue estimates for three values of the standard deviation in (19), $\sigma=.1h^{2}$, $.2h^{2}$ and $.3h^{2}$ with $h=1/32$. These are shown in the figure, for $\nu=1/175$ on the left and $\nu=1/185.6$ on the right. For both values of $\nu$, the rightmost unperturbed eigenvalue (center of the set of perturbations) has negative real part, showing that the associated steady solution is stable, but for the close-to-critical value $\nu=1/185.6$, some of the perturbed eigenvalues have a positive real part. As for the step problem, it is readily seen that the magnitude of the perturbations does not depend on $\nu$ and varies linearly with $\sigma$.

4 Unsteady flow simulations

In this section, we explore the connection between time integration of the Navier–Stokes equations and the eigenvalue perturbation results in the previous section. We will do this by computing time-accurate solutions of the Navier–Stokes equations using the adaptive (stabilized) Trapezoidal Rule (sTR) time stepping methodology built into ifiss. The suitability of sTR for long-time integration is discussed in [14]. Full details of the ifiss  implementation of sTR can be found in section 10.2.3 of [10]. (Stabilization is based on time step averaging, which prevents the “ringing” to which TR is susceptible for stiff systems.) In what follows, we present results obtained from a nonlinear version of the integrator, denoted (sTR$k$), where a fixed number ($k=1$ or $k=2$) of Picard corrections are performed at every time step. We present results for the benchmark problems of Sections 3.1 and 3.2. Our objective is to test the sensitivity of the reference flow with respect to instantaneous spatial perturbations, loosely simulating a laboratory experiment where a reference steady flow is subject to an external disturbance, and the flow is monitored to see if it returns to the steady state.

4.1 Evolution of expansion flow around a symmetric step

Motivated by the eigenvalue calculations shown in Section 3.1, we consider the three distinct values of the viscosity parameter $\nu=1/210$ (linearly stable), $\nu=1/220$ (close to critical) and $\nu=1/250$ (unstable).

We model the laboratory scenario computationally via a two-stage process.

We start from a quiescent state and a tiny time step (1e-9). The inflow profile is smoothly ramped up to a fully developed flow using an exponential startup. The sTR2 integration is then carried out for 330 time steps with a relatively tight accuracy tolerance (i.e., a bound on an estimate of local truncation error), 3e-5. The number of steps taken is arbitrary but needs to be chosen large enough so that the reference flow is visually steady. More precisely, when this phase is complete, the instantaneous acceleration $a(t)$, defined in terms of the flow velocity $\!\vec{\,u}_{h}(\cdot,t)$ at time $t$ by $a(t)=\sqrt{\int_{\mathcal{D}}\left(\frac{\partial\vec{u}_{h}}{\partial t}\right)^{2}}\,,$ should be around $10^{-2}$ or even less.

At the point in time, $T$ say, where the first stage is completed, the integration is interrupted and a perturbation is added to the flow field $\vec{u}_{h}(\cdot,T)$. The perturbation is of the form $\!\vec{\,\delta}_{h}$ specified in Claim 2.1 where the associated scalar field $\phi_{h}$ derives from (19)–(20). We construct $\phi_{h}$ with $\sigma=.3(h^{2})$.333It is also necessary to scale the perturbation so as not to “shock” the transient simulation — the perturbation field is thus scaled by a factor of 1e-5. This ensures that the magnitude of the perturbation is comparable to the time accuracy used for the simulation.

• 2.

The time integration is then restarted without reducing the time step, using sTR1 in place of sTR2 (because it is marginally less dissipative). The restarted integration is continued for a fixed number (typically, 200 or 700) time steps, stopping prematurely only if the time reaches $T^{*}$=1e14 — which we interpret as reaching a “computational steady-state” — at which point the adaptive time-stepping routine is taking very large time steps, see Figure 12 below, and the acceleration $a(T^{*})$ will almost certainly be smaller than unit roundoff.

The mean vorticity $\omega(t)$, or the average vertical velocity at the outflow,

[TABLE]

provides a convenient way of assessing the degree of departure from the reflectionally symmetric flow (for which $\omega$=0). At the conclusion of the first phase of the time integration, a pseudo-steady flow is computed for each of the three values of $\nu$. The evolution of the mean vorticity and the acceleration visualized in Figure 8 shows that a symmetric flow is established for $\nu=1/250$ before the interruption is made after 330 time steps; this corresponds to time $T\approx 32$.

Moving on to the second stage, we show results for the subcritical cases $\nu=1/210$ and $1/220$, for three representative flow perturbations, each of which derives from a particular collocation point $\boldsymbol{\xi}^{(k)}$ used in (6). For the first of these, no perturbation is made (this corresponds to the point $\boldsymbol{\xi}\equiv{\bf 0}$) and the integration simply continues from the first-stage stopping point $T$. The other two are representatives of a “benign” perturbation and a “lively” perturbation and have the spatial structure shown in Figure 9.

The evolution of the flow after the restart for $\nu=1/210$ is depicted in Figure 10, where mean vorticity is shown in both actual and logarithmic scales. The unperturbed flow is perfectly stable; the sTR1 integrator reaches the end time ($T^{*}$=1e14) at time step 396, 66 steps after the restart. The distinctive jumps in the acceleration are associated with the stabilization of the integrator, which has the effect of periodically injecting a small amount of dissipation into the flow. The benign perturbation, which respects the reflectional symmetry, has no effect on the long-term flow evolution. In contrast, the lively perturbation excites visible instability at about time step 390, 60 steps after the restart. But (as seen in particular from the acceleration), the size of the perturbation is not big enough to stop the long-term evolution to the symmetric flow at the designated end time $T^{*}$. The growth in vorticity for the stable examples toward the end of the simulation (about step 385) is a roundoff effect caused by allowing the simulation to proceed after changes in the steady solution are near machine precision.

The evolution of flow for the intermediate viscosity parameter $\nu$ = 1/220 is shown in Figure 11. The unperturbed case is just about stable: the sTR1 integrator reaches the end time 75 steps after the restart. A virtually identical evolution is evident when the perturbation to the flow is benign. (As in the previous example, roundoff effects lead to some growth in the vorticity after a steady solution is obtained.) The time evolution for the lively perturbation is noticeably different, however. In this case, the sTR1 integrator rejects time step 415 (85 steps after the restart) and the computational flow evolves to a numerically noisy solution where the magnitude of the oscillation is of the order of the time-stepping accuracy.

These observations are substantiated in Figure 12, which shows the history of the time step sizes chosen by the adaptive integrator. For each of the plots in this figure, the switch from the first to the second stage is identified by a vertical dotted line. When either no perturbation or a benign perturbation is made, the time step sizes rapidly increase because the integration goes to a steady state for the subcritical values of $\nu$. This behavior can also be seen for the lively perturbation and $\nu=1/210$. In contrast, the integrator behaves differently for $\nu=1/220$ — here the time step size is cut back at around 70 time steps after the perturbation is made in order to resolve the nonstationary solution shown at the bottom of Figure 11. Computing solutions when so close to the stability limit is a delicate business.

Results for the super-critical viscosity parameter $\nu$ = 1/250 are in Figure 13. In this instance, no perturbations are needed to excite instability. The time step history of the complete flow evolution from $t$=0 to $t$=1e14 is presented at the top of the figure. Note the scale on the vertical axis — this is a pretty demanding computational exercise! The evolution of mean vorticity and acceleration after the interrupt is shown in the two plots at the bottom of Figure 13 and should be contrasted with the results for the subcritical viscosity shown in Figure 11. Just when the symmetric flow looks to be steady (400 time steps; 70 after the restart) the time step is cut back to O(1) and after a transient the flow goes to a computational steady state that does not have reflectional symmetry. This is evident from the flow snapshots plotted at/after the interrupt shown in Figure 14; the particular steady-state solution (top eddy longer than the bottom one) is solely determined by the build-up of roundoff error. The two “cups” between 400 and 600 in the time step history shown at the top of Figure 13 suggest that the sTR1 algorithm needed two attempts to fix on the specific stationary solution — it is instructive to contrast this with the evolution that results when vigorously perturbing the flow close to the critical viscosity, which is shown at the bottom of Figure 12.

4.2 Evolution of flow around an obstacle

Motivated by the eigenvalue calculations discussed in Section 3.2, we now consider three distinct values of the viscosity parameter for the obstacle problem: $\nu=1/175$ (subcritical), $\nu=1/185.6$ (close to critical) and $\nu=1/200$ (unstable). We consider $\nu=1/200$ first. The same two-stage process described above gives the results shown in Figure 15. These results should be compared with those in Figure 13. The difference is that instead of going to a nonsymmetric steady state solution, the computational flow evolves to a periodic (vortex-shedding) solution, at which point the time step becomes essentially constant. The vortex-shedding solution is persistent — it is unchanged when we run the solver for another 10,000 time steps. The same long-term behavior is obtained if a perturbation is added at the interrupt point. The different outcomes for the two benchmark problems are representative of the difference between a pitchfork bifurcation (for the step) and a Hopf bifurcation (for the obstacle) [5],[10, p. 343].

To study the flow breakdown mechanism in detail, the second phase of the time integration is computed with a very small accuracy tolerance (1e-9) using the unstabilised TR1 integrator.444Stabilization of TR is not appropriate when the accuracy tolerance is so small. In all cases discussed below the time integrator is run for 2500 steps after the interrupt. Figure 16 shows the evolution of the mean vorticity and the acceleration using this refined strategy for each value of the viscosity parameter, when no perturbation is done. In the super-critical case of $\nu=1/200$ (bottom), there is a fast breakdown to the vortex-shedding solution. (Note that the evolution is plotted against physical time in this figure.) For both the subcritical ($\nu=1/175$) and near-critical ($\nu=1/185.6$) cases, there are long delays (until $t\approx 6e7$ and $t\approx 1.1e4$, respectively) after the interrupt, after which numerical instability kicks in and (as in the preceding section) generates a numerically noisy solution. The onset of instability is dramatically later for the subcritical case.

We explore the breakdowns in more depth in Figure 17, which shows magnified images of the noisy solution measures at the time they become unsteady. These images show that the magnitudes of the numerical oscillations (of order 1e-8 in the sub-critical case and 1e-7 in the near-critical case) are comparable to the time-stepping accuracy. Even when no explicit perturbation is done, time accuracy plays a role in long-term simulation to compute steady solutions in near-critical regimes.

Continuing this exploration of subcritical cases, we now consider the effects of perturbations at the interrupt. As in the previous section, we look at one perturbation that respects the reflectional symmetry of the flow solution in Figure 5 and is expected to be “benign” and one that breaks the reflectional symmetry and so is expected to be “lively”. The results for $\nu=1/175$ are shown in Figure 18 and those for $\nu=1/185.6$ are in Figure 19. In these figures, the vertical scaling for the mean vorticities are now set to be equal in order to discern differences for the two viscosity values. These images should be compared with those corresponding to analogous experiments with no perturbation in Figures 16–17. In particular, for the sub-critical viscosity $\nu=1/175$ with either type of perturbation, after a long delay, the solution moves away from a steady state. This is not surprising, since the same phenomenon occurs when no perturbation is done. The onset of periodic behavior for the perturbed data is slightly earlier than for no perturbation (and earlier still for the lively perturbation), but the magnitude of the oscillations is small. The results for $\nu=1/185.6$ (Figure 19) bear some similarity to these — most notably, the behavior for the benign perturbation is virtually identical to that for no perturbation (middle of Figure 16). The structure of the oscillations for the near-critical viscosity is more like that for the super-critical viscosity (for both perturbations as well as without perturbation, compare the images for lively perturbation in Figure 19 with the images in Figures 16–17). In contrast, for the sub-critical viscosity $\nu=1/175$, the structure of the oscillations is is more like that arising when no perturbation is done. But for the lively perturbation and $\nu=185.6$, the onset of unstable behavior is significantly earlier (bottom of Figure 19).

Finally, when we check to see what happens when the perturbation is significantly larger (of the order of the perturbation made in computing the pseudo-eigenvalues in Figure 6) we observe that there is a big difference in the time-stepping behavior in any case where the perturbation is not benign. This is illustrated by the results shown in the bottom plot in Figure 20. In this case the size of the lively perturbation is big enough to destabilize the integrator and a noisy periodic solution is computed. This mirrors the vortex-shedding solution that is computed in the unstable case but has an amplitude that is too small to be seen when plotted.

5 Concluding remarks

Our aims in this study were twofold. First, we have developed a new approach to assess the stability of dynamical systems by constructing perturbed systems based on collocation methods. This is reminiscent of methods for computing pseudospectra, but it has the advantage that the process of sampling (approximate) spectra is significantly less costly. Second, we compared the results of such assessments with the performance of time-stepping computations for a nontrivial application, the incompressible Navier-Stokes equations. In particular, for two benchmark problems, we examined the behavior of a stable integration scheme for simulating transient behavior for values of the viscosity in the system that are “sub-critical”, nearly critical (very slightly smaller than the critical value), and super-critical.

In general, we found that the predictions of instability made by the collocation method were consistent with the behavior of integrators: in the nearly critical regime (of parameter values, viscosity in this case), there is more sensitivity to perturbation than in the sub-critical regime, and outcomes are qualitatively like those for super-critical parameters. We also note that making such assessments is complicated somewhat by the delicate nature of computations in regimes at or near stability limits. Eigenvalues and pseudoeigenvalues are not the sole determining factor affecting stability; the form of the perturbation also plays a significant role.

**Acknowledgements: We thank Mark Embree and an anonymous referee for very constructive comments. **