<figure><img src="image_1.png" /></figure>

HAL
open science

# A numerical method for the computation of bifurcation points in fluids mechanics

Jean-Marc Cadou, Michel Potier-Ferry, Bruno Cochelin

► To cite this version:

Jean-Marc Cadou, Michel Potier-Ferry, Bruno Cochelin. A numerical method for the computation of bifurcation points in fluids mechanics. European Journal of Mechanics - B/Fluids, 2006, 25, pp.234-254. 10.1016/j.euromechflu.2005.07.002. hal-00089043

HAL Id: hal-00089043

https://hal.science/hal-00089043v1

Submitted on 9 Aug 2006

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
A numerical method for the computation of bifurcation points
in fluid mechanics

J.M. Cadou<sup>a</sup>, M. Potier-Ferry<sup>b</sup>, B. Cochelin<sup>c</sup>

<sup>a</sup> Laboratoire Génie Mécanique et Matériaux, Université de Bretagne Sud, Rue de Saint Maudé, BP 92116, 56321 Lorient cedex, France

<sup>b</sup> Laboratoire de Physique et Mécanique des Matériaux, ISGMP, Université de Metz, Ile du Saulcy, 57045, Metz, France

<sup>c</sup> Laboratoire de Mécanique et Acoustique, Ecole Supérieure de Mécanique de Marseille, ESM2, IMT, Technopôle de Château Gombert, 13451, Marseille, France

Two original algorithms are proposed for the computation of bifurcation points in fluid mechanics. These algorithms consist of finding the zero values of a specific indicator. To compute this indicator a perturbation method is used which leads to an analytical expression of this indicator. Two kinds of instability are considered: stationary and Hopf bifurcations. To prove the efficiency and advantages of such numerical methods several numerical tests are discussed.

Keywords: Bifurcation indicators; Hopf bifurcation; Asymptotic numerical method

# 1. Introduction

In this paper, we are interested in numerical investigations of two types of classical instability: stationary and Hopf bifurcations. A stationary bifurcation corresponds to a transition from a steady flow (usually symmetric) to another flow (in most cases, with nonsymmetric solutions), whereas a Hopf bifurcation indicates the appearance of a time-periodic solution from a steady branch. The conditions necessary for encountering a Hopf bifurcation are as follows: [1,2] a stationary solution exists for which two of the eigenvalues of the Jacobian matrix cross the imaginary axis. This means that if $\zeta_k$ are the eigenvalues of the Jacobian matrix, then at a Hopf bifurcation point, two eigenvalues are purely imaginary (i.e. $\zeta_1 = \bar{\zeta}_2 = i\omega$, Re($\zeta_k$) $\neq$ 0 for $k \ge 3$).

Common methods for computing Hopf bifurcation points precisely are generally divided into two families: the “indirect” and the “direct” method. The indirect method consists of finding the solutions of the equation Re($\zeta$) = 0 along the stationary solution branches. This method requires, for each steady state solution, the computation of the eigenvalues of the Jacobian matrix [3]. In some cases, when the number of unknowns of the system is too high, it is sometimes impossible or at least highly prohibitive to compute all eigenvalues. In such cases, only part of the spectrum is calculated [4] using, for example, an Arnoldi algorithm [5].
The direct method consists of solving an augmented system whose solutions are Hopf bifurcation points [6,7]. This method does not need to follow the stationary solution branch, but it requires good initial values (a good initial Reynolds number and a good approximation of the initial eigenvalue) to obtain an acceptable convergence.

In this paper, we propose another method, which can also be included in the family of indirect methods using a different function (not, Re(ζ) = 0) to characterize the bifurcation points. This method is based on the introduction of indicators with the distinctive feature of having zero values at singular points. These indicators were first developed in structural stability mechanics. Boutyour et al. [8] have defined a stationary indicator. Bensaadi [9] has proposed theoretical work on the determination of Hopf bifurcation and Tri et al. [10] have applied this theory in structural problems. We use these indicators and adapt them for fluid instability problems.

The basic idea is to perturb the stationary solution (written $U_λ$) by a load vector $μf$, where $μ$ is the intensity and $f$ is a given vector force. The linearized and perturbed problem arising from this perturbation can then be written in the following form:

$$ \mathcal{E}_T(U_\lambda, \omega) \cdot \Delta\mathcal{V} = \mu f \quad (1) $$

where $\Delta\mathcal{V}$ is the fluctuation in the velocity resulting from the perturbation force. $\mathcal{E}_T(U_\lambda, \omega)$ is an operator which depends on the Reynolds number (via the fundamental stationary solution $U_\lambda$) and, for the Hopf bifurcation, of the angular frequency $\omega$. This operator depends on whether it is a stationary or a Hopf bifurcation that we are looking for. In Eq. (1) one can see that when $μ$ is null, then (1) is equivalent to an eigenvalue problem. So the quantity $μ$ is our indicator of bifurcation (as will be shown later): computation of the bifurcation points consists of finding the point of the fundamental solution branch where $μ$ is equal to zero.

When stationary bifurcations are considered, the indicator only depends on the fundamental point, $U_λ$. This stationary solution is computed using a perturbation method combined with the finite element method (Asymptotic Numerical Method [11–14], “ANM”). Thus, $U_λ$ is an explicit analytical function of a path following parameter “a”. Considering the previous idea, the bifurcation indicator only depends on this path parameter. This property has been used by Boutyour et al. [8] and Tri et al. [15] to determine the stationary bifurcation points in structural and fluid mechanics. They used a perturbation method to compute both the indicator $μ$ and the vector $\Delta\mathcal{V}$. The advantages of such a method are as follows: On one hand, in order to determine the quantities $μ$ and $\Delta\mathcal{V}$ the operator $\mathcal{E}_T(U_λ)$ is the same as that used to compute the stationary solution branches. This therefore leads to a very small increase in CPU time to evaluate the indicator as compared with the CPU time required to compute the stationary solution. On the other hand, when solving a nonlinear problem with ANM, the solutions are continuously known, as is the indicator. The indicator is then determined along the whole solution branch and not just for some points as the quantity Re(ζ) with the indirect method mentioned above.

In the case of a Hopf bifurcation, the indicator depends on the path following parameter “a” and also on the angular frequency $\omega$. So Tri et al. [10] have used a perturbation method with these two parameters ($a$ and $\omega$). However this method is very difficult to apply and it is almost impossible to build an automatic method for detection of bifurcation points in this way. To overcome these two drawbacks, we propose another way here, which consists of setting the Reynolds number and taking the angular frequency as the perturbation parameter. In this way we can easily use a perturbation method and also a continuation method (Cochelin [11]) to compute the indicator for each value of the angular frequency.

The first part of this paper is devoted to the determination of stationary bifurcation points on stationary branches, which are solutions of the Navier–Stokes equations. The second part proposes the determination of the first Hopf bifurcation on the stationary solution branches of the Navier–Stokes equations.

## 2. Bifurcation indicators for fluid mechanics

We consider a viscous incompressible fluid whose motion is governed by the stationary Navier–Stokes equations:

$$ \begin{cases} -\nu u_{i,jj} + u_j u_{i,j} + \frac{1}{\rho} p_{,i} = 0 & \text{in } \Omega, \\ u_{i,i} = 0 & \text{in } \Omega \end{cases} \quad (2) $$
where $v$ is the kinematic viscosity, $\rho$ is the density, and where the velocity $u$ and the pressure $p$ are the unknowns.
The boundary conditions are:

$$ u = \lambda u_d \quad \text{on } \partial_u \Omega \qquad (3) $$

where $u_d$ is the prescribed velocity on the boundary $\partial_u \Omega$ of the fluid domain $\Omega$. The parameter $\lambda$ can be identified as the Reynolds number. It should be noted that only Dirichlet boundary conditions are considered here. However Neumann boundary conditions could easily be used. The problem is written in the following operator form:

$$ \begin{cases} L(U) + Q(U, U) = 0, \\ u = \lambda u_d \quad \text{on } \partial_u \Omega \end{cases} \qquad (4) $$

where $L(U)$ represents a linear operator, function of $U$ (Laplacian and divergence operator), $Q(U, U)$ a quadratic operator (convection term), and $U$ a mixed unknown vector defined as:

$$ U = \begin{Bmatrix} u \\ p \end{Bmatrix}. \qquad (5) $$

## 2.1. Indicator for stationary bifurcation

In this section, we present a numerical method for determining stationary bifurcation points. A knowledge of the fundamental solution of the Navier-Stokes equations, $U_\lambda$, is assumed. This steady solution can be determined, for example, with an ANM [12,13] or with the incremental-iterative Newton-Raphson method. ANM consists of applying a perturbation technique in a stepwise manner [11]. The solution is represented by a sequence of local polynomial approximations. The first method is used in our approach since less CPU time is needed to evaluate the solution branch than with the Newton-Raphson method [13]. Moreover, the ANM solution curves are determined continuously and not point by point as with the incremental-iterative method.

The fundamental solution $U_\lambda$ verifies the following equations:

$$ \begin{cases} L(U_\lambda) + Q(U_\lambda, U_\lambda) = 0, \\ u_f = \lambda u_d \quad \text{on } \partial_u \Omega. \end{cases} \qquad (6) $$

This fundamental solution is perturbed by a load vector $\mu.f$ where $\mu$ is an unknown scalar and $f$ is a known random vector. Thus a fluctuation in the velocity field $\Delta V$ takes place. Let us denote $U = U_\lambda + \Delta V$. Inserting this expression in the stationary Navier-Stokes equations (4) and neglecting second-order terms in $\Delta V$, we obtain the following new problem:

$$ \begin{cases} L_T(\Delta V) = \mu f, \\ \Delta V = 0 \quad \text{on } \partial_u \Omega. \end{cases} \qquad (7) $$

Here $L_T(\bullet)$ is the tangent operator defined by the expression: $L_T(\bullet) = L(\bullet) + Q(\bullet, U_\lambda) + Q(U_\lambda, \bullet)$.

A supplementary condition has to be imposed to obtain a well-formulated problem. Here we require that the increment $\Delta V - \Delta V_0$ is orthogonal to $\Delta V_0$, which gives the following relationship:

$$ \langle \Delta V - \Delta V_0, \Delta V_0 \rangle = 0 \qquad (8) $$

with $\langle \bullet, \bullet \rangle$ being the Euclidian scalar product, and $\Delta V_0$ the solution of the linear system (7) with $\mu = 1$:

$$ \begin{cases} L_T(\Delta V_0) = f, \\ \Delta V_0 = 0 \quad \text{on } \partial_u \Omega. \end{cases} \qquad (9) $$

Finally, the system to be solved is written as:

$$ \begin{cases} L_T(\Delta V) = \mu f, \\ \langle \Delta V - \Delta V_0, \Delta V_0 \rangle = 0, \\ \Delta V = 0 \quad \text{on } \partial_u \Omega. \end{cases} \qquad (10) $$
The scalar $\mu$ and the vector $\Delta V$ are computed using:

$$ \left\{ \begin{array}{l} \mu = \frac{\langle \Delta V_0, \Delta V_0 \rangle}{\langle [K_T(U_\lambda)]^{-1} \cdot \{f\}, \Delta V_0 \rangle}, \\ \Delta V = \mu [K_T(U_\lambda)]^{-1} \cdot \{f\}, \\ \Delta V = 0 \quad \text{on } \partial_u \Omega \end{array} \right. \qquad (11) $$

where $[K_T(U_\lambda)]$ is a discrete matrix corresponding to the tangent operator $L_T(\Delta V)$.

Note that the resolution of system (11) demands, the triangulation of the matrix $K_T(U_\lambda)$ for each Reynolds number. So that a considerable amount of CPU time is needed to evaluate it, when there is a high number of unknowns. A perturbation method can be used to reduce this CPU time. This consists of using asymptotic expansions of both indicator $\mu$ and the vector $\Delta V$. This method has been applied in structural mechanics [8,16] and fluid mechanics [15] and leads to a very small increase in computing time because the operator needed to evaluate $\mu$ and $\Delta V$ is the same as that needed to compute the fundamental solution. Moreover, if the fundamental solution is computed using polynomial approximations, it is then an easy and economical way (in terms of CPU time) of using power series of the path parameter “a” to determine the scalar $\mu$ and the vector $\Delta V$ (see [8,16]).

The scalar $\mu$ represents our stationary bifurcation indicator and is defined at each point of the fundamental solution of Eq. (4). When the operator $L_T(\Delta V)$ is singular, the indicator is equal to zero. In elastic structural problems, this indicator can be identified as a stiffness measure. A stationary bifurcation point then corresponds to a disappearing stiffness (Boutyour et al. [8]).

## 2.2. Indicator for Hopf bifurcation

We now consider the nonsteady Navier-Stokes equations which are also written in terms of operators here:

$$ \left\{ \begin{array}{l} M(\dot{U}) + L(U) + Q(U, U) = 0, \\ U = \lambda U_d \quad \text{on } \partial_u \Omega \end{array} \right. \qquad (12) $$

where $M$ is the mass matrix. A perturbation $\Delta V$ of the stationary solution $U_\lambda$ is inserted in Eq. (12). Disregarding second-order terms in $\Delta V$, Eq. (12) becomes:

$$ \left\{ \begin{array}{l} M(\dot{\Delta V}) + L_T(\Delta V) = 0, \\ \Delta V = 0 \quad \text{on } \partial_u \Omega \end{array} \right. \qquad (13) $$

with $L_T(\Delta V)$ being the tangent operator defined at the fundamental point $U_\lambda$ (see Section 2.1). In the case of Hopf bifurcation, we are looking for points where the tangent operator $L_T(\Delta V)$ has a complex pair of eigenvalues which cross the imaginary axis. Let us introduce:

$$ \Delta V(t) = \Delta V \cdot e^{i\omega t} \qquad (14) $$

where $\omega$ is the angular frequency. Relationship (14) is inserted in Eq. (13), and gives:

$$ \left\{ \begin{array}{l} L_T(\Delta V) + i\omega M \Delta V = 0, \\ \Delta V = 0 \quad \text{on } \partial_u \Omega \end{array} \right. \qquad (15) $$

where $\Delta V$ is a complex vector whose real part is written as $\Delta V^a$ and its imaginary part as $\Delta V^b$. We change the right-hand side of Eq. (15) by $\mu f$:

$$ \left\{ \begin{array}{l} L_T(\Delta V) + i\omega M \Delta V = \mu f, \\ \Delta V = 0 \quad \text{on } \partial_u \Omega \end{array} \right. \qquad (16) $$

where $\mu$ is still a scalar, and $f$ is a random load vector, assumed to be real. The scalar $\mu$ is our bifurcation indicator, which has the properties of always being positive and null at a bifurcating point (Bensaadi [9], Cadou [17]). System (15) is written in the real domain:

$$ \begin{bmatrix} K_T(U_\lambda) & -\omega M \\ \omega M & K_T(U_\lambda) \end{bmatrix} \begin{bmatrix} \Delta V^a \\ \Delta V^b \end{bmatrix} = \begin{bmatrix} \mu f \\ 0 \end{bmatrix}. \qquad (17) $$
Finally, to compute the Hopf bifurcation points, the following system has to be solved:

$$ \begin{cases} \ell_T(U_\lambda, \omega) \cdot \Delta V = \mu f, \\ \Delta V = 0 \quad \text{on } \partial\Omega_u \end{cases} \qquad (18) $$

with $\ell_T(U_\lambda, \omega)$ being an operator defined at a point on the stationary solution branch ($U_\lambda$) and for a given value of the pulsation ($\omega$). This operator is twice the size of the stationary equations. Note that Eqs. (18) and (7) are quite similar.

The next step is to find where in the fundamental branch $U_\lambda$, the scalar $\mu$ is null. System (18) is equivalent to a linear eigenvalue problem when $\mu$ is equal to zero. In this case, the angular frequency $\omega$ is the imaginary part of the critical eigenvalues and $\Delta V$ is the corresponding eigenvector.

To ensure the uniqueness of the solution, an additional condition is needed. This condition is introduced here by specifying a value for the norm of the vector $\Delta V$:

$$ \|\Delta V\|^2 = \|\Delta V_0\|^2. \qquad (19) $$

The vector $\Delta V_0$ is a solution of the problem:

$$ \Delta V_0 = \ell_T^{-1}(U_\lambda, 0) \cdot f. \qquad (20) $$

Condition (19) leads to the following definition of the bifurcating indicator:

$$ \mu = \sqrt{\frac{\langle \Delta V_0, \Delta V_0 \rangle}{\langle U^*, U^* \rangle}} \qquad (21) $$

where

$$ U^* = \ell_T^{-1}(U_\lambda, \omega) \cdot f. \qquad (22) $$

Finally, the complex vector $\Delta V$, solution of Eq. (18), is given by:

$$ \Delta V = \mu U^*. \qquad (23) $$

It can be noted that condition (8) in Section 2.1 can also be used to define the indicator $\mu$:

$$ \langle \Delta V - \Delta V_0, \Delta V_0 \rangle = 0. \qquad (24) $$

In this case the bifurcating indicator is defined by the following expression:

$$ \mu = \frac{\langle \Delta V_0, \Delta V_0 \rangle}{\langle U^*, \Delta V_0 \rangle}, \quad \text{with condition (24).} \qquad (25) $$

This second condition is very interesting because it leads to a linear expression of the bifurcating indicator $\mu$. Unfortunately, the numerical results presented in the following sections show the difficulties of using such an expression with a perturbation method.

### 3. Numerical implementation

In this section, we will first discuss the numerical method we use to compute the stationary and Hopf indicator. We then discuss the best choice for the additional conditions (expressions (19) and (24)) in order to compute the Hopf bifurcation indicator with a perturbation technique. Finally, we propose a continuation method for computing the indicator of Hopf bifurcation for all values of the angular frequency. Discussion and results for the stationary indicator will be given in Section 4.1.

#### 3.1. Spatial discretization

The numerical tests presented here are obtained using the finite element method. A quadrilateral element is used (Q9/3D [18]) with nine nodes for the velocity and three for the pressure. A penalty method [18] is used in the discrete continuity equation. Velocities are interpolated by biquadratic functions, whereas discontinuous linear functions are used for the pressure. The use of the finite element method leads to a discrete mass matrix $M$ and a tangent matrix $K_T(U_\lambda)$ (resulting from the discretization of the tangent operator $L_T(U_\lambda)$).
3.2. *Comparison of the two additional conditions (19) and (24)*

In this section, the two indicators defined in Eqs. (21) and (25) are compared. The numerical test is performed considering the flow around a cylinder. The configuration and the boundary conditions are presented in Fig. 1. The Reynolds number of the flow is determined from the inlet velocity and the diameter of the cylinder. We first propose a “direct computation” of the indicators and of the corresponding vectors $\Delta V$. For a given value of the angular frequency $\omega$, systems (25) and (21) are solved. A given $\omega$ together with the triangulation of the operator $\mathcal{E}_T(U_\lambda, \omega)$ provides the corresponding values of the Hopf indicators $\mu$ and the complex vector $\Delta V$. The results for each indicator and for a Reynolds number in the region of 52 are represented as a function of the angular frequency in Fig. 2(a) for system (25) and in Fig. 2(b) for system (21). The algorithm for the direct computation can be summarized in the following steps:

① Given initial values of $\omega_0$ (usually equal to zero) and $U_\lambda$
② Compute $\Delta V_0 = \mathcal{E}_T^{-1}(U_\lambda, 0) \cdot f$
③ Evaluate $\omega = \omega + \Delta\omega$, where $\Delta\omega$ is a known increment
④ Compute $\Delta V$ and $\mu$ by solving system (21) or (25)

5a If $\omega = \omega_{\max}$
then
– change $U_\lambda$ and goto ②
– or stop algorithm

5b Else goto ③

Note that this procedure is an incremental method and the solution is then known point by point. A small increment $\Delta\omega$ gives an excellent description of the curves $(\omega, \mu)$ but involves some additional CPU time. The most time-consuming step in the direct computation is step 4 because it requires the triangulation of the operator $\mathcal{E}_T(U_\lambda, \omega)$.

Another difficulty may arise when definition (25) is used to evaluate the indicator $\mu$. For some values of the angular frequency $\omega$, the vector $U^*$ (with $U^* = \mathcal{E}_T^{-1}(U_\lambda, \omega) \cdot f$) can be orthogonal to the initial vector $\Delta V_0$. In this case, which remains an exception, the scalar $\mu$ is infinite (see Eq. (25)). In Fig. 2(a) such a point is reached for $\omega$ in the region of 5. Clearly, for this latter angular frequency, an infinite value of indicator $\mu$ is computed with expression (25), resulting from the zero value of the denominator of the fractional expression (25). On the other hand, expression (21) does not lead to such singular points. Thus, with this latter expression the indicator $\mu$ is well-defined for any angular frequencies.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Flow around a cylinder, description of the geometry.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 2. Indicator versus angular frequency, flow around a cylinder with Reynolds number equal to 52.</figcaption></figure>

The initial $\omega_0$ value is generally chosen to be zero. However, for some examples, when the critical angular frequency $\omega_c$ is known, it is more appropriate to take an initial value $\omega_0$ near this critical number. For the flow around the cylinder, the critical values of the Reynolds number and the angular frequency (through the Strouhal number) are now well known ($Re_c \approx 46–47$, $\omega_c \approx 4$). In this special test, it is helpful to reduce the total computing time of the direct computation by choosing an initial value $\omega_0$ close to this critical value. However, in almost all studies, these critical flow parameters are not known. So the use of an incremental method is not the best choice. Moreover, for solving problem (21)–(23) which only depends on one parameter $\omega$, a perturbation method [19,12] is better than any iterative method, as for example, the Newton–Raphson method. Indeed, with this method, only one triangulation of the operator $\mathcal{L}_T(U_\lambda, \omega)$ is needed in order to compute an analytical part of the $(\mu, \omega)$ curve. In the next paragraphs we will review the perturbation and continuation techniques for calculating this $(\mu, \omega)$ curve.

### 3.3. Perturbation technique

The basis of the perturbation method is the performance of a Taylor expansion in terms of the pulsation $\omega$ of the two unknowns $(\mu, \Delta V)$ near to a known solution $(\mu_0, \Delta V_0)$:

$$ \begin{cases} \Delta V = \Delta V_0 + \omega \Delta V_1 + \omega^2 \Delta V_2 + \cdots + \omega^p \Delta V_p, \\ \mu = \mu_0 + \omega \mu_1 + \omega^2 \mu_2 + \cdots + \omega^p \mu_p. \end{cases} \qquad (26) $$

These asymptotic expansions are then inserted in Eqss. (16) and (19). This latter condition (19) is chosen instead of (24) so as to ensure the continuity of the $(\mu, \omega)$ curves. This property is compulsory when a perturbation technique is used. With the help of balancing terms with identical powers for $\omega$, we obtain the following set of equations:

**Order 0 in $\omega$:**

$$ \begin{cases} K_T(U_\lambda) \cdot \Delta V_0 = f, \\ \mu_0 = 1, \\ \Delta V_0 = 0 \quad \text{on } \partial \Omega_u. \end{cases} \qquad (27) $$

**Order 1 in $\omega$:**

$$ \begin{cases} K_T(U_\lambda) \cdot \Delta V_1 + iM \cdot \Delta V_0 = \mu_1 f, \\ \langle \Delta V_1, \Delta V_0 \rangle = 0, \\ \Delta V_1 = 0 \quad \text{on } \partial \Omega_u. \end{cases} \qquad (28) $$
Order p in $\omega$:

$$ egin{cases} K_T(U_\lambda) \cdot \Delta V_p + iM \cdot \Delta V_{(p-1)} = \mu_p f, \\ 2 \cdot \langle \Delta V_p, \Delta V_0 \rangle + \sum_{r=1}^{p-1} \langle \Delta V(p-r), \Delta V_r \rangle = 0, \\ \Delta V_p = 0 \quad \text{on } \partial\Omega_u. \end{cases} \tag{29} $$

Finally, at each order, the indicator and the complex vector could be computed from the following expressions:

$$ egin{cases} \Delta V_0 = [K_T(U_\lambda)]^{-1} \cdot f, \\ \mu_0 = 1, \\ \Delta V_0 = 0 \quad \text{on } \partial\Omega_u. \end{cases} \tag{30} $$

Since $f$ has been chosen as being real, the $\Delta V_0$ vector is also real. This condition on the $f$ random vector leads to some remarkable properties concerning the indicator and the complex vector at each order. For instance, by inserting expression (30) at order 0 in Eqs. (28), the indicator and complex vector at order 1 become:

$$ egin{cases} \Delta V_1 = -i[K_T(U_\lambda)]^{-1} \cdot M \cdot \Delta V_0, \\ \mu_1 = 0. \end{cases} \tag{31} $$

At order 1, the indicator is equal to zero and the corresponding vector is purely imaginary.

$$ egin{cases} \Delta V_2 = \mu_2 \Delta V_0 + \Delta V_2^*, \\ \mu_2 = \frac{-\langle \Delta V_1, \Delta V_1 \rangle - 2\langle \Delta V_2^*, \Delta V_0 \rangle}{2 \cdot \langle \Delta V_0, \Delta V_0 \rangle}, \\ \text{with } \Delta V_2^* = [K_T(U_\lambda)]^{-1} \cdot M \cdot \Delta V_1. \end{cases} \tag{32} $$

At order 2, the indicator is computed with the use of expressions (32) and the complex vector at this order is real. Finally, a general expression for both indicator and complex vector can be made for any order $p$:

Order $p$, with $p$ even:

$$ egin{cases} \Delta V_p = \mu_p \Delta V_0 + \Delta V_p^*, \\ \mu_p = \frac{-2\langle \Delta V_p^*, \Delta V_0 \rangle - \sum_{r=1}^{p-1} \langle \Delta V_{(p-r)}, \Delta V_r \rangle}{2 \cdot \langle \Delta V_0, \Delta V_0 \rangle}, \\ \Delta V_p^* = [K_T(U_\lambda)]^{-1} \cdot M \cdot \Delta V_{(p-1)}. \end{cases} \tag{33} $$

Order $p$ with $p$ odd:

$$ egin{cases} \Delta V_p = i \Delta V_p^*, \\ \mu_p = 0, \\ \Delta V_p^* = -[K_T(U_\lambda)]^{-1} \cdot M \cdot \Delta V_{(p-1)}. \end{cases} \tag{34} $$

Previous expressions (33) and (34) are used to build series (26) which are written in the following forms:

$$ egin{cases} \Delta V = \{\Delta V_0\} + \omega\{i\Delta V_1^*\} + \omega^2\{i\Delta V_2^*\} + \cdots + \omega^{(p-1)}\{i\Delta V_{(p-1)}^*\} + \omega^p\{\Delta V_p\}, \\ \mu = \mu_0 + \omega^2\mu_2 + \omega^4\mu_4 + \cdots + \omega^p\mu^p \quad \text{with } p \text{ even.} \end{cases} \tag{35} $$

This method has been applied to compute the indicator and the complex vector for the flow around the cylinder. The results for several truncation orders of asymptotic expansions and for a Reynolds number equal to 46 are plotted and compared with the direct method (Section 3.2) in Fig. 3. It is shown that the series have a finite radius of convergence. The asymptotic representation of the solution is only valid up to a certain value of the angular frequency $\omega$. For the highest truncation order (order 22), the asymptotic expansion seems to be valid up to a value of $\omega$ equal to 1. Since the “area of interest” of this test is located near $\omega \approx 4$, this first step of the perturbation technique is not sufficient to study the stability of this example. In order to determine the indicator for any angular frequency values, we propose applying the perturbation method step by step. This method, called the continuation technique [11], has been successfully applied for the computation of the complete solution branches of nonlinear problems (nonlinear elasticity [11] or Navier-Stokes equations [13]).
<figure><img src="image_5.png" /><figcaption>Fig. 3. Bifurcation indicator versus angular frequency—comparison of the polynomial approximations at several orders of truncation and the reference (direct computation with Eq. (21)), for flow around a cylinder at Reynolds number equal to 46.</figcaption></figure>

### 3.4. Continuation technique

The continuation technique consists of applying the perturbation technique from a new starting point ($\mu_0$, $\Delta V_0$). First the range of validity of approximation (35) is computed automatically using the following criterion [11]:

$$ \omega_m = \left( \varepsilon \frac{\| \Delta V_1 \|}{\| \Delta V_p \|} \right)^{1/(p-1)} \qquad (36) $$

where $\varepsilon$ is a small number and $p$ the truncation order of the series expansion. The end of the solution path ($\omega, \mu(\omega), \Delta V(\omega) = (\omega_m, \mu_m, \Delta V_{\omega_m})$ is defined by criterion (36) and expansions (35) and this point is taken as the new starting point ($\omega_0, \mu_0, \Delta V_{\omega_m}$). The angular frequency then becomes:

$$ \omega = \omega_0 + \hat{\omega} \qquad (37) $$

where $\hat{\omega}$ is the new perturbation parameter. The unknowns ($\mu, \Delta V$) are again expanded into power series and balancing identical powers in $\hat{\omega}$ gives a new set of linear equations:

Order 0 in $\hat{\omega}$:

$$ \left\{ \begin{aligned} \mu_0 &= \sqrt{\frac{\langle \Delta V_0^i, \Delta V_0^i \rangle}{\langle \Delta V_0^*, \Delta V_0^* \rangle}}, \\ \Delta V_0^* &= \xi_T (U_\lambda, \omega_0)^{-1} \cdot f, \\ \Delta V_0 &= \mu_0 \Delta V_0^* \end{aligned} \right. \qquad (38) $$

where $\Delta V_0^i$ is the solution of the first step of the method (i.e. $\Delta V_0^i$ is solution of Eq. (30)).

Order 1 in $\hat{\omega}$:

$$ \left\{ \begin{aligned} \mu_1 &= -\frac{\langle \Delta V_1^*, \Delta V_0 \rangle}{\langle \Delta V_0^*, \Delta V_0 \rangle}, \\ \Delta V_1^* &= i\xi_T (U_\lambda, \omega_0)^{-1} \cdot M \cdot \Delta V_0, \\ \Delta V_1 &= \mu_1 \Delta V_0^* + \Delta V_1^*. \end{aligned} \right. \qquad (39) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Bifurcation indicator versus angular frequency—comparison of the continuation technique solution and the reference (direct computation with Eq. (21)), for flow around a cylinder at a Reynolds number of 52.</figcaption></figure>

Order p in $\hat{\omega}$:

$$ egin{cases} 
\mu_p = \frac{-2 \cdot \langle \Delta V_p^*, \Delta V_0 \rangle - \sum_{r=1}^{p-1} \langle \Delta V(p-r), \Delta V_r \rangle}{2 \cdot \langle \Delta V_0^*, \Delta V_0 \rangle}, \\ 
\Delta V_p^* = i \xi_T (U_\lambda, \omega_0)^{-1} M \cdot \Delta V_{p-1}, \\ 
\Delta V_p = \mu_p \Delta V_0^* + \Delta V_p^* 
\end{cases} 
\quad (40) $$

The linear systems (38)–(40) have the same tangent operator $\xi_T(U_\lambda, \omega_0)$. Only one triangulation of this operator (at order 0, Eq. (38)) leads to the determination of an analytical part of the curve $(\mu, \omega)$. The following problems (40) need only backward and forward substitutions. It should be noted that problems (38)–(40) can also be solved using an iterative solver (like GMRES). But in this case, the CPU times required to solve all these problems are the same. Finally, the total CPU times for computing all the terms of the asymptotic expansions (35) can be greater with an iterative solver than with a direct solver even if the problem involves a large number of unknowns.

Two different ways of calculating the starting point $(\mu_0, \Delta V_0)$ at each step of the continuation technique can be deduced from the property of the common tangent operator $\xi_T(U_\lambda, \omega_0)$. The first method has been discussed above (by inserting the optimal value $\omega_m$ in the asymptotic expansions (35)). The second way is to compute values $(\mu_0, \Delta V_0)$ for the new value of the angular frequency $\omega_0$ by solving system (38) for the order 0. This second method enables us to avoid accumulations of errors due to successive polynomial approximations (with method one) and has been chosen in the numerical tests presented in this paper. The continuation technique has been applied to the example of the flow around a cylinder. A comparison of this method and direct computation is given in Fig. 4 for a Reynolds number of 52.

The final algorithm for precise determination of the Hopf bifurcating point is as follows:

0. Set the parameters of the method:

(a) The area of interest for the Reynolds number: $Re^c \in [Re^1, Re^2]$.

(b) The order of truncation of the asymptotic expansions: $p$.

(c) The number of steps in the continuation technique.

(d) The minimum value of the indicator: $\epsilon$.

1. Compute the stationary solution $U_\lambda$ and the corresponding Reynolds number, $Re$, with ANM [13].

2. If $Re \in [Re^1, Re^2]$, then:

2.1 Computation of the $\mu_i$ and $\Delta V_i$:
(i) first step with the expressions (30), (33) and (34).

(ii) following steps with the expressions (38) and (40).

2.2 Evaluation of the range of validity of asymptotic expansions with the criterion (36).

2.3 Building of series with the expression (35), for the first step, or with the expression (26) for the following steps.

2.5 For this Reynolds number, determine the value of the angular frequency ω where the indicator μ is minimum, ($\omega_{mini}^{Re}, \mu_{mini}^{Re}$):

2.6.1 First step: ($\omega_{mini}, \mu_{mini}$) = ($\omega_{mini}^{Re}, \mu_{mini}^{Re}$)

2.6.2 Following steps: If $\mu_{mini}^{Re} < \mu_{mini}$ Then ($\omega_{mini}, \mu_{mini}$) = ($\omega_{mini}^{Re}, \mu_{mini}^{Re}$)

3.1 If $\mu_{mini} < \epsilon$ then

3.1 $Re^c = Re, \omega^c = \omega_{mini}$

3.2 Stop

3.2 Else goto 1.

# 4. Applications and discussion

In this section, we will consider several traditional fluid stability tests to assess the ability of the methodology that has been presented to predict stationary and periodic bifurcations.

## 4.1. Computation of stationary bifurcation points in fluid flows

We will first apply the stationary indicator defined in Section 2.1 to the example of a two-dimensional symmetric flow through a channel with a symmetrical expansion about its centerline. Shapira et al. [20] have proposed a linear stability analysis with computation of eigenvalues for several channel width values. In this paper, we consider a configuration that Shapira has also investigated: a channel expansion ratio of 2 and an expansion angle of $\alpha = 90^{\circ}$ (see Fig. 5). The configuration and the boundary conditions are given in Fig. 5. A parabolic velocity profile is imposed at the inlet channel (in the x direction). The Reynolds number is based on the channel diameter (d), the maximum velocity at the inlet and the kinematic viscosity $\nu$:

$$Re = \frac{U_{\max} \cdot d}{\nu} \qquad (41)$$

The mesh used in the numerical computations is made of 642 elements (Q9/3D) (roughly 6000 degrees of freedom for velocity). The stationary solution is computed with an Asymptotic Numerical Method [13]. This method is an automatic procedure which allows a nonlinear solution branch to be followed. Such a method can lead to a step accumulation in the continuation procedure near a turning point or a bifurcating point. To illustrate this idea, the velocity $U_x$ measured at the location P1 (see Fig. 5) is plotted versus the Reynolds number in Fig. 6(a). In this plot, a step accumulation is observed for a Reynolds number close to 200 (each step in the continuation procedure is represented by a square symbol on curve Fig. 6(a). We then apply the procedure for detection of stationary bifurcation

<figure><img src="image_7.png" /><figcaption>Fig. 5. Configuration and boundary conditions for flow in a channel.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Study of stability of flow through a channel with a symmetrical expansion.</figcaption></figure>

points given in Section 2.1 for a Reynolds number between 200 and 250. In Fig. 6(b) the indicator is plotted versus the Reynolds number. This indicator is equal to zero for a Reynolds number value between 214 and 215. For the same example, Shapira et al. [20] have found that the eigenvalue which corresponds to the least stable mode is real and changes sign for a Reynolds number between 212.2 and 216 (the eigenvalue is equal to -0.637E-3 for $Re = 212.4$ and 0.1E-3 for $Re = 216$).

This indicates that for a Reynolds number greater than 216, several other steady solutions exist. Thus, to compute these stationary solution branches after the bifurcating point, we performed computations with several parameters of the ANM. The streamlines of the flow for three different computations with the ANM are given in Fig. 7 for a Reynolds number equal to 600. Another interesting feature of the ANM is the ability to jump near bifurcation points to different bifurcation branches (see Cochelin [11], Boutyour et al. [8] for further explanations) by modifying only some parameters of the continuation method.

Our indicator therefore gives accurate values for the critical Reynolds number where a bifurcation appears. The eigenmode corresponding to the bifurcation can be computed either by applying linear interpolation or using an automatic procedure (Vannucci et al. [16]). In this paper, an easy automatic algorithm, based on ANM, is proposed for calculating and switching on stationary branches emanating from a bifurcating point in solid mechanics problems.

It should be noted that this indicator has been applied for studying the Coanda effect which occurs in several flow configurations [21] (sudden expansion with or without divergent, open cavity). This work shows a high degree of consistency between the results obtained with the stationary indicator and the results found in the literature.

## 4.2. Computation of Hopf bifurcation points in fluid flows

### 4.2.1. Flow around a cylinder

In this section, we will apply our Hopf bifurcation indicator to a traditional test: flow around a cylinder. The configuration and the boundary conditions are given in Fig. 1 and are the same as in Jackson's study [7]. The Reynolds number is calculated taking the cylinder diameter, the kinematic viscosity and the inlet velocity into consideration. The mesh used for all computations has more than 3000 velocity nodes. A Galerkin weighting was used in all numerical tests.

In the case of stationary bifurcation, for each Reynolds number, an indicator value was computed. In the case of Hopf bifurcation, for each Reynolds number, we compute an indicator versus the angular frequency curve, written $G_{Re}(\mu, \omega)$. We know for this example that the critical Reynolds number is between 40 and 50. So computations have been performed in this range of Reynolds numbers. The results are given in Fig. 8 where the indicator is plotted versus the angular frequency for several values of Reynolds numbers. A Hopf bifurcation is characterized by a zero indicator value. In fact, looking for some critical parameters (Reynolds number and angular frequency) whose indicator is zero, consists of detection of the minimum of function $G_{Re}(\mu, \omega)$. In Table 1 the minima of this function are given for five
<figure><img src="image_9.png" /><figcaption>Fig. 7. Flow in a channel—streamlines at $Re = 600$.</figcaption></figure>

<table><caption>Table 1
Evolution of the minimum of function $G_{Re}(\mu, \omega)$ versus Reynolds number, for flow around a cylinder</caption><thead><tr><th>Re</th><th>38.00</th><th>43.83</th><th>45.98</th><th>46.13</th><th>47.53</th></tr></thead><tbody><tr><td>μ<sub>mini</sub></td><td>0.903 × 10<sup>-1</sup></td><td>0.184 × 10<sup>-1</sup></td><td>1 × 10<sup>-8</sup></td><td>0.112 × 10<sup>-2</sup></td><td>0.11 × 10<sup>-1</sup></td></tr><tr><td>ω(μ<sub>mini</sub>)</td><td>3.195</td><td>3.769</td><td>3.972</td><td>3.986</td><td>4.117</td></tr><tr><td>Number of steps</td><td>8</td><td>10</td><td>30</td><td>18</td><td>10</td></tr></tbody></table>

Reynolds numbers. The smallest minimum for these five Reynolds numbers occurs for $Re = 45.98$ and $\omega = 3.972$ ($G_{Re}=45.98$ ($\mu = 1 \times 10^{-8}$, $\omega = 3.972$)). Since the critical angular frequency is found, it is then easy to compute the eigenvalue and the corresponding eigenmode for this point of bifurcation. Indeed, when $\mu$ is zero then (16) becomes:

$$ \begin{cases} L_T(\Delta V) + i\omega M \Delta V = 0, \\ \Delta V = 0 \quad \text{on } \partial_{\mu}\Omega. \end{cases} \qquad (42) $$

Therefore, $i\omega$ is the eigenvalue and $\Delta V$ is the eigenmode. In this test, the eigenvalues corresponding to the Hopf bifurcation are $\zeta_1 = \bar{\zeta}_2 = 3.972i$. The corresponding eigenmode is simply computed with Eq. (26) with $\omega = \omega_c = 3.972$:

$$ \Delta V = \Delta V_0 + \omega_c \Delta V_1 + \omega_c^2 \Delta V_2 + \cdots + \omega_c^p \Delta V_p. \qquad (43) $$

With these asymptotic expansions (43), there is no singular matrix at a bifurcating point because the operator $\mathcal{L}_T(U_\lambda, \omega)$ is triangulated far from the singular point. In Fig. 9, streamlines for real 9(a) and imaginary parts 9(b)
<figure><img src="image_10.png" /><figcaption>(a) Indicator versus angular frequency for several Reynolds numbers</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Zoom of Fig. (8(a))</figcaption></figure>

<caption>Fig. 8. Computations of indicator versus angular frequency for different Reynolds numbers, for flow around a cylinder.</caption>

of the bifurcating eigenvector are drawn. The Strouhal number corresponding to this bifurcating point is equal to $St_c = 0.137$. This number is computed from the expression:

$$St_c = \frac{D \cdot \omega_c}{2\pi U x_c} \qquad (44)$$

where the subscript c indicates critical parameters. Jackson [7] found the following critical parameters for this test: $Re_c = 46.184$, $St_c = 0.138$. These are similar to the values we found.

### 4.2.2. Flow around a square cylinder

We now consider the flow around a square cylinder. The geometry and the boundary conditions are given in Fig. 10. The mesh used leads to 3098 velocity nodes for 725 elements Q9/3D. This test is the same as that proposed by
<figure><img src="image_12.png" /><figcaption>Fig. 9. Streamlines of the complex bifurcating eigenvector at $Re = 45.98$ and $\omega = 3.972$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Geometry and boundary conditions, for flow around a square cylinder.</figcaption></figure>

<table><caption>Table 2
Evolution of the minimum of function $G_{Re}(\mu, \omega)$ versus Reynolds number, flow around a square cylinder</caption><thead><tr><th>Re</th><th>41.50</th><th>50.64</th><th>54.70</th><th>55.00</th><th>56.53</th></tr></thead><tbody><tr><th>μ<sub>mini</sub></th><td>0.19</td><td>0.47 × 10<sup>-1</sup></td><td>2.65 × 10<sup>-7</sup></td><td>0.37 × 10<sup>-2</sup></td><td>0.18 × 10<sup>-1</sup></td></tr><tr><th>ω(μ<sub>mini</sub>)</th><td>3.5</td><td>4.30</td><td>4.65</td><td>4.68</td><td>4.8</td></tr><tr><th>Number of steps</th><td>7</td><td>8</td><td>16</td><td>13</td><td>8</td></tr></tbody></table>

Kelkar and Patankar [22]. The Hopf bifurcation indicator has been computed for a Reynolds number between 40 and 60. This Reynolds number is based on the length of the side of the square cylinder and also on the uniform inlet velocity. In Table 2, the minima of function $G_{Re}(\mu, \omega)$ are given for five Reynolds numbers. The indicator is almost null ($\mu = 2.65 \times 10^{-7}$) for a Reynolds number equal to 54.70 and an angular frequency equal to 4.65. The critical Strouhal number of this flow is computed using expression (44) in which $D$ is the length of the side of the square cylinder and the critical velocity is in the inlet velocity for a Reynolds number equal to 54.70. The critical Strouhal number is then $St_c = 0.135$. This value is consistent with the numerical results reported by Kelkar [22] ($St_c = 0.126$) and the experimental results of Okajima [23], close to 0.1 for a square cylinder.

### 4.2.3. Flow in a channel with a symmetric expansion

The two first tests are standard tests used in fluid stability analysis. With these tests our results have been compared with the results found in the literature. We will now study the stability of the flow in a channel with a symmetric expansion. This example has been studied in Section 4.1 and a stationary bifurcation has been found for a Reynolds number close to 215. We are now looking for a Hopf bifurcation from the bifurcated solutions that are displayed in Fig. 7(a). Computations have been made for some Reynolds numbers between 10 and 900. The Hopf indicator is then plotted versus the angular frequency for four Reynolds numbers in Fig. 11. With these four values, we can isolate
<figure><img src="image_14.png" /><figcaption>(a) Re = 70</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Re = 270</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Re = 550</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) Re = 830</figcaption></figure>

Fig. 11. Computations of the indicator versus the angular frequency for different Reynolds numbers, for flow in a channel with a symmetric expansion.

<table><caption>Table 3
Evolution of the minimum of function $\mathcal{G}_{Re}(\mu, \omega)$ versus the Reynolds number, for flow in a channel with a symmetric expansion</caption><thead><tr><th>Re</th><th>606</th><th>607</th><th>609</th><th>610</th></tr></thead><tbody><tr><td>μ<sub>mini</sub></td><td>5.017 × 10<sup>-3</sup></td><td>2.949 × 10<sup>-3</sup></td><td>8.888 × 10<sup>-4</sup></td><td>1.118 × 10<sup>-3</sup></td></tr><tr><td>ω(μ<sub>mini</sub>)</td><td>0.536</td><td>0.537</td><td>0.538</td><td>0.539</td></tr><tr><td>Number of steps</td><td>68</td><td>69</td><td>71</td><td>71</td></tr></tbody></table>

the zone of Reynolds numbers for which Hopf bifurcation can occur. Indeed, the indicator remains far from zero while the Reynolds number is less than 270. No Hopf bifurcation can be found in this area. On the other hand, the minima for the indicator tend to be zero for Reynolds numbers greater than 550. A first zone, for a Reynolds number between 606 and 610, is then isolated. The Hopf indicator for this zone becomes zero for an angular frequency value close to 0.5396, thus indicating a Hopf bifurcation for these critical parameter values ($Re_c = 609$ and $ω_c = 0.538$ see Table 3). These critical numbers lead to a Strouhal number of 0.14 (computed with expression (44) in which the length $D$ is the channel diameter and $U_{xc}$ is the maximum velocity at the inlet).

When the Reynolds number is increased, a second zone is found where the indicator is null. For Reynolds numbers between 617 and 621 (see Table 4), the indicator is then null for a critical Reynolds number of 620 and a critical angular frequency of 0.6667. The Strouhal number for these critical parameters is then $St_c = 0.17$.
<table><caption>Table 4
Evolution of the minimum of function $G_{Re}(\mu, \omega)$ versus the Reynolds number, for flow in a channel with a symmetric expansion</caption><thead><tr><th>Re</th><th>617</th><th>618</th><th>619</th><th>620</th><th>621</th></tr></thead><tbody><tr><td>μ<sub>mini</sub></td><td>8.203 × 10<sup>-3</sup></td><td>4.884 × 10<sup>-3</sup></td><td>1.595 × 10<sup>-3</sup></td><td>1.588 × 10<sup>-4</sup></td><td>1.672 × 10<sup>-3</sup></td></tr><tr><td>ω(μ<sub>mini</sub>)</td><td>0.6631</td><td>0.6646</td><td>0.6661</td><td>0.6667</td><td>0.6676</td></tr><tr><td>Number of steps</td><td>87</td><td>88</td><td>90</td><td>133</td><td>96</td></tr></tbody></table>

<table><caption>Table 5
First Hopf bifurcation in the 2-D lid-driven cavity. Comparison of the critical values obtained in the literature</caption><thead><tr><th>Authors</th><th></th><th>Re<sub>c</sub></th><th>Critical frequency</th></tr></thead><tbody><tr><td>Proposed method</td><td>13 122 d.o.f. mesh b Fig. 14(b)</td><td>7890</td><td>0.44</td></tr><tr><td>Poliashenko and Audin [26]</td><td></td><td>7763</td><td>0.45</td></tr><tr><td rowspan="4">Fortin et al. [4]</td><td>1 4082 d.o.f.</td><td>7745</td><td>0.45</td></tr><tr><td>18 678 d.o.f.</td><td>7937</td><td>0.45</td></tr><tr><td>25 814 d.o.f.</td><td>7998.5</td><td>0.45</td></tr><tr><td>31 922 d.o.f.</td><td>8000</td><td>0.45</td></tr><tr><td>Cazemier et al. [24,27]</td><td>DNS</td><td>7922</td><td>0.44</td></tr><tr><td>Cazemier et al. [24,27]</td><td>POD</td><td>7819</td><td>0.60</td></tr><tr><td>Tiensinga et al. [28]</td><td></td><td>8375</td><td>0.43</td></tr><tr><td rowspan="3">Peng et al. [25]</td><td>100 × 100 grid points</td><td>7402</td><td>0.59</td></tr><tr><td>150 × 150 grid points</td><td>7694</td><td>—</td></tr><tr><td>200 × 200 grid points</td><td>7704</td><td>—</td></tr><tr><td rowspan="3">Sahin and Owens [29]</td><td>129 × 129 grid points</td><td>8244</td><td>0.45</td></tr><tr><td>193 × 193 grid points</td><td>8109</td><td>0.45</td></tr><tr><td>257 × 257 grid points</td><td>8069</td><td>0.69</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>Fig. 12. Geometry for the 2-D lid-driven cavity flow.</figcaption></figure>

### 4.2.4. The 2-D lid-driven cavity flow

We will now consider the flow in a 2-D lid-driven cavity. The domain and boundary conditions of this problem are shown in Fig. 12. In the past decade there have been a large number of papers in the literature on the study of the stability of the 2-D lid-driven cavity. We limit ourselves here to the first Hopf bifurcation. A detailed description of the whole dynamic phenomenon that occurs in this example can be found in Cazemier [24] or Peng et al. [25].

In Table 5 we have given the critical Reynolds number and frequency for which the steady flow becomes periodic. The results given by the proposed method are compared with the results found in the literature. A large number of authors have studied this problem using different methods to characterize the critical parameters. For example, Poliashenko and Audin [26] use a direct method, Fortin et al. [4] compute the smallest eigenvalues of the matrix
resulting from the discretization of the Navier-Stokes equations. Cazemier et al. [27] propose a direct numerical simulation (DNS) for checking the results from a proper orthogonal decomposition (POD) coupled with an eigenvalue computation of the resulting low dimensional matrix. Tiesenga et al. [28] have used a Newton-Picard method which results in the determination of the eigenvalues. The results of Peng et al. [25] are obtained by a direct numerical simulation while those of Sahin and Owens [29] rely on Arnoldi's method.

For our computations, we have used a mesh with 80 × 80 elements, that leads to 13.122 d.o.f. This corresponds to the maximum size problem that can be processed with our computer resources. The corresponding mesh is drawn in Fig. 14(b). Despite this coarser mesh, when compared with those used by Fortin et al. [4] the critical Reynolds number we obtain is very similar to those in the literature ($Re_c \approx 8000$). The critical frequency result based on our indicator is about the same as those in the literature (0.45).

For this example, several authors have computed the eigenvalues of the Jacobian matrix for a Reynolds number close to 8000 (see for example Tiesenga et al. [28] or Fortin et al. [4]). These results are given in Table 6. For this Reynolds number, $Re_c = 7890$, the evolution of the indicator versus the frequency is plotted in Fig. 13. We have also given all the minima of function $G_{Re=7890}(\mu, \omega)$ in this figure. If we compare these values with the results presented in Table 6, it is easy to see that these minima correspond to eigenvalues near the imaginary axis.

<figure><img src="image_19.png" /><figcaption>Fig. 13. The 2-D lid-driven cavity. Indicator versus frequency, $Re = 7890$.</figcaption></figure>

<table><caption>Table 6
Eigenvalues near the imaginary axis from Refs. [28,4]</caption><caption>Eigenvalues and corresponding frequencies from</caption><thead><tr><th colspan="3">Tiesenga et al. [28]</th><th colspan="3">Fortin et al. [4]</th></tr><tr><th>λ</th><th>St</th><th>Re</th><th>λ</th><th>St</th><th>Re</th></tr></thead><tbody><tr><td>-0.0144±0.9392i</td><td>0.1495</td><td>8000</td><td>≈ 2.83786i</td><td>≈ 0.45</td><td>8000</td></tr><tr><td>-0.0268±1.8667i</td><td>0.2971</td><td>8000</td><td>≈ 1.95i</td><td>≈ 0.31</td><td>8000</td></tr><tr><td>0.0009±2.764i</td><td>0.4399</td><td>8375</td><td>≈ 1</td><td>≈ 0.16</td><td>8000</td></tr><tr><td>-0.0312±3.7956i</td><td>0.6041</td><td>8375</td><td></td><td></td><td></td></tr><tr><td>-0.0151±3.3252i</td><td>0.5292</td><td>8875</td><td></td><td></td><td></td></tr><tr><td>-0.0348±4.4752i</td><td>0.7123</td><td>8875</td><td></td><td></td><td></td></tr></tbody></table>
<table><caption>Table 7
Computations of first Hopf bifurcation point—critical numbers for several mesh size</caption><thead><tr><th>Numerical tests</th><th>mesh</th><th>Re<sub>c</sub></th><th>S<sub>t</sub></th></tr></thead><tbody><tr><td rowspan="3">Flow around a cylinder</td><td>2800 d.o.f.</td><td>52.716</td><td>0.167</td></tr><tr><td>5000 d.o.f.</td><td>49.33</td><td>0.171</td></tr><tr><td>6400 d.o.f.</td><td>45.61</td><td>0.163</td></tr><tr><td rowspan="5">2-D lid driven cavity</td><td>3362 d.o.f., reg.-mesh</td><td>9224</td><td>0.19</td></tr><tr><td>3362 d.o.f., irr.-mesh</td><td>9151</td><td>0.48</td></tr><tr><td>7442 d.o.f., irr.-mesh</td><td>7119</td><td>0.45</td></tr><tr><td>13 122 d.o.f., reg.-mesh</td><td>7673</td><td>0.45</td></tr><tr><td>13 122 d.o.f., irr.-mesh</td><td>7890</td><td>0.44</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>Fig. 14. Meshes used for the 2-D lid-driven cavity.</figcaption></figure>

### 4.2.5. Influence of the grid size and discretization scheme

We will now present a short examination of the influence of the mesh size for characterizing the Hopf bifurcation. This study concerns the flow around a cylinder and the 2-D lid-driven cavity. For the flow around a cylinder, we use three different meshes and the corresponding numbers of d.o.f. are given in Table 7. Five meshes are used for the 2-D lid-driven cavity. One can distinguish two kinds of mesh: regular mesh (Fig. 14(a), denoted by reg.-mesh in Table 7) and irregular mesh (Fig. 14(b), denoted by irr.-mesh in Table 7). Clearly the critical Reynolds number and the Strouhal number converge with the fineness of the mesh. Our final results are about $Re_c = 45.5–46$, $St = 0.16–0.17$ for the flow around the cylinder, $Re_c = 7600–7900$, $St = 0.44–0.45$ for the driven cavity. As noted by Fortin et al. [4], the critical frequency is not highly mesh sensitive. Moreover, in the second example, the accuracy of the results seems to be more depending on the total number of d.o.f., rather than on a local fineness near the wall. Nevertheless, meshes used in this study seem to be not enough fine to give relevant conclusions.

The weighting method is another important key for finite element solution of the Navier-Stokes equations. For example, a Petrov-Galerkin [18,30] weighting can be chosen. The advantage of this formulation is the elimination of the spatial oscillations that occur with a Galerkin weighting [18]. The difference between these formulations lies in the discretization of the convection term and the mass matrix [18,30]. Note that all previous results have been obtained with a Galerkin method. Details concerning the computation of the fundamental solutions and the bifurcating points with the Asymptotic Numerical Method and a Petrov-Galerkin weighting are presented in Refs. [13,17]. In Ref. [17], it has been stated that the Petrov-Galerkin formulation often leads to incorrect estimations of the Hopf bifurcation. This conclusion is consistent with the point of view of Gresho [31] and Sahin and Owens [29], who do not recommend the use of upwind methods (such as the Petrov-Galerkin method) for computing fluid flow and instabilities.
## 5. Conclusions

According to the computational tests presented here, the numerical methods we propose in this paper give critical values that are consistent with the results in the literature, for both stationary bifurcations and Hopf bifurcations.

In the case of the stationary indicator, the numerical method presented here can be improved using asymptotic expansions (see Boutyour et al. [8] and Tri et al. [15]) and an automatic algorithm could be used to switch on the bifurcating branches (see Vannucci et al. [16]).

In the case of Hopf bifurcation, the method gives accurate values of the critical parameters. Moreover, this method enables ready computation of the complex bifurcating eigenvector. Note that with a direct method [7] the complex bifurcating eigenvector is also determined and the critical parameters are also very accurate. With the algorithm presented in this paper, several critical points can be determined (see Section 4.2.3 for the flow in a channel with a symmetric expansion). With this method, the convergence to a Hopf bifurcation point is not secure and depends on the choice of the initial values. However, the detection of a Hopf bifurcation point is difficult to automate with the proposed method, because the indicator does not change its sign at the critical point. Thus, the method presented in this paper could be combined with a direct method, for example by providing a good initial guess for the exact computation.

We have not discussed the CPU time needed with the proposed algorithm in this paper. In fact the largest amount of time is for the triangulation of the operator $\mathcal{E}_T(U_\lambda, \omega_0)$. This is nearly ten times that required to triangulate the tangent matrix $L_T$ of the stationary nonlinear problem due to an appropriate storage of the operator $\mathcal{E}_T(U_\lambda, \omega_0)$. On the other hand, the total CPU time is governed by the number of steps in the continuation method (this means the number of triangulations of the operator $\mathcal{E}_T(U_\lambda, \omega_0)$). This number of steps can be reduced by using a continuation technique based on Padé approximants [32] instead of polynomial approximations. This continuation method was recently proposed by Elhage-Hussein et al. [33] and leads to a reduction of nearly 50 per cent in the number of steps within the asymptotic numerical method. Such a method can easily be applied to the algorithm presented here and the same reduction in the number of steps in the continuation is expected.

## References

[1] J.E. Marsden, M. McCracken, The Hopf Bifurcation and its Applications, Springer, New York, 1976.

[2] D. Roosevelt, An algorithm for the computation of Hopf bifurcation points in comparison with other methods, J. Comput. Appl. Math. 12–13 (1985) 517–529.

[3] B. Hassard, N. Kazarinoff, Y.-H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981.

[4] A. Fortin, M. Jardak, J.J. Gervais, R. Pierre, Localization of Hopf bifurcations in fluid flow problems, Int. J. Numer. Methods Fluids 24 (1997) 1185–1210.

[5] N. Nayar, J.M. Ortega, Computation of selected eigenvalues of generalized eigenvalue problems, J. Comput. Phys. 108 (1993) 8–14.

[6] A.D. Jepson, Numerical Hopf bifurcation, Thesis, California Institute of Technology, 1981.

[7] C.P. Jackson, A finite-element study of the onset of vortex shedding in flow past variously shaped body, J. Fluid Mech. 182 (1987) 23–45.

[8] E.H. Boutyour, B. Cochelin, M. Potier-Ferry, Calculs des points de bifurcation par une méthode asymptotique numérique, in: Proceedings 1st Congrès National de Mécanique au Maroc, 1993, pp. 371–378.

[9] M.E.H. Bensaadi, Méthode asymptotique-numérique pour le calcul de bifurcations de Hopf et de solutions périodiques, Thesis, Université de Metz, France, 1995.

[10] A. Tri, Méthodes asymptotiques-numériques pour les fluides visqueux incompressibles et la détection de la bifurcation de Hopf, Thesis, Université de Metz, France, 1996.

[11] B. Cochelin, A path-following technique via an asymptotic-numerical method, Comput. Struct. 53 (5) (1994) 1181–1192.

[12] N. Damil, M. Potier-Ferry, A new method to compute perturbed bifurcations: application to the buckling of imperfect elastic structures, Int. J. Engrg. Sci. 28 (9) (1990) 943–957.

[13] J.M. Cadou, B. Cochelin, N. Damil, M. Potier-Ferry, Asymptotic numerical method for stationary Navier–Stokes equations and with Petrov–Galerkin formulation, Int. J. Numer. Methods Engrg. 50 (2001) 825–845.

[14] H. Lahmam, J.M. Cadou, H. Zahrouni, N. Damil, M. Potier-Ferry, High-order predictor-corrector algorithms, Int. J. Numer. Methods Engrg. 55 (2002) 685–704.

[15] A. Tri, B. Cochelin, M. Potier-Ferry, Résolution des équations de Navier–Stokes et détection des bifurcations stationnaires par une méthode asymptotique numérique, Revue Européenne des éléments finis 5 (1996) 415–442.

[16] P. Vannucci, B. Cochelin, N. Damil, M. Potier-Ferry, An asymptotic-numerical method to compute bifurcating branches, Int. J. Numer. Methods Engrg. 41 (1998) 1365–1389.

[17] J.M. Cadou, Méthode Asymptotique-Numérique pour le calcul des branches solutions et des instabilités dans les fluides et pour les problèmes d'interaction fluide-structure, Thesis, Université de Metz, France, 1997.

[18] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, vol. 2, fourth ed., McGraw-Hill Book Company, 1991 (Section 12.13).
[19] J.M.T. Thompson, A.C. Walker, The non-linear perturbative analysis of discrete structural systems, Int. J. Solids Struct. 4 (1968) 757-768.

[20] M. Shapira, D. Degani, D. Weihs, Stability and existence of multiple solutions for viscous flow in suddenly enlarged channels, Computers & Fluids 18 (1990) 239-258.

[21] C. Allery, J.M. Cadou, A. Hamdouni, D. Razafindralandy, Application of the asymptotic numerical method to the Coanda effect study, Revue Européenne des éléments finis 13 (1-2) (2004) 57-77.

[22] K.M. Kelkar, S.V. Patankar, Numerical prediction of vortex shedding behind a square cylinder, Int. J. Numer. Methods Fluids 14 (1992) 8-14.

[23] A. Okajima, Strouhal numbers of rectangular cylinders, J. Fluids Mech. 123 (1982) 379-398.

[24] W. Cazemier, Proper orthogonal decomposition and low-dimensional models for turbulent flows, Thesis, University of Groningen, The Netherlands, 1997.

[25] Y.F. Peng, Y.H. Shiau, R.R. Hwang, Transition in a 2-d lid-driven cavity flow, Computers & Fluids 32 (2003) 337-352.

[26] M. Poliashenko, C.K. Aidun, A direct method for computation of simple bifurcations, J. Comput. Phys. 121 (1995) 246-260.

[27] W. Cazemier, R.W.C.P. Verstappen, A.E.P. Veldam, Proper orthogonal decomposition and low-dimensional models for driven cavity flows, Phys. Fluids 10 (7) (1998) 1685-1699.

[28] G. Tiensenga, F.W. Wubs, A.E.P. Veldam, Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method, J. Comput. Appl. Math. 140 (2002) 751-772.

[29] M. Sahin, R.G. Owens, A novel fully-implicit finite volume method applied to lid-driven cavity problem. Part II. Linear stability analysis, Int. J. Numer. Methods Fluids 42 (2003) 79-88.

[30] T.J.R. Hughes, A. Brooks, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199-259.

[31] P.M. Gresho, R.L. Lee, Don't suppress the wiggles-they're telling you something, Computers & Fluid 9 (1981) 223-253.

[32] G.A. Baker, P. Graves-Morris, Padé Approximants, Part I: Basic Theory, Encyclopedia of Math. Appl., vol. 13, Addison-Wesley, Reading, MA, 1981.

[33] A. Elhage-Hussein, M. Potier-Ferry, N. Damil, A numerical continuation method based on padé approximants, Int. J. Solids Struct. 37 (2000) 6981-7001.