Spatial stability of similarity solutions for viscous ﬂows in channels with porous walls

The spatial stability of similarity solutions for an incompressible ﬂuid ﬂowing along a channel with porous walls and driven by constant uniform suction along the walls is analyzed. This work extends the results of Durlofsky and Brady @ Phys. Fluids 27 , 1068 ~ 1984 !# to a wider class of similarity solutions, and examines the spatial stability of small amplitude perturbations of arbitrary shape, generated at the entrance of the channel. It is found that antisymmetric perturbations are the best candidates to destabilize the solutions. Temporally stable asymmetric solutions with ﬂow reversal presented by Zaturska, Drazin, and Banks @ Fluid Dyn. Res. 4 , 151 ~ 1988 !# are found to be spatially unstable. The perturbed similarity solutions are also compared with fully bidimensional ones obtained with a ﬁnite difference code. The results conﬁrm the importance of similarity solutions and the validity of the stability analysis in a region whose distance to the center of the channel is more than three times the channel half-width. © 2000 American Institute of Physics. @ S1070-6631 ~ 00 ! 02003-1 #


I. INTRODUCTION
The Navier-Stokes equation for an incompressible viscous flow along a channel with porous walls, driven by constant uniform suction, admits a similarity solution. This solution was studied in 1953 by Berman, 1 who reduced the bidimensional Navier-Stokes equation to a fourth order nonlinear ordinary differential equation with two boundary conditions at each wall. This equation depends on a sole nondimensional parameter, the transversal Reynolds number, R, defined in terms of the channel width and the suction velocity. Similarity solutions were tested in early experiments on laminar flows in ducts with wall injection, where evidence of good agreement between measurements and theoretical predictions was reported. 2 Many authors have studied Berman's equation ͑see, for example, Refs. 3-12͒ and found a very rich structure of solutions. Zaturska et al., 8 however, proved that most of these solutions are temporally unstable to linear perturbations. On the other hand, the spatial stability of the solutions has not been completely analyzed. In 1984 Durlofsky and Brady 13 studied the spatial stability of symmetric solutions under linear symmetric perturbations.
This work presents the spatial stability analysis for both symmetric and asymmetric solutions under arbitrary small amplitude perturbations. The basic equations of the problem are shown in the next section and the stability analysis is developed in Sec. III. A comparison between the spatially perturbed similarity solutions and fully bidimensional ones is presented in Sec. IV. Section V is devoted to the conclusions.

II. BASIC EQUATIONS
The coordinate system of Fig. 1 is considered, with the origin in the center of the channel and the x-axis parallel to the walls. The parameters of the problem are the suction velocity V, the channel width 2h, and the ͑constant͒ density and viscosity of the fluid, respectively. Steady state solutions for the nondimensional velocity and pressure fields u គ (x,y), p(x,y) satisfy the continuity and Navier-Stokes equations. Thus, div͑u គ ͒ϭ0, ͑1͒ u គ Ϫgrad͑u គ ͒ϭϪgrad͑ p ͒ϩdiv͑ grad͑u គ ͒͒/R, with u គ ϭϮŷ on the walls. The unit vectors in the x and y directions are denoted by x and ŷ , respectively. The transversal Reynolds number is RϭVh/. The variables are expressed in units of the length h, the velocity V and the pressure V 2 .
A similarity solution of the form, ͑2͒ p͑x,y ͒ϭ͑ y ͒ϩAx 2 /2, leads to the following equation for f (y): where Ј denotes differentiation along the transversal coordinate. The boundary conditions for f (y) are Once Eq. ͑3͒ is solved, the pressure field corresponding to Eq. ͑2͒ can be obtained from the following expressions: where two stable asymmetric solutions ͑types I a and I a Ј) are generated in a pitchfork bifurcation. These solutions remain stable below R Ia ϭ12.963, where a Hopf bifurcation takes place. A saddle-node bifurcation at R II ϭ12.165 gives birth to two temporally unstable symmetric solutions ͑types II and III͒ one of which undergoes a pitchfork bifurcation at R III ϭ15.415, generating two new asymmetric branches ͑types III a and III a Ј). As R→ϱ, solutions of type I and II differ by exponentially small terms only, 7 and they coalesce in a cusp type bifurcation.
Solutions of type II with RϽR II Јϭ13.119 and those of type III, III a , and III a Ј have regions of flow reversal in the center of the channel. Solutions of type I a and I a Ј have a region of flow reversal near one wall if RϾR Ia Ј ϭ6.557.

III. SPATIAL STABILITY ANALYSIS
The spatial stability analysis of Eq. ͑3͒ was performed taking solutions of the form f(x,y)ϭx f (y)ϩx H (y) for the stream function. The term x H (y) is regarded as a perturbation to the similarity solution. The linearized approximation for H (y), valid in the region not close to the center of the channel, leads to the following characteristic value problem: with homogeneous boundary conditions,

͑5͒
For a given Reynolds number, the existence of characteristic values with Re͑͒Ͻ1, implies the instability of the similarity solution. Under these conditions certain perturbations at the mouth of the channel can grow downstream, destroying the similarity solution in the interior.
This kind of analysis ͑for symmetric solutions͒ has been performed by Durlofsky and Brady. 13 These authors, however, considered boundary conditions of the form, and thus restricted the analysis to symmetric perturbations.

A. Symmetric solutions
For symmetric solutions of type I with boundary conditions given by Eq. ͑6͒, Durlofsky and Brady found that Ͼ1 for all eigenvalues and 0ϽRϽϱ. Consequently, these solutions are stable under symmetric perturbations. However, this stability is marginal when R→ϱ since the lowest eigenvalue tends to one. When general conditions ͑5͒ are chosen, antisymmetric perturbations are included in the analysis. In Fig. 3, eigenvalues corresponding to symmetric ͑dotted lines͒ and antisymmetric ͑full lines͒ eigenfunctions of Eq. ͑4͒ are plotted vs the Reynolds number. These results were obtained with a code based on spectral methods following the procedure used by Orszag 15 to analyze the temporal stability of  Poiseuille flow. Half of the branches are those described in Ref. 13. The others are due to the antisymmetric perturbations. The lowest branch, and thus the one that determines the instability of the solution, corresponds to an antisymmetric eigenfunction. This branch crosses the line ϭ1 at R ϭR I , which means that these solutions lose their temporal and spatial stability at the same value of the Reynolds number. This fact becomes apparent when the eigenvalue equation arising from the temporal stability analysis is considered ͑see Ref. 8͒. If temporal perturbations of the form G(y)exp(st) are proposed, the following characteristic value problem for the eigenfunctions G and the eigenvalues s is obtained,

͑7͒
For this kind of perturbations, there is an exchange of stability if sϭ0. A comparison between Eqs. ͑4͒ and ͑7͒ shows that, at a given R, a temporal mode with sϭ0 is also a spatial marginal mode ͑ϭ1͒.
Consequently, the spatial stability of the symmetric solutions of type I is determined by the antisymmetric eigenvalues, and the range of Reynolds numbers for which these solutions are spatially stable coincides with the range of temporal stability.
The real part of the eigenvalues for symmetric ͑dotted lines͒ and antisymmetric ͑full lines͒ eigenfunctions of the other two symmetric solutions ͑types II and III͒ are plotted in Figs. 4͑a͒ and 4͑b͒. Figure 4͑a͒ shows the eigenvalues with minimum magnitude of Re͑͒ corresponding to symmetric and antisymmetric eigenfunctions of solutions of types II and III. For type II solutions there is a branch of real eigenvalues with 0Ͻр1, so that these solutions are unstable. Nevertheless, for large R, the instability of solutions of type II with respect to symmetric perturbations becomes marginal. In fact, all branches of type II eigenvalues coalesce with the corresponding ones of type I as R→ϱ. This behavior, which is linked to the presence of the cusp bifurcation at Rϭϱ mentioned in Sec. II, can be observed comparing eigenvalues at Rϭ20 in Figs. 3 and 4͑a͒.
As in the case of type I solutions, the antisymmetric perturbations prove to be more dangerous for stability since the minimum magnitude eigenvalues are negative in this case. It is interesting to note in this figure that the structure of the spectrum is somewhat different from the one shown in Ref. 13 since complex eigenvalues, not considered in that paper, are included in Fig. 4 by plotting their real parts. This is exhibited in Fig. 4͑a͒ as the merging of two real branches at the lower right hand corner, Rϳ19.5, and at Rϳ12.2. Figure 4͑b͒ presents branches with negative eigenvalues of larger magnitude. Note that the scale of the vertical axis is different from that of Fig. 4͑a͒. The branches in the upper right hand corner of the figure correspond to type III solutions. They merge with type II branches with flow reversal at the minimum value of RϭR II ϭ12.165. These branches of type II solutions tend to Ϫϱ as the Reynolds number approaches the value R II Јϭ13.119, which is the largest value of R for solutions of type II with flow reversal. The uppermost branch in this figure was already calculated in Ref. 13. The overall structure of this part of the spectrum for type II solutions is quite different from that of Fig. 4͑a͒.

B. Asymmetric solutions
When the symmetric solution of type I loses its temporal stability, two asymmetric solutions which are mirror images of each other in the center line of the channel ͑types I a and I a Ј) are generated in a pitchfork bifurcation. The temporal stability analysis for these solutions shows that they are stable for RϽR Ia ϭ12.963. 8 The asymmetric solutions are characterized by a displacement of the stagnation point towards one of the walls. The axial velocity increases near this wall and decreases near the other one. This asymmetry effect increases with the Reynolds number. When RϾR Ia Ј ϭ6.557 there is a region close to one of the channel walls where the flow reverses. Figure 5 shows the spectrum for both symmetric ͑type I, It is interesting to note that eigenvalues with a very large negative real part as those presented in Fig. 6 are found not only in solutions of type I a , but also in symmetric solutions of type II near the point where the flow reverses (RϽR II Ј), as was shown in Fig. 4͑b͒. These eigenvalues and their corresponding eigenfunctions were investigated analytically. 16 The analysis showed that the structure of the characteristic value problem depends strongly on the sign of f Љ(y) at the points where f Ј(y) vanishes. In Fig. 7 an eigenfunction of the solution of type I a for Rϭ6.6 is plotted. This value of the Reynolds number is close to the critical value R Ia Ј ϭ6.557, where flow reversal takes place. The modulus of the associated eigenvalue is as large as 10 6 . The inset shows in more detail the structure of the boundary layer in the region of flow reversal close to the wall. Inside the channel the eigenfunction is proportional to f Ј(y), with boundary layers at the walls, a property that can be derived from analytical considerations. 16 The asymmetric solutions of types III a and III a Ј ͑the mirror image of III a ) were also investigated. They were found to be unstable under both symmetric and antisymmetric pertur- bations. The results are shown in Figs. 8͑a͒ and 8͑b͒, where the real part of eigenvalues of type II and III solutions ͓already described in Figs. 4͑a͒ and 4͑b͔͒ are plotted in dotted lines, and those of type III a and III a Ј in full lines. The minimum magnitude eigenvalues are presented in Fig. 8͑a͒. As in Fig. 4, a different scale for the vertical axis is used in Fig.  8͑b͒ to plot the eigenvalues with larger negative real part.

IV. FULLY BIDIMENSIONAL SOLUTIONS
In order to compare the spatially perturbed similarity solutions with fully bidimensional steady state ones, Eqs. ͑1͒ with no similarity assumptions were integrated numerically. A code in finite differences which solves the equations iteratively was developed. Multigrid and local relaxation techniques were employed to improve convergence. The problem has a new nondimensional parameter, L, the channel length given in terms of its half-width.
To complete the characterization of the problem, boundary conditions at both edges of the channel should be imposed. By assuming symmetry with respect to the xϭ0 axis, the domain is reduced to the xϾ0 half of the channel. The stream function and the vorticity are equal to zero at xϭ0. The conditions at the entrance, xϭL, depend on the flow conditions outside the channel.
The vorticity-stream function representation was used and the system was discretized using centered finite differ-ences over a uniform grid with 513 nodes in the x direction and 65 nodes in the y direction.
Boundary conditions of the form f c (y)ϭ f (y) ϩ␣H (y), with ␣Ӷ1, were selected at the entrance of the channel. The functions f (y) and H (y) are solutions of Eqs. ͑3͒ and ͑4͒, respectively.
In order to analyze the difference between a bidimensional solution and the similarity solution, the function Normalized functions ␦(x,y) are presented together with their amplitudes ⌬(x)ϭmax y ␦(x,y), for Lϭ20 and Rϭ3, 5, 7, and 9 in Fig. 9. The conditions imposed at the mouth of the channel correspond to solutions of type I perturbed with their first eigenfunction. The two curves plotted on the y ϭ1 plane are ⌬(x) and a straight line that represents the linear behavior in x of the similarity solution. The surface plots show that for every value of x the perturbation conserves its profile. The amplitudes ⌬(x) decay faster than the solution in the first two cases. This agrees with the fact that those cases correspond to stable solutions ͑Ͼ1͒. On the other hand, the last two graphs show a spatial attenuation that is smaller than the linear one of the similarity solution. This is to be expected since for Rϭ7 and 9 the similarity solutions are unstable.
In Fig. 10 the local attenuation rate of the solution, defined as ⌫(x)ϭd(ln(⌬)/d(ln(x)) ͑full lines͒, and the characteristic values ͑dotted lines͒ provided by the linear analysis are plotted. An excellent agreement is observed when the distance to the origin is more than three times the half-width of the channel. These results, commented in detail in Ref. 16, confirm the validity of the stability analysis.

V. CONCLUSIONS
Solutions of type I and I a are the only branches of solutions of Berman's equation which are spatially stable for a certain range of Reynolds numbers. For solutions of type I the range of spatial stability is the same as the range of temporal stability. The instability threshold is determined by an antisymmetric perturbation. On the other hand, the range of spatial stability of the type I a solutions is considerably shorter than the range of temporal stability. For this type of solutions the spatial stability is restricted to those solutions which have no flow reversal.
A singularity in the spectrum of solutions with flow reversal is present. Branches of eigenvalues with negative real part diverge as the Reynolds number approaches the critical value where solutions with flow reversal first appear.
Two-dimensional simulations of the flow in the channel were developed. These simulations were carried out with the aid of a finite difference code written in a vorticity-stream function formulation of the problem. The results given by this code support the validity of the linear stability analysis in the region not close to the center of the channel.
In conclusion, we have shown that the potential presence of antisymmetric perturbations considerably restricts the range of Reynolds numbers where similarity solutions are stable. There are stable similarity solutions for RϽ6.557 only. These solutions are symmetric if RϽ6.001 and asymmetric if RϾ6.001. However, it is important to keep in mind that in practical cases some unstable modes may not develop due to the finite length of the channel.

ACKNOWLEDGMENTS
This work was supported in part by Contract No. CTI*CT 91-0944 of the European Union, by the University of Buenos Aires ͑Grants Nos. EX088 and TX32͒, and by CONICET ͑Consejo Nacional de Investigaciones Cientificas y Técnicas͒.