A Study of Bifurcation of Kolmogorov Flows with

an Emphasis on the Singular Limit

Hisashi Okamoto^{1}

Abstract.

We consider a family of stationary Navier-Stokes flows in 2D flat tori. The flow is driven by an outer force which is of the form (sin y, 0). Varying the Reynolds number and the aspect ratio of the torus, we numerically compute bifurcating solutions by a path-continuation method. Folds and cusps are obtained in the range where the Reynolds number is < 100. Some solutions are computed up until the Reynolds number becomes 10,000. Asymptotic properties as the Reynolds number tends to infinity are discussed. Also given is an analysis as the aspect ratio of the torus tends to zero.

1991 Mathematics Subject Classification: Primary 76D30; Secondary 76C05, 35Q30, 35Q35. Keywords and Phrases: Kolmogorov flows in 2D tori, incompressible fluid, bifurcation, singular perturbation, internal layer, inviscid limit.

1 Introduction

The Navier-Stokes equations have attracted very much attention of both mathematicians and physicists; accordingly scientific papers on them are almost innumerable. Nonetheless, many difficult problems remain to be analyzed; this is especially true when the Reynolds number is large ( see [7] and [4] ). One of the purposes of the present paper is to point out that something new can be found even if we restrict ourselves to steady-states.

1Supported by the Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan, # 09304023, # 09554003

Documenta Mathematica · Extra Volume ICM 1998 · III · 513–522

We compute numerically a family of stationary motions of incompressible viscous fluid, which is governed by the following Navier-Stokes equations:

u\frac{∂u}{∂x} + v \frac{∂u}{∂y} = \frac{1}{R}△u − \frac{∂p}{∂x} + \frac{1}{R} sin y,  (1)
u\frac{∂v}{∂x} + v \frac{∂v}{∂y} = \frac{1}{R}△v − \frac{∂p}{∂y} ,  (2)
\frac{∂u}{∂x} + \frac{∂v}{∂y} = 0,  (3)

where R is the Reynolds number, (u, v) the velocity, and p the pressure. Those equations are satisfied in (x, y) ∈ [−π/α, π/α] × [−π, π], with periodic boundary condition. Namely we consider the Navier-Stokes flows in a two-dimensional flat torus and α is its aspect ratio. See [5] and [11].

If the driving force is (sin ky, 0) with k ≥ 2, then interesting phenomena are already known for nonstationary motions as well as steady-states ( see [2], and the references in [11] for instance ). In our problem, which comes from Iudovich [5], the flow is driven by an outer force (R^{−1} sin y, 0). This simplifies the problem greatly. We would like to refer the reader to [3] and [11], where motives of investigation and historical comments are found. The purpose of the present paper is to report that we can observe many interesting phenomena when we change the aspect ratio α as well as the Reynolds number R. Varying α and R, we numerically compute bifurcating solutions ( [11] ). Folds and cusps are obtained in the range where R < 100. Some solutions are computed until R becomes 10,000 ( [13] ). We hope that such a list of solutions serves as raw materials for future study of the Navier-Stokes equations. In particular, we would like to obtain a bifurcation diagram which is global in the sense that solutions of all the parameters (α, R) are listed in the diagram. Such global bifurcation diagrams are computed in many onedimensional systems, notably reaction-diffusion systems ( [10] ). However, global diagram for the Navier-Stokes equations are substantially more difficult to obtain and we are forced to be content with a partial answer, which we are going to present in the present paper. The following study of Kolmogorov flows is motivated by A. Majda’s pioneering works on incompressible fluid motions ( see, e.g., [7] ) and Nishiura’s analysis of reaction-diffusion systems [10].

From our numerical computation, we can guess some interesting asymptotic behavior as the Reynolds number tends to infinity. Since the Navier-Stokes equations are defined in a 2D torus, there can not be a boundary layer. However, we can observe some internal layers, which we will explain in the next section. In order to analyze those internal layers, we apply in section 3 a singular perturbation method in the range where (α, R) ≈ (1, ∞). The solution obtained by the perturbation method shows a good agreement with the numerical solution. Analysis when α → 0 is given in section 4.

Figure 1: Neutral curves of mode n ( n = 1, 2, ⋯) (left). Nontrivial solutions bifurcate from the points on the neutral curves. The curve of mode n starts from (α, R) = (0, \sqrt{2}) and ends at (α, R) = (1/n, ∞). Schematic bifurcation diagram of solutions of mode 1 (right). The point A represents (α, R, ψ) = (1.0, ∞, −(cos x + cos y)/2). Only the upper half of the bifurcating solutions are drawn.

2 Global picture of solutions and inviscid limit

We first recall some numerical facts reported in [11]. We discretize (1)–(3) by the spectral method. The resulting nonlinear equations are solved by the pathcontinuation method ( see [6], for instance ). One easily notice that (u, v, p) = (sin y, 0, 0) solves the equations and the boundary conditions for all R > 0 and all α > 0. We call it a trivial solution. It is Meshalkin and Sinai [9] which proves that any bifurcation from the trivial solution occurs by steady-states. Namely, the Hopf bifurcation from the trivial solution is prohibited. Iudovich [5] showed that there are bifurcation from the trivial solution if 0 < α < 1 and that there is none if 1 ≤ α < ∞. See Figure 1 (left). Bifurcating solutions are classified by a positive number called a mode. Roughly speaking, solutions of mode n contains n pairs of eddies in the rectangle (−π/α, π/α) × (−π, π). See [11]. When 0 < α < 0.966 ⋯, then with R as a bifurcation parameter, there exists a pitchfork of bifurcating solutions. There is no secondary bifurcation in the class of those solutions satisfying ψ(x, y) = ψ(−x, −y). When 0.966 ⋯ < α < 1, the branch of nontrivial solutions possesses two turning points ( = limit points ) but still there is no secondary bifurcation in the same function class ( [11] ). We recently re-computed the stability of those solutions in the function space where we do not assume any symmetry. We have found that the solutions are stable even in this general setting. Therefore the nontrivial solutions on the pitchfork is stable however large the Reynolds number may be. The global view of the solutions of mode one is given in Figure 1 (right). See [11] for detail. This suggests that a possibility that the global attractor for 1/2 < α < 1 is of one-dimension however large the Reynolds number may be. Such low-dimensionality is reported in a different context in Afendikov and Babenko [1] and Chen and Price [2].

The pitchfork bifurcations are supercritical for all α ∈ (0, 1). Namely, the nontrivial solutions in a neighborhood of the bifurcation point lie in the right hand side, where the Reynolds number is greater than the critical Reynolds number. This was shown numerically in [11] but a mathematical proof was not available there, although the supercriticality for sufficiently small α is proved by Afendikov and Babenko [1], and independently by M. Yamada. See [11]. Recently Matsuda and Miyatake [8] gave a proof of supercriticality when α is close to one.

We now consider the asymptotic behavior of the solutions as R → ∞. When R increases with a fixed α, the nontrivial solution seems to converge on a certain function. The numerical experiments in [13] suggests that the vorticity as a limit of R → ∞ is at most C^{1}. See Figure 2, which shows that (−△)^{−3/2}ψ ( ψ is the stream function ) loses smoothness along certain curves. We call these curves internal layers. The layer yields an energy spectrum of k^{−r}, where k is the wave number and r is between −7 and −4, depending on the aspect ratio. See [13] for detail. We remark that our “singularity” is much weaker than those found in [7].

Since the knowledge of solutions with large R may help us understand the turbulent motion of fluid, it would be of practical importance to mathematically analyze an asymptotic behavior of steady-states as R → ∞. In the present case, asymptotic analysis seems to be very difficult for a general α. First of all, we encounter the following problem: The Euler equations, which are obtained from (1)–(3) by setting R = ∞, possess an infinite number of solutions. In fact, they have a continuum of steady-states. On the other hand, the Navier-Stokes equations have finitely many steady-states for a fixed R < ∞. Therefore the vast majority of the stationary Euler flows are disconnected with the Navier-Stokes flows. Hence we would like to know how the Euler flows which are connected with the Navier-Stokes flows are distinguished from those which are disconnected. Some partial answers are given in [12, 13] but we do not know the real mechanism for it. We will show in the next section that a certain heuristic analysis is possible for those solutions which lie in the neighborhood of the point A in Figure 1 (right).

3 Asymptotic analysis as (R, α) → (∞, 1)

The equations (1)–(3) are written by the stream function as follows ( see [11] ):

\frac{△^{2}ψ + cos y}{R} + ψ_{x}△ψ_{y} − ψ_{y}△ψ_{x} = 0,  (4)

where the subscript means differentiation. Substituting x = x′/α, y = y′ and dropping the primes, we obtain

0 = (α^{2}∂_{x}^{2} + ∂_{y}^{2})^{2} ψ + cos y + Rα [ψ_{x} (α^{2}∂_{x}^{2} + ∂_{y}^{2}) ψ_{y} − ψ_{y} (α^{2}∂_{x}^{2} + ∂_{y}^{2}) ψ_{x}] , (5)

Figure 2: Graphs of (−△)^{−3/2}ψ. Bird’s-eye views (left) and slices of the graphs along the line y = αx (right). α is 0.7 (top), 0.984 (center), and 0.999 (bottom), respectively. The equations are discretized by the Fourier-Galerkin method. The resulting nonlinear equations, which contains 544 to more than 5,000 independent variables depending on the Reynolds number R, are solved by the pathcontinuation method. which is satisfied in −π < x < π, −π < y < π. Let δ = 1/(Rα) and γ = 1 − α^{2}. Then, by defining J(f, g) = f_{x}g_{y} − f_{y}g_{x}, we obtain

0 = δ (△^{2}ψ + cos y − 2γ△ψ_{xx} + γ^{2}ψ_{xxxx}) + J(ψ, △ψ) − γJ(ψ, ψ_{xx}). (6)

We now consider those solutions which are close to the point A in Figure 1 (right). We expand γ ∈ R and ψ as follows:

ψ = ∑_{j,k=0}^{∞} ϵ^{j}δ^{k}ψ^{j,k}(x, y),  γ = ∑_{j,k=0}^{∞} ϵ^{j}δ^{k}γ(j, k),  (7)

where ϵ is an artificial parameter. It is taken along the vertical tangent of surface of solution set at (α, R) = (1, ∞). See Figure 1(right) and Figure 3 (left). γ(0, 0) = 0 is assumed so as to comply with the numerical results. Substituting (7) into (6), we compute coefficients of ϵ^{j}δ^{k}. Then we first obtain J(ψ^{0,0}, △ψ^{0,0}) = 0. We already know from the numerical results in [11, 13] that ψ^{0,0} = −(cos y ± cos x)/2. Since both cases are dealt with in the same way, we choose ψ^{0,0} = −(cos y + cos x)/2. From the coefficient of ϵ^{0}δ^{1}, we obtain

△^{2}ψ^{0,0} + cos y + J(ψ^{0,0}, △ψ^{0,1}) + J(ψ^{0,1}, △ψ^{0,0}) − γ(0, 1)J(ψ^{0,0}, ψ_{xx}^{0,0}) = 0,

which is written as

\frac{cos y − cos x}{2} − γ(0, 1)\frac{sin x sin y}{4} + \frac{1}{2} (sin x∂_{y} − sin y∂_{x}) (I + △)ψ^{0,1} = 0.

Figure 3: Two coordinates near the turning pint A (left). Plots of αR(−1)^{n}n(1 − 2n^{2})a(n, n). α = 0.9999 and R = 3000, 5000, 8000. (right)

Let us define the operator K by K = \frac{1}{2} (sin x∂_{y} − sin y∂_{x}). Note that

K (log |\begin{matrix}cos & \frac{x−y}{2} \\ cos & \frac{x+y}{2} \end{matrix} | ) = \frac{1}{2}(cos y − cos x)  and  K(cos x − cos y) = sin x sin y.

Note also that a function u = u(x, y) satisfies Ku = 0 if and only if there exists a function of one variable f such that u(x, y) = f(cos x + cos y). These facts lead us to

log |\begin{matrix}cos & \frac{x−y}{2} \\ cos & \frac{x+y}{2} \end{matrix} | + \frac{1}{4}γ(0, 1)(cos y − cos x) + (I + △)ψ^{0,1} + f(cos x + cos y) = 0

with some function f. Multiplying this equation with cos x − cos y, we integrate it on [−π, π]^{2}. Then we obtain γ(0, 1) = 0, whence

log |\begin{matrix}cos & \frac{x−y}{2} \\ cos & \frac{x+y}{2} \end{matrix} | + (I + △)ψ^{0,1} + f(cos x + cos y) = 0.  (8)

If we further assume that f ≡ 0, then we obtain

ψ^{0,1}(x, y) = − ∑_{n=1}^{∞} \frac{2(−1)^{n−1}}{n(1 − 2n^{2})} sin nx sin ny + c_{1} cos x + c_{2} cos y,  (9)

where c_{j}’s are constant.

Because of the limitation of the paper size we do not compute other coefficients in the present paper. Even with our incomplete computation, we can guess a certain interesting asymptotic behavior as (α, R) → (1, ∞). In fact,

ψ = −\frac{cos x + cos y}{2} + ϵψ^{1,0}
+\frac{1}{αR} ∑_{n=1}^{∞} \frac{2(−1)^{n}}{n(1 − 2n^{2})} sin nx sin ny + smooth function + ⋯ (10)

shows that the nonsmooth function appearing in (9) is a dominant factor for a large ( but not too large ) wave number range when (R, α) → (∞, 1). Figure 3 (right) shows the plot of αRn(1 − 2n^{2})a(n, n). This figure indicates that

ψ = −\frac{cos x + cos y}{2} + \frac{2}{αR} ∑_{n=1}^{∞} \frac{(−1)^{n}}{n(1 − 2n^{2})} sin nx sin ny + ⋯  (11)

is a good approximation to the solutions on the turning points in Figure 1 (right) in an intermediate wave number space.

4 Kolmogorov flows of small aspect ratio

The purpose of the present section is to consider the asymptotic behavior of the solutions of (5) as α → 0 with a fixed R. The stationary solution of the Navier-Stokes equations are expanded into the Fourier series: ψ(x, y) = ∑_{(m,n)≠(0,0)} a(m, n) exp(iαmx + iny). Then the Fourier coefficients satisfy the following equations:

\frac{1}{R} (α^{2}j^{2} + k^{2})^{2} a(j, k) + \frac{1}{2R}δ_{k,±1}δ_{j,0}

− ∑_{p=−∞}^{+∞} ∑_{q=−∞}^{+∞} αa(p, q)a(j − p, k − q)(kp − qj)(α^{2}jp + kq) = 0, (12)

where Kronecker’s delta is used. In particular we obtain

a(j, 0) + R ∑_{p,q} a(p, q)a(j − p, −q)\frac{pq}{αj^{2}} = 0.  (13)

This suggests the following asymptotic relation:

a(j, 0) = O(α^{−1})  as  α → 0,

which we assume from now on. Also we assume that

a(j, k) = O(1)  as  α → 0  ( k ≠ 0).

These asymptotic relations are compatible with our numerical experiment, which we can not show because of the page limitation. We now define b(j, k) as follows: b(j, k) = lim_{α→0} a(j, k) for k ≠ 0 and b(j, 0) = lim_{α→0} αa(j, 0). Then, we have the following equations:

\frac{1}{R}k^{4}b(j, k) + \frac{1}{2R}δ_{k,±1}δ_{j,0} − ∑_{p≠j} b(p, k)k^{3}(p − j)b(j − p, 0) = 0, (k ≠ 0) (14)
b(j, 0) + R ∑_{p≠0,q≠0} b(p, q)b(j − p, −q)pqj^{−2} = 0.  (15)

After some computation, with the aid of symmetry b(−j, k) = b(j, −k) = (−1)^{j+k−1}b(j, k), we can prove that b(j, k) = 0 if |k| > 1. Then we can rewrite (14) by means of {b(j, 1)}_{j=−∞}^{+∞} only:

b(j, 1) + \frac{1}{2}δ_{j,0} − R^{2} ∑_{p−j:odd} ∑_{s≠0} b(p, 1)\frac{2s(−1)^{s−1}}{j − p} b(s, 1)b(j − p − s, 1) = 0. (16)

This equation is supplemented by the equation (15) with q = ±1.

We have seen that the Navier-Stokes equations (12), which are written by the two-dimensional array {a(m, n)}, are reduced to a system of nonlinear equations of a one-dimensional array {b(j, 1)}. We have thus achieved a substantial reduction. However, the reduced equation contains a new difficulty. In fact, the equation (16) has a trivial solution such that b(0, 1) = −1/2, b(j, 1) = 0 (j ≠ 0). Linearizing (16) at this trivial solution, we can easily see that the system (16) possesses one and only one bifurcation point, which degenerates with infinite multiplicity. Consequently, the set of solutions of (16) near the trivial solution would look like Figure 4 (left). Note, however, that this figure is based on a naive guess from the linear analysis and the truth may well be different.

Figure 4: Infinite number of pitchforks at (α, R) = (0, \sqrt{2}) (left). Graph of the stream function of mode 1. α = 0.02, R = 100.0. Since x ranges from −50π to 50π, it is rescaled to a scale similar to y (right).

Now let us come back to the equation (4). Figure 4 (right) is the graph of the numerical stream function when α is small. This and the Fourier analysis above suggest that

ψ ∼ \frac{f(αx)}{α} + g(αx) cos y + h(αx) sin y + O(α)  as α → 0,  (17)

where f, g, and h are functions of one variable. Figure 4 (right) shows that f(ξ) = μ(R) (|ξ| − \frac{π}{2} ), where μ(R) is a constant depending on R. Substituting (17) into (5), we obtain

g(ξ) = \frac{−1}{1 + R^{2}(f′(ξ))^{2}} , h(ξ) = \frac{Rf′(ξ)}{1 + R^{2}(f′(ξ))^{2}} .

Thus, the solutions are calculated up to order O(1) as α → 0. However, further expansion accompanies a substantial difficulty.

The analysis of the singular perturbation problems of this and the preceding sections will be left to the future work.

References

[1] A.L. Afendikov and K.I. Babenko, Bifurcations in the presence of a symmetry group and loss of stability of some plane flows of a viscous fluid, Soviet Math. Dokl., vol. 33 (1986), pp. 742–747.

[2] D. Armbruster et al., Symmetry and dynamics for 2-D Navier-Stokes flows, Physica D, vol. 95 (1996), pp. 81–93. M.S. Jolly, Bifurcation computations on an approximate inertial manifold for the 2D Navier-Stokes equations, ibid, vol. 63 (1993), pp. 8-20. Z.-M. Chen and W.G. Price, Long-time behavior of Navier-Stokes flow on a two-dimensional torus excited by an external sinusoidal force, J. Stat. Phys., vol. 86 (1997), pp. 301–335.

[3] V. I. Arnold, Kolmogorov’s hydrodynamic attractors, Proc. R. Soc. Lond. A, vol. 434 (1991), pp. 19–22.

[4] P. Constantin, Some mathematical problems of fluid mechanics, in Proceedings if ICM 1994, ed. S. G. Chatterji, Birkhäuser (1995), pp. 1086–1095.

[5] V. I. Iudovich, Example of the generation of a secondary stationary or periodic flow when there is loss of stability of the laminar flow of a viscous incompressible fluid, J. Appl. Math. Mech., vol. 29 (1965), pp. 527–544.

[6] H.B. Keller, Lectures on Numerical Methods in Bifurcation Theory ( Tata Institute of Fundamental Research No. 79 ), Springer Verlag (1987).

[7] A.J. Majda, The interaction of nonlinear analysis and modern applied mathematics, Proceedings of ICM 1990, ed. I. Satake, Springer Verlag, (1991), pp. 175–191.

[8] M. Matsuda and S. Miyatake, Bifurcation curves on Kolmogorov flow, preprint.

[9] L.D. Meshalkin and Y.G. Sinai, Investigation of the stability of a stationary solution of a system of equations for the plane movement of an incompressible viscous liquid, J. Appl. Math. Mech., vol. 25 (1962), pp. 1700–1705.

[10] Y. Nishiura, Global structure of bifurcating solutions of some reactiondiffusion systems, SIAM J. Math. Anal., vol. 13 (1982), pp. 555–593. Y. Nishiura and H. Fujii, Stability of singularly perturbed solutions to systems of reaction-diffusion equations, ibid., vol. 18 (1987), pp. 1726–1770.

[11] H. Okamoto and M. Shōji, Bifurcation diagrams in the problem of incompressible viscous fluid flows in 2D tori, Japan J. Indus. Appl. Math., vol. 10 (1993), pp. 191–218.

[12] H. Okamoto, A variational problem arising in the two dimensional Navier-Stokes equations with vanishing viscosity, Appl. Math. Lett., vol. 7 (1994), pp. 29–33.

[13] H. Okamoto, Nearly singular two-dimensional Kolmogorov flows for large Reynolds number, J. Dynamics and Diff. Eqns., vol. 8 (1996), pp. 203–220.

Hisashi Okamoto

Research Institute for Mathematical Sciences, Kyoto University, Kyoto, 606-8502 Japan, okamoto@kurims.kyoto-u.ac.jp