On Spectral Approximations for the Stability Analysis of Boundary Layers

Abstract

This paper devotes to the approximation of spectral and boundary-value problems arising in the stability analysis of incompressible boundary layers. As an alternative to the collocation method with mappings, the Galerkin–collocation method based on Laguerre functions is adopted. A robust numerical implementation of the latter method is discussed. The methods are compared within the stability analysis of the Blasius and Ekman layers. The Galerkin–collocation method demonstrates an exponential convergence rate for scalar stability characteristics, and has a number of advantages.

INTRODUCTION

Three-dimensional small-amplitude disturbances against a main laminar flow are of interest in numerical studies of boundary-layer instabilities. Equations governing the evolution of such disturbances are considered on the half-line , where is the wall-normal coordinate, with the boundary co-nditions

1.1

for the velocity components , and (see, e.g., [1–3]). The boundary conditions (1.1) represent the no-slip condition at the flow-exposed surface and decaying of disturbances at far distance from the surface. This paper devotes to approximations on of such problems.

Spectral methods, including collocation and Galerkin–collocation methods, are a good choice for the approximation of governing equations, since the equations are linear while sought disturbances are smooth functions of . Within the collocation method, the solution is approximated by a series of infinitely-differentiable functions with a non-finite supply; and the expansion coefficients are determined by requiring the equations to be satisfied at given grid nodes called collocation points. Within the Galerkin–collocation method, the equations are approximated in the weak form, Lagrange interpolation functions associated with some grid are used as trial and test functions, and the inner products are computed by a high-order quadrature formula associated with the same grid. Note that the Galerkin–collocation method is often called the Galerkin method with numerical integration [4]. For problems considered on a finite interval, these spectral methods are discussed, for example, in [4–6]; and procedures from the well-known software packages [7, 8] can be used for the numerical implementation of these methods.

There are three main approaches for approximating problems considered on the half-line under the boundary conditions (1.1). The first one introduces an artificial boundary at finite but large distance from the surface. Then the equations are considered on the interval under some (e.g., zero or asymptotic [9]) boundary conditions at instead of those at infinity. Coupled with a spectral method for the finite interval, this approach is widely used in numerical studies of boundary-layer instabilities (see [3] and references therein). This approach requires choosing the sufficient value of for a particular problem, and might depend on flow parameters. Note that boundary conditions for the velocity components at any might allow for solutions that do not decay as . For boundary-layer stability problems, solutions with such a behavior are known; these solutions correspond to the modes of continuous spectrum [1, 3], with their physical relevance being still an open question.

The second approach uses a mapping that transform a system of functions with well-established approximation properties on a finite interval (for example, the Lagrange interpolation polynomials associated with the Chebyshev points) into that on a half-line [4, 10]. The approximation properties of such mapped systems of functions are discussed in [11]. From [12, 13] onwards, various mappings are compared for model problems. This approach is used in hydrodynamic and aerodynamic applications [10, 14, 15].

The third approach [11, 16] uses the Laguerre functions , where is the Laguerre polynomial of degree . As to our knowledge, spectral methods based on the Laguerre functions have not been previously used for studying boundary-layer instabilities.

In [11] the convergence of the spectral Galerkin method based on either mapped systems of polynomials or the Laguerre functions is studied theoretically. Upper bounds on the approximation errors are obtained for both type of methods and then verified on model elliptic equations. Note that these bounds are obtained in different norms. For the method based on the Laguerre functions, the norm is the usual (i.e., with the unit weight function) -norm over the half-line; this norm has a clear physical interpretation in the study of boundary-layer instabilities—it is the disturbance kinetic energy density. For the method based on mapped systems of polynomials, that is the weighted norm determined by the mapping. The work [11] provides a number of examples, where either the first method converges faster than the second one, or the second one converges faster than the first one, or both methods show close results. It is of interest to compare these methods for boundary-layer stability problems.

The present work is organized as follows. In Section 2, we describe the approximation by the Galerkin–collocation method based on the Laguerre functions of the equations governing evolution of small-amplitude disturbances of viscous incompressible boundary layers. Section 3 devotes to a robust numerical implementation of this method. Section 4 compares the proposed method with the collocation method with mappings for the stability analysis of the Blasius and Ekman layers. Section 5 summarizes the results.

Throughout this paper, denotes the -norm for vectors and matrices, the superscripts and denote the symbols of transposition and conjugate transposition respectively, and denotes the Kronecker delta.

APPROXIMATION OF PROBLEMS ARISING WITHIN THE STABILITY ANALYSIS OF BOUNDARY LAYERS

In the Cartesian coordinates, (streamwise), (wall-normal) and (spanwise), consider a flow of a viscous incompressible fluid over the flat surface . Against the background of a main laminar flow, we consider three-dimensional small-amplitude time-dependent disturbances which are represented as follows

2.1

where , , , and are the complex-valued amplitudes of the streamwise, wall-normal and spanwise velocity components, and the pressure, respectively. The amplitudes depend only on and . Here is the time, is the streamwise wavenumber, and is the spanwise wavenumber.

Two problems are considered in this paper to present and compare approximation methods in the wall-normal direction .

The first problem is the temporal stability analysis of the Blasius layer under the local-parallel assumption. In this case, it is assumed that the linear dimensionless equations governing the evolution of small-amplitude disturbances are as follows

2.2

where is the Reynolds number, and . The streamwise velocity of the main flow depends only on ; and satisfies the Blasius equationwhich can be solved by standard numerical methods (see, e.g., references in [3]). A physical interpretation of Eqs. (2.2) as well as the definition of the Reynolds number can be found in [2, 3].

The second problem is the temporal stability analysis of the Ekman layer. In this case, it is assumed that the linear dimensionless equations governing the evolution of small-amplitude disturbances are as follows

2.3

where is the Reynolds number, and is the Rossby number. The streamwise velocity and spanwise velocity of the main flow depend only on and are known analytically. A physical interpretation of the Eqs. (2.3) as well as the definition of the Reynolds and Rossby numbers can be found in [14].

Within the stability analysis, the velocity components satisfy the boundary conditions (1.1) for both considered problems. In addition, we consider the disturbance kinetic energy density

2.4

as a physically-relevant measure of disturbance magnitude.

Approximation by the Galerkin–Collocation Method Based on the Laguerre Functions

Let us consider Eqs. (2.2) under the boundary conditions (1.1). Suppose is the space of complex-valued functions square-integrable over the half-line . This space is equipped with the inner product and the norm that is similar to the energy functional (2.4). Suppose is the space, whose elements satisfy zero boundary condition at and belong to together with their first derivatives. Let us multiply the momentum equations by and the continuity equation by ; and integrate these equations over the half-line , using the integration by parts. Thus, we obtain the weak form of Eqs. (2.2). We seek for and (at any fixed ) such that the weak form of Eqs. (2.2) is valid for any and .

Let be the Laguerre polynomial of degree . Suppose is the Laguerre–Gauss–Radau grid, whose non-zero nodes are the roots of the derivative of . The Laguerre–Gauss–Radau quadrature formulaassociated with this grid is exact for any polynomial of degree or less [17]. Then, the following quadrature formula is valid

2.5

Suppose are the Lagrange interpolation polynomials for the Laguerre–Gauss–Radau grid. Likewise, are those for the grid . It is easy to see that

2.6

In the sequel, functions of the form

2.7

are called the Laguerre interpolation functions. The functions at equal zero at . These functions are used as trial functions for the velocity components, and as test functions for the momentum equations. The functions are used as trial functions for the pressure, and as test functions for the continuity equation. The quadrature formula (2.5) is used for computing the inner products.

Let us point out the approximation of some operators in the weak form of (2.2). Let be the approximation of a function from by (), and be the column whose elements are the corresponding expansion coefficients. Likewise, let be the approximation of a function from by (), and be the column whose elements are the corresponding expansion coefficients. Then, the following equalities are valid

2.8

2.9

where the Euclidean inner product is denoted by the braces, is the diagonal matrix of order whose entries are the quadrature weights (2.5), is the matrix of size whose entries are the derivatives of () at the Laguerre–Gauss–Radau nodes, and is the matrix of size whose entries are the values of at the Laguerre–Gauss–Radau nodes. Note that Eqs. (2.8), (2.9) hold since the quadrature (2.5) is exact for any function of the form , where is a polynomial of degree or less. The matrix is called the differentiation matrix. The matrix is called the projection matrix; only have to be computed since at the interior nodes by definition. The computation of the matrices and is discussed in Section 3.

As a result of the described approximation of Eqs. (2.2), we obtain a system of ordinary differential and algebraic equations of the form

2.10

where is the -component column, whose elements represent the values of the velocity components at the interior grid nodes. In (2.10), is additionally scaled such that correspond to the energy functional (2.4). Here , , and are matrices of size , , and , respectively. The ordinary differential equations in (2.10) correspond to the momentum equations, while the algebraic equations correspond to the continuity equation. From (2.8) and (2.9), it is easy to see that the discrete analogue of the Laplace operator is a symmetric negative-definite matrix as well as equality is valid.

Similar to the polynomial interpolation, approximation properties of the functions (2.7) are determined by the Lebesgue function and the Lebesgue constant

2.11

At given , the function is equal to at and decays exponentially at . Figure 1 shows the function at ; it is qualitatively the same at other . As for the polynomial interpolation, the Lebesgue constant increases with . Figure 1 shows that the increase of is at logarithmic rate. In addition, the increase of is compared to the increase of the Lebesgue constant for polynomial interpolation at the Chebyshev points. It is shown that the values of are slightly smaller, while the growth rate is similar.

[See PDF for image]

Fig. 1.

On the left: the Lebesgue function (2.11) at . The Laguerre–Gauss–Radau nodes are marked with green dots. On the right: the increase of (green solid) with , and that of the Lebesgue constant for the polynomial interpolation for the Chebyshev grid (black dashed).

Scaling for the Stability Analysis of Boundary Layers

Within the boundary-layer stability analysis, there are two characteristic wall-normal length scales – the thickness of the laminar boundary layer , and the finite height such that disturbances might be regarded as negligible at . The value of can be found before the stability analysis, using only the main flow data (see introduction in [18] and references therein). In contrast, the value of can be found only within the stability analysis by studying the convergence of the sought disturbances with increasing . We also note that some disturbances (e.g., Tollmien–Schlichting waves) can extend significantly above the boundary layer, i.e., can be much larger than . This physical discussion leads to the following requirements. The grid nodes should be separated such that both the boundary-layer domain () and its outside are covered. In addition, as increases, both the number of nodes inside and outside the boundary layer must increase.

Therefore, the Laguerre–Gauss–Rado grid should be scaled, since these nodes are distributed along the entire half-line , with increasing with . For example, one can satisfy these requirements by the scaling

2.12

where is a scaling factor, and is the integer part of . This scaling means that the half of the grid nodes lies inside the boundary layer. An advantage of (2.12) is that this scaling does not depend on . Note that this is not the only possible way of scaling. For example, one can additionally adjust the parameter for a particular problem.

Figure 2 shows the Laguerre–Gauss–Radau nodes and weights under the scaling (2.12) at various with being fixed. Note that increases slowly with and decreases with under the scaling (2.12).

[See PDF for image]

Fig. 2.

On the left: the Laguerre–Gauss–Radau nodes under the scaling (2.12) at fixed and at , where . The independent variable is stretched along the vertical axis. On the right: the quadrature weights (2.5) under the scaling (2.12) at the same and .

Approximation by the Collocation Method with Mappings

Let us briefly describe the approximation of Eqs. (2.2) under the boundary conditions (1.1) by collocation method.

Let be the Chebyshev points, i.e., . Suppose are the Lagrange interpolation polynomials for this grid, and are those for the grid . Let () be a smooth monotonic function that ensures an one-to-one mapping between the interval and the half-line such that , , and . Such a mapping gu-arantees that the half of the grid nodes lies inside the boundary layer, similarly to the scaling (2.12). Then we use the functions as basis functions for the velocity components and the functions as basis functions for the pressure. The functions at satisfy zero boundary conditions at and . The approximation properties of such functions are discussed in [11].

For computing the energy functional (2.4), we use the quadrature formulawhere , and are the weights of the Clenshaw–Curtis quadrature [6]. This formula is exact for any functions of the form , where is a polynomial of degree or less. The values of the derivatives of and at the grid nodes can be computed by the procedures from [7] or [8], coupled with the chain rule.

As a result of the described approximation of Eqs. (2.2), we obtain a system of ordinary differential and algebraic equations of the form (2.10). Note that the collocation method does not ensure that the discrete analogue of the Laplace operator is a symmetric negative-definite matrix. In addition, the equality does not hold, in general.

As a mapping, the following ones are used in the present paper

2.13

2.14

2.15

The algebraic (2.14) and exponential (2.15) mappings are known [4]; and the mapping with the tangent function (2.13) is currently implemented in software package [15], which is designed for predicting an onset of laminar–turbulent transition in industrial applications.

Figure 3 shows the streamwise velocity of the Blasius layer and the streamwise and spanwise velocities of the Ekman layer. For both considered main flows, the typical values of (see, [9, 14]) are marked, while the values of correspond to the upper limits of the subfigures. In addition, Fig. 3 shows the Laguerre–Gauss–Radau nodes under the scaling (2.12), and the grid nodes obtained by either the mapping (2.13), (2.14), or (2.15).

[See PDF for image]

Fig. 3.

The independent variable is along the vertical axis. On the left: the Blasius velocity profile (, ). On the right: the Ekman velocity profiles, and (, ). The Laguerre–Gauss–Radau nodes under the scaling (2.12) (green dots) and the Chebyshev points under the mappings (2.13)–(2.15) (red, blue and pink dots, respectively) at .

NUMERICAL IMPLEMENTATION OF THE GALERKIN–COLLOCATION METHOD

To implement the approximation method described in Section 2.1, one has to compute the nodes and weights of the quadrature formula (2.5), the derivatives of (2.7) at the grid nodes, and the values of (2.7) at . In this section, we provide an algorithm for computing these quantities; the proposed algorithm is stable, including the case of large .

By definition [19], the Laguerre polynomials are orthogonal in the inner productwith the exponential weight function. They satisfy the following three-term relations

3.1

and can be represented as

3.2

In addition, the Laguerre polynomials satisfy the relations

3.3

3.4

whose derivation from (3.1) and (3.2) is straightforward, see [20].

It is known [17] that the Laguerre–Gauss–Radau nodes are eigenvalues of the symmetric tridiagonal matrix

3.5

This allows for the robust computation of .

Some Results

The quadrature weights (2.5) are determined by , which can be computed by the three-term relations (3.1). At large and , the values of appear to be very large (up to the machine infinity), and the values of appear to be very small (up to the machine zero). Therefore, the stable computation of is an issue. At the same time, the values of are bounded from below since the Laguerre functions are bounded, , at any and [16].

We propose to compute the weights (2.5) by

3.6

with an additional scaling at computing by (3.1). If we have at some , where is a given threshold parameter, then we divide the previously computed and by and keep using the three-term relation. This scaling by is done whenever the result exceeds in absolute value. As a result, we have , where is the number of fractions done, and is the value obtained by the procedure. Numerical experiments show the overall procedure is robust to round-off errors at large and ; the computed values of up to are shown in Fig. 2.

Let us consider the numerical interpolation by the Laguerre interpolation functions (2.7). The following statement is valid.

Lemma 1.Suppose the functionis equal toat some grid. Then, the interpolantconstructed with the functions of the formwhereare the Lagrange interpolation polynomials for this grid, is represented as

3.7

where

Proof. It is straightforward that the following representation

3.8

is valid, where . Then, (3.8) andend the proof.

As for the polynomial interpolation (see, e.g., [6]), the representation (3.7) is called the barycentric form of the interpolant, while are called the barycentric weights. The barycentric weights can be computed in a robust way due to the following statement.

Theorem 1.Supposeare the Laguerre–Gauss–Radau nodes, andare the quadrature weights in (2.5). Then, the barycentric weightsfor this grid, and thosefor the grid, are represented as

3.9

3.10

where

Proof. Let us prove (3.10) first. The polynomial has the same roots as ; therefore, these polynomials differ only by a multiplicative factor. Using that the leading coefficient of is equal to , we obtain that

3.11

For the interior Laguerre–Gauss–Radau nodes, it is valid that

3.12

where the left equality is obtained by taking the first derivative of (3.4), and the right one follows from (3.3). Substituting (3.12) into (3.11) and using the definition of (2.5), we obtain the statement (3.10) up to a sign. To end the proof, note that have to change the sign, with .

To prove (3.9), note thatand therefore

3.13

For the interior grid nodes the first term in the denominator of (3.13) equals , while the second term is expressed in terms of ; those lead to (3.9) at . For the boundary node , we have . To end the proof, note that the Laguerre polynomials satisfy , and therefore .

Note that the statement similar to (3.9) is proven in [21] for the polynomial interpolation for the Laguerre–Gauss–Radau grid.

The barycentric weights and contain the multiplicative factor , which decays at very large rate with increasing . At some , this factor becomes smaller than the machine zero. However, there is no need to compute for the interpolant representation (3.7), since the barycentric weights are involved both in the nominator and denominator. Thus, the interpolant representation (3.7) with the barycentric weights computed by Theorem 1 allow for the numerical interpolation from the Laguerre–Gauss–Radau grid to another grid. In addition, substituting (3.10) into (3.7) at the point , the following corollary is obtained.

Corollary 1.Supposeare the interior Laguerre–Gauss–Radau nodes,are the quadrature weights in (2.5) corresponding to these nodes, andare the Laguerre interpolation functions (2.7) for this grid. Then,

One can also derive the explicit formula for the derivatives of the interpolation Laguerre functions (2.7) at the Laguerre–Gauss–Radau nodes. Note that such formula is given in [16], Eq. (3.17), but without a proof.

Theorem 2.Supposeare the Laguerre–Gauss–Radau nodes,are the quadrature weights in (2.5), andare the Laguerre interpolating functions (2.7) for this grid. Then

Proof. The derivative of the function (2.7) is as follows

3.14

where are the Lagrange interpolation polynomials. By (2.6), it is valid that

3.15

By taking the derivative of (3.15), we obtain thatat . Substituting this expression to (3.14) and using the result for the barycentric weights (3.9), we prove the theorem at .

By taking the second derivative of (3.15), we obtain thatAt , it is valid that and . At other nodes , it is valid that (3.3). Thus,

3.16

The substitution of (3.16) into (3.14) ends the proof at .

To sum up, the Galerkin–collocation method based on the Laguerre functions can be implemented as follows. The Laguerre–Gauss–Radau nodes are computed as the eigenvalues of (3.5). The quadrature weights associated with this grid are computed by (3.6) with an additional scaling while using the three-term relations (3.1). Then, the values and derivatives of the functions (2.7) at the grid nodes are computed, see Corollary 1 and Theorem 2. For the interpolation of a grid function given at the nodes to another grid, the barycentric formula (3.7) is used, where the barycentric weights are computed by Theorem ; the multiplicative factor is common for all barycentric weights at given , therefore there is no need to compute it. The proposed algorithm performs robustly, including the case of large .

NUMERICAL EXPERIMENTS

As a result of the approximation of either the system (2.2) or (2.3) under the boundary conditions (1.1) by either the Galerkin–collocation method from Section 2.1 or the collocation method from Section 2.3, we obtain a differential-algebraic system of the form (2.10).

Note that lies in the kernel of . Let be a rectangular matrix, whose columns form an orthonormal basis in the kernel of . Likewise, let be a rectangular matrix, whose columns form an orthonormal basis in the kernel of . Under an additional assumption that both , and are of full rank, the system (2.10) is equivalent to the system of ordinary differential equations

4.17

where , and . The detailed justification of such a reduction of differential-algebraic systems is given in [22]; the assumption made is valid for the considered problems. Note that the approximation by the Galerkin–collocation method ensures that , and therefore . It is also worth noting that is the discrete analogue of the energy functional (2.4).

Within the numerical stability analysis, the eigenvalue of with the largest real part is of the most interest [2, 3]. This eigenvalue is called the leading eigenvalue, and the corresponding eigenvector is called the leading eigenvector; we denote the leading eigenvalue by . The another quantity of physical interest [2, 3] iswhich is called the maximum energy amplification. The quantity represents the maximum possible growth of the dusturbance kinetic energy density at given time period . In case the matrix is non-normal, the value of might significantly exceed the growth of the leading eigenvector [2, 3]; in this case, the initial disturbance at which is achieved (which is called the optimal disturbance [2, 3]) usually differs from the leading eigenvector. The maximum energy amplification can be computed by the efficient matrix algorithm [23] based on a low-rank approximation; this algorithm guarantees the result with a given accuracy.

To compare the approximation methods, we study the convergence of both scalar characteristics and . Comparing the methods by the convergence of vector characteristics, namely either the leading eigenvector or the optimal disturbance, could not be done quantitatively due to additional errors from interpolation from one grid to another.

For the considered test problems, the stability analysis can be performed only numerically. To establish the referential values and , we set the artificial boundary with zero boundary conditions for the velocity components at and then approximate the equations by the Galerkin–collocation method with the Lagrange interpolation polynomials for the Gauss–Lobatto grid as trial and test functions. The approximation properties of these basis functions are well-established [4], while the method was widely used for hydrodynamic stability problems considered on finite domains. Therefore, this method is reliable that is the most important for obtaining referential values. For boundary-layer stability problems, this method is certainly inefficient, since the Gauss–Lobatto nodes are refined both to and , while the value of has to be tuned manually. Tracking the convergence of the referential solution by increasing and , we achieve the convergence of and up to a desired precision.

As the first test problem, we perform temporal stability analysis of the Ekman layer (2.3). The Ekman layer could be considered as the simplest model of atmospheric boundary layers, with its stability being studied in detail (see [14] and references therein). Using the results of [14], we choose the parameter values as Re = 500, , , , and , where , . The referential values of the leading eigenvalue and the maximum energy amplification are computed at and .

Figure 4 demonstrates the relative errors at computing and for various approximation methods. All methods show an exponential convergence rate of with increasing . However, only the Galerkin–collocation method based on the Laguerre functions shows an exponential convergence rate of . Note that converges slightly faster than for this method. It is also worth noting that there are no significant differences between collocation methods at various mappings. Additional experiments not presented in this paper show that these findings remain qualitatively the same at increasing or decreasing the Reynolds number.

[See PDF for image]

Fig. 4.

The relative error at computing (on the left) and (on the right) with increasing for the Ekman layer. Results for the Galerkin–collocation method based on the Laguerre functions with scaling (2.12) are marked with green. Those for the method used to obtain referential values are marked with black. Those for the collocation methods with mappings (2.13)–(2.15) are marked with red, blue and pink, respectively.

As the second test problem, we perform temporal stability analysis of the Blasius layer (2.2). This main flow could be considered as the simplest model of aerodynamic boundary layers, with its stability being studied in detail (see [2, 3] and references therein). Using the results of [9], we choose the parameter values as , , and for computing ; and the parameter values , , , for computing . The typical boundary-layer thickness is . The referential values of the leading eigenvalue and the maximum energy amplification are computed at and .

Figure 5 demonstrates the relative errors at computing and for various approximation methods. The collocation method shows an exponential convergence rate of with increasing ; and there are no significant differences between mappings used. In contrast, for the Galerkin–collocation method based on the Laguerre functions, the accuracy that can be achieved is limited; nevertheless, the obtained accuracy might be more than enough in applications. This issue is remedied by choosing a larger , that is also shown in Fig. 5. One can also improve the scaling (2.12) by increasing the share of the nodes outside the boundary layer.

[See PDF for image]

Fig. 5.

The relative error at computing (on the left) and (on the right) with increasing for the Blasius layer. The colors mean the same as for Fig. 4. In addition, results for the Galerkin–collocation method based on the Laguerre functions with scaling (2.12) are shown at larger (green dotted).

As for the Ekman layer, the Galerkin–collocation method based on the Laguerre functions shows an exponential convergence rate of , while the collocation method leads to a slower convergence of . This disadvantage of the collocation method appears to be irremediable; and the reason is an ill approximation of the operator , which leads to the presence of slowly damped unphysical solutions. The basis for this conclusion are as follows. First, it is observed that convergence to the referential value for the collocation method is from above, i.e., . Second, let us consider the discrete analogue of the operator on the half-line under zero boundary conditions at and ; denote this matrix by . The collocation method, in general, does not ensure the symmetry of , although in this case it provides negative definiteness. Nevertheless, the matrix resulted as the approximation by the collocation method has a large condition number (e.g., of order at at the mapping (2.13)), and the spectrum of contains several very small in absolute value negative eigenvalues corresponding to strongly oscillating eigenvectors. For comparison, has the conditional number of order at , when approximated by the Galerkin–collocation method based on the Laguerre functions.

SUMMARY

This paper proposes the Galerkin–collocation method based on the Laguerre functions for approximating spectral and boundary-value problems arising in studying boundary-layer instabilities. These problems considered on the half-line , where is the wall-normal coordinate. The robust numerical implementation of this method is proposed (see Section 3), including the procedure for computing the weights of the Laguerre–Gauss–Radau quadrature formula, the explicit expressions for values and derivatives of the Laguerre interpolation functions at the grid nodes, and the procedure for numerical interpolation from the Laguerre–Gauss–Radau grid to another grid.

Within temporal stability analysis of the Blasius and Ekman layers, the proposed method is compared to the collocation method with mappings; the latter method is often used for numerical analysis of boundary-layer instabilities. The comparison is made at computing both the leading eigenvalue and the maximum energy amplification . It is shown that both type of methods show an exponential convergence rate of ; and differences between the methods are insignificant. However, the Galerkin–collocation method based on the Laguerre functions shows an exponential convergence rate of , while the collocation method leads to a slower convergence of this quantity. It is shown that converges faster than for the Galerkin–collocation method based on the Laguerre functions.

The Galerkin–collocation method based on the Laguerre functions might be successfully applied as the wall-normal approximation for more complex boundary-layer stability problems (see [24]).

FUNDING

The work is supported by the Russian Science Foundation (grant no. 22-11-00025).

CONFLICT OF INTEREST

The author of this work declares that he has no conflicts of interest.

Translated by the author

Publisher’s Note.

Pleiades Publishing remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

AI tools may have been used in the translation or editing of this article.

REFERENCES

1 Zhigulev, V. N.; Tumin, A. M. Origin of Turbulence: Dynamic Theory of Excitation and Evolution of Instabilities in Boundary Layers; 1987; Novosibirsk, Nauka: 0633.76001 [in Russian]

2 Schmid, P. J.; Henningson, D. S. Stability and Transition in Shear Flows; 2000; New York, Springer: 0966.76003

3 Boiko, A. V.; Dovgal, A. V.; Grek, G. R.; Kozlov, V. V. Physics of Transitional Shear Flows: Instability and Laminar–Turbulent Transition in Incompressible Near-Wall Shear Layers; 2011; Berlin, Springer: 1096.76001

4 Canuto, C.; Hussaini, M. Y.; Quarteroni, A.; Zang, T. Spectral Methods: Fundamentals in Single Domains; 2006; Berlin, Springer: [DOI: https://dx.doi.org/10.1007/978-3-540-30726-6] 1093.76002

5 Boyd, J. P. Chebyshev and Fourier Spectral Methods; 2000; 0994.65128

6 Trefethen, L. N. Approximation Theory and Approximation Practice; 2019; Philadelphia, SIAM: [DOI: https://dx.doi.org/10.1137/1.9781611975949] 1264.41001

7 Weideman, J. A. C.; Reddy, S. C. A MATLAB differentiation matrix suite. ACM Trans. Math. Software; 2000; 26, pp. 465-519.1939962 [DOI: https://dx.doi.org/10.1145/365723.365727] 1381.86002

8 T. A. Driscoll, N. Hale, and L. N. Trefethen, Chebfun Guide (2014).

9 Boiko, A. V.; Demyanko, K. V.; Nechepurenko, Yu. M. Asymptotic boundary conditions for the analysis of hydrodynamic stability of flows in regions with open boundaries. Russ. J. Numer. Anal. Math. Model.; 2019; 34, pp. 15-29.3909801 [DOI: https://dx.doi.org/10.1515/rnam-2019-0002] 1428.65072

10 Spalart, P. R.; Moser, R. D.; Rogers, M. M. Spectral methods for the Navier–Stokes equations with one infinite and two periodic directions. J. Comput. Phys.; 1991; 96, pp. 297-324.1128225 [DOI: https://dx.doi.org/10.1016/0021-9991(91)90238-G] 0726.76074

11 Shen, J.; Wang, L.-L. Some recent advances on spectral methods for unbounded domains. Commun. Comput. Phys.; 2009; 5, pp. 195-241.25136851364.65265

12 Grosch, C. E.; Orszag, S. A. Numerical solution of problems in unbounded regions: Coordinate transforms. J. Comput. Phys.; 1977; 25, pp. 273-296.488870 [DOI: https://dx.doi.org/10.1016/0021-9991(77)90102-4] 0403.65050

13 Boyd, J. P. The optimization of convergence for Chebyshev polynomial methods in an unbounded domain. J. Comput. Phys.; 1982; 45, pp. 43-79.650425 [DOI: https://dx.doi.org/10.1016/0021-9991(82)90102-4] 0488.65035

14 Foster, R. C. Structure and energetics of optimal Ekman layer perturbations. J. Fluid Mech.; 1997; 333, pp. 97-123. [DOI: https://dx.doi.org/10.1017/S0022112096004107] 0959.76026

15 Boiko, A. V.; Demyanko, K. V.; Kirilovskiy, S. V.; Nechepurenko, Yu. M.; Poplavskaya, T. V. Modeling of transonic transitional three-dimensional flows for aerodynamic applications. AIAA J.; 2021; 59, pp. 3598-3610. [DOI: https://dx.doi.org/10.2514/1.J060174] 1501.76039

16 Shen, J. Stable and efficient spectral methods in unbounded domains using Laguerre functions. SIAM J. N-umer. Anal.; 2000; 28, pp. 1113-1133.1786133 [DOI: https://dx.doi.org/10.1137/S0036142999362936] 0979.65105

17 Gautschi, W. Gauss–Radau formulae for Jacobi and Laguerre weight functions. Math. Comput. Simul.; 2000; 54, pp. 403-412.1808697 [DOI: https://dx.doi.org/10.1016/S0378-4754(00)00179-8] 0981.41017

18 Boiko, A. V.; Demyanko, K. V.; Nechepurenko, Yu. M.; Zasko, G. V. On the use of probability-based methods for estimating the aerodynamic boundary-layer thickness. Fluids; 2021; 6, 267. [DOI: https://dx.doi.org/10.3390/fluids6080267] 1536.76043

19 Szegö, G. Orthogonal Polynomials; 1975; Providence, R.I., Am. Math. Soc.: 0305.42011

20 Shen, J.; Tang, T. Spectral and High-Order Methods with Applications; 2006; Beijing, Science: 1234.65005

21 Wang, P.; Huybrechs, D.; Vandewalle, S. Explicit barycentric weights for polynomial interpolation in the roots or extrema of classical orthogonal polynomials. Math. Comput.; 2014; 83, pp. 2893-2914.3246814 [DOI: https://dx.doi.org/10.1090/S0025-5718-2014-02821-4] 1297.41001

22 Nechepurenko, Yu. M. On the dimension reduction of linear differential-algebraic control systems. Dokl. Math.; 2012; 86, pp. 457-459.3027170 [DOI: https://dx.doi.org/10.1134/S1064562412040059] 1264.93033

23 Nechepurenko, Yu. M.; Sadkane, M. A low-rank approximation for computing the matrix exponential norm. SIAM J. Matrix Anal. Appl.; 2011; 32, pp. 349-363.2817492 [DOI: https://dx.doi.org/10.1137/100789774] 1298.65085

24 Zasko, G. V.; Boiko, A. V.; Demyanko, K. V.; Nechepurenko, Yu. M. Simulating the propagation of boundary-layer disturbances by solving boundary-value and initial-value problems. Russ. J. Numer. Anal. Math. Model.; 2024; 39, pp. 47-59.4703259 [DOI: https://dx.doi.org/10.1515/rnam-2024-0005] 1536.76043