2.1. Steady-state solution

At small Reynolds numbers, the flow is steady and symmetric with respect to the $x$ -axis, though it becomes unstable at larger Reynolds numbers than a critical value. The solution corresponding to the steady symmetric flow, say $(\overline{{\it\psi}},\overline{{\it\omega}})$ , satisfies basic (2.2) and (2.3) under boundary conditions (2.5)–(2.8) even for larger Reynolds numbers above the critical value, which we call the main flow and investigate its linear and weakly nonlinear stability to elucidate the nonlinear interaction between the two wakes. The stream function and the vorticity $(\overline{{\it\psi}},\overline{{\it\omega}})$ are odd functions of $y$ , i.e.  $\overline{{\it\psi}}(x,-y)=-\overline{{\it\psi}}(x,y)$ and $\overline{{\it\omega}}(x,-y)=-\overline{{\it\omega}}(x,y)$ , and they are obtained numerically by solving the steady-state vorticity transport equation:

together with the Poisson equation:

under the same boundary conditions for $({\it\psi},{\it\omega})$ including (2.5)–(2.8) and (2.10).

2.2. Linear and weakly nonlinear stability analyses

In order to investigate the stability of the steady symmetric flow, we consider a disturbance $(\widehat{{\it\psi}},\widehat{{\it\omega}})$ added to the main flow $(\overline{{\it\psi}},\overline{{\it\omega}})$ and express the vorticity and the stream function as

respectively. Substituting these expressions into (2.2) and subtracting (2.11) from the resultant equation, we obtain the disturbance equation as

Similarly, the Poisson equation for $\widehat{{\it\psi}}$ is obtained as

The boundary conditions for $(\widehat{{\it\psi}},\widehat{{\it\omega}})$ on the upstream and side boundaries are given as

The outlet condition on the downstream boundary (CD) is taken as the Sommerfeld radiation condition

where $c$ is taken to be $\partial \overline{{\it\psi}}/\partial y$ at each point on CD. The non-slip boundary condition is taken on the surfaces of the square cylinders, which is written as

On the cylinder surfaces, the values of $\widehat{{\it\psi}}\ (=\widehat{{\it\psi}}_{\text{B}1},\ \widehat{{\it\psi}}_{\text{B}2})$ are determined in the same way as in the numerical simulations.

In the linear stability analysis, we neglect the nonlinear term of the disturbance $(\widehat{{\it\psi}},\widehat{{\it\omega}})$ in (2.14), and assume the time dependence of the disturbance as $\widehat{{\it\psi}}=\widetilde{{\it\psi}}(x,y)\text{e}^{{\it\lambda}t}$ and $\widehat{{\it\omega}}=\widetilde{{\it\omega}}(x,y)\text{e}^{{\it\lambda}t}$ . Then we obtain a linearized disturbance equation as

which constitutes the eigenvalue problem for the linear stability analysis together with the Poisson equation

under the appropriate boundary conditions. The boundary conditions for $(\widetilde{{\it\psi}},\widetilde{{\it\omega}})$ are the same as for $(\widehat{{\it\psi}},\widehat{{\it\omega}})$ except that the condition

is used instead of (2.17). The solution of the eigenvalue problem gives an infinite number of eigenvalues, ${\it\lambda}$ and eigenfunctions ( $\widetilde{{\it\psi}},\widetilde{{\it\omega}}$ ), which are classified into symmetric and antisymmetric modes, $(\widetilde{{\it\psi}}_{s},\widetilde{{\it\omega}}_{s})$ and $(\widetilde{{\it\psi}}_{a},\widetilde{{\it\omega}}_{a})$ , according to the symmetry of the corresponding flow field with respect to the $x$ -axis. The eigenfunctions for the symmetric and antisymmetric modes are odd and even functions, respectively, so that the relations $\widetilde{{\it\psi}}_{s}(x,-y)=-\widetilde{{\it\psi}}_{s}(x,y)$ , $\widetilde{{\it\omega}}_{s}(x,-y)=-\widetilde{{\it\omega}}_{s}(x,y)$ and ( $\widetilde{{\it\psi}}_{a}(x,-y)=\widetilde{{\it\psi}}_{a}(x,y)$ $\widetilde{{\it\omega}}_{a}(x,-y)=\widetilde{{\it\omega}}_{a}(x,y)$ ) hold. The eigenvalues, ${\it\lambda}$ ( $={\it\sigma}-\text{i}{\it\mu}$ ), are complex in general whose real and imaginary parts, ${\it\sigma}$ and ${\it\mu}$ , show the growth rate and the frequency (angular frequency) of the disturbance, respectively. Among infinite number of eigenmodes, the mode having the largest growth rate for each Reynolds number and gap ratio is called the most unstable mode. The Reynolds numbers at which the growth rate ${\it\sigma}$ of each mode vanishes are called neutral Reynolds numbers, above which the mode of disturbance grows, otherwise it decays, and the smallest neutral Reynolds number gives the critical state so that it is called a critical Reynolds number. As will be shown later, the antisymmetric mode determines the critical Reynolds number for smaller gap ratios than a threshold value, and the symmetric mode for larger gap ratios. There are two kinds of instability, one of which is the case of ${\it\mu}=0$ at the criticality where ${\it\sigma}$ vanishes, and in the other case ${\it\mu}\neq 0$ as in the present flow. In the former case, the growing disturbance is stationary and as a result of growth of stationary disturbance, another steady flow emerges from the main flow $(\overline{{\it\psi}},\overline{{\it\omega}})$ , i.e. the number of steady solution changes, which is called a pitchfork bifurcation when $(\overline{{\it\psi}},\overline{{\it\omega}})$ has a symmetry and the symmetry is broken due to the bifurcation. While in the latter case, the growing disturbance is oscillatory and the resultant flow due to the instability is oscillatory flow with the superposition of the oscillatory disturbance on the main flow. The emergence of oscillatory flow from the main flow is called a Hopf bifurcation in dynamical systems theory (see Drazin Reference Drazin1992, §1.6).

In the case where the symmetric steady flow is unstable only to one of the symmetric or antisymmetric modes of disturbance, the disturbance is expected to tend to an equilibrium state via nonlinearity of itself. On the other hand, if both modes are growing, they make nonlinear interactions with each other as well as with their own to lead to an equilibrium state. There are three possibilities for the equilibrium state, in which only one of the two modes survives overwhelming the other, or both modes have finite amplitudes. In considering the symmetries of the solution resulting from the nonlinear interactions, it is helpful to recall that ${\it\psi}(x,y,t)$ has the same symmetry both in $x$ and $y$ with ${\it\omega}(x,y,t)$ because the operator ${\it\Delta}$ does not change these symmetries. However, the nonlinear term $N({\it\psi},\ {\it\omega})$ in (2.2) yields a component of an odd function of $y$ if both ${\it\psi}(x,y,t)$ and ${\it\omega}(x,y,t)$ are even functions of $y$ , whereas it keeps the parity of function if both of them are odd. Note that, although it is confusing enough, the functions $\widetilde{{\it\psi}}_{s}$ and $\widetilde{{\it\omega}}_{s}$ for the symmetric mode are odd functions of $y$ , and their flow fields are symmetric with respect to the $x$ -axis. Since the symmetric mode of disturbance $\widetilde{{\it\psi}}_{s}$ is an odd function of $y$ so that $\widetilde{{\it\psi}}_{s}$ induces only symmetric higher harmonics in nonlinear states after interacting with the steady symmetric flow as well as with itself. However, as the antisymmetric mode is an even function of $y$ , it induces symmetric component as the first harmonic mode through nonlinear interaction of itself. Hence, in the case where only the symmetric mode survives the resultant equilibrium state displays a symmetric flow pattern, being superposed on the steady symmetric flow. On the other hand, even if only the antisymmetric mode survives the flow field is not antisymmetric, but asymmetric with respect to the $x$ -axis because it is a superposition of the antisymmetric disturbance on the steady symmetric flow. This flow property makes it difficult to evaluate the amplitude of the antisymmetric mode of disturbance in the nonlinear state from numerical simulation data. In order to resolve this difficulty, we need to appropriately define the amplitudes of the two modes.

We use Watson’s amplitude expansion method to derive the dynamical equations for the amplitudes of the two modes of disturbance, $A_{s}$ and $A_{a}$ , together with Herbert’s method to define the amplitudes of disturbance uniquely, taking consideration into the symmetry of the first harmonic component. We choose two points $\text{P}_{1}$ ( $(x,y)=(20,({\it\Gamma}+1)/2)$ ) and $\text{P}_{2}\ ((x,y)=(20,-({\it\Gamma}+1)/2))$ shown in figure 1 in defining the real-valued amplitudes of the antisymmetric and symmetric components of disturbance as $A_{s}=-(\partial \widehat{{\it\psi}}_{s}/\partial x)(\text{P}_{1})$ and $A_{a}=-(\partial \widehat{{\it\psi}}_{a}/\partial x)(\text{P}_{1})$ , following Herbert (Reference Herbert1983). With these definitions, the $y$ -components of velocity $v_{1}$ and $v_{2}$ at $\text{P}_{1}$ and $\text{P}_{2}$ are expressed as $v_{1}\!=-(\partial \widehat{{\it\psi}}_{s}/\partial x)(\text{P}_{1})\!-(\partial \widehat{{\it\psi}}_{a}/\partial x)(\text{P}_{1})$ and $v_{2}=(\partial \widehat{{\it\psi}}_{s}/\partial x)(\text{P}_{2})-(\partial \widehat{{\it\psi}}_{a}/\partial x)(\text{P}_{2})$ as long as we consider oscillation having only Fourier modes with the frequencies ${\it\mu}_{s}$ and ${\it\mu}_{a}$ .

We expand $\widehat{{\it\psi}}_{s}$ and $\widehat{{\it\psi}}_{a}$ in powers of $A_{s}$ and $A_{a}$ as well as in a Fourier series as

The expansion coefficient functions satisfy the relations

because of the reality of $\widehat{{\it\psi}}$ , where the asterisk denotes the complex conjugate hereafter. We also expand nonlinear growth rates ${\it\lambda}_{s}$ and ${\it\lambda}_{a}$ in powers of $A_{s}$ and $A_{a}$ , which gives

As the complex growth rates are denoted as ${\it\lambda}_{s}={\it\sigma}_{s}-\text{i}{\it\mu}_{s}$ and ${\it\lambda}_{a}={\it\sigma}_{a}-\text{i}{\it\mu}_{a}$ , equations (2.24) are equivalent to the expansions,

Note that the real parts of the complex growth rates, ${\it\sigma}_{s}$ and ${\it\sigma}_{a}$ , are expressed in terms of $A_{s}$ and $A_{a}$ as

respectively.

The basic equations governing the time development of the disturbance, $\widehat{{\it\psi}}$ and $\widehat{{\it\omega}}$ , are (2.14) and (2.15), which we express symbolically in a simple form by using only $\widehat{{\it\psi}}$ as

where $\mathscr{L}\widehat{{\it\psi}}$ and $\mathscr{N}(\widehat{{\it\psi}}_{1},\widehat{{\it\psi}}_{2})$ are defined as

Substituting (2.22), (2.24) and (2.26) into (2.27), and equating the same Fourier component and the same power of the amplitude like $A_{s}^{2m+n}A_{a}^{2\ell }\text{e}^{-\text{i}n{\it\mu}_{s}t}$ or $A_{a}^{2m+n}A_{s}^{2\ell }\text{e}^{-\text{i}n{\it\mu}_{a}t}$ , we obtain a set of equations for the expansion coefficient functions ${\it\phi}_{snm\ell }$ and ${\it\phi}_{anm\ell }$ . Among the resultant infinite sets of equations, the lowest-order equations of $O(A_{s})$ or $O(A_{a})$ are nothing but the linear stability (2.19) and (2.20), which are expressed simply as

Since (2.29) for ${\it\phi}_{s100}$ and ${\it\phi}_{a100}$ are linear and homogeneous, and their boundary conditions are as well, we normalize them to be $-(\partial {\it\phi}_{s100}/\partial x)(\text{P}_{1})=1$ and $-(\partial {\it\phi}_{a100}/\partial x)(\text{P}_{1})=1$ , being consistent with the definitions of $A_{s}$ and $A_{a}$ .

At the second order of $O(A_{s}^{2})$ or $O(A_{a}^{2})$ , we obtain the equations for the mean flow distortions ( ${\it\phi}_{s010}$ , ${\it\phi}_{a010}$ ) and the first harmonic modes ( ${\it\phi}_{s200}$ , ${\it\phi}_{a200}$ ) as

and

respectively. As (2.30) and (2.31) are also linear, but inhomogeneous, we can solve these equations under the homogeneous boundary conditions similar to (2.16)–(2.18).

At the third order of $O(A_{s}^{3})$ , $O(A_{a}^{3})$ , $O(A_{s}A_{a}^{2})$ or $O(A_{s}^{2}A_{a})$ , we obtain the equations for ${\it\phi}_{s110}$ , ${\it\phi}_{a110}$ , ${\it\phi}_{s101}$ and ${\it\phi}_{a101}$ as

Equations (2.32) for ${\it\phi}_{s110}$ , ${\it\phi}_{a110}$ , ${\it\phi}_{s101}$ and ${\it\phi}_{a101}$ include unknown constants ${\it\lambda}_{s10}$ , ${\it\lambda}_{a10}$ , ${\it\lambda}_{s01}$ and ${\it\lambda}_{a01}$ , so that we need some technique in solving these equations. We will follow Watson (Reference Watson1960) to solve the equations and express the solutions as

Substituting (2.33) into (2.32) and equating the terms including the factors ${\it\lambda}_{s10}$ , ${\it\lambda}_{a10}$ , ${\it\lambda}_{s01}$ and ${\it\lambda}_{a01}$ , we obtain the equations for ${\it\chi}_{s010}$ , ${\it\chi}_{a010}$ , ${\it\chi}_{s001}$ and ${\it\chi}_{a001}$ , which are easily solved as

On the other hand, the equations for ${\it\chi}_{s110}$ , ${\it\chi}_{a110}$ , ${\it\chi}_{s101}$ and ${\it\chi}_{a101}$ consist of terms which are left in (2.32) after subtracting the terms including the factors ${\it\lambda}_{s10}$ , ${\it\lambda}_{a10}$ , ${\it\lambda}_{s01}$ and  ${\it\lambda}_{a01}$ , and are expressed as

The functions ${\it\chi}_{s1m\ell }$ and ${\it\chi}_{a1m\ell }$ with $\ell +m>0$ are obtained without any arbitrary constants by solving inhomogeneous equations (2.36) under the boundary conditions similar to (2.16)–(2.18).

The expansion coefficients ${\it\lambda}_{s10}$ , ${\it\lambda}_{a10}$ , ${\it\lambda}_{s01}$ and ${\it\lambda}_{a01}$ of ${\it\lambda}_{s}$ and ${\it\lambda}_{a}$ in (2.24) are expressed in terms of ${\it\chi}_{s1m\ell }$ and ${\it\chi}_{a1m\ell }\ (\ell +m>0)$ as

where we used the relations $(\partial {\it\phi}_{s1m\ell }/\partial x)(\text{P}_{1})=0$ and $(\partial {\it\phi}_{a1m\ell }/\partial x)(\text{P}_{1})=0$ , for ( $\ell +m>0$ ) which come from the unique definitions of $A_{s}$ and $A_{a}$ at $\text{P}_{1}$ . Thus, we have determined all of the expansion coefficients in (2.24). Defining two complex amplitudes $\widetilde{A}_{s}$ and $\widetilde{A}_{a}$ by $\widetilde{A}_{s}=A_{s}\text{e}^{-\text{i}{\it\mu}_{s}t}$ and $\widetilde{A}_{a}=A_{a}\text{e}^{-\text{i}{\it\mu}_{a}t}$ , we finally obtain a set of amplitude equations for the complex amplitudes $\widetilde{A}_{s}$ and $\widetilde{A}_{a}$ as

3.1. Numerical methods and accuracy assessments

In numerical simulations, numerical domain ABCD is discretized with a non-uniform rectangular mesh having variable differences ${\rm\Delta}x$ and ${\rm\Delta}y$ in the $x$ and $y$ directions. The mesh sizes ${\rm\Delta}x$ and ${\rm\Delta}y$ are smallest near the cylinder surface, and become larger going away from it, finally reaching constant values beyond prescribed distances. The spatial derivatives are approximated with the central finite difference of fourth-order accuracy. The vorticity transport equation (2.2) is solved under the appropriate boundary and initial conditions by approximating the time derivative with the fourth-order Runge–Kutta method taking a finite difference ${\rm\Delta}t$ . The Poisson equation (2.3) for ${\it\psi}$ is solved under the appropriate boundary conditions with the SOR (successive over-relaxation) method.

The extent of the numerical domain is crucial to assure the accuracy of numerical results. In order to study the effect of the extent of the numerical domain, we obtain maximum ( $v_{1max}$ , $v_{2max}$ ) and minimum ( $v_{1min}$ , $v_{2min}$ ) values and amplitudes ( $v_{1max}-v_{1min}$ , $v_{2max}-v_{2min}$ ) of oscillation in the velocities $v_{1}$ and $v_{2}$ at $\text{P}_{1}$ ( $(x,y)=(20,{\it\Gamma}/2+1/2)$ ) and $\text{P}_{2}$ ( $(x,y)=(20,-{\it\Gamma}/2-1/2)$ ). We compare these values between numerical domain I ( $L_{1}=9.96$ , $L_{2}=140.1$ , $W=40$ ) and domain II ( $L_{1}=19.4$ , $L_{2}=280.2$ , $W=80$ ) for $\mathit{Re}=60$ and ${\it\Gamma}=2.0$ in table 1. The relative errors are smaller than $0.6\,\%$ . Hence, we adopt numerical domain I in our numerical simulations.

Table 1.

Numerical domain assessment, $\mathit{Re}=60$ , ${\it\Gamma}=2.0$ . Domain I: $L_{1}=9.96$ , $L_{2}=140.1$ , $W=40$ . Domain II: $L_{1}=19.4$ , $L_{2}=280.2$ , $W=80$ .

We carried out assessments of appropriateness of the number of mesh points, where we compare the values ( $v_{1max}$ , $v_{2max}$ , $v_{1min}$ , $v_{2min}$ , $v_{1max}-v_{1min}$ , $v_{2max}-v_{2min}$ ) by taking two different numbers of mesh points, i.e. I ( $1569\times 201$ ) and II ( $3137\times 201$ ). These values are tabulated in table 2. At $\mathit{Re}=60$ and ${\it\Gamma}=2.0$ , the relative errors between I and II are smaller than $0.4\,\%$ , so that we adopt I ( $1569\times 201$ ) as the number of mesh points in our numerical simulations, in which 20 mesh points are taken along each side of the square cylinders.

Table 2.

Assessment of the number of mesh points, $\mathit{Re}=60$ , ${\it\Gamma}=2.0$ .

The magnitude of finite difference ${\rm\Delta}t$ in time discretization is also assessed by evaluating the oscillation amplitudes $v_{1max}-v_{1min}$ and $v_{2max}-v_{2min}$ for two finite differences ${\rm\Delta}t=0.001$ and 0.002. The relative errors between the two cases are smaller than $0.3\,\%$ for $\mathit{Re}=60$ and ${\it\Gamma}=2.0$ as shown in table 3, and we adopt ${\rm\Delta}t=0.002$ .

Table 3.

Assessment of the finite differences in the temporal discretization, ${\rm\Delta}t$ . Discretization I: ${\rm\Delta}t=0.002$ . Discretization II: ${\rm\Delta}t=0.001$ .

3.2. Linear and weakly nonlinear stability analyses

In numerical calculations of the steady symmetric flow and linear and weakly nonlinear stability analyses, we use the same non-uniform rectangular mesh adopted in numerical simulation, and approximate the spatial derivatives with the central finite difference of fourth-order accuracy. The eigenvalue problem composed of the linear disturbance equation (2.29) and the appropriate boundary conditions is solved with the double QR method. Inhomogeneous equations, (2.30), (2.31) and (2.36), are solved under the appropriate boundary conditions with the SOR iteration method, where convergence is determined when the maximum differences between two consecutive iteration values of ${\it\psi}$ and ${\it\omega}$ over the entire numerical domain becomes less than $1\times 10^{-7}$ .

4.1. Linear stability analysis of the steady symmetric flow

We have obtained the steady symmetric flow numerically by solving (2.11) and (2.12) under the boundary conditions (2.5)–(2.8) and (2.10), and analysed its linear stability by solving the eigenvalue problem consisting of (2.19), (2.20) and the boundary conditions for $\widetilde{{\it\psi}}$ and $\widetilde{{\it\omega}}$ for each Reynolds number and each gap ratio in the range of $1.0\leqslant {\it\Gamma}\leqslant 3.5$ . The complex linear growth ratio ${\it\lambda}\ (={\it\sigma}-\text{i}{\it\mu})$ thus obtained is depicted in figure 2(a) for ${\it\Gamma}=2.0$ , for instance, where the solid and broken lines indicate the real part ${\it\sigma}$ and imaginary part ${\it\mu}$ , respectively. The antisymmetric mode of disturbance has a larger linear growth rate, say ${\it\sigma}_{a00}$ , than the symmetric mode ( ${\it\sigma}_{s00}$ ) throughout the range of Reynolds number depicted in this figure. The critical and neutral Reynolds numbers are evaluated as $\mathit{Re}_{a}=48.17$ and $\mathit{Re}_{s}=54.41$ for the antisymmetric and symmetric modes, respectively, and hence the flow becomes unstable to the antisymmetric mode first for ${\it\Gamma}=2.0$ . The angular velocities, say ${\it\mu}_{s00}$ and ${\it\mu}_{a00}$ , increase very slowly with $\mathit{Re}$ being almost constant, and the symmetric mode has a larger value of ${\it\mu}_{s00}$ than the antisymmetric mode ( ${\it\mu}_{a00}$ ). The Strouhal numbers $\mathit{St}$ are evaluated for both modes as $\mathit{St}_{s}={\it\mu}_{s00}/(2{\rm\pi})=0.1584$ , $\mathit{St}_{a}={\it\mu}_{a00}/(2{\rm\pi})=0.1418$ at the criticality, respectively. On the other hand, in the case of ${\it\Gamma}=2.5$ the symmetric mode gives the critical condition as seen in figure 2(b), having a larger linear growth rate ${\it\sigma}_{s00}$ than the antisymmetric mode ( ${\it\sigma}_{a00}$ ). The critical and neutral Reynolds numbers are $\mathit{Re}_{s}=50.31$ and $\mathit{Re}_{a}=52.19$ , and the Strouhal numbers are $\mathit{St}_{s}=0.1640$ , $\mathit{St}_{a}=0.1490$ at each critical state.

Figure 2.

Linear growth rate ${\it\lambda}$ : (a) ${\it\Gamma}=2.0$ ; (b) ${\it\Gamma}=2.5$ . Solid line: real part ${\it\sigma}$ . Broken line: imaginary part ${\it\mu}$ .

As no appropriate experimental or numerical data are found to compare with our linear stability results at similar Reynolds numbers, we compare them with numerical simulation data performed at rather larger Reynolds number of $\mathit{Re}=100$ by Kang (Reference Kang2003), where Strouhal numbers $\mathit{St}_{s}=0.177$ and $\mathit{St}_{a}=0.168$ are reported for circular cylinder case with ${\it\Gamma}=2.0$ . Our values of $\mathit{St}_{s}$ and $\mathit{St}_{a}$ are a little smaller than Kang’s numerical values, because of the smaller value of the Reynolds number at the criticality in the present analysis and the difference of the cylinder shape.

We show the steady symmetric flow and the flow fields of the eigenfunctions corresponding to the two most unstable modes in the linear stability analysis in figure 3 for ${\it\Gamma}=2.0$ and $\mathit{Re}=60$ as typical examples. The flow pattern of the steady symmetric flow is plain as shown in figure 3(a), which almost keeps its shape even for different Reynolds numbers, with only the length of the attached vortex being longer for larger Reynolds numbers. The real part of the eigenfunction $\widetilde{{\it\psi}}_{s}$ for the symmetric mode of disturbance is an odd function of $y$ . The flow field of $\widetilde{{\it\psi}}_{s}$ is depicted in figure 3(b) by contour lines, which consists of two lines of vortices, where the dark and light lines indicate contours of positive and negative values of $\widetilde{{\it\psi}}_{s}$ , respectively. Each couple of adjacent and facing vortices have opposite senses of rotation. This mode of disturbance yields a symmetric flow field characterized by double vortex streets. Note that the contour lines of $\widetilde{{\it\psi}}_{s}$ extend over the whole range of this figure up to $x=100$ , while the eigenfunction $\widetilde{{\it\psi}}_{a}$ for the antisymmetric mode diminishes in magnitude in a shorter distance than $x\sim 50$ as shown in figure 3(c). The rapid decay of the eigenfunctions is responsible to early annihilation of the two vortex streets antisymmetrically arranged behind the cylinders, and the asymmetric flow induced by the antisymmetric mode of disturbance loses its initial vortex configuration after advected only in a short distance behind the cylinders, as will be explained in § 5. The imaginary parts of the eigenfunctions are similar to the corresponding real parts except phase differences, so that they are omitted to be shown to save space.

Figure 3.

Steady symmetric flow and the eigenfunctions corresponding to the two most unstable modes in the linear stability analysis (contour lines of stream function), $\mathit{Re}=60$ , ${\it\Gamma}=2.0$ : (a) steady symmetric flow (main flow); (b) symmetric mode of disturbance; (c) antisymmetric mode of disturbance.

The neutral Reynolds numbers, at which each of the growth rates ${\it\sigma}_{s00}$ and ${\it\sigma}_{s00}$ of symmetric and antisymmetric modes of disturbance vanishes, are shown in figure 4, where the solid and broken lines indicate the symmetric and antisymmetric modes, respectively. The two lines intersect at a point $({\it\Gamma}_{c},Re)=(2.25,51.3)$ . We call the gap ratio ${\it\Gamma}_{c}=2.25$ the threshold gap ratio. For smaller gap ratios than this value ( ${\it\Gamma}<{\it\Gamma}_{c}$ ), an antisymmetric mode of disturbance is the most unstable mode, while for ${\it\Gamma}>{\it\Gamma}_{c}$ the most unstable mode is a symmetric mode. The parameter domain of ( ${\it\Gamma}$ , $\mathit{Re}$ ) is classified by the symmetry of the growing mode(s) of disturbance into four regions I, II, III and IV as shown in figure 4. In region I, the steady symmetric flow is stable to both modes of disturbance ( ${\it\sigma}_{s00}<0$ , ${\it\sigma}_{a00}<0$ ). Regions II and III shows that the steady symmetric flow is unstable only to the antisymmetric ( ${\it\sigma}_{a00}>0$ , ${\it\sigma}_{s00}<0$ ) or the symmetric mode ( ${\it\sigma}_{s00}>0$ , ${\it\sigma}_{a00}<0$ ), respectively. Region IV ( ${\it\sigma}_{s00}>0$ , ${\it\sigma}_{a00}>0$ ) is the parameter region that we are now concerned with, in which both modes are growing. When the both modes grow, they compete with each other through nonlinear interactions between them. As a result of the competition, one of the symmetric and antisymmetric modes of disturbance survive overwhelming the other, or both modes remains having finite amplitudes. The resultant equilibrium state of the total flow is established as a superposition of the surviving mode(s) of disturbance on the steady symmetric flow.

Figure 4.

Critical Reynolds number as a function of the gap ratio: solid line, symmetric mode; broken line, antisymmetric mode; I, steady symmetric flow; II, unstable to the antisymmetric mode of disturbance; III, unstable to the symmetric mode of disturbance; IV, unstable to the both antisymmetric and symmetric modes of disturbance.

It is noted here that the linear stability of the flow behind two circular cylinders shows similar behaviour to the present case of square cylinders, although differences between the flows around square and circular cylinders are often overemphasized. The critical gap ratio is ${\it\Gamma}_{c}=2.34$ at $\mathit{Re}_{c}=45.2$ for the circular cylinders, across which the most unstable mode changes from the antisymmetric to the symmetric mode similarly to the square cylinder case as shown by Akinaga & Mizushima (Reference Akinaga and Mizushima2005).

4.2. Weakly nonlinear stability analysis

We have evaluated all of the coefficients included in the coupled amplitude equations (2.38) by numerically solving (2.30), (2.31) and (2.36) under the appropriate boundary conditions for parameter values in region IV. Among numerous sets of numerical results, we show only two cases of ${\it\Gamma}=2.0$ and 2.5 at the same Reynolds number $\mathit{Re}=60$ in table 4, as typical examples. The two gap ratios ${\it\Gamma}=2.0$ and 2.5 are chosen as representative cases where the antisymmetric or the symmetric mode of disturbance gives the critical condition, respectively. Needless to say, the real parts of ${\it\lambda}_{s00}$ and ${\it\lambda}_{a00}$ are positive in these tables because we are considering the cases where both modes are growing, only for which case the amplitude equations are applicable. It is noticeable that the real parts of ${\it\lambda}_{s10}$ , ${\it\lambda}_{s01}$ , ${\it\lambda}_{a10}$ and ${\it\lambda}_{a01}$ are all negative, which we assume for other values of ${\it\Gamma}$ in the following discussion on the existence and stability of the equilibrium solutions of the amplitude equations.

Table 4.

Coefficients of the coupled amplitude equations, $\mathit{Re}=60$ .

In order to explore nonlinear behaviors of the complex amplitudes we solved the coupled amplitude equations for ( $\widetilde{A}_{s}(t)$ , $\widetilde{A}_{a}(t)$ ) as an initial value problem and plotted their traces in phase diagrams in figure 5 for the two cases ${\it\Gamma}=2.0$ and 2.5 at $\mathit{Re}=60$ . In these figures, the horizontal and vertical axes stand for the absolute values $|\widetilde{A}_{s}|$ and $|\widetilde{A}_{a}|$ of the amplitudes of the symmetric and the antisymmetric modes. The open circles indicate the initial state ( $|\widetilde{A}_{s}(t=0)|$ , $|\widetilde{A}_{a}(t=0)|$ ) which can be chosen arbitrarily. The time-dependent solution proceeds along each line in the direction of the arrow and finally approaches to the equilibrium state indicated by filled circles. In figure 5(a), the flow converges to one of the two single-mode solutions, ( $|\widetilde{A}_{s}|\neq 0,|\widetilde{A}_{a}|=0$ ) or ( $|\widetilde{A}_{s}|=0,|\widetilde{A}_{a}|\neq 0$ ), depending on the initial condition for ${\it\Gamma}=2.0$ , which tells us that the equilibrium state is not determined uniquely, but depends on the initial state. Although there is no orbit converging to a mixed-mode solution ( $|\widetilde{A}_{s}|\neq 0,|\widetilde{A}_{a}|\neq 0$ ) because it is unstable, we can evaluate the amplitudes of the mixed-mode solution from (2.38), which is indicated by the cross symbol in figure 5(a). In the case of ${\it\Gamma}=2.5$ (figure 5 b), the solution converges to the mixed-mode solution indicated by the filled circle and the single-mode solutions (crosses in figure 6 b) are unstable so could not be attained. Thus, the equilibrium solution is uniquely determined independently of the initial condition in this case.

Figure 5.

Phase diagram of the solution obtained by the coupled Landau equations, $\mathit{Re}=60$ : solid line, orbit; ○, initial condition; ●, equilibrium amplitude; $\times$ , unsteady equilibrium solution; (a) ${\it\Gamma}=2.0$ ; (b) ${\it\Gamma}=2.5$ .

As we are interested in the equilibrium flow resulting from nonlinear interactions between the two modes of disturbance, we evaluate the equilibrium amplitudes of the two modes and investigate their stability. The equilibrium solutions ( $\overline{A}_{s}$ , $\overline{A}_{a}$ ) of the coupled amplitude equations (2.38) are expressed as $\overline{A}_{s}=|\overline{A}_{s}|\exp (\text{i}{\it\mu}_{s}t)$ and $\overline{A}_{a}=|\overline{A}_{a}|\exp (\text{i}{\it\mu}_{a}t)$ with $\text{d}|\overline{A}_{s}|/\text{d}t=0$ and $\text{d}|\overline{A}_{a}|/\text{d}t=0$ . There are four possible equilibrium solutions, i.e. the null solution $(|\overline{A}_{s}|=0,\ |\overline{A}_{a}|=0)$ , a symmetric single-mode solution $(|\overline{A}_{s}|\neq 0,\ |\overline{A}_{a}|=0)$ , an antisymmetric single-mode solution $(|\overline{A}_{s}|=0,\ |\overline{A}_{a}|\neq 0)$ , and a mixed solution $(|\overline{A}_{s}|\neq 0,\ |\overline{A}_{a}|\neq 0)$ .

The stability analysis of the equilibrium solutions is straightforward although it is lengthy and caution is required for the treatment of trivial eigenvalues, and so we will briefly summarize the procedure. Adding disturbances $\widehat{A}_{s}$ and $\widehat{A}_{a}$ to the equilibrium solutions $\overline{A}_{s}$ and $\overline{A}_{a}$ , we denote $\widetilde{A}_{s}=\overline{A}_{s}+\widehat{A}_{s}$ and $\widetilde{A}_{a}=\overline{A}_{a}+\widehat{A}_{a}$ . Substituting these expressions into (2.38), subtracting the equations for $\overline{A}_{s}$ and $\overline{A}_{a}$ from the resultant equations, and neglecting the nonlinear terms of $\widehat{A}_{s}$ and $\widehat{A}_{a}$ , we obtain linearized disturbance equations for $\widehat{A}_{s}$ and $\widehat{A}_{a}$ . Then, substituting $\widehat{A}_{s}=A_{s}^{\prime }\exp (({\it\zeta}+\text{i}{\it\mu}_{s})t)$ and $\widehat{A}_{a}=A_{a}^{\prime }\exp (({\it\zeta}+\text{i}{\it\mu}_{a})t)$ into the linearized disturbance equations, we obtain an eigenvalue problem for ${\it\zeta}$ . If real parts of all of the eigenvalues, ${\it\zeta}$ , are negative, then the equilibrium solution is stable; otherwise unstable.

The equilibrium amplitudes and the frequencies are easily obtained from (2.38) for each solution. The null solution is simply expressed as

which exists unconditionally and is apparently unstable because we have assumed that both modes are growing, i.e.  ${\it\sigma}_{s}>0$ and ${\it\sigma}_{a}>0$ . In fact, the eigenvalues for this solution are ${\it\zeta}={\it\sigma}_{s}$ and ${\it\sigma}_{a}$ both of which are double roots.

The symmetric single-mode solution is obtained as

which is possible only if ${\it\sigma}_{s10}<0$ . This condition ${\it\sigma}_{s10}<0$ means that the Hopf bifurcation is supercritical at the first bifurcation point where the steady symmetric flow becomes unstable only to the symmetric mode of disturbance. The single symmetric-mode solution exists in the whole range in figure 6, whose linear stability analysis gives the eigenvalues as ${\it\zeta}=-2{\it\sigma}_{s00}$ , ${\it\sigma}_{a00}+|\overline{A}_{s}|^{2}{\it\sigma}_{a01}$ (double root), and a trivial one ${\it\zeta}=0$ . Hence, the symmetric single-mode solution is stable if ${\it\sigma}_{a01}<-|\overline{A}_{s}|^{-2}{\it\sigma}_{a00}={\it\sigma}_{a00}{\it\sigma}_{s10}/{\it\sigma}_{s00}$ (region $\text{F}_{1}\text{F}_{2}\text{H}_{2}\text{H}_{1}$ in figure 6), otherwise unstable. Similarly, the antisymmetric single-mode solution

exists in the whole region and stable in $\text{F}_{1}\text{V}_{2}\text{V}_{1}\text{F}_{3}$ .

Figure 6.

Expected stable equilibrium solutions: $\text{S}_{s}$ , symmetric single-mode solution; $\text{S}_{a}$ , antisymmetric single-mode solution; $\text{S}_{m}$ , mixed-mode solution.

The amplitudes and frequencies are obtained for the mixed-mode solution as

The mixed mode exists if $({\it\sigma}_{s00}{\it\sigma}_{a10}-{\it\sigma}_{s01}{\it\sigma}_{a00})({\it\sigma}_{s10}{\it\sigma}_{a10}-{\it\sigma}_{s01}{\it\sigma}_{a01})<0$ and if $({\it\sigma}_{s10}{\it\sigma}_{a00}-{\it\sigma}_{s00}{\it\sigma}_{a01})({\it\sigma}_{s10}{\it\sigma}_{a10}-{\it\sigma}_{s01}{\it\sigma}_{a01})<0$ i.e. in $\text{F}_{1}\text{C}_{4}\text{G}_{1}\text{C}_{3}$ and $\text{G}_{1}\text{H}_{2}\text{C}_{2}\text{C}_{1}\text{V}_{1}$ in figure 6. Two eigenvalues for the mixed-mode solution are obtained as

and the other two are trivial. As ${\it\sigma}_{s01}<0$ and ${\it\sigma}_{a01}<0$ , the mixed solution is stable in $\text{G}_{1}\text{H}_{2}\text{C}_{2}\text{C}_{1}\text{V}_{1}$ and unstable in $\text{F}_{1}\text{C}_{4}\text{G}_{1}\text{G}_{3}$ .

To summarize, if we confine our discussion to the cases of ${\it\sigma}_{s01}<0$ and ${\it\sigma}_{a01}<0$ as well as ${\it\sigma}_{s10}<0$ and ${\it\sigma}_{a10}<0$ as the present numerical calculations show, the mixed-mode solution is stable and the two single-mode solutions are unstable in the region $\text{G}_{1}\text{G}_{2}\text{O}\text{G}_{3}$ , and vice versa in $\text{F}_{1}\text{V}_{2}\text{G}_{1}\text{H}_{1}$ . Furthermore, the trace of the coefficients moves along a line such as indicated by ${\it\alpha}{\it\beta}$ from $\text{G}_{1}\text{G}_{2}\text{O}\text{G}_{3}$ to $\text{F}_{1}\text{V}_{2}\text{G}_{1}\text{H}_{1}$ , or from $\text{F}_{1}\text{V}_{2}\text{G}_{1}\text{H}_{1}$ to $\text{G}_{1}\text{G}_{2}\text{O}\text{G}_{3}$ inevitably going through the point $\text{G}_{1}$ , though the line ${\it\alpha}{\it\beta}$ in figure 6 is schematic. Thus we arrive at a conclusion that the stability of the mixed-mode solution exchanges with both the two single-mode solutions simultaneously. This conclusion is valid to all cases of two-mode nonlinear interactions as far as the two amplitudes are governed by (2.38) with coefficients of ${\it\sigma}_{s01}<0$ and ${\it\sigma}_{a01}<0$ as well as ${\it\sigma}_{s10}<0$ and ${\it\sigma}_{a10}<0$ , i.e. to all cases of double Hopf bifurcations.

We demonstrate the changes of stability of the mixed mode and the two single modes of solutions in figure 7 for $\mathit{Re}=60$ where the real parts of ${\it\zeta}$ , $\text{Re}[{\it\zeta}]$ , for the three modes are depicted as functions of ${\it\Gamma}$ , where the solid and broken lines indicate the growth rate of disturbance for the symmetric and antisymmetric single modes solution, respectively, and the dotted line for the mixed-mode solution. At $\mathit{Re}=60$ , both the symmetric and antisymmetric modes of disturbance grow due to instability of the steady symmetric flow for ${\it\Gamma}\geqslant 1.61$ , and hence we depicted $\text{Re}[{\it\zeta}]$ only in the range of $1.61\leqslant {\it\Gamma}\leqslant 3.5$ . In this figure, the mixed mode equilibrium solution is stable and the two single-mode solutions are unstable in the range of $2.10\leqslant {\it\Gamma}\leqslant 3.41$ , while the mixed-mode solution is unstable and the two single-mode solutions are stable outside the range. Thus, we can confirm our conclusion derived above that the stability of the mixed-mode solution exchanges with both the two single-mode solutions simultaneously.

Figure 7.

Growth rate ${\it\zeta}_{r}$ of equilibrium solutions depending on the gap ratio, $\mathit{Re}=60$ : solid line, symmetric single-mode solution; broken line, antisymmetric single-mode solution; dotted line, mixed-mode solution.

Thus we have evaluated the equilibrium amplitude and examined the stability for the three modes of solutions and subdivided the region IV of figure 4 into IV-1 and IV-2 as in figure 8. In subregion IV-1, the two single-mode solutions are stable and the mixed-mode solution is unstable, where the equilibrium state is not uniquely determined. On the other hand, the region IV-2 indicates that only the mixed-mode solution is possible and the equilibrium state is uniquely determined through the nonlinear interactions between the two modes of disturbance.

Figure 8.

Expected stable equilibrium solutions: solid line, symmetric mode; broken line, antisymmetric mode; dotted line, boundary between IV-1 and IV-2; I, steady symmetric flow; II, antisymmetric single-mode solution; III, symmetric single-mode solution; IV-1, symmetric or antisymmetric single-mode solution; IV-2, mixed-mode solution.