On the three-dimensional stability of Poiseuille flow in a finite-length duct
OON THE THREE-DIMENSIONAL STABILITY OF POISEUILLEFLOW IN A FINITE-LENGTH DUCT
LEI XU AND ZVI RUSAK Abstract.
The stability of a three-dimensional, incompressible, viscous flowthrough a finite-length duct is studied. A divergence-free basis technique isused to formulate the weak form of the problem. A SUPG(streaminglineupwind Petrov-Galerkin) based scheme for eigenvalue problems is proposedto stabilize the solution. With proper boundary condtions, the least-stableeigenmodes and decay rates are computed. It is again found that the flowsare 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 thePoiseuille base flow does not exhibits such structure. At these Reynolds num-bers, the eigenmodes are dominant in the vicinity of the duct wall and areconvected downstream. The boundary-layer-structure brings singularity tothe modes at high Re with unbounded perturbation gradient. It is shown thatdue to the singular structure of the least-stable eigenmodes, the linear Navier-Stoker operator tends to have pseudospectrua and the nonlinear mechanismkicks in when the perturbation energy is still small at high Re . The decreasingstable region as Re increases is a result of both the decreasing decay rate andthe singular structure of the least-stable modes. The result demonstrated thatat very high Re , linearization of Navier-Stokes equation for duct flow may notbe a good model problem with physical disturbances. Introduction
The stability of pipe or duct Poiseuille flow and its transition to turbulenceremains an active topic. Reynolds [1] first demonstrated that the pipe flow usuallymakes a transistion to turbulence at Reynolds number near 2000. It was shownby Huerre and Rossi [2] that laminar flow can be achieved up to Reynolds number10 by carefully reducing external disturbances. Various numerical simulations havebeen 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]. Ductflows with various cross-section shapes were studied by Delplace [11] and similarbehavior as pipe flow was found. It seems such flows are linearly stable for allReynolds numbers, yet they tend to be unstable above certain Reynolds number inexperiments, see [12], [13].One way to explain such phenomena is to analyze the nonlinear dynamics andidentify the disturbance threshold that can keep the laminar flow. The inverserelationship between the level of flow perturbation and the Reynolds number upto which the laminar flow can sustain has been studied by [14], [15]. A questionis therefore raised to identify the scaling factor γ , where (cid:15) = O ( Re γ ), with (cid:15) themimial amplitude of all perturbations that can trigger instability. Hof [14] showed Department of Mechanical, Aerospace and Nuclear Engineering, Rensselaer PolytechnicInstitute. a r X i v : . [ phy s i c s . f l u - dyn ] A ug * by experiments that the factor γ = −
1, while other estimates ranges from − − /
4, see [16], [17]. Another way to understand this sensitivity is to look at thepseudospectra of the linearized problem, see also [10], [15]. The linear operatorresulting from Navier-Stokes equations exhibits strong sensitivity at high Re withrespect to pipe geometry fluctuations and thus instability can occur with finite Re .In this paper we computed the stability of three-dimensional Poiseuille flow ina finite-length duct, with physical boundary conditions that fixes the inlet veloc-ity and the normal stresses at the outlet. The three-dimensional structure of theleast-stable modes are resolved and the eigenvalues are computed. It is again shownthat the duct Poiseuille flow is linearly stable up to Re = 2500. The study usesa SUPG(Streamline upwind Petrov-Galerkin) based finite element method to com-pute a generalized matrix eigenvalue problem. The SUPG scheme, proposed byHughes [18], [19], [20], is well-known to provide the answer of finding a higher-order accurate method for the advection-diffusion equations without deteriorationof convergence rate. The numerical scheme also employs a weakly-divergence-freebasis technique to overcome the singular mass matrix difficulty brought by theincompressibility condition. The divergence-free basis numercial scheme has beenstudied by various articles. Griffths [21], [22], [23] proposed several element-baseddivergence-free basis on triangular and quadrilateral elements. Fortin [24] also in-vestigated the discrete divergence-free subspace with piecewise constant pressurespace and various discrete velocity spaces. Ye and Hall [25] showed another discretedivergence-free basis with continuous discrete pressure space. The main advantageof weakly-divergence-free basis technique is that the degrees of freedom is greatlyreduced when solving a large scale matrix problem. Apparantly, this method givesadded benefits when analyzing the eigenvalue problem with incompressibility con-straint.The computation result shows a boundary-layer-structure of the least-stableeigenmodes at high Re . As Re increases, the modes are dominant close to theduct wall, forming a thin boundary layer with large gradients. The interestingstructure, however, is not accompanied by a boundary layer of the base flow. Theunbounded mode gradient is crucial in understanding the nonlinear effects gener-ated by small perturbations. In high Re flows, even physically small disturbancescan induce non-negligible nonlinear effects, due to the large value of u · ∇ u where u denotes the velocity perturbation. The boundary-layer-structure also qualitativelyexplains the sensitivity of the linearized Navier-Stokes problem in the pseudospec-tra theory. Due to the eigenmodes dominance in the vicinity of the duct wall, smallfluctuations of duct geometry may significantly alter the least-stable mode shapes.This may cause large deviation of eigenvalues from ideal with small linear operatorperturbations. This result gives another point of view regarding studies in [10], [15].2. Mathematical Model
Unsteady Navier-Stokes equations.
We consider a visous, incompressibleNewtonian 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-dimensionalizedduct length is L , see figure 1. Rectangular coordinates are used, with x - y planelying on the duct cross-section and z axis extending along the duct length. The flowdomain is given by Ω = { ( x, y, z ) | ≤ x ≤ , ≤ y ≤ , ≤ z ≤ L } . The velocity N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT3
Figure 1.
Duct geometrycomponents are scaled with the characteristic speed U , such that the nondimen-sionalized flux through the duct is unity. The fluid density ρ and the viscosity µ are constants. The pressure is scaled with the dynamic pressure ρU . The time isscaled with dU . As a result, the nondimensional viscosity is ν = Re = µρUL . Thenondimensionalized Navier-Stokes equations are described as follows, ∂ u ∂t + u · ∇ u = −∇ p + ν ∇ u , in Ω(1) ∇ · u = 0 , in Ω , with boundary conditions, inlet condition : u = u ( x, y ) , z = 0 , f or all t ≥ wall conidtion : u = , x = 0 , x = 1 , y = 0 , y = 1 , f or all t ≥ outlet condtion : − p n + ν ∂ u ∂ n = 0 , z = L, f or all t ≥ , (4)and initial condition is, u = f ( x, y, z ) , when t = 0 , f satisf ies above B.C.. Here the inlet velocity is fixed to the inlet velocity profile u ( x, y ), which will bedescribed in detail shortly after. The duct wall imposes a no-slip condition. Theoutlet is subjected to a natural Neumann condition. The vector n denotes theunit outward normal vector of the duct outlet at z = L . This outlet conditiondecouples the velocity field u and the pressure field p in the weak form. At highReynolds number, the outlet boundary condition approaches a constant pressureoutlet condition. The initial condition can be chosen arbitrarily.It is known that there exists a steady state solution u ( x, y, z ) = [0 , , w ( x, y )] T , p ( x, y, z ) = p ( z ) to the above model problem, satisfying, ∂ w ∂x + ∂ w ∂y = 1 ν ∂p ∂z . (5) * Here ν ∂p ∂z is a negative scaled constant such that the constraint of total flux 1 issatisfied. It can be verified that the solution of eqation (5) is,(6) w ( x, y ) = K (cid:20) −
14 ( y − y ) −
14 ( x − x )+12 ∞ (cid:88) m =1 sin ( mπx )( A m e mπy + B m e − mπy )+12 ∞ (cid:88) m =1 sin ( mπy )( C m e mπx + D m e − mπx ) (cid:35) , ≤ x ≤ , ≤ y ≤ A m = − mπ ) [( − m − e − mπ − e mπ − e − mπ ,B m = 2( mπ ) [( − m − e mπ − e mπ − e − mπ ,C m = A m ,D m = B m ,K = 28 . . Given the steady state solution (6), we consider the associated linear stabilityproblem. We denote the perturbation to the base flow u as u (cid:48) (cid:15) . The velocity fieldcan be expressed as u = u + u (cid:48) (cid:15) . Similarly the pressure field p can be expressedas p = p + p (cid:48) (cid:15) . By neglecting higher oder terms when the perturbations are small,the linearized Navier-Stokes equations are, ∂ u (cid:48) (cid:15) ∂t + u · ∇ u (cid:48) (cid:15) + u (cid:48) (cid:15) · ∇ u = −∇ p (cid:48) (cid:15) + ν ∇ u (cid:48) (cid:15) , in Ω(7) ∇ · u (cid:48) (cid:15) = 0 , in Ω . Boundary conditions are, inlet condition : u (cid:48) (cid:15) = , z = 0 , f or all t ≥ wall conidtion : u (cid:48) (cid:15) = , x = 0 , x = 1 , y = 0 , y = 1 , f or all t ≥ outlet condtion : − p (cid:48) (cid:15) n + ν ∂ u (cid:48) (cid:15) ∂ n = 0 , z = L, f or all t ≥ , (10)and initial condition is, u (cid:48) (cid:15) = f (cid:48) (cid:15) ( x, y, z ) , when t = 0 , f (cid:15) satisf ies above B.C.. Eigenvalue problem formulation.
We now turn to study the spectral de-composition of the above linear system. The mode can be expressed as, u (cid:48) (cid:15) = e λt u (cid:48) ( x, y, z ) , (11) p (cid:48) (cid:15) = e λt p (cid:48) ( x, y, z ) . For simplicity, we will omit the prime (cid:48) of the velocity perturbation mode u (cid:48) and thepressure perturbation mode p (cid:48) from now on. The equations of the corresponding N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT5 eigenvalue problem after substituting the velocity and pressure fields are, λ u = − u · ∇ u − u · ∇ u − ∇ p + ν ∇ u (12) ∇ · u = 0Here 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 λ determines whether each mode is asymptotically stableor 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 isdecaying, or when the flow is unstable as the perturbation grows exponentially.2.3. Divergence-free subspace and weak form.
The classical approach in find-ing the solutions of equations (12) is to solve the velocity mode u in functional space V = (cid:8) u ∈ [ H (Ω)] | u satisf ies boundary condition (8) − (10) (cid:9) and the pressuremode in Q = (cid:8) p | p ∈ L (Ω) / R (cid:9) . We now try to find the solution of u in its sub-space V = { u ∈ V |∇ · u = 0 } . When limiting the solution space in the subspace V , the incompressiblity constraint is no longer needed. Thus a weak form can beexpressed from equation (12), f ind u ∈ V , λ ∈ C s.t,λM ( u , w ) = A ( u , w ) + C ( u , w ) , ∀ w ∈ V , (13) M ( u , w ) = (cid:10) u i , w i (cid:11) ,A ( u , w ) = − ν (cid:10) ∇ k u i , ∇ k w i (cid:11) ,C ( u , w ) = (cid:10) − u k ∇ k u i − u k ∇ k u i , w i (cid:11) , || u || L (Ω) = 1 . The inner-product is defined as the integral over Ω, i.e., (cid:104) f, g (cid:105) = (cid:82) Ω f gd Ω. Notethe pressure term and boundary integrals are eliminated as a result of the specifiedboundary conditions (8)-(10). The eigenfunctions are normalized such that its L norm is unity.2.4. Finite element formulations. A Q − P mixed form is used in this paper.It is well-known that the trilinear velocity-constant pressure element suffers from acheckerboard pressure mode on regular meshes. The spurious pressure mode is dueto the fact that this element does not satisfy the Babˇuska-Brezzi(B.B.) condition.However, if the velocity field is the only interested variable, the spurious pressuremode is irrelevant, since pressure term is eliminated in equation (12). Boffi etal. ([26]) have shown that such mixed form is valid in analyzing eigenvalue problemwith the Galerkin method.The finite element formulation first constructs a weakly-divergence-free basiswith Q elements. Then a discretized weak form from equation (13) is formulated.Let V h and Q h be the finite dimensional subspace of V and Q . The weaklydivergence-free subspace follows as, V h = { u h ∈ V h | (cid:104)∇ · u h , q h (cid:105) = 0 , ∀ q h ∈ Q h } (14) = Span ( Φ , . . . , Φ N )Here the divergence-free basis function Φ i belongs to the Q trilinear element. Thetrial function q h is constant in each element. The discrete subspace is equivalent as * Figure 2.
Weakly-divergence-free basis
Figure 3.
Reference elementa subspace imposing the conservation of mass condition on each element ¯Ω e , whereΩ e denotes the interior of the eth element. (cid:90) ¯Ω e ∇ · u h d Ω = (cid:90) ∂ ¯Ω e u h · n dσ = (cid:88) i =1 (cid:90) F i u h · n dσ = 0 , ∀ ¯Ω e , (15)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 Φ i can be ex-pressed as vortices lying on each face of the element, shown in figure 2. The sixfaces of a cubic element is associated with six localized vortices. The basis can bephysically viewed as a basis for a vorticity field ω = ∇ × u , in the discretized sense.The direction of the ’effective vorticity’ is normal to each element face. Note thatfor each element Ω e , the six basis functions are linearly dependent and the degreeof freedom is five. N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT7
One of the weakly-divergence-free basis can be expressed in the reference domainas, ˆΦ ( ξ, η, ζ ) = [0 , ˆ φ + ˆ φ − ˆ φ − ˆ φ , − ˆ φ + ˆ φ − ˆ φ + ˆ φ ] T , (16)where ˆ φ i ( ξ, η, ζ ) is the typical trilinear Q shape function, with nodal value 1 oneach of the element vertex. The reference element domain is shown in figure 3, thedomain ˆ T = { ( ξ, η, ζ ) | ξ ∈ [ − , , η ∈ [ − , , ζ ∈ [ − , } . Shape function index isassociated with node index described in figure 3. The other weakly-divergence-freebasis ˆΦ i ( ξ, η, ζ ) on the reference domain can be constructed similarly.The finite element uses a simplest uniform cubic element discretizing domain Ωthroughout this paper. The cube has identical edge length h in each direction. Asa result, the coordinate transformation between ( x, y, z ) and ( ξ, η, ζ ) has a simplelinear relationship, x = h ξ + x ,y = h η + y ,z = h ζ + z , where x , y , z are coordinates of the element center. With the above assump-tions, it can be verified that the weakly-divergence-free condition is satified foreach Φ i ( x, y, z ) = ˆΦ i ( ξ, η, ζ ), (cid:90) ∂ ¯Ω e Φ i · n dσ = 0 , ∀ ¯Ω e . (17)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’ elementdescribed above, see figure 4. One of the ’half vortex’ element on the referencedomain is, ˆΦ ( ξ, η, ζ ) = [0 , ˆ φ − ˆ φ , − ˆ φ − ˆ φ ] T , (18)The index 7 is labelled after the 6 ’full vortices’ on each of the element faces. Theother basis functions at the duct outlet can be expressed in a similar way. Thetransformation to physical domain is straightforward. It can be verified that the’half vortex’ basis satisfies the weakly-divergence-free condition. There is also alinear dependency of four ’half vortices’ ˆΦ to ˆΦ and a ’full vortex’ ˆΦ on theelement outlet surface. Thus global basis Φ i mapped from ˆΦ must be excludedfrom 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 numberof elements in x , y , z directions, the total degree of freedom in the finite elementformulation is, N dof =( N x − N y − N z −
2) + ( N x − N y − N z − N x − N y −
1) + 2( N x − N y − . This is a result of the linear dependency of the six basis functions on each elementand the inclusion of boundary vortex elements. It should be noted that all vorticeswith ’effective vorticity’ along the z direction are linearly dependent with the restof the vortices, thus they must be excluded from the divergence-free basis set. Thebasis functions Φ i are constructed such that the ’effective vorticity’ has positivecomponent in each direction. This ensures the consistency when Φ i is evaluated onnearby elements. * Figure 4.
Weakly-divergence-free basis at outletAlthough the construction of the finite element scheme enforces stric cubic ele-ments, it is perhaps the simplest way to analyze the problem with minimal degreesof freedom.2.5.
Streamline upwind Petrov-Galerkin(SUPG) method.
Hughes [20] in-troduced SUPG scheme for computing time-dependent and steady state fluid prob-lems, and enjoyed great success. When element Peclet number, defined as α = h | u | / (2 ν ), where | u | is the maximum magnitude of u in each element, is large,the Galerkin method suffers from wild oscillations that contaminate the numeri-cal soltuions. In eigenvalue analysis, the same numerical instability prevales whenthe element Peclet number is large. Thus a SUPG formulation for the eigenvalueproblem is proposed to eliminate the instability.Similiar to the classical SUPG method, we take the inner-product of both sidesof equation (12) with w h and with τ u · ∇ w h in element internal domain ˜Ω = ∪ Ω e ,where Ω e is the internal of the eth element. It follows that the discretized weakform is, f ind u h ∈ V h , λ ∈ C s.t,λM ( u h , w h ) = A ( u h , w h ) + C ( u h , w h ) , ∀ w h ∈ V h , (20) M ( u h , w h ) = (cid:10) u ih , w ih (cid:11) + (cid:10) u i , τ u k ∇ k w ih (cid:11) ˜Ω ,A ( u h , w h ) = − ν (cid:10) ∇ k u ih , ∇ k w ih (cid:11) + ν (cid:10) ∇ u ih , τ u k ∇ k w ih (cid:11) ˜Ω ,C ( u h , w h ) = (cid:10) − u k ∇ k u ih − u kh ∇ k u i , w ih (cid:11) + (cid:10) − u k ∇ k u ih − u kh ∇ k u i , τ u k ∇ k w ih (cid:11) ˜Ω , || u h || L (Ω) =1 . The parameter τ can be chosen as τ = min ( α/ , h/ (2 | u | ), see [27].Appendix gives the convergence analysis of this new scheme. Convergence rateis shown to be O ( h ) for eigenvalue and O ( h ) for eigenmode energy norm. N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT9
Figure 5.
Convergence plot of leading eigenvaluesThe weak form (20) is equivalent to a generalized eigenvalue problem, H c = λM c , (21) H ij = A ( Φ j , Φ i ) + C ( Φ j , Φ i ) ,M ij = M ( Φ j , Φ i ) . The eigenmode u to be solved is u = (cid:80) i c i Φ i .The implicitly-restarted Arnoldi method, provided by ARPACK, or MATLB eigsfunction 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 − L = 3. Regular mesh size of N x = N y = N z/L is adopted.Figure 5 shows the mesh convergence plot. Computed least-stable eigenvlaue λ , λ and the third eigenvalue λ tend to converge as the degrees of freedom in-creases. The following analyasis will all adopt the case with N x = 40, which is thelargest degree of freedom computed so far.The computed eigenvalues show that the largest eigenvalue has multiplicity of2. Acutually the first two eigenmodes are identical by interchanging x, y coordi-nates and u, v components. The third eigenmode has a multiplicity of 1 and thiseigenmode possesses x - y symmetry.Figure 6 shows the three least-stable eigenvalues as a function of Re . It is againfound that the duct Poiseuille flow is linearly stable for all Re up to 2500, whichagrees with previous pipe or duct theoretical studies.It is interesting to investigate the leading-order eigenmode structures. Figure 7shows the normalized eigenmode corresponding to the first eigenvalue λ = − . Re = 2000. Figure 8 shows the normalized eigenmode corresponding to thethird eigenvalue λ = − . Re . Mesh size with N x × N y × N z =40 × ×
120 are used to compute the solutions. The cut planes are chosen atplaces where the mode is dominant. The second mode is not plotted because it isa symmetric mode of the first one. Although complicated in the structures of theshape, 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 walland are convected downstream. When the pipe inlet profile is fixed to the base flow Figure 6.
Least-stable eigenvalues as a function of Re Figure 7.
Eigenmode corresponding to the first eigenvalue λ = − . Re = 2000, nondimensional duct length is L = 3 withmesh N x × N y × N z = 40 × × || u || = 1. First row, eigenmode u, v, w in x - y cut plane at z = 2 . u, v, w in y - z cut plane at x = 0 . u, v, w in x - z cut plane at y = 0 . u , least-stable perturbations tends to dominate near the outlet. It can alsobe observed that these modes do not have a simple normal mode structure, e.g., u ( x, y, z ) = Ψ ( x, y ) e ikz , (22)where k is the wave number and Ψ is a common shape function. Therefore thissimple theoretical assumption may not be the best way to analyze pipe or duct flowstability.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 9shows the side view of the first eigenmode w ( x, y, z = 2 .
9) corresponding to λ at N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT11
Figure 8.
Eigenmode corresponding to the third eigenvalue λ = − . Re = 2000, nondimensional duct length is L = 3 withmesh N x × N y × N z = 40 × × || u || = 1. First row, eigenmode u, v, w in x - y cut plane at z = 2 . u, v, w in y - z cut plane at x = 0 . u, v, w in x - z cut plane at y = 0 . Re = 500 , , x ∼ . Re = 500 to x ∼ . Re = 2500. The physical consequences are discussed in thenext section. It seems true that ’instability’ tends to be generated near the walland 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 Poiseuilleflow is studied.It has been investigated in [15] that the instability of pipe Poiseuille is attributedto the sensitivity of the linear operator L from the linearization of the Navier-Stokesequations, where L u = − u · ∇ u − u · ∇ u − ∇ p + ν ∇ u . It was shown that when Re increase, a small perturbation to the linear operator L results in a relativelylarge eigenvalue deviation. As a result, pipe Poiseuille flow may become unstableat finite Re although in ideal case it is linearly stable. We connect this claim tothe 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 expectedthat a small fluctuation of the duct wall configuration will alter the eigenmodeshape in a significant way. The base flow field, however, is not affected greatly bywall fluctuations since its velocity field is dominant at the bulk. With this physicalpicture in mind, a perturbation to L induced by small fluctuation of duct wall will Figure 9.
Side view of the first eigenmode w ( x, y, z = 2 . Re = 500. Central figure, mode at Re = 1500.Bottom figure, mode at Re = 2500. Position x where the modepeaks are labelled.change the eigenvalue noticeablely as a result of the singularity of mode shapes inthe limit of high Re . This gives a physical point of view regarding the pseudospectratheory.To investigate the nonlinear effects, we assume the initial condition to be one ofthe least-stable modes (cid:15) u , where || u || = 1 corresponding to λ . The instananeousflow energy growth is investigated. Physically it can be regarded as the study ofthe nonlinear dynamics for a short period of time. Multiplying u into equation (1)and integrate on both sides gives,12 ∂∂t (cid:104) u , u (cid:105) = A ( u , u ) + C ( u , u ) − (cid:104) u · ∇ u , u (cid:105) , (23) N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT13 where A , C are defined in equation (13). Since the initial condition is chosen tobe the least-stable eigenmode (cid:15) u , it follows,12 ∂∂t (cid:104) u , u (cid:105) (cid:12)(cid:12)(cid:12)(cid:12) t =0 = (cid:15) λ (cid:104) u , u (cid:105) − (cid:15) (cid:104) u · ∇ u , u (cid:105) , (24)where λ is the least-stable eigenvalue. Thus the total flow energy (cid:104) u , u (cid:105) growthrate at instant t = 0 is driven by the stabilizing linear effect and the nonlineareffect. By switching the sign of the perturbation amplitude (cid:15) , the nonlinear energyproduction rate can always be positive, we therefore consider the nonlinear term tobe destabilizing in the worst case.If the nonlinear effect is negligible, the flow energy decays. The assumption totake the eigenmode as the initial condition physically makes sense when flow isattracted to the Poiseuille base flow branch. On the other hand, when (cid:15) is large,the nonlinear energy production rate may exceed the linear stabilizing effect andthe flow energy starts to grow at t = 0. The linearization assumption breaksdown in this situation. Thus we study the critical amplitude (cid:15) at various Re .Perturbation with above the critical value may not result in a final transistion intoturbulence since nonlinear dynamics after t = 0 will be complicated and is notanalyzed. However, it is a sign that linear stability problem fails to give a propermathematical model and the flow tends to be dominated by nonlinear effects. Withthe difficulty in quantifying subsequent nonlinear dynamics, the analysis only gives arough estimate on how small perturbations may generate non-negligible destablizingeffects. It should be noted that a eigenmode initial perturbation is different fromthe other studies to induce impulsive or periodic distrubances and this is not verylikely to be physically induced.The critical amplitude is determined by a differential balance equation, (cid:15) = inf | λ u i u i d Ω || u k ∇ k u i u i d Ω | . (25)When the energy production term | u k ∇ k u i u i d Ω | exceeds the linear stablizing term λ u i u i d Ω in any infinitesimal volume, the linearzation assumption is no longervalid. Thus the critical amplitude of the flow energy can be computed numericallyby finite element integration on each hexahedral element. Figure 10 shows thecritical perturbation amplitude (cid:15) vs Re , the effective nonlinear energy production η = sup | u k ∇ k u i u i d Ω || λ u i u i d Ω | vs Re , and the first eigenvalue λ vs Re . The threshold isseen to be decreasing greatly as Re increases. Any reasonable perturbation at veryhigh Re tend to have non-negligible nonlinear effect. A fit to the (cid:15) − Re curvegives approximately (cid:15) ∼ O ( Re − . ). This number is different from the resultsreported by [14], [16] and [17], mainly because the initial condition given here isnot applicable to the other studies. It may also attribute to the difference betweenpipe and duct flow. The significance of the nonlinear analysis demonstrates thatthe linear stability assumption fails to be a good model problem as Re becomeshigher. Two reasons contribute to the breakdown. First, the decreasing absolutevalue of first eigenvalue as Re increass. Second, the increasing nonlinear energyproduction rate as a result of the growing eigenfunction gradient near the ductwall, as it is found that the maximum nonlinear energy production is near the peakof the boundary layer. Figure 10. (a). Critical perturbation amplitude (cid:15) vs. Re . (b).Effective nonliner energy production rate η vs. Re . (c). Leading-order eigenvalue λ vs. Re .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 flowis 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 forthe sensitivity of the linearized Navier-Stokes problem with respect to geometryvariation as well as the sensitivity of nonlinear effects.6. Appendix
The convergence analysis in this paper assumes that the domain boundary ∂ Ωis smooth enough, and also that the base flow u , the eigenmode u and the corrre-sponding pressure field p are smooth enough, for simplicity.Space V = (cid:8) u ∈ [ H (Ω)] | u satisf ies boundary condition (8) − (10) (cid:9) and Q = (cid:8) p | p ∈ L (Ω) / R (cid:9) are now assumed to be complex spaces. V h ⊂ V , Q h ⊂ Q , V = { u ∈ V |∇ · u = 0 } and V h = { u h ∈ V h | (cid:104)∇ · u h , q h (cid:105) = 0 , ∀ q h ∈ Q h } also ex-tend to complex spaces accordingly. Lemma 6.1.
Let a ( u , w ) = − A ( u , w ) + γ (cid:104) u , w (cid:105) c ( u , w ) = − C ( u , w ) , where u , w ∈ V . A , C are defined in equation (13). There exists a sufficientlylarge γ > such that, ∃ α > , a ( u , u ) + c ( u , u ) ≥ α || u || Proof.
First we have,(26) − A ( u, u ) = ν (cid:10) ∇ k u i , ∇ k u i (cid:11) > . Next (cid:104) u · ∇ u , u (cid:105) is estimated, (cid:104) u · ∇ u , u (cid:105) = (cid:90) Ω u k ∇ k u i u i d Ω(27) = (cid:90) Ω ∇ k ( u k u i u i ) d Ω − (cid:90) Ω u k u i ∇ k u i d Ω . N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT15
Here ∇ k u k = 0 is implied. From equation (27), (cid:104) u · ∇ u , u (cid:105) = 12 (cid:90) Ω ∇ k ( u k u i u i ) d Ω(28) = 12 (cid:90) ∂ Ω u k u i u i n k dσ. Here n is the unit outward normal vector. Since the problem assumes outflow atthe outlet, namely u k n k | outlet >
0, and at the inlet and wall u i = 0, it follows that, (cid:104) u · ∇ u , u (cid:105) > . (29)Then (cid:104) u · ∇ u , u (cid:105) is estimated. Since the base flow field u and its gradient isbounded, there exists a constant tensor N ki such that, |∇ k u i | ≤ N ki . (30)Then, |(cid:104) u · ∇ u , u (cid:105)| = (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) Ω u k ∇ k u i u i d Ω (cid:12)(cid:12)(cid:12)(cid:12) (31) ≤ (cid:90) Ω N ki | u k || u i | d Ω < γ || u || , Here γ is a sufficiently large positive number. It is apparant that when γ > γI − N is positive definite.It then follows, ∃ α > , (cid:104) u · ∇ u , u (cid:105) + γ (cid:104) u , u (cid:105) > α (cid:104) u , u (cid:105) . (32)Combining equations (26), (27), (32), ∃ α > a ( u , u ) + c ( u , u ) ≥ α || u || . (cid:3) Lemma 6.2.
Let Ω e be the interior of the eth element, ˜Ω = ∪ Ω e , and, a supg ( u h , w h ) + c supg ( u h , w h )(33) = (cid:104) ν ∇ u h , ∇ w h (cid:105) + (cid:104) u · ∇ u h , w h (cid:105) + (cid:104) u h · ∇ u , w h (cid:105) + (cid:104) γ u h , w h (cid:105)− (cid:10) ν ∇ u h , τ u · ∇ w h (cid:11) ˜Ω + (cid:104) u · ∇ u h , τ u · ∇ w h (cid:105) ˜Ω + (cid:104) u h · ∇ u , τ u · ∇ w h (cid:105) ˜Ω + (cid:104) γ u h , τ u · ∇ w h (cid:105) ˜Ω . where u h , w h ∈ V h . V h is the discretized continuous piecewise polynomial solutionspace that satifies boundary conditions (8)-(10). τ = h | u | ζ ( α ) ,α = h | u | ν . Here α is the element Peclet number. ζ ( α ) satisfies m = sup α ζ ( α ) α ≤ C , where C is the inverse estimate constant that satifies, || ∆ v h || ˜Ω ≤ Ch − ||∇ v h || Ω . (34) h is the mesh parameter. In the regular cubic element case, it is the edge length. | u | denotes the maximum of base flow velocity magnitude within each element. τ varies across each element. The subscript ˜Ω means that the inner-product is takenin the interior of each element and then summing the elementwise integration up.Then ∃ γ > defined in lemma 6.1, and ∃ δ > , when τ < δ , the followingcoercive condition holds, ∃ α > , a supg ( u h , u h ) + c supg ( u h , u h ) ≥ α || u h || + 16 || τ / u · ∇ u h || . Proof. a supg ( u h , u h ) + c supg ( u h , u h )(35) = (cid:104) ν ∇ u h , ∇ u h (cid:105) + (cid:104) u · ∇ u h , u h (cid:105) + (cid:104) u h · ∇ u , u h (cid:105) + (cid:104) γ u h , u h (cid:105)− (cid:10) ν ∇ u h , τ u · ∇ u h (cid:11) ˜Ω + (cid:104) u · ∇ u h , τ u · ∇ u h (cid:105) ˜Ω + (cid:104) u h · ∇ u , τ u · ∇ u h (cid:105) ˜Ω + (cid:104) γ u h , τ u · ∇ u h (cid:105) ˜Ω . First we recall the inverse estimate relationship in equation (34), ντ = νh | u | ζ ( α ) = 14 h ζ ( α ) α ≤ h m ≤ h C (36) ⇒ ν || τ / ∇ u h || ≤ ν ||∇ u h || With this inequality, − (cid:10) ν ∇ u h , τ u · ∇ u h (cid:11) ˜Ω is estimated, | − (cid:10) ν ∇ u h , τ u · ∇ u h (cid:11) ˜Ω | (37) ≤ ν || τ / ∇ u h || + 13 || τ / u · ∇ u h || ≤ ν ||∇ u h || + 13 || τ / u · ∇ u h || . Then (cid:104) u h · ∇ u , τ u · ∇ u h (cid:105) ˜Ω is estimated, | (cid:104) u h · ∇ u , τ u · ∇ u h (cid:105) ˜Ω | (38) ≤|| τ / u h · ∇ u || + 14 || τ / u · ∇ u h || ≤|| N τ / u h || + 14 || τ / u · ∇ u h || . Here N is a constant due to the boundedness of ∇ u .Next we estimate (cid:104) γ u , τ u · ∇ u h (cid:105) ˜Ω , | (cid:104) γ u h , τ u · ∇ u h (cid:105) ˜Ω | (39) ≤|| γτ / u h || + 14 || τ / u · ∇ u h || . Subsituting all the above estimations into equation (35) and use the results inlemma 6.1, a supg ( u h , u h ) + c supg ( u h , u h )(40) ≥ ν ||∇ u h || + 16 || τ / u · ∇ u h || + (cid:10) ( α − N τ − γ τ ) u h , u h (cid:11) ˜Ω . N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT17
Here α is the estimation constant in lemma 6.1. It follows that when τ < α N + γ = δ, (41) ∃ α > , a supg ( u h , u h ) + c supg ( u h , u h ) ≥ α || u h || + 16 || τ / u · ∇ u h || . (cid:3) Remark 6.0.1.
For typical flow problem α , γ can be chosen on the order of N ∼|∇ u | , so the condition for τ is not restrictive. Remark 6.0.2.
The introduction of ’ γ ’ terms ensures the coercivity of the linearproblem. The corresponding eigenvalue problem is, a supg ( u , w ) + c supg ( u , w ) = m supg (˜ λ u , w ) , (42) m supg ( u , w ) = − M ( u , w ) . Here M ( u , w ) is a bilinear form defined in equationThe discretized problem is equivalent to a matrix eigenvalue problem, H c − γM c = ˜ λ h M c . (43) Here H , M are matrices defined in equation (20). It is apparant that the eigenvalueis shifted by γ from the original problem 21, and that the eigenfunctions remainunchanged. Lemma 6.3.
Let a ( u , w ) = (cid:10) ∇ k u i , ∇ k w i (cid:11) , ∀ u , w ∈ V , (44) b ( w , p ) = (cid:104)∇ · w , p (cid:105) , ∀ w ∈ V , p ∈ Q. Define a linear system in its weak form, f ∈ V , a ( u , w ) + b ( w , p ) = (cid:104) f , w (cid:105) , ∀ w ∈ V , (45) b ( u , q ) = 0 , ∀ q ∈ Q, Here the binlinear form b ( u , q ) satisfies the B.B. condition such that, ∃ β > , inf q ∈ Q,q (cid:54) =0 sup u ∈ V , u (cid:54) = | b ( u , q ) ||| u || V || q || Q ≥ β. The corresponding discretized weak form is, a ( u h , w h ) + b ( w h , p h ) = (cid:104) f h , w h (cid:105) , ∀ w h ∈ V h , (46) b ( u h , q h ) = 0 , ∀ q h ∈ Q h . Then whether or not the B.B. condition is satisfied or not in equations (46), thereis a unique u h that satifies the discretized equations and, || u − u h || ≤ C inf u I ∈ V h ,p I ∈ Q h ( || u − u I || + || p − p I || ) . (47) Here C is a constant that depends on u and p .Proof. It is well-known that equation (45) has a unique solution ( u , p ) ∈ V × Q ,due to the coercivity of a ( · , · ) and the satisfaction of the B.B condition in b ( · , · ).The solution u h of equation (46) is in subspace V h = Ker ( B h ), where(48) (cid:104) B h u h , q h (cid:105) = b ( u h , q h ) , ∀ w h ∈ V h , q h ∈ Q h . Thus by restricting the trial function w h ∈ V h , a ( u h , w h ) = (cid:104) f , w h (cid:105) , ∀ w h ∈ V h . (49)Then there exists a unique u h ∈ V h that satisfies equation (46). Note p h may nothave unique solution.Next we slightly tweak the discretized problem such that space Q h is replacedby ˆ Q h = Q h /Ker ( B Th ), , a ( u h , w h ) + b ( w h , ˆ p h ) = (cid:104) f , w h (cid:105) , ∀ w h ∈ V h , (50) b ( u h , ˆ q h ) = 0 , ∀ ˆ q h ∈ ˆ Q h , It is straight-forward to check that the B.B. condition is satisfied in equation (50).Therefore, there exists a unique pair ( u h , ˆ p h ) ∈ V h × ˆ Q h as the solution of equa-tion (50). Apparantly u h is also the unique solution of equation (46).It is evident that the following relationship also holds, a ( u , w h ) + b ( w h , p ) = (cid:104) f , w h (cid:105) , ∀ w h ∈ V h , (51) b ( u , ˆ q h ) = 0 , ∀ ˆ q h ∈ ˆ Q h , Subtracting equation (50) from equation (52), we have, a ( u h − u I , w h ) + b ( w h , p h − p I ) = a ( u − u I , w h ) + b ( w h , p − p I ) , ∀ w h ∈ V h , (52) b ( u h − u I , ˆ q h ) = b ( u − u I , ˆ q h ) , ∀ ˆ q h ∈ ˆ Q h . Here u I ∈ V h , p I ∈ Q h are any interpolations of u and p . Note p I is not in thetruncated subspace ˆ Q h . Becuase the B.B. condition holds in equation (52), ∃ u ∈ V h , ˆ u ∈ Ker ( ˆ B h ) , (53) u h − u I = u + ˆ u ,b ( u , ˆ q h ) = b ( u − u I , ˆ q h ) , ∀ ˆ q h ∈ ˆ Q h . Here ˆ B h is defined as, (cid:68) ˆ B h u h , ˆ q h (cid:69) = b ( u h , ˆ q h ) , ∀ u h ∈ V h , ˆ q h ∈ ˆ Q h . If the B.B condition does not hold in equation (46),
Ker ( B Th ) (cid:54) = Φ, therefore any p I ∈ Q h can be expressed as, p I = ˆ p + ˜ p, ˆ p ∈ ˆ Q h , ˜ p ∈ Ker ( B Th ) . (54)We now demonstrate that Ker ( ˆ B h ) = Ker ( B h ).First, if u h ∈ Ker ( B h ), (cid:104) B h u h , q h (cid:105) = 0 , ∀ q h ∈ Q h . (55)Because ˆ Q h ⊆ Q h , it follows Ker ( B h ) ⊆ Ker ( ˆ B h ).Next, if u h ∈ Ker ( ˆ B h ), (cid:68) ˆ B h u h , ˆ q h (cid:69) = 0 , ∀ ˆ q h ∈ ˆ Q h . (56) N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT19
Pick any q h ∈ Q h , q h = ˆ q h + ˜ q h , ˆ q h ∈ ˆ Q h , ˜ q h ∈ Ker ( B Th ) , (cid:104) B h u h , q h (cid:105) = (cid:104) B h u h , ˆ q h + ˜ q h (cid:105) = 0 , ∀ q h ∈ Q h . (57)Then we have Ker ( ˆ B h ) ⊆ Ker ( B h ).Therefore Ker ( ˆ B h ) = Ker ( B h ).From the above discussion we have, a (ˆ u , w h ) + b ( w h , p h − p I ) =(58) a ( u − u I , w h ) + b ( w h , p − p I ) − a ( u , w h ) , ∀ w h ∈ V h , u I ∈ V h , p I ∈ Q h . Let w h = ˆ u and notice b (ˆ u , p h − p I ) = 0, since ˆ u ∈ Ker ( ˆ B h ) = Ker ( B h ). It follows, a (ˆ u , ˆ u ) = a ( u − u I , ˆ u ) + b (ˆ u , p − p I ) − a ( u , ˆ u )(59)It is known that B.B. condition holds in equation (50), thus we have ∃ β > , || u || ≤ β − || ˆ B h u || . (60)Therefore, from equation (54), || u || ≤ β − sup ˆ q h ∈ ˆ Q h | b ( u , ˆ q h ) ||| ˆ q h || = β − sup ˆ q h ∈ ˆ Q h | b ( u − u I , ˆ q h ) ||| ˆ q h || ≤ Cβ − || u − u I || (61)Then, α || ˆ u || ≤ a (ˆ u , ˆ u ) ≤|| a || · || u − u I || · || ˆ u || + C || p − p I || · || ˆ u || +(62) Cβ − || a || · || u − u I || · || ˆ u || ⇒ || ˆ u || ≤ C ( || u − u I || + || p − p I || )Combining equation (61), || u h − u I || ≤ C ( || u − u I || + || p − p I || ) . (63)Finally, || u h − u || ≤ C inf u I ∈ V h ,p I ∈ Q h ( || u − u I || + || p − p I || ) . (64)Here constant C may be of different value in each equation for simplicity. (cid:3) Corollary 6.1.
For any u ∈ V , where V is the divergence-free subspace of V ,there exists an interpolation ˜ u I ∈ V h such that, || ˜ u I − u || ≤ C inf u I ∈ V h ,p I ∈ Q h ( || u − u I || + || p − p I || ) . (65)This corollary shows that the divergence-free vector field u can be approximatedby a weakly-divergence-free interpolation with the above error bounds. Whetheror not the B.B condition is satisfied or not does not matter. Theorem 6.1.
Given f ∈ V , there exists a unique pair ( u , p ) ∈ V × Q such that, a ( u , w ) + c ( u , w ) + b ( w , p ) = (cid:104) f , w (cid:105) , ∀ w ∈ V , (66) b ( u , q ) = 0 , ∀ q ∈ Q. Here the binlinear form b ( u , q ) satisfies the B.B. condition such that, ∃ β > , inf q ∈ Q,q (cid:54) =0 sup u ∈ V , u (cid:54) = | b ( u , q ) ||| u || V || q || Q ≥ β. Then the following discretized linear system also has a unique solution u h ∈ V h , a supg ( u h , w h ) + c supg ( u h , w h ) = (cid:104) f , w h (cid:105) + (cid:104) f , τ u · ∇ w h (cid:105) ˜Ω , ∀ w h ∈ V h . (67) Here the discretized binlinear form b ( u h , q h ) may not satisfy the B.B condition.If we assume V h to be continous piecewise polynomials of degree k and Q h be ofdegree k − , the following esitmation holds, || u − u h || ≤ C u,p h . (68) Here C u,p is a function of u and p .Proof. It is well-known that when B.B. condition is satisfied, equation (66) has aunqiue solution because the coercivity holds true.The unique solution ( u , p ) ∈ V × Q also satisfies, a supg ( u , w h ) + c supg ( u , w h ) + b ( w h , p )+(69) (cid:104) τ u · ∇ w h , ∇ p (cid:105) ˜Ω = (cid:104) f , w h (cid:105) + (cid:104) f , τ u · ∇ w h (cid:105) ˜Ω , ∀ w h ∈ V h ,b ( u , q h ) = 0 , ∀ q h ∈ Q h . Combining equation (67) and (69), a supg ( u h , w h ) + c supg ( u h , w h ) =(70) a supg ( u , w h ) + c supg ( u , w h ) + b ( w h , p ) + (cid:104) τ u · ∇ w h , ∇ p (cid:105) ˜Ω , ∀ w h ∈ V h . Recall V h is the discrete divergence-free subspace. Let u I ∈ V h , be any inter-polation of u , and p I ∈ Q h be any interpolation of p . Note b ( w h , p I ) = 0. Wehave, a supg ( u h − u I , w h ) + c supg ( u h − u I , w h ) = a supg ( u − u I , w h ) + c supg ( u − u I , w h )(71) + b ( w h , p − p I ) + (cid:104) τ u · ∇ w h , ∇ p (cid:105) ˜Ω , ∀ w h ∈ V h . N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT21
Let w h = u h − u I ∈ V h , from lemma 6.2, α || u h − u I || + 16 || τ / u · ∇ ( u h − u I ) || (72) ≤ (cid:10) ν ∇ ( u − u I ) , ∇ ( u h − u I ) (cid:11) + (cid:10) u · ∇ ( u − u I ) , u h − u I (cid:11) + (cid:10) ( u − u I ) · ∇ u , u h − u I (cid:11) + (cid:10) γ ( u − u I ) , u h − u I (cid:11) − (cid:10) ν ∇ ( u − u I ) , τ u · ∇ ( u h − u I ) (cid:11) ˜Ω + (cid:10) u · ∇ ( u − u I ) , τ u · ∇ ( u h − u I ) (cid:11) ˜Ω + (cid:10) ( u − u I ) · ∇ u , τ u · ∇ ( u h − u I ) (cid:11) ˜Ω + (cid:10) γ ( u − u I ) , τ u · ∇ ( u h − u I ) (cid:11) ˜Ω . + b ( u h − u I , p − p I ) + (cid:10) τ u · ∇ ( u h − u I ) , ∇ p (cid:11) ˜Ω ≤ α || u h − u I || + C || u − u I || + 112 || τ / u · ∇ ( u h − u I ) || + C ν || τ / ∇ ( u − u I ) || + C || p − p I || + C || τ / ∇ p || ⇒ α || u h − u I || ≤ C || u − u I || + C ν || τ / ∇ ( u − u I ) || + C || p − p I || + C || τ / ∇ p || Therefore, ∃ C > , (73) || u − u I || ≤ C ( || u − u I || + ν || τ / ∇ ( u − u I ) || + || p − p I || + || τ / ∇ p || ) . Since u ∈ V and u I ∈ V h , from Corollary 6.1, ∃ u I ∈ V h , || u − u I || ≤ C inf ˜ u I ∈ V h ,p I ∈ Q h ( || u − ˜ u I || + || p − p I || ) ∼ O ( h k ) , (74)From the choice of τ = h u min ( α , τ = O ( h | u | ) , P eclet number α is large, (75) τ = O ( h ν ) , P eclet number α is small. Remember α = | u | h ν in each element.We therefore have, (cid:40) ν τ = ν h | u | = h | u | α ≤ h | u | ≤ Ch , α large, | u | bounded,ν τ = νh ≤ Ch , α small,ν || τ / ∇ ( u − u I ) || ≤ C u h l , (76) 2 l = (cid:40) k + 1 , α large, k, α small, where C ’s are constant numbers in each estimation and C u is a function of u .We also have, (cid:40) τ = h | u | = h να ≤ Ch ν , α large,τ = h ν ≤ Ch ν , α small, (77) || τ / ∇ p || ≤ C p h . Combining all results above and notice the leading order error is O ( h ) when k ≥ || u − u I || ≤ C u,p h . (78)Note all the C u , C p are functions of u and p . Since given f , the pair ( u , p ) is unique.We can also write that, || u − u I || ≤ C f h . (79)The leading-order error comes from (cid:10) τ u · ∇ ( u h − u I ) , ∇ p (cid:11) ˜Ω . However, thisscheme does not degrade the accuracy of Q − P element. It will be interesting tofind a way to increase the order of convergence in the future. (cid:3) Corollary 6.2.
Given the linear systems, a ( T f , w ) + c ( T f , w ) = (cid:104) f , w (cid:105) , ∀ w ∈ V , (80) and, a supg ( T h f , w h ) + c supg ( T h f , w h ) = (cid:104) f , w h (cid:105) + (cid:104) f , τ u · ∇ w h (cid:105) ˜Ω , ∀ w h ∈ V h , (81) The operator T : V → V , T h : V → V h is well-defined, with the norm bound, || ( T − T h ) | E ( µ ) || L ( V ) ≤ Ch. (82)
Here T | E ( µ ) means the restriction to the generalized eigenspace corresponding toeigenvalue µ of operator T . Here µ = ˜ λ − , where ˜ λ is defined in equation (42),and, ˜ λT u = u . (83)The convergence of the flow eigenvalue problem is a direct result of Babˇuska andOsborn [28], Theorem 6.2. (Babˇuska and Osborn) The distance of eigenspace E ( µ ) and thediscretized eigenspace E h ( µ ) satisfies the following bound, ˆ δ ( E ( µ ) , E h ( µ )) ≤ C || ( T − T h ) | E ( µ ) || L ( V ) . (84) Here ˆ δ ( E ( µ ) , E h ( µ )) = max ( δ ( E ( µ ) , E h ( µ )) , δ ( E h ( µ ) , E ( µ ))) , (85) δ ( E , F ) = sup u ∈ E, || u || =1 inf v ∈ F || u − v || Theorem 6.3. (Babˇuska and Osborn) Let µ be a non-zero eigenvalue of T withalgebraic multiplicity equal to m and let ˆ µ h denote the arithmetic mean of the m discrete eigenvalues of T h converging towards µ . Let φ , · · · , φ m be a basis ofgeneralized eigenvectors in E ( µ ) and let φ ∗ , · · · , φ ∗ m be a dual basis of generalizedeigenvectors in E ∗ ( µ ) . Then, | µ − ˆ µ h | ≤ m m (cid:88) i =1 | (cid:104) ( T − T h ) φ i , φ ∗ i (cid:105) | + C || T − T h || L ( V ) || T ∗ − T ∗ h || L ( V ∗ ) . (86) N THE THREE-DIMENSIONAL STABILITY OF POISEUILLE FLOW IN A FINITE-LENGTH DUCT23
Theorem 6.4.
The eigenvalue problem has an error estimate for the i th eigenvalue, | ˜ λ i − ( 1 m m (cid:88) i k =1 ˜ λ − i k h ) − | ≤ C i h , (87) where m is the algebraic multiplicity of ˜ λ i and ˜ λ i k h are the m discrete eigenvaluesconverging towards ˜ λ .The estimate for the i th eigenmode is, || u i − u ih || ≤ C i h. (88) Proof.
The proof follows from Boffi [29] Theorem 10 . (cid:3) References [1] O. Reynolds,
An experimental investigation of the circumastances which determine whetherthe motion of water shall be direct or sinuous and of the law of resistance in parallel channels ,Phil. Trans. R. Soc. Lond. (1883).[2] P. Huerre, M. Rossi,
Hydrodynamic instabilities in openflows , In Hydrodynamics and non-linear instabilities(ed. C. Godreche & P. Manneville) (1998).[3] L. Boberg, U. Brosa,
Onset of turbulence in a pipe , Z. Naturforschung (1988)[4] P. L. O’Sullivan, K. S. Breuer,
Transient growth in a circular pipe. I. Linear distrubances ,Phys. Fluids (1994).[5] V. G. Priymak, T. Miyazaki,
Accurate Navier-Stokes investigation of transitional and tur-bulent flows in a circular pipe , J. Comp. Phys. (1998).[6] O. Yu. Zikanov,
On the instability of pipe Poiseuille flow , Phys. Fluids (1996).[7] L. Bergstr¨om,
Optimal growth of small disturbances in pipe Poiseuille flow , Phys. Fluids A(1993).[8] P. J. Schmid, D. S. Henningson,
Optimal energy growth in Hagen-Poiseuille flow , J. FluidMech (1994).[9] A. E. Trefethen, L. N. Trefethen, P. J. Schmid,
Spectra and pseudospectra for pipe Poiseuilleflow , Comp. Metho. Appl. Mech. Engr. (1999).[10] ´A. Meseguer, L. N. Trefethen,
Linearized pipe flow to Reynolds number , J. Comp. Phys.(2003).[11] F. Delplace, G. Delaplace, S. Lefebvre, et al. Friction curves for the flow of Newtonianand non-Newtonian liquids in ducts of complex crosssectional shape.
Proc. 7th Internationalcongress on engineering and food. (1997)[12] A. G. Darbyshire, T. Mullin,
Transition to turbulence in constant mass-flux pipe flow , J.Fluid Mech. (1995).[13] I. J. Wygnanski, F. H. Champagne,
On transition in a pipe. Part 1. The origin of puffs andslugs and the flow in a turbulent slug , J. Fluid Mech (1973).[14] B. Hof, A. Juel, T. Mullin,
Scaling of the turbulence transition threshold in a pipe , Phys.Rev. Lett. (2003)[15] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, T. A. Driscoll,
Hydrodynamic stability withouteigenvalues , Science (1993).[16] ´A. Meseguer
Streak breakdown instability in pipe Poiseuille flow , Phys. Fluids (2003).[17] B. Eckhardt, A. Mersmann,
Transition to turbulence in a shear flow , Phys. Rev. E (1999).[18] T. J. R. Hughes, A. N. Brooks,
A multidimensional upwind scheme with no crosswind dif-fusion , in: T.J.R. Hughes, ed., Finite Element Methods for Convection Dominated Flows(ASME, New York, 1979).[19] T. J. R. Hughes, A. N. Brooks,
A theoretical framework for Petrov-Galerkin methods withdiscontinuous weighting functions: application to the streamline upwind procedure , in: R. H.Gallagher, D. H. Norrie, J. T. Oden and O. C. Zienkiewicz, Finite Elements in Fluids, Vol.IV (Wiley, London, 1982).[20] T.J.R. Hughes,
Recent progress in the development and understanding of SUPG methodswith special reference to the compressible Euler and Navier-Stokes equations , Internat. J.Numer. Methods Fluids (1987). [21] D. F. Griffiths,
Finite element for incompressible flow , Math. Methods Appl. Sci. (1979).[22] D. F. Griffiths,
The construction of approximately divergence-free finite element , in: TheMathematics of Finite Element and its Applications, Vol. 3, ed. J.R. Whiteman (AcademicPress, New York, 1979).[23] D. F. Griffiths,
An approximately divergence-free 9-node velocity element for incompressibleflows , Internat. J. Numer. Methods Fluids (1981).[24] M. Fortin,
Old and new finite elements for incompressible flows , Internat. J. Numer. MethodsFluids (1981).[25] X. Ye, C. A. Hall,
A discrete divergence-free basis for finite element methods , NumericalAlgorithms (1997).[26] D. Boffi, F. Brezzi and L. Gastaldi,
On the convergence of eigenvalues for mixed formulations ,Ann. Scuola Norm. Sup. Pisa Cl. Sci. (1997)[27] T.J.R. Hughes, M. Mallet and A. Mizukami,
A new finite element formulation for computa-tional fluid dynamics: II. Beyond SUPG , Comput. Meths. Appl. Mech.(1986)[28] I. Babˇuska and J. Osborn,
Eigenvalue problems , In Handbook of Numerical Analysis, Vol. II,North-Holland, Amsterdam (1991).[29] D. Boffi,