The linearized pressure Poisson equation for global instability analysis of incompressible flows

The linearized pressure Poisson equation (LPPE) is used in two and three spatial dimensions in the respective matrix-forming solution of the BiGlobal and TriGlobal eigenvalue problem in primitive variables on collocated grids. It provides a disturbance pressure boundary condition which is compatible with the recovery of perturbation velocity components that satisfy exactly the linearized continuity equation. The LPPE is employed to analyze instability in wall-bounded flows and in the prototype open Blasius boundary layer flow. In the closed flows, excellent agreement is shown between results of the LPPE and those of global linear instability analyses based on the time-stepping nektar++, Semtex and nek5000 codes, as well as with those obtained from the FreeFEM++ matrix-forming code. In the flat plate boundary layer, solutions extracted from the two-dimensional LPPE eigenvector at constant streamwise locations are found to be in very good agreement with profiles delivered by the NOLOT/PSE space marching code. Benchmark eigenvalue data are provided in all flows analyzed. The performance of the LPPE is seen to be superior to that of the commonly used pressure compatibility (PC) boundary condition: at any given resolution, the discrete part of the LPPE eigenspectrum contains converged and not converged, but physically correct, eigenvalues. By contrast, the PC boundary closure delivers some of the LPPE eigenvalues and, in addition, physically wrong eigenmodes. It is concluded that the LPPE should be used in place of the PC pressure boundary closure, when BiGlobal or TriGlobal eigenvalue problems are solved in primitive variables by the matrix-forming approach on collocated grids.

with a new configuration or parameter range, since it ensures access to a large part of the eigenspectrum, computation of which scales linearly with the size of the Krylov subspace dimension chosen. This advantage may outweigh the shared-memory requirements and parallelization issues arising from the discretization and storage of large matrices. The present contribution is devoted to one specific, but important, aspect of the matrix-forming approach, namely the boundary condition to be imposed on the pressure perturbations in the incompressible limit, if the primitive-variable formulation is used and both perturbation velocity components and pressure perturbation are collocated on the same grid.
The issue of boundary conditions for a pressure perturbation that ensures flow incompressibility is typically absent in time-stepping methods, since it has been dealt with in the direct numerical simulation code used, e.g., by time-advancing velocity fields that satisfy the continuity equation, as done in the nektar++ [10,30], nek5000 [13,16] or Semtex [7] spectral element codes. The analogous procedure of weak formulation of the equations of motion, followed in the finite-element package FreeFEM++ [25], ensures that no issues exist regarding appropriate boundary conditions for the pressure perturbation, when this package is adapted to perform global linear instability analysis [12,17,32,53].
In a matrix-forming context, the problem of determining correct pressure boundary conditions can be circumvented when the flow to be analyzed comprises only a base flow velocity component along the homogeneous spatial direction. In that case cross-differentiation of the governing equations permits elimination of the pressure perturbation, at the cost of increasing the order of derivatives to be discretized numerically. The two-dimensional analog of the Orr-Sommerfeld equation arises, first derived and solved by Tatsumi and Yoshimura [54]. In this relatively simple geometry spectral methods are optimally employed to discretize the governing equations using basis recombination techniques which exactly satisfy the Dirichlet and Neumann type of perturbation velocity boundary conditions [8,54,67] at the expense of a relatively large number of nodes needed for the accurate representation of the resulting third-and fourth-order derivatives. When base velocity components additional to that along the homogeneous spatial direction are present the advantage of eliminating pressure perturbation is lost and the linearized Navier-Stokes equations in primitive variables must be closed with a priori unknown boundary conditions for the pressure perturbation. In this situation, another means of circumventing the need for pressure boundary conditions has been discussed in the literature on incompressible direct numerical simulation using high-order methods [28,29], where a staggered grid for pressure is used. Using spectral collocation on staggered grids, the momentum equations are collocated and solved on two-dimensional tensor-product grids based on extremum (e.g., Gauss-Lobatto) gridpoints, continuity is solved on interior (e.g., Gauss) gridpoints and spectral interpolation operators are used in order to transfer velocity from the extremum onto the interior grid and vice versa [40]. This approach has been used successfully in the context of global instability analysis by Theofilis and Colonius [62] to obtain eigenvalue problem results in a compressible flat plate boundary layer as limiting validation cases of the algorithm used for the solution of the BiGlobal eigenvalue problem in compressible open cavity flow 1 When a collocated scheme for the disturbance velocity and pressure eigenfunctions is used, boundary conditions for the latter quantity must be provided at the (solid or open) domain boundary. Some guidance regarding the appropriate choice of these boundary conditions may be sought in the persisting discussion in the incompressible direct numerical simulation community [21,47] where solutions are being proposed for the analogous problem of pressure boundary conditions that ensure satisfaction of continuity; a full discussion can be found in the excellent review article of Rempfer [47]. Two of the key ideas in the latter work, namely that: ...it is illegal to write down the momentum equation taken at the boundary and derive a pressure boundary condition from it by simply projecting the result on the wall-normal coordinate... and, further, ...to close a differential problem, the boundary conditions need to provide some additional information that is not already contained in the field equations... are relevant to the discussion that follows. It should, however, be stressed that in the present matrix-forming context the issues arising from the splitting formulation of the incompressible equations of motion, where inappropriate use of pressure boundary conditions in the time-integration algorithm may violate mass conservation, are not directly related to the boundary conditions to close the linearized Navier-Stokes system. Pressure compatibility (PC) boundary conditions, derived from the momentum equations and collocated at the wall, have been shown to perform well in classic linear stability analysis [41] and their extension in two spatial dimensions is the natural candidate to provide the sought closure of the system of linearized equations of motion; indeed, these conditions have been used successfully in the global instability analysis of a number of well-studied incompressible flows [35,36,57,63,64].
On the other hand, various forms of the pressure Poisson equation have been used in order to provide boundary conditions consistent with mass conservation in algorithms for the solution of the incompressible Navier-Stokes equations. This equation could be considered as providing the additional information in the sense of [47]. This paper discusses the derivation and implementation of an appropriate linearized pressure Poisson equation (LPPE) to provide boundary conditions for the pressure perturbation in closed and open flows with two or three inhomogeneous spatial directions. Section 2 presents the derivation of the LPPE in the context of incompressible BiGlobal and TriGlobal analysis. In Sect. 3, instability analysis results of the LPPE boundary closure are compared in terms of eigenspectra and amplitude functions with solutions of the twoand three-dimensional eigenvalue problem obtained herein using the nektar++, Semtex and nek5000 timestepping codes, as well as the FreeFEM++ matrix-forming code, as applied to wall-bounded flow examples in which only discrete eigenspectra exist. In addition, results of the two-dimensional eigenvalue problem subject to the LPPE and PC boundary closures are obtained in the incompressible flat plate boundary layer, which comprises both a discrete and a continuous eigenspectrum; here results extracted from the two-dimensional eigenfunctions are compared with profiles predicted by the NOLOT/PSE code. All results obtained establish the LPPE boundary closure as a reliable means of closing the BiGlobal or TriGlobal EVP and demonstrate its superior performance over the PC boundary closure. A short discussion summarizes the findings in Sect. 4.

Modal linear BiGlobal analysis
Linear modal global instability analysis of incompressible flows developing in domains in which one spatial direction, z in Cartesian coordinates, can be considered homogeneous, can be performed by substituting the decomposition q(x, y, z, t) =q(x, y) + εq(x, y)e i(βz−ωt) , ε 1 ( 1 ) into the incompressible Navier-Stokes and continuity equations. In the most general case, the base flow comprises all three velocity components and pressure,q = (ū,v,w,p) T , and the linearized system of equations governing modal BiGlobal instability reads [60]û Solution of the linearized Navier-Stokes equations (LNSE) defined by system (2)(3)(4)(5) for the determination of the two-dimensional velocity and pressure amplitude functions,q = (û,v,ŵ,p) T , may be sought in either a temporal (β ∈ R, ω ∈ C) or spatial (β ∈ C, ω ∈ R) framework. In either case, conditions for both the perturbation velocity field and the perturbation pressurep must be provided at the domain boundary. Boundary conditions for the velocity depend on the type of domain boundary considered. At solid walls, no-slip is straightforwardly imposed through homogeneous Dirichlet conditions on all components of the perturbation velocity vector. The same type of conditions are also imposed at the inflow boundary of open domains, when the analysis aims at excluding perturbations from entering the computational domain [65]. Vanishing of linear perturbations is also imposed at far-field boundaries, if the latter are taken far away from solid surfaces, although it should be stressed that, much like the situation in classic linear stability theory, this choice misrepresents or altogether eliminates the continuous branch of perturbations oscillatory to infinity [49,56]. At open outflow boundaries, it is in principle unclear what form the perturbation velocity may assume, although linear extrapolation from the interior has been found to not only work well, but also to have little effect on the form of linear perturbations in the interior of the domain, when the latter is analytically known [64]. Assuming that the three velocity components provide boundary conditions for three equations in the system (2)(3)(4)(5), a boundary condition for the pressure perturbation, consistent with the role of this Lagrangian multiplier in enforcing incompressibility of the disturbance velocity field, needs to be supplied in order to close the system.
Substitution of (6) into the incompressible Navier-Stokes and continuity equations leads to the TriGlobal eigenvalue problemû x +v y +ŵ z = 0, The real TriGlobal eigenvalue problem defined by (7-10) can be solved for the determination of the eigenvalues ω ∈ C and the three-dimensional amplitude functionsq = (û,v,ŵ,p)T. While boundary conditions for the perturbation velocity components can be determined in a relatively straightforward manner as discussed in the previous section, a boundary condition for the pressure perturbation is again unknown.

The linearized pressure Poisson equation (LPPE)
The information to be used for the pressure perturbation at the boundary, which is independent of the field equations [47], is provided by the pressure Poisson equation which is obtained from the incompressible Navier-Stokes equations after taking the divergence of the momentum equation and applying continuity to eliminate the time derivative and the viscous term. Substituting the linear decomposition (1) into (11) and assuming that the terms at O(1) are identically satisfied by the base flow, at O(ε) the linearized pressure Poisson equation in two inhomogeneous spatial directions is obtained, The same equation is obtained if (3) and (4) are differentiated w.r.t. x and y, respectively, the result is added to (5) multiplied by the factor iβ and use of the base flow continuity equation as well as (2) is made. In a TriGlobal analysis context, decomposition (6) is introduced into equation (11) leading to the threedimensional linearized pressure Poisson equation ∂ x x + ∂ yy + ∂ zz p + 2 ū xûx +ū yvx +ū zŵx +v xû y +v yvy +v zŵy +w xûz +w yvz +w zŵz = 0, after subtraction of the base flow continuity terms and neglecting O(ε 2 ) terms. Both of equations (12) and (13) will be referred to in what follows as the linearized pressure Poisson equation (LPPE) since there is no risk of confusion between the respective BiGlobal and TriGlobal contexts in which they are used.
Key features of the LPPE are, firstly, the exact nature of this equation, which contrasts with approximate conditions based on expansions of a variable near a boundary, where the normal vector needs to be approximated usually by a low-order scheme, secondly, the strong coupling between pressure perturbation, the disturbance velocity and the underlying base flow velocity components and, thirdly, the presence of only (better-conditioned) first-order base flow derivatives to be discretized numerically. Perhaps the most interesting characteristic of the LPPE in both forms (12) and (13) is that it is independent of the flow Reynolds number, since the diffusion terms of both the basic flow and the linear perturbations are eliminated by virtue of the respective continuity equations. A noteworthy aspect of (12) in the context of BiGlobal analysis is the fact that the perturbation velocity component along the homogeneous spatial direction does not appear explicitly and also that all three basic flow velocity components appear, making the LPPE applicable to the limiting cases of the wavenumber vector being parallel or normal to the plane on which the base flow develops. Another point worth making in the case of BiGlobal analysis is the fact that in the derivation of the LPPE no use has been made of the modal temporal Ansatz (1). This implies that this equation could in principle also be used to perform non-modal BiGlobal instability analysis, although to date no such work has appeared in the literature.
The LPPE-based boundary closures (12) and (13) are collocated at the (grid-conforming) boundary as part of the spatial discretization procedure employed to solve the BiGlobal (2)(3)(4)(5) or TriGlobal (7-10) eigenvalue problem, respectively. Unlike typical splitting schemes of direct numerical simulation, either explicit or semi-implicit in time, in substeps of which the nonlinear coupling between pressure and velocity components is routinely eliminated, the LPPE closure retains this coupling and ensures conservation of mass at the level of linear flow perturbations. Depending on the type of the boundary at which it is imposed, the LPPE may be simplified on account of tangential or normal base flow velocity components vanishing at the boundary. This is particularly important at corners of solid walls, where no discontinuity arises. Analogously, at a uniform free-stream, both forms of the LPPE reduce to a Helmholtz equation for the pressure perturbation.

Results
In what follows instability analysis results in four examples of wall-bounded and open flows will be presented, in which the temporal BiGlobal and TriGlobal eigenvalue problems will be solved using the LPPE boundary closure. The eigenspectra in all cases analyzed are well understood from a physical point of view, which facilitates assessment of the quality of results obtained.
Solutions to the real BiGlobal EVP which can be derived from (2-5) using a simple transformation [60] have been obtained subject to the two-dimensional LPPE boundary closure (12) in a spanwise homogeneous square lid-driven cavity [4,58]. On the other hand, flow driven by a constant pressure gradient in a duct of square cross-section [54,63] is an example of wall-bounded flow, instability of which is governed by a complex (two-dimensional) BiGlobal EVP. The open flow addressed is the classic Blasius boundary layer, in which again a real two-dimensional EVP may be solved. In the author's opinion, this flow is hardly appropriate as a test bed for application of global stability analysis methods, owing to the convective nature of its instability [27]. Nevertheless, the BiGlobal EVP has been solved by several authors with mixed degrees of success and it is interesting to contribute here to the related discussion in the literature, by separating issues arising due to the pressure boundary condition imposed at the wall boundary from those related with the open boundary treatment. Finally, an example of flow with three inhomogeneous spatial directions, which serves as demonstrator of applicability of the three-dimensional LPPE (13), is the cubic lid-driven cavity [15], instability of which is governed by the (three-dimensional) real TriGlobal EVP.
Eigenspectra obtained by the LPPE closure are compared with those delivered by two well-validated incompressible DNS codes, nek5000 [16] and nektar++ [10], both of which feature spectral element spatial discretization and time-stepping procedures for global instability analysis. In the case of the two-dimensional lid-driven cavity, the matrix-forming finite-element code FreeFEM++ is also employed to solve the BiGlobal EVP. In the Blasius flow, the nonlinear Parabolized Stability Equations (PSE) code NOLOT [24,26], a European industry-standard for the analysis of laminar-turbulent transition in boundary layer flows [52,55], is used to obtain amplitude function predictions, against which those obtained by the LPPE and the pressure compatibility boundary conditions are compared.

Matrix-forming solution of the real BiGlobal EVP on collocated grids
Linear instability of spanwise homogeneous flow in the square lid-driven cavity has been addressed by a number of authors [14,45]. Theofilis [58] and Albensoeder et al. [4] were the first to independently discover unstable eigenvalue branches subcritical to those known up to that time and present the critical conditions of the square [4,58] and rectangular [63] cavity. As mentioned, from a numerical point of view, in this wall-bounded flow the base flow develops on a plane perpendicular to the wavenumber vector, which permits reducing the complex global eigenvalue problem into one with real coefficients and thus halve the memory requirements for the storage of the discretized matrix [60]. From the point of view of the boundary closure based on the LPPE, the lid-driven cavity flow includes gradients of the base flow velocity componentsū(x, y) andv(x, y), which will be absent in the example to be considered in the next section, but their inclusion was found to be essential for the success of the present analysis. Instability in the square two-dimensional lid-driven cavity flow has been recently addressed by Taira and co-workers (personal communication), who obtained base flows using modified versions of the incompressible Cliff and the compressible CharLES 2 codes [9,22,23,31], while their stability analysis was performed by an in-house matrix-forming code [42]. Table 1 presents comparison of results already shown in reference [63], where a fully spectral algorithm was used for both the base flow computation and the solution of the eigenvalue problem. In addition, this table presents recently obtained (unpublished) results from the Taira group and it may be seen that the agreement with the reference results [63] is remarkable.
While the above results have been obtained using the singular version of the cavity so as to facilitate comparisons with earlier work, flow resulting from a regularized version of the lid motion is discussed next. In the regularized lid-driven cavity, the uniform lid motion of its standard (singular) counterpart is replaced by the two-dimensional version of the analytic function proposed by Leriche et al. [33], using p = 16. This choice delivers a base flow which is close to that of the singular configuration, while it facilitates exponential convergence of the spectral collocation method used to solve the two-dimensional equations of motion, the latter discussed in [60]. The PC equations (16) need to be modified to account for the nonzero base flow lid-velocity and gradients and becomê The real system which results from the linearized Navier-Stokes equations (2)(3)(4)(5) after the transformation discussed in [60] is used has been collocated on a two-dimensional grid of mapped Chebyshev-Gauss-Lobatto spectral nodes. Either of the LPPE (12) or the PC (15) was used to provide boundary closure. A reasonably large Krylov subspace dimension of 100 was used in a hard-coded version of the Arnoldi algorithm, while at the low Reynolds number considered, resolutions of 40 points per spatial direction were employed for both of the LPPE and PC results that will be discussed in Sect. 3.1.3.  Table 2 3

.1.2 Matrix-forming using FreeFEM++ and time-stepping using nektar++ and Semtex
In this section, instability analysis results have been obtained using two alternative methodologies, one based on the finite-element package FreeFEM++ [12,25] and one using time-stepping, as implemented in the (independent) direct numerical simulation codes nektar++ [10] and Semtex. All three of these codes have the capability to compute the base flow using their respective spatial discretization techniques, which are also employed for the solution of the eigenvalue problem. The latter is solved in FreeFEM++ by a matrix-forming approach, while in nektar++ and Semtex [7] the time-stepping algorithm discussed by Barkley et al. [5] is used. The FreeFEM++ based code used for the present instability analysis has been presented by Tammisola et al. [53] and was further validated by Lashgari et al. [32]. It uses an unstructured mesh comprising a total of 29,132 triangles and 14,787 vertices has been used to spatially discretize both the base flow of the regularized cavity (14) and the corresponding global eigenvalue problem. This implies a total number of degrees of freedom (and leading dimension of the matrix discretizing the EVP) of close to 2 × 10 5 . The matrix is stored in compressed format, and linear algebra operations are performed using the UMFPACK sparse matrix library [12]. The base flow is obtained using a Newton technique, while the eigenvalue problem is solved by calls to the ARPACK library using a Krylov subspace dimension of 30.
In nektar++, spatial discretization employed a structured mesh comprising 18 2 elements, clustered toward the cavity walls. In each direction within an element, a polynomial of degree 5 ≤ p ≤ 11 was used to ensure spatial convergence. The steady laminar base flow was computed by time-marching the equations of motion until convergence in time was obtained at t ≈ 150. Instability analysis was performed using the HalfMode option, corresponding to a real BiGlobal EVP. A hard-coded version of the Arnoldi algorithm with a Krylov subspace dimension of 80 was used for the iterative solution of the EVP, in which the unsteady linearized equations of motion were marched in time with a time-step of t = 10 −2 within 0 ≤ t ≤ 1, until a prescribed tolerance of τ = 10 −5 was reached in the desired eigenvalues. Figure 1 shows the leading members of the eigenspectra obtained using the LPPE, FreeFEM++, nektar++, Semtex and PC discussed in the previous two sections; the practically identical results obtained by nektar++ and Semtex have been represented by a single symbol. The most striking feature of this figure is the very  good agreement between the results of the first three methodologies. The respective algorithms underlying the solution of the real BiGlobal EVP use spectral collocation, second-order finite-elements, or spectral/ hp element spatial discretization and enforce continuity either strongly (LPPE) or weakly (FreeFEM++, nektar++ and Semtex). All approaches employ Arnoldi iteration to resolve the eigenspectrum in the neighborhood of ω = (0, 0) but use different sizes of the Krylov subspaces that they, respectively, construct. Despite these algorithmic differences, the three approaches produce an identical number of eigenvalues within any given radius, both at the present conditions and others not shown here for brevity. The leading five eigenvalues indicated as Mode I-V in Fig. 1 are presented in Table 2.

Comparisons
By contrast, the eigenspectrum obtained by imposition of the pressure compatibility equations (15) is seen to deviate substantially from the result of the other three approaches. While the leading part of the LPPE eigenspectrum is converged at the 40 2 resolution utilized, its PC counterpart is still sensitive to resolution. Visual inspection of the amplitude functions of the respective leading eigenmodes, shown in Fig. 2, verifies the non-physical nature of the least damped member of the PC eigenspectrum and provides a hint as to the origin of the discrepancy. Oscillations near the top and downstream cavity walls are visible in thep amplitude function delivered by the PC closure, although the overall field shape of this eigenfunction is similar to the pressure perturbation delivered by the LPPE (12). The latter amplitude function is smooth at the interior and, more importantly, at all of the stationary and moving cavity walls.
It can also be seen in Fig. 1 that some eigenvalues obtained by the PC closure are close (Modes II and V) or even agree very well (Mode III) with the correct results obtained by the other three methodologies. If resolution is increased, agreement improves further, and it could be speculated that the origin of the numerical artifacts appearing in the PC eigenspectrum is related with the strong gradients of the base flow at the corners of the moving lid. However, analysis of the PC boundary closure is beyond the scope of the present work, the aim of which is to establish the LPPE closure for the solution of the global eigenvalue problem on collocated grids. The agreement of the LPPE results for the real BiGlobal EVP with those of FreeFEM++, nektar++ and Semtex certainly points in that direction; however, the performance of (12) in solving the complex BiGlobal EVP or that in open flow domains also needs to be assessed; this is done next.
3.2 Instability of flow in a square duct

Matrix-forming solution of the complex BiGlobal EVP
Pressure-gradient driven flow in a rectangular duct is representative of the class of flows in which the basic state comprises one steady velocity component,w(x, y), along the homogeneous spatial direction, aligned with the wavenumber vector βe z . The eigenvalue problem governing stability of this flow is complex and was first solved in the classic work of Tatsumi and Yoshimura [54] as a generalization of the Orr-Sommerfeld equation in two inhomogeneous spatial dimensions.
Here, the two-dimensional eigenvalue problem is solved by collocating the linearized Navier-Stokes and continuity equations (2-5) on a non-staggered grid and imposing the two-dimensional LPPE (12) to provide boundary conditions for pressure. Results are compared with those obtained by imposition of the pressure compatibility (PC) boundary conditions which are obtained by taking both basic flow and perturbation velocity components to vanish at the wall, such that (3-5) reduce tô Equations (16) are collocated at the domain boundaries x = ±1 and y = ±1, respectively, as discussed by Theofilis et al. [63]. In both matrix-forming solutions, a spectral collocation discretization using upwards of 30 points per spatial direction is sufficient to converge the respective results.

DNS and time-stepping solution of the BiGlobal EVP
In order to aid identification of physically relevant EVP results from numerical artifacts introduced by the boundary conditions, direct numerical simulations have been performed for the square duct base flow and subsequently time-stepping solutions of the BiGlobal EVP have been obtained. Numerically obtained base flows were compared with the analytically known [63], Reynolds number independent, spanwise component of the steady laminar base flow, in x ∈ [−A, A] × y ∈ [−1, 1]. Choosing a periodicity length of z ∈ [0, 2π] along the spanwise direction, the associated constant pressure gradient is also known in closed form Direct numerical simulations have been performed at Re = 100 using both the nektar++ [10] and nek5000 [16] incompressible spectral element codes in a three-dimensional duct of square cross-section, A = 1. The length of the domain, L z = 2π, along the homogeneous spatial direction, z, has been chosen so as to match the fundamental wavenumber, β = 2π/L z = 1, at which instability analysis has been performed. In the nektar++ code, a pressure gradient along the spanwise direction, ∂P/∂z = 2/Re, has been imposed, while in the nek5000 code a constant volume flow rate Q was used to drive the flow. The linear relationship between Q andw has been verified numerically by Tatsumi and Yoshimura [54], who also presented the value Q = 1.1264 corresponding tow(x = 0, y = 0) ≈ 0.58937 that (17) delivers at A = 1. Both of the present simulations were validated against these results, in terms of recovery of the midpoint spanwise velocity value, an inflow pressure ofP inf = 0.12566 that is consistent with (18) when an outflow pressureP out = 0 is imposed, and a perfectly linear pressure decay along z. Subsequently, simulations were set up such that w(x = 0, y = 0) = 1, consistently with earlier analyses [54,63]. To this end, a volume flow rate Q = 1.90816 was used in nek5000, while in nektar the above-mentioned midpoint velocity value was used to scalew to unity, leading toP inf = 0.21321. In both cases, simulations started with fluid at rest and the unsteady equations Instability analysis in a time-stepping context did not employ nek5000, since Fourier expansions are not implemented in this fully three-dimensional code, which solves the present BiGlobal as a three-dimensional eigenvalue problem. By contrast, in the nektar++ package that was used for the present analysis the homogeneous spatial direction, z, is expanded using a small number of 4 − 8 Fourier modes and the stability analysis is performed using the FullMode option, which solves the complex BiGlobal EVP. Additional parameters utilized in the time-stepping context at this low Re have been a rather short integration time, T = 1, and large value of the time-step, t = 0.025, and a Krylov subspace dimension of 40. Convergence of the results was achieved in less that 500 iterations.
The eigenspectra obtained by closing the governing equations with either of the LPPE or the PC boundary conditions in a matrix-forming context, as well as the leading eigenvalues obtained by nektar++ close to the origin, are shown in Fig. 4. In all three sets of results, all but the strongest damped members of the respective eigenspectra are converged in their eigenvalues. Of interest in this figure are the following aspects. First, excellent agreement can be seen between the results of the LPPE matrix-forming and the nektar++ timestepping methodologies for the solution of the complex BiGlobal EVP. The eigenvalues marked as Modes I-V on this Figure are found in Table 3. Second, the results of the PC closure are in agreement with those of LPPE and nektar++ for most of the eigenvalues, including the least damped eigenmode recovered by the LPPE and nektar++, indicated as Mode I. However, what appears to be a branch exclusive to the PC results is also visible and actually contains eigenvalues that are less damped than Mode I of Table 3; in the absence of independent confirmation of the results of LPPE, the least damped member of the PC eigenspectrum could have been mistaken as the least damped flow eigenmode.
Visual inspection of the amplitude functions corresponding to the least damped eigenmodes obtained by the LPPE and PC closure methodologies, shown in Fig. 5, reveals that the "additional" eigenmode delivered by the pressure compatibility condition is spurious; the same holds for all other modes found in this eigenspectrum but not present in that delivered by the linearized pressure Poisson equation (12) or nektar++. By contrast, eigenmodes delivered by both of the PC and the LPPE closures are well resolved and correct. Consequently, the PC boundary closure could be used (and has been used in the literature, e.g., [63]) but care must be exercised in analyzing its results, since there exists no a priori means of distinguishing between physically relevant and numerically erroneous eigenmodes. In practice, one could inspect all of the physically most relevant-most amplified/least damped members of the eigenspectrum delivered by the PC condition; however, this task, besides not being satisfactory from a theoretical point of view, is also tedious to perform at each different set of parameters analyzed; use of the LPPE to close the eigenvalue problem is free from such additional post-processing, as it delivers only physically relevant eigenvalues.

The flat plate boundary layer [11]
An open flow example which is used to evaluate the performance of the LPPE boundary closure for global instability analysis is the classic laminar incompressible flow over a flat plate. From the outset, it should be noted that the author's point of view is that, from a physical perspective, instability analysis of the flat plate boundary layer can optimally be performed by local and non-local linear theory tools or direct numerical simulation. Nevertheless, the global eigenspectrum of the Blasius boundary layer has appeared several times in the literature, on occasion together with imposition of the PC boundary closure alongside claims of excellent agreement between two-dimensional eigenfunctions and results of the Orr-Sommerfeld equation or the linear Parabolized Stability Equations (PSE) [6,24,26]. It was thus found interesting to highlight the differences between the results of the LPPE and PC boundary closures in the context of this open flow configuration. Lengths are scaled with the (finite) plate half-width, and the base flow has been obtained by incompressible two-dimensional direct numerical simulation [11]. The two-dimensional (β = 0) temporal BiGlobal eigenvalue problem is solved in a rectangular domain defined by x ∈ [100, 650]× y ∈ [1,36], corresponding to a Reynolds number Re = 2400 at the beginning of the domain. The matrix-forming approach discussed by Paredes et al. [44] is used, in which equations (2)(3)(4)(5) are discretized by a 16th-order FD-q finite-difference method using N x = 601 and N y = 101 nodes along the x− and y−directions, respectively. The eigenspectra delivered by the LPPE and the PC closures are shown in Fig. 6, where it can be seen that, as in the closed flow cases previously discussed, substantial differences exist. The PC spectrum is actually more unstable than that of the LPPE and also contains unstable branches than are altogether absent in the LPPE result. As is well known for this flow [2,48], these eigenspectra are not converged. Instead, they are both composed of what is commonly referred to in the literature as box modes, the location of which shifts in parameter space, depending on resolution and size of the discretized domain [48]. The reader is warned that both in matrixforming and in time-stepping solutions of the two-dimensional eigenvalue problem, spurious modes related with insufficient domain length may also be present in the eigenspectrum [5]; this is not an issue in the present example, in which a large number of Tollmien-Schlichting wave periods has been included in the calculation domain. As is also known from non-modal analyses of the Blasius boundary layer [39,46,51,66], complete description of the physics of linear instability in this flow requires adopting a non-modal approach, which has been introduced in a global context by Abdessemed et al. [1]. This approach will not be followed here, since the present focus is on performance of the boundary conditions to close the eigenvalue problem system, a First, the amplitude functions of two members of the LPPE and PC spectra which correspond to two close-by eigenvalues, ω LPPE = 0.121308 + 0.00739111i and ω PC = 0.117821 + 0.00903062i are monitored in terms of the overall shape of each component of the respective eigenvector, as well as in terms of profiles at selected locations. The former result is shown in Fig. 7 where one may notice that, although the results are qualitatively analogous, the LPPE eigenfunction persists at substantially larger distances along the plate; also, oscillations are present in the far-field part of the amplitude functions delivered by the PC boundary closure, as have been seen in both closed flow examples analyzed in the previous sections.
Second, cuts through the LPPE and PC amplitude functions at an arbitrarily chosen location along the plate, x = 370, are presented as profiles in Fig. 8. Shown are the normalized streamwise,û(y), and wall-normal, v(y), profiles close to the wall, where noticeable differences can be seen; these may be explained by the slightly different eigenvalues ω LPPE and ω PC . However, as the free-stream is approached (not shown) the LPPE profile decays smoothly while the PC result shows oscillations that increase in amplitude, as can be appreciated in the results of Fig. 7. The question that these comparisons pose is: which global instability analysis result is closer to that of Blasius theory?
This question is answered by the third kind of comparison performed, between predictions of the LPPE boundary closures against results of linear PSE computations. The reasoning behind this comparison lies in the results of the celebrated work of Bertolotti et al. [6], in which PSE predictions of the spatial evolution of Tollmien-Schlichting (TS) waves developing in the incompressible flat plate boundary layer were shown to be in excellent agreement with results delivered by independently performed spatial direct numerical simulation of laminar-turbulent transition, in both of the linear and nonlinear regimes. In the PSE work performed here, the European aeronautics industry-standard NOLOT [24,26] code has been utilized.
The PSE space marching procedure is started at x = 50, using the (complex) eigenvalue ω LPPE = 0.121308 + 0.00739111i that the LPPE analysis has delivered. The equations are marched until an arbitrary streamwise location, x = 370, where the streamwise perturbation velocity of the TS wave predicted by PSE is compared withû(x = 370, y) extracted from the LPPE eigenvector shown in Fig. 7; excellent agreement between the moduli of the respective (complex) streamwise perturbation velocity components is It has to be stressed at this point that the intention of the manuscript is not to claim that use of the LPPE can deliver converged spectra in open flows in general, or in the Blasius boundary layer in particular. As stated at the beginning of this section, the author does not view solution of a two-dimensional eigenvalue problem as the appropriate methodology to perform analysis of instability in the flat plate boundary layer. The pitfalls of solving a two-dimensional BiGlobal EVP in boundary layer flows have been discussed by Rodríguez and Theofilis [48], who demonstrated the sensitivity of results of 2D EVP solutions to inflow/outflow boundary conditions. However, in literature relating to instability analysis of the Blasius boundary layer claims can be found of excellent agreement between results delivered by the Orr-Sommerfeld equation, (classic) PSE and the two-dimensional BiGlobal EVP. Such claims can be misleading, especially when no mention is made of the wall-boundary conditions imposed. The work presented in this manuscript demonstrates by comparison with independently performed PSE analyses, that, as far as the wall-boundary condition is concerned, results obtained by the PC are inferior to those delivered by the LPPE.
A further comment that can be made on open flow eigenspectrum convergence is that either the PSE, as an initial value problem, or the BiGlobal EVP, as one in which inflow and outflow boundary conditions need to be specified, will deliver results as good (or as bad) as the inflow (for the PSE) or as the inflow and outflow (for the 2D BiGlobal EVP) boundary conditions imposed. It is these conditions that lead, or not, to convergence of the spectra in open flows. The LPPE wall-boundary condition, subject of the present manuscript, deals with the ability to obtain accurate solutions of the LNSE when a wall is present in the flow and not with the ability to converge an open flow eigenspectrum. Of course, when all inflow/outflow/far-field/wall-boundary conditions have been set appropriately, converged spectra should be expected, as has been demonstrated in the closed flows discussed in this manuscript.

The 3d lid-driven cavity
The cubic lid-driven cavity permits assessing the performance of the three-dimensional LPPE (13) for modal TriGlobal linear stability analysis. Albensoeder and Kuhlmann [3] have provided accurate predictions for the steady three-dimensional base flow while linear instability has been addressed successfully numerically by Giannetti et al. [18], Feldman and Gelfgat [15] and Gómez et al. [19], among others, and experimentally by Liberzon et al. [34]. In the present work, the domain is defined in = [0, 1] 3 and flow is set up by horizontal motion of the wall y = 1 along the positive Ox axis direction. In order to preserve spectral accuracy, in a manner analogous to (14), a regularization of the lid motion using p = 16 has been implemented, [33] u(x, y = 1, z) Base flow computations have been performed with nek5000 [16] to time-march the flow to steady state, obtained after ≈ 100 and 250 time units at Re = 200 and 1000, respectively. The cubic domain was discretized by 10 elements per spatial direction, each of which was resolved with polynomials of degrees varying from 5 to 11. Figure 10 displays the convergence of residuals, defined by The three-dimensional eigenvalue problem is solved by the matrix-forming approach discussed in Gómez et al. [20]. First, using a modified version of the nek5000 internal routine int_tp, the base flow is interpolated spectrally from the spectral element grid onto a cubic domain discretized by 101 mapped Chebyshev Gauss Lobatto points in each spatial direction. Subsequently, the interpolated base flow is again interpolated spectrally onto the mesh used for the eigenvalue problem solutions, the latter built using 4th-and 6th-order FD-q methods [44] and comprising from 31 3 to 61 3 discretization nodes. The system (7-10) is then solved subject to homogeneous Dirichlet boundary conditions for the disturbance velocity components on all boundaries and pressure perturbation boundary conditions provided by the three-dimensional LPPE (13).   [18] ±0.284 −0.1356 Gómez et al. [19] ±0.285 −0.1360 Liu [37] ±0.299 −0.1353 Shown are also results obtained by nek5000 in the regularized cavity, as well as available literature data in the singular configuration Table 4 presents the convergence history of the leading (damped) eigenmodes at Re = 200 and Re = 1000, where it can be seen that convergence of five significant places in most of the eigenvalues has been reached in the LPPE results. One also observes that the monotonic decay of the perturbation shown at Re = 200 in Fig. 10 corresponds to the least damped stationary mode recovered by the LPPE boundary closure, ω i = −0.45073 at this Reynolds number. The relative difference of this eigenvalue and the slope of the straight line computed by post-processing the nek5000 DNS results is lower than 0.02%. Table 4 also presents the converged results of time-stepping solution of the TriGlobal EVP, obtained by nek5000. The base flow has been interpolated on grids featuring between 4 3 and 8 3 elements, which are resolved by polynomials of degrees 5 ≤ p ≤ 9. The At the higher Reynolds number value examined, Re = 1000, the results delivered by the LPPE boundary closure, new results obtained by the nek5000 code in the regularized cavity, as well as results obtained by Giannetti et al. [18], Gómez et al. [19] and Liu [37,38] are also presented in Table 4. Here it may be observed that the perturbation being damped in the nek5000 DNS results, also shown in Fig. 10, is a traveling mode and has a damping rate that can be predicted by the slope of the residual and the stability analysis as ω i ≈ −0.1402. On the other hand, despite the fact that work performed in [18,19,37,38] has considered the singular configuration, the distances in parameter space between eigenvalues predicted by any of the three previous works and the present analysis are of the same order of magnitude for all three modes. Perhaps the only qualitative difference between the instability analysis results in the singular and the regularized lid-driven cavity at Re = 1000 is the fact that the leading eigenmode in the latter configuration is a stationary disturbance, as opposed to the traveling mode being the least damped in the singular cavity.
The amplitude functions of the leading stationary mode at Re = 1000 obtained by the LPPE closure are shown in Fig. 11 where, besides the rather well-known velocity perturbations, also the amplitude function of the pressure perturbations is shown. Isosurfaces very close to zero have been included in the presentation of thep(x, y, z) amplitude function in order to highlight the degree of complexity expected in this component of the eigenvector as the Reynolds number increases. This is a consequence of the gradients of the base flow at Re = 1000 in conjunction with the coupled base flow derivatives appearing in the right-hand side of (13) and implies that the base flow monitored should be sufficiently well represented in order for reliable results of the eigenvalue problem (7-10) to be obtained. As a matter of fact, in all analyses performed the pressure amplitude function has been found to be a reliable diagnostic tool to determine the quality of the base flow analyzed; converged instability analysis results can only be obtained when this component of the eigenvector is well resolved.

Summary
Two forms of the linearized pressure Poisson equation (LPPE), one pertinent to BiGlobal and one to TriGlobal linear instability analysis, have been derived. These equations have been used to provide boundary closure for the pressure perturbation in the respective systems of linearized Navier-Stokes equations, the latter solved as an eigenvalue problem in a matrix-forming context on collocated meshes. Steady laminar flows in closed domains have been analyzed in two spatial dimensions, where results of the LPPE boundary closure have been compared with those obtained by the FreeFEM++ matrix-forming and the nektar++ and Semtex timestepping codes; excellent agreement was obtained in all three classes of flow instability problems considered. The superior performance of the LPPE over the more commonly used PC boundary closure was demonstrated in terms of the absence from the LPPE results of spurious modes found in the PC eigenspectrum. The structure of the three-dimensional pressure perturbation amplitude function was seen to provide a diagnostic tool to assess the quality of the results obtained by the pressure compatibility closure. In the flat plate boundary layer, stability of which is discussed from a physical point of view by Saric et al. [50], results of the two-dimensional eigenvalue problem obtained by the LPPE boundary closure were seen to agree very well with PSE space marching solutions of the linearized perturbation equations pertinent to this non-parallel base flow. Finally, the regularized cubic lid-driven cavity was analyzed using the LPPE boundary closure and results were obtained in relatively close agreement with existing literature on the singular counterpart of this flow as well as with those obtained herein by the time-stepping module of the nek5000 code. The results presented, and others not shown here for brevity, demonstrate that the LPPE provides reliable pressure perturbation boundary conditions for the solution of the BiGlobal and TriGlobal eigenvalue problems in a matrix-forming context on collocated grids.
Acknowledgements Prof. K. Taira and his co-workers have provided the comparison results included in Table 1. The LPPE and PC boundary closures have been implemented in the codes described in the Ph.D. Thesis of Dr. P. Paredes [43]. Results on the flat plate, section 3.  Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.