Linear stability analysis#

Linear stability analysis is a method used to determine whether small disturbances to a steady base state grow or decay over time. The system is linearized around the base flow, and solutions are sought in the form of exponentially growing or decaying modes. By solving an eigenvalue problem, one identifies the growth rates and shapes of these modes. If any mode grows over time, the base flow is considered unstable. This approach provides insight into the natural tendencies of the system to amplify disturbances without external forcing, focusing purely on the system’s internal dynamics.

References:

Definition of the linear operator#

We start with a general nonlinear equation, written in compact form as

\[ \mathcal{B}\frac{\mathrm{d}\mathbf{q}}{\mathrm{d}t}=\mathcal{N}(\mathbf{q}), \]

where \(\mathbf{q} = [\mathbf{u}, p, \rho, ...]^{T} \in \mathbb{R}^N \) is the state vector (conservative variables), \(\mathcal{N}\) denotes the nonlinear operator of the system (e.g. Navier-Stokes), \(\mathcal{B}\) is a restriction operator; zero if the time derivative is not considered, else one.

We decompose the flow field into a time-invariant base flow, \(\overline{\mathbf{q}} \in \mathbb{R}^N\), representing either a time-averaged state or a fixed-point solution, and a perturbation, \(\mathbf{q}' \in \mathbb{R}^N\), such that

\[ \mathbf{q}(\mathbf{x},t) = \overline{\mathbf{q}}(\mathbf{x})+\epsilon\mathbf{q}'(\mathbf{x},t)\ . \]

In contrast to the resolvent analysis, we assume the perturbation to be small, \(\epsilon \ll 1\). Upon inserting this ansatz into the governing equations and linearisation, we arrive at the homogeneous equation for a linear time-invariant dynamical system describing the perturbation:

\[ \frac{\mathrm{d}\mathbf{q}'}{\mathrm{d}t}=\mathcal{L}(\mathbf{q}'), \]

where

\[ \mathcal{L}\equiv\nabla_{\mathbf{q}}\mathcal{N}\big|_{\overline{\mathbf{q}}}\in \mathbb{R}^{N\times N} \]

represents the Jacobian evaluated at the base state \(\overline{\mathbf{q}}\).

Eigenvalue analysis of the linear operator#

We assume the perturbation takes the form of normal modes,

\[ \mathbf{q}'=\hat{\mathbf{q}}\mathrm{e^{-\mathrm{j}\omega t}}+c.c., \]
where \(\hat{\mathbf{q}}\) denotes the complex amplitude, and the complex frequency is given by \(\omega = \omega_r + \mathrm{j}\omega_i\). Inserting this in the linearized perturbation equation leads to the generalized eigenvalue problem,

\[ \mathbf{A} \hat{\mathbf{q}} = -\mathrm{j} \omega \mathbf{B} \hat{\mathbf{q}}, \]

where \(\omega\) is the eigenvalue, \(\hat{\mathbf{q}}\) is the eigenvector and \(\mathbf{A}, \mathbf{B}\) are the discrete versions of the Jacobian \(\mathcal{L}\) and \(\mathcal{B}\), respectively.

Solving the eigenvalue problem yields complex eigenvalues, whose components have the following interpretations:

  • \(\omega_i > 0\): exponential temporal growth

  • \(\omega_i < 0\): exponential temporal decay

  • \(\omega_r\): oscillation frequency of the mode

  • \(\hat{\mathbf{q}}\): spatial structure (mode shape)

Additionally to the direct generalized eigenvalue problem above, one can formulate its adjoint problem with

\[ \mathbf{A}^\dagger\,\hat{\mathbf{q}}^{\dagger} = -\,\mathrm{j}\,\omega^\dagger\,\mathbf{B}\,\hat{\mathbf{q}}^{\dagger}, \]

Here, the dagger superscript denotes the Hermitian transpose for matrices and the complex conjugate for scalar quantities. The spectrum of the generalized adjoint eigenvalue problem is the complex-conjugate transpose of the direct spectrum, so that for each direct eigenmode, \(\hat{\mathbf{q}}_i\), there exists a corresponding adjoint eigenmode, \(\hat{\mathbf{q}}^\dagger_i\). The adjoint mode provides the spatial distribution of the flow’s receptivity, indicating where the corresponding direct mode is most sensitive to external forcing in a linear sense.

FELiCS implementation#

The discretized eigenvalue problem is solved using a SLEPc-based solver implemented in FELiCS. The solver computes a set of eigenvalues \(\omega_k\) (\(k=1,2,...,n\)) located nearest to a defined guess value \(\omega_{\text{guess}}\). For each eigenvalue \(\omega_k\), the corresponding eigenvector \(\hat{\mathbf{q}}_k\) is obtained.