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.

1. 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 [4–6] 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. 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 [7]. Further insight into these results was obtained by other combined numerical and analytical techniques [8–10].

However experimental results particularly for shear driven flows, 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 [11–13]. 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 [14, 15] 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 [16]. A third approach is to understand the stability properties of the perturbed base flow [17, 18]. Further developments in the 1990s showed that the instability can be attributed to the “non-normality” of the linear stability operator [19, 20]. The eigenvectors of the linear operator are not orthogonal and this allows for the possibility of transient growth before the eventual asymptotic decay implied by the negative eigenvalues. The non-linear term can then amplify this growth [21–23]; see Grossmann [24] and Schmid [25] 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 (grid 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 (for example La Porta et al. [26, 27]). 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 (see for example [28] and references therein). In the context of geophysical flows, the Lagrangian picture has been used extensively to understand the backreaction effect of non-linear perturbations on the mean background flow following the formulation by Andrews and McIntyre [29]. Formal work regarding mathematical properties of the Lagrangian trajectories has been also performed [30] recently. 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). As a simple test case, we focus only on 2D perturbations of the incompressible plane Couette flow but the formalism is general and can be extended to other flows. One of the motivations 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. Thus, the investigation of Lagrangian stability can provide an 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. It is not expected to model the turbulent regime where orbit crossing is likely to occur.

LPT has been used in other branches of physics, most notably in cosmology, to model the growth of non-linear structure in the universe. The statistical theories of homogeneous and isotropic turbulence and the growth of non-linear cosmological large scale structure share many common features. The velocity field in the case of turbulence and the density field in the case of cosmology are both modeled as random fluctuations in a homogenous and isotropic background. In turbulence, the convective term in Navier-Stokes is the main source of non-linearity. Higher order velocity correlation functions are the main quantities of interest and one is interested in their scaling properties. In cosmology, the non-linearity arises both from the convective term in the Euler equation and gravity. Density correlations are of importance and they are used to constrain cosmological parameters. In the past, many perturbative techniques from the theory of turbulence have been successfully applied to analytically model cosmological structure. For example: Taruya and Hiramatsu [31] use the Direct Interaction Approximation method of Kraichnan [32–34] to address the “closure problem” of the hierarchy of moment equations. Crocce and Scoccimarro [35] use the techniques discussed in Wyld [36] and L'vov and Proccacia [37] to formulate a renormalized perturbation theory. The adhesion approximation [38], used to deal with Lagrangian particle crossings is essentially the model of 3D Burgers turbulence, see work by Frisch and collaborators [39, 40], Gaite [41].

We attempt to do the reverse: use a technique from the theory of large scale structure to understand the transition to turbulence. In this paper, we present the first step in applying LPT to the analysis of the incompressible Navier-Stokes equation. For simplicity we restrict to 2D perturbations and compute the first order solution of the scheme. The linear analysis using LPT analytically confirms the linear Eulerian stability result that the plane Couette flow is asymptotically stable at all Reynolds numbers. In addition it recovers the feature of transient growth, which in the Eulerian case is attributed to the non-normality of the linear stability operator. To the best of our knowledge a perturbative analysis in the Lagrangian frame has been performed in the past by Pierson [42] in the context of geophysical flows. But our work differs from Pierson's because the order counting and flow geometries are different giving rise to a different set of equations and solutions. Pierson uses 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 section 2 sets up the equations in the Lagrangian frame. Section 3 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. Section 4 discusses the procedure to recover the Eulerian velocity from the Lagrangian velocity. Section 5 provides a discussion and summary.

2. Equations in the Lagrangian Frame

The incompressible fluid is described by the system of equations:

$\begin{array}{ll}\frac{d\text{v}}{dt}=-\frac{\nabla P}{\rho }+\nu {\nabla }^2\text{v},& (1)\end{array}$

$\begin{array}{ll}\nabla ·\text{v}=0,& (2)\end{array}$

where d/dt is the usual convective derivative $\frac{d}{dt}=\frac{\partial }{\partial t}+(\text{v}·\nabla )$ and P, ρ and ν denote the pressure, density and kinematic viscosity 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. (see Figure 1). 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 2D perturbations and hence the fluid displacements are confined only to the xy-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). Thus,

$\begin{array}{ll}\text{R}=\text{r}(\text{R},0).& (3)\end{array}$

FIGURE 1

Figure 1. Schematic representation of the 2D plane Couette flow. Two parallel plates are a distance h apart; the steady laminar flow has a linear profile w.r.t. the y coordinate i.e., v_s.s. = {cy, 0, 0}. For the semi-bounded case, h → ∞.

The physical position at any later time is:

$\begin{array}{ll}\text{r}\equiv \text{r}(\text{R},t)& (4)\end{array}$

and the corresponding Lagrangian velocity is:

$\begin{array}{ll}{\text{v}}_L=\dot{\text{r}}(\text{R},t),& (5)\end{array}$

where the subscript “L” denotes Lagrangian and “dot” denotes the time derivative d/dt. Note that in the Lagrangian frame the spatial variable R does not change with time. Thus, d/dt is not represented in terms of a partial time derivative and the convective term as is done in the Eulerian frame; instead it is the total derivative and is also denoted as D/Dt in the literature.

Taking the curl of Equation (1) and making r the fundamental variable, the system of Equations (1) and (2) can be equivalently expressed as:

$\begin{array}{ll}{\nabla }_r\times \ddot{\text{r}}=\nu {\nabla }_r^2({\nabla }_r\times \dot{\text{r}}),& (6)\end{array}$

$\begin{array}{ll}{\nabla }_r·\dot{\text{r}}=0.& (7)\end{array}$

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 6.1). This yields the system:

$\begin{array}{l}{ϵ}_{ijk}{ϵ}_{jmn}{ϵ}_{l{m}^{\prime }{n}^{\prime }}{\ddot{r}}_{k,l}{r}_{m,m^{\prime }}{r}_{n,n^{\prime }}=\nu {ϵ}_{ijk}{ϵ}_{jfg}{ϵ}_{m{f}^{\prime }{g}^{\prime }}{ϵ}_{lde}{ϵ}_{q{d}^{\prime }{e}^{\prime }}{ϵ}_{lab}{ϵ}_{p{a}^{\prime }{b}^{\prime }}{r}_{d,d^{\prime }}{r}_{e,e^{\prime }}\\ {\left\{\frac{1}{2J}{r}_{a,a^{\prime }}{r}_{b,b^{\prime }}{(\frac{1}{2J}{ṙ}_{k,m}{r}_{f,f^{\prime }}{r}_{g,g^{\prime }})}_{,p}\right\}}_{,q}& (8)\end{array}$

$\begin{array}{ll}{ϵ}_{ilm}{ϵ}_{j{l}^{\prime }{m}^{\prime }}{ṙ}_{i,j}{r}_{l,l^{\prime }}{r}_{m,m^{\prime }}=0,& (9)\end{array}$

where the commas denote derivatives with respect to the Lagrangian spatial coordinate and

$\begin{array}{ll}J=\text{Det}\left(\frac{\partial r_i}{\partial R_j}\right)=\frac{1}{6}{ϵ}_{ijk}{ϵ}_{lmn}{r}_{i,l}{r}_{j,m}{r}_{k,n}& (10)\end{array}$

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.

3. Lagrangian Perturbation Theory (LPT)

In the Lagrangian perturbative scheme, r is formally expanded as:

$\begin{array}{ll}\text{r}(\text{R},t)={\displaystyle \sum _{n=0}^{\infty }}{\text{p}}^{(n)}(\text{R},t){\Delta }^n& (11)\end{array}$

where p⁽ⁿ⁾ is the n-th order term and Δ is just a formal book-keeping parameter related to the strength of the initial velocity perturbation. The definition of the Lagrangian coordinate (Equation (3)) is used to set the initial values of the displacement vectors at each order. We choose:

$\begin{array}{ll}{\text{p}}^{(0)}(\text{R},0)=\text{R},& (12)\end{array}$

$\begin{array}{ll}{\text{p}}^{(n)}(\text{R},0)=0\forall n>0.& (13)\end{array}$

3.1. Zeroth Order

The background steady state solution for the plane Couette flow is given by ${\text{v}}_E^{s.s.}=\left\{cy,0,0\right\}$. 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.,

$\begin{array}{ll}{\text{p}}^{(0)}(\text{R},t)=\text{R}+{\text{v}}_L(\text{R},0)t.& (14)\end{array}$

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 6.2.

3.2. First Order Equations

Substitute the ansatz of Equation (11) in the system of Equations (8) and (9) and keep terms up to first order. Using the zeroth order solution given by Equation (14) (see Appendix 6.3 for details), the equation for the first order displacement is:

$\begin{array}{ll}{\nabla }_R·{\dot{\text{p}}}^{(1)}=c\left(p_{Y,X}^{(1)}+t{ṗ}_{Y,X}^{(1)}\right)& (15)\end{array}$

$\begin{array}{l}{\ddot{p}}_{Y,X}^{(1)}-{\ddot{p}}_{X,Y}^{(1)}+ct{\ddot{p}}_{X,X}^{(1)}=\nu \left({\nabla }_R^2-2ct\frac{{\partial }^2}{\partial X\partial Y}+c^2{t}^2\frac{\partial }{\partial X^2}\right)\\ \left({ṗ}_{Y,X}^{(1)}-{ṗ}_{X,Y}^{(1)}+ct{ṗ}_{X,X}^{(1)}-c{p}_{X,X}^{(1)}\right).& (16)\end{array}$

Here $p_{Y,X}^{(1)}$ denotes the derivative of the Y component of p⁽¹⁾ w.r.t. X (partial spatial derivative). In the usual Eulerian perturbation theory the perturbed velocity also remains divergence-free at first (and higher) orders; this condition arises because the flow is incompressible. This condition translated into the Lagrangian frame at first order gives Equation (15). Note that the divergence of the perturbed velocity in the Lagrangian frame is not zero at any order; the non-zero terms arise from the transformation between the Eulerian and Lagrangian coordinate. In order to satisfy Equation (15), we assume p⁽¹⁾ to have the form

$\begin{array}{ll}{\text{p}}^{(1)}(X,Y,t)=\left\{\frac{\partial \psi }{\partial Y}-ct\frac{\partial \psi }{\partial X},-\frac{\partial \psi }{\partial X},0\right\},& (17)\end{array}$

where, ψ is a scalar function of Lagrangian coordinates and time i.e., ψ ≡ ψ(X, Y, t). This gives:

$\begin{array}{ll}{\dot{\text{p}}}^{(1)}(X,Y,t)=\left\{\frac{\partial \dot{\psi }}{\partial Y}-c\frac{\partial \psi }{\partial X}-ct\frac{\partial \dot{\psi }}{\partial X},-\frac{\partial \dot{\psi }}{\partial X},0\right\}.& (18)\end{array}$

ψ is analogous to a stream-function, but not the same as one encountered in the usual Orr-Sommerfeld analysis. Substituting the form of Equation (17) in Equation (16) and simplifying gives:

$\begin{array}{l}\frac{d}{dt}\left[(1+c^2{t}^2)\frac{{\partial }^2}{\partial X^2}+\frac{{\partial }^2}{\partial Y^2}-2ct\frac{{\partial }^2}{\partial X\partial Y}\right]\dot{\psi }\\ =\nu {\left[(1+c^2{t}^2)\frac{{\partial }^2}{\partial X^2}+\frac{{\partial }^2}{\partial Y^2}-2ct\frac{{\partial }^2}{\partial X\partial Y}\right]}^2\dot{\psi }.& (19)\end{array}$

This system can be recast as:

$\begin{array}{ll}{A}\dot{\psi }=\varphi ,& (20)\end{array}$

$\begin{array}{ll}\dot{\varphi }=\nu {A}\varphi ,& (21)\end{array}$

where the operator ${A}$

$\begin{array}{ll}{A}=(1+c^2{t}^2)\frac{{\partial }^2}{\partial X^2}+\frac{{\partial }^2}{\partial Y^2}-2ct\frac{{\partial }^2}{\partial X\partial Y}.& (22)\end{array}$

The solution is obtained by first solving Equation (21) for ϕ and then solving Equation (20) for $\dot{\psi }$. Integrate to get ψ.

3.3. Initial and Boundary Conditions

The symmetry of the underlying flow implies periodic boundary conditions along the X-axis for solving both ϕ and $\dot{\psi }$. 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 defined by Equations (20) and (21) is second order in time and fourth order in the spatial variable Y. Accordingly two temporal initial conditions and four boundary conditions (two for $\dot{\psi }$ and two for ϕ) 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₀(X, Y) formally denote this initial perturbation.

$\begin{array}{ll}{\dot{\text{p}}}^{(1)}(X,Y,t=0)={\text{v}}_0(X,Y).& (23)\end{array}$

The definition of the Lagrangian coordinate provides the other initial condition. Equation (13) for n = 1 is:

$\begin{array}{ll}{\text{p}}^{(1)}(X,Y,t=0)=0\forall X,Y.& (24)\end{array}$

The boundary conditions are more involved. The symmetry of the underlying flow implies periodic boundary conditions along the X-axis for solving both ϕ and $\dot{\psi }$. The no slip condition imposed on the wall at Y = 0 means:

$\begin{array}{ll}{\dot{\text{p}}}^{(1)}(X,Y=0,t)=0\forall t.& (25)\end{array}$

These are two conditions corresponding to the X and Y coordinate. Note that $\dot{\text{p}}$ depends explicitly on ψ and $\dot{\psi }$ 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 $\dot{\psi }$. Using the definition of ${\dot{\text{p}}}^{(1)}$ from Equation (18) in Equation (25) gives:

$\begin{array}{ll}\frac{\partial \dot{\psi }}{\partial Y}{|}_{Y=0}=c\frac{\partial \psi }{\partial X}{|}_{Y=0}\forall t,& (26)\end{array}$

$\begin{array}{ll}\frac{\partial \dot{\psi }}{\partial X}{|}_{Y=0}=0\forall t.& (27)\end{array}$

Using Equation (27) and using the fact that the time derivative commutes with the spatial derivative gives $\frac{d}{dt}(\frac{\partial \psi }{\partial X}{|}_{Y=0})=0$, at all times t. Evaluating the Y-coordinate of Equation (24) at Y = 0 using Equation (17) gives $\frac{\partial \psi }{\partial X}{|}_{Y=0}=0$ at the initial time t = 0. Since the time derivative is zero, $\frac{\partial \psi }{\partial X}{|}_{Y=0}$ stays zero at all times. Thus the r.h.s. of Equation (26) is zero at all times and it follows that:

$\begin{array}{ll}\frac{\partial \dot{\psi }}{\partial Y}{|}_{Y=0}=0\forall t.& (28)\end{array}$

We assume periodic b.c. (Fourier decomposition) along the Y axis for ϕ. Equations (27) and (28) provide the boundary conditions for $\dot{\psi }$.

3.4. 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:

$\begin{array}{ll}\varphi (X,Y,t)={\displaystyle \sum _{k_x,k_y}}\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}{e}^{i{k}_{y}Y}f(t).& (29)\end{array}$

Substituting in Equation (21), gives:

$\begin{array}{ll}\frac{df}{dt}=\nu \left(-k_x^2(1+c^2{t}^2)-k_y^2+2ct{k}_x{k}_y\right)f.& (30)\end{array}$

The solution is:

$\begin{array}{ll}f(t)=f(0)e^{\nu \left(-k^{2}t+k_x{k}_{y}c{t}^2-\frac{k_x^2{c}^2{t}^3}{3}\right)}& (31)\end{array}$

where $k^2=k_x^2+k_y^2$ and f(0) is the integration constant to be set later. The solution for ϕ is:

$\begin{array}{ll}\varphi (X,Y,t)=f(0){\displaystyle \sum _{k_x,k_y}}\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}{e}^{i{k}_{y}Y}{e}^{\nu \left(-k^{2}t+k_x{k}_{y}c{t}^2-\frac{k_x^2{c}^2{t}^3}{3}\right)}.& (32)\end{array}$

One can now solve Equation (20) for $\dot{\psi }$. From the structure of the equation one can assume that the purely temporal functions are the same for both ϕ and $\dot{\psi }$. Any additional time dependence introduced via the operator ${A}$ is necessarily also a function of Y. This gives the ansatz:

$\begin{array}{ll}\dot{\psi }(X,Y,t)={\displaystyle \sum _{k_x,k_y}}\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}{g}_{k_y}(Y,t)f(t).& (33)\end{array}$

The subscript k_y denotes the dependence of g on the parameter k_y. We will drop it in subsequent evaluations. Substituting in Equation (20) and using Equation (32) gives:

$\begin{array}{ll}-k_x^2(1+c^2{t}^2)g(Y,t)-2i{k}_{x}ct{g}^{\prime }(Y,t)+g^{″}(Y,t)=e^{i{k}_{y}Y},& (34)\end{array}$

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: g(Y, t) = g_homo.(Y, t) + g_part.(Y, t) where

$\begin{array}{ll}{g}_{homo.}(Y,t)=C_1(t)e^{(-k_x+k_{x}ict)Y}+C_2(t)e^{(k_x+k_{x}ict)Y}& (35)\end{array}$

and

$\begin{array}{ll}{g}_{part.}(Y,t)=C_3(t)e^{i{k}_{y}Y},& (36)\end{array}$

where

$\begin{array}{ll}{C}_3^{-1}(t)=-k_x^2(1+c^2{t}^2)-k_y^2+2k_x{k}_{y}ct.& (37)\end{array}$

Satisfying the boundary conditions represented by Equations (27) and (28) fixes C₁(t) and C₂(t).

$\begin{array}{l}{C}_1^{-1}(t)=2k_x(kx(1-ict)+i{k}_y),\\ C_2^{-1}(t)=2k_x(kx(1+ict)-i{k}_y).& (38)\end{array}$

The net solution for g is

$\begin{array}{ll}g(Y,t)=C_1(t)e^{(-k_x+ict{k}_x)Y}+C_2(t)e^{(k_x+ict{k}_x)Y}+C_3(t)e^{i{k}_{y}Y}.& (39)\end{array}$

ψ is obtained by integrating $\dot{\psi }$ w.r.t. time:

$\begin{array}{ll}\psi ={\displaystyle \sum _{k_x,k_y}}\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}{\displaystyle \int }g(Y,t^{\prime })f(t^{\prime })d{t}^{\prime }+C(X,Y),& (40)\end{array}$

where C(X, Y) is the constant of integration. This can be re-written as:

$\begin{array}{ll}\psi ={\displaystyle \sum _{k_x,k_y}}f(0)\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}h(Y,t)+C(X,Y),& (41)\end{array}$

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 $C(X,Y)=-f(0)\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}h(Y,t=0)$ giving:

$\begin{array}{ll}\psi ={\displaystyle \sum _{k_x,k_y}}f(0)\tilde{\varphi }(k_x,k_y)e^{i{k}_{x}X}\left[h(Y,t)-h(Y,t=0)\right].& (42)\end{array}$

The only quantities remaining to be determined are f(0) and $\tilde{\varphi }(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 behavior: 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>>{(c{k}_x)}^{-1}$ they integrate to zero leaving:

$\begin{array}{ll}\psi ~e^{i{k}_{x}X}{e}^{i{k}_{y}Y}{\displaystyle \int }{C}_3(t^{\prime })f(t^{\prime })d{t}^{\prime }& (43)\end{array}$

Substituting for C₃ and f from Equations (31) and (37) gives:

$\begin{array}{ll}\begin{array}{lll}\psi \hfill & ~\hfill & \tilde{\psi }(k_x,k_y)f(0)e^{i{k}_{x}X}{e}^{i{k}_{y}Y}{{\displaystyle \int }}^{}-\frac{e^{-\frac{\nu k_x^2{c}^{2}t{\prime }^3}{3}}}{k_x^2{c}^{2}t{\prime }^2}dt\prime \hfill \end{array}& (44)\end{array}$

$\begin{array}{lll}~& \frac{f(0)}{k_x^2{c}^{2}t}\left[e^{-\frac{\nu k_x^2{c}^2{t}^3}{3}}+t{\left(\frac{\nu k_x^2{c}^2}{3}\right)}^{1/3}\gamma \left(\frac{2}{3},\frac{\nu k_x^2{c}^2{t}^3}{3}\right)\right]+C(X,Y)& (45)\end{array}$

where $\gamma (s,z)={\int }_o^z{e}^{-z\prime }z{\prime }^{s-1}dz\prime$ is the lower incomplete gamma function.

Thus ψ and hence the velocity perturbation evolves as $~e^{-\frac{\nu k_x^2{c}^2{t}^3}{3}}$ at late times and the flow is linearly stable in the Lagrangian frame. This late time behavior is in exact agreement with that of Press and Marcus [10] 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. One has to choose appropriate basis functions which will satisfy them. Although we do not provide a complete numerical solution for this case, we briefly sketch its form in Appendix 6.4. We note that the temporal function f(t) has a similar exponential form; velocity perturbation evolves as $~e^{-\frac{\nu k_x^2{c}^2{t}^3}{3}}$ at late times. Defining the Reynolds number for the bounded flow of height h as Re = v^s.s.(h)·h/ν = ch²/ν, the perturbation evolves as $~exp(-\frac{h^2{k}_x^2{c}^3{t}^3}{3Re})$.

3.5. Discussion

The main temporal dependence of the Lagrangian solution is given by the exponential term in Equation (31). All three terms in the exponent arise from the viscous component as is indicated by the factor ν multiplying them. The linear term −νk²t in the exponent also arises in standard Eulerian linear perturbation theory, but quadratic and cubic terms are new in the Lagrangian frame. It is also interesting to note that the solution for a small time interval can grow as $e^{\nu k_x{k}_{y}c{t}^2}$ before eventually falling off as $e^{-\nu k_x^2{c}^2{t}^3/3}$ at late times. For sufficiently large k_y, this term can dominate the dynamics. This hints at the phenomenon of transient growth which is believed to be responsible for the instabilities observed in shear flow experiments. In Eulerian perturbation theory this transient growth has been attributed to the fact that the linear stability operator is non-normal. The stability in these cases is not governed just by the spectrum, but by its pseudospectrum [19, 20]. However, the linear transient growth explanation has also found critics [43]. In particular, the effect of the non-linearity on the mean background flow plays an important role in the transition to turbulence. Recently, Fukumoto et al. [44] used Lagrangian approach to study the weakly non-linear stability in the case of elliptical flows. Understanding the transient growth and weakly non-linear regime for the plane Couette flow using LPT is left for future work.

Stability analyses in the Eulerian and Lagrangian frame differ in one fundamental aspect. In the former, non-linear convective terms such as (v·∇)v, which are second order in the velocity perturbation are ignored whereas in the Lagrangian frame, they get absorbed into the time derivative operator. This suggests a potential advantage of the Lagrangian frame over the Eulerian frame. However, other than the fact that we get exponents that are quadratic and cubic in t, we do not get qualitatively different results at linear order in the two frames. In Appendix 6.5 we estimate the effect of the convective term in Eulerian perturbation theory using an interative procedure. It shows qualitatively a similar behavior as compared to the first order Lagrangian solution, but there is no quantitative agreement nor does it provide any additional quantitative insights into the nature of transient growth. Finally, we must remind ourselves that there is never a clear matching between orders in the Eulerian and Lagrangian frame. A more quantitative comparison of the linear LPT solution to that obtained in the Eulerian picture can be done only when the Lagrangian to Eulerian map is inverted as is described in the next section (see section 4). This requires further numerical work and is beyond the scope of this paper.

4. Transforming Back to the Eulerian Frame

The LPT scheme outlined in the previous section solves 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). In recent work, Nadkarni-Ghosh and Chernoff [45] 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 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), but outline its solution. Let r_phys denote the solution in the physical frame that satisfies the original set and r_LPT denote the solution in the calculational frame obtained by the perturbative treatment discussed in the previous section. The two are related as:

$\begin{array}{ll}\text{r}={\text{r}}_{LPT}+\Delta \text{r}(t).& (46)\end{array}$

By substituting Equation (47) 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 $\Delta \dot{\text{r}}(0)=0$. The net physical solution and velocity:

$\begin{array}{ll}\text{r}(\text{R},t)=\text{R}+{\text{p}}^{(0)}(\text{R},t)+{\text{p}}^{(1)}(\text{R},t)+\Delta \text{r}(t),& (47)\end{array}$

$\begin{array}{ll}{\text{v}}_L(\text{R},t)={\dot{\text{p}}}^{(0)}(\text{R},t)+{\dot{\text{p}}}^{(1)}(\text{R},t)+\Delta \dot{\text{r}}(t).& (48)\end{array}$

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 $\text{r}={F}(\text{R},t)$ then the Eulerian velocity at the coordinate point r is:

$\begin{array}{ll}{\text{v}}_E(\text{r},t)={\text{v}}_L({{F}}^{-1}(\text{r},t),t).& (49)\end{array}$

5. Conclusion

The main motivation behind this paper was to outline the formalism of LPT, a technique successful in other branches of physics, and apply it to the problem of flow stability. LPT has been used to model non-linear growth in cosmology for over four decades. Early work was done in the 1970's by Zeldovich [46] with linear order perturbation theory. In the 1990s, others including Buchert and collaborators [47] developed higher orders and recently Nadkarni-Ghosh and Chernoff [45, 48] addressed issues of convergence of the theory prior to orbit crossing. To the best of our knowledge, such techniques have been sparingly used to investigate the stability of fluid flows. As was outlined in the introduction, the fields of turbulence and large scale growth of structure in the universe share a common kind of complexity. In both, linear Eulerian stability theory is usually employed to get analytical results and the non-linear regime is often modeled using numerical simulations. It has often been the case that techniques developed in one branch of science later found utility in other seemingly unrelated branches of science. New techniques bring new insights and new ideas come from mapping concepts of one set of things to another. Higher order perturbation theory is always complicated, whether in Eulerian or Lagrangian frames, but is nevertheless useful to give analytical insights and can sometimes be used in conjunction with simulations to improve their efficiency.

Here, we presented a first step in this direction. We focussed on the simplest shear flow: the plane Couette flow and restricted to 2D perturbations. This greatly simplified the expressions. In the Eulerian analysis, it is often enough to consider 2D perturbations thanks to Squire's theorem [49], which states that an unstable 3D eigenmode for some Reynolds number implies a unstable 2D eigenmode for a lower Reynolds number. Unfortunately, Squire's theorem, which is based on the Orr-Sommerfeld analysis, need not apply to the Lagrangian analysis so one may not be able to make conclusions based on 2D perturbations. Furthermore, it is known that transient growth is usually weaker in 2D than in 3D perturbations [23]. For 3D perturbations or other types of flows the framework remains the same; of course the form of the equations is more complicated and their analysis will require numerics. This complexity is to be expected, but is not formally prohibitive.

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 [45, 48]). 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. Recent work [30, 50] also addresses the issue of analyticity of Lagrangian particle trajectories from a more formal perspective. The level of complexity of higher order LPT is perhaps comparable to higher order perturbation theory in the Eulerian frame. Each mode may have an associated time scale of evolution and their interaction may possibly complicate convergence. This phenomenon has been demonstrated on a two time-scale problem, where the scales are well separated: one fast and one slow (see Berry and Shukla [51] and references therein). It may be possible to apply LPT recursively to get 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. Nevertheless, Lagrangian perturbative methods provide an alternate way to analyze the stability of fluid flows and the solutions could either be used in isolation or could be potentially useful in making educated guesses for starting numerical simulations that aim to understand the transition to turbulence.

Author Contributions

SN-G and JB formulated the problem. SN-G performed the analytic calculations and wrote the paper. JB supervised the project.

Funding

SN-G would like to thank the Science and Engineering Research Board (SERB) of the Department of Science and Technology (DST), Govt. of India for the grant YSS/2014/000526.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

SN-G would like to thank Mahendra Verma, Stefano Anselmi, and Rama Govindrajan for valuable discussions and suggesting useful references. Thanks are due also to Sagar Chakraborty for discussions and comments on the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2018.00037/full#supplementary-material

Footnotes

1. ^We have chosen dimensional units since for the semi-bounded flow, there is no inherent length scale to define the Reynolds number. Thus, p is a displacement and ψ has dimensions of m²; g~m²; f~t⁻¹.

References

1. Orr WMF. The stability or instability of the steady motions of a liquid. Part I. Proc Irish Acad. (1907) 27:9–68.

2. Sommerfeld A. Ein Beitrag zur hydrodynamische Erklärung der turbulenten Flüssigkeitsbewegungen. Proceedings of 4th International Congress on Mathematics Rome (1908). p. 116–24.

3. Bayly BJ, Orszag SA, Herbert T. Instability mechanisms in shear-flow transition. Annu Rev Fluid Mech. (1988) 20:359–91. doi: 10.1146/annurev.fl.20.010188.002043

CrossRef Full Text | Google Scholar

4. Joseph DD. Eigenvalue bounds for the Orr-Sommerfeld equation. J Fluid Mech. (1968) 33:617–21. doi: 10.1017/S0022112068001552

CrossRef Full Text | Google Scholar

5. Case KM. Stability of inviscid plane couette flow. Phys Fluids (1960) 3:432–5. doi: 10.1063/1.1706055

CrossRef Full Text | Google Scholar

6. Lin CC. Some mathematical problems in the theory of the stability of parallel flows. J Fluid Mech. (1961) 10:430–8. doi: 10.1017/S0022112061001025

CrossRef Full Text | Google Scholar

7. Orszag SA. Accurate solution of the Orr-Sommerfeld stability equation. J Fluid Mech. (1971) 50:689–703. doi: 10.1017/S0022112071002842

CrossRef Full Text | Google Scholar

8. Hughes TH. Variable mesh numerical method for solving the orr-sommerfeld equation. Phys Fluids (1972) 15:725–8. doi: 10.1063/1.1693974

CrossRef Full Text | Google Scholar

9. Davey A. On the stability of plane Couette flow to infinitesimal disturbances. J Fluid Mech. (1973) 57:369–80. doi: 10.1017/S0022112073001217

CrossRef Full Text | Google Scholar

10. Marcus PS, Press WH. On Green's functions for small disturbances of plane Couette flow. J Fluid Mech. (1977) 79:525–34.

Google Scholar

11. Tillmark N, Alfredsson PH. Experiments on transition in plane Couette flow. J Fluid Mech. (1992) 235:89–102. doi: 10.1017/S0022112092001046

CrossRef Full Text | Google Scholar

12. Daviaud F, Hegseth J, Bergé P. Subcritical transition to turbulence in plane Couette flow. Phys Rev Lett. (1992) 69:2511–4.

PubMed Abstract | Google Scholar

13. Lemoult G, Aider J-L, Wesfreid JE. Experimental scaling law for the subcritical transition to turbulence in plane Poiseuille flow. Phys Rev E (2012) 85:025303. doi: 10.1103/PhysRevE.85.025303

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Orszag SA, Patera AT. Subcritical transition to turbulence in plane channel flows. Phys Rev Lett. (1980) 45:989–93. doi: 10.1103/PhysRevLett.45.989

CrossRef Full Text | Google Scholar

15. Orszag SA, Kells LC. Transition to turbulence in plane Poiseuille and plane Couette flow. J Fluid Mech. (1980) 96:159–205. doi: 10.1017/S0022112080002066

CrossRef Full Text | Google Scholar

16. Cherhabili A, Ehrenstein U. Finite-amplitude equilibrium states in plane Couette flow. J Fluid Mech. (1997) 342:159–77. doi: 10.1017/S0022112097005661

CrossRef Full Text | Google Scholar

17. Barkley D, Tuckerman LS. Stability analysis of perturbed plane Couette flow. Phys Fluids (1999) 11:1187–95. doi: 10.1063/1.869987

CrossRef Full Text | Google Scholar

18. Dauchot O, Daviaud F. Streamwise vortices in plane Couette flow. Phys Fluids (1995) 7:901–3. doi: 10.1063/1.868565

CrossRef Full Text | Google Scholar

19. Schmid PJ, Henningson DS, Khorrami MR, Malik MR. A study of eigenvalue sensitivity for hydrodynamic stability operators. Theor Comput Fluid Dyn. (1993) 4:227. doi: 10.1007/BF00417929

CrossRef Full Text | Google Scholar

20. Trefethen LN. Pseudospectra of linear operators. SIAM Rev. (1997) 39:383–406. doi: 10.1137/S0036144595295284

CrossRef Full Text | Google Scholar

21. Butler KM, Farrell BF. Three-dimensional optimal perturbations in viscous shear flow. Phys Fluids A (1992) 4:1637–50. doi: 10.1063/1.858386

CrossRef Full Text | Google Scholar

22. Reddy SC, Henningson DS. Energy growth in viscous channel flows. J Fluid Mech. (1993) 252:209–38. doi: 10.1017/S0022112093003738

CrossRef Full Text | Google Scholar

23. Trefethen LN, Trefethen AE, Reddy SC, Driscoll TA. Hydrodynamic stability without eigenvalues. Science (1993) 261:578–84. doi: 10.1126/science.261.5121.578

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Grossmann S. The onset of shear flow turbulence. Rev Mod Phys. (2000) 72:603–18. doi: 10.1103/RevModPhys.72.603

CrossRef Full Text | Google Scholar

25. Schmid PJ. Nonmodal stability theory. Annu Rev Fluid Mech. (2007) 39:129–62. doi: 10.1146/annurev.fluid.38.050304.092139

CrossRef Full Text | Google Scholar

26. La-Porta A, Voth GA, Crawford AM, Alexander J, Bodenschatz E. Fluid particle accelerations in fully developed turbulence. Nature (2001) 409:1017–9. doi: 10.1038/35059027

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Mordant N, Metz P, Michel O, Pinton JF. Measurement of lagrangian velocity in fully developed turbulence. Phys Rev Lett. (2001) 87:214501. doi: 10.1103/PhysRevLett.87.214501

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Falkovich G, Haitao X, Alain P, Bodenschatz E, Biferale L, Boffetta G, et al. On Lagrangian single-particle statistics. Phys Fluids (2012) 24:055102-8. doi: 10.1063/1.4711397

CrossRef Full Text | Google Scholar

29. Andrews DG, McIntyre ME. An exact theory of nonlinear waves on a Lagrangian-mean flow. J Fluid Mech (1978) 89:609–46. doi: 10.1017/S0022112078002773

CrossRef Full Text | Google Scholar

30. Zheligovsky V, Frisch U. Time-analyticity of Lagrangian particle trajectories in ideal fluid flow. J Fluid Mech. (2014) 749:407. doi: 10.1017/jfm.2014.221

CrossRef Full Text | Google Scholar

31. Taruya A, Hiramatsu T. A closure theory for nonlinear evolution of cosmological power spectra. Ap J. (2008) 674:617–35. doi: 10.1086/526515

CrossRef Full Text | Google Scholar

32. Kraichnan RH. The structure of isotropic turbulence at very high Reynolds numbers J Fluid Mech. (1959) 5:497–543. doi: 10.1017/S0022112059000362

CrossRef Full Text | Google Scholar

33. Kida S, Goto S. A Lagrangian direct-interaction approximation for homogeneous isotropic turbulence J Fluid Mech. (1997) 345:307–45. doi: 10.1017/S0022112097006289

CrossRef Full Text | Google Scholar

34. Goto S, Kida S. Direct-interaction approximation and Reynolds-number reversed expansion for a dynamical system. Phys D (1998) 117:191–214. doi: 10.1016/S0167-2789(97)00314-X

CrossRef Full Text | Google Scholar

35. Crocce M, Scoccimarro R. Renormalized cosmological perturbation theory. Phys Rev D (2006) 73:63519. doi: 10.1103/PhysRevD.73.063519

CrossRef Full Text | Google Scholar

36. Wyld HW. Formulation of the theory of turbulence in an incompressible fluid. Ann Phys. (1961) 14:143–65. doi: 10.1016/0003-4916(61)90056-2

CrossRef Full Text | Google Scholar

37. L'vov V, Procaccia I. Exact resummations in the theory of hydrodynamic turbulence. I. The ball of locality and normal scaling. Phys Rev E (1995) 52:3840–57.

38. Gurbatov SN, Saichev IA, Shandarin SF. The large-scale structure of the universe in the frame of the model equation of non-linear diffusion Mon Not R Astron Soc. (1989) 236:385–402. doi: 10.1093/mnras/236.2.385

CrossRef Full Text | Google Scholar

39. Frisch U, Bec J. Burgulence. In: Lesieur M, Yaglom A, David F, editors New Trends in Turbulence (Berlin; Heidelberg: EDP Sciences, Springer-Verlag) (2001). p. 341–83.

Google Scholar

40. Frisch U, Bec J, Villone B. Singularities and the distribution of density in the Burgers/adhesion model Phys Rev D (2001). 73:63519. doi: 10.1016/S0167-2789(01)00195-6

CrossRef Full Text | Google Scholar

41. Gaite J. A non-perturbative Kolmogorov turbulence approach to the cosmic web structure. Europhys Lett. (2012) 98:49002. doi: 10.1209/0295-5075/98/49002

CrossRef Full Text | Google Scholar

42. Pierson WJ. Perturbation analysis of the Navier-Stokes equations in lagrangian form with selected linear solutions. J Geophys Res. (1962) 67:3151. doi: 10.1029/JZ067i008p03151

CrossRef Full Text | Google Scholar

43. Waleffe F. Transition in shear flows. Nonlinear normality versus non-normal linearity. Phys Fluids (1995) 7:3060–6. doi: 10.1063/1.868682

CrossRef Full Text | Google Scholar

44. Fukumoto Y, Hirota M, Mie Y. Lagrangian approach to weakly nonlinear stability of elliptical flow. Phys Scripta (2010) 142:014049. doi: 10.1088/0031-8949/2010/T142/014049

CrossRef Full Text | Google Scholar

45. Nadkarni-Ghosh S, Chernoff DF. Modelling non-linear evolution using Lagrangian perturbation theory re-expansions. Mon Not R Astron Soc. (2013) 431:799–823. doi: 10.1093/mnras/stt217

CrossRef Full Text

46. Zel'dovich YB. Gravitational instability: an approximate theory for large density perturbations. Astron Astrophys. (1970) 5:84–9.

47. Ehlers J, Buchert T. Newtonian cosmology in Lagrangian formulation: foundations and perturbation theory. Gen Relativ Gravit. (1997) 29:733–64. doi: 10.1023/A:1018885922682

CrossRef Full Text | Google Scholar

48. Nadkarni-Ghosh S, Chernoff DF. Extending the domain of validity of the Lagrangian approximation. Mon Not R Astron Soc. (2011) 410:1454–88. doi: 10.1111/j.1365-2966.2010.17529.x

CrossRef Full Text | Google Scholar

49. Squire HB. On the stability for three-dimensional disturbances of viscous fluid flow between parallel walls. Proc R Soc Lond Series A (1933) 142:621–8. doi: 10.1098/rspa.1933.0193

CrossRef Full Text | Google Scholar

50. Rampf C, Villone B, Frisch U. How smooth are particle trajectories in a ∧ CDM Universe? Mon Not R Astron Soc. (2015) 452:1421–36. doi: 10.1093/mnras/stv1365

CrossRef Full Text | Google Scholar

51. Berry MV, Shukla P. Typical weak and superweak values. J Phys A Math Theor. (2010) 43:045102. doi: 10.1088/1751-8113/43/4/045102

CrossRef Full Text | Google Scholar