Stability analysis of fluid flows using Lagrangian Perturbation Theory (LPT): application to the plane Couette flow

We present a new application of Lagrangian Perturbation Theory (LPT): the stability analysis of fluid flows. As a test case that demonstrates the framework we focus on the plane Couette flow. The incompressible Navier-Stokes equation is recast such that the particle position is the fundamental variable, expressed as a function of Lagrangian coordinates. The displacement due to the steady state flow is taken to be the zeroth order solution and the position is formally expanded in terms of a small parameter (generally, the strength of the initial perturbation). The resulting hierarchy of equations is solved analytically at first order. We find that we recover the standard result in the Eulerian frame: the plane Couette flow is asymptotically stable for all Reynolds numbers. However, it is also well established that experiments contradict this prediction. In the Eulerian picture, one of the proposed explanations is the phenomenon of `transient growth' which is related to the non-normal nature of the linear stability operator. The first order solution in the Lagrangian frame also shows this feature, albeit qualitatively. As a first step, and for the purposes of analytic manipulation, we consider only linear stability of 2D perturbations but the framework presented is general and can be extended to higher orders, other flows and/or 3D perturbations.


I. INTRODUCTION
Understanding the transition of fluid flows from the stable to turbulent regime is one of the central questions in the studies of turbulence. In usual linear stability analysis one formally expands the Navier-Stokes equation about a steady state flow assuming that the perturbations to the background flow are small. Study of linear stability of laminar flows, such as the plane Couette flow, commenced over a hundred years ago with the seminal work by Orr [1] and Sommerfeld [2]. The resulting Orr-Sommerfeld equation has been studied extensively over the decades (see Bayly et al. [3] for a review). Purely analytic investigations were not conclusive (e.g., Joseph [4]) and many efforts were devoted to obtaining a numerical solution of this equation. The basic idea was the same: expand the velocity in terms of an orthogonal set of basis functions and recast the system as an eigenvalue problem. But the analyses differed in their choices of basis functions, which resulted in different convergence rates. For example, Gallagher and Mercer [5] used the Chandrasekhar-Reid functions [6], Grosch and Salwen [7] used a modified version of the same and Orszag [8] used Chebyshev polynomials, which have the best convergence rate. The numerical results indicated that the plane Couette flow is linearly stable for all Reynolds numbers whereas the plane Poiseuille flow exhibits a transition to turbulence at Re = 5772.22 [8]. Further insight into these results was obtained by other combined numerical and analytical techniques [9][10][11].
However, experimental results show a discrepancy with linear predictions. For example, the plane Couette flow shows a transition to turbulence when none is expected while the plane Poiseuille transitions to turbulence at a Reynolds number much lower than the linear estimate e.g., [12][13][14]. Various approaches have been employed to explain this transition. One way is to computationally investigate the full non-linear Navier-Stokes equation as was done by Orszag and collaborators [15,16] who showed that the energy growth in the system corresponds to a sub-critical bifurcation. Another is to look for finite amplitude equilibrium states near the transition and examine their stability against two and three dimensional perturbations [17]. A third approach is to understand the stability properties of the perturbed base flow [18,19]. Some authors have also explored mechanism of transient growth due non-normality of the linear operator [20][21][22]; see Grossmann [23] and Schmid [24] for recent reviews.
Majority of the analytical stability analysis has been carried out in the Eulerian frame. In this frame the velocity is the fundamental variable and is expressed as a function of a fixed Eulerian coordinate system (gird coordinates). On the other hand in the Lagrangian frame, the particle position is fundamental and is expressed as a function of a Lagrangian coordinate (usually the initial position) and time. Eulerian measurements are easier, whereas Lagrangian methods require sophisticated 3D particle tracking techniques e.g., [25,26]. The choice of frame depends on the problem at hand. Much analytic and numerical work has been done in the Lagrangian frame in terms of analyzing particle statistics, predicting scaling laws, structure of correlation functions etc. However, linear stability analysis in this frame has been relatively rare.
In this paper, we examine the stability of laminar flows using a perturbative scheme in the Lagrangian frame i.e., Lagrangian Perturbation Theory (LPT). We focus only on 2D perturbations of the incompressible plane Couette flow but the formalism is general and can be extended to other flows. The main motivation to use this method is that the Lagrangian derivative includes by definition the non-linear term (v · ∇)v and hence it ought to be able to better estimate the non-linear effect. Furthermore, a flow which is unstable in the Eulerian frame is also unstable in the Lagrangian frame. This means that stability in the Lagrangian frame implies stability in the Eulerian frame. Thus, the investigation of Lagrangian stability can provide a independent confirmation of Eulerian stability.
The main drawback of this scheme is that relies on the one-to-oneness of the map between the Eulerian and the Lagrangian frame and fails when particles cross. Thus, it is not expected to model the turbulent regime, but it could potentially be useful to model the transition to turbulence. LPT has been employed in the past in astrophysics and cosmology (for e.g., [27,28]) to study the non-linear evolution of gravitating matter in an expanding universe. In recent work, Nadkarni-Ghosh and Chernoff [29,30] extended this theory by introducing the new feature LPT re-expansions. They showed that by an iterative application of the perturbative scheme, it was possible to achieve convergence to the exact non-linear answer, prior to orbit crossing. Multiple steps with even a first order linear scheme can successfully model the non-linear regime. In this paper we present the first order solution of the scheme for a single step. In doing this, we are able to analyze linear stability, but we do not achieve the full non-linear solution.
To the best of our knowledge a perturbative analysis in the Lagrangian frame has only been performed in the past by Pierson [31] in the context of geophysical flows. But our work differs from theirs because the order counting and flow geometries are different giving rise to a different set of equations and solutions. They use the Lagrangian particle labels as the zeroth order solution whereas we use the displacement due to the base flow as the zeroth order solution. The latter approach allows one to more easily track the dependence of the base flow making it easier to generalize.
The paper is organized as follows §II sets up the equations in the Lagrangian frame. §III outlines the perturbative solution. The full solution needs to be computed numerically, but we show analytically that at late times the perturbations decay, confirming linear stability. §IV discusses the procedure to recover the Eulerian velocity from the Lagrangian velocity. §V provides a discussion and summary.

II. EQUATIONS IN THE LAGRANGIAN FRAME
The incompressible fluid is described by the dimensionless system of equations where d/dt is the usual convective derivative d dt = ∂ ∂t + (v · ∇)v and P , ρ and Re denote the pressure, density and Reynolds number respectively. In this paper we will focus only on laminar flows, in particular on the plane Couette flow which consists of two parallel plates moving with respect to each other with a steady state velocity v s.s. . The velocity profile is invariant along the flow direction (defined to be the x-axis) and varies only in the direction perpendicular to the flow (defined to be the y-axis). Let r = {x, y, z} denote the physical position of a fluid element. We will restrict to 2-d perturbations and hence the fluid displacements are confined only to the x − y plane. The system given by equations (1) and (2) is Eulerian: the velocity is the fundamental quantity and is usually solved in terms of the fixed coordinate system (grid coordinates) i.e. v = v E (r), where the subscript 'E' denotes Eulerian. On the other hand, in the Lagrangian framework, the particle position is the fundamental quantity and is usually solved in terms of some fixed Lagrangian coordinate and time. We choose the Lagrangian coordinate R = {X, Y, Z} to be the physical position at the initial time (t = 0), The position at any later time is and the corresponding Lagrangian velocity is where the 'dot' denotes the total time derivative or Lagrangian derivative d/dt (this is also sometimes denoted in the literature as D/Dt) and the subscript 'L' denotes Lagrangian.
Making r the fundamental variable, the system of equations (1) and (2) can be equivalently expressed as where we have replaced equation (1) with its curl. The spatial derivatives in the above equation are with respect to the physical variable r. These have to be transformed to derivatives with respect to the Lagrangian coordinate (see Appendix VII A). This yields the system where the commas denote derivatives with respect to the Lagrangian spatial coordinate and is the determinant of the Jacobian of the transformation between the Eulerian and Lagrangian coordinates. Note that the time derivative commutes with the spatial Lagrangian differentiation i.e., 'dots' and 'commas' commute.

III. LAGRANGIAN PERTURBATION THEORY (LPT)
In the Lagrangian perturbative scheme, r is formally expanded as where p (n) is the n-th order term and ǫ is just a formal book-keeping parameter. The definition of the Lagrangian coordinate (equation (3)) is used to set the initial values of the displacement vectors at each order. We choose A. Zeroth order The background steady state solution for the plane Couette flow is given by v s.s. E = {cy, 0, 0}. By definition, y = Y at the initial time. Then the initial velocity in the Lagrangian frame is v L (R, 0) = {cY, 0, 0}. Thus the particle at initial position R at t = 0 is at R + v L (R, 0)t after time t. We take this to be the zeroth order solution for the position vector i.e.
In the Eulerian framework, stability of the flow is examined by perturbing around the steady state velocity in Eulerian coordinates. The zeroth order solution in that case is always an exact solution of the incompressible Navier-Stokes system. It is necessary to check that the transformation to Lagrangian coordinates preserves this property of the background solution. This is shown in Appendix VII B.

B. First order equations
Substitute the ansatz of equation (11) in the system of equations (8) and (9) and keep terms upto first order. Using the zeroth order solution given by equation (14) (see Appendix VII C for details), the equation for the first order displacement is In the usual Eulerian perturbation theory the perturbed velocity also remains divergence-free at first (and higher) orders. However, the divergence in the Lagrangian frame is not zero at any order. This is not surprising since the two coordinate systems are distinct and the non-zero divergence arises from the transformation between them.
In order to satisfy equation (15), we assume p (1) to have the form where, ψ is a scalar function of Lagrangian coordinates and time i.e., ψ ≡ ψ(X, Y, t). This giveṡ ψ is analogous to a stream-function, but not the same as one encountered in the usual Orr-Sommerfeld analysis. Substituting the form in equation (17) in equation (16) and simplifying gives This system can be recast as where the operator A The solution is obtained by first solving equation (21) for φ and then solving equation (20) forψ. Integrate to get ψ.

C. Initial and boundary conditions
The symmetry of the underlying flow implies periodic boundary conditions along the X-axis for solving both φ anḋ ψ. We assume that the X-dependent part of the solution can be separated from the rest and represent it by a Fourier series expansion. With this ansatz, the net system of equations for the stream function ψ is second order in time and fourth order in spatial variable Y . Accordingly two temporal initial conditions and four boundary conditions are needed.
The perturbation velocity profile at t = 0 is specified initially. In numerical simulations this is done by exciting specific modes or specifying an initial power spectrum of the perturbation. Let v 0 (X, Y ) formally denote this initial perturbation.ṗ (1) The definition of the Lagrangian coordinate provides the other initial condition. Equation (13) for n = 1 is The boundary conditions are more involved. The symmetry of the underlying flow implies periodic boundary conditions along the X-axis for solving both φ andψ. The no slip condition imposed on the wall at Y = 0 meanṡ These are two conditions corresponding to the X and Y coordinate. Note thatṗ depends explicitly on ψ andψ but only indirectly on φ. For simplicity we will assume that the flow is semi-bounded i.e., wall is placed only at Y = 0. This allows us to relate equation (25) to conditions onψ. We assume a Fourier decomposition along the Y axis for φ.
Using the definition ofṗ (1) from equation (18) in equation (25) gives Evaluating the Y -coordinate of equation (24) at Y = 0 using equation (17) (26) is zero at all times and it follows that

D. First order solution
We now solve for φ and ψ subject to the above conditions. Equation (21) for φ has spatial and temporal derivatives appearing on different sides of the equation and hence is separable. Following the arguments in the earlier section, the ansatz for φ is Substituting in equation (21), gives The solution is where k 2 = k 2 x + k 2 y and f (0) is the integration constant to be set later. The solution for φ is One can now solve equation (20) forψ. From the structure of the equation one can assume that the purely temporal functions are the same for both φ andψ. Any additional time dependence introduced via the operator A is necessarily also a function of Y . This gives the ansatż The function G(Y, t) is parametrized by k y . For convenience of notation we will henceforth drop the subscript. Substituting in equation (20) and using equation (32) gives where the primes denote differentiation w.r.t Y . The solution for G(Y, t) can be split into a homogenous part and a particular solution: and where Satisfying the boundary conditions represented by equations (27) and (28) fixes C 1 (t) and C 2 (t).
The net solution for g is ψ is obtained by integratingψ w.r.t. time: where C(X, Y ) is the constant of integration. This can be re-written as where h(Y, t) = g(Y, t)f (t)/f (0)dt. The initial condition equation (24) is satisfied if ψ = 0 at t = 0. This fixes The only quantities remaining to be determined are f (0) andφ(k x , k y ). This is set by the initial velocity perturbation. The individual components φ(k x , k y ) and characterize its shape and f (0) sets the overall initial amplitude. Late time behaviour: The above prescription completely specifies the initial conditions and the numerical solution for ψ can be computed at any later time. Simplified analytic expressions can be obtained at late times. Note that the terms arising from the homogenous part of the solution are damped and oscillatory. If the integration is over a sufficiently large time interval t >> (ck x ) −1 they integrate to zero leaving Substituting for C 3 and f from equations (31) and (37) gives where γ(s, z) = z o e − z ′ z ′s−1 dz ′ is the lower incomplete gamma function. Thus ψ and hence the velocity perturbation evolves as ∼ e − k 2 x c 2 t 3 3Re at late times and the flow is linearly stable in the Lagrangian frame. This late time behavior is in exact agreement with that in Press & Marcus [11] which was obtained using symmetry arguments for the unbounded Couette flow. Once ψ is known, the full first order displacement and velocity can be computed from equations (17) and (18).
In the above analysis we used Fourier basis functions for both the X and Y axis in solving for φ. This allowed us to solve for φ and ψ in succession. If the flow is bounded at some finite Y then the boundary conditions are more complicated. φ and ψ then get coupled through the boundary conditions and it is not possible to obtain analytic expressions. Alternatively, one may choose another set of basis functions which will satisfy the correct boundary conditions.

IV. TRANSFORMING BACK TO THE EULERIAN FRAME
The perturbative scheme and the LPT solution outlined in the previous section solve equations (6) and (7) for the variable r. However, the original system whose solution we seek is given by equations (1) and (2). Equation (6) is obtained by taking spatial derivatives of equation (1) and hence the LPT solution is insensitive to any spatially homogenous time dependent transformation ∆r(t). This distinction leads to the definition of two separate frames of reference: the physical frame and the calculational frame. Let r phys denote the solution in the physical frame that satisfies the original set and r LP T denote the solution in the calculational frame obtained by the perturbative treatment discussed in the previous section. The two are related as r = r LP T + ∆r(t). (46) By substituting equation (46) in equation (1) (with v = dr/dt), one obtains a differential equation for ∆r(t), where the source terms are determined by the LPT solution. The initial conditions for ∆r are specified by the transformation between the physical frame and the calculational frame at the initial time. In the simple case of the plane Couette flow, we can assume ∆r(0) = 0 and ∆ṙ(0) = 0. In recent work, Nadkarni-Ghosh & Chernoff [30] showed that convergence properties of the perturbative solution crucially depend on fixing this degree of freedom, although the work was in the context of a different physical system, namely dark matter fluid gravitating an an expanding universe. Since our aim in this paper is to merely examine the stability, we do not explicitly calculate the exact form of ∆r(t). The net physical solution and velocity Note however that the v is known as a function of the Lagrangian coordinate. In order to obtain the Eulerian velocity v E (r, t) one has to solve for the initial R of the fluid element which is located at r at time t i.e. if r = F (R, t) then the Eulerian velocity at the coordinate point r is

V. DISCUSSION AND CONCLUSION
In this paper we have performed a stability analysis of the plane Couette flow in the Lagrangian frame. The zeroth order solution is taken to be the displacement due to the steady state Couette flow and perturbations to this solution are examined at linear order. The analysis presented here differs from the usual (Eulerian) linear stability analysis in one fundamental aspect. In the latter, terms such as (v · ∇)v, which are second order in the velocity perturbation are ignored and the equations are recast as a linear eigenvalue problem ∂v/∂t = Sv, where the matrix S depends purely on spatial terms and its eigenvalues determine the stability of the system. The resulting solution is exponential, the exponent is always linear in the time coordinate t. For laminar flows such as the shear-driven plane Couette flow and the pressure-driven plane Poiseuille flow in a pipe, the eigenvalues are always negative and the flows are linearly stable. However, it is observed in experiments and simulations that these flows do undergo a transition to turbulence. One mechanism proposed to explain this is the 'non-normal' growth of transients due to a coupling to the non-linear term (see the review by Grossmann [23]). We note that such 'transient' growth is also observed in the linear Lagrangian solution. The form of f (t) in equation (31) has a term that grows as e kxkyct 2 /Re before eventually falling off as e −k 2 x c 2 t 3 /3Re at late times. For sufficiently large k y , this term can dominate the dynamics. This growth may potentially be attributed to the fact that the (v · ∇)v term is not entirely left out of the Lagrangian analysis. This makes linear Lagrangian perturbation theory more powerful in some ways than the standard linear Eulerian perturbation theory.
We have focussed on the simplest shear flow: the plane Couette flow and restricted to 2D perturbations. This greatly simplified the expressions (see Appendix VII C). For other flows or 3D perturbations the form of the equations will be more involved, but the framework remains the same. It is also possible to extend to the perturbative formalism to the non-linear regime by keeping terms to higher order in the displacement field in equation (11). Alternatively, it is possible to model the non-linear regime by repeated expansions of the linear PT (Nadkarni-Ghosh and Chernoff [29,30]). This technique was initially developed in order to overcome the fact that, independent of orbit crossing, the Lagrangian series solution has a finite time range of validity. The basic idea is that LPT can be thought of as a numerical finite difference scheme with an associated time stepping criterion. It may be possible to apply similar ideas to model the non-linear solution of the Navier-Stokes equation, however, issues of convergence, numerical stability etc. will need to be considered carefully to get meaningful results.

A. Mathematical Transformations
We start with the expression for divergence of the velocity.
Einstein's repeated summation convention is followed. The inverse transformation from R-space to r-space is given as where and ǫ jlm is the usual Levi-Civita symbol (Arfken & Weber [32]). The incompressibility condition reduces to where commas denote spatial derivatives with respect to the Lagrangian coordinate. Note that this can also be written as Consider the curl of the Navier-Stokes equation. The i-th component of the l.h.s. is where the last equality follows from equation (51). The r.h.s. of the Navier-Stokes is proportional to ∇ 2 r (∇ r ×ṙ). For any scalar f i , ∇ r converted to Lagrangian coordinates is Using equation (51) gives Substituting f i = ǫ ijk ∂ṙ k ∂rj and again using equation (51) to transform derivatives gives Thus the curl of the Navier-Stokes equation in Lagrangian coordinates reduces to ǫ ijk ǫ jmn ǫ lm ′ n ′r k,l r m,m ′ r n,n ′ = 1 Re ǫ ijk ǫ jf g ǫ mf ′ g ′ ǫ lde ǫ qd ′ e ′ ǫ lab ǫ pa ′ b ′ r d,d ′ r e,e ′ 1 2J r a,a ′ r b,b ′ 1 2Jṙ k,m r f,f ′ r g,g ′ We check that this is an exact solution of the incompressible Navier-Stokes system given by equations (6) and (7). The transformation between the Lagrangian and Eulerian frame is Here i is the row-wise index, j is column-wise. To check if equation (6) and (7), it suffices only to consider the x-component since the others are trivially zero. The divergence-less condition given by equation (7) is In equation (6), the l.h.s. is zero since there is no t dependence inṙ and Lagrangian derivative is just the total time derivative acting onṙ. So it remains to prove that r.h.s =0. The x-component is Applying the change of derivatives again and using the fact that ∂Y /∂x = 0, ∂ṙ x /∂X = 0 and ∂ 2ṙ x /∂Y 2 = 0, all terms become zero. Thus both Navier-Stokes and the incompressibility are satisfied when the base flow is expressed in Lagrangian coordinates. It is a natural candidate for the zeroth order particle position and one can examine the stability of the system by perturbing about this steady state solution.