2.1. Taylor coefficients for the base flow

We now show how the Taylor coefficients, i.e. the higher-order sensitivities of the base flow with respect to $a$, can be computed. In this subsection, the scalar parameter $a$ will be kept generic. We start with a detailed discussion on the first- and second-order coefficients, before presenting a general formulation for coefficients of arbitrary order.

The first-order sensitivity of $\boldsymbol {N}$ towards $a$ is

(2.7)\begin{align} \boldsymbol{N}_{1}\equiv\left.\frac{\text{d}\boldsymbol{N}(\boldsymbol{q}(a),a)}{\text{d}a}\right|_{a=a_{0}}=\left.\frac{\partial\boldsymbol{N}(\boldsymbol{q}(a),a_{0})}{\partial\boldsymbol{q}}\right|_{\boldsymbol{q}(a)=\boldsymbol{q}_{0}}\left.\frac{\text{d}\boldsymbol{q}(a)}{\text{d}a}\right|_{a=a_{0}}+\left.\frac{\partial\boldsymbol{N}(\boldsymbol{q}_{0},a)}{\partial a}\right|_{a=a_{0}}, \end{align}

where $\boldsymbol {q}_{0}$ is the base flow that solves (2.3) when $a=a_{0}$. For reasons of brevity, we replace the partial derivatives with subscripts from now on, and rewrite (2.7) as

(2.8)\begin{equation} \boldsymbol{N}_{1}\equiv\boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{1}+\boldsymbol{N}_{a}=\boldsymbol 0. \end{equation}

Here, ${\boldsymbol {N}_{\boldsymbol {q}}}$ is the linearization of the NS equations with respect to $\boldsymbol {q}$ at $\boldsymbol {q}_{0}$, thus a linear operator or a tensor of order two, which explicitly reads

(2.9)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{k}=\left(\begin{array}{@{}c@{}} (\boldsymbol{u}_{0}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_{k}+(\boldsymbol{u}_{k}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_{0}+\boldsymbol{\nabla} p_{k}-\dfrac{1}{{{Re}}}\,\nabla^{2}\boldsymbol{u}_{k}\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{k} \end{array}\right) \end{equation}

when applied on any $\boldsymbol {q}_{k}$. The explicit expression for $\boldsymbol {N}_a$ depends on the parameter chosen for the expansion, as discussed later.

The first-order sensitivity of $\boldsymbol {q}$ towards $a$ can now be obtained by solving the linear system

(2.10)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{1}=-\boldsymbol{N}_{a}. \end{equation}

The first-order Taylor polynomial then reads

(2.11)\begin{equation} \boldsymbol{q}(a)\approx\boldsymbol{q}^{(1)}(a_{0}+\epsilon)=\boldsymbol{q}_{0}+\boldsymbol{q}_{1} \epsilon. \end{equation}

Following an analogous procedure, the equation for the second-order expansion reads

(2.12)\begin{equation} \boldsymbol{N}_{2}=\boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{2}+\boldsymbol{q}_{1}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{1}+2\boldsymbol{N}_{\boldsymbol{q}a}\boldsymbol{q}_{1}+\boldsymbol{N}_{aa}=\boldsymbol 0,\end{equation}

where $\boldsymbol {N}_{\boldsymbol {qq}}$ is a symmetric bilinear form or a tensor of order three that reads

(2.13)\begin{equation} \boldsymbol{q}_{j}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{k}=\left(\begin{array}{@{}c@{}} (\boldsymbol{u}_{j}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_{k}+(\boldsymbol{u}_{k}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}_{j}\\[3pt] 0 \end{array}\right) \end{equation}

when applied on any $\boldsymbol {q}_{j}$ and $\boldsymbol {q}_{k}$. We will provide explicit expressions for $\boldsymbol {N}_{\boldsymbol {q}a}$ and $\boldsymbol {N}_{aa}$ later when we specify $a$. We obtain the second-order sensitivity of $\boldsymbol {q}$ with respect to $a$ by solving the linear equation system

(2.14)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{2}=-\boldsymbol{q}_{1}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{1}-2\boldsymbol{N}_{\boldsymbol{q}a}\boldsymbol{q}_{1}-\boldsymbol{N}_{aa}, \end{equation}

and the Taylor polynomial of order 2 is complete with

(2.15)\begin{equation} \boldsymbol{q}(a)\approx\boldsymbol{q}^{(2)}(a_{0}+\epsilon)=\boldsymbol{q}_{0}+\boldsymbol{q}_{1} \epsilon+\tfrac{1}{2}\,\boldsymbol{q}_{2} \epsilon^{2}. \end{equation}

This procedure can be pursued until arbitrary order. To obtain a general expression, it comes in handy that the higher-order derivatives of the NS equations with respect to the base flow vanish, i.e. that we have

(2.16)\begin{equation} \boldsymbol{N}_{\boldsymbol{qqq}}=\boldsymbol{0}=\boldsymbol{N}_{\boldsymbol{qqq}a}=\boldsymbol{N}_{\boldsymbol{qqqq}}=\cdots. \end{equation}

A closed formula for the computation of $\boldsymbol {q}_{k}$ for any scalar parameter $a$ then reads

(2.17) \begin{align} \boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{k} & =-\sum_{i=1}^{k-2} \sum_{j=1}^{k-i-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\left(\begin{array}{@{}c@{}} k-i\\ j \end{array}\right)\frac{1}{2}\,\boldsymbol{q}_{j}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}a^{i}}\boldsymbol{q}_{k-i-j}\ \nonumber\\ & \quad-\sum_{i=1}^{k-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\left[\frac{1}{2}\,\boldsymbol{q}_{i}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{k-i}+\boldsymbol{N}_{\boldsymbol{q}a^{i}}\boldsymbol{q}_{k-i}\right]-\boldsymbol{N}_{a^{k}}, \end{align}

where we adopted the notation $\boldsymbol {N}_{aa}=\boldsymbol {N}_{a^{2}}$ and so on for higher-order derivatives. The individual $\boldsymbol {q}_{k}$ have to be calculated one after the other from $k=1$ until $k=n$, the order of the desired Taylor polynomial. Note that the left-hand side of each linear equation system always contains the same linear operator $\boldsymbol {N}_{\boldsymbol {q}}$, which is advantageous for a fast computation of several higher orders.

2.2. Expansion of the equations with respect to the parameter $1/{{Re}}$

We will now specify the parameter $a$ as the inverse of the Reynolds number

(2.18)\begin{equation} a=\frac{1}{{{Re}}}, \end{equation}

and define explicit expressions for the operators defined in § 2.1. This choice for the parameter $a$ is convenient because $1/{{Re}}$ enters linearly in the NS equation. The case $a={{Re}}$ is also briefly discussed in this study and can be found in Appendix A. The first partial derivative of the NS equations (2.1) with respect to $a$ is

(2.19)\begin{equation} \boldsymbol{N}_{a}=\left(\begin{array}{@{}c@{}} -\nabla^{2}\boldsymbol{u}\\[2pt] 0 \end{array}\right), \end{equation}

and the higher order derivatives follow with

(2.20)\begin{equation} \boldsymbol{N}_{aa}=\boldsymbol{N}_{aaa}=\cdots=\boldsymbol{0}. \end{equation}

Similarly, the first partial derivative of the linear operator (2.9) with respect to $a$ is

(2.21)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}a}\boldsymbol{q}_{k}=\left(\begin{array}{@{}c@{}} -\nabla^{2}\boldsymbol{u}_{k}\\[2pt] 0 \end{array}\right) \end{equation}

and for the higher orders,

(2.22)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}aa}=\boldsymbol{N}_{\boldsymbol{q}aaa}=\cdots=\boldsymbol{0}. \end{equation}

Furthermore, it holds that

(2.23)\begin{equation} \boldsymbol{N}_{\boldsymbol{qq}a}=\boldsymbol{N}_{\boldsymbol{qq}aa}=\cdots=\boldsymbol{0}. \end{equation}

For the case $a=1/{{Re}}$, (2.17) thus simplifies to

(2.24)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{k}=\left\{ \begin{array}{@{}ll} -\boldsymbol{N}_{a}, & \text{for}\ k=1,\\ -\left(\begin{array}{@{}c@{}} k\\ 1 \end{array}\right)\boldsymbol{N}_{\boldsymbol{q}a}\boldsymbol{q}_{k-1} -\dfrac{1}{2}\displaystyle\sum_{i=1}^{k-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\boldsymbol{q}_{i}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}} \boldsymbol{q}_{k-i}, & \text{for}\ k>1 . \end{array}\right. \end{equation}

3.1. Taylor coefficients for the linear operator

In the following, we will express the higher-order sensitivities of $\boldsymbol{\mathsf{L}}$ through partial derivatives of $\boldsymbol {N}$ that were defined in § 2. Because we have defined $\boldsymbol{\mathsf{L}}\equiv \boldsymbol {N}_{\boldsymbol {q}}$, it follows that $\boldsymbol{\mathsf{L}}_{q}=\boldsymbol {N}_{\boldsymbol {qq}}$, $\boldsymbol{\mathsf{L}}_{a}=\boldsymbol {N}_{\boldsymbol {q}a}$, and so on. Thus it can be found that at first order,

(3.10)\begin{equation} \boldsymbol{\mathsf{L}}_{1}\boldsymbol{q}_{j}=\boldsymbol{q}_{1}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{j}+\boldsymbol{N}_{\boldsymbol{q}a}\boldsymbol{q}_{j},\end{equation}

and for second order,

(3.11)\begin{equation} \boldsymbol{\mathsf{L}}_{2}\boldsymbol{q}_{j}=\boldsymbol{q}_{2}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{j}+2 \boldsymbol{q}_{1}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}a}\boldsymbol{q}_{j}+\boldsymbol{N}_{\boldsymbol{q}aa}\boldsymbol{q}_{j} .\end{equation}

The arbitrary-order expression takes the form

(3.12)\begin{equation} \boldsymbol{\mathsf{L}}_{k}\boldsymbol{q}_{j}=\boldsymbol{q}_{k}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{j}+\sum_{i=1}^{k-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\boldsymbol{q}_{k-i}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}a^{i}}\boldsymbol{q}_{j}+\boldsymbol{N}_{qa^{k}}\boldsymbol{q}_{j}.\end{equation}

3.2. Taylor coefficients for the eigenvalue

The sensitivities of $\lambda (a)$ are derived by calculating the sensitivities towards $a$ of (3.1). The first-order eigenvalue sensitivity towards $a$, $\lambda _{1}$, is given by the sesquilinear scalar product

(3.13)\begin{equation} \lambda_{1}=\frac{1}{c}\left\langle\hat{\boldsymbol{q}}_{0}^{*},\boldsymbol{\mathsf{L}}_{1}\hat{\boldsymbol{q}}_{0}\right\rangle,\quad c\equiv\langle\hat{\boldsymbol{q}}_{0}^{*} ,\boldsymbol{\mathsf{M}}\hat{\boldsymbol{q}}_{0}\rangle .\end{equation}

This expression can be derived by calculating the first-order sensitivity towards $a$ at $a_{0}$ of the eigenproblem (3.1),

(3.14)\begin{equation} \left(\boldsymbol{\mathsf{L}}_{0}-\lambda_{0}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{1}+\left(\boldsymbol{\mathsf{L}}_{1}-\lambda_{1}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{0}=\boldsymbol{0},\end{equation}

and then taking the scalar product with the corresponding zeroth-order adjoint eigenvector; see e.g. Luchini & Bottaro (Reference Luchini and Bottaro2014b). From (3.3), it holds that

(3.15)\begin{equation} \left\langle\hat{\boldsymbol{q}}_{0}^{*},\left(\boldsymbol{\mathsf{L}}_{0}-\lambda_{0}\boldsymbol{\mathsf{L}}\right)\hat{\boldsymbol{q}}_{k}\right\rangle= 0 \end{equation}

for any $k$, therefore $\lambda _{1}$ can be calculated without knowing $\hat {\boldsymbol {q}}_{1}$.

The second-order sensitivity is derived in the exact same manner, i.e. by calculating the second-order sensitivity of (3.1) and then taking the scalar product with the corresponding zeroth-order adjoint eigenvector (Orchini et al. Reference Orchini, Magri, Silva, Mensah and Moeck2020; Boujo Reference Boujo2021):

(3.16)\begin{equation} \lambda_{2}=\frac{1}{c}\left\langle\hat{\boldsymbol{q}}_{0}^{*} ,\boldsymbol{\mathsf{L}}_{2}\hat{\boldsymbol{q}}_{0}\right\rangle+\frac{2}{c}\left\langle\hat{\boldsymbol{q}}_{0}^{*} ,\left(\boldsymbol{\mathsf{L}}_{1}-\lambda_{1}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{1}\right\rangle .\end{equation}

Here, the first-order coefficient of the eigenvector, $\hat {\boldsymbol {q}}_{1}$, is needed, which will be derived in the next subsection. This procedure can be followed until arbitrary order, arriving at

(3.17)\begin{equation} \lambda_{k}=\frac{1}{c}\left\langle\hat{\boldsymbol{q}}_{0}^{*} ,\boldsymbol{\mathsf{L}}_{k}\hat{\boldsymbol{q}}_{0}\right\rangle+\frac{1}{c}\sum_{i=1}^{k-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\left\langle\hat{\boldsymbol{q}}_{0}^{*} ,\left(\boldsymbol{\mathsf{L}}_{i}-\lambda_{i}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{(k-i)}\right\rangle . \end{equation}

Note that the respective lower-order sensitivities of the eigenvalue and the eigenvector have to be known beforehand, which means that higher-order eigenvalue coefficients have to be calculated one after the other.

3.3. Taylor coefficients for the eigenvector

For the first-order sensitivity $\hat {\boldsymbol {q}}_{1}$, the linear equation system (3.14) is solved:

(3.18)\begin{equation} \left(\boldsymbol{\mathsf{L}}_{0}-\lambda_{0}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{1}=-\left(\boldsymbol{\mathsf{L}}_{1}-\lambda_{1}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{0} .\end{equation}

The first-order coefficient of the eigenvalue, $\lambda _{1}$, has to be calculated first (see (3.13)). The operator $(\boldsymbol{\mathsf{L}}_{0}-\lambda _{0}\boldsymbol{\mathsf{M}})$ is singular by definition, but the equation system is solvable because of the Fredholm alternative (Oden Reference Oden1979): if the right-hand side of the equation system is orthogonal to the adjoint eigenvector $\hat {\boldsymbol {q}}_{0}^{*}$, then a solution can be found. This condition is met because $\lambda _{1}$ in (3.13) was calculated by imposing orthogonality of the right-hand side to $\hat {\boldsymbol {q}}_{0}^{*}$.

Note that the solution $\hat {\boldsymbol {q}}_{1}$ is not defined uniquely, but consists of a component orthogonal to $\hat {\boldsymbol {q}}_{0}$, which is fully determined, and a component parallel to $\hat {\boldsymbol {q}}_{0}$, which is undetermined. As shown e.g. in Mensah et al. (Reference Mensah, Orchini and Moeck2020) and Boujo (Reference Boujo2021), the latter has no influence on the value of $\lambda _{2}$, defined in (3.16). More generally, it can be shown that the component parallel to $\hat {\boldsymbol {q}}_{0}$ has no influence on the value of any higher-order eigenvalue coefficient. Once a choice has been made on the definition of $\hat {\boldsymbol {q}}_{1}$, the same choice has to be maintained for the calculation of the higher-order coefficients, since it has an influence on the definition of the higher-order eigenvector sensitivities.

The second-order sensitivity $\hat {\boldsymbol {q}}_{2}$ is obtained by solving the equation system

(3.19)\begin{equation} \left(\boldsymbol{\mathsf{L}}_{0}-\lambda_{0}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{2}=-\left(\boldsymbol{\mathsf{L}}_{2}-\lambda_{2}\right)\hat{\boldsymbol{q}}_{0}-2\left(\boldsymbol{\mathsf{L}}_{1}-\lambda_{1}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{1} .\end{equation}

As for the first order, despite the operator on the left-hand side being singular, the equation system is guaranteed to be solvable thanks to the Fredholm alternative condition imposed by (3.16). Analogously to the first order, also the component of $\hat {\boldsymbol {q}}_{2}$ parallel to $\hat {\boldsymbol {q}}_{0}$ is not defined uniquely. It was shown in Mensah et al. (Reference Mensah, Orchini and Moeck2020), and generalized to arbitrary high order, that any solution of (3.19) will lead to the same value of $\lambda _{3}$ or any subsequent eigenvalue coefficient. Similarly to the first-order case, once defined, it is important that $\hat {\boldsymbol {q}}_{2}$ remains the same for all higher-order calculations.

A general expression for the equation system of an arbitrary $\hat {\boldsymbol {q}}_{k}$ is

(3.20)\begin{equation} \left(\boldsymbol{\mathsf{L}}_{0}-\lambda_{0}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{k}=-\left(\boldsymbol{\mathsf{L}}_{k}-\lambda_{k}\right)\hat{\boldsymbol{q}}_{0}-\sum_{i=1}^{k-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\left(\boldsymbol{\mathsf{L}}_{i}-\lambda_{i}\boldsymbol{\mathsf{M}}\right)\hat{\boldsymbol{q}}_{(k-i)}.\end{equation}

The considerations concerning solvability and uniqueness of the Taylor coefficients can be extended to arbitrary order.

5.1. Taylor expansion at ${{Re}}=47$

We first choose the expansion point to be Reynolds number 47, so that $a_0=\tfrac {1}{47}$, at the onset of periodic vortex shedding, and predict with our method the base flow and the leading eigenvalue of its eigenproblem.

For ${{Re}}\le 47$, solving the NS base flow equations (2.1) is particularly easy, because the linear operator has only eigenvalues with negative real part, meaning that several iterative methods can be applied when solving the nonlinear equation. Here, the base flow equations were first solved at ${{Re}}_{0}=47$ with a Newton method. Then the coefficients for the Taylor polynomials were calculated until order 40 for the base flow, (2.17), as well as for the eigenvalues and eigenvectors, (3.17) and (3.20).

In figure 1, the predicted profiles of the velocity in the $x$-direction, $u$, and the velocity in the $y$-direction, $v$, are plotted at $\boldsymbol {x}=(2,0)$. The lightest (darkest) curve is the prediction of the order 1 (order 40) Taylor polynomial. The exact base flow, calculated with the method described in § 5.2, is depicted with black circles. It can be seen that for ${{Re}}=100$, the prediction is highly accurate, already at order $n\approx 7$ no difference to the true base flow is visible. The same goes for ${{Re}}=150$: at order $n\approx 12$, the difference with respect to the true base flow is negligible. At ${{Re}}=200$, however, the Taylor series diverges. The divergence can be seen on the profile of $v$. The prediction becomes worse at increasing order of the polynomials. The same behaviour can be observed in the whole computational domain. Further details on convergence will be discussed in § 5.3. Between ${{Re}}=150$ and ${{Re}}\approx 180$, the base flow is predicted less and less accurately, but the prediction improves with each order, until it reaches its convergence limit at ${{Re}}\approx 180$, where the prediction becomes worse when the order increases. Note that enforcing the transversal velocity $v$ to be zero at the symmetry axis does not change the convergence behaviour.

Figure 1.

Taylor series expansion of the base flow along $a=1/{{Re}}$ conducted at $a_{0}=1/47$. Velocity profiles of the predicted base flows at $\boldsymbol {x}=(2,y)$ are plotted for Reynolds numbers 100 (a,b), 150 (c,d) and 200 (e,f). Directly calculated velocity profiles are plotted as black circles for comparison.

When expanding the leading eigenvalue into a Taylor polynomial, a similar convergence behaviour can be observed. In figure 2, the real and imaginary parts of the eigenvalue are plotted over the Reynolds number. Again, the lightest curve is the Taylor polynomial of order 1, the darkest of order 40. Directly calculated eigenvalues are plotted as black circles for comparison. Until ${{Re}}=150$, the eigenvalue can be predicted with the Taylor expansion. From then on, the prediction becomes less accurate, until the Taylor polynomials start to diverge. A difference between the convergence quality of the eigenvalue and the base flow cannot be observed. This suggests that at least in this case, the convergence behaviour of the eigenvalue is dictated by the convergence of the base flow. This is possible because the base flow coefficients $\boldsymbol {q}_{k}$ feed into the eigenvalue corrections (3.17).

Figure 2.

Taylor series expansion of the eigenproblem along $a=1/{{Re}}$ conducted at $a_{0}=1/47$. The predicted real and imaginary parts of the eigenvalue are plotted over the Reynolds number for orders $n$ from 1 to 40. Directly calculated eigenvalues are plotted as black circles for comparison.

In figure 3, the computation time is plotted for the calculations of the various higher-order sensitivities. The calculation time is measured on one CPU and scaled by the time that is needed to solve one linear equation system with only real entries. This scaled CPU time is plotted on a logarithmic scale over the order of the Taylor polynomial. For the base flow sensitivities, computation time increases approximately exponentially. However, the calculations of base flow sensitivities up to order 40 is still faster than solving the nonlinear base flow equations (dashed line, calculated with a Newton solver in 7 iterations). For the eigenvalue sensitivities, a similar behaviour can be observed. These computation times are larger than those of the base flow sensitivities, because for each coefficient, a linear equation system with complex entries has to be solved, which has twice as many degrees of freedom. Again, calculating all coefficients until order 40 is faster than solving the eigenproblem directly for one eigenvalue and eigenvector (calculated as a shift-inverse problem with the Python package SciPy). For both eigenvalue and base flow coefficients, the computation time for approximately order 16 is twice as long as the time for the first solution of the linear equation system, i.e. the time to calculate the first-order coefficients. Thus, once the first-order coefficients are computed, higher orders follow at comparatively low computational costs.

Figure 3.

Computation time for the higher-order sensitivities of the base flow and the eigenvalue, plotted over the order $n$. Dashed lines refer to the computation time for solving the nonlinear base flow equations (red) and the eigenproblem (blue). Note the logarithmic scale.

5.2. Subsequent computation of the base flow and its leading eigenvalue up to ${{Re}}=1000$

Even though the convergence of the base flow polynomial is finite, we can extend our prediction range by starting a new Taylor expansion at a new $a_{0}$ inside the convergence radius. Here, the following procedure was applied. Starting at ${{Re}}_{0}=30$, so that the expansion point is $a_0=\tfrac {1}{30}$, the base flow was expanded into a Taylor series until Taylor order 40. Then the Taylor approximation was evaluated at ${{Re}}_{0}=40$, thus yielding the base flow solution at a new Reynolds number inside the converged region. From there, a new Taylor series was expanded around $a_0=\tfrac {1}{40}$, and so on. At each step, the Reynolds number was increased by 10. By this method, the base flow was computed successfully until ${{Re}}=1000$. We note that the small step size used here and the high expansion order chosen may not have been strictly necessary to obtain an accurate solution of the NS equations as per (2.1), and to keep its residual (numerically) zero. However, a detailed investigation on the optimal performance of our algorithm in terms of step sizes and perturbation orders is beyond the scope of our work and will not be discussed here.

In figure 4, the calculated base flow is shown for Reynolds numbers 100, 500 and 1000. It can be seen that the gradients get sharper for higher Reynolds numbers. Accordingly, the width of the recirculation zone in the wake increases, while the length varies non-monotonically as the Reynolds number goes up. Near the symmetry axis, the back flow velocity increases significantly.

Figure 4.

Base flow velocity components (a,c,e) $u(\boldsymbol {x})$ and (b,d,f) $v(\boldsymbol {x})$, for (a,b) ${{Re}}=100$, (c,d) ${{Re}}=500$, and (e,f) ${{Re}}=1000$. Computed with subsequent Taylor polynomial prediction, starting point was ${{Re}}=30$.

Next, the leading eigenvalue was traced until Reynolds number 1000. To do so, the eigenproblem was solved for several base flow solutions that were obtained with the method described previously, and the eigenvalue predictions conducted with the perturbation method were overlaid to construct a continuous curve. In figure 5(a), the real and imaginary parts of the eigenvalue versus ${{Re}}$ are shown. The dotted vertical lines indicate where the eigenproblem was solved. In the area between Reynolds numbers 200 and 420, the Reynolds number interval had to be decreased due to the poor convergence of the Taylor polynomials. This coincides with the area in which the convergence of the base flow is poor, thus indicating that the convergence of the eigenvalue is determined by the base flow in the whole Reynolds number range that was considered. Figure 5(b) shows part of the eigenvalue spectrum at ${{Re}}=1000$ (black dots) and the evolution of the leading eigenvalue with increasing Reynolds number. While the imaginary part decreases almost monotonically over the Reynolds number range, the real part first increases very quickly, and then increases and decreases again with a slowing rate. At ${{Re}}=1000$ it is still the leading eigenvalue of the entire spectrum. Figure 6 displays the eigenvectors of the traced eigenvalue for Reynolds numbers 100, 500 and 1000. Similar to the base flow itself, the gradients get sharper when the Reynolds number increases.

Figure 5.

(a) Real and imaginary part of the traced eigenvalue over the Reynolds number. (b) Detail of the eigenvalue spectrum of the base flow at ${{Re}}=1000$. The traced eigenvalue is shown in varying colours for different Reynolds numbers.

Figure 6.

Eigenvectors for (a,b) ${{Re}}=100$, (c,d) ${{Re}}=500$, and (e,f) ${{Re}}=1000$.

5.3. Remarks on convergence of the Taylor polynomials

In this subsection, we will take a closer look at the convergence behaviour of the Taylor polynomials that were computed in the previous subsection.

To quantify the convergence of the predicted base flow, we look at the $L_{2}$-norm of the NS base flow equations. We define

(5.1)\begin{equation} L_{2}(n,a,a_{0}):=\frac{1}{\sqrt{A}}\left\Vert \boldsymbol{N}\left(\boldsymbol{q}^{(n)}(a), a\right)\right\Vert _{2}=\frac{1}{\sqrt{A}}\left\langle \boldsymbol{N}\left(\boldsymbol{q}^{(n)}(a), a\right)\!,\boldsymbol{N}\left(\boldsymbol{q}^{(n)}(a), a\right)\right\rangle ^{{1}/{2}}\end{equation}

with the surface area of the computational domain $A$, and equations (2.1), (2.3) and (2.4). It holds that

(5.2)\begin{equation} \lim_{n\rightarrow \infty} L_{2}(n,a,a_{0}) = 0 \end{equation}

if

(5.3)\begin{equation} |\epsilon|\equiv|a-a_{0}|< r \end{equation}

for a certain convergence radius $r$. In a looser sense, one may state that the series truncated at a finite order $n$ provides a good approximation to the actual solution if $L_{2}(n,a,a_{0}) < \delta$, where $\delta$ is a small error threshold. This means that the NS base flow equation should be (approximately) fulfilled for the predicted base flow in a certain region around $a_{0}$, and the smaller $L_{2}$, the better the prediction.

First, we look at the convergence of the base flow prediction for the parameter $a={{Re}}$. The parameter-specific details of the calculations can be found in Appendix A. Figure 7(a) shows the $L_{2}$-norm of the NS equations for the predicted base flows when expanded around $a_{0}={{Re}}_{0}=30$. The lightest curve depicts the first-order prediction, $n=1$, and the darkest curve depicts the prediction of order $40$. The norm is shown on a logarithmic scale.

Figure 7.

Taylor expansion of the base flow; the $L_{2}$-norm of the NS equations is plotted over the expansion parameter and serves as an error measurement. (a) Expansion parameter $a={{Re}}$, conducted at $a_{0}=30$, diverges for $a<0$ and $a>2a_{0}=60$ (vertical dashed line). (b) Expansion parameter $a=1/{{Re}}$, conducted at $a_{0}=1/10$, diverges for $a<0$ and $a>2a_{0}=1/5$. Note the logarithmic scale of the $L_{2}$-norm. Mesh A (AF).

It can be seen that the Taylor polynomial converges if $|a-a_{0}|< a_{0}$, and diverges otherwise. This behaviour can be explained by the manner in which the Reynolds number appears in the NS equations. Just like the Taylor series for the function $f(x)=1/x$, the NS equations have a singularity at ${{Re}}=0$. This singularity limits the convergence of any power series related to the equations, and the base flow convergence radius is therefore $a_{0}={{Re}}_{0}$. Therefore, the base flow cannot be predicted for ${{Re}}>2{{Re}}_{0}$ if the chosen parameter is the Reynolds number itself, because the related series diverges.

Using the parameter $a=1/{{Re}}$ for the perturbation expansion yields a qualitatively similar but quantitatively different convergence result. In figure 7(b), the $L_{2}$-norm of the NS equations is depicted for the predicted base flows, when expanded around $a_{0}=1/{{Re}}_{0}=1/10$. It is known that the NS equations have a singular point at $a=0$, i.e. ${{Re}}=\infty$. This divergence is due to the fact that with no viscous term, the NS equations turn into the Euler equations, which need a different boundary condition at the cylinder wall. Thus also in this case, $a_{0}$ defines an upper limit for the convergence radius. Nonetheless, there may exist other singularities that would reduce the radius of convergence further.

While the convergence behaviour when expanding around small Reynolds numbers is defined by the singular point at ${{Re}}=0$, the behaviour when starting at higher Reynolds numbers is not described as easily. In figure 8, the convergence behaviour of the base flow Taylor polynomials is depicted for the parameter $a=1/{{Re}}$ and $a_{0}=1/47$. The $L_{2}$-norm of the NS equations against the $a$ parameter is plotted in figure 8(a). For convenience, the same results are shown in figure 8(b) by using the parameter $1/a = {{Re}}$ on the horizontal axis. In figure 8(a), it can be seen that the Taylor polynomials converge in a symmetric region around the expansion point. This symmetry is lost in figure 8(b), because the scaling $1/a$ stretches the results in different ways for small/large Reynolds numbers. It can be seen that a good approximation is obtained up to ${{Re}}\approx 180$. We remark that if we had chosen to use $a={{Re}}$ as the expansion parameter, then an upper limit for the convergence radius would be given by ${{Re}}_0 = 47$, because of the singularity of the equations at ${{Re}}=0$. Thus the result would surely diverge for ${{Re}}>94$. This comparison shows that the choice of the expansion parameter has an influence on the largest achievable Reynolds number, and that in the case considered here, using the parameter $a=1/{{Re}}$ allows us to obtain accurate results for larger values of ${{Re}}$.

Figure 8.

Taylor series expansion of the base flow along $a=1/{{Re}}$; a norm of the NS equations is plotted over (a) the inverse of the Reynolds number, and (b) the Reynolds number itself. Starting point is $a_{0}=1/47$. Mesh A (AF).

6.1. Model equations

Following Hill (Reference Hill1992), Marquet et al. (Reference Marquet, Sipp and Jacquin2008) and Boujo (Reference Boujo2021), we model the small control cylinder of diameter $d$ located in $\boldsymbol {x}_c$ as a steady forcing acting on the base flow. Its magnitude is defined by the drag force that would act on the small cylinder in uniform flow, and the forcing is distributed along a two-dimensional Gaussian function. In particular, we use the same drag coefficient as Boujo (Reference Boujo2021), but add a model for the variance of the Gaussian function, to make the predictions more accurate and to avoid convergence issues of the Taylor polynomial. The forcing is thus defined as

(6.1)\begin{equation} \boldsymbol{F}(\boldsymbol{x}_c,d)=\tfrac{1}{2}\,d\,C_{d}(\boldsymbol{x}_{c})\,\Vert\boldsymbol{u}_{0}(\boldsymbol{x}_{c})\Vert\,\boldsymbol{u}_{0}(\boldsymbol{x}_{c})\, \boldsymbol{\chi}(\boldsymbol{x}-\boldsymbol{x}_{c},\sigma^{2}), \end{equation}

with the local Reynolds number

(6.2)\begin{equation} {{Re}}_{d}(\boldsymbol{x}_{c})\equiv\frac{d\,\Vert\boldsymbol{u}_{0}(\boldsymbol{x}_{c})\Vert}{\nu}, \end{equation}

and the drag coefficient $C_{d}(\boldsymbol {x}_{c})=0.8558+10.05({{Re}}_{d}(\boldsymbol {x}_{c}))^{-0.7004}$. Pivotal is the velocity of the unperturbed base flow $\boldsymbol {u}_0$ at the centre $\boldsymbol {x}_c$ of the small cylinder. Here, $\boldsymbol {\chi }(\boldsymbol {x}-\boldsymbol {x}_{c},\sigma ^{2})$ is a two-dimensional Gaussian function that is distributed symmetrically around $\boldsymbol {x}_{c}$ with variance $\sigma ^{2}$. We add a linear dependence of the standard deviation with regard to the magnitude of the forcing to improve results in comparison to the small control cylinder incorporated directly into the grid. The standard deviation is thus

(6.3)\begin{equation} \sigma = \max\left[\sigma_{min}, d\left(\tfrac{1}{2}\,C_d\,\Vert\boldsymbol{u}_0(\boldsymbol{x}_c)\Vert^2-1 \right)\right]\!, \end{equation}

where $\sigma _{{min}}$ ensures that the Gaussian function can be depicted accurately enough with the given grid resolution. Here, we choose $\sigma _{{min}}=0.025$ in combination with grid B (see table 1). We recall that only dimensionless variables appear in (6.3), and that the velocity and diameter of the small control cylinder has to be made dimensionless for (6.3) to hold.

In the following, we consider a small control cylinder of diameter $d=0.1$. The expansion parameter $a$ is chosen to be the amplitude of the steady forcing such that the base flow equations become

(6.4) \begin{equation} \boldsymbol{N}(\boldsymbol{q},a)=\left(\begin{array}{@{}c@{}} (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u}+\boldsymbol{\nabla} p-\dfrac{1}{Re}\,\nabla^{2}\boldsymbol{u}\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} \end{array}\right)+a\left(\begin{array}{@{}c@{}}\boldsymbol{F}\\0\end{array}\right)= \boldsymbol 0 . \end{equation}

Thus $a$ is a scalar factor that is expanded from $a_{0}=0$. It is related to the diameter of the small control cylinder in the above model via

(6.5)\begin{equation} a(d,\boldsymbol{x}_{c})\equiv\frac{d\,C_{d}(\boldsymbol{x}_{c})}{0.1\,C_{0.1}(\boldsymbol{x}_{c})} \end{equation}

such that $a(d=0.1,\boldsymbol {x}_{c})=1$.

The first-order partial derivative of the steady NS equations is

(6.6)\begin{equation} \boldsymbol{N}_{a}=\left(\begin{array}{@{}c@{}} \boldsymbol F\\0 \end{array}\right), \end{equation}

and (2.17) simplifies to

(6.7)\begin{equation} \boldsymbol{N}_{\boldsymbol{q}}\boldsymbol{q}_{k}=\left\{ \begin{array}{@{}ll} -\left(\begin{array}{@{}c@{}} \boldsymbol F\\0 \end{array}\right)\!, & \text{for}\ k=1,\\ -\dfrac{1}{2}\displaystyle\sum_{i=1}^{k-1}\left(\begin{array}{@{}c@{}} k\\ i \end{array}\right)\boldsymbol{q}_{i}^{\rm T}\boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{k-i}, & \text{for}\ k>1 . \end{array}\right. \end{equation}

The coefficients for the Taylor expansion of the linear operator, needed in the prediction of the eigenvalue, are (see (3.12))

(6.8)\begin{equation} \boldsymbol{\mathsf{L}}{}_{k}\boldsymbol{q}_{j}=\boldsymbol{q}_{k}^{\rm T} \boldsymbol{N}_{\boldsymbol{qq}}\boldsymbol{q}_{j}. \end{equation}

6.2. Higher-order sensitivity maps of the control cylinder

The higher-order sensitivity maps are calculated by computing the Taylor coefficients of the model described above up to order 10 at 0.1 increments in the $x$- and $y$-directions. Thus for each Reynolds number, a total of 496 control cylinder positions were computed. As discussed before, the preconditioning of the linear equation systems can be exploited for the whole map once the undisturbed base flow and eigenproblem are given, making the computation very efficient. Computing one map until Taylor order 10, with the implementation described here, cost the equivalent of approximately 30 base flow computations and 26 eigenvalue computations (on one CPU, non-parallel). We recall that with this method, the base flows and leading eigenvalues were predicted not only for 496 discrete cylinder positions, but also for a continuous range of control cylinder diameters within the validity of the used model. Also, as will be seen later, calculating the coefficients up to order 10 is not needed for sufficient convergence at most positions. Because the computation time increases exponentially with every order (see figure 3), taking into account the convergence behaviour would reduce the computation time even further.

In figure 9 the regions in which a small control cylinder stabilizes the flow are drawn for Reynolds numbers 60, 70, 80 and 90, and Taylor polynomial orders from 1 to 10. The drawn lines are B-spline interpolations of the calculated grid data. For comparison, the eigenproblem was also solved with the small control cylinder inserted directly into the grid, with finer grid resolution in the vicinity of the wall, and Dirichlet boundary conditions. To identify the value at which the growth rate of the eigenvalue is zero, for each point, two control cylinder positions were calculated at distance 0.06, and the resulting position was computed by first-order interpolation. When necessary, one or two additional positions were calculated. This procedure was possible because the approximate coordinates of the resulting points were already known from the high-order perturbation expansion results.

Figure 9.

Predicted stabilizing regions until order 10, using a steady forcing to model the small control cylinder with diameter 0.1, for Reynolds numbers 60, 70, 80 and 90. Black circles are values calculated from incorporating the small control cylinder directly into the grid. Small crosses are experimental data from Strykowski & Sreenivasan (Reference Strykowski and Sreenivasan1990). Note that the minimum visible order varies: for ${{Re}}=60$, order 1 already depicts a stabilizing region; for ${{Re}}=70$, order 2 is needed; for ${{Re}}=80$, order 3; and for ${{Re}}=90$, order 5 (the lighter the contours, the higher the expansion order).

As can be seen in figure 9, the higher-order Taylor predictions of the steady forcing model are generally in very good agreement with the values from a direct computation using the embedded control cylinder. The deviation for ${{Re}}=60$ and 70, visible at the streamwise end of the stabilizing region, stems from the grid resolution: with a finer grid, a smaller $\sigma _{min}$ in (6.3) would have been possible. This was verified by performing calculations on the finer BF grid; see table 1. Since in these regions the respective base flow shows a comparatively low absolute velocity and thus a comparatively low forcing value is needed, the borders of the stability regions would have been more accurate with a steeper Gaussian function. At higher Reynolds numbers, this problem does not occur.

Note that the higher the Reynolds number, the higher the expansion order needed to depict the stable regions accurately. While for ${{Re}}=60$ the second-order prediction is already fairly accurate, for ${{Re}}=70$ it is vital for detecting a stabilizing area at all, since the first-order expansion predicts no stabilization in the whole computational domain. These results compare well to the second-order results from Boujo (Reference Boujo2021). Subsequently, for ${{Re}}=80$, expansion order 3 is needed to detect the stabilizing area, and for ${{Re}}=90$, the region appears only for orders $n\ge 5$. Experimental data from Strykowski & Sreenivasan (Reference Strykowski and Sreenivasan1990) for Reynolds numbers 60 and 70 have been added for completeness. They are slightly more unstable than the numerical data. The authors themselves discuss that, due to the experimental set-up, three-dimensional effects occur, and thus the two-dimensionality assumption is not completely accurate. In this light, the results compare well to the numerical data depicted here.