On the three-dimensional stability of Poiseuille flow in a finite-length duct

TL;DR

This paper investigates the three-dimensional stability of Poiseuille flow in a finite-length duct, revealing boundary-layer structures in eigenmodes at high Reynolds numbers and questioning the linear model's validity under such conditions.

Contribution

It introduces a SUPG-based scheme for eigenvalue problems and uncovers the boundary-layer structure of least-stable eigenmodes at high Reynolds numbers.

Findings

Flows are asymptotically stable up to Re=2500.

Eigenmodes develop boundary-layer structures at high Re.

Linear models may be inadequate at high Re due to singular eigenmodes.

Abstract

The stability of a three-dimensional, incompressible, viscous flow through a finite-length duct is studied. A divergence-free basis technique is used to formulate the weak form of the problem. A SUPG (streamingline upwind Petrov-Galerkin) based scheme for eigenvalue problems is proposed to stabilize the solution. With proper boundary condtions, the least-stable eigenmodes and decay rates are computed. It is again found that the flows are asymptotically stable for all $Re$ up to $2500$. It is discovered that the least-stable eigenmodes have a boundary-layer-structure at high $Re$, although the Poiseuille base flow does not exhibits such structure. At these Reynolds numbers, the eigenmodes are dominant in the vicinity of the duct wall and are convected downstream. The boundary-layer-structure brings singularity to the modes at high $Re$ with unbounded perturbation gradient. It is shown that due to the singular structure of the least-stable eigenmodes, the linear Navier-Stoker operator tends to have pseudospectrua and the nonlinear mechanism kicks in when the perturbation energy is still small at high $Re$. The decreasing stable region as $Re$ increases is a result of both the decreasing decay rate and the singular structure of the least-stable modes. The result demonstrated that at very high $Re$, linearization of Navier-Stokes equation for duct flow may not be a good model problem with physical disturbances.

Equations

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$inletcondition$$

$$wallconidtion$$

$$outletcondtion$$

$$\mathbf{u}={\mathbf{f}}_0(x,y,z),whent=0,{\mathbf{f}}_{0}satisfiesaboveB.C..$$

$$\frac{{\partial }^2{w}_0}{\partial x^2}+\frac{{\partial }^2{w}_0}{\partial y^2}=\frac{1}{\nu }\frac{\partial p_0}{\partial z}.$$

$$w_0(x,y)=$$

$$\frac{1}{2}\sum _{m=1}^{\infty }sin(m\pi x)(A_m{e}^{m\pi y}+B_m{e}^{-m\pi y})+$$

$$\frac{1}{2}\sum _{m=1}^{\infty }sin(m\pi y)(C_m{e}^{m\pi x}+D_m{e}^{-m\pi x})],$$

$$0\le x\le 1,0\le y\le 1$$

$$A_m$$

$$B_m$$

$$C_m$$

$$D_m$$

$$K_0$$

$$\frac{\partial {\mathbf{u}}_{ϵ}^{\prime }}{\partial t}+{\mathbf{u}}_0\cdot \nabla {\mathbf{u}}_{ϵ}^{\prime }+{\mathbf{u}}_{ϵ}^{\prime }\cdot \nabla {\mathbf{u}}_0$$

$$\nabla \cdot {\mathbf{u}}_{ϵ}^{\prime }$$

$$inletcondition$$

$$wallconidtion$$

$$outletcondtion$$

$${\mathbf{u}}_{ϵ}^{\prime }={{\mathbf{f}}^{\prime }}_{ϵ0}(x,y,z),whent=0,{\mathbf{f}}_{ϵ0}satisfiesaboveB.C..$$

$${\mathbf{u}}_{ϵ}^{\prime }$$

$$p_{ϵ}^{\prime }$$

$$\lambda \mathbf{u}=-{\mathbf{u}}_0\cdot \nabla \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$find$$

$$\lambda M_0(\mathbf{u},\mathbf{w})$$

$$M_0(\mathbf{u},\mathbf{w})$$

$$A_0(\mathbf{u},\mathbf{w})$$

$$C_0(\mathbf{u},\mathbf{w})$$

$$||\mathbf{u}|{|}_{L^2(\Omega )}$$

$${\mathbf{V}}_h^0$$

$$=Span({\mathbf{\Phi }}_1,\dots ,{\mathbf{\Phi }}_N)$$

$${\int }_{{\bar{\Omega }}^e}\nabla \cdot {\mathbf{u}}_{h}d\Omega ={\int }_{\partial {\bar{\Omega }}^e}{\mathbf{u}}_h\cdot \mathbf{n}d\sigma =\sum _{i=1}^6{\int }_{F_i}{\mathbf{u}}_h\cdot \mathbf{n}d\sigma =0,\forall {\bar{\Omega }}^e,$$

$${\hat{\mathbf{\Phi }}}_1(\xi ,\eta ,\zeta )=[0,{\hat{\varphi }}_2+{\hat{\varphi }}_3-{\hat{\varphi }}_6-{\hat{\varphi }}_7,-{\hat{\varphi }}_2+{\hat{\varphi }}_3-{\hat{\varphi }}_6+{\hat{\varphi }}_7{]}^T,$$

$$\{\begin{array}{l}x=\frac{h}{2}\xi +x_0,\\ y=\frac{h}{2}\eta +y_0,\\ z=\frac{h}{2}\zeta +z_0,\end{array}$$

$${\int }_{\partial {\bar{\Omega }}^e}{\mathbf{\Phi }}_{\mathbf{i}}\cdot \mathbf{n}d\sigma =0,\forall {\bar{\Omega }}^e.$$

$${\hat{\mathbf{\Phi }}}_7(\xi ,\eta ,\zeta )=[0,{\hat{\varphi }}_3-{\hat{\varphi }}_7,-{\hat{\varphi }}_3-{\hat{\varphi }}_7{]}^T,$$

$$N_{dof}=$$

$$+2(N_x-2)(N_y-1)+2(N_x-1)(N_y-2).$$

Full text

Department of Mechanical, Aerospace and Nuclear Engineering, Rensselaer Polytechnic Institute

1. Introduction

The stability of pipe or duct Poiseuille flow and its transition to turbulence remains an active topic. Reynolds [1] first demonstrated that the pipe flow usually makes a transistion to turbulence at Reynolds number near $2000$. It was shown by Huerre and Rossi [2] that laminar flow can be achieved up to Reynolds number $~{}10^{5}$ by carefully reducing external disturbances. Various numerical simulations have been studied to analyze pipe flow stability, including [3], [4], [5], [6]. In addition, the linearized pipe Poiseuille flow has been investigated by [7], [8], [9], [10]. Duct flows with various cross-section shapes were studied by Delplace [11] and similar behavior as pipe flow was found. It seems such flows are linearly stable for all Reynolds numbers, yet they tend to be unstable above certain Reynolds number in experiments, see [12], [13].

One way to explain such phenomena is to analyze the nonlinear dynamics and identify the disturbance threshold that can keep the laminar flow. The inverse relationship between the level of flow perturbation and the Reynolds number up to which the laminar flow can sustain has been studied by [14], [15]. A question is therefore raised to identify the scaling factor $\gamma$, where $\epsilon=O(Re^{\gamma})$, with $\epsilon$ the mimial amplitude of all perturbations that can trigger instability. Hof [14] showed by experiments that the factor $\gamma=-1$, while other estimates ranges from $-1$ to $-7/4$, see [16], [17]. Another way to understand this sensitivity is to look at the pseudospectra of the linearized problem, see also [10], [15]. The linear operator resulting from Navier-Stokes equations exhibits strong sensitivity at high $Re$ with respect to pipe geometry fluctuations and thus instability can occur with finite $Re$.

In this paper we computed the stability of three-dimensional Poiseuille flow in a finite-length duct, with physical boundary conditions that fixes the inlet velocity and the normal stresses at the outlet. The three-dimensional structure of the least-stable modes are resolved and the eigenvalues are computed. It is again shown that the duct Poiseuille flow is linearly stable up to $Re=2500$. The study uses a SUPG(Streamline upwind Petrov-Galerkin) based finite element method to compute a generalized matrix eigenvalue problem. The SUPG scheme, proposed by Hughes [18], [19], [20], is well-known to provide the answer of finding a higher-order accurate method for the advection-diffusion equations without deterioration of convergence rate. The numerical scheme also employs a weakly-divergence-free basis technique to overcome the singular mass matrix difficulty brought by the incompressibility condition. The divergence-free basis numercial scheme has been studied by various articles. Griffths [21], [22], [23] proposed several element-based divergence-free basis on triangular and quadrilateral elements. Fortin [24] also investigated the discrete divergence-free subspace with piecewise constant pressure space and various discrete velocity spaces. Ye and Hall [25] showed another discrete divergence-free basis with continuous discrete pressure space. The main advantage of weakly-divergence-free basis technique is that the degrees of freedom is greatly reduced when solving a large scale matrix problem. Apparantly, this method gives added benefits when analyzing the eigenvalue problem with incompressibility constraint.

The computation result shows a boundary-layer-structure of the least-stable eigenmodes at high $Re$. As $Re$ increases, the modes are dominant close to the duct wall, forming a thin boundary layer with large gradients. The interesting structure, however, is not accompanied by a boundary layer of the base flow. The unbounded mode gradient is crucial in understanding the nonlinear effects generated by small perturbations. In high $Re$ flows, even physically small disturbances can induce non-negligible nonlinear effects, due to the large value of $\mathbf{u}\cdot\nabla\mathbf{u}$ where $\mathbf{u}$ denotes the velocity perturbation. The boundary-layer-structure also qualitatively explains the sensitivity of the linearized Navier-Stokes problem in the pseudospectra theory. Due to the eigenmodes dominance in the vicinity of the duct wall, small fluctuations of duct geometry may significantly alter the least-stable mode shapes. This may cause large deviation of eigenvalues from ideal with small linear operator perturbations. This result gives another point of view regarding studies in [10], [15].

2. Mathematical Model

2.1. Unsteady Navier-Stokes equations

We consider a visous, incompressible Newtonian flow through a three-dimensional duct. The duct has a sqaure cross-section. The distances are scaled with the duct width $d$ and the non-dimensionalized duct length is $L$, see figure 1. Rectangular coordinates are used, with $x$-$y$ plane lying on the duct cross-section and $z$ axis extending along the duct length. The flow domain is given by $\Omega=\{(x,y,z)|0\leq x\leq 1,0\leq y\leq 1,0\leq z\leq L\}$. The velocity components are scaled with the characteristic speed $U$, such that the nondimensionalized flux through the duct is unity. The fluid density $\rho$ and the viscosity $\mu$ are constants. The pressure is scaled with the dynamic pressure $\rho U^{2}$. The time is scaled with $\frac{d}{U}$. As a result, the nondimensional viscosity is $\nu=\frac{1}{Re}=\frac{\mu}{\rho UL}$. The nondimensionalized Navier-Stokes equations are described as follows,

[TABLE]

with boundary conditions,

[TABLE]

and initial condition is,

[TABLE]

Here the inlet velocity is fixed to the inlet velocity profile $\mathbf{u}_{0}(x,y)$, which will be described in detail shortly after. The duct wall imposes a no-slip condition. The outlet is subjected to a natural Neumann condition. The vector $\mathbf{n}$ denotes the unit outward normal vector of the duct outlet at $z=L$. This outlet condition decouples the velocity field $\mathbf{u}$ and the pressure field $p$ in the weak form. At high Reynolds number, the outlet boundary condition approaches a constant pressure outlet condition. The initial condition can be chosen arbitrarily.

It is known that there exists a steady state solution $\mathbf{u}_{0}(x,y,z)=[0,0,w_{0}(x,y)]^{T}$, $p(x,y,z)=p_{0}(z)$ to the above model problem, satisfying,

[TABLE]

Here $\frac{1}{\nu}\frac{\partial p_{0}}{\partial z}$ is a negative scaled constant such that the constraint of total flux $1$ is satisfied. It can be verified that the solution of eqation (5) is,

[TABLE]

[TABLE]

Given the steady state solution (6), we consider the associated linear stability problem. We denote the perturbation to the base flow $\mathbf{u}_{0}$ as $\mathbf{u^{\prime}_{\epsilon}}$. The velocity field can be expressed as $\mathbf{u}=\mathbf{u}_{0}+\mathbf{u^{\prime}_{\epsilon}}$. Similarly the pressure field $p$ can be expressed as $p=p_{0}+p^{\prime}_{\epsilon}$. By neglecting higher oder terms when the perturbations are small, the linearized Navier-Stokes equations are,

[TABLE]

Boundary conditions are,

[TABLE]

and initial condition is,

[TABLE]

2.2. Eigenvalue problem formulation

We now turn to study the spectral decomposition of the above linear system. The mode can be expressed as,

[TABLE]

For simplicity, we will omit the prime ′ of the velocity perturbation mode $\mathbf{u^{\prime}}$ and the pressure perturbation mode $p^{\prime}$ from now on. The equations of the corresponding eigenvalue problem after substituting the velocity and pressure fields are,

[TABLE]

Here $\mathbf{u}$ and $p$ are mode shapes and they are functions of space coordinates $(x,y,z)$ only. Equations (12) are subjected to boundary conditions (8), (9), (10).

The sign of eigenvalue $\lambda$ determines whether each mode is asymptotically stable or not. It is crucial to find the least-stable eigenvalues with the greatest real part. These most dangerous modes are dominant either when the flow perturbation is decaying, or when the flow is unstable as the perturbation grows exponentially.

2.3. Divergence-free subspace and weak form

The classical approach in finding the solutions of equations (12) is to solve the velocity mode $\mathbf{u}$ in functional space $\mathbf{V}=\left\{\mathbf{u}\in[H^{1}(\Omega)]^{3}|~{}\mathbf{u}~{}satisfies~{}boundary~{}condition~{}(\ref{E:eigenBCinlet})-(\ref{E:eigenBCoutlet})\right\}$ and the pressure mode in $Q=\left\{p|p\in L^{2}(\Omega)/\mathbb{R}\right\}$. We now try to find the solution of $\mathbf{u}$ in its subspace $\mathbf{V}^{0}=\left\{\mathbf{u}\in\mathbf{V}|\nabla\cdot\mathbf{u}=0\right\}$. When limiting the solution space in the subspace $\mathbf{V}^{0}$, the incompressiblity constraint is no longer needed. Thus a weak form can be expressed from equation (12),

[TABLE]

The inner-product is defined as the integral over $\Omega$, i.e., $\left\langle f,g\right\rangle=\int_{\Omega}fgd\Omega$. Note the pressure term and boundary integrals are eliminated as a result of the specified boundary conditions  (8)-(10). The eigenfunctions are normalized such that its $L_{2}$ norm is unity.

2.4. Finite element formulations

A $Q_{1}-P_{0}$ mixed form is used in this paper. It is well-known that the trilinear velocity-constant pressure element suffers from a checkerboard pressure mode on regular meshes. The spurious pressure mode is due to the fact that this element does not satisfy the Babǔska-Brezzi(B.B.) condition. However, if the velocity field is the only interested variable, the spurious pressure mode is irrelevant, since pressure term is eliminated in equation (12). Boffi et al. ([26]) have shown that such mixed form is valid in analyzing eigenvalue problem with the Galerkin method.

The finite element formulation first constructs a weakly-divergence-free basis with $Q_{1}$ elements. Then a discretized weak form from equation (2.3) is formulated.

Let $\mathbf{V}_{h}$ and $Q_{h}$ be the finite dimensional subspace of $\mathbf{V}$ and $Q$. The weakly divergence-free subspace follows as,

[TABLE]

Here the divergence-free basis function $\mathbf{\Phi}_{i}$ belongs to the $Q_{1}$ trilinear element. The trial function $q_{h}$ is constant in each element. The discrete subspace is equivalent as a subspace imposing the conservation of mass condition on each element $\bar{\Omega}^{e}$, where $\Omega^{e}$ denotes the interior of the eth element.

[TABLE]

where $F_{i}$ are the six faces of a cubic element. Fortin [24] showed that in three-dimensional cubic regular element, the weakly-divergence-free basis $\mathbf{\Phi}_{i}$ can be expressed as vortices lying on each face of the element, shown in figure 2. The six faces of a cubic element is associated with six localized vortices. The basis can be physically viewed as a basis for a vorticity field $\mathbf{\omega}=\nabla\times\mathbf{u}$, in the discretized sense. The direction of the ’effective vorticity’ is normal to each element face. Note that for each element $\Omega^{e}$, the six basis functions are linearly dependent and the degree of freedom is five.

One of the weakly-divergence-free basis can be expressed in the reference domain as,

[TABLE]

where $\hat{\phi}_{i}(\xi,\eta,\zeta)$ is the typical trilinear $Q_{1}$ shape function, with nodal value $1$ on each of the element vertex. The reference element domain is shown in figure 3, the domain $\hat{T}=\left\{(\xi,\eta,\zeta)|~{}\xi\in[-1,1],\eta\in[-1,1],\zeta\in[-1,1]\right\}$. Shape function index is associated with node index described in figure 3. The other weakly-divergence-free basis $\mathbf{\hat{\Phi}}_{i}(\xi,\eta,\zeta)$ on the reference domain can be constructed similarly.

The finite element uses a simplest uniform cubic element discretizing domain $\Omega$ throughout this paper. The cube has identical edge length $h$ in each direction. As a result, the coordinate transformation between $(x,y,z)$ and $(\xi,\eta,\zeta)$ has a simple linear relationship,

[TABLE]

where $x_{0},y_{0},z_{0}$ are coordinates of the element center. With the above assumptions, it can be verified that the weakly-divergence-free condition is satified for each $\mathbf{\Phi_{i}}(x,y,z)=\mathbf{\hat{\Phi}}_{i}(\xi,\eta,\zeta)$,

[TABLE]

It should be pointed out that due to the non-zero outlet boundary condition (10), ’half vortex’ elements must be included in addition to the ’full vortex’ element described above, see figure 4. One of the ’half vortex’ element on the reference domain is,

[TABLE]

The index $7$ is labelled after the $6$ ’full vortices’ on each of the element faces. The other basis functions at the duct outlet can be expressed in a similar way. The transformation to physical domain is straightforward. It can be verified that the ’half vortex’ basis satisfies the weakly-divergence-free condition. There is also a linear dependency of four ’half vortices’ $\mathbf{\hat{\Phi}}_{7}$ to $\mathbf{\hat{\Phi}}_{10}$ and a ’full vortex’ $\mathbf{\hat{\Phi}}_{6}$ on the element outlet surface. Thus global basis $\mathbf{\Phi}_{i}$ mapped from $\mathbf{\hat{\Phi}}_{6}$ must be excluded from the final basis set.

With $N_{x}N_{y}N_{z}$ elements in the flow domain, where $N_{x}$, $N_{y}$, $N_{z}$ are the number of elements in $x$, $y$, $z$ directions, the total degree of freedom in the finite element formulation is,

[TABLE]

This is a result of the linear dependency of the six basis functions on each element and the inclusion of boundary vortex elements. It should be noted that all vortices with ’effective vorticity’ along the $z$ direction are linearly dependent with the rest of the vortices, thus they must be excluded from the divergence-free basis set. The basis functions $\mathbf{\Phi}_{i}$ are constructed such that the ’effective vorticity’ has positive component in each direction. This ensures the consistency when $\mathbf{\Phi}_{i}$ is evaluated on nearby elements.

Although the construction of the finite element scheme enforces stric cubic elements, it is perhaps the simplest way to analyze the problem with minimal degrees of freedom.

2.5. Streamline upwind Petrov-Galerkin(SUPG) method

Hughes [20] introduced SUPG scheme for computing time-dependent and steady state fluid problems, and enjoyed great success. When element Peclet number, defined as $\alpha=h|\mathbf{u}_{0}|/(2\nu)$, where $|\mathbf{u}_{0}|$ is the maximum magnitude of $\mathbf{u}_{0}$ in each element, is large, the Galerkin method suffers from wild oscillations that contaminate the numerical soltuions. In eigenvalue analysis, the same numerical instability prevales when the element Peclet number is large. Thus a SUPG formulation for the eigenvalue problem is proposed to eliminate the instability.

Similiar to the classical SUPG method, we take the inner-product of both sides of equation (12) with $\mathbf{w}_{h}$ and with $\tau\mathbf{u}_{0}\cdot\nabla\mathbf{w}_{h}$ in element internal domain $\tilde{\Omega}=\cup\Omega^{e}$, where $\Omega^{e}$ is the internal of the eth element. It follows that the discretized weak form is,

[TABLE]

The parameter $\tau$ can be chosen as $\tau=min(\alpha/3,1)h/(2|\mathbf{u}_{0}|)$, see [27].

Appendix gives the convergence analysis of this new scheme. Convergence rate is shown to be $O(h^{2})$ for eigenvalue and $O(h)$ for eigenmode energy norm.

The weak form (2.5) is equivalent to a generalized eigenvalue problem,

[TABLE]

The eigenmode $\mathbf{u}$ to be solved is $\mathbf{u}=\sum_{i}c_{i}\mathbf{\Phi}_{i}$.

The implicitly-restarted Arnoldi method, provided by ARPACK, or MATLB eigs function is capable of finding the least-stable eigenvalues and eigenmodes.

3. Numerical Results

We now show some of the computed results of the new numerical scheme.

$Re$ from $500-2500$ are computed for a finite-length duct with nondimensional duct length $L=3$. Regular mesh size of $N_{x}=N_{y}=Nz/L$ is adopted.

Figure 5 shows the mesh convergence plot. Computed least-stable eigenvlaue $\lambda_{1},\lambda_{2}$ and the third eigenvalue $\lambda_{3}$ tend to converge as the degrees of freedom increases. The following analyasis will all adopt the case with $N_{x}=40$, which is the largest degree of freedom computed so far.

The computed eigenvalues show that the largest eigenvalue has multiplicity of $2$. Acutually the first two eigenmodes are identical by interchanging $x,y$ coordinates and $u,v$ components. The third eigenmode has a multiplicity of $1$ and this eigenmode possesses $x$-$y$ symmetry.

Figure 6 shows the three least-stable eigenvalues as a function of $Re$. It is again found that the duct Poiseuille flow is linearly stable for all $Re$ up to $2500$, which agrees with previous pipe or duct theoretical studies.

It is interesting to investigate the leading-order eigenmode structures. Figure 7 shows the normalized eigenmode corresponding to the first eigenvalue $\lambda_{1}=-0.4448$ at $Re=2000$. Figure 8 shows the normalized eigenmode corresponding to the third eigenvalue $\lambda_{3}=-0.4698$ at the same $Re$. Mesh size with $N_{x}\times N_{y}\times N_{z}=40\times 40\times 120$ are used to compute the solutions. The cut planes are chosen at places where the mode is dominant. The second mode is not plotted because it is a symmetric mode of the first one. Although complicated in the structures of the shape, there is a characteristic boundary-layer-structure that prevails at all high $Re$ leading-order modes shapes. The modes are dominant in the vicinity of the wall and are convected downstream. When the pipe inlet profile is fixed to the base flow profile $\mathbf{u}_{0}$, least-stable perturbations tends to dominate near the outlet. It can also be observed that these modes do not have a simple normal mode structure, e.g.,

[TABLE]

where $k$ is the wave number and $\Psi$ is a common shape function. Therefore this simple theoretical assumption may not be the best way to analyze pipe or duct flow stability.

The boundary layer gets thinner and closer to the duct wall as $Re$ increases. The gradient of the normalized modes become unbounded as $Re$ grows. Figure 9 shows the side view of the first eigenmode $w(x,y,z=2.9)$ corresponding to $\lambda_{1}$ at $Re=500,1000,2500$. The position where the mode function peaks are labelled. It can be seen the peak position from the duct wall decreases from $x\sim 0.1375$ at $Re=500$ to $x\sim 0.1$ at $Re=2500$. The physical consequences are discussed in the next section. It seems true that ’instability’ tends to be generated near the wall and the outlet in the very beginning phase of transition.

4. Relationship with pseudospectra theory and nonlinear Effects

The relationship with pseduospectra theory and nonlinear effects of Poiseuille flow is studied.

It has been investigated in [15] that the instability of pipe Poiseuille is attributed to the sensitivity of the linear operator $\mathcal{L}$ from the linearization of the Navier-Stokes equations, where $\mathcal{L}\mathbf{u}=-\mathbf{u_{0}}\cdot\nabla\mathbf{u}-\mathbf{u}\cdot\nabla\mathbf{u_{0}}-\nabla p+\nu\nabla^{2}\mathbf{u}$. It was shown that when $Re$ increase, a small perturbation to the linear operator $\mathcal{L}$ results in a relatively large eigenvalue deviation. As a result, pipe Poiseuille flow may become unstable at finite $Re$ although in ideal case it is linearly stable. We connect this claim to the boundary-layer-structure of the least-stable eigenmode, in a qualitative way. The boundary-layer-structure has most of its active part near the duct wall. As $Re$ becomes higher, the layer becomes thinner and closer to the wall. It is expected that a small fluctuation of the duct wall configuration will alter the eigenmode shape in a significant way. The base flow field, however, is not affected greatly by wall fluctuations since its velocity field is dominant at the bulk. With this physical picture in mind, a perturbation to $\mathcal{L}$ induced by small fluctuation of duct wall will change the eigenvalue noticeablely as a result of the singularity of mode shapes in the limit of high $Re$. This gives a physical point of view regarding the pseudospectra theory.

To investigate the nonlinear effects, we assume the initial condition to be one of the least-stable modes $\epsilon\mathbf{u}_{1}$, where $||\mathbf{u}_{1}||=1$ corresponding to $\lambda_{1}$. The instananeous flow energy growth is investigated. Physically it can be regarded as the study of the nonlinear dynamics for a short period of time. Multiplying $\mathbf{u}$ into equation (1) and integrate on both sides gives,

[TABLE]

where $A_{0}$, $C_{0}$ are defined in equation (2.3). Since the initial condition is chosen to be the least-stable eigenmode $\epsilon\mathbf{u}_{1}$, it follows,

[TABLE]

where $\lambda_{1}$ is the least-stable eigenvalue. Thus the total flow energy $\frac{1}{2}\left<\mathbf{u},\mathbf{u}\right>$ growth rate at instant $t=0$ is driven by the stabilizing linear effect and the nonlinear effect. By switching the sign of the perturbation amplitude $\epsilon$, the nonlinear energy production rate can always be positive, we therefore consider the nonlinear term to be destabilizing in the worst case.

If the nonlinear effect is negligible, the flow energy decays. The assumption to take the eigenmode as the initial condition physically makes sense when flow is attracted to the Poiseuille base flow branch. On the other hand, when $\epsilon$ is large, the nonlinear energy production rate may exceed the linear stabilizing effect and the flow energy starts to grow at $t=0$. The linearization assumption breaks down in this situation. Thus we study the critical amplitude $\epsilon$ at various $Re$. Perturbation with above the critical value may not result in a final transistion into turbulence since nonlinear dynamics after $t=0$ will be complicated and is not analyzed. However, it is a sign that linear stability problem fails to give a proper mathematical model and the flow tends to be dominated by nonlinear effects. With the difficulty in quantifying subsequent nonlinear dynamics, the analysis only gives a rough estimate on how small perturbations may generate non-negligible destablizing effects. It should be noted that a eigenmode initial perturbation is different from the other studies to induce impulsive or periodic distrubances and this is not very likely to be physically induced.

The critical amplitude is determined by a differential balance equation,

[TABLE]

When the energy production term $|u_{1}^{k}\nabla_{k}u_{1}^{i}u_{1}^{i}d\Omega|$ exceeds the linear stablizing term $\lambda_{1}u_{1}^{i}u_{1}^{i}d\Omega$ in any infinitesimal volume, the linearzation assumption is no longer valid. Thus the critical amplitude of the flow energy can be computed numerically by finite element integration on each hexahedral element. Figure 10 shows the critical perturbation amplitude $\epsilon$ vs $Re$, the effective nonlinear energy production $\eta=\sup\frac{|u_{1}^{k}\nabla_{k}u_{1}^{i}u_{1}^{i}d\Omega|}{|\lambda_{1}u_{1}^{i}u_{1}^{i}d\Omega|}$ vs $Re$, and the first eigenvalue $\lambda_{1}$ vs $Re$. The threshold is seen to be decreasing greatly as $Re$ increases. Any reasonable perturbation at very high $Re$ tend to have non-negligible nonlinear effect. A fit to the $\epsilon-Re$ curve gives approximately $\epsilon\sim O(Re^{-0.6})$. This number is different from the results reported by [14], [16] and [17], mainly because the initial condition given here is not applicable to the other studies. It may also attribute to the difference between pipe and duct flow. The significance of the nonlinear analysis demonstrates that the linear stability assumption fails to be a good model problem as $Re$ becomes higher. Two reasons contribute to the breakdown. First, the decreasing absolute value of first eigenvalue as $Re$ increass. Second, the increasing nonlinear energy production rate as a result of the growing eigenfunction gradient near the duct wall, as it is found that the maximum nonlinear energy production is near the peak of the boundary layer.

5. Conclusions

In this paper have presented a SUPG based finite element method with divergence-free-basis technique to compute the Poiseuille flow eigenvalue problem. The flow is found to be linearly stable for $Re$ up to $2500$ and the modes show a boundary-layer-structure at high $Re$. This characteristic structure at high $Re$ accounts for the sensitivity of the linearized Navier-Stokes problem with respect to geometry variation as well as the sensitivity of nonlinear effects.

6. Appendix

The convergence analysis in this paper assumes that the domain boundary $\partial{\Omega}$ is smooth enough, and also that the base flow $\mathbf{u}_{0}$, the eigenmode $\mathbf{u}$ and the corrresponding pressure field $p$ are smooth enough, for simplicity.

Space $\mathbf{V}=\left\{\mathbf{u}\in[H^{1}(\Omega)]^{3}|~{}\mathbf{u}~{}satisfies~{}boundary~{}condition~{}(\ref{E:eigenBCinlet})-(\ref{E:eigenBCoutlet})\right\}$ and $Q=\left\{p|p\in L^{2}(\Omega)/\mathbb{R}\right\}$ are now assumed to be complex spaces. $\mathbf{V}_{h}\subset\mathbf{V}$, $Q_{h}\subset Q$, $\mathbf{V}^{0}=\left\{\mathbf{u}\in\mathbf{V}|\nabla\cdot\mathbf{u}=0\right\}$ and $\mathbf{V}_{h}^{0}=\left\{\mathbf{u}_{h}\in\mathbf{V}_{h}|\left\langle\nabla\cdot\mathbf{u}_{h},q_{h}\right\rangle=0,~{}\forall q_{h}\in Q_{h}\right\}$ also extend to complex spaces accordingly.

Lemma 6.1**.**

Let

[TABLE]

where $\mathbf{u},\mathbf{w}\in\mathbf{V}$. $A_{0}$, $C_{0}$ are defined in equation (2.3). There exists a sufficiently large $\gamma>0$ such that,

[TABLE]

Proof.

First we have,

[TABLE]

Next $\left<\mathbf{u}_{0}\cdot\nabla\mathbf{u},\mathbf{u}\right>$ is estimated,

[TABLE]

Here $\nabla_{k}u_{0}^{k}=0$ is implied. From equation (27),

[TABLE]

Here $\mathbf{n}$ is the unit outward normal vector. Since the problem assumes outflow at the outlet, namely $u_{0}^{k}n^{k}|_{outlet}>0$, and at the inlet and wall $u^{i}=0$, it follows that,

[TABLE]

Then $\left<\mathbf{u}\cdot\nabla\mathbf{u}_{0},\mathbf{u}\right>$ is estimated. Since the base flow field $\mathbf{u}_{0}$ and its gradient is bounded, there exists a constant tensor $N_{ki}$ such that,

[TABLE]

Then,

[TABLE]

Here $\gamma$ is a sufficiently large positive number. It is apparant that when $\gamma>0$ is large enough, matrix $\gamma I-N$ is positive definite.

It then follows,

[TABLE]

Combining equations (26), (27), (32), $\exists~{}\alpha_{0}>0$, $a(\mathbf{u},\mathbf{u})+c(\mathbf{u},\mathbf{u})\geq\alpha_{0}||\mathbf{u}||_{1}^{2}$.

∎

Lemma 6.2**.**

Let $\Omega^{e}$ be the interior of the eth element, $\tilde{\Omega}=\cup\Omega^{e}$, and,

[TABLE]

where $\mathbf{u}_{h},\mathbf{w}_{h}\in\mathbf{V}_{h}$. $\mathbf{V}_{h}$ is the discretized continuous piecewise polynomial solution space that satifies boundary conditions (8)-(10).

[TABLE]

Here $\alpha$ is the element Peclet number. $\zeta(\alpha)$ satisfies

[TABLE]

where $C$ is the inverse estimate constant that satifies,

[TABLE]

$h$* is the mesh parameter. In the regular cubic element case, it is the edge length. $|\mathbf{u}_{0}|$ denotes the maximum of base flow velocity magnitude within each element. $\tau$ varies across each element. The subscript $\tilde{\Omega}$ means that the inner-product is taken in the interior of each element and then summing the elementwise integration up.*

Then $\exists~{}\gamma>0$ defined in lemma 6.1, and $\exists~{}\delta>0$, when $\tau<\delta$, the following coercive condition holds,

[TABLE]

Proof.

[TABLE]

First we recall the inverse estimate relationship in equation (34),

[TABLE]

With this inequality, $-\left<\nu\nabla^{2}\mathbf{u}_{h},\tau\mathbf{u}_{0}\cdot\nabla\mathbf{u}_{h}\right>_{\tilde{\Omega}}$ is estimated,

[TABLE]

Then $\left<\mathbf{u}_{h}\cdot\nabla\mathbf{u}_{0},\tau\mathbf{u}_{0}\cdot\nabla\mathbf{u}_{h}\right>_{\tilde{\Omega}}$ is estimated,

[TABLE]

Here $N$ is a constant due to the boundedness of $\nabla\mathbf{u}_{0}$.

Next we estimate $\left<\gamma\mathbf{u},\tau\mathbf{u}_{0}\cdot\nabla\mathbf{u}_{h}\right>_{\tilde{\Omega}}$,

[TABLE]

Subsituting all the above estimations into equation (35) and use the results in lemma 6.1,

[TABLE]

Here $\alpha_{1}$ is the estimation constant in lemma 6.1. It follows that when

[TABLE]

∎

Remark 6.0.1**.**

For typical flow problem $\alpha_{1},\gamma$ can be chosen on the order of $N\sim|\nabla\mathbf{u}_{0}|$, so the condition for $\tau$ is not restrictive.

Remark 6.0.2**.**

The introduction of ’$\gamma$’ terms ensures the coercivity of the linear problem. The corresponding eigenvalue problem is,

[TABLE]

Here $M(\mathbf{u},\mathbf{w})$ is a bilinear form defined in equation

The discretized problem is equivalent to a matrix eigenvalue problem,

[TABLE]

Here $H$, $M$ are matrices defined in equation (2.5). It is apparant that the eigenvalue is shifted by $\gamma$ from the original problem 21, and that the eigenfunctions remain unchanged.

Lemma 6.3**.**

Let

[TABLE]

Define a linear system in its weak form, $\mathbf{f}\in\mathbf{V}$,

[TABLE]

Here the binlinear form $b(\mathbf{u},q)$ satisfies the B.B. condition such that,

[TABLE]

The corresponding discretized weak form is,

[TABLE]

Then whether or not the B.B. condition is satisfied or not in equations (46), there is a unique $\mathbf{u}_{h}$ that satifies the discretized equations and,

[TABLE]

Here $C$ is a constant that depends on $\mathbf{u}$ and $p$.

Proof.

It is well-known that equation (45) has a unique solution $(\mathbf{u},p)\in\mathbf{V}\times Q$, due to the coercivity of $a_{0}(\cdot,\cdot)$ and the satisfaction of the B.B condition in $b(\cdot,\cdot)$.

The solution $\mathbf{u}_{h}$ of equation (46) is in subspace $\mathbf{V}_{h}^{0}=Ker(B_{h})$, where

[TABLE]

Thus by restricting the trial function $\mathbf{w}_{h}\in\mathbf{V}_{h}^{0}$,

[TABLE]

Then there exists a unique $\mathbf{u}_{h}\in\mathbf{V}_{h}$ that satisfies equation (46). Note $p_{h}$ may not have unique solution.

Next we slightly tweak the discretized problem such that space $Q_{h}$ is replaced by $\hat{Q}_{h}=Q_{h}/Ker(B_{h}^{T})$, ,

[TABLE]

It is straight-forward to check that the B.B. condition is satisfied in equation (50). Therefore, there exists a unique pair $(\mathbf{u}_{h},\hat{p}_{h})\in\mathbf{V}_{h}\times\hat{Q}_{h}$ as the solution of equation (50). Apparantly $\mathbf{u}_{h}$ is also the unique solution of equation (46).

It is evident that the following relationship also holds,

[TABLE]

Subtracting equation (50) from equation (52), we have,

[TABLE]

Here $\mathbf{u}^{I}\in\mathbf{V}_{h}$, $p^{I}\in Q_{h}$ are any interpolations of $\mathbf{u}$ and $p$. Note $p^{I}$ is not in the truncated subspace $\hat{Q}_{h}$. Becuase the B.B. condition holds in equation (52),

[TABLE]

Here $\hat{B}_{h}$ is defined as,

[TABLE]

If the B.B condition does not hold in equation (46), $Ker(B_{h}^{T})\neq\Phi$, therefore any $p^{I}\in Q_{h}$ can be expressed as,

[TABLE]

We now demonstrate that $Ker(\hat{B}_{h})=Ker(B_{h})$.

First, if ${\mathbf{u}}_{h}\in Ker(B_{h})$,

[TABLE]

Because $\hat{Q}_{h}\subseteq Q_{h}$, it follows $Ker(B_{h})\subseteq Ker(\hat{B}_{h})$.

Next, if ${\mathbf{u}}_{h}\in Ker(\hat{B}_{h})$,

[TABLE]

Pick any $q_{h}\in Q_{h}$, $q_{h}=\hat{q}_{h}+\tilde{q}_{h}$, $\hat{q}_{h}\in\hat{Q}_{h},\tilde{q}_{h}\in Ker(B_{h}^{T}),$

[TABLE]

Then we have $Ker(\hat{B}_{h})\subseteq Ker({B}_{h})$.

Therefore $Ker(\hat{B}_{h})=Ker({B}_{h})$.

From the above discussion we have,

[TABLE]

Let $\mathbf{w}_{h}=\hat{\mathbf{u}}$ and notice $b(\hat{\mathbf{u}},p_{h}-p^{I})=0$, since $\hat{\mathbf{u}}\in Ker(\hat{B}_{h})=Ker(B_{h})$. It follows,

[TABLE]

It is known that B.B. condition holds in equation (50), thus we have

[TABLE]

Therefore, from equation (6),

[TABLE]

Then,

[TABLE]

[TABLE]

Combining equation (61),

[TABLE]

Finally,

[TABLE]

Here constant $C$ may be of different value in each equation for simplicity. ∎

Corollary 6.1**.**

For any $\mathbf{u}\in\mathbf{V}^{0}$, where $\mathbf{V}^{0}$ is the divergence-free subspace of $\mathbf{V}$, there exists an interpolation $\tilde{\mathbf{u}}^{I}\in\mathbf{V}_{h}^{0}$ such that,

[TABLE]

This corollary shows that the divergence-free vector field $\mathbf{u}$ can be approximated by a weakly-divergence-free interpolation with the above error bounds. Whether or not the B.B condition is satisfied or not does not matter.

Theorem 6.1**.**

Given $\mathbf{f}\in\mathbf{V}$, there exists a unique pair $(\mathbf{u},p)\in\mathbf{V}\times Q$ such that,

[TABLE]

Here the binlinear form $b(\mathbf{u},q)$ satisfies the B.B. condition such that,

[TABLE]

Then the following discretized linear system also has a unique solution $\mathbf{u}_{h}\in\mathbf{V}_{h}^{0}$,

[TABLE]

Here the discretized binlinear form $b(\mathbf{u}_{h},q_{h})$ may not satisfy the B.B condition.

If we assume $\mathbf{V}_{h}$ to be continous piecewise polynomials of degree $k$ and $Q_{h}$ be of degree $k-1$, the following esitmation holds,

[TABLE]

Here $C_{u,p}$ is a function of $\mathbf{u}$ and $p$.

Proof.

It is well-known that when B.B. condition is satisfied, equation (66) has a unqiue solution because the coercivity holds true.

The unique solution $(\mathbf{u},p)\in\mathbf{V}\times Q$ also satisfies,

[TABLE]

Combining equation (67) and (69),

[TABLE]

Recall $\mathbf{V}_{h}^{0}$ is the discrete divergence-free subspace. Let $\mathbf{u}^{I}\in\mathbf{V}_{h}^{0}$, be any interpolation of $\mathbf{u}$, and $p^{I}\in Q_{h}$ be any interpolation of $p$. Note $b(\mathbf{w}_{h},p^{I})=0$. We have,

[TABLE]

Let $\mathbf{w}_{h}=\mathbf{u}_{h}-\mathbf{u}^{I}\in\mathbf{V}_{h}^{0}$, from lemma 6.2,

[TABLE]

Therefore,

[TABLE]

Since $\mathbf{u}\in\mathbf{V}^{0}$ and $\mathbf{u}^{I}\in\mathbf{V}^{0}_{h}$, from Corollary 6.1,

[TABLE]

From the choice of $\tau=\frac{h}{2\mathbf{u}_{0}}min(\frac{\alpha}{3},1)$,

[TABLE]

Remember $\alpha=\frac{|\mathbf{u}_{0}|h}{2\nu}$ in each element.

We therefore have,

[TABLE]

where $C$’s are constant numbers in each estimation and $C_{u}$ is a function of $\mathbf{u}$.

We also have,

[TABLE]

Combining all results above and notice the leading order error is $O(h^{2})$ when $k\geq 1$,

[TABLE]

Note all the $C_{u}$, $C_{p}$ are functions of $\mathbf{u}$ and $p$. Since given $\mathbf{f}$, the pair $(\mathbf{u},p)$ is unique. We can also write that,

[TABLE]

The leading-order error comes from $\left<\tau\mathbf{u}_{0}\cdot\nabla(\mathbf{u}_{h}-\mathbf{u}^{I}),\nabla p\right>_{\tilde{\Omega}}$. However, this scheme does not degrade the accuracy of $Q_{1}-P_{0}$ element. It will be interesting to find a way to increase the order of convergence in the future. ∎

Corollary 6.2**.**

Given the linear systems,

[TABLE]

and,

[TABLE]

The operator $\mathbf{T}:\mathbf{V}\rightarrow\mathbf{V}$, $\mathbf{T}_{h}:\mathbf{V}\rightarrow\mathbf{V}_{h}$ is well-defined, with the norm bound,

[TABLE]

Here $T|_{\mathbf{E}(\mu)}$ means the restriction to the generalized eigenspace corresponding to eigenvalue $\mu$ of operator $T$. Here $\mu=\tilde{\lambda}^{-1}$, where $\tilde{\lambda}$ is defined in equation (42), and,

[TABLE]

The convergence of the flow eigenvalue problem is a direct result of Babǔska and Osborn [28],

Theorem 6.2**.**

(Babǔska and Osborn) The distance of eigenspace $\mathbf{E}(\mu)$ and the discretized eigenspace $\mathbf{E}_{h}(\mu)$ satisfies the following bound,

[TABLE]

Here

[TABLE]

Theorem 6.3**.**

(Babǔska and Osborn) Let $\mu$ be a non-zero eigenvalue of $T$ with algebraic multiplicity equal to $m$ and let $\hat{\mu}_{h}$ denote the arithmetic mean of the $m$ discrete eigenvalues of $T_{h}$ converging towards $\mu$. Let $\phi_{1},\cdots,\phi_{m}$ be a basis of generalized eigenvectors in $E(\mu)$ and let $\phi^{*}_{1},\cdots,\phi^{*}_{m}$ be a dual basis of generalized eigenvectors in $E^{*}(\mu)$. Then,

[TABLE]

Theorem 6.4**.**

The eigenvalue problem has an error estimate for the $i$th eigenvalue,

[TABLE]

where $m$ is the algebraic multiplicity of $\tilde{\lambda}_{i}$ and $\tilde{\lambda}_{i_{k}h}$ are the $m$ discrete eigenvalues converging towards $\tilde{\lambda}$.

The estimate for the $i$th eigenmode is,

[TABLE]

Proof.

The proof follows from Boffi [29] Theorem $10.4$ and the above Babǔska and Osborn theory. ∎