Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 45, 3, pp. 659-683, Warsaw 2007 INSTABILITY OF VISCOUS INCOMPRESSIBLE FLOW IN A CHANNEL WITH TRANSVERSELY CORRUGATED WALLS Jacek Szumbarski Warsaw University of Technology, Institute of Aeronautics and Applied Mechanics, Poland e-mail: jasz@meil.pw.edu.pl Destabilization of the laminar flow in a channel with transversly wavy walls is investigated. The basic flow is determined by a semi-analytical methodbased on the concept of immersedboundaries.The stability equ- ations are discretized by the spectrally accurateChebyshev-taumethod. It is demonstrated that the transversalwallwaviness canbeused inorder to achieve destabilization of the laminar flow at lowReynolds numbers. Possible applications include the improvement of mixing and heat/mass transfer in various devices used in heat technology, biotechnology and medicine. Key words: wavy channel, laminar flow, linear stability 1. Introduction The enhancement of mixing is of crucial importance for efficient operation of such devices as compact heat exchangers, cooling systems in microelectro- nics, small bioreactors, oxygenators, dialisators and similar. Flows inside such devices are usually characterized by low or – at most – moderate Reynolds numbers (say, of the order of hundreds), thus they tend to be laminar with poor transport properties. The obvious way to improve themixing is flow tur- bulization, however it can be hard to achieve. In some situations – like flows of suspensions of shear-sensitive biological material or blood – turbulization may be even undesirable. Then, improvement of themixingmust be achieved by an appropriate flow modification of the flow kinematical structure in the laminar (or transitional) regime. It seems that the most widely analyzed method for the enhancement of mixing of internal laminar flows is the application of various variants of uni- directional wall waviness (taken in a broad sense, i.e. including such configu- rations like ribs or grooves) oriented parallel to the main flow direction. Such 660 J. Szumbarski wall waviness is referred here as the longitudinal one. In particular, in wall waviness presented in Fig.1 would be considered longitudinal if the direction of the main flow was parallel to the 0x axis of the reference frame. Fig. 1. A channel with sinusoidal wall waviness oriented transversely to the flow. The corrugation amplitude of both walls is the same while the phase is opposite. Thus, the flow domain is symmetric with respect to the plane y=0 During the last three decades, many variants of such geometrical modi- fications have been investigated both numerically and experimentally. First comprehensive experimental investigations of kinematical structures andmass transfer in the flow through a symmetric divergent-convergent channel with arc-shaped walls were carried out by Goldstein and Sparrow (1977). Numeri- cal investigations of two-dimensional flow in adivergent-convergent symmetric channel (in two dimensions) with sinusoidal walls were performed by Sobey (1980) and verified experimentally by Stephanoff et al. (1980). Systematic in- vestigations of flows in corrugated channels have been carried out since the middle 80s in Japan, mostly by Nishimura and his co-workers. They inve- stigated a variety of kinematical structures, which appear in flows through channels with sinusoidal walls driven by steady (Nishimura et al., 1990), pul- sating (Nishimuraet al., 1989;NishimuraandKawamura, 1995) andoscillatory pressure gradient (Nishimura et al., 1987; Nishimura andKojima, 1995). The basic conclusion from these works is that the significant improvement of mi- xing properties can be achieved when the natural flow separation phenomena cooperate with forced flow oscillations. Later, in numerical works by Ghad- dar and El-Hajj (2000), Niceno and Nobile (2001), Wang and Vanka (1995) it was shown that a significant improvement of heat or mass transfer occurs only when the wall corrugation induces flow instability and the velocity field becomes oscillatory. For sinusoidalwavywalls with large amplitudes (reaching 27% of the average channel height), the oscillatory mode of the flow appears for the Reynolds number over 200 and more than threefold heat transfer en- hancement (defined as the ratio between average Nuselt numbers for the flow Instability of viscous incomressible flow... 661 in the corrugated and plane channels) is achieved. A serious ”penalty” for this improvement is the significantly increased hydraulic resistance, which can be even by order of magnitude larger than in the reference Poiseuille flow. This effect is especially severe when the corrugated walls have sharp edges (like in the case of the arc-shaped waviness). Several authors have recently studied stability properties of laminar inter- nal flows in the presence of a longitudinal wall corrugation. In the work by Cho et al. (1998), the linear stability of the laminar flow in the plane channel with sinusoidal walls is analyzed. Using a numerical solution to the lineari- zed equations for small disturbances, the authors have shown the existence of two modes of instability: oscillatory and non-oscillatory. The first mode is attributed to traveling wave-like disturbances, and in fact it corresponds to the fundamental Tollmien-Schlichting (TS) wave modified by the presence of the wall waviness. The non-oscillatory mode corresponds to stationary (not traveling) disturbances. The critical Reynolds numbers (based on the avera- ge velocity and average half-distance between wavy walls) for both unstable modes depend of the corrugation amplitude and the spanwise wave number. For large corrugation amplitudes, two-dimensional TS waves are destabilized first, and the critical Reynolds number can drop to values as low as 100. The stationary disturbances are three-dimensional and for lower corrugation am- plitudes their critical Reynolds number is smaller than for 2D TSwaves. The effect of instability of TS waves on the effectiveness of heat transfer in the flow through a divergent-convergent channel was investigated numerically by Blancher et al. (1998). Recently, the same authors have carried out numerical analysis of spatially developing nonlinear instabilities in such a flow (Blancher et al., 2004). Someotherworkswere originallymotivated by the interest in the transitio- nal phenomena in flows over roughwalls rather than themixing enhancement. For that reason, the amplitudes of the corrugation considered there are rather small and thus application of the perturbationmethods is justified. Examples of such an approach can be found in Floryan et al. (2002), Selvarajan et al. (1999) or, recently,Wu and Luo (2006). Another possibility to investigate the stability offlowsover slightly corrugatedwalls is tomodel thegeometrymodifi- cation bywall-distributed transpirations (blowing and suction). Using this ap- proach, Floryan (1997) demonstrated that even small-amplitude, periodically distributed wall transpiration can significantly reduce the stability thresholds and that the unstable modes have a form of stationary streamwise-oriented vortices. Later, Szumbarski (2002c) showed that there exist two unstable mo- des and also the explained their origin. A computationally convenient setting for the investigation of a small or moderate wall corrugation is provided by the spectrally accuratemethod of immersed boundaries proposed originally in Szumbarski and Floryan (1999) and later extended by Szumbarski (2002b) to 662 J. Szumbarski linear stability equations. Thismethodwas successfully applied in the further works investigating the instability of the wavy-channel flows (Cabal et al., 2002; Floryan, 2003, 2005; Szumbarski, 2002b) and, recently, also the pheno- menon of transient growth of disturbances (Szumbarski and Floryan, 2006). The results of these works confirmed large potential of the small-amplitude wall waviness for flow destabilization. It has been established that wall cor- rugation with the amplitude equal only 2% of the average channel height can reduce the critical Reynolds number from 5772 (for the reference Poiseuille flow) to about 1000. They also allow one to properly interpret some results of the previous research. In particular to identify the 3-dimensional stationary modes discovered byCho et al. (1998) as certain Squiremodesmodifiedby the presence of the wall waviness. High sensitivity of the critical conditions of the channel flow to small modifications such as a wall corrugation/transpiration or weak pulsations of the pressure gradient has been recently demonstrated byWu and Luo (2006). Surprisingly enough, the number of works investigating the stability of la- minar flows in transversely corrugated channels is very limited. Such a type of geometricalmodifications attracted someattention in the late 80’s,mostly due to relations to surface elements known as ”riblets”. Most works on riblets are focused on the effect of reduction of the friction drag in turbulent boundary layers, see Viswanath (2002) for the recent review of the research in this field. Investigations of the laminar-turbulent transition in the boundary layer over riblet-mounted surfaces were carried in 90’s, notably by Luchini (1995), Lu- chini andTrombetta (1995), Grek et al. (1996), Ehrenstein (1996). It has been established that the riblets affect the natural transition in the laminar boun- dary layer in the following manner: they accelerate the initial growth of the Tollmien-Schlichtingwaves but further nonlinear stage of transition (TSwaves secondary instability, formation of lambda vortices and turbulent spots) is de- layed and occurs farther downstream than in the case of a smooth surface. In the paper byEhrenstein (1996) it was shown that the transversal short-period wall waviness destabilizes the originally two-dimensional Tollmien-Schlichting waves in the Poiseuille flow. For larger corrugation amplitudes, the critical Reynolds number of these disturbances can drop below 2600. This result is not very spectacular and actually seems to suggest that the transversal cor- rugation is by no means competitive to the longitudinal one as the method of the mixing enhancement. Indeed, similar reduction of the critical Reynolds number can be achieved with the use of the longitudinal waviness with the amplitude smaller by the order of magnitude and essentially without any no- ticeable increase in the flow resistance. It is true indeed that the transversal waviness is much less efficient in destabilizing Tollmien-Schlichting waves that the longitudinal one. However, existence of someother formsof disturbances,whichmaybedestabilizedmuch Instability of viscous incomressible flow... 663 stronger than TS waves, cannot be excluded. We will show that such distur- bances really exist and dramatic reduction of the critical Reynolds number of the flow can be achieved. The paper is organized as follows. In Section 2, the problem of determina- tion of the basic (undisturbed) flow is formulated and the semi-analytical me- thod based on the concept of immersed boundaries is described. In Section 3, the linear stability of the basic flow is investigated. An unstable normal mo- de is identified and its parametric variation is studied numerically. Kinematic structures in the disturbance velocity field corresponding to the unstable mo- de are also presented. Finally, in Section 4, a short summary is provided and some conclusions are formulated. 2. Basic flow 2.1. Formulation of the problem Let us first consider the reference case: laminar flow of a viscous liquid in a straight channel, i.e. in a region between two parallel flat walls. The coordinate system is chosen in such a way that the walls correspond to the planes Y =−H and Y =H. Let GP < 0 denote a constant pressure gradient which drives fluidmotion in the positive direction of the Z axis. The velocity of the flow inside the channel takes the form V 0 = [0,0,W0(Y )] = [ 0,0,Wmax ( 1− Y 2 H2 )] (2.1) where µ denotes the dynamic viscosity and Wmax is the maximum velocity given by the formula Wmax =− GPH 2 2µ It is convenient to apply adimensionless description.To this end,wedefine the length scale H, the velocity scale Wmax and the pressure scale ρW 2 max, where ρ denotes the density of the fluid. Then, velocity field (2.1) is transformed to the standard Poiseuille flow w0(y)= 1−y 2 (2.2) where the dimensionless coordinate y belongs to the interval [−1,1]. Thenon- dimensional pressure gradient can be now expressed by the Reynolds number Re=WmaxH/ν as gP =GP H ρW2max =− 2 Re 664 J. Szumbarski Letus consider now the laminar flow in the transverselywavy channel depicted in Fig.1. In the dimensionless coordinates, the shape of the wavy walls is described by the following x-periodic functions yD(x)=−1+S cos(αx) yG(x)= 1−S cos(αx) (2.3) where S ­ 0 is the amplitude of wall corrugation. Thus, the channel is sym- metricwith respect to the plane y=0.Thegeometric period of thewall shape is λX =2π/α. The average height of the wavy channel is equal 2, i.e. it is the same as the wall distance in the reference straight channel. The basic flow inside the wavy channel is assumed unidirectional and its velocity field can be expressed as vB = [0,0,W(x,y)] (2.4) Clearly, the continuity equation divvB = 0 is satisfied. We assume that the flow is driven by the same pressure gradient gP as the reference Poiseuille flow. The Navier-Stokes equation for the basic flow is reduced to the following Poisson equation −gP + 1 Re ∆W =0 ⇒ ∆W =−2 (2.5) where ∆= ∂xx+∂yy. Since the function W is x-periodic, equation (2.5) should be satisfied in the domain Ω = {(x,y) : x ∈ (0,λX), yD(x) < y < yG(x)}. The boundary conditions imposed at the channel walls are W(x,yD(x))= 0 W(x,yG(x))= 0 (2.6) 2.2. Solution to the boundary value problem In order to determine the basic flow,we take upa semi-analytical approach based on the concept of immersed boundaries. To this end, we define the extended computational domain Ωext = {(x,y) : x∈R, y∈ [−1−S,1+S]}. Obviously, the inclusion Ω⊂Ωext holds. Next, we assume that the basic flow is defined for every point in the exten- ded domain Ωext. The velocity of this flow can be written as the following sum W(x,y)=w0(y)+W ′(x,y) (2.7) Expression (2.7) is admissible since formula (2.2) for the reference flow can be trivially extended to the interval [−1−S,1+S]. The component W ′(x,y) describes a modification of the reference flow due to the wall corrugation. Instability of viscous incomressible flow... 665 Equation (2.5) will be satisfied providing that the function W ′(x,y) is x- periodic and harmonic, i.e. it fulfills the Laplace equation ∆W ′ =0 in Ωext. The harmonic function W ′ can be expressed in the form of the Fourier series W ′ = ∞∑ n=−∞ Wn(y)e inαx W−n =W ∗ n (2.8) where theamplitude functions Wn are the solutions to theordinarydifferential equations (D2−n2α2)Wn =0 n=0,±1,±2, . . . (2.9) and can be expressed analytically as follows W0 =C0+S0y (2.10) Wn =Cncosh(nαy)+Sn sinh(nαy) Since the function W ′ is real, the coefficients C0 and S0 are real numbers and C−n =C ∗ n, S−n =−S ∗ n for all n 6= 0. After insertion of (2.10) into (2.8), we get the formula W ′(x,y)=C0+S0y+ ∑ n 6=0 [Cncosh(nαy)+Sn sinh(nαy)]e inαx (2.11) The coefficients {Cn} and {Sn} should be such that the boundary conditions (2.6) are satisfied. Note that the wall contours y= yD(x) and y= yG(x) are ”immersed” in the extended domain Ωext, thus conditions (2.6) have internal rather than boundary character. Themethod used to enforce these conditions is similar to the immerse boundary technique described by Szumbarski and Floryan (1999). The basic idea is to calculate explicitly theFourier coefficients of the velocity distributions at the wavy walls and then tomake them vanish. Practically, only a certain number of the leading Fourier modes can be elimi- nated. In the current study, the least squares formulation of the elimination procedure is applied. Consider the bottom wall of the channel y = yD(x). From (2.11), the velocity distribution at the bottom wall can be expressed as follows W(x,yD(x))=w0(yD(x))+C0+S0yD(x)+ (2.12) + ∑ n 6=0 [Cncosh(nαyD(x))+Sn sinh(nαyD(x))]e inαx Since both the velocity field and the wall shape are periodic with respect to the x coordinate, the above formula describes the x-periodic function εD. 666 J. Szumbarski This function can be written in the form of the Fourier expansion εD(x)≡W(x,yD(x)) = ∞∑ m=−∞ EDme imαx (2.13) Weneed toderive formulae for the coefficients {EDm}.To this end, the following Fourier expansions of all x-periodic functions appearing in equation (2.12) have to be calculated w0(yD(x)) = ∞∑ m=−∞ BDme imαx ζDm(x)≡ cosh(mαyD(x))= ∞∑ n=−∞ GDm,ne imαx (2.14) γDm(x)≡ sinh(mαyD(x)) = ∞∑ n=−∞ HDm,ne imαx Insertion of formulae (2.14) into (2.12) yields after some algebra the following expression for the function εD εD(x)=C0+ ∞∑ m=−∞ ( BDm+S0A D m+ ∑ n 6=0 CnG D n,m−n+ ∑ n 6=0 SnH D n,m−n ) eimαx (2.15) Thus, for the bottom wavy wall we have ED0 =C0+B D 0 +S0A D 0 + ∑ n 6=0 Cn(G D n,n) ∗+ ∑ n 6=0 Sn(H D n,n) ∗ (2.16) EDm =B D m+S0A D m+ ∑ n 6=0 CnG D n,m−n+ ∑ n 6=0 SnH D n,m−n m 6=0 Analogous relations can be derived for the upper wall of the channel. Boundary conditions (2.6) are enforced by setting EDm =0 E G 0 =wG E G m =0 m=±1,±2, . . . (2.17) In practice, infiniteFourier expansions have to be truncated to a finite number of terms. Assume that, the summation in (2.8) is restricted to the Fourier modes with |n| ¬N, while conditions (2.17) are imposed for |m| ¬M. It is convenient to introduce 4N +2 real unknowns: the real coefficients C0 and S0 and the real and imaginary parts of the remaining complex coefficients CRm =Re(Cm) C I m = Im(Cm) SRm =Re(Sm) S I m = Im(Sm) (2.18) Instability of viscous incomressible flow... 667 for n= 1,2, . . . ,N. After separation of the real and imaginary parts of equ- ations (2.17) for m> 0, the following set of 4M+2 linear algebraic equation is obtained (m=1,2, . . . ,M) 1 2 C0+ 1 2 A D,G 0 S0+ N∑ n=1 [Re(GD,Gn,m)C R n +Im(G D,G n,m)C I n]+ + N∑ n=1 [Re(HD,Gn,m )S R n +Im(H D,G n,m )S I n] =− 1 2 B D,G 0 Re(AD,Gm )S0+ N∑ n=1 [Re(G D,G n,m+n+G D,G n,m−n)C R n +Im(G D,G n,m+n−G D,G n,m−n)C I n]+ + N∑ n=1 [Re(H D,G n,m+n+H D,G n,m−n)S R n +Im(H D,G n,m+n−H D,G n,m−n)S I n] =−Re(B D,G m ) (2.19) Im(AD,Gm )S0+ N∑ n=1 [Im(G D,G n,m+n+G D,G n,m−n)C R n −Re(G D,G n,m+n−G D,G n,m−n)C I n]+ + N∑ n=1 [Im(H D,G n,m+n+H D,G n,m−n)S R n −Re(H D,G n,m+n−H D,G n,m−n)S I n] =−Im(B D,G m ) where the upper indices D and G correspond to the bottom and upper wall, respectively. After system (2.19) is solved, the velocity field can be evaluated using formulae (2.7) and (2.11). Theproperties of theabovemethodwereanalyzed indetails bySzumbarski (2007). It has been established that themethodperformsbetterwhen M >N and overdetermined linear system (2.19) is solved in the sense of the least squares (Golub and van Loan, 1996). If this system is written shortly as Pc= r (2.20) where c= [C0,C R 1 ,C I 1, . . . ,C R M,C I M,S0,S R 1 ,S I 1, . . . ,S R M,S I M] ⊤, then the least squares solution to (2.20) can be found by solving the following algebraic problem P ⊤ Pc=P⊤r (2.21) with the symmetric and positive definite matrix Z = P⊤P. Linear system (2.21) can be solved bymeans of the QR-decomposition of thematrix P, i.e. P=QR, where Q denotes the orthogonalmatrix (i.e. Q⊤Q= I) and R is the square upper-triangularmatrix. Using the orthogonality of thematrix Q, one gets from (2.21) the linear system Rc=Q⊤r (2.22) 668 J. Szumbarski which can be solved by the elementary back-substitution technique. The QR- decomposition of the matrix P can be computed effectively by the modified Gram-Schmidt or Householder methods (Golub and van Loan, 1996). 2.3. Numerical tests We will present some results of the numerical tests of the method described above. We will be mostly interested in the accuracy of the enforcement of boundary conditions (2.6). For symmetry reasons, the norm of the boundary error can be defined as follows ε= sup x∈[0,λX] |W(x,yD(x))|= sup x∈[0,λX] |W(x,yG(x))| (2.23) Fig. 2. The norm of the boundary error ε plotted versus the number N of Fourier modes used for aproximation of the basic (undisturbed) flow, computed for different values of the amplitude S; the wave number is α=1 In all tests presented here it has been set M =4N. In Fig.2, the norm ε is plotted as a function of the number N of the Fourier modes, computed for different values of the corrugation amplitude S and the wave number α=1. It can be seen that the boundary error diminishes at an exponential rate, so the method is spectrally convergent. Obviously, the rate of convergence for larger corrugation amplitudes is slower. Figure 3 shows analogous results, but this time the amplitude S is fixed (S = 0.3) and the parameter is the wave number α. Again, the norm ε decreases asymptotically at the exponential rate, which, however, drops when the wave number α increases. The contour map of the velocity field of the basis flow computed for the amplitude S=0.4 and the wave number α=1 is shown in Fig.4. Themaxi- mum velocity (about 1.5) is achieved at the locations of the largest distance between walls, while its value at the narrowest part of the channel is about 3 times smaller. Thus, we observe strong spanwise modulation of the velocity, especially in the plane of symmetry y=0. Instability of viscous incomressible flow... 669 Fig. 3. The norm of the boundary error ε plotted versus the number N of Fourier modes used for aproximation of the basic (undisturbed) flow, computed for different values of the amplitude number α; the corrugation amplitude is S=0.3 Fig. 4. A contour plot of the streamwise velocity of the basic flow in channel with symmetric transverse wall waviness. The amplitude is S=0.4 and the wave number is α=1. Note that the basic flow is unidirectional, i.e. the velocity components in the planes z= const vanish identically In is important to evaluate the hydraulic resistance of the flow in the wavy channel. Note that the cross-sectional segmentswith the spanwise length equal to one geometric period λX taken from thewavy-wall and the reference channels have the same areas. In addition, the pressure drop driving both the reference and basic flow is the same. Thus, an appropriate way of evaluation of the flow resistance is to compare the volumetric flow rates computed per one geometric period. The results of such calculations, obtained for different values of S and α, are presented in Fig.5. The remarkable feature of the plotted lines is that they all intersect for α ≈ 1.2, where the ratio between the volumetric rates is very close to 1. It is also observed that the long-wave transversal corrugation leads to some reduction of the flow resistance, while for the shorter-wave corrugation the hydraulic losses increasewith both α and the amplitude S. 670 J. Szumbarski Fig. 5. The ratio between the volumetric flow rates in the way and plain channel flow plotted versus the wave number α for different values of the corrugation amplitude S 3. Analysis of linear stability 3.1. Derivation of the stability equations In this Section, we derive equations of the linear stability theory for the flow in the transversely wavy channel defined in Section 2. The velocity of the basic flow described in the previous Section can be written as vB =W(x,y)ez, while the pressure pB is a linear function of the z coordinate. We introduce time-dependent disturbances v′ = v′(t,x,y,z) and p′ = p′(t,x,y,z), i.e. the velocity and pressure of the disturbed flow can be written as v=vB +v ′ p= pB +p ′ (3.1) Themathematical description of the spatio-temporal evolution of flow distur- bances can be obtained by inserting (3.1) to theNavier-Stokes and continuity equations. Since the flow disturbances are assumed small, the nonlinear term in theNavier-Stokes equation can be neglected. Eventually, the following equ- ations are obtained ∂tu+W∂zu=−∂xp ′+ 1 Re ∆u ∂tv+W∂zv=−∂yp ′+ 1 Re ∆v (3.2) ∂tw+W∂zw+u∂xW +v∂yW =−∂zp ′+ 1 Re ∆w ∂xu+∂yv+∂zw=0 Equations (3.2) contain four unknown scalar fields: three Cartesian compo- nents of the disturbance velocity field and the field of pressure disturbances. Instability of viscous incomressible flow... 671 The coefficients of these equations depend on two space coordinates: x (perio- dically) and y. We will consider a special class of solutions to equations (3.2) known as the normalmodes. These solutions are crucial for the determination of the asymptotic behavior of small disturbances in the flow. For the case con- sidered in this study, the velocity and pressure fields of the normal modes are defined by the following expressions [u,v,w](t,x,y,z) = ei(δx+βz−σt) ∞∑ m=−∞ [gmu ,g m v ,g m w ](y)e imαx+C.C. (3.3) p′(t,x,y,z) = ei(δx+βz−σt) ∞∑ m=−∞ qm(y)eimαx+C.C. where the symbol C.C. stands for the complex conjugate terms. In the above, the symbol β denotes the streamwise wave number, δ is the Floquet parameter and the number σ=σR+iσI is the complex frequency of the normal mode. Thus, formulae (3.3) describe disturbances, which are pe- riodic in the z variable with the period λZ =2π/β. The type of dependence with respect to the x variable is related to the Floquet parameter, or – more precisely – to the value of the fractional part of the ratio δ/α. If this number is irrational, then the disturbed flow is quasi-periodic in the x direction. Other- wise, the disturbance field is x-periodic, however the period can be equal to some multiplicity of the geometric period λX = 2π/α. In fact, it is sufficient to consider δ from the range [0,α/2]. The time variation of the normal mode is determined by its complex fre- quency σ. If σI is negative, then the mode is attenuated or stable. If σI is negative, then the mode is amplified or unstable. In σI = 0 we say that the mode is neutrally stable or critical. As a rule, all normal modes are stable if the Reynolds number is sufficiently small. In the linear stability analysis, we usually look for the smallest value of the Reynolds (denoted by ReL) for which at at least one neutrally stable normal mode exists. If Re>ReL, some disturbances will expotentially grow in time no matter how small they are initially. We say that for Re>ReL the basic (or undisturbed) flow becomes unstable with respect to disturbances with an arbitrary small amplitude. The real part of the complex frequency determines kinematic character of the disturbance field. If σR is different from zero, the disturbances have a form of the traveling wave (the speed of this wave in the streamwise direction is equal σR/β). Such disturbances are also referred to as the oscillatory ones because at any fixedpoint in space one observes a time-periodicmodulation of thedisturbanceamplitude (superimposedon theexponential decay orgrowth). If σR =0, then time variation of the amplitude of disturbances at any space location ismonotonic (nonoscillatory); such disturbances are sometimes called 672 J. Szumbarski stationary. Substitution of expressions (3.3) into equations (3.2) leads to a countable set of ordinary differential equations for the amplitude functions gmu , g m v , g m w and qm, (m=0,±1,±2, . . .).Themathematical descriptionof thedisturbance dynamics can be simplified in a similar manner as in the case of parallel flows (see Schmid andHenningson (2001) formore details) by eliminating the pressure disturbance field and introducing the y-component of vorticity with the formula η= ∞∑ m=−∞ θm(y)ei(tmx+βz−σt)+C.C. (3.4) where θm = i(βgmu − tmg m w ) m=0,±1,±2, . . . (3.5) and tm = δ+mα. Using relations (3.5) together with the formulae implied by the continuity equations itmg m u +∂yg m v +iβg m w =0 m=0,±1,±2, . . . (3.6) one can express the amplitude functions of the velocity components u and w by the velocity component v and the vorticity components η, namely gmu = i tm∂yg m v −βθ m) k2m gmw = i β∂yg m v + tmθ m) k2m (3.7) In the above, we have introduced the real numbers k2m = t 2 m+β 2, whichmust be different from zero for all integer indices m. The latter condition is always satisfied if the streamwisewave number β 6=0,which is assumed in this study. The special case β=0 (meaning that the disturbance field does not actually depend on the streamwise coordinate z) requires a separate treatment and it will not be considered here. After some rather laborious algebra, the following set of equations can be derived (m=0,±1,±2, . . .) −iσLmgmv +S mgmv + + ∑ n>0 (H m,n V gm−nv + Ĥ m,n V gm+nv +H m,n θ θm−n+ Ĥ m,n θ θm+n)= 0 (3.8) −iσθm+Qmθm− itmDF 0 Wg m v + + ∑ n>0 (E m,n V gm−nv + Ê m,n V gm+nv +E m,n θ θm−n+ Ê m,n θ θm+n)= 0 Instability of viscous incomressible flow... 673 The operators appearing in the above equations are defined as follows Lm = ∂yy −k 2 m Q m = iβF0W − 1 Re Lm Sm = iβ(F0WL m−D2F0W)− 1 Re (Lm)2 H m,n V = iβ [β2+ tm−ntm+n k2m−n FnW∂yy − (D 2+k2m)F n W + 2nαtm−n k2m−n DFnW∂y ] Ĥ m,n V = iβ [β2+ tm−ntm+n k2m+n (FnW) ∗∂yy − (D 2+k2m)(F n W) ∗+ − 2nαtm+n k2m+n (DFnW) ∗∂y ] H m,n θ =−i 2nαβ2 k2m−n (FnW∂y +DF n W) (3.9) Ĥ m,n θ = i 2nαβ2 k2m+n [(FnW) ∗∂y +(DF n W) ∗] E m,n V = i nα k2m−n (β2+ tmtm−n)F n W∂y − itmDF n W Ê m,n V =−i nα k2m+n (β2+ tmtm+n)(F n W) ∗∂y − itm(DF n W) ∗ E m,n θ = iβ β2+ tmtm−2n k2m−n FnW Ê m,n θ = iβ β2+ tmtm+2n k2m+n (FnW) ∗ The question arises how to formulate appropriate boundary conditions for equations (3.8). Since the physical domain Ω has been extended to the com- putational domain Ωext, these equations are supposed to be satisfied in the interval −1−S < y < 1+S. However, the ”classical” boundary conditions for the amplitude functions cannot be formulated. Instead, we have to derive some conditions for these functions which would be equivalent to the physical constrains imposed on the velocity disturbance field, namely u(t,x,yD(x),z) =u(t,x,yG(x),z) = 0 v(t,x,yD(x),z) = v(t,x,yG(x),z) = 0 (3.10) w(t,x,yD(x),z) =w(t,x,yG(x),z) = 0 Such conditions can be obtained bymeans of the immerse boundary approach developed by the author (Szumbarski, 2002a). The general idea is to derive explicitly formulae for the Fourier coefficients of the x-periodic distributions of the velocity components at both channel walls and then to demand that 674 J. Szumbarski these coefficients identically vanish. Such a formulation is actually equivalent to the following denumerable set of integral conditions ∞∑ m=−∞ λX∫ 0 [gmu ,g m v ,g m w ](t,yD(x))e i(m−k)αx dx=0 (3.11) ∞∑ m=−∞ λX∫ 0 [gmu ,g m v ,g m w ](t,yG(x))e i(m−k)αx dx=0 k=0,±1,±2, . . . With the use of relations (3.7), conditions (3.11) can be expressed in terms of the amplitude functions gmv and θ m (m=0,±1,±2, . . .). 3.2. Spectral discretization of the stability equations and the boundary conditions In order to obtain a numerically tractable problem, we have to discretize dif- ferential equations (3.8) and to derive a finite dimensional approximation of integral conditions (3.11). To this end, we introduce truncated Chebyshev expansions gnv(y)≈ KV∑ k=0 ΓnkTk(y) θ n(y)≈ Kθ∑ k=0 ΘnkTk(y) (3.12) where the basic functions are defined as Tk(y)= tk[y/(1+S)] and tk denotes the ”standard” Chebyshev polynomial of the kth order. For a detailed de- scription of the properties of these polynomials, the reader may refer to Boyd (2001). Next, finite Chebyshev expansions (3.12) are inserted into equations (3.8)1. Then each of these equations is multiplied by the polynomial Tj (j = 0,1, . . . ,KV − 4) and integrated in the interval [−1− S,1+S] with the modified Chebyshev weight function ω(y)= 1+S √ (1+S)2−y2 (3.13) Equivalently, it is demanded that the residua of all 4th order equations (3.8)1 are ω-orthogonal to thefinite dimensional space of thepolynomials of an order not higher than KV −4.As a result, a set of KV −3 linear algebraic equations is obtained for each integer number m. Analogous procedure applied to the 2nd order differential equations (3.8)2 produces Kθ−1 algebraic equations for each integer index m. This way, the number of algebraic equations obtained due to the Chebyshev-Galerkin projection for the mth Fourier mode of the Instability of viscous incomressible flow... 675 disturbance field is equal to KV +Kθ − 4, while the number of unknowns corresponding to this mode is KV +Kθ +2. Thus, for each integer m one is free to add six algebraic equations in order to enforce boundary constrains (3.11). Wewill skipwriting an explicit formof the algebraic equations obtained in effect of the discretization procedure because it has been described in details elsewhere (Szumbarski, 2002a, 2007). For practical calculations, the Fourier representation of the disturbance field has to be truncated to a finite num- ber of modes. To this end, we will assume that all Fourier modes with the index m outside the range {−M,.. . ,M} are omitted. This way we obtain (2MS+1)(KV +Kθ+2) linear algebraic equations with the unknownCheby- shev coefficients {Γmk , k = 0, . . . ,KV , |m| ¬MS} and {Θ m k , k = 0, . . . ,Kθ, |m| ¬MS}. If we define the vector z= [ {Γn0 ,Γ n 1 , . . . ,Γ n KV−1 ,ΓnKV ;Θ n 0 ,Θ n 1 , . . . ,Θ n Kθ−1 ,ΘnKθ},n= (3.14) =−MS, . . . ,0, . . . ,MS ]⊤ then the algebraic system can be conveniently written in amatrix-vector form as follows Pz= iσQz Bz=0 (3.15) Thus, the normal modes are approximated by the eigenvectors of generali- zed eigenvalue problem (3.15) and the complex frequencies are simply the corresponding eigenvalues. In the current study, we are interested in flow de- stabilization, so we will look for the normal mode with the largest imaginary part σI and investigate its parametric variation. 3.3. Numerical analysis of an unstable normal mode Investigation of unstable modes can be a tedious numerical task because of the number of parameters involved (wave numbers, Floquet parameter, amplitude of the wall corrugation and the Reynolds number) and the size of eigenvalue problem (3.15). The exact number of equations in (3.15) cannot be established a priori (the convergence analysis can be done only after an appropriate eigensolution has been identified), but on the basis of the author’s previous experience a size of 2000-3000 can be anticipated. These expectation has been actually confirmed, and some results of the convergence study are presented in the Appendix. The searching for unstablemodes can be perform efficiently done by para- metric continuation of the selected eigensolution of the Orr-Sommerfeld and Squire equations corresponding to the referencePoiseuille flow.Anatural con- tinuation parameter is the amplitude of the wave corrugation S. The numeri- 676 J. Szumbarski cal tools for such continuations are the methods of inverse or subspace itera- tions (Golub and van Loan, 1996; Saad, 1992). The key problem is to identify the eigenmodes of the Poiseuille flow, which are most sensitive to destabili- zation, while the corrugation amplitude increases. Some of these modes are already known: Ehrenstein (1996) showed that a short-wave transversal cor- rugation (with α around 12) can destabilize the least attenuated, originally two-dimensionalOrr-Sommerfeldmode (or theTollmien-Schlichtingwave) and the critical Reynolds number ReL of the flow can be reduced to 2600, appro- ximately. Here, however, we are looking for a possibility of effective flow de- stabilization at theReynolds numberswhich aremuch lower. It turns out that such a possibility appears in the range of parameters that has not apparently been investigated before. In addition, the normal mode being destabilized is quite different than themodes investigated in the previous works. This mode turns out to be the least attenuated Squire mode denoted here as Sq1. Origi- nally (i.e. in the Poiseuille flow), the wave vector of this mode can be written as κ≡ [κx,κy,κz] = [0,0,β]. It can be shown that the corresponding veloci- ty field has only a transversal component (here, the x-component), which is periodic in the streamwise direction z. Moreover, in the Poiseuille flow, the mode Sq1 is damped for all Reynolds numbers. The situation becomes quite different when the transversal waviness is in- troduced. If the corrugation wave number α assumes values in some range around 1.2 then the imaginary part of the complex frequency σI of themode Sq1 grows quickly with the increasing amplitude S, and the mode becomes unstable even though the Reynolds number is rather low. The variations of the imaginary part of the complex frequency of the mode Sq1 with the stre- amwise wave number β computed for the wall waviness with the amplitude S =0.3 and different various wave numbers α are shown in Fig.6. The Rey- nolds number corresponding to these results is only 100. The destabilization effect is most pronounced for the waviness with the wave number around 1.2, i.e. for the corrugation period which is about 3 times larger than the average wall distance (or the channel height) and about 10 times larger than for the corrugation considered in Ehrenstein (1996). The obtained results demonstra- te the possibility of radical reduction of the critical Reynolds number, even much below the value of 100. It is also very important that such a spectacular effect is achieved without introducing any additional hydraulic drag. In fact, the volumetric flow rate of the flow in awavy channel is slightly larger than in the case of the reference Poiseuille flow. Since the driving pressure drop and the area of the cross-section in both cases are the same, we conclude that the most destabilizing wall waviness actually leads also to some reduction of the hydraulic resistance. This feature is quite opposite to what is observed in the case of the longitudinal waviness with a similar amplitude, where significant increase of the hydraulic resistance is inevitable. Instability of viscous incomressible flow... 677 Fig. 6. The imaginary part of the complex frequency of the unstable normal mode Sq1 plotted versus the streamwise wave number β, computed for different values of the wave number α. The amplitude of the wall waviness is S=0.3 and the Reynolds number is Re=100 The destabilization effect of the transversal waviness of the channel walls canbenicely illustratedbyplotting thecurveof neutral stability, i.e. the line in theRe-β plane alongwhich σI =0. Such curves, calculated for three different values of the corrugation amplitude S and for the wave number α = 1, are presented in Fig.7. One can see that the critical Reynolds number ReL for the amplitude S = 0.4 is approximately equal to 58, hence it is smaller by two orders of magnitude (!) than ReL for the reference Poiseuille flow. Fig. 7. Neutral stability curves of the mode Sq1 computed for different values of the corrugation amplitude S. The wave number α=1 It is interesting to investigate the kinematical structure of the unstablemo- de Sq1 is presented in Fig.8 and Fig.9. The geometric parameters of the wall waviness are S=0.4 and α=1, the streamwisewave number β=0.4 and the Reynolds number is Re = 100. The disturbance velocity field has been nor- 678 J. Szumbarski Fig. 8. The disturbance velocity fields of the unstable normal mode Sq1 computed at two channel sections: (a) z=0, (b) z=λ z /4 for the Reynolds number is Re=100, corrugation amplitude S=0.4 and wave number α=1. The streamwise wave number is β=0.4. The contour maps show the streamwise velocity component w (dashed lines correspond to negative values). The velocity field is normalized so that its maximummagnitude in the flow domain is attained in the plane z=0 and equals 1 Instability of viscous incomressible flow... 679 malized in such away that themaximumof the velocitymagnitude inside the channel is attained at the plane z=0 and is equal to 1. The spanwise struc- ture of the velocity field at the plane z=0 is presented in Fig.8a. The upper plot shows a contour map of the streamwise velocity component w, while the bottom one presents the velocity vectors projected on the plane of the section. Figure 8b shows analogous results computed for the plane z =λz/4= π/2β. The structure of the velocity field in the symmetry plane y = 0 is depic- ted in Fig.9. The space-periodic system of counter-rotating vortices is clearly seen. All presented pictures give an idea of complicated three-dimensional ki- nematics of flow disturbances related to the unstable mode Sq1. It should be emphasized that the real part of thismode is not zero, so the flowdisturbances have a form of a traveling wave. The velocity of this wave in the streamwi- se direction is equal to the ratio σR/β, and in the case considered above, it is equal 0.854, which is close to the average velocity of the basic flow in the symmetry plane y=0. Fig. 9. The disturbance velocity fields of the instable mode Sq1 computed in the plane y=0 (the central plane of the wavy channel). All parameters and normalization are the same as in Fig.8 4. Summary and conclusions In this study, the laminar viscous flow in a channel with transversely oriented wall waviness has been determined using a semi-analytical method of immer- sed boundaries which avoids the domain transformation. This method has been proved to be spectrally convergent, and if the geometric period of the wall waviness is large enough, it gives good accuracy even for relatively large corrugation amplitudes. 680 J. Szumbarski The linear stability analysis of the basic flow has been performed and the possibility of flowdestabilization at theReynoldsnumberas lowas 60hasbeen established. The wave number of the most destabilizing wall corrugations is about 1.2, which corresponds to the geometric period being about 3 times larger than the average channel height. The unstable normal mode originates from the fundamental Squiremode of the reference Poiseuille flow. Originally, this mode is attenuated and has a simple kinematical structure (its velocity has only a transversal component). The presence of wall waviness strongly destabilizes this mode, and the disturbance velocity field contains a system of counter-rotating vortices. The unstable disturbance field has a form of the wave which travels downstream with a velocity close to the average velocity of the undisturbed flow in the symmetry plane y=0 of the channel. It can be expected that for the Reynolds numbers slightly higher than the critical value ReL, the nonlinear saturation process of the unstable modewill occur.Thenewnonlinear state of flow, emergingdue to the saturation process, wouldhave a complexvortical structure, goodmixingproperties andhydraulic resistance comparable with the reference Poiseuille flow. One of possible ways to investigate the saturation process is to derive the low dimensional dynami- cal system out of the full Navier-Stokes equation using the Galerkin method and the normalmodes of stability equations derived in the role of divergence- free basic functions. The obtained system of nonlinear ordinary differential equations can be solved numerically, and bifurcation behavior can be investi- gated. In further perspective, the properties ofmass and heat transport of the saturated flowwill also be investigated. Appendix Herewe present the results of convergence study for the spectral discretization method applied to the stability equations. The test calculations have been carried for thewavenumber α=1, theamplitude of thewallwaviness S=0.4 and the Reynolds number Re = 100. For such parameters, the fundamental Squiremode Sq1 is unstable. InTable 1,we show the real and imaginary parts of the complex frequency σ of this mode as well as the value of the boundary error ε, i.e. the maximum magnitude of the disturbance velocity field at the channelwalls.Thevelocity field of themodehasbeennormalized in suchaway that the maximum velocity magnitude in the flow domain is 1. The number of Chebyshev polynomials used for approximation of the amplitude functions is KV =Kθ =50. Rapid convergence of the complex frequency is observedwith the increasing number of the Fourier modes N. Note that high accuracy of σ is achieved in Instability of viscous incomressible flow... 681 spite of the relatively large boundary error. The latter decreases with N in the manner that proves spectral convergence of the approximation. Table 1 N σR σI ·10 2 ε 8 0.34174012 1.37161397 0.10651771 10 0.34171914 1.37262042 0.05299126 12 0.34171812 1.37276618 0.02218982 14 0.34171807 1.37277909 0.00994565 16 0.34171807 1.37278009 0.00471424 References 1. Blancher S., Creff R., Le Quere P., 1998, Effect of Tollmien-Schlichting wave on convective heat transfer in a wavy channel. Part I: Linear analysis, International Journal of Heat and Fluid Flow, 19, 39-48 2. Blancher S., Creff R., Le Quere P., 2004, Analysis of convective hydro- dynamic instabilities in a symmetric wavy channel, Physics of Fluids, 16, 10, 3726-3737 3. Boyd J.P., 2001, Chebyshev and Fourier Spectral Methods, 2nd Ed. Dover Publication Inc., Minneola NY 4. Cabal A., Szumbarski J., Floryan J.M., 2002, Stability of flow in a wavy channel, Journal of Fluid Mechanics, 457, 191-212 5. Cho K.J., Kim M., Shin H.D., 1998, Linear stability of two-dimensional steady flow in wavy-walled channel,Fluid Dynamic Research, 23, 349-370 6. Ehrenstein U., 1996, On the linear stability of channel flow over riblets, Physics of Fluids, 8, 11, 3194-3196 7. Floryan J.M., Szumbarski J., Wu X., 2002, Stability of flow in a channel with vibrating walls,Physics of Fluids, 14, 11, 3927-3936 8. Floryan J.M., 1997, Stability of wall-bounded shear layers in the presence of simulated distributed surface roughness, J. Fluid Mech., 335, 29-55 9. Floryan J.M., 2005, Two-dimensional instability of flow in a rough channel, Physics of Fluids, 17, 044101-1 (electronic version) 10. Floryan J.M., 2003, Vortex instability in a diverging-converging channel, J. Fluid. Mech., 482, 17-50 11. Ghaddar N., El-Hajj A., 2000, Numerical study of heat transfer augmen- tation of viscous flow in corrugated channels, Heat Transfer Engineering, 21, 35-46 682 J. Szumbarski 12. Goldstein J.L., SparrowE.M., 1977,Heat-mass transfer characteristics for flow in a corrugated wall channel, Journal of Heat Transfer, 99, 187-195 13. Golub G.H., Van Loan Ch.F., 1996, Matrix Computations, 3rd Ed., The John Hopkins University Press, Baltimore 14. Grek G.R., Kozlov V.V., Titarenko S.V., 1996, An experimental study of the influence of riblets on transition, Journal of FluidMechanics, 315, 31-49 15. Luchini P., Trombetta G., 1995, Effects of riblets upon flow stability, Ap- plied Scientific Research, 54, 4, 313-321 16. Luchini P., 1995, Asymptotic analysis of laminar boundary-layer flow over finely grooved surfaces, European Journal of Mechanics B Fluids, 14, 2, 169- 195 17. NicenoB.,NobileE., 2001,Numerical analysis of fluidflowandheat transfer in periodic wavy channel, International Journal of Heat and Fluid Flow, 22, 156-167 18. NishimuraT.,Murakami S., ArakawaS.,KawamuraY., 1990,Flow ob- servations and mass transfer characteristics in symmetrical wavy-walled chan- nels at moderate Reynolds number for steady flow, Int. J. Heat and Mass Transfer, 33, 835-844 19. Nishimura T., Arakawa S., Murakami S., Kawamura Y., 1989, Oscil- latory viscous flow in symmetric wavy-wall channels, Chemical Engineering Science, 44, 2211-2224 20. NishimuraT., KawamuraY., 1995,Three-dimensionality of oscillatory flow in a two-dimensional symmetric sinusoidal wavy-walled channel, Exp. Therm. Fluid Sci., 10, 62-73 21. Nishimura T., Yoshino T., Kawamura Y., 1987, Numerical flow analysis of pulsatile flow in a channel with symmetric wavywalls atmoderateReynolds numbers, J. Chem. Eng. Jpn., 20, 479-485 22. Nishimura T., Kojima N., 1995, Mass transfer enhancement in a sinusoidal wavy-walled channel for pulsatile flow, Int. Journal of Heat andMass Transfer, 38, 1719-1731 23. Saad Y., 1992,Numerical methods for large eigenvalue problems, Manchester University Press 24. SchmidP.J., HenningsonD.S., 2001, Stability and transition in shear flows, Applied Mathematical Sciences, 142, Springer, NewYork 25. SelvarajanS., TulapurkaraE.G.,VasantaRamV., 1999,Stability cha- racteristics of wavy channel flows,Physics of Fluids, 11, 3, 579-589 26. Sobey I.J., 1980, On flow through furrowed channels. Part 1: Calculated flow patterns, Journal of Fluid Mechanics, 96, 1-26 27. Stephanoff K.D., Sobey I.J., Bellhouse B.J., 1980, On flow through furrowed channels. Part 2:Observedflowpatterns, Journal of FluidMechanics, 96, 27-32 Instability of viscous incomressible flow... 683 28. Szumbarski J., Floryan J.M., 1999, A direct spectral method for determi- nation of flows over corrugated boundaries, Journal of Computational Physics, 156, 378-402 29. Szumbarski J., 2002c,On the origin of unstablemodes in viscous channel flow subject to periodically distributed surface suction, Journal of Theoretical and Applied Mechanics, 40, 4, 847-871 30. Szumbarski J., 2002a, Immersedboundary approach to stability equations for a spatially periodic viscous flow,Archives of Mechanics, 54, 3, 199-222 31. Szumbarski J., 2002b,LinearAnalysis of the disturbance dynamics in spatial- ly periodic channel flow. Poster presentation [in Polish],XVPolish Conference of Fluid Mechanics, Augustów 2002 32. Szumbarski J., Floryan J.M., 2006, Transient disturbance growth in a cor- rugated channel, Journal of Fluid Mechanics, 568, 243-272 33. Szumbarski J., 2007, Instability of the viscous liquid flow in the wavy-wall channel [in Polish], Habilitation Dissertation, In: Prace Naukowe Politechniki Warszawskiej, Mechanika, Warszawa 34. Viswanath P.R., 2002,Aircraft viscous drag reduction using riblets,Progress in Aerospace Sciences, 38, 571-600 35. Wang G., Vanka S.P., 1995, Convective heat transfer in periodic wavy pas- sages, Int. Journal of Heat and Mass Transfer, 38, 17, 3219-3230 36. WuX., Luo J., 2006, Influence of small imperfections on the stability of plane Poiseuille flow and the limitation of Squire’s theorem, Physics of Fluids, 18, 044104 (electronic form) Zjawisko destabilizacji przepływu laminarnego w kanale z poprzecznie pofalowanymi ścianami Streszczenie Pracapoświęcona jest zjawiskudestabilizacjiprzepływu laminarnegowkanalewy- wołanej poprzecznym pofalowaniem ścian. Przepływ niezaburzony wyznaczono przy pomocy semi-analitycznejmetody zanurzonych granic.Wyprowadzono równania sta- teczności, które następnie poddano dyskretyzacji przy użyciuwielomianówCzebysze- wa. Pokazano, że poprzeczne pofalowanie ścian o odpowiedniej geometrii prowadzi do destabilizacji przepływu przy niskich liczbach Reynoldsa i omówiono własności pola zaburzeń.Wydaje się, że opisany efektmożemieć znaczenie dla intensyfikacjimiesza- nia i poprawy efektywności procesów transportu w technice cieplnej, biotechnologii i medycynie. Manuscript received December 18, 2006; accepted for print March 12, 2007