## Abstract

Extensional self-similar flows in a channel are explored numerically for arbitrary stretching–shrinking rates of the confining parallel walls. The present analysis embraces time integrations, and continuations of steady and periodic solutions unfolded in the parameter space. Previous studies focused on the analysis of branches of steady solutions for particular stretching–shrinking rates, although recent studies focused also on the dynamical aspects of the problems. We have adopted a dynamical systems perspective, analysing the instabilities and bifurcations the base state undergoes when increasing the Reynolds number. It has been found that the base state becomes unstable for small Reynolds numbers, and a transitional region including complex dynamics takes place at intermediate Reynolds numbers, depending on the wall acceleration values. The base flow instabilities are constitutive parts of different codimension-two bifurcations that control the dynamics in parameter space. For large Reynolds numbers, the restriction to self-similarity results in simple flows with no realistic behaviour, but the flows obtained in the transition region can be a valuable tool for the understanding of the dynamics of realistic Navier–Stokes solutions.

## 1. Introduction

We study the behaviour of a viscous fluid confined within a two-dimensional infinitely long channel whose parallel walls are stretching or shrinking independently, i.e. the walls spatially accelerate or decelerate linearly and independently along their length. This problem has been frequently used as a paradigm for self-similar flows, i.e. exact solutions of the Navier–Stokes equations that satisfy spatial scale invariance properties in unbounded domains. Historically, self-similar solutions have been introduced to simplify the functional dependence of the flow on its unbounded coordinate. Multiple examples of self-similar flows can be found in the literature, ranging from the classical stagnation point flows by Hiemenz or modelizations of channel flows with porous walls with either suction or injection [1]. For a comprehensive review, we refer the reader to the monograph by Drazin & Riley [2]. Besides its inherent mathematical interest, self-similar solutions of extensional channels are being used to model a wide variety of multidisciplinary phenomena such as physiological blood flows and oxygen transport in cardiac systems [3]. Special cases of self-similar extensional flows in channels and pipes have been studied by many authors in the past addressing particular geometries and boundary conditions. For example, steady flows generated between parallel plates with the same stretching rates were studied by Brady & Acrivos [4]. Steady and time-periodic symmetry-breaking flow generalizations of similar solutions within the context of suction through porous walls were studied by Cox [5], Watson *et al.* [6] and King & Cox [7], who identified their primary instabilities as a result of pitchfork and Hopf bifurcations that eventually lead to period doubling cascades and chaos. On similar grounds, Cox [8] also concluded that the chaotic dynamics was extremely sensitive to symmetry breaking. Wang & Wu [9] addressed the existence of steady extensional flows with one of the channel walls stretching and the other stationary, thus breaking the symmetry beforehand.

Quite recently, Espín & Papageorgiou [10] have studied the range of validity of the aforementioned self-similar solutions by direct numerical simulation of the Navier–Stokes equations on a long bidimensional domain and with symmetric boundary conditions, i.e. with the plates shrinking or stretching simultaneously at the same rate. Espín & Papageorgiou [10] concluded that whereas the self-similar solution for symmetric stretching of the plates is valid for moderately large Reynolds numbers, its validity is extremely limited when the plates are shrinking due to strong instabilities exhibited by the flow.

In order to shed light on the underlying mechanisms responsible for the sensitivity of chaotic dynamics reported by Cox [8] or the lack of validity of the self-similar solution structure confirmed by Espín & Papageorgiou [10] we need to explore the space of steady or time-periodic solutions (stable or unstable) admissible by arbitrary shrinking–stretching boundary conditions. In other words, we need to unfold the state space by parametrizing the admissible solutions on the upper and lower wall acceleration or deceleration rates. Understanding qualitative changes in the dynamics and predicting flows arising from instabilities must be carried out using bifurcation theory. Bifurcation phenomena (i.e. the unexpected change in the qualitative dynamical behaviour of physical systems when external conditions are slightly varied) are ubiquitous in nature. An example of these changes are hydrodynamic instabilities in fluid flows such as the sudden oscillatory or turbulent motion exhibited by the wake left by a ship or an airplane when their speed exceeds some threshold value. Dynamical systems theory provides the mathematical building blocks to understand and predict in advance these bifurcations, therefore playing a very important role in many areas of physics and engineering. These building blocks are the generic bifurcations (fold, Hopf, etc.) that physical systems experience under parameter variation, regardless of the mechanisms underlying the dynamics. If the system is controlled by a single parameter, codimension-one bifurcations are expected but higher codimension bifurcations appear if the system depends on two or more parameters.

The paper is structured as follows. The mathematical formulation of the problem as well as the details of the spatio-temporal discretization of the Proudman–Johnson equation are given in §2. Section 3 describes the base state and primary instabilities it undergoes and localizes three cusp bifurcations that organize the subsequent bifurcations. Sections 4–6 describe the complex dynamics arising into the three cusp regions, and §7 presents a discussion and concluding remarks.

## 2. Governing equations

Let us consider an infinitely long two-dimensional channel with parallel extensional walls separated by distance 2ℓ, filled with an incompressible fluid of kinematic viscosity *ν*. All variables are rendered non-dimensional using ℓ as unit of length and the viscous time ℓ^{2}/*ν* as unit of time. Figure 1 shows a schematic representation of the channel geometry, using non-dimensional variables. The extensional channel walls have velocities proportional to the coordinate *x* along the channel. The two-dimensional Navier–Stokes equations for this problem are
**v**(*x*,*y*,*t*)=(*u*,*v*) is the two-dimensional velocity field in Cartesian coordinates. The channel domain is *y*=±1. The no-slip boundary conditions for the velocity field are
*σ*_{±} are the non-dimensional stretching rates of the walls. These are the two non-dimensional parameters of the problem. They can be positive or negative, depending on the extensional or compressional character of the wall motion. It is convenient to introduce two alternative parameters,
*α*∈[0,2*π*] measures the ratio between the stretching rates of the walls, *α* in terms of dimensional variables, indicated using an asterisk
*ψ*(*x*,*y*,*t*) as **v**=(*ψ*_{y},−*ψ*_{x}), where subindices *x*, *y* indicate partial derivatives. The equation for the streamfunction is
*ψ*: *ψ*_{y}(*x*,±1,*t*)=*σ*_{±}*x*, and *ψ*_{x}(*x*,±1,*t*)=0. The streamfunction equation is of fourth order, and we have four boundary conditions, but they are Neumann boundary conditions, so *ψ* is defined up to an additive function of time. In order to fix this constant, let us compute the mass flux *ϕ* in the streamwise direction *x*:
*ϕ* independent of *x*, depending only on time. The mass flux *ϕ*(*t*) may be imposed externally or even appear as a result of the channel flow dynamics. The difference between *ψ* at the channel walls is fixed by the imposed mass flux. We can use the gauge freedom in the election of *ψ* to impose the condition *ψ*(1)+*ψ*(−1)=0, that fixes the gauge. The final boundary conditions for *ψ* (that replace *ψ*_{x}(*x*,±1,*t*)=0) are
*ϕ*=0 will be considered in this study.

The form of the boundary condition (2.7) suggests the self-similar ansatz (for *ϕ*=0)
*y*-derivative. The equation for the self-similar function *f* is
*β*(*t*) resulting from the governing equation (2.9) measures the growth rate and has been extensively used in earlier papers on the subject. The last equation in (2.10) including *β* is known as the Proudman–Johnson equation [12].

The channel geometry is invariant under reflections on the coordinate axes, *K*_{x} and *K*_{y}. Translations along the streamwise direction also leave the channel invariant, but the stretching completely destroys the translational symmetry. The self-similar formulation (2.8) is equivariant under the reflection symmetry *K*_{x}, i.e. we are working in the *K*_{x}-invariant subspace, therefore only the *K*_{y} symmetry plays a role. The action of this symmetry on the function *f*, the parameters of the problem and the physical variables *u*, *v*, *p* and *ψ* is
*K*_{y}-equivariant only when *σ*_{+}=*σ*_{−}, i.e. for two values of *α*, *π*/4 and 5*π*/4.

Even when the governing equations have the *K*_{y} symmetry, the solutions may not have it. The base state (unique for sufficiently small *R*) is *K*_{y}-symmetric, but this symmetry can be broken at bifurcations, resulting in states with no symmetry at all. In the classical paper, Brady & Acrivos [4], the analysis was limited to solutions having all the symmetries of the problem (*α*=*π*/4), and as we will see, there are additional solutions breaking the symmetry. Moreover, the presence of symmetry changes the kind of bifurcations a dynamical system may undergo, as we have observed in this problem.

### (a) Numerical methods

For a uniform and simple treatment, it is convenient to separate *f* into two parts, one satisfying the non-homogeneous boundary conditions, and another satisfying homogeneous boundary conditions. So we introduce the split
*f*_{0} is the third-order polynomial
*f*_{1} is
*L* and *N* are the linear and nonlinear parts of the equation, respectively. For the time evolution, *f*_{1} is approximated via a spectral expansion in modulated Legendre polynomials of the form
*Φ*_{m}(±1)=*Φ*′_{m}(±1)=0, where *P*_{m}(*y*) are the Legendre polynomials. Then equation (2.15) is projected on the same basis via the standard Hermitian inner product in [−1,1]
*M* nonlinear ordinary differential equations (ODEs),
_{ℓm}, L_{ℓm} and N_{ℓ} are the linear and nonlinear operators
*h* is the time step and the superindex refers to the discrete times where the solution is obtained, *t*_{n}=*nh*.

For the computation and continuation in parameter space of the steady solutions, the fourth-order equation
*f*_{2}=*f*_{1}′′, we obtain
*D* denotes *y*-derivative. We expand *f*_{1} and *f*_{2} in basis functions satisfying the boundary conditions (two for *f*_{1} and none for *f*_{2}):
*f*_{1} and *f*_{2} are projected onto the Legendre basis (a Petrov–Galerkin method). The resulting nonlinear equations can be solved using the Newton method, and a pseudo-arclength continuation scheme has been implemented along any straight line in parameter space (*R*,*α*).

A pseudo-arclength Newton–Krylov–Poincaré continuation method for computing periodic solutions has been implemented [13,14]. Stable and unstable periodic orbits are obtained using this method. Finally, linear stability analyses of the fixed points and periodic solutions are carried out, focusing on the most unstable eigenvalues using Krylov–Arnoldi methods [15].

## 3. Base state and primary instabilities

For *R* small, and specified stretching–shrinking rates *σ*_{±}, the self-similar equations have a unique solution, that we call the base state. Figure 2 shows the base state at *R*=20, for 16 equispaced values of *α* in [0,2*π*], illustrating the different stretching–shrinking possibilities. Figure 3*a* shows the variation of the wall-normal velocity at the channel centre *V* =−*f*(0) versus *α* at *R*=20. The circles in figure 3*b* correspond to the 16 cases shown in figure 2. The symmetry about the line *α*=*π*/4, 5*π*/4 is easily recognisable from the diagrams shown in figure 2; for example, the diagram for *α*=0 is the horizontal reflection about *y*=0 of the *α*=*π*/2 diagram. The exchange between extensional and shrinking walls, corresponding to the rotation *α*→*π*+*α*, is not a symmetry of the problem, and the corresponding flows are quite different, as illustrated with the cases *α*=*π*/2 and *α*=3*π*/2.

When increasing the Reynolds number above *R*≈25, more than one steady state is obtained for certain *α*-values. Four cusp bifurcations appear at parameter values given in table 1. Two of them, *C*_{1} and *C*_{3}, take place on the symmetry line *α*=*π*/4 and 5*π*/4; the other two, *C*_{2} and *α*=*π*/2 and 2*π*. Figure 4*a* shows the variation of the wall-normal velocity at the channel centre *V* versus *α* at *R*=300. Figure 4*b*–*d* are close-ups of the three cusps *C*_{1}, *C*_{2} and *C*_{3}, respectively. The cusp *C*_{3}, is not shown. These curves have been obtained by continuation of steady solutions of the governing self-similar equation (2.9). At each of the saddle-node bifurcations in figure 4*b*–*d*, a real eigenvalue crosses the imaginary axis. Figure 5 shows the structure of the flow using streamlines, for the three solutions coexisting inside the cusp bifurcation; the middle branch solution is always unstable, while the upper and lower branches are stable close enough to the corresponding cusp point. For the cusps *C*_{1} and *C*_{3}, the upper and lower branches are symmetrical, and the corresponding flows are *K*_{y} symmetric, as shown in the figure.

When increasing the Reynolds number *R*, the basic state undergoes a variety of bifurcations depending on the *α* value. A comprehensive exploration of the parameter space (*R*,*α*)∈[0,500]×[0,2*π*) has been carried out, and a variety of Hopf bifurcations and codimension-two bifurcations, in addition to the four cusps, also appear. They are summarized in figure 6 and will be discussed in detail in the next sections. Some of the bifurcations and instabilities in figure 6 have already been reported in previous studies. For example, the *α*=5*π*/4 line of figure 6 corresponds to the particular case of symmetric shrinking (i.e. the case with the two walls decelerating at the same spatial rate) studied by many authors in the past. The cusp *C*_{1} is a parametric unfolding of the pitchfork bifurcation reported, for example, in Watson *et al.* [6] or in Espín & Papageorgiou [10]. Similarly, the line *α*=*π*/2 (or its symmetric counterpart *α*=0) corresponds to the particular case of one stretching wall and the other stationary, reported by Wang & Wu [9], that also leads to a pitchfork bifurcation at cusp *C*_{2}. Finally, line *α*=*π*/4 corresponds to the case of symmetric stretching addressed by Brady & Acrivos [11] and is associated with the pitchfork bifurcation cusp *C*_{3}.

### (a) Limit cycle solutions

In the *α* ranges between the cusp bifurcations, the base state becomes unstable for moderate *R*-values by means of Hopf bifurcations, where a periodic orbit solution is born. Between the cusps *C*_{1} and *C*_{2}, the Hopf bifurcation *H*_{1} takes place at small values of *R*. We describe in detail the bifurcated solutions for *α*=*π*, corresponding to a shrinking upper wall, while the other wall is at rest. To our knowledge, this case has not been explored before in the context of extensional flows. However, similar results in the case of flows driven by suction through porous walls were found by King & Cox [7].

Figure 7*a* displays the bifurcated stable periodic orbit for *α*=*π* and various Reynolds numbers; the figure shows phase portraits of the limit cycles projected on the plane (*U*,*V*)=(*f*′(0),−*f*(0)), the streamwise and wall-normal velocities at the channel centre, for *x*=1. The base state becomes unstable at *R*=46.34. When increasing *R* the limit cycle grows substantially in size: the maximum value of |*U*| increases from |*U*|=263.8 at *R*=200 to |*U*|=8111 at *R*=500, i.e. a factor of 30 for an increase in *R* of a factor 2.5. The time series develops a sharp slow/fast behaviour when increasing *R*, as shown in figure 7*b* for *R*=200; the slowing down takes place near the upper right corner of the limit cycle. This dynamical behaviour is typical of homoclinic or saddle-node in cycle (SNIC) bifurcations of limit cycles, where either the periodic orbit approaches a homoclinic connection of a saddle equilibrium, or where the saddle is born within the cycle, respectively. However, the two aforementioned scenarios have been discarded by monitoring the dependence of the period of the limit cycles as a function of *R*. Figure 7*c* reveals that the period of the limit cycles grows moderately when increasing *R* and apparently remains bounded and without any indication of divergence, contrary to what would be expected in SNIC or homoclinic scenarios. When increasing *R*, the slow–fast dynamics of the limit cycles become more pronounced and distinguishable along the periodic orbits with fast transient flow stages developing narrow boundary layers. As shown in figure 7*c*, this increase in *R* also leads to an exponential growth of the *β* factor that modulates *f*′′′ in equation (2.10) and that is responsible for the quadratic growth of the pressure in the streamwise *x* coordinate of the channel. The high values of *β* and the associated very thin boundary layers for high *R* are a clear indication that in real flows other phenomena will appear, destroying the self-similar character of the solutions, as already reported by Espín & Papageorgiou [10] for the case of decelerating wall flows (*α*=5*π*/4). Decelerating or shrinking walls occur in the interval *α*∈(*π*/2,2*π*), where the mentioned limit cycles originate along the Hopf curves *H*_{1} and *C*_{1} around *α*=5*π*/4, the limit cycles developing very narrow boundary layers are the only stable solutions of the self-similar equations. Away from the *C*_{1} cusp region, no additional fixed points (apart from the unstable base state) have been identified. These facts, and the difficulties in the numerical solution of the governing equations due to the narrow boundary layers, are the most likely explanation of why the two generic shrinking cases for *α*=3*π*/4 (shrinking–stretching walls) and *α*=*π* (shrinking–stationary walls) have not been studied previously. For high values of *R*, the exponential growth of *β* along with the presence of very thin boundary layers makes the numerical computations really challenging. All spectral computations reported in this study have used a minimum value of *M* to ensure a decay of at least 10 orders of magnitude in the spectral Legendre coefficients *a*_{m} of the expansion (2.16) and a suitable decrease of the time-step *h* to avoid numerical instabilities or noticeable changes of the dynamics when it is reduced. For most of the computations, *M*=70 Legendre polynomials and *h*=2×10^{−3}*R*^{−1} were sufficient to satisfy the numerical accuracy and stability requirements. However, the spatial and temporal resolutions needed to be increased in order to properly solve the aforementioned narrow boundary layers associated with the slow–fast dynamics of the limit cycles generated above curves *H*_{1} and

Figure 7*d*–*i* show the structure of the limit cycle for *α*=*π* and *R*=200. The solution along the cycle presents two very different behaviours. For nearly half of the period, the flow exhibits four recirculating cells, as illustrated with the streamlines in figure 7*d*. The profiles of the streamwise *u* in black (blue online) and wall-normal *v* in grey (orange online) velocities in figure 7*f* are almost symmetrical, and the *K*_{y} symmetry is only broken inside the thin boundary layers near the walls (the boundary conditions are not *K*_{y}-symmetric, except for *α*=*π*/4 and 5*π*/4). We call this solution symmetric (Sym in the figure). There is also a small fraction of the period when the solution has only two recirculating cells, as illustrated with the streamlines in figure 7*e*. The profiles of *u* and *v* in figure 7*g* are clearly not *K*_{y}-symmetric. We call this solution asymmetric (asym in the figure). There are also two intervals of transition between the sym and asym phases of the limit cycle, indicated in figure 7*h* for the phase-portrait, and in figure 7*i* for the pressure factor *β*. The transition and asymmetric phases become a smaller fraction of the limit cycle period when increasing the Reynolds number *R*. The sudden and large increase in *β* takes place in the asymmetric phase.

The other regions where the primary instability is due to a Hopf bifurcation are within the intervals *α*∈(0,*π*/4) between cusps *C*_{3}, and *α*∈(*π*/4,*π*/2), between *C*_{3} and *C*_{2}; see figure 6. In these two ranges, the Hopf bifurcations take place at much higher Reynolds numbers (typically ranging between 300 and 500) along curves *H*_{2} and its symmetric counterpart *α*=*π*/8 the Hopf bifurcation takes place at *R*=465.5, leading to time periodic regimes whose main features are outlined in figure 8. These periodic orbits also exhibit slow–fast dynamics when increasing *R*, as can be easily seen in figure 8*b*. However, the pressure factor *β* attains moderate values when compared with the decelerating walls case previously analysed. The orbital period exhibits an almost linear growth as a function of *R*, without any trace or signal of potential blow up, thus discarding SNIC or homoclinic secondary bifurcations of the period orbits. Within the range of explored values we have not identified secondary instabilities of these limit cycles.

Figure 8*d* shows how the Hopf frequency (the frequency of the limit cycle at the Hopf bifurcation point) varies along the two Hopf branches *H*_{1} and *C*_{2} end of the *H*_{1} branch, and at both ends of the

## 4. Secondary bifurcations inside cusp *C*_{1}

The interaction of the Hopf bifurcation curves and the cusp *C*_{1} is shown in figure 9*a*. The *H*_{1} Hopf curve and the fold curve *C*_{1} become tangent at the codimension-two fold-Hopf bifurcation point *FH*^{S}, with parameter values (*α*,*R*)_{FHS}≈(4.0207,113.14). The fold-Hopf bifurcation admits different scenarios, depending on the specifics of the dynamical system considered. In order to ascertain which scenario corresponds to our problem, we have computed phase portraits of the solutions at *R*=125, for different *α* values, as indicated by square symbols in figure 9*b*. The phase portraits corresponding to points a–c shown in the first row of figure 10 exhibit the two steady points emerging from the saddle-node bifurcation curve *a*. At figure 10*b*, we can see the limit cycle emerging from the saddle point (middle fixed point in the figure), at the Hopf curve *H*_{1}. It is an unstable limit cycle; for *R*<*R*_{FHS}, the limit cycle is stable, as has been discussed in the previous section. The unstable limit cycle becomes stable in a Neimark–Sacker bifurcation curve *NS*_{−}, as shown in figure 10*c*. These dynamics correspond to the third scenario of the fold-Hopf bifurcation, as described in Kuznetsov [16]. The two-torus that is born at *NS*_{−} is unstable, and it is in fact the boundary between the basins of attraction of the stable lower and upper limit cycles shown in figure 10*c*. This boundary set is an invariant set that is inherently unstable and may have a very complex structure. It has been studied in some detail in other fluid problems, as in pipe flow [17,18], where it has been shown to have a fractal structure, and there are unstable states in the boundary set, called edge states, that play a significant dynamical role. In our case, the boundary set is an unstable two-torus close to the *NS*_{−} bifurcation curve, and we have not explored in detail its subsequent evolution, although some considerations will be made about it when exploring the codimension-two *PN*_{−} point.

When decreasing the parameter *α* at *R*=125, the stable limit cycle in figure 10*c*, born at the *NS*_{−} bifurcation curve (the lower limit cycle in the figures), undergoes a period doubling cascade (figure 10*d*) starting at the bifurcation curve *PD*_{−}. The cascade results in the strange attractor shown in figure 10*e*. A further decrease in *α* results in the disappearance of the strange attractor; the only stable solution that remains is the upper limit cycle, as shown in figure 10*f*. Some of the unstable periodic orbits associated with the strange attractor survive, as the unstable limit cycle depicted in figure 10*f*. The disappearance of the strange attractor happens along the curve *Cr*_{−}. This is a global bifurcation called a boundary crisis: the strange attractor collides with the boundary of its basin of attraction, becoming unstable [19]. At the collision the boundary set (that was a two-torus near the *NS*_{−} bifurcation curve) is destroyed; the basin of attraction of the strange attractor is absorbed by the basin of attraction of the upper stable limit cycle.

The period doubling cascade ending at the boundary crisis happens in a very narrow parameter range. Continuation of the bifurcation curves *NS*_{−}, *PD*_{−} and *Cr*_{−} at larger Reynolds numbers *R*>125 results in their merging in a codimension-two bifurcation point *PN*_{−} at (*α*,*R*)_{PN−}≈(3.990,136.95). This is a codimension-two bifurcation of the lower limit cycle, undergoing simultaneously a period-doubling and a Neimark–Sacker bifurcation. It can be understood as a codimension-two bifurcation of the Poincaré map associated with the lower limit cycle, where three eigenvalues cross the unit cycle: an eigenvalue −1, corresponding to the period-doubling, and a couple of complex-conjugate eigenvalues of modulus one, corresponding to the Neimark–Sacker bifurcation [16]. This complicated bifurcation has not been studied in full detail yet, although some particular scenarios in another system have been recently analysed [20].

The bifurcation curves *NS*_{−}, *PD*_{−} and *Cr*_{−} have also been continued at lower *R* and *α* parameter values. The period doubling cascade interval widens when approaching the symmetry line at *α*=5*π*/4. Around this line, we have two period doubling cascades, associated with the upper and lower limit cycles, starting at the curves *PD*_{+} and *PD*_{−} respectively, as shown in figure 9. Around the symmetry line at *α*=5*π*/4 both period-doubling cascades overlap partially, as shown in figure 9*b*. Phase portraits illustrating the period-doubling cascade *PD*_{+} at *R*=125 are shown in figure 10*g*–*l*. On the symmetry line *α*=5*π*/4 both strange attractors coexist, as shown in figure 10*k*; in fact the two strange attractors are *K*_{y} symmetric. When we move slightly away from the *α*=5*π*/4 line, one of the attractors undergoes a boundary crisis and disappears (see figures 10*j* and *l*). The boundary crisis bifurcation curves *Cr*_{+} and *Cr*_{−} are asymptotic to the *α*=5*π*/4 line. Therefore, although the strange attractors exist for large *R* values on the symmetry line (corresponding to both walls shrinking at the same rate), a very small *K*_{y}-symmetry breaking results in dramatic changes in the flow. The period-doubling curves *PD*_{+} and *PD*_{−} also approach the *α*=5*π*/4 line when increasing *R*. Therefore all the complex dynamics becomes confined in a very narrow *α* interval when increasing *R* and can only be observed in detail at moderate *R* values (between about *R*=105 and *R*=125). It is also worth mentioning that the strange attractors depicted in figure 10*k* exhibit the slow down characteristic of the limit cycles analysed in §3 (see also figure 7*a*). The slow down in this case corresponds to the strange attractors approaching the saddle fixed point along the line *V* =0 in figure 10*k*. In this figure, we can see that the strange attractor, moving away from the saddle, closely follows the unstable manifold of the saddle.

The precise values of the succession of period doublings have been computed for *R*=125 (varying *α*) and on the symmetry line *α*=5*π*/4 (varying *R*). The Newton–Krylov–Poincare method combined with Arnoldi stability analysis allows one to monitor the Floquet exponents of the limit cycles, and a precise determination of the critical parameters of the period doubling cascade. The results are shown in table 2. P2, P4, P8… refer to the succession of period doubled orbits; *δ*_{n} measures the ratio of each bifurcation interval to the next
*μ* refers to *R* or *α*, the parameter varied along the period doubling cascade. The first period doubling *n*=1 results in P2, *n*=2 results in P4 and so on. We have been able to compute up to five period doublings using continuation. We observe that in both cases analysed, the parameter *δ*_{n} for large *n* tends to the Feigenbaum constant *δ*=4.6692016…, a universal constant in many period doubling processes [21,22]. The period doubling cascade has been observed previously in the particular case *α*=5*π*/4 (the symmetric case) by Watson *et al.* [6]. In this work, two successive period doublings were observed, and we have included the corresponding values in the last column of the table; the critical parameters *R*_{n} were computed with an accuracy of at least three figures. The remaining values shown in table 2 have been computed with five figure accuracy; the first three period doublings are illustrated in figure 10*g*–*i*. We have also been able to find limit cycles of period three (P3) and five (P5), and their windows of stability; the presence of an unstable cycle of period three is a clear indication of the presence of a strange attractor [23,24], with unstable limit cycles of all periods multiple of the basic limit cycle LC. These P3 and P5 cycles are born in cyclic-fold bifurcations (a saddle-node bifurcation of cycles), at Reynolds numbers *R*_{P3}=114.10093 and *R*_{P5}=113.26726. In these bifurcations two limit cycles, one stable and the other unstable, are born, as shown in figure 11. Although the cycles exist for *R*>*R*_{P3} and *R*>*R*_{P5}, respectively, the interval of stability *δR* of the initially stable branch is very narrow: *δR*_{P3}=0.27775, and *δR*_{P5}=0.03787.

The results obtained in the present study can be quantitatively compared with those obtained by Watson *et al.* [6] for the self-similar flow at *α*=5*π*/4, i.e. for the case of parallel plates with the same shrinking rates. In this work, the critical Reynolds numbers for the cusp bifurcation *C*_{1} and the Hopf bifurcation *H*_{1} are given as *R**_{C1}=17.30715 and *R**_{H1}=55.77. The Reynolds number *R** used in Watson *et al.* [6] is related with *R* used in our computations according to *R*_{C1}=24.4760 and *R*_{H1}=78.8380 and converted to *R** result in 17.3071 and 55.7469, respectively. The agreement is very good.

## 5. Secondary bifurcations inside cusp *C*_{2}

Figure 12*a* shows the secondary bifurcations and codimension-two points associated with the asymmetric cusp *C*_{2}. The Hopf curve *TB*_{1}, as was suggested by figure 8*d*. The parameter values of *TB*_{1} are (*α*,*R*)_{TB1}≈(1.7181,150.52). In contrast, the Hopf curve *H*_{1} does not interact with the bifurcation curves associated with *C*_{2}, and figure 12*a* suggests that it approaches *F*_{2} for large values or *R*, outside the scope of the present study.

There is another Takens–Bogdanov point *TB*_{2} on the fold curve *F*_{2}, at parameter values (*α*,*R*)_{TB2}≈(1.5823,367.70). The numerical exploration shows that it is a standard Takens–Bogdanov bifurcations that does not interact with the other branches inside the cusp *C*_{2} and takes place in a very narrow region in parameter space, as shown in the inset in figure 12*b*. Figure 13 shows three phase portraits of the solutions crossing the bifurcations associated with *TB*_{2}, corresponding to the points a, b and c in figure 12*b*. Figure 13*a* shows the two unstable fixed points born at the fold bifurcation *F*_{2}. In figure 13*b*, one of the new fixed points becomes stable in a Hopf bifurcation *H*_{3}, and an unstable limit cycle is born. This unstable limit cycle disappears in an homoclinic collision *Hom*_{2} with the other fixed point born at *F*_{2}, in an infinite period bifurcation. After the homoclinic connection, typical of the Takens–Bogdanov bifurcation [16], only the two fixed points remain, as illustrated in figure 13*c*. The three bifurcation curves *F*_{2}, *H*_{3} and *Hom*_{2} are tangent at *TB*_{2}.

The bifurcation curves emerging from *TB*_{1}, in contrast with *TB*_{2}, cover a wide region inside the cusp and exhibit complex dynamics. The Hopf bifurcation curve *Hom*_{1} emerging from *TB*_{1} interact via two new codimension-two bifurcations, joined by the curve of cyclic-folds *CF*_{1}. These new bifurcations are shown in figure 12*c*.

Figure 14 illustrates the phase portraits corresponding to the points a, b, c and d for *R*=200 in figure 12*c*, close to the Takens–Bogdanov bifurcation *TB*_{1}. Figure 14*a* for *α*=1.7 shows the upper and lower stable fixed points along with the middle saddle unstable one. Orbits emerging from the saddle that surround the upper fixed point grow in size when *α* is increased, eventually leading to the homoclinic connection shown in figure 14*b*. The homoclinic connection for larger *α* values results in an unstable periodic orbit, as shown in figure 14*c*. This unstable periodic orbit collapses into the upper fixed point, in the Hopf bifurcation *d*.

For *R*=200 near the Takens–Bogdanov point *TB*_{1}, the Hopf bifurcation curve is subcritical, in the sense that the stable point is surrounded by an unstable limit cycle that bounds its basin of attraction. In §2 it was shown that the Hopf bifurcations along the curve *α*∈(*π*/4,*π*/2) are supercritical, in the sense that the limit cycle is stable and does not coexist with the stable fixed point. The simplest mechanism responsible for the subcritical to supercritical change of the Hopf bifurcation along *Ba*_{1}, located at (*α*,*R*)_{TB1}≈(1.7114,230.75), a bifurcation curve *CF*_{1} emerges. Along this curve, the stable and unstable limit cycles, born at the subcritical and supercritical Hopfs, merge in a fold bifurcation of cycles (a cyclic-fold) and disappear.

Figure 15 shows phase portraits at the different regions delimited by the bifurcation curves above the Bautin bifurcation, corresponding to the points a–f in figure 12*c* at *R*=250. The first two phase portraits, figure 15*a* and *b*, show the crossing of the Hopf curve *H*_{1}, where the stable limit cycle around the lower fixed point is born; as we mentioned before, the Hopf curve *H*_{1} does not interact with the dynamics associated with the Takens–Bogdanov point *TB*_{1} and the Bautin point *Ba*_{1}, and indeed the dynamics on the remaining phase portraits in figure 15 take place in the upper part of the phase space, around the middle and upper fixed points. By increasing *α* the supercritical Hopf curve *H*_{2} is crossed and a stable limit cycle around the upper fixed point appears (see figure 15*c*). Figure 15*d* shows the formation of the homoclinic connection with the middle saddle point, on the curve *Hom*_{1}; a further increase in *α* results in the break-up of the homoclinic loop, giving rise to an unstable limit cycle, as shown in figure 15*e*. These two limit cycles approach each other when *α* is increased, merging on the bifurcation curve *CF*_{1} and disappearing, as shown in figure 15*f*.

Increasing the Reynolds number above *R*=250, the bifurcation curves *CF*_{1} and *Hom*_{1} approach each other, and they meet at a codimension-two bifurcation *CFH*_{1} located at (*α*,*R*)_{CFH1}≈(1.7053,300.0). This is a codimension-two global bifurcation of limit cycles: the cyclic-fold where the stable and unstable limit cycles are born happens when both cycles collide homoclinically with the middle saddle point. For *R* above the *CFH*_{1} point, the cyclic-fold *CF*_{1} curve no longer exists, and the dynamics simplifies: only the homoclinic collision of the stable upper limit cycle with the middle saddle remains. Figure 16 shows the crossing of the *Hom*_{1} curve at *R*=320, and the three phase portraits correspond to the points a, b, c at *R*=320 in figure 12*c*.

## 6. Secondary bifurcations inside cusp *C*_{3}

The third and last cusp *C*_{3} shown in figure 6 unfolds the flow dynamics generated between parallel plates with the same stretching rates (*α*=*π*/4) originally studied by Brady & Acrivos [4]. The bifurcation diagram in a neighbourhood of the cusp is shown in figure 17*a*, where the dynamics is governed by a Takens–Bogdanov bifurcation *TB*_{3} located at (*α*,*R*)_{TB3}=(0.67910,401.45). Similarly to what has been just described in cusp *C*_{2}, the *TB*_{3} point is located at the confluence of the Hopf (*F*_{3} bifurcation curves, and an homoclinic (*Hom*_{3}) bifurcation curve characteristic of the Takens–Bogdanov bifurcation emerges from the same point. Also in this case, close to *TB*_{3} the Hopf bifurcation *α*>*π*/4 as was shown in §2) is supercritical, and the change in behaviour takes place at a Bautin codimension-2 point *Ba*_{3} located at (*α*,*R*)_{Ba3}=(0.72289,494.97). From the Bautin point *Ba*_{3} emerges a cyclic-fold transition curve *CF*_{3} that presumably should merge with *Hom*_{3} for larger values of *R*. The aforementioned bifurcation scenarios and corresponding phase portraits are qualitatively identical to the ones described in previous section for the cusp *C*_{2}, in particular those addressed in figure 12*c* along path *a*–*d* for *R*=200 (depicted in figure 14), path *a*–*f* for *R*=250 (depicted in figure 15) and path *a*–*c* for *R*=320 (depicted in figure 16). The bifurcation curves *Hom*_{3} and *CF*_{3} and their symmetric counterparts approach the symmetry line at *α*=*π*/4 for large *R* values. Outside this narrow region of complex dynamics, there only remains one stable solution, a limit cycle around the upper or lower unstable fixed point, depending on the side of the symmetry line considered. The complex dynamics becomes confined at intermediate *R* values (between about *R*=400 and *R*=550) and for larger *R* it exists only in a very narrow *α* interval around the symmetry line *α*=*π*/4.

A quantitative comparison of the results obtained in the present study with those obtained by Espín & Papageorgiou [10] is shown in table 3, for the self-similar flow at *α*=*π*/4, i.e. for the case of parallel plates with the same stretching rates. The critical Reynolds numbers agree with five figures, and the period of the stable limit cycle solution at *R*=523.26 differs in less than 0.5%, so the agreement is fully satisfactory. The Reynolds number and the time scale used in [10] are different from those used in the present study, and in table 3 we have converted our values to theirs.

In the classic work of Brady & Acrivos [4] on extensional channel flow, for the case *α*=*π*/4, *K*_{y}-symmetric steady solutions were computed up to *R*≈127 600, and three different branches of solutions were computed. Branch I in the mentioned paper corresponds to the base state in the present study and to the mid-branch after the cusp *C*_{3} bifurcation. Branches II and III were disconnected from branch I and appeared at a fold bifurcation at about *R*≈438. A disconnected branch of solutions cannot be computed by continuation methods, neither by time evolution if they are unstable, so finding disconnected branches is a difficult problem, and one never knows if there are additional branches not yet found. In this case, we used the fact that branches I and II become very close at high *R* to compute them; this is how branches II and III were obtained in the mentioned reference. The result is shown in figure 17*b*, fold curve *F*_{4}; on the symmetry line *α*=*π*/4, the fold takes place at *R*=434.13, very close to the value reported in Brady & Acrivos [4]. Apart from the *K*_{y}-symmetric steady solutions, there are other asymmetric fixed points. In Watson *et al.* [6], the two asymmetric branches born at the cusp *C*_{3} were found, and also the Hopf bifurcation *H*_{2} (for *α*=*π*/4). We have performed a detailed exploration of the fixed points up to *R*=850 and *α*∈[0,*π*/2], and the result is shown in figure 17*b*. A new cusp bifurcation of one of the branches born at the fold *F*_{4} has been found at *R*=744.4, with the corresponding fold bifurcation curves *F*_{5} and

Figure 18*a*,*c* shows the variation of *V* with *α* for *R*=480 and 800, for the branches of fixed points. The new branch born at the fold *F*_{4} is the closed curve, limited by two fold bifurcations. The open curve shows the solutions associated with the cusp *C*_{3}. At *R*=800 it can observed that the closed branch has undergone a cusp bifurcation, and the closed curve exhibits a couple of new fold points very close to the *α*=*π*/4 line. The new fixed points are shown in figure 18*b*,*d* which are the corresponding phase portraits for *R*=480 and 800 at *α*=*π*/4. All the new fixed points are unstable and are away from the region where the complex dynamics associated with *C*_{3} takes place. The phase portraits in figure 18 also show the unstable manifolds associated with the new fixed points that evolve towards the stable solutions (fixed points and limit cycles associated with *C*_{3}). Therefore, the presence of these new fixed points does not modify the dynamics already discussed.

Figure 19(1–7) shows the streamlines corresponding to the seven fixed points in figure 18*d*. The first and third rows are the three solutions associated with the cusps *C*_{3} and *C*_{4}. The central column shows the three *K*_{y}-symmetric solutions, being 5 and 7 the solutions born at the fold bifurcation *F*_{4} (branches II and III in Brady & Acrivos [4]). Solutions 2 and 4 in the central column are very similar, corresponding to the mentioned approach of branches I and II in Brady & Acrivos [4] (also shown in figure 18*d*), at high Reynolds numbers. Figure 19*a*,*b* outlines the variation of the pressure parameter *β*/*R* and the wall-normal velocity *V* as a function of the Reynolds number *R* for two fixed values of *α*. These figures clearly show the successive bifurcations that the base state undergoes when increasing *R*, and also the unfolding of the pitchfork bifurcations resulting from the *K*_{y}-symmetry breaking for a slightly increased value *α*=0.82. In particular, figure 19*a* allows a direct comparison with figure 2*a* of Brady & Acrivos [4].

## 7. Discussion and conclusion

A comprehensive exploration of the dynamics of the self-similar extensional channel flow has been carried out. Historically, the self-similar studies of this flow focused on the computation of steady solutions and their asymptotic properties at large Reynolds numbers [4,9,12]; however, these large *R* solutions are usually unstable, and more recent works have focused on dynamical aspects [6,10]. All these previous studies focused on specific values on the stretching ratio of the walls *α*; only the cases *α*=*π*/4, 5*π*/4 and *π*/2 have been previously studied, corresponding to both walls stretching or shrinking at the same rate, and one wall stretching and the other at rest. What we have found is that the interesting dynamics takes place in *α* intervals centred around the *α*=*π*/4 and 5*π*/4 cases, and also close to but not including the *α*=*π*/2 case. For *α* values outside these three regions, corresponding to the three cusps *C*_{1}, *C*_{2} and *C*_{3}, the base state undergoes a Hopf bifurcation resulting in a periodic solution, and there are no additional bifurcations when increasing *R*.

It has also been found that the base state becomes unstable (in the self-similar subspace) for small Reynolds numbers, and a transitional region including complex dynamics takes place at intermediate *R* values, depending on the value of *α*. For *α*∈(0.57*π*,1.93*π*), where at least one of the walls is shrinking, the transition takes place in the approximate range *R*∈(25,140). For *α*∈(0,*π*/2), where both walls are stretching, the transition takes place in the range *R*∈(200,600). In the remaining region, associated with cusp *C*_{2}, the transition takes place in the range *R*∈(100,450). For Reynolds numbers above the transition region, the only stable solution is a limit cycle, whose structure has been discussed in detail in §3a. This periodic solution develops slow–fast dynamics when increasing *R*, and the boundary layers become very thin with a large pressure parameter *β*. The limit cycle is not unique but appears in pairs in a narrow region around the symmetry lines *α*=*π*/4 and 5*π*/4, where a pair of symmetrically related limit cycles are present, and in the case *α*=5*π*/4 a strange attractor persists at high *R*. When the symmetry is broken, one of the limit cycles (and the strange attractor) disappears, leaving only a single periodic solution; the flow is very sensitive to symmetry-breaking.

The transition region is very rich in a variety of codimension one and two bifurcations and is an excellent arena for applying dynamical systems theory. In particular, the complex time-dependent flows emerge from, and are locally governed by, the explored codimension-two bifurcation points.

A natural question arises, about the feasibility of these flows in real problems. The flows obtained in this study are exact solutions of the Navier–Stokes equations, but restricted to the self-similar subspace. Do these solutions persist, are they observable and do they play a dynamical role when solving the full Navier–Stokes equations without restrictions? The answer is already partially known. Espín & Papageorgiou [10] compared the self similar solutions with DNS of the Navier–Stokes equations in large but finite domains, in the case where both walls have the same stretching rates (the cases *α*=*π*/4 and 5*π*/4). In this study, it was established that *all branches of the self-similar solutions, including time-periodic ones, can occur in finite channels provided the Reynolds numbers are not too high: approximately R=500 for accelerating and R=33 for decelerating wall flows, respectively. At sufficiently high Reynolds numbers, the self-similar inflow conditions are incapable of producing the exact solution in the interior, but instead new steady or time-periodic states emerge.* Other authors have also pointed out that the self-similar solutions at high

*R*are not feasible. Brady & Acrivos [11] explicitly say that

*the results suggest that similarity solutions should be viewed with caution because they may not represent a real flow once a critical Reynolds number is exceeded.*These results and considerations strongly suggest that the solutions found below and into the transition region could be observed in real flows, while the solutions above the transition region do not represent any real flow. In fact, the general behaviour of the solutions of the Navier–Stokes equations is laminar flow at low

*R*, a transition region where complex spatio-temporal flows emerge at intermediate

*R*, and turbulent flow at high

*R*. Restricting the computations to a self-similar subspace suppresses the turbulent phase, and in fact above the transition region only simple flows with no realistic behaviour are found (periodic solutions with narrow boundary layers and very large pressure parameters). However, the analysis of the basic flow and part of the transition region (for moderate Reynolds numbers) in the self-similar subspace can be a valuable tool for the understanding of the dynamics of realistic Navier–Stokes solutions.

## Data accessibility

The datasets and codes supporting this article have been uploaded as supplementary material.

## Authors' contributions

F.Ma. identified the bifurcation scenarios, did most of the computations and wrote the paper with A.M.; A.M. provided most of the numerical schemes, computed steady solutions and wrote the paper with F.Ma.; F.Me. provided the continuation schemes and linear stability analysis for periodic orbits and computed some of them; P.W. provided the extensional formulation and helped with the writing. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by the Spanish Government grant nos. FIS2013-40880-P and FIS2016-77849-R.

## Acknowledgements

The Catalan Government support of the authors in UPC, members of the research group SGR 1515; F.Me. is an Associate Professor of the Serra Húnter Programme.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3789979.

- Received March 3, 2017.
- Accepted May 5, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.