## Abstract

With the goal of providing the first example of application of a recently proposed method, thus demonstrating its ability to give results in principle, global stability of a version of the rotating Couette flow is examined. The flow depends on the Reynolds number and a parameter characterizing the magnitude of the Coriolis force. By converting the original Navier–Stokes equations to a finite-dimensional uncertain dynamical system using a partial Galerkin expansion, high-degree polynomial Lyapunov functionals were found by sum-of-squares of polynomials optimization. It is demonstrated that the proposed method allows obtaining the exact global stability limit for this flow in a range of values of the parameter characterizing the Coriolis force. Outside this range a lower bound for the global stability limit was obtained, which is still better than the energy stability limit. In the course of the study, several results meaningful in the context of the method used were also obtained. Overall, the results obtained demonstrate the applicability of the recently proposed approach to global stability of the fluid flows. To the best of our knowledge, it is the first case in which global stability of a fluid flow has been proved by a generic method for the value of a Reynolds number greater than that which could be achieved with the energy stability approach.

## 1. Introduction

Hydrodynamic stability is the field that investigates the transient effects of an initial perturbation of a known steady flow. The area has attracted the attention of many researchers and is closely related to the study of transition to turbulence [1,2].

Using Lyapunov stability theory, a steady flow can be proved to be stable with respect to perturbations of arbitrary amplitude by constructing a Lyapunov functional *V* [**u**], which is a positive-definite functional of the velocity perturbation **u** that decays monotonically on any non-zero solution **u**(*t*,**x**) of the Navier–Stokes equations [3].

Usually, the Lyapunov functional is chosen to be the perturbation energy *E*=∥**u**∥^{2}/2, where the norm is defined as an integral of |**u**|^{2} over the flow domain. Then the problem of proving that *E* is a Lyapunov functional reduces to a tractable linear eigenvalue problem [4,5]. However, the resulting estimates of the global stability range, usually expressed by the energy stability limit Reynolds number *Re*_{E}, could be very conservative in the sense that *Re*_{E} is generally far below the maximum *Re* for which the flow is globally stable.

Recently, a method was proposed by Goulart & Chernyshenko [6] for exploiting the sum-of-squares (SOS) decomposition [7,8] to construct polynomial Lyapunov functionals differing from *E*, thus extending the range of *Re* in which the flow can be proved to be globally stable. In this approach, the Navier–Stokes equations are first reduced to a finite-dimensional uncertain dynamical system, that is a system of ordinary differential equations (ODEs) with right-hand side containing terms for which only bounds, but not exact expressions, are available. For incompressible flows both the right-hand side of the ODEs and the bounds are polynomial. The corresponding Lyapunov stability condition can then be reduced to a condition of positive definiteness of other certain polynomials [6]. Noting that the condition of a polynomial being positive-definite can be replaced by a stronger, but more tractable numerically, condition of the polynomial being a SOS of other polynomials, an admissible Lyapunov functionals can be found using the polynomial SOS optimization approach [7,9].

The SOS technique has been applied in stability analysis of constrained ODEs [10], hybrid systems [8] or time-delay systems [11], but there are few results for partial differential equations. The relevant publications are [12] and [13], where SOS-based algorithmic methodologies are presented for the analysis of systems described by certain types of parabolic partial differential equations. It was shown how certain Lyapunov structures could be constructed to prove stability using transformations defined through integration by parts. It is worth noting that in both [12,13], the partial differential equations are considered directly, which is different from the construction process of a Lyapunov functional in [6].

Application of SOS of polynomials to fluid flows goes beyond nonlinear stability. After an overview of nonlinear stability applications including a short announcement of a part of the results of this work, applications for deriving rigorous bounds on time-averaged characteristics of turbulent flows and applications to flow control are discussed in [14].

While [6] provides a full theoretical description of the new approach, it was applied there only to a truncated Galerkin approximation rather than to the full Navier–Stokes equations, thus leaving open the question of whether there is at least one fluid flow for which this method will work. Such an example is given in this paper. To the best of our knowledge, this study is the first case in which global stability of a fluid flow has been proved by a generic method for the value of a Reynolds number greater than that which could be achieved with the energy stability approach.

## 2. Problem formulation

Our goal is to demonstrate that the method proposed in [6] can actually be used to prove the stability of a fluid flow for the values of the Reynolds number above the energy stability limit. This section describes briefly the method and the flow to which it is to be applied for achieving this goal.

### (a) The method

An unsteady flow of incompressible viscous fluid in a given domain with time-independent boundary conditions is considered. The perturbation velocity **u**(*t*,**x**), defined as the deviation of the instantaneous flow velocity from the steady solution *A* depends on the flow in question. It is convenient to leave it in a compact general form here. This formulation is a little more general than in [6], where *A* had a particular form, but all the results of Goulart & Chernyshenko [6] apply with obvious minor modifications, which are made without further comments in the summary of the method given in this section. The perturbation velocity **u** is subject to homogeneous boundary conditions. The steady flow can be proved to be stable with respect to perturbations of arbitrary amplitude by constructing a Lyapunov functional *V* [**u**].

#### (i) The Lyapunov functional

For this purpose, in [6] the perturbation velocity is represented as
**u**_{i}, *i*=1,…,*k*, are an orthonormal set of solenoidal vector fields with the inner product 〈**w**_{1},**w**_{2}〉 defined as the integral of **w**_{1}⋅**w**_{2} over the flow domain **u**_{s} is solenoidal and orthogonal to all the **u**_{i}, and both **u**_{i} and **u**_{s} satisfy the homogeneous boundary conditions. The Lyapunov functional is sought in the form^{1} *V* [**u**]=*V* (**a**,*q*^{2}), where **a**=(*a*_{1},…,*a*_{k}), and *q*^{2}=∥**u**_{s}∥^{2}/2=〈**u**_{s},**u**_{s}〉/2. For *V* [**u**] to be a Lyapunov functional, the function *V* (**a**,*q*^{2})>0 and d*V*/d*t*<0, for all (**a**,*q*^{2})≠0. With the use of (2.1) and (2.2), the latter condition can be rewritten as [6]
**f** are
*Γ* and *χ* are^{2}
**Θ**(**u**_{s},**a**)=**Θ**_{a}(**u**_{s})+**Θ**_{b}(**u**_{s},**a**)+**Θ**_{c}(**u**_{s}) has the components
**h**_{ij}=**u**_{j}⋅∇**u**_{i}−∇**u**_{j}⋅**u**_{i}. The notation used can be clarified by the Einstein equivalent of the formula for **h**_{i}: *x*^{m} are the *m*th components of the vectors **h**_{i}, **u**_{i} and **x**, respectively.

Constructing *V* satisfying (2.3) and *V* >0 for all (**a**,*q*^{2})≠0 would prove the global stability of the flow under consideration. However, (2.3) involves **u**_{s} while *V* is a function of **a** and *q*^{2} only. Hence the next step is required [6].

#### (ii) The bounds

The terms in (2.3) dependent on **u**_{s} can be bounded by functions of **a** and *q*^{2}. Note that *χ*, *Θ*_{ai}, and *Θ*_{bi} are linear functionals of **u**_{s}, while *Γ* and *Θ*_{ci} are quadratic functionals of **u**_{s}.

For *Θ*_{ai} defined by (2.7) the Cauchy–Schwarz inequality gives **h**_{i} onto the solenoidal subspace and a small modification to account for the boundary conditions [6]. Other linear functionals can be bounded similarly.

The linear functional *χ* is a special case. If the basis **u**_{j} is chosen to consist of the eigenfunctions of the classical energy stability problem [15] for *χ*≡0 [6]. This can become obvious if one notices the similarity between the energy stability operator [15] and (2.6). We will use such a basis in this study. Accordingly, even though we refer to the following formulae as obtained in [6], in fact they are simplified versions that we derived with an additional assumption *χ*=0. In most cases, this can be done by simply omitting some terms. The readers wishing to use the method with the basis **u**_{i} in which *χ*≠0 should refer to [6] rather than to the formulae below.

The bounds on the quadratic functionals can be obtained by maximizing the functionals subject to a constraint *q*^{2}=1. This reduces to a linear eigenvalue problem in a fairly standard way. The resulting bound for *Γ* has the form [6]
*κ*_{s} reduces to the energy stability problem with an additional constraint 〈**u**_{s},**u**_{i}〉=0 ∀ *i*. When the linear eigenvalue problem has a discrete set of eigenvalues, only a finite number of them can be positive [15] and, hence, *κ*_{s} can be made negative by selecting a suitable basis **u**_{i}.

Rather than using the generic approach to bounding *Θ*_{ci} as proposed in [6], in appendix A we derive an explicit tight bound for it.

Putting together all the above gives the following set of bounds
*κ*_{s} and the coefficients of the quadratic polynomial *p*(**a**,*q*^{2}) depend on the particular flow in question and on the selection of the particular eigenfunctions of the corresponding energy stability problem to be used as the finite basis **u**_{i}.

This allows formulating the stability analysis problem as a SOS optimization problem.

#### (iii) Reduction to a polynomial sum-of-squares optimization

Since the term *Γ*(**u**_{s}) in (2.3) is upper-bounded by a negative-definite function *κ*_{s}*q*^{2} but is not lower-bounded, one has to impose an additional requirement that the candidate function *V* satisfy

The condition (2.10) ensures that the term (∂*V* /∂(*q*^{2}))*Γ*(**u**_{s}) makes a negative contribution to the left-hand side of the Lyapunov condition (2.3). Combining the condition (2.10) with the bounds (2.9) gives that (2.3) holds if
*V* and *p* enter it nonlinearly. This can be circumvented via introduction of additional variables. As shown in [6], (2.11) is equivalent to
*V* simultaneously satisfying the three conditions *V* >0, (2.10) and (2.12), for all (**a**,*q*^{2}≠0) then the flow in question is globally stable.

If *V* is sought for in the form of a polynomial then checking each of these three conditions amounts to checking the global non-negativity of a polynomial function, which is known to be NP-hard; see [9]. However, a sufficient, and numerically tractable, condition for global non-negativity of a polynomial is that it can be written as a SOS of other polynomials. Accordingly, a sufficient condition for joint satisfaction of our three conditions is that
*Σ*_{l} represents the set of all SOS polynomials in **z**≠0 and/or (**a**,*q*^{2})≠0.

Existence of a function *V* satisfying the conditions (2.14) can be checked via the polynomial SOS approach of [7,9], which amounts to solving a convex optimization problem in the form of a large semi-definite program. For a prescribed degree of the candidate polynomial *V* , the solution to that problem either provides an explicit expression for *V* or states that such *V* does not exist. The SOS approach is well known and will not be described here.

#### (iv) Uncertain system interpretation

Substituting the partial Galerkin expansion (2.2) into the Navier–Stokes equations (2.1) and projecting onto **u**_{i} gives
*q*^{2} of the residual field **u**_{s} is
**Θ**, *Γ* and *χ* are defined in §2a(i) as functionals of **u**_{s} and functions of **a**. One can, however, allow **Θ**, *Γ* and *χ* in (2.15) and (2.16) to assume any values as far as they satisfy the set of known bounds defined in terms of **u**_{s} and **a**, such as (2.9). In this sense, ((2.9), (2.15), (2.16)) is an uncertain dynamical system. The solution of this system is therefore not unique. However, if all the solutions of ((2.9), (2.15), (2.16)) tend to zero as time tends to infinity, then the solution of the Navier–Stokes system also tends to zero. The stability conditions described above are in fact the stability condition for this uncertain system. This turns out to be quite useful in understanding and interpretation of the further results.

Note that *V* (**a**,*q*^{2}) is a Lyapunov function for the uncertain system and at the same time a Lyapunov functional for the full Navier–Stokes equations.

If one takes *V* =∥**a**∥^{2}/2+*q*^{2} then (2.11) reduces to the classic energy stability problem. Hence, this approach [6] is guaranteed to give at least as good results as the energy stability approach. In order to demonstrate that it is capable of providing a better result, it has to be applied to a particular flow.

### (b) Double-periodic rotating Couette flow

Given the novelty and complexity, and the relatively demanding computational requirements for solving the semi-definite programming problems stemming from polynomial SOS optimization, for the first application of the stability analysis method of Goulart & Chernyshenko [6], it is desirable to select as simple a flow as possible. The particular flow we select is a version of the famous rotating Couette flow between two co-axial cylinders.

#### (i) Governing equations

The gap between the cylinders is assumed to be much smaller than the cylinder radius. A local Cartesian coordinate system **x**=(*x*,*y*,*z*) is oriented such that the axis of rotation is parallel to the *z*-axis, while the circumferential direction corresponds to the *x*-axis. Only flows independent of *x* are considered. The flow velocity is represented as (*y*+*u*,*v*,*w*), so that **u**=(*u*,*v*,*w*) is the velocity perturbation and *a*
*b*
*c*
*d*where *Ω* is a non-dimensional parameter characterizing the Coriolis force, *Re* is the Reynolds number and *p* is pressure. More compactly, (2.17) can be written in the vector form (2.1) with

In §2a, *A* was assumed to contain only the terms depending on the base flow *Ω* are present. Conveniently, their presence does not affect any of the formulae in §2a, mostly because these terms correspond to the Coriolis force and thus do not enter the energy equation.

For simplicity, the flow is assumed to be 2*π*-periodic in *y* and *z*, *u* and *v* are assumed to be odd in *y* and even in *z*, while *w* is assumed odd in *z* and even in *y*:

#### (ii) Stability properties of the flow

We first apply the well-known energy stability approach. Setting the Lyapunov functional as the perturbation energy *E*=∥**u**∥^{2}/2 leads to a linear eigenvalue problem [15]. For (2.17)–(2.18), the resulting eigenfunctions **e**_{n,m}(**x**), as can be verified by direct substitution into that eigenvalue problem, are:
*n*=1,2,…,*m*=0,±1,±2,…. The corresponding eigenvalues are

Note that for this flow, conveniently, neither the eigenfunctions nor the eigenvectors depend on *Ω*. The eigenvalue λ_{1,−1} is positive for _{1,−1}. Hence, the energy stability limit is *Ω*<1 and
*Re*_{L}=*Re*_{E} for

For the classical case of no-slip conditions at the wall, it has been proved [18] that the linear stability and the global stability limits coincide. For the double-periodic flow considered here the same is true, too. This can be proved by the same method as in [18], which amounts to selecting the Lyapunov functional in the form
*V* =*V* (**a**,*q*^{2}) might be too restrictive to obtain the same results as with (2.22).

Hence, for the flow we are considering both the energy stability limit and the actual global stability limit are known, giving the framework for considering the performance of the method [6].

## 3. Application of the method to the double-periodic rotating Couette flow

### (a) Selection of the Galerkin basis fields **u**_{i}

In solving the SOS feasibility problem (2.14), the uncertain parts of the system (2.15)–(2.16) are characterized only by their bounds. Further, those bounds are linked directly to the Galerkin basis fields **u**_{i}, *i*=1,…,*k*, as it is evident from (2.5)–(2.7). The choice of the Galerkin basis fields is therefore crucial. If *k* is too small, then the dynamic model of the system becomes over-simplified and does not adequately capture the salient features of the flow. Consequently, it might be difficult to achieve a better stability result than *Re*_{E}. If *k* is too large, the computational cost in SOS analysis will be prohibitively high. We, therefore, aim to select a limited number of Galerkin modes from amongst the eigenfunctions (2.19) of the system (2.1) in such a way that the best stability bound can be obtained.

It is shown in appendix B that for *k*=4 with modes **e**_{1,±1},**e**_{1,±2} there is no polynomial Lyapunov function for the uncertain system (2.15)–(2.16) at *Re*>*Re*_{E}. Note that in this case **f** is linear, while the proof of the existence of a polynomial Lyapunov functional given in [6] explicitly requires **f** to be quadratic, which of course can always be achieved by taking more Galerkin modes. For our particular flow, we therefore consider *k*≥6.

It is difficult to foresee the effect of mode selection on the system stability via the changes in **f**. It is also difficult to foresee the effect of mode selection on *p*(**a**,*q*^{2}). However, minimizing *κ*_{s} is clearly beneficial. According to the definition (2.8), *κ*_{s} can be minimized simply by selecting the *k* modes of the finite basis to consist of the eigenfunctions with the largest eigenvalues λ_{n,m}. Hence, we select the Galerkin modes by following such a criterion.

We define three sets of eigenfunctions {**e**_{n,m}}:
*Re*∈(*Re*_{E},*Re*_{E}+3.771), where *Re*_{E}+3.771 is the value of *Re*_{L} at *Ω*=0.1. The mode selection result for *Re*≥*Re*_{E}+3.771 can be derived similarly if needed.

### (b) Global stability of the flow

Proving that the flow is globally stable at a particular value of *Re* does not prove that it is globally stable for any other *Re*. However, for the flow in question (see §2b(ii)) there exists a global stability limit *Re*_{G} (equal to the linear stability limit *Re*_{L} for this flow) such that the flow is globally stable for all *Re*<*Re*_{G}. Hence, any *Re* for which the flow is globally stable gives a lower bound for *Re*_{G}. In particular, the SOS stability bound for the uncertain system *Re*_{SOS}≤*Re*_{G}.

The largest possible value of *Re*_{SOS} was obtained by trial end error. For a given *Re*, we try to find *V* satisfying the SOS constraints (2.14). If this is successful, we increase *Re* by *δRe* and repeat the trial. To do this, we assume some partially fixed structure of the candidate function *V*, that is the degree and the values of a part of the coefficients, and consider the remaining coefficients as the decision variables. We then tune the decision variables to satisfy (2.14), using the SOS package YALMIP [20], under which the decision variables were found using the semi-definite program (SDP) solver MOSEK [21]. Prior to solving the SDP, the SOS problem was pre-processed by the linear program solver GUROBI [22] in order to reduce and simplify the SOS program [23].

### (c) The best bound

The best result was achieved using the Galerkin mode set *K*_{1}. For *K*_{1}, (2.4) gives
**Θ**_{a}=0. Further, following the bound evaluation procedure introduced in §2a(ii), we have

Unlike the simple calculation of *κ*_{s}, estimating **Θ**_{b} and **Θ**_{c}, while not very complicated, does involve certain amount of numerical calculations, such as optimization of finite-dimensional polynomials. The validity of the bounds (3.1) and (3.2) was partially verified by solving the Navier–Stokes equations (2.1) with boundary conditions (2.18) for five sets of initial conditions for the velocity. In each of the cases, calculations were performed with *Re*=2*Re*_{L},5*Re*_{L},10*Re*_{L} and *Ω*=0.1,0.25,0.5,0.75,0.9, giving 75 combinations in total. In each case, (3.1) and (3.2) hold true. The tightness of bound evaluation for each component of **Θ**_{b} and **Θ**_{c} was verified independently by maximizing them numerically over **u**_{s} under the constraint of ∥**u**_{s}∥=1 for a set of fixed values **a**. All the tests gave satisfactory results.

The best stability bound was obtained with the following 4th-degree candidate function candidate:
*α*≥0 and *P*(**a**,**c**_{i}) is a general *i*th degree polynomial in **a**, with **c**_{i} denoting the coefficient vector. Noting that the constant term and the linear terms are obviously redundant in *V* (**a**,*q*^{2}) owing to the first SOS constraint in (2.14), they are eliminated in advance. The Lyapunov functional with the structure of (3.4) can be adjusted by tuning the decision variables **c**_{2},**c**_{4} and *α*. Since *ϵ*_{ij} required to construct the functions ℓ_{i} in (2.14) were fixed to

We then performed a trial-and-error procedure with *δRe*=0.01. Figure 1 shows the result. In the range *Ω*∈(0.2529,0.7471), *Re*_{SOS} coincides with *Re*_{L}. Hence, in that range the method [6] produces the exact result. Outside this range the method gives only the lower bound for the true global stability limit *Re*_{G}=*Re*_{L}. This lower bound is still better than the bound *Re*_{E} obtained by the standard energy stability analysis.

The obtained bound is of the form

The Lyapunov functionals obtained depend on *Ω*. For example, for *Ω*=0.1 it is
*Ω* the monomials themselves and their coefficients show some variation.

Note that the first term is proportional to the square of the total energy of the perturbation **u**. This was not imposed by us but was obtained as a result of the calculations. The likely form of the Lyapunov functional is discussed in [6].

These results demonstrate the feasibility of the method [6], thus achieving the primary goal of the present study.

### (d) Further observations

We remind that the stability of the flow at a given *Re* does not automatically imply stability at all smaller *Re*, even though such behaviour is typical. To simplify the exposition, we assume that our system does possess this typical behaviour. In fact, our results were obtained for a discrete, although reasonably dense, set of *Re* values.

The SOS stability limits obtained for the flow in question always turned out to be of the form *Re* can be used as the measure of the quality of the bound. The stability bound *Re*_{SOS,T} of the truncated system (that is the system (2.15) with **Θ**=0) is also presented in figure 1. It is
*Re*_{SOS,T}=5.52 in (3.7) to Δ*Re*_{SOS}=0.85 given by (3.5). Another immediate observation is that over a certain range of *Ω* the SOS bound for the stability margin of the truncated system is less than the true global stability margin of the full system, which coincides with *Re*_{L}(*Ω*). These differences can be due to truncation, uncertainty and/or the use of SOS relaxation, in particular, the use of a not high-enough degree of the candidate Lyapunov function. In this section, we summarize the related observations.

#### (i) Stability limits of the full system and the truncated system

We observed that a better SOS stability limit of the truncated system does not necessarily imply a better SOS stability limit of the uncertain system.

For instance, we considered another set of six Galerkin modes: *K*_{1}′={**e**_{1,0},**e**_{1,±1},**e**_{2,0},**e**_{2,±1}}, which differs from the Galerkin set *K*_{1} in that the mode **e**_{1,−2} is replaced by **e**_{2,1}. The SOS stability limit for the truncated *K*_{1}′ system was found to be *Re*_{SOS,T}=*Re*_{L}(*Ω*), while it is *K*_{1}. On the other hand, the SOS stability limits for the corresponding uncertain systems were found to be *K*_{1}′ and *K*_{1}.

This means that the selection of the Galerkin modes cannot be made on the basis of analysing the truncated system alone.

#### (ii) More Galerkin modes

Intuitively, it seems that increasing the number of Galerkin modes explicitly taken into account should improve the resulting SOS stability limit. The observations show a different picture.

For the Galerkin mode set *K*_{2}, which includes the mode set *K*_{1} and two additional modes **e**_{2,−2} and **e**_{2,1}, the uncertain fluid dynamical system is of the 8th order. In this case,
**f** is not shown here due to its large size. For a comparison with the 6-mode case, we consider the same 4th-degree Lyapunov function candidate (3.2). Owing to the increase in the number of the modes, a direct solution of the SOS optimization problem, (2.14) requires substantially greater computational effort than in the 6-mode case. For *K*_{2}, there are 532 parametric variables and 18 independent variables in the SOS optimization. Solving the SOS problem (2.14) for *K*_{2} gives the SOS stability limit
*Re*_{SOS,K2}=0.50, which is less than Δ*Re*_{SOS,K1}=0.85.

For 10 Galerkin modes set *K*_{3} no feasible Lyapunov function was found for the SOS optimization problem (2.14), even after we decreased *Re* to *Re*_{E}. In other words, Δ*Re*_{SOS,K3}=0.

To understand why increasing the number of Galerkin modes does not necessarily improve the SOS stability limit *Re*_{SOS}, we revisit the stability condition for the uncertain system:
*κ*_{s} becomes larger in magnitude, thus increasing the potential for (3.12) to be satisfied. However, using more Galerkin modes changes unfavourably the bounds of the uncertainties **Θ**_{b} and **Θ**_{c}, thus increasing *p*(**a**,*q*^{2}). As a result, the potential for (3.12) to be satisfied is reduced.

This means that there might be a finite optimum number of Galerkin modes to be included explicitly into the analysis by the method of Goulart & Chernyshenko [6].

#### (iii) Conservativeness analysis

Everywhere in this section we consider the *K*_{1} system, which gave the best SOS limits we obtained.

The global stability of the uncertain system implies the global stability of the Navier–Stokes system, but not the vice versa. The double-periodic rotating Couette flow (2.17)–(2.18) is globally stable for *Re*<*Re*_{L}(*Ω*), while the best SOS stability limit we could obtain implies global stability for *Ω* when *Re*_{SOS}<*Re*_{L}(*Ω*).

To demonstrate this, note that the obtained limits *Re*_{SOS} and *Re*_{SOS,T} are independent of *Ω* for sufficiently small or large *Ω*, as shown in figure 1. First, we consider *Re*_{SOS,T.} It is easy to show analytically that the truncated system has five steady solutions:
*Re*>*Re*_{L}. Hence, the truncated system is not globally stable when *δRe*.

For the case of the uncertain system we could not obtain an analytic solution. Increasing the degree of the Lyapunov function candidate to six by taking
**c**_{2},**c**_{4},**c**_{6} and *α* are decision variables, gave the same SOS stability limit as the 4th-degree Lyapunov function. Taking even higher degree polynomial led to so high computational costs that the calculations had to be abandoned.

Fortunately, we can demonstrate numerically that the SOS stability limit of the uncertain system obtained with the Lyapunov function of 4th degree is very close to the actual global stability limit of the uncertain system. For this reason, increasing the degree of Lyapunov function candidates is not helpful.

The idea is to evaluate the lower bound for such *Re* that there exist a steady non-zero solution of (2.15)–(2.16). This is similar to what we did for the truncated system, but we use numerical rather than analytic solutions. Naturally, the existence of steady non-zero solutions implies that the zero solution is not globally stable. Starting from any pre-specified Reynolds number for which there exists a non-zero equilibrium of the uncertain system, we decrease *Re* gradually. The smallest *Re* for which the uncertain system still has a non-zero equilibrium will be an upper bound of the global stability limit *Re*_{U} of the uncertain system.

Consider the following nonlinear optimization problem:
*F*_{i} are

The constraints *F*_{1} and *F*_{2} ensure that *F*_{3} and *F*_{4} are the bound for the uncertain terms in (2.15)–(2.16). In other words, we minimize *Re* subject to a constraint that there exist **a**,**Θ**,*Γ*,*q* satisfying a steady version of (2.15)–(2.16) and not coinciding with **a**=0,*q*=0 solution, the stability of which we are investigating. The variable *Γ* in the program can be eliminated by combining the constraints *F*_{2} and *F*_{3} as **a****Θ**−*κ*_{s}*q*^{2}≤0. Since the constraints in (3.13) include equalities and inequalities and are highly nonlinear in *Re*, instead of solving (3.13) directly we consider the following optimization problem for a given *Re*:
*ε*>0 is introduced to smooth the objective function. Note that
*ϕ*(*ω*_{1},*ω*_{2},*ϵ*) implies *ω*_{1}<0,*ω*_{2}<0. As such, all the constraints in (3.13) would be satisfied if *ϕ*(*ω*_{1},*ω*_{2},*ϵ*)<0 in (3.14). Now, the lower bound for *Re* that leads to the non-global stability of system (2.15)–(2.16) can be obtained by decreasing *Re* and solving (3.14) repeatedly. This procedure is stopped once the minimum of *ϕ*(*ω*_{1},*ω*_{2},*ϵ*) is no longer negative.

Let *ϵ*=0.1 in (3.14). Tables 2 and 3 show the optimization results for *Ω*=0 and *Ω*=0.1 when *Re* decreases from *Re*_{E}+1.00 to *Re*_{E}+0.85. The initial searching point in each trial is always set as the stopping point in the previous trial. The sequential quadratic programming method [24] associated with the function NLPSolve in MAPLE optimization toolbox is used to solve (3.14). From the tables, we can see that all the constraints in (3.13) can be satisfied when *Re*≥*Re*_{E}+0.86 in the sense that the residual error |**f**+**Θ**|^{2} is negligible. This implies that the global stability limit for the uncertain system is less than *Re*_{E}+0.86. Recalling that *Re*_{SOS}=*Re*_{E}+0.85 for *Ω*=0 and *Ω*=0.1, we conclude that at least for these values of *Ω* the actual global stability limit for the uncertain system is between *Re*_{SOS} and *Re*_{SOS}+0.01. The verification for other values of *Ω* can be conducted similarly.

Overall, the analysis of this section shows that for the flow in question further improvement of the SOS stability limit can be achieved only by improving the uncertainty bounds, since increasing the number of Galerkin modes explicitly taken into account gives no improvement, while the SOS stability limit for the uncertain system is already tight.

## 4. Discussion and conclusion

The uncertainty bounds obtained by the methods proposed in [6] and appendix A are tight for each of the components of **Θ**_{b} and **Θ**_{c}. However, they are attained at different **a** and **u**_{s}. Hence, the bound on |**Θ**| used to obtain the stability limit is not tight. This provides a potential for improving the results.

For all practical purposes, the question of a global stability of a particular flow can often be answered by physical experiment or numerical calculation, at least up to the existence of unstable stationary points or orbits. Theoretical studies of global stability, however, can provide deeper insight into various aspects of the flow and of the methods used. For example, one could hope to gain such insight from the particular form of the Lyapunov functional. It appears, however, that the method [6] generates Lyapunov functionals of a rather complicated form, which is difficult to interpret. In addition, the obtained Lyapunov functional depends on computational parameters such as the number of the Galerkin modes taken explicitly into account and the degree and the form of the polynomial representing the candidate Lyapunov function.

Reduction of the Navier–Stokes equation to an uncertain system, which forms the basis of the method of [6], turns out to be useful for a number of problems, which might be even of larger interest than the global stability problem. For example, it might prove useful in studies of bounds for long-time averages of the characteristics of turbulent flows and other infinite-dimensional systems with complicated behaviour, and in designing control for such systems [14]. For such studies, the observations of the properties of the method of [6] made in the present work might constitute even a greater interest than the central result of this paper.

It remains to give a summary of the obtained results and observations.

A new expression for one of the uncertainty bounds required by the method of [6] was obtained (see appendix A). It allows a considerable reduction of the computational cost as compared to the approach proposed in [6].

A systematic approach to selecting the particular Galerkin modes for the uncertain system was proposed.

It was shown that the selection of Galerkin modes leading to a better stability result for the truncated Galerkin system does not necessarily lead to a better stability result for the uncertain system.

It was also shown that increasing the dimension of the uncertain system does not necessarily improve the stability bound obtained for the full system, and that increasing the degree of the candidate polynomial Lyapunov function also does not necessarily improve the bounds. This suggests that further progress in this problem is more likely to be achieved by improving the bounds in the uncertain system than by improving SOS optimization.

For the particular version of the double-periodic rotating Couette flow we considered, we demonstrated that

(i) for

*Ω*∈(0.2529,0.7471), the SOS stability limit coincides with the actual global stability limit.(ii) for

*Ω*∈(0,0.2529]∪[0.7491,1), the SOS stability limit*Re*_{SOS}≥*Re*_{E}+0.85,

where *Re*_{E} is the energy stability limit. This demonstrates the feasibility of the method proposed in [6].

To the best of our knowledge, this is the first case in which a systematic method applicable in principle to any fluid flow was used successfully to prove global stability of a fluid flow for the value of the Reynolds number greater than that which could be achieved with the energy stability approach.

## Authors' contributions

S.C. and P.G. conceived the idea of the paper. S.C. performed the reduction to the uncertain system for the flow in question using Mathematica package and D.H. did the same using Maple, thus confirming independently the correctness of the uncertain system. D.H. performed all the sum-of-squares calculations, thus obtaining the main results. D.L. and O.T. did all the Navier–Stokes calculations. P.G. gave advice to D.H. on the use of YALMIP. S.C. conceived the idea of appendix A, and F.F. implemented it in the form presented, and wrote the appendix. All authors regularly discussed the progress during the entire work, which was truly collaborative, so that each author contributed to some extent to each of the results. D.H. and S.C. wrote the first draft of the paper. All authors took part in bringing the draft to the final form.

## Competing interests

We have no competing interests.

## Funding

This work was supported by the UK Engineering and Physical Sciences Research Council grants EP/J011126/1, EP/J010537/1 and EP/J010073/1 and received support-in-kind from Airbus Operation Ltd, ETH Zurich (Automatic Control Laboratory), University of Michigan (Department of Mathematics) and University of California, Santa Barbara (Department of Mechanical Engineering). D.H. was partially supported by the State Key Laboratory of Industrial Control Technology grant (ICT1304), Zhejiang University, China.

## Appendix A. Bound evaluation for **Θ**_{c}(**u**_{s})

Note that **u**_{s}⋅∇**u**_{i}⋅**u**_{s}=**u**_{s}⋅*D*_{i}(**x**)⋅**u**_{s}, where *D*_{i}(**x**)=(∇**u**_{i}+∇^{T} **u**_{i})/2 are the rate of strain tensors of the Galerkin basis vector fields. Hence, for *i*=1,…,*k*,
*D*_{i}(**x**), it is immediate that *D*_{i}(**x**), which as usual is denoted by *ρ*(*D*_{i}(**x**)). As such, **u**_{s} with compact support similar to a delta function centred around the point **x**_{m}, where *ρ*(*D*_{i}(**x**)) attains its global maximum. The idea here is similar to the one used in [25], while the difference is that **u**_{s} is a vector and not a scalar, and that it has to be solenoidal, and orthogonal to all the basis field functions **u**_{i},*i*=1,…,*k*.

Now, by (A3), this yields the tight bound

## Appendix B. Non-existence of polynomial Lyapunov functions

An example is given to illustrate that an excessively low-order uncertain system may not be suitable in the sense that there does not exist a polynomial Lyapunov function that guarantees its stability when *Re*>*Re*_{E}.

In this case, the uncertain system (2.15)–(2.16) involves only four modes. More precisely, let **u**_{1}=**e**_{1,1},**u**_{2}=**e**_{1,−1},**u**_{3}=**e**_{1,2} and **u**_{4}=**e**_{1,−2}. It is easy to calculate that
**Θ**(**u**_{s},**a**) and *Γ*(**u**_{s}) are, respectively,
*V* =*V* (**a**,*q*^{2}) satisfying ∂*V* /∂*q*^{2}≥0, a sufficient condition for stability is (2.11), i.e.
*Re*<*Re*_{E}, the total energy *E*_{t}=|**a**|^{2}/2+*q*^{2} is the Lyapunov function as |∂*E*_{t}/∂**a**−∂*E*_{t}/∂(*q*^{2})**a**|=0 and (∂*E*_{t}/∂**a**)**f**+∂*E*_{t}/∂(*q*^{2})⋅*κ*_{s}*q*^{2}<0. However, this is not the case for *Re*>*Re*_{E}, namely, there does not exist a polynomial Lyapunov function for the uncertain system (2.15)–(2.16), if only the modes **e**_{1,±1},**e**_{1,±2} are considered.

We first prove that the total energy *E*_{t} is not a Lyapunov function. Note that *Re* increases and exceeds *Re*_{E}, λ_{1,−1} will be the first to become positive, implying that *E*_{t} is not a Lyapunov function.

Moreover, as *V* cannot be proportional to *E*_{t}. Indeed, let **f**(**a**)=*O*(|**a**|), *a*_{2} and *a*_{1}=*a*_{3}=*a*_{4}=*q*=0, (∂*V* /∂**a**)**f**>0, so that *V* is not a Lyapunov function.

On the other hand, in view of (B1), *p*(**a**,*q*^{2})=*O*(*q*^{2}(|**a**|^{2}+*q*^{2})). If the main term of *V* is not a function of *E*_{t}, then the order of the right-hand side of (B2) is greater than the order of its left-hand side, and the inequality (B2) cannot hold true, so that *V* cannot be a Lyapunov function in this case either.

Overall, when *Re*>*Re*_{E}, there does not exist a polynomial Lyapunov function for the uncertain system (2.15)–(2.16), if only the modes **e**_{1,±1},**e**_{1,±2} are considered.

## Footnotes

- Received September 3, 2015.
- Accepted October 26, 2015.

- © 2015 The Authors.