- Email: [email protected]
Hydrodynamic instability of fiber suspensions in channel flows
Fluid Dynamics Research 34 (2004) 251 – 271
Hydrodynamic instability of ber suspensions in channel "ows You Zhenjianga , Lin Jianzhonga;∗ , Yu Zhaoshengb a
State Key Laboratory of Fluid Power Transmission and Control, Department of Mechanics, Zhejiang University, Yugu Road, Hangzhou 310027, PR China b Department of Mechanical and Mechatronic Engineering, University of Sydney, NSW 2006, Australia Received 24 May 2003; received in revised form 10 October 2003; accepted 8 January 2004 Communicated by Z.-S. She
Abstract A linear stability analysis of channel "ow in the presence of ber suspensions is performed. Based on slender-body theory, a modi ed stability equation is derived by using the natural closure approximation to determine the ber orientation. It is found that the "ow instability of ber suspensions is governed by two parameters: H , a ratio between the axial stretching resistance of bers and the inertial force of the "uid, and the orientational di;usivity coe
1. Introduction The behavior of bers in a "ow a;ects the rheology and light scattering properties of suspensions that are of interest in many areas of industry. In the use of ber additives as drag-reducing agents, ∗
Corresponding author. Tel.: +86-571-87952882; fax: +86-571-87951464. E-mail address: [email protected] (L. Jianzhong).
c 2004 Published by The Japan Society of Fluid Mechanics and Elsevier B.V. 0169-5983/$30.00 All rights reserved. doi:10.1016/j."uiddyn.2004.01.002
252
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
for example, Arranaga (1970) gave the friction reduction characteristics of brous and colloidal substances. Mansour (1985) presented an explicit approximation to the relationship that describes velocity pro les and drag reduction in turbulent ber suspension "ow. Nsom (1994) presented the results of a theoretical study of dilute and semi-concentrated suspensions of sti; bers, in which the drag reduction is directly computed with respect to the cylinder radii ratio, the ber concentration and density. In addition, it is reported that a minimum of highly conducting bers enhances the thermal and electrical conductivity of composites (Mackaplow et al., 1994; Sundararajakumar and Koch, 1999). When the suspended particles are slender, their orientations strongly a;ect the rheological properties of the "ow. In fact, the "ow can induce a preferred orientation of the bers. The rheological properties of ber suspensions depend on the volume fraction and the orientations of the bers. Je;ery (1922) was among the rst to investigate the motion of a single rigid ellipsoidal particle immersed in a viscous Newtonian liquid. Folgar and Tucker (1984) developed an evolution equation for concentrated ber suspensions, where the ber– ber interactions are taken into account by adding a di;usion term of Je;ery’s equation. Rahnama et al. (1995) considered long-range hydrodynamic interactions between bers, including an anisotropic rotary di;usivity, which is consistent with experimental observations in semi-dilute suspensions (Stover et al., 1992). A series of tensors are used by Advani and Tucker (1987) to describe ber orientation. Although these tensors o;er a concise representation of the ber orientation, and are convenient for numerical simulation, some additional information is required to completely determine the "ow-induced orientation. A variety of closure approximations have been proposed to obtain a closed set of equations for orientation tensors (Lipscomb et al., 1988; Advani and Tucker, 1990; Verleye and Dupret, 1993; Cintra and Tucker, 1995). Batchelor (1971) applied slender-body theory and developed a constitutive model to determine the stress of large-aspect-ratio bers in dilute and semi-dilute suspensions. Later Batchelor’s model was adopted to study the rheological properties of ber suspensions in consideration of inter- ber interactions (Shaqfeh and Fredrickson, 1990; Mackaplow and Shaqfeh, 1996). However, only limited attention has been devoted to the instability of ber suspensions. It is very important to understand the e;ect of bers on the "ow instability. Pilipenko et al. (1981) discussed the linear instability of Couette "ow of ber suspensions between coaxial cylinders. Galdi and Reddy (1999) studied the well-posedness of the equations governing the "ow of ber suspensions, through the use of second- and fourth-order orientation tensors to account for the presence of bers. They demonstrated that the linear closure relation leads to anomalous behavior, in that the rest state of the "uid is unstable for certain ranges of the ber particle number. Azaiez (2000a,b) presented the results of a linear instability analysis of a mixing layer at high Reynolds numbers, and included the e;ects of the presence of rigid bers on the temporal instability of the "ow. Gupta et al. (2002) analyzed the linear stability of the Taylor–Couette "ow of semi-dilute non-Brownian suspension. They found that the ber additives suppress the centrifugal instability, which is attributed to the fact that the suspension develops negative rst and second normal stresses in the TC "ow. Munganga and Reddy (2002) established the existence and uniqueness of solutions to the governing equations for ber suspensions, and the existence of a unique classical solution, local in time, is proved for the cases of both linear and quadratic closure rules. Linear stability analysis has been widely applied to the hydrodynamic instability problems of single- or multi-phase "ow. Yamamoto et al. (1995, 1998) studied the stability of the "ow in a helical tube of small curvature for a wide range of the torsion parameter experimentally and theoretically.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
253
Both results showed that the torsion had a destabilizing e;ect on the "ow. Wie and Malik (1998) investigated the e;ect of surface waviness on boundary layer transition in two-dimensional "ow and concluded that the wall waviness destabilized disturbance waves. McAlpine and Drazin (1998) presented a linear stability analysis of two-dimensional perturbations of the steady "ow driven between inclined plane walls by a line source at the intersection of the walls. Vladimir et al. (2001) discussed the e;ects of wall compliance on the hydrodynamic stability and transition delay of a wall-jet "ow. Yanase et al. (2002) studied the laminar "ow through a curved rectangular duct over a wide range of the aspect ratio. Five branches of steady solutions were found and linear stability characteristics were investigated for all these solutions. Hirayama and Takaki (1995) analyzed a two-dimensional uniform "uidized bed. They found that the inertia term due to the average "uid velocity was responsible for the instability, while the particle di;usion and the e;ective particle viscosity suppressed the growth of disturbances. Hirayama and Takaki (2000) gave a linear stability analysis of one-dimensional "uidized bed containing two kinds of solid particles. This system in some cases was found less unstable than a "uidized bed with one kind of particle. While the instability of Newtonian channel "ow, a class of wall-bounded "ow, has been broadly studied (Pekeris and Shkoller, 1967; Orszag, 1971; Itoh, 1974), but the instability of ber suspensions in channel "ows is still to be explored. Hence, in the present paper, we examine the e;ects of bers on the "ow instability in channel "ows through a parametric analysis and interpret the mechanisms behind them. 2. Mathematical model A Schematic of the channel "ow is shown in Fig. 1, where U0 and V0 represent the streamwise and transverse velocities, respectively. The distance between two plates is d which is taken to be much larger than the dimension of bers. The background "ow is assumed to be laminar and is driven by streamwise pressure gradient. The maximum of U0 is U on the x-axis. The Reynolds number is de ned as Re = U (d=2)= , where is the "uid density and the dynamic viscosity of the "uid. The basic equations that describe the "ow are the continuity and momentum equations ∇ · u = 0;
(2.1)
9u 1 + u · ∇u = − ∇p + ∇ · ; 9t
(2.2)
Fig. 1. Schematic diagram of the channel "ow.
254
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
where p is the isotropic pressure, and the stress tensor. The boundary conditions are u |y=±1 = 0 :
(2.3)
In the current work, the "ow of the ber suspensions is considered a continuum, and the apparent stress is composed of contributions from both the Newtonian suspending "uid and the anisotropic particles of bers. Hence, the stress tensor can be divided into two terms = s + f :
(2.4)
The rst term corresponds to the contribution of the suspending "uid 1 ”; Re
(2.5)
” = ∇u + ∇uT :
(2.6)
s =
Here ” is the rate of strain tensor. The second term on the right side of (2.4) f represents the contribution of the bers. Once this stress is expressed in terms of the "ow variables, the "ow behavior of the whole system is then determined.
3. Model of ber stress In order to determine the ber stress f , we need solve for the orientation distribution of the bers in the "ow. The bers are regarded as rigid and axisymmetric cylinders. Assume the spatial distribution of bers to be uniform. A unit vector p, as shown in Fig. 2, is used to denote the orientation of the ber. The orientation state of bers at a point in space can be described by a probability distribution function, (p). This function is de ned so that the probability of nding a ber between p and (p + dp) is given by (p) dp. If we assume that bers move with the bulk motion of the "uid then (p) may be regarded as a convected quantity. The continuity condition is then D + ∇ · (p˙ ) = 0: Dt
Fig. 2. Orientation of a single ber described by a unit vector p.
(3.1)
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
For the motion of a single ber, we will use the equation (Cintra and Tucker, 1995) 1 1 p˙ = − ! · p + ( · p − : ppp) − Dr · ∇ ; 2 2 ! = ∇uT − ∇u
255
(3.2) (3.3)
where ! is the vorticity tensor, r is the aspect ratio of bers, and = (r 2 − 1)=(r 2 + 1). Dr is the rotary di;usivity resulting from inter- ber hydrodynamic interactions. Folgar and Tucker (1984) suggested a simple empirical expression where Dr is assumed to be isotropic and can be replaced by CI |”|I , where ” is de ned in (2.6). The coe
r2 : 3ln( 2 = )
(3.9)
256
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
For brevity, we de ne H = A =Re and F = a4 : ”. The second term on the right-hand side of Eq. (3.8) is absorbed into the pressure term because of its isotropy. Then the ber stress is written as f = HF:
(3.10)
Substituting Eqs. (3.10) and (2.4)–(2.6) into (2.2), we have the momentum equation 1 2 9u + u · ∇u = − ∇p + ∇ u + ∇ · (HF): 9t Re
(3.11)
4. Derivation and solution of stability equation We examine the linear instability of the "ow by disturbing the base "ow in nitesimally, and limit the analysis to the case of a two-dimensional "ow. For the channel "ow, we take the base "ow as U0 (y) = 1 − y2 ; V0 (y) = 0; !0 (y) = 2y; y3 ; (4.1) 3 where U0 (y); V0 (y) are the base-state streamwise and transverse velocities, respectively, !0 (y) the spanwise vorticity and 0 (y) the corresponding streamfunction expressed in dimensionless form. The components of the base state tensors a02 (y) and F0 (y) are then derived (Appendix A). The streamfunction (y), vorticity !(y), velocity, ber orientation tensor a2 (y) and F(y) are represented by the base-state pro le plus a small perturbation: 0 (y)
=y−
(x; y; t) =
0 (y)
+ ’ (y)ei"(x−ct) ;
u(x; y; t) = U0 (y) + u (y)ei"(x−ct) ; v(x; y; t) = V0 (y) + v (y)ei"(x−ct) ; !(x; y; t) = !0 (y) + ! (y)ei"(x−ct) ; a2 (x; y; t) = a02 (y) + a2 (y)ei"(x−ct) ; F(x; y; t) = F0 (y) + F (y)ei"(x−ct) ;
(4.2)
We now examine the temporal instability, and consider normal modes of the perturbation, where the wavenumber " is treated as real and c = cr + ici is complex. The growth rate of the disturbance is given by &i = "ci , while cr is the propagation velocity of disturbance. Substituting formula (4.2) into Eq. (3.11), we have the linear stability equation governing the channel "ow of ber suspensions (Appendix B) 1 2 i"[(U0 − c)(D2 − "2 ) − D2 U0 ]’ − (D − "2 )2 ’ Re = H [(D2 + "2 )F12 + i"D(F11 − F12 )]
(4.3)
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
257
with the boundary conditions ’ (±1) = 0;
D’ (±1) = 0;
(4.4)
where D = d=dy and the prime denotes derivation of the base-state velocity with respect to y. One should note that the special case of H = 0 represents conditions of Newtonian "ow, and Eq. (4.3) degenerates into standard Orr–Sommerfeld (O–S) equation. The components of the tensor F in the modi ed O–S equation (4.3) can be expressed as the functions of ’ (y), and (4.3) is rewritten as follows: 4
Ji Di ’ (y) = 0;
(4.5)
i=0
where coe
258
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 3. The neutral instability curve of Newtonian channel "ow.——: present results, - - - - -: Pekeris and Shkoller’s results. Nishioka, Iida and Ichikawa’s experimental results: : damped; : nearly neutral; ♦: ampli ed.
Fig. 4. Neutral instability curves of channel "ow.——–: H = 0; –·–·–: H = 0:0005, CI = 0:005; - - - -: H = 0:0005, CI = 0:01.
1=r 2 6 6 1=r, and the Reynolds number we set is the order of 104 . It is seen that if r increases to 103 , H has the order of magnitude 10−5 –10−3 . The higher r rises, the larger the range of H is. Without loss of generality, we shall restrict the investigation to conditions for which H is in the order of 10−4 . Fig. 4 shows the changes of the neutral instability curve considering the e;ect of the ber parameter H and the interaction coe
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
259
Fig. 5. Growth rate of small disturbance as a function of the wavenumber (Re = 10 000): (a) with CI = 0:01, ——: H = 0, - - - - -: H = 0:0001,–· –· –: H = 0:0003, ––––: H = 0:0005; (b) with H = 0:0001, ———: CI = 0, - - - - -: CI = 0:005: –· –· –: CI = 0:01, ––––: CI = 0:02.
Reynolds number and a reduction of the unstable region of small disturbances, therefore reinforcing the "ow stability. It is evident that the larger the Reynolds number, the clearer the stabilization provided by bers. Additionally, the results are similar for di;erent CI , keeping H constant. The variation of the growth rate of small disturbances versus the wavenumber for di;erent values of H and CI is illustrated in Fig. 5. As H is increased, it is clear that the "ow instability is altered. The range of unstable disturbance waves is diminished from that of the Newtonian "ow, and the largest growth rate that governs the instability of the "ow is observably reduced. Moreover, the entire unstable spectrum is shifted towards shorter waves. An increase of CI attenuates the instability similarly to that of H but the e;ect is not so remarkable. Fig. 6(a) and (b) depicts the distribution of the streamwise and transverse velocity disturbances vertical to the wall, respectively. To present clearly the behavior of variables close to the walls, only one half of the channel width is shown here. It can be seen that in the case of the Newtonian "ow, the contribution of the velocity disturbance is mainly concentrated in the near-wall region. Fiber additives cause the maximum of the velocity to reduce and move away from the walls. It is noticeable that inter- ber hydrodynamic interactions reduce the peak value of the disturbance, and hence attenuate the instability of the "ow. In Fig. 7, each component of the base-state ber orientation tensor as a function of CI is illustrated. The base-state orientation tensor a02 is the stable solution for the undisturbed "ow. It is worthwhile to note that a02 has no concern with the spatial coordinates or ber parameter H , as it is only a function of the interaction coe
260
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 6. Variation of the velocity disturbances: (a) the streamwise velocity disturbances; (b) the transverse velocity disturbances. —–: H = 0; –· –· –: H = 0:0005, CI = 0:01; - - - - - : H = 0:0005, CI = 0:02.
Fig. 7. Variation of the three components of the ber orientation tensor of the base-state "ow.
Fig. 8 represents the variation of the rst component of the orientation tensor disturbance a11 , from which the orientation distribution of bers in the "ow is derived. The value a11 re"ects the disturbance’s e;ect on the degree of ber alignment with the "ow direction. If the interaction coe
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
261
Fig. 8. Variation of the orientation of bers: (a) the component of the disturbance of the rst-order orientation tensor; (b) the decline angle of ber. ——: H = 0:0001, CI = 0:01: –· –· –: H = 0:0005, CI = 0:01; - - - - - : H = 0:0001, CI = 0:005.
In what follows, the contours of di;erent terms appearing in the vorticity equation (5.1) are presented. The results are given for one wavelength in the streamwise direction, 0 ¡ x ¡ 2 =", but for only a small fraction of the extent in the direction normal to the wall, 0:7 ¡ y ¡ 1. Examinations of the full range [ − 1; 1] in the transverse direction show very minute variations of all the terms outside the near-wall region. In order to present clearly the features of variables in the near-wall domain, only a fraction of the range is considered in the analysis below. For comparison in the following gures, the two graphs side by side are plotted with the same range and interval of contour levels (the wall is on the top boundary). The solid lines stand for positive values and the dashed lines stand for negative values of a function: 2 9 9 1 2 9 92 92 + U0 ! = −v D!0 + ∇ ! +H (F − F11 ) : (5.1) F12 + − 9t 9x Re 9x2 9y2 9x9y 22 Fig. 9 shows the rate of production of vorticity disturbance for ber suspensions. There are three distinctions between Newtonian "uid and ber suspensions. First, in the region close to the wall, the disturbance magnitude of the vorticity production rate is greater in the Newtonian "uid, which indicates that the disturbance energy obtained from the wall action is weakened because of bers. Second, there is a streamwise phase shift out of the near-wall region in the ber suspension "ow, leading to a better alignment of the positive contours (or negative contours) in the two regions of the "ow. Third, the gradient of contours in the vertical direction is smaller in the ber suspension "ow, as compared to that in the Newtonian "ow. This is another proof of the e;ect of bers on restraining the growth of the vorticity production rate disturbance near the wall. Fig. 10 represents the convective term (−v D!0 ) for ber suspensions. Because the derivative of base-state vorticity is a constant, this term varies in the same way as the transverse velocity disturbance. As a positive contribution to the vorticity production disturbance, the convective term is slightly reduced near the wall owing to the presence of bers in the "ow. There is also a streamwise phase shift in the ber suspension "ow, compared with Newtonian "ow. The region of
262
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 9. Contours of the rate of production of vorticity disturbance at the maximum growth rate: " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01), contour level range: min = −0:1, max = 0:1, interval: 0.008.
Fig. 10. Contours of the convective term at the maximum growth rate: " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01), contour level range: min = −0:15, max = 0:15, interval: 0.005.
the positive contours (or negative contours) is nearly superposed to the lower part of the vorticity production rate contours, indicating that the convective term is the main source for the vorticity disturbance. Fig. 11 represents the viscous term for ber suspensions. Since this term is in direct proportion to the second-order derivative of the vorticity disturbance, the smaller variation gradient of the vorticity disturbance near the wall leads to the smaller magnitude of the viscous term, compared to the Newtonian "ow. There is a major di;erence in the distribution of the positive and negative contours, contrasted with the lower part of contours of vorticity production rate. This result suggests that the viscous term depletes the vorticity disturbance. Relatively, the area of negative values in ber suspensions is larger than that in the Newtonian "ow. The contributions of shear and normal stress are identically zero in the Newtonian "uid, so the two cases of ber suspensions with di;erent H are presented for comparison. In Fig. 12, the augment of the e;ect of shear stress term is obvious with increasing H . Similar to the viscous term, the contribution of the shear stress term is negative to the vorticity disturbance. Fig. 13 illustrates the contribution of the rst normal stress di;erence in the two types of ber suspension "ow. The distinct e;ect of bers on the normal stress term can be found. Since the extents
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
263
Fig. 11. Contours of the viscous term at the maximum growth rate: " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01), contour level range: min = −0:05, max = 0:05, interval: 0.004.
Fig. 12. Contours of the shear stress term at the maximum growth rate: (a) " = 0:96, &i = 0:0036 (H = 0:0001, CI = 0:01), (b) " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01). Contour level range: min = −0:02, max = 0:02, interval: 0.002.
Fig. 13. Contours of the normal stress term at the maximum growth rate: (a) " = 0:96, &i = 0:0036 (H = 0:0001, CI = 0:01), (b) " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01). Contour level range: min = −0:002, max = 0:002, interval: 5e − 5.
of positive and negative contours are close in the region where contours of vorticity production rate are positive (or negative), the rst normal stress di;erence is deduced to contribute negligibly to the vorticity disturbance, as compared to the other terms.
264
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Table 1 Relative values of di;erent terms in the equation for the rate of change of enstrophy H
CI
"
&i
E1
E2 E1
E3 E1
E4 E1
E5 E1
0 0.0001 0.0005 0.0005
0 0.01 0.01 0.02
0.95 0.96 1.00 1.01
0.0042 0.0036 0.0012 0.0003
1:47e − 4 1:14e − 4 6:72e − 5 1:80e − 5
1.14 1.27 2.80 4.67
−0:14 −0:25 −1:19 −2:25
0 −0:02 −0:71 −2:03
0 2:8e − 4 0.10 0.61
The analysis of the energy equation for disturbances can reveal the mechanism behind the enhancement or attenuation of the "ow instability. As mentioned before, the vorticity has great e;ect on the "ow instability and attracts much attention. So in this paper we adopt the energy analysis based on the enstrophy, which is de ned as the integral of the square of the vorticity (Azaiez, 2000a,b). The enstrophy is considered the basic energy quantity for wall-bounded shear "ows. To obtain the enstrophy budget, take the scalar product of the equation for vorticity disturbance (5.1) with the vorticity disturbance ! , then integrate over the wavelength of the linear disturbance in the streamwise direction and the interval (0:7; 1) in the transverse direction. The resulting disturbance enstrophy budget is given by E1 = E2 + E 3 + E 4 + E 5 ; where
2 =" 1 2 =" 1 1 d 2 E1 = ! dy d x; E2 = −v ! D!0 dy d x; 2 dt 0 0:7 0 0:7 2 =" 1 2 =" 1 2 9 1 92 2 E3 = ! ∇ ! dy d x; E4 = H ! − 2 F12 dy d x; 2 Re 0 9x 9y 0:7 0 0:7 2 =" 1 92 ! − F11 ) dy d x: E5 = H (F22 9x9y 0 0:7
(5.2)
(5.3)
It is signi cant to present the relative contribution of each term on the right-hand side of Eq. (5.2) to the accumulation rate of enstrophy. In Table 1, the results of several types of ber suspensions with di;erent H and CI are listed. The case of the Newtonian "ow is also given for comparison. After the parameters H and CI are known, the most unstable disturbance can be found, the maximum growth rate &i , the corresponding wavenumber " and the value of each term of Eq. (5.2) are then obtained. The Reynolds number is xed to 104 . It is obvious that both the growth rate of disturbance &i and the accumulation rate of enstrophy E1 are reduced as H or CI is increased. For the Newtonian "ow, the accumulation rate of enstrophy grows mainly out of E2 , the transport of base-state vorticity by the vertical velocity disturbance. The viscous dissipation of enstrophy E3 has a negative e;ect on the production of enstrophy, but the e;ect appears very small. For the "ow of ber suspensions, E3 is still negative but the proportion of the relative contribution to the enstrophy production becomes larger as H or CI is increased. E4 and E5 , which are special for the "ow of ber suspensions, represent the production of enstrophy by the shear and normal stress disturbances respectively.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
265
E4 always has a negative e;ect and is much larger than E5 , which always has a positive e;ect and plays a destabilizing role. Physical di;erences in the stability behavior between the wall-bounded shear "ows and the free shear "ows may be found in the analysis of vorticity and energy production of disturbances. For free shear "ows such as the mixing layer, the e;ect of Reynolds number on the vorticity production is negligible and the corresponding viscous term in the vorticity equation can be ignored. While in the wall-bounded channel "ow, the viscous e;ect plays an important role in suppressing the vorticity production and "ow instability. In the mixing layer, the energy sink term E4 is dominant. The negative rate of work by shear stresses can nearly balance the positive contributions from E5 , the normal stress disturbance and E2 , the transfer of energy through Reynolds stress. But in the channel "ow, the e;ect of ber shear stress disturbance is not so prominent, even less than the contribution of viscous term E3 . These two terms resist jointly the energy source term E2 . In order to con rm the stability analysis of ber suspensions, we also design a water cycle system, made of smooth Plexiglas, to conduct the experiment. The "ow is maintained by a constant head tank, which is provided with a continuously adjustable over"ow device. The discharged water "ows through a recti er and enters a horizontal channel, whose width, depth and length are, respectively, 10.5, 1.9 and 150 cm. At the end of the channel is a valve, adjusting the "ux of "ow to the reservoir below. A screw pump lifts the water at the low level to the head tank. A type of synthetic nylon is selected for this experiment, because its density is near that of water (1:04 × 103 kg=m3 ), it absorbs water to only a slight degree (7% at 90% R.H.) and it has a circular cross-section (20 m in diameter). The bers are cut short to keep the range of aspect ratio 20–100. Electron microscope examination of samples in the semi-dilute regime indicates uniform cross-section, smooth surface and nearly straight axis of single ber. It may be considered a proof of the rationality of slender-body theory applied to the stability analysis of ber suspensions in this case. Qualitative and quantitative measurements of disturbances in channel suspensions are carried out using dye emission techniques. To avoid the inlet section, the dye injector nozzle, 1:5 mm in diameter, is located on the centerline and 47 cm from the beginning of the channel. A valve is set to adjust the velocity of dye emission close to the main "ow. During the course of the experiments, measurements in pure water "ow are performed rst, followed by ber suspension experiments. Each experiment starts by adjusting the main "ow valve to set the "ow rate, then adjusting the dye "ow valve carefully to match the injection velocity with the passing "uid. When the initial traces of the dye keep clear and steady, the properties of natural disturbance waves can then be recorded by camera. To perform the ber suspension experiments under the same conditions, a small quantity of bers are added into the reservoir with gentle stirring. After several "ow loops, the suspension becomes homogeneous and the photorecording can be performed. The value of parameter H is obtained based on Re, concentration and aspect ratio of bers. The experimental results include two groups of data; each one comprises a type of water "ow and two types of ber suspensions. Fig. 14 shows the typical photorecord of each case. For each type of "ow, the number of records in each sample is N (¿ 10). Every sample of photo is scanned and magni ed in order to mark the starting point of any observed sinusoidal disturbances and the point of loss of stability. The position of the average starting location is de ned as P1 , and the average ending location of sinusoidal disturbance from the beginning of the channel is de ned as P2 . The stable length is L = P1 − P2 . Results of measurements are listed in Table 2.
266
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 14. Typical dye traces of each type of "ow.
The results show that the stable length of Newtonian "ow is shorter than that of ber suspensions, keeping Re constant. It indicates that ber additives have a stabilizing e;ect on the "ow, consistent with the conclusion of theoretical analysis. For ber suspensions, L increases with an augment of H . Speci cally the enhancement in H leads to signi cant stabilization, which is much clearer for higher values of Re. This property veri es the analysis of neutral instability curves of ber suspensions. We also nd the drag reduction with ber suspensions in the transition stage. The amount of drag reduction depends on parameter H . The theoretical results are qualitatively consistent with the experiments. The mechanism of drag reduction by ber additives is revealed through the variation of velocity pro le of suspension and the decrease of wall shear stress.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
267
Table 2 Results of stable length of streak line Case
Re
H
N
P1 (cm)
P2 (cm)
L (cm)
(a) (b) (c) (d) (e) (f)
1039 1039 1039 1667 1667 1667
0 1:2 × 10−4 3:8 × 10−4 0 0:8 × 10−4 2:5 × 10−4
15 12 13 15 14 23
48.94 48.91 48.94 48.95 49.14 49.03
65.13 65.51 69.45 58.73 64.36 74.33
16.19 16.60 20.51 9.78 15.22 25.30
6. Summary A complete analysis for the linear instability of two-dimensional channel "ow with ber suspensions is conducted in this paper. The planar orientation state of bers is determined using the natural closure approximation and the model of the ber contribution to the total stress is established based on the slender-body theory. The linear stability equation governing the channel "ow of ber suspensions is derived and solved with an e
268
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
ber suspensions. The results of experiments tell the stabilizing e;ect of ber additives on the "ow and the clearer e;ect with larger Re, which are in agreement with the conclusions of the theoretical analysis. Acknowledgements The nancial support for this work from the National Natural Science Foundation of China for Distinguished Scholars through Grant No. 19925210 is gratefully acknowledged. The authors would like to thank Dr. Ken Kiger at University of Maryland for making revisions of our manuscript. Appendix A. Derivation of base-state tensors for the undisturbed (ow In this appendix, we derive the expression of the components of the base-state tensors a02 and F0 in the case of undisturbed ber suspensions. The ber is considered as a slender body (r 1 ; → 1). For the steady two-dimensional channel "ow of ber suspensions, the equation for the orientation tensor a2 is turned into u0 · ∇a02 − ∇u0T · a02 − a02 · ∇u0 = −0 : a04 + 2CI |”0 |(I − 2a02 );
(A.1)
where u02 = (U0 V0 ) and ”0 = ∇u0 + ∇u0T . Substituting Eq. (3.7) into (A.1), we get 1 1−1 CI (41 − 2) + 13 CI ; a02 = (A.2) 1 3 CI (41 − 2) + 1 1 CI where 1 = CI 4CI2 + 1 − 2CI2 . As a consequence, the tensor F0 re"ecting the base-state tensor is obtained (a011 a022 + a012 2 ) DU0 2a011 a012 DU0 ; (A.3) F0 = (a011 a022 + a012 2 ) DU0 2a012 a022 DU0 where D = d=dy. Appendix B. Derivation of coe*cients of modied O–S equation The momentum equation of ber suspensions (3.11) is expanded using the expressions (4.2), with nonlinear terms of the disturbance ignored: 1 3 9p i"(U0 − c)D’ − i"’ DU0 = − + (D ’ − "2 D’ ) + H (i"F11 + DF12 ); 9x Re "2 (U0 − c)’ = −
1 9p + (i"3 ’ − i"D2 ’ ) + H (i"F12 + DF22 ): 9y Re
Eliminating the pressure p from (B.1), we obtain the linear stability equation (4.3).
(B.1)
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
269
Substituting the closure formulation (3.7) into the orientation tensor equation (3.6), we have the orientation equations in the component form: A11 a11 − A12 a12 = (A2 D2 + i"A1 D + "2 A0 )’ ; B12 a12 − B11 a11 = (B2 D2 + i"B1 D + "2 B0 )’ ; a11 + a22 = 0;
(B.2)
where A11 = B12 = i"(U0 − c) + (4CI + 2a012 )DU0 ; A12 = −B11 = 2(1 − a011 )DU0 ; A2 = 0;
A1 = 2(a012 2 − a011 2 ) + 2a011 ;
B2 = 0;
B1 = 2a012 (1 − 2a011 );
A0 = −2a012 ;
B0 = 2a011 − 1:
(B.3)
The solution of Eqs. (B.2) gives a11 = (i"E1 D + "2 E0 )’ ; a12 = (i"G1 D + "2 G0 )’ ; a22 = −(i"E1 D + "2 E0 )’ ;
(B.4)
where A12 B1 + B12 A1 ; det A11 B1 + B11 A1 G1 = ; det det = A11 B12 − B11 A12 : E1 =
A12 B0 + B12 A0 ; det A11 B0 + B11 A0 G0 = ; det
E0 =
(B.5)
Since the orientation of bers is determined, the tensor F can be obtained and the two terms we need are as follows: F11 − F12 = (M2 D2 + i"M1 D + "2 M0 )’ ; F12 = (N2 D2 + i"N1 D + "2 N0 )’ ;
(B.6)
where M2 = KL;
M1 = 4S + 2(KG1 + LE1 )DU0 ;
N2 = −S;
N1 = KL + (LG1 − KE1 )DU0 ;
K = 2a011 − 1;
L = 2a012 ;
M0 = KL + 2(KG0 + LE0 )DU0 ; N0 = −S + (LG0 − KE0 )DU0 ;
S = a2011 − a2012 :
Substituting the expressions of (B.6) into Eq. (4.3), we obtain the O–S equation (4.5). The coe
(B.7) nal form of the modi ed
270
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
J2 = −i"(U0 − c) −
2"2 + H (D2 N2 + i"D(2N1 + M2 ) + "2 (N2 + N0 − M1 )); Re
J3 = H (2DN2 + i"(M2 + N1 )); J4 =
1 + HN2 : Re
(B.8)
References Advani, S.G., Tucker, C.L., 1987. The use of tensors to describe and predict ber orientation in short ber composites. J. Rheol. 31, 751–784. Advani, S.G., Tucker, C.L., 1990. Closure approximations for three-dimensional structure tensors. J. Rheol. 34, 367–386. Arranaga, A.B., 1970. Friction reduction characteristics of brous and colloidal substances. Nature 225, 447–449. Azaiez, J., 2000a. Linear stability of free shear "ows of ber suspensions. J. Fluid Mech. 404, 179–209. Azaiez, J., 2000b. Reduction of free shear "ows instability: e;ects of polymer versus ber additives. J. Non-Newtonian Fluid Mech. 91, 233–254. Batchelor, G.K., 1971. The stress generated in a non-dilute suspension of elongated particles by pure straining motion. J. Fluid Mech. 46, 813–829. Cintra, J.S., Tucker, C.L., 1995. Orthotropic closure approximations for "ow-induced ber orientation. J. Rheol. 39, 1095–1122. Folgar, F.P., Tucker, C.L., 1984. Orientation behaviour of bers in concentrated suspensions. J. Reinf. Plast. Comp. 3, 98–119. Galdi, G.P., Reddy, B.D., 1999. Well-posedness of the problem of ber suspension "ows. J. Non-Newtonian Fluid Mech. 83, 205–230. Gupta, V.K., Sureshkumar, R., Khomami, B., Azaiez, J., 2002. Centrifugal instability of semidilute non-Brownian ber suspensions. Phys. Fluids 14, 1958–1971. Hirayama, O., Takaki, R., 1995. Linear stability of a two-dimensional uniform "uidized bed. Fluid Dyn. Res. 15, 425–435. Hirayama, O., Takaki, R., 2000. Linear stability of one-dimensional "uidized bed containing two kinds of solid particles. Fluid Dyn. Res. 27, 199–216. Itoh, N., 1974. A power series method for the numerical treatment of the Orr–Sommerfeld equation. Trans. Jpn. Soc. Aero Space Sci. 17, 65–75. Je;ery, G.B., 1922. The motion of ellipsoidal particles immersed in a viscous "uid. Proc. R. Soc. London A 102, 161–179. Koch, D.L., 1995. A model for orientational di;usion in ber suspensions. Phys. Fluids 7, 2086–2088. Lipscomb, G.G., Denn, M.M., Hur, D.U., Boger, D.V., 1988. The "ow of ber suspensions in complex geometries. J. Non-Newtonian Fluid Mech. 26, 297–325. Mackaplow, M.B., Shaqfeh, E.S.G., 1996. A numerical study of the rheological properties of suspensions of rigid, non-Brownian bers. J. Fluid Mech. 329, 155–186. Mackaplow, M.B., Shaqfeh, E.S.G., Schiek, R.L., 1994. A numerical study of heat and mass transport in ber suspensions. Proc. R. Soc. London A 447, 77–110. Mansour, A.R., 1985. Simple and explicit equations for the friction factor in turbulent ber suspension "ow. J. Non-Newtonian Fluid Mech. 17, 245–247. McAlpine, A., Drazin, P.G., 1998. On the spatio-temporal development of small perturbations of Je;ery–Hamel "ows. Fluid Dyn. Res. 22, 123–138. Munganga, J.M.W., Reddy, B.D., 2002. Local and global existence of solutions to equations for "ows of ber suspensions. Math. Models Methods Appl. Sci. 12 (8), 1177–1203. Nishioka, M., Iida, S., Ichikawa, Y., 1975. An experimental investigation of the stability of plane Poiseuille "ow. J. Fluid Mech. 72, 731–751. Nsom, B., 1994. Computation of drag reduction in ber suspensions. Fluid Dyn. Res. 14, 275–288.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
271
Orszag, S.A., 1971. Accurate solution of the Orr–Sommerfeld stability equation. J. Fluid Mech. 50, 689–703. Osborne, M.R., 1967. Numerical methods for hydrodynamic stability problems. SIAM J. Appl. Math. 15, 539–557. Pekeris, C.L., Shkoller, B., 1967. Stability of plane Poiseuille "ow to periodic disturbances of nite amplitude in the vicinity of the neutral curve. J. Fluid Mech. 29, 31–38. Pilipenko, V.N., Kalinichenko, N.M., Lemak, A.S., 1981. Stability of the "ow of a ber suspension in the gap between coaxial cylinders. Sov. Phys. Dokl. 26, 646–648. Rahnama, M., Koch, D.L., Shaqfeh, E.S.G., 1995. The e;ect of hydrodynamic interactions on the orientation distribution in a ber suspension subject to simple shear "ow. Phys. Fluids 7, 487–506. Shaqfeh, E.S.G., Fredrickson, G.H., 1990. The hydrodynamic stress in a suspension of rods. Phys. Fluids 2, 7–24. Stover, C.A., Koch, D.L., Cohen, C., 1992. Observations of ber orientation in simple shear "ow of a semi-dilute suspension. J. Fluid Mech. 238, 277–296. Sundararajakumar, R.R., Koch, D.L., 1999. Electrical conductivity of isotropic ber suspensions. Proc. R. Soc. London A 455, 1923–1930. Thomas, L.H., 1953. The stability of plane Poiseuille "ow. Phys. Rev. 91, 780–783. Verleye, V., Dupret, F., 1993. Prediction of ber orientation in complex injection molded parts. ASME Appl. Mech. Division 175, 139–163. Vladimir, L., David, E.L., Jacob, C., 2001. E;ects of wall compliance on the hydrodynamic stability and transition delay of a wall-jet "ow. Fluid Dyn. Res. 29, 115–134. Wie, Y., Malik, M., 1998. E;ect of surface waviness on boundary-layer transition in two-dimensional "ow. Comput. Fluids 27, 157–181. Yamamoto, K., Akita, T., Ikeuchi, H., Kita, Y., 1995. Experimental study of the "ow in a helical circular tube. Fluid Dyn. Res. 16, 237–249. Yamamoto, K., Yanase, S., Jiang, R., 1998. Stability of the "ow in a helical tube. Fluid Dyn. Res. 22, 153–170. Yanase, S., Kaga, Y., Daikai, R., 2002. Laminar "ows through a curved rectangular duct over a wide range of the aspect ratio. Fluid Dyn. Res. 31, 151–183. You, Z.J., Wang, C.B., 2002. Analysis of the numerical methods for the linear stability of channel "ow. J. Zhejiang Inst. Sci. Technol. 19, 207–211.
Hydrodynamic instability of ber suspensions in channel "ows You Zhenjianga , Lin Jianzhonga;∗ , Yu Zhaoshengb a
State Key Laboratory of Fluid Power Transmission and Control, Department of Mechanics, Zhejiang University, Yugu Road, Hangzhou 310027, PR China b Department of Mechanical and Mechatronic Engineering, University of Sydney, NSW 2006, Australia Received 24 May 2003; received in revised form 10 October 2003; accepted 8 January 2004 Communicated by Z.-S. She
Abstract A linear stability analysis of channel "ow in the presence of ber suspensions is performed. Based on slender-body theory, a modi ed stability equation is derived by using the natural closure approximation to determine the ber orientation. It is found that the "ow instability of ber suspensions is governed by two parameters: H , a ratio between the axial stretching resistance of bers and the inertial force of the "uid, and the orientational di;usivity coe
1. Introduction The behavior of bers in a "ow a;ects the rheology and light scattering properties of suspensions that are of interest in many areas of industry. In the use of ber additives as drag-reducing agents, ∗
Corresponding author. Tel.: +86-571-87952882; fax: +86-571-87951464. E-mail address: [email protected] (L. Jianzhong).
c 2004 Published by The Japan Society of Fluid Mechanics and Elsevier B.V. 0169-5983/$30.00 All rights reserved. doi:10.1016/j."uiddyn.2004.01.002
252
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
for example, Arranaga (1970) gave the friction reduction characteristics of brous and colloidal substances. Mansour (1985) presented an explicit approximation to the relationship that describes velocity pro les and drag reduction in turbulent ber suspension "ow. Nsom (1994) presented the results of a theoretical study of dilute and semi-concentrated suspensions of sti; bers, in which the drag reduction is directly computed with respect to the cylinder radii ratio, the ber concentration and density. In addition, it is reported that a minimum of highly conducting bers enhances the thermal and electrical conductivity of composites (Mackaplow et al., 1994; Sundararajakumar and Koch, 1999). When the suspended particles are slender, their orientations strongly a;ect the rheological properties of the "ow. In fact, the "ow can induce a preferred orientation of the bers. The rheological properties of ber suspensions depend on the volume fraction and the orientations of the bers. Je;ery (1922) was among the rst to investigate the motion of a single rigid ellipsoidal particle immersed in a viscous Newtonian liquid. Folgar and Tucker (1984) developed an evolution equation for concentrated ber suspensions, where the ber– ber interactions are taken into account by adding a di;usion term of Je;ery’s equation. Rahnama et al. (1995) considered long-range hydrodynamic interactions between bers, including an anisotropic rotary di;usivity, which is consistent with experimental observations in semi-dilute suspensions (Stover et al., 1992). A series of tensors are used by Advani and Tucker (1987) to describe ber orientation. Although these tensors o;er a concise representation of the ber orientation, and are convenient for numerical simulation, some additional information is required to completely determine the "ow-induced orientation. A variety of closure approximations have been proposed to obtain a closed set of equations for orientation tensors (Lipscomb et al., 1988; Advani and Tucker, 1990; Verleye and Dupret, 1993; Cintra and Tucker, 1995). Batchelor (1971) applied slender-body theory and developed a constitutive model to determine the stress of large-aspect-ratio bers in dilute and semi-dilute suspensions. Later Batchelor’s model was adopted to study the rheological properties of ber suspensions in consideration of inter- ber interactions (Shaqfeh and Fredrickson, 1990; Mackaplow and Shaqfeh, 1996). However, only limited attention has been devoted to the instability of ber suspensions. It is very important to understand the e;ect of bers on the "ow instability. Pilipenko et al. (1981) discussed the linear instability of Couette "ow of ber suspensions between coaxial cylinders. Galdi and Reddy (1999) studied the well-posedness of the equations governing the "ow of ber suspensions, through the use of second- and fourth-order orientation tensors to account for the presence of bers. They demonstrated that the linear closure relation leads to anomalous behavior, in that the rest state of the "uid is unstable for certain ranges of the ber particle number. Azaiez (2000a,b) presented the results of a linear instability analysis of a mixing layer at high Reynolds numbers, and included the e;ects of the presence of rigid bers on the temporal instability of the "ow. Gupta et al. (2002) analyzed the linear stability of the Taylor–Couette "ow of semi-dilute non-Brownian suspension. They found that the ber additives suppress the centrifugal instability, which is attributed to the fact that the suspension develops negative rst and second normal stresses in the TC "ow. Munganga and Reddy (2002) established the existence and uniqueness of solutions to the governing equations for ber suspensions, and the existence of a unique classical solution, local in time, is proved for the cases of both linear and quadratic closure rules. Linear stability analysis has been widely applied to the hydrodynamic instability problems of single- or multi-phase "ow. Yamamoto et al. (1995, 1998) studied the stability of the "ow in a helical tube of small curvature for a wide range of the torsion parameter experimentally and theoretically.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
253
Both results showed that the torsion had a destabilizing e;ect on the "ow. Wie and Malik (1998) investigated the e;ect of surface waviness on boundary layer transition in two-dimensional "ow and concluded that the wall waviness destabilized disturbance waves. McAlpine and Drazin (1998) presented a linear stability analysis of two-dimensional perturbations of the steady "ow driven between inclined plane walls by a line source at the intersection of the walls. Vladimir et al. (2001) discussed the e;ects of wall compliance on the hydrodynamic stability and transition delay of a wall-jet "ow. Yanase et al. (2002) studied the laminar "ow through a curved rectangular duct over a wide range of the aspect ratio. Five branches of steady solutions were found and linear stability characteristics were investigated for all these solutions. Hirayama and Takaki (1995) analyzed a two-dimensional uniform "uidized bed. They found that the inertia term due to the average "uid velocity was responsible for the instability, while the particle di;usion and the e;ective particle viscosity suppressed the growth of disturbances. Hirayama and Takaki (2000) gave a linear stability analysis of one-dimensional "uidized bed containing two kinds of solid particles. This system in some cases was found less unstable than a "uidized bed with one kind of particle. While the instability of Newtonian channel "ow, a class of wall-bounded "ow, has been broadly studied (Pekeris and Shkoller, 1967; Orszag, 1971; Itoh, 1974), but the instability of ber suspensions in channel "ows is still to be explored. Hence, in the present paper, we examine the e;ects of bers on the "ow instability in channel "ows through a parametric analysis and interpret the mechanisms behind them. 2. Mathematical model A Schematic of the channel "ow is shown in Fig. 1, where U0 and V0 represent the streamwise and transverse velocities, respectively. The distance between two plates is d which is taken to be much larger than the dimension of bers. The background "ow is assumed to be laminar and is driven by streamwise pressure gradient. The maximum of U0 is U on the x-axis. The Reynolds number is de ned as Re = U (d=2)= , where is the "uid density and the dynamic viscosity of the "uid. The basic equations that describe the "ow are the continuity and momentum equations ∇ · u = 0;
(2.1)
9u 1 + u · ∇u = − ∇p + ∇ · ; 9t
(2.2)
Fig. 1. Schematic diagram of the channel "ow.
254
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
where p is the isotropic pressure, and the stress tensor. The boundary conditions are u |y=±1 = 0 :
(2.3)
In the current work, the "ow of the ber suspensions is considered a continuum, and the apparent stress is composed of contributions from both the Newtonian suspending "uid and the anisotropic particles of bers. Hence, the stress tensor can be divided into two terms = s + f :
(2.4)
The rst term corresponds to the contribution of the suspending "uid 1 ”; Re
(2.5)
” = ∇u + ∇uT :
(2.6)
s =
Here ” is the rate of strain tensor. The second term on the right side of (2.4) f represents the contribution of the bers. Once this stress is expressed in terms of the "ow variables, the "ow behavior of the whole system is then determined.
3. Model of ber stress In order to determine the ber stress f , we need solve for the orientation distribution of the bers in the "ow. The bers are regarded as rigid and axisymmetric cylinders. Assume the spatial distribution of bers to be uniform. A unit vector p, as shown in Fig. 2, is used to denote the orientation of the ber. The orientation state of bers at a point in space can be described by a probability distribution function, (p). This function is de ned so that the probability of nding a ber between p and (p + dp) is given by (p) dp. If we assume that bers move with the bulk motion of the "uid then (p) may be regarded as a convected quantity. The continuity condition is then D + ∇ · (p˙ ) = 0: Dt
Fig. 2. Orientation of a single ber described by a unit vector p.
(3.1)
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
For the motion of a single ber, we will use the equation (Cintra and Tucker, 1995) 1 1 p˙ = − ! · p + ( · p − : ppp) − Dr · ∇ ; 2 2 ! = ∇uT − ∇u
255
(3.2) (3.3)
where ! is the vorticity tensor, r is the aspect ratio of bers, and = (r 2 − 1)=(r 2 + 1). Dr is the rotary di;usivity resulting from inter- ber hydrodynamic interactions. Folgar and Tucker (1984) suggested a simple empirical expression where Dr is assumed to be isotropic and can be replaced by CI |”|I , where ” is de ned in (2.6). The coe
r2 : 3ln( 2 = )
(3.9)
256
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
For brevity, we de ne H = A =Re and F = a4 : ”. The second term on the right-hand side of Eq. (3.8) is absorbed into the pressure term because of its isotropy. Then the ber stress is written as f = HF:
(3.10)
Substituting Eqs. (3.10) and (2.4)–(2.6) into (2.2), we have the momentum equation 1 2 9u + u · ∇u = − ∇p + ∇ u + ∇ · (HF): 9t Re
(3.11)
4. Derivation and solution of stability equation We examine the linear instability of the "ow by disturbing the base "ow in nitesimally, and limit the analysis to the case of a two-dimensional "ow. For the channel "ow, we take the base "ow as U0 (y) = 1 − y2 ; V0 (y) = 0; !0 (y) = 2y; y3 ; (4.1) 3 where U0 (y); V0 (y) are the base-state streamwise and transverse velocities, respectively, !0 (y) the spanwise vorticity and 0 (y) the corresponding streamfunction expressed in dimensionless form. The components of the base state tensors a02 (y) and F0 (y) are then derived (Appendix A). The streamfunction (y), vorticity !(y), velocity, ber orientation tensor a2 (y) and F(y) are represented by the base-state pro le plus a small perturbation: 0 (y)
=y−
(x; y; t) =
0 (y)
+ ’ (y)ei"(x−ct) ;
u(x; y; t) = U0 (y) + u (y)ei"(x−ct) ; v(x; y; t) = V0 (y) + v (y)ei"(x−ct) ; !(x; y; t) = !0 (y) + ! (y)ei"(x−ct) ; a2 (x; y; t) = a02 (y) + a2 (y)ei"(x−ct) ; F(x; y; t) = F0 (y) + F (y)ei"(x−ct) ;
(4.2)
We now examine the temporal instability, and consider normal modes of the perturbation, where the wavenumber " is treated as real and c = cr + ici is complex. The growth rate of the disturbance is given by &i = "ci , while cr is the propagation velocity of disturbance. Substituting formula (4.2) into Eq. (3.11), we have the linear stability equation governing the channel "ow of ber suspensions (Appendix B) 1 2 i"[(U0 − c)(D2 − "2 ) − D2 U0 ]’ − (D − "2 )2 ’ Re = H [(D2 + "2 )F12 + i"D(F11 − F12 )]
(4.3)
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
257
with the boundary conditions ’ (±1) = 0;
D’ (±1) = 0;
(4.4)
where D = d=dy and the prime denotes derivation of the base-state velocity with respect to y. One should note that the special case of H = 0 represents conditions of Newtonian "ow, and Eq. (4.3) degenerates into standard Orr–Sommerfeld (O–S) equation. The components of the tensor F in the modi ed O–S equation (4.3) can be expressed as the functions of ’ (y), and (4.3) is rewritten as follows: 4
Ji Di ’ (y) = 0;
(4.5)
i=0
where coe
258
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 3. The neutral instability curve of Newtonian channel "ow.——: present results, - - - - -: Pekeris and Shkoller’s results. Nishioka, Iida and Ichikawa’s experimental results: : damped; : nearly neutral; ♦: ampli ed.
Fig. 4. Neutral instability curves of channel "ow.——–: H = 0; –·–·–: H = 0:0005, CI = 0:005; - - - -: H = 0:0005, CI = 0:01.
1=r 2 6 6 1=r, and the Reynolds number we set is the order of 104 . It is seen that if r increases to 103 , H has the order of magnitude 10−5 –10−3 . The higher r rises, the larger the range of H is. Without loss of generality, we shall restrict the investigation to conditions for which H is in the order of 10−4 . Fig. 4 shows the changes of the neutral instability curve considering the e;ect of the ber parameter H and the interaction coe
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
259
Fig. 5. Growth rate of small disturbance as a function of the wavenumber (Re = 10 000): (a) with CI = 0:01, ——: H = 0, - - - - -: H = 0:0001,–· –· –: H = 0:0003, ––––: H = 0:0005; (b) with H = 0:0001, ———: CI = 0, - - - - -: CI = 0:005: –· –· –: CI = 0:01, ––––: CI = 0:02.
Reynolds number and a reduction of the unstable region of small disturbances, therefore reinforcing the "ow stability. It is evident that the larger the Reynolds number, the clearer the stabilization provided by bers. Additionally, the results are similar for di;erent CI , keeping H constant. The variation of the growth rate of small disturbances versus the wavenumber for di;erent values of H and CI is illustrated in Fig. 5. As H is increased, it is clear that the "ow instability is altered. The range of unstable disturbance waves is diminished from that of the Newtonian "ow, and the largest growth rate that governs the instability of the "ow is observably reduced. Moreover, the entire unstable spectrum is shifted towards shorter waves. An increase of CI attenuates the instability similarly to that of H but the e;ect is not so remarkable. Fig. 6(a) and (b) depicts the distribution of the streamwise and transverse velocity disturbances vertical to the wall, respectively. To present clearly the behavior of variables close to the walls, only one half of the channel width is shown here. It can be seen that in the case of the Newtonian "ow, the contribution of the velocity disturbance is mainly concentrated in the near-wall region. Fiber additives cause the maximum of the velocity to reduce and move away from the walls. It is noticeable that inter- ber hydrodynamic interactions reduce the peak value of the disturbance, and hence attenuate the instability of the "ow. In Fig. 7, each component of the base-state ber orientation tensor as a function of CI is illustrated. The base-state orientation tensor a02 is the stable solution for the undisturbed "ow. It is worthwhile to note that a02 has no concern with the spatial coordinates or ber parameter H , as it is only a function of the interaction coe
260
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 6. Variation of the velocity disturbances: (a) the streamwise velocity disturbances; (b) the transverse velocity disturbances. —–: H = 0; –· –· –: H = 0:0005, CI = 0:01; - - - - - : H = 0:0005, CI = 0:02.
Fig. 7. Variation of the three components of the ber orientation tensor of the base-state "ow.
Fig. 8 represents the variation of the rst component of the orientation tensor disturbance a11 , from which the orientation distribution of bers in the "ow is derived. The value a11 re"ects the disturbance’s e;ect on the degree of ber alignment with the "ow direction. If the interaction coe
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
261
Fig. 8. Variation of the orientation of bers: (a) the component of the disturbance of the rst-order orientation tensor; (b) the decline angle of ber. ——: H = 0:0001, CI = 0:01: –· –· –: H = 0:0005, CI = 0:01; - - - - - : H = 0:0001, CI = 0:005.
In what follows, the contours of di;erent terms appearing in the vorticity equation (5.1) are presented. The results are given for one wavelength in the streamwise direction, 0 ¡ x ¡ 2 =", but for only a small fraction of the extent in the direction normal to the wall, 0:7 ¡ y ¡ 1. Examinations of the full range [ − 1; 1] in the transverse direction show very minute variations of all the terms outside the near-wall region. In order to present clearly the features of variables in the near-wall domain, only a fraction of the range is considered in the analysis below. For comparison in the following gures, the two graphs side by side are plotted with the same range and interval of contour levels (the wall is on the top boundary). The solid lines stand for positive values and the dashed lines stand for negative values of a function: 2 9 9 1 2 9 92 92 + U0 ! = −v D!0 + ∇ ! +H (F − F11 ) : (5.1) F12 + − 9t 9x Re 9x2 9y2 9x9y 22 Fig. 9 shows the rate of production of vorticity disturbance for ber suspensions. There are three distinctions between Newtonian "uid and ber suspensions. First, in the region close to the wall, the disturbance magnitude of the vorticity production rate is greater in the Newtonian "uid, which indicates that the disturbance energy obtained from the wall action is weakened because of bers. Second, there is a streamwise phase shift out of the near-wall region in the ber suspension "ow, leading to a better alignment of the positive contours (or negative contours) in the two regions of the "ow. Third, the gradient of contours in the vertical direction is smaller in the ber suspension "ow, as compared to that in the Newtonian "ow. This is another proof of the e;ect of bers on restraining the growth of the vorticity production rate disturbance near the wall. Fig. 10 represents the convective term (−v D!0 ) for ber suspensions. Because the derivative of base-state vorticity is a constant, this term varies in the same way as the transverse velocity disturbance. As a positive contribution to the vorticity production disturbance, the convective term is slightly reduced near the wall owing to the presence of bers in the "ow. There is also a streamwise phase shift in the ber suspension "ow, compared with Newtonian "ow. The region of
262
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 9. Contours of the rate of production of vorticity disturbance at the maximum growth rate: " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01), contour level range: min = −0:1, max = 0:1, interval: 0.008.
Fig. 10. Contours of the convective term at the maximum growth rate: " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01), contour level range: min = −0:15, max = 0:15, interval: 0.005.
the positive contours (or negative contours) is nearly superposed to the lower part of the vorticity production rate contours, indicating that the convective term is the main source for the vorticity disturbance. Fig. 11 represents the viscous term for ber suspensions. Since this term is in direct proportion to the second-order derivative of the vorticity disturbance, the smaller variation gradient of the vorticity disturbance near the wall leads to the smaller magnitude of the viscous term, compared to the Newtonian "ow. There is a major di;erence in the distribution of the positive and negative contours, contrasted with the lower part of contours of vorticity production rate. This result suggests that the viscous term depletes the vorticity disturbance. Relatively, the area of negative values in ber suspensions is larger than that in the Newtonian "ow. The contributions of shear and normal stress are identically zero in the Newtonian "uid, so the two cases of ber suspensions with di;erent H are presented for comparison. In Fig. 12, the augment of the e;ect of shear stress term is obvious with increasing H . Similar to the viscous term, the contribution of the shear stress term is negative to the vorticity disturbance. Fig. 13 illustrates the contribution of the rst normal stress di;erence in the two types of ber suspension "ow. The distinct e;ect of bers on the normal stress term can be found. Since the extents
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
263
Fig. 11. Contours of the viscous term at the maximum growth rate: " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01), contour level range: min = −0:05, max = 0:05, interval: 0.004.
Fig. 12. Contours of the shear stress term at the maximum growth rate: (a) " = 0:96, &i = 0:0036 (H = 0:0001, CI = 0:01), (b) " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01). Contour level range: min = −0:02, max = 0:02, interval: 0.002.
Fig. 13. Contours of the normal stress term at the maximum growth rate: (a) " = 0:96, &i = 0:0036 (H = 0:0001, CI = 0:01), (b) " = 1:00, &i = 0:0012 (H = 0:0005, CI = 0:01). Contour level range: min = −0:002, max = 0:002, interval: 5e − 5.
of positive and negative contours are close in the region where contours of vorticity production rate are positive (or negative), the rst normal stress di;erence is deduced to contribute negligibly to the vorticity disturbance, as compared to the other terms.
264
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Table 1 Relative values of di;erent terms in the equation for the rate of change of enstrophy H
CI
"
&i
E1
E2 E1
E3 E1
E4 E1
E5 E1
0 0.0001 0.0005 0.0005
0 0.01 0.01 0.02
0.95 0.96 1.00 1.01
0.0042 0.0036 0.0012 0.0003
1:47e − 4 1:14e − 4 6:72e − 5 1:80e − 5
1.14 1.27 2.80 4.67
−0:14 −0:25 −1:19 −2:25
0 −0:02 −0:71 −2:03
0 2:8e − 4 0.10 0.61
The analysis of the energy equation for disturbances can reveal the mechanism behind the enhancement or attenuation of the "ow instability. As mentioned before, the vorticity has great e;ect on the "ow instability and attracts much attention. So in this paper we adopt the energy analysis based on the enstrophy, which is de ned as the integral of the square of the vorticity (Azaiez, 2000a,b). The enstrophy is considered the basic energy quantity for wall-bounded shear "ows. To obtain the enstrophy budget, take the scalar product of the equation for vorticity disturbance (5.1) with the vorticity disturbance ! , then integrate over the wavelength of the linear disturbance in the streamwise direction and the interval (0:7; 1) in the transverse direction. The resulting disturbance enstrophy budget is given by E1 = E2 + E 3 + E 4 + E 5 ; where
2 =" 1 2 =" 1 1 d 2 E1 = ! dy d x; E2 = −v ! D!0 dy d x; 2 dt 0 0:7 0 0:7 2 =" 1 2 =" 1 2 9 1 92 2 E3 = ! ∇ ! dy d x; E4 = H ! − 2 F12 dy d x; 2 Re 0 9x 9y 0:7 0 0:7 2 =" 1 92 ! − F11 ) dy d x: E5 = H (F22 9x9y 0 0:7
(5.2)
(5.3)
It is signi cant to present the relative contribution of each term on the right-hand side of Eq. (5.2) to the accumulation rate of enstrophy. In Table 1, the results of several types of ber suspensions with di;erent H and CI are listed. The case of the Newtonian "ow is also given for comparison. After the parameters H and CI are known, the most unstable disturbance can be found, the maximum growth rate &i , the corresponding wavenumber " and the value of each term of Eq. (5.2) are then obtained. The Reynolds number is xed to 104 . It is obvious that both the growth rate of disturbance &i and the accumulation rate of enstrophy E1 are reduced as H or CI is increased. For the Newtonian "ow, the accumulation rate of enstrophy grows mainly out of E2 , the transport of base-state vorticity by the vertical velocity disturbance. The viscous dissipation of enstrophy E3 has a negative e;ect on the production of enstrophy, but the e;ect appears very small. For the "ow of ber suspensions, E3 is still negative but the proportion of the relative contribution to the enstrophy production becomes larger as H or CI is increased. E4 and E5 , which are special for the "ow of ber suspensions, represent the production of enstrophy by the shear and normal stress disturbances respectively.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
265
E4 always has a negative e;ect and is much larger than E5 , which always has a positive e;ect and plays a destabilizing role. Physical di;erences in the stability behavior between the wall-bounded shear "ows and the free shear "ows may be found in the analysis of vorticity and energy production of disturbances. For free shear "ows such as the mixing layer, the e;ect of Reynolds number on the vorticity production is negligible and the corresponding viscous term in the vorticity equation can be ignored. While in the wall-bounded channel "ow, the viscous e;ect plays an important role in suppressing the vorticity production and "ow instability. In the mixing layer, the energy sink term E4 is dominant. The negative rate of work by shear stresses can nearly balance the positive contributions from E5 , the normal stress disturbance and E2 , the transfer of energy through Reynolds stress. But in the channel "ow, the e;ect of ber shear stress disturbance is not so prominent, even less than the contribution of viscous term E3 . These two terms resist jointly the energy source term E2 . In order to con rm the stability analysis of ber suspensions, we also design a water cycle system, made of smooth Plexiglas, to conduct the experiment. The "ow is maintained by a constant head tank, which is provided with a continuously adjustable over"ow device. The discharged water "ows through a recti er and enters a horizontal channel, whose width, depth and length are, respectively, 10.5, 1.9 and 150 cm. At the end of the channel is a valve, adjusting the "ux of "ow to the reservoir below. A screw pump lifts the water at the low level to the head tank. A type of synthetic nylon is selected for this experiment, because its density is near that of water (1:04 × 103 kg=m3 ), it absorbs water to only a slight degree (7% at 90% R.H.) and it has a circular cross-section (20 m in diameter). The bers are cut short to keep the range of aspect ratio 20–100. Electron microscope examination of samples in the semi-dilute regime indicates uniform cross-section, smooth surface and nearly straight axis of single ber. It may be considered a proof of the rationality of slender-body theory applied to the stability analysis of ber suspensions in this case. Qualitative and quantitative measurements of disturbances in channel suspensions are carried out using dye emission techniques. To avoid the inlet section, the dye injector nozzle, 1:5 mm in diameter, is located on the centerline and 47 cm from the beginning of the channel. A valve is set to adjust the velocity of dye emission close to the main "ow. During the course of the experiments, measurements in pure water "ow are performed rst, followed by ber suspension experiments. Each experiment starts by adjusting the main "ow valve to set the "ow rate, then adjusting the dye "ow valve carefully to match the injection velocity with the passing "uid. When the initial traces of the dye keep clear and steady, the properties of natural disturbance waves can then be recorded by camera. To perform the ber suspension experiments under the same conditions, a small quantity of bers are added into the reservoir with gentle stirring. After several "ow loops, the suspension becomes homogeneous and the photorecording can be performed. The value of parameter H is obtained based on Re, concentration and aspect ratio of bers. The experimental results include two groups of data; each one comprises a type of water "ow and two types of ber suspensions. Fig. 14 shows the typical photorecord of each case. For each type of "ow, the number of records in each sample is N (¿ 10). Every sample of photo is scanned and magni ed in order to mark the starting point of any observed sinusoidal disturbances and the point of loss of stability. The position of the average starting location is de ned as P1 , and the average ending location of sinusoidal disturbance from the beginning of the channel is de ned as P2 . The stable length is L = P1 − P2 . Results of measurements are listed in Table 2.
266
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
Fig. 14. Typical dye traces of each type of "ow.
The results show that the stable length of Newtonian "ow is shorter than that of ber suspensions, keeping Re constant. It indicates that ber additives have a stabilizing e;ect on the "ow, consistent with the conclusion of theoretical analysis. For ber suspensions, L increases with an augment of H . Speci cally the enhancement in H leads to signi cant stabilization, which is much clearer for higher values of Re. This property veri es the analysis of neutral instability curves of ber suspensions. We also nd the drag reduction with ber suspensions in the transition stage. The amount of drag reduction depends on parameter H . The theoretical results are qualitatively consistent with the experiments. The mechanism of drag reduction by ber additives is revealed through the variation of velocity pro le of suspension and the decrease of wall shear stress.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
267
Table 2 Results of stable length of streak line Case
Re
H
N
P1 (cm)
P2 (cm)
L (cm)
(a) (b) (c) (d) (e) (f)
1039 1039 1039 1667 1667 1667
0 1:2 × 10−4 3:8 × 10−4 0 0:8 × 10−4 2:5 × 10−4
15 12 13 15 14 23
48.94 48.91 48.94 48.95 49.14 49.03
65.13 65.51 69.45 58.73 64.36 74.33
16.19 16.60 20.51 9.78 15.22 25.30
6. Summary A complete analysis for the linear instability of two-dimensional channel "ow with ber suspensions is conducted in this paper. The planar orientation state of bers is determined using the natural closure approximation and the model of the ber contribution to the total stress is established based on the slender-body theory. The linear stability equation governing the channel "ow of ber suspensions is derived and solved with an e
268
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
ber suspensions. The results of experiments tell the stabilizing e;ect of ber additives on the "ow and the clearer e;ect with larger Re, which are in agreement with the conclusions of the theoretical analysis. Acknowledgements The nancial support for this work from the National Natural Science Foundation of China for Distinguished Scholars through Grant No. 19925210 is gratefully acknowledged. The authors would like to thank Dr. Ken Kiger at University of Maryland for making revisions of our manuscript. Appendix A. Derivation of base-state tensors for the undisturbed (ow In this appendix, we derive the expression of the components of the base-state tensors a02 and F0 in the case of undisturbed ber suspensions. The ber is considered as a slender body (r 1 ; → 1). For the steady two-dimensional channel "ow of ber suspensions, the equation for the orientation tensor a2 is turned into u0 · ∇a02 − ∇u0T · a02 − a02 · ∇u0 = −0 : a04 + 2CI |”0 |(I − 2a02 );
(A.1)
where u02 = (U0 V0 ) and ”0 = ∇u0 + ∇u0T . Substituting Eq. (3.7) into (A.1), we get 1 1−1 CI (41 − 2) + 13 CI ; a02 = (A.2) 1 3 CI (41 − 2) + 1 1 CI where 1 = CI 4CI2 + 1 − 2CI2 . As a consequence, the tensor F0 re"ecting the base-state tensor is obtained (a011 a022 + a012 2 ) DU0 2a011 a012 DU0 ; (A.3) F0 = (a011 a022 + a012 2 ) DU0 2a012 a022 DU0 where D = d=dy. Appendix B. Derivation of coe*cients of modied O–S equation The momentum equation of ber suspensions (3.11) is expanded using the expressions (4.2), with nonlinear terms of the disturbance ignored: 1 3 9p i"(U0 − c)D’ − i"’ DU0 = − + (D ’ − "2 D’ ) + H (i"F11 + DF12 ); 9x Re "2 (U0 − c)’ = −
1 9p + (i"3 ’ − i"D2 ’ ) + H (i"F12 + DF22 ): 9y Re
Eliminating the pressure p from (B.1), we obtain the linear stability equation (4.3).
(B.1)
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
269
Substituting the closure formulation (3.7) into the orientation tensor equation (3.6), we have the orientation equations in the component form: A11 a11 − A12 a12 = (A2 D2 + i"A1 D + "2 A0 )’ ; B12 a12 − B11 a11 = (B2 D2 + i"B1 D + "2 B0 )’ ; a11 + a22 = 0;
(B.2)
where A11 = B12 = i"(U0 − c) + (4CI + 2a012 )DU0 ; A12 = −B11 = 2(1 − a011 )DU0 ; A2 = 0;
A1 = 2(a012 2 − a011 2 ) + 2a011 ;
B2 = 0;
B1 = 2a012 (1 − 2a011 );
A0 = −2a012 ;
B0 = 2a011 − 1:
(B.3)
The solution of Eqs. (B.2) gives a11 = (i"E1 D + "2 E0 )’ ; a12 = (i"G1 D + "2 G0 )’ ; a22 = −(i"E1 D + "2 E0 )’ ;
(B.4)
where A12 B1 + B12 A1 ; det A11 B1 + B11 A1 G1 = ; det det = A11 B12 − B11 A12 : E1 =
A12 B0 + B12 A0 ; det A11 B0 + B11 A0 G0 = ; det
E0 =
(B.5)
Since the orientation of bers is determined, the tensor F can be obtained and the two terms we need are as follows: F11 − F12 = (M2 D2 + i"M1 D + "2 M0 )’ ; F12 = (N2 D2 + i"N1 D + "2 N0 )’ ;
(B.6)
where M2 = KL;
M1 = 4S + 2(KG1 + LE1 )DU0 ;
N2 = −S;
N1 = KL + (LG1 − KE1 )DU0 ;
K = 2a011 − 1;
L = 2a012 ;
M0 = KL + 2(KG0 + LE0 )DU0 ; N0 = −S + (LG0 − KE0 )DU0 ;
S = a2011 − a2012 :
Substituting the expressions of (B.6) into Eq. (4.3), we obtain the O–S equation (4.5). The coe
(B.7) nal form of the modi ed
270
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
J2 = −i"(U0 − c) −
2"2 + H (D2 N2 + i"D(2N1 + M2 ) + "2 (N2 + N0 − M1 )); Re
J3 = H (2DN2 + i"(M2 + N1 )); J4 =
1 + HN2 : Re
(B.8)
References Advani, S.G., Tucker, C.L., 1987. The use of tensors to describe and predict ber orientation in short ber composites. J. Rheol. 31, 751–784. Advani, S.G., Tucker, C.L., 1990. Closure approximations for three-dimensional structure tensors. J. Rheol. 34, 367–386. Arranaga, A.B., 1970. Friction reduction characteristics of brous and colloidal substances. Nature 225, 447–449. Azaiez, J., 2000a. Linear stability of free shear "ows of ber suspensions. J. Fluid Mech. 404, 179–209. Azaiez, J., 2000b. Reduction of free shear "ows instability: e;ects of polymer versus ber additives. J. Non-Newtonian Fluid Mech. 91, 233–254. Batchelor, G.K., 1971. The stress generated in a non-dilute suspension of elongated particles by pure straining motion. J. Fluid Mech. 46, 813–829. Cintra, J.S., Tucker, C.L., 1995. Orthotropic closure approximations for "ow-induced ber orientation. J. Rheol. 39, 1095–1122. Folgar, F.P., Tucker, C.L., 1984. Orientation behaviour of bers in concentrated suspensions. J. Reinf. Plast. Comp. 3, 98–119. Galdi, G.P., Reddy, B.D., 1999. Well-posedness of the problem of ber suspension "ows. J. Non-Newtonian Fluid Mech. 83, 205–230. Gupta, V.K., Sureshkumar, R., Khomami, B., Azaiez, J., 2002. Centrifugal instability of semidilute non-Brownian ber suspensions. Phys. Fluids 14, 1958–1971. Hirayama, O., Takaki, R., 1995. Linear stability of a two-dimensional uniform "uidized bed. Fluid Dyn. Res. 15, 425–435. Hirayama, O., Takaki, R., 2000. Linear stability of one-dimensional "uidized bed containing two kinds of solid particles. Fluid Dyn. Res. 27, 199–216. Itoh, N., 1974. A power series method for the numerical treatment of the Orr–Sommerfeld equation. Trans. Jpn. Soc. Aero Space Sci. 17, 65–75. Je;ery, G.B., 1922. The motion of ellipsoidal particles immersed in a viscous "uid. Proc. R. Soc. London A 102, 161–179. Koch, D.L., 1995. A model for orientational di;usion in ber suspensions. Phys. Fluids 7, 2086–2088. Lipscomb, G.G., Denn, M.M., Hur, D.U., Boger, D.V., 1988. The "ow of ber suspensions in complex geometries. J. Non-Newtonian Fluid Mech. 26, 297–325. Mackaplow, M.B., Shaqfeh, E.S.G., 1996. A numerical study of the rheological properties of suspensions of rigid, non-Brownian bers. J. Fluid Mech. 329, 155–186. Mackaplow, M.B., Shaqfeh, E.S.G., Schiek, R.L., 1994. A numerical study of heat and mass transport in ber suspensions. Proc. R. Soc. London A 447, 77–110. Mansour, A.R., 1985. Simple and explicit equations for the friction factor in turbulent ber suspension "ow. J. Non-Newtonian Fluid Mech. 17, 245–247. McAlpine, A., Drazin, P.G., 1998. On the spatio-temporal development of small perturbations of Je;ery–Hamel "ows. Fluid Dyn. Res. 22, 123–138. Munganga, J.M.W., Reddy, B.D., 2002. Local and global existence of solutions to equations for "ows of ber suspensions. Math. Models Methods Appl. Sci. 12 (8), 1177–1203. Nishioka, M., Iida, S., Ichikawa, Y., 1975. An experimental investigation of the stability of plane Poiseuille "ow. J. Fluid Mech. 72, 731–751. Nsom, B., 1994. Computation of drag reduction in ber suspensions. Fluid Dyn. Res. 14, 275–288.
Y. Zhenjiang et al. / Fluid Dynamics Research 34 (2004) 251 – 271
271
Orszag, S.A., 1971. Accurate solution of the Orr–Sommerfeld stability equation. J. Fluid Mech. 50, 689–703. Osborne, M.R., 1967. Numerical methods for hydrodynamic stability problems. SIAM J. Appl. Math. 15, 539–557. Pekeris, C.L., Shkoller, B., 1967. Stability of plane Poiseuille "ow to periodic disturbances of nite amplitude in the vicinity of the neutral curve. J. Fluid Mech. 29, 31–38. Pilipenko, V.N., Kalinichenko, N.M., Lemak, A.S., 1981. Stability of the "ow of a ber suspension in the gap between coaxial cylinders. Sov. Phys. Dokl. 26, 646–648. Rahnama, M., Koch, D.L., Shaqfeh, E.S.G., 1995. The e;ect of hydrodynamic interactions on the orientation distribution in a ber suspension subject to simple shear "ow. Phys. Fluids 7, 487–506. Shaqfeh, E.S.G., Fredrickson, G.H., 1990. The hydrodynamic stress in a suspension of rods. Phys. Fluids 2, 7–24. Stover, C.A., Koch, D.L., Cohen, C., 1992. Observations of ber orientation in simple shear "ow of a semi-dilute suspension. J. Fluid Mech. 238, 277–296. Sundararajakumar, R.R., Koch, D.L., 1999. Electrical conductivity of isotropic ber suspensions. Proc. R. Soc. London A 455, 1923–1930. Thomas, L.H., 1953. The stability of plane Poiseuille "ow. Phys. Rev. 91, 780–783. Verleye, V., Dupret, F., 1993. Prediction of ber orientation in complex injection molded parts. ASME Appl. Mech. Division 175, 139–163. Vladimir, L., David, E.L., Jacob, C., 2001. E;ects of wall compliance on the hydrodynamic stability and transition delay of a wall-jet "ow. Fluid Dyn. Res. 29, 115–134. Wie, Y., Malik, M., 1998. E;ect of surface waviness on boundary-layer transition in two-dimensional "ow. Comput. Fluids 27, 157–181. Yamamoto, K., Akita, T., Ikeuchi, H., Kita, Y., 1995. Experimental study of the "ow in a helical circular tube. Fluid Dyn. Res. 16, 237–249. Yamamoto, K., Yanase, S., Jiang, R., 1998. Stability of the "ow in a helical tube. Fluid Dyn. Res. 22, 153–170. Yanase, S., Kaga, Y., Daikai, R., 2002. Laminar "ows through a curved rectangular duct over a wide range of the aspect ratio. Fluid Dyn. Res. 31, 151–183. You, Z.J., Wang, C.B., 2002. Analysis of the numerical methods for the linear stability of channel "ow. J. Zhejiang Inst. Sci. Technol. 19, 207–211.