On the non-linear instability of fiber suspensions in a Poiseuille flow

International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / n l m

On the non-linear instability of fiber suspensions in a Poiseuille flow Zhanhong Wan a , Jianzhong Lin a,b,∗ , Hongbing Xiong a a b

State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, Hangzhou, Zhejiang 310027, PR China China Jiliang University, Hangzhou 310018, PR China

A R T I C L E

I N F O

Article history: Received 29 December 2005 Received in revised form 7 June 2008 Accepted 8 June 2008

Keywords: Non-linear instability Fiber suspensions Poiseuille flow

A B S T R A C T

The non-linear instability characteristics of fiber suspensions in a plane Poiseuille flow are investigated. The evolution equation of the perturbation amplitude analogous to Landau equation is formulated and solved numerically for different fiber parameters. It is found that the equilibrium amplitude decreases with the increase of the fiber aspect ratio and volume fraction, i.e. the addition of the fibers reduces the amplitude of the perturbation, and leads to the reduction of the flow instability. This phenomenon becomes significant for large volume concentration and aspect ratio. The mechanism of the reduction of the flow instability is also analyzed by taking into account of the modification of the mean flow and the energy transfer from the mean flow to the perturbation flow. © 2008 Elsevier Ltd. All rights reserved.

1. Introduction It is reported that fiber suspensions have more and more applications in industry, such as composite manufacturing, paper making, fiber spinning [1,2] and so on. The characterization of the flow or the fiber orientation in a suspension is of particular interest in the industrial processing research. When fibers form a flow-induced orientation state, fiber suspension exhibits anisotropy. Therefore, the addition of fibers to a Newtonian fluid can drastically change the flow kinematics, even at very low concentrations. Besides, the flow can induce a preferred orientation of the fibers which, upon solidification, influences the mechanical properties of the resulting fiberreinforced composite. The composite is stiffer and stronger in the direction of greatest orientation, and weaker and more compliant in the direction of least orientation. Drag reduction in turbulent flow, affected by suspending a very small concentration of fibers, is another application. The mode of flow instability is responsible for the characterization of the flow as well as the spatial and orientational distributions of fibers, therefore, it is significant to understand the effects of fibers on the flow instability. Actually, the hydrodynamic instability in flows has been a critical problem in the theory of fluid motion for over a century. For Newtonian fluid, the recent monograph of Schmid and Henningson [3] gave a good summary on the hydrodynamic instability. For non-Newtonian fluid such as suspensions, the effects of

∗ Corresponding author at. State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, Hangzhou, Zhejiang 310027, PR China. Tel.: +86 571 87952882; fax: +86 571 87951464. E-mail address: [email protected] (J. Lin).

0020-7462/$ - see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2008.06.006

fibers on the mechanisms of flow instability have not been fully understood yet. The mechanisms of instability occurred in the flows of non-Newtonian fluid is different from that of Newtonian fluid because the non-Newtonian fluid shows some particular characteristics such as stress relaxation [4], pressure and shear dependent viscosities [5] and shear thinning [6]. There are some studies devoted to the hydrodynamic instability in the flows of non-Newtonian fluid. Azaiez and Homsy [7] investigated the influences of viscoelasticity on the instability in the free shear flow by examining three different rheological models. They derived a modified Rayleigh equation, the solution of which shows that viscoelasticity reduces the instability of the flow but does not suppress it. Frigaard et al. [8] obtained the relation between the minimum Reynolds number and the Bingham number by analyzing the instability in the plane Poiseuille flow of a Bingham fluid with 2-D disturbances introduced, and it is found that the minimum Reynolds number for linear instability increases almost linearly with increasing Bingham number. Sureshkumar and Beris [9] pointed out that the viscoelastic stresses first destabilize the plane Poiseuille flow and then stabilize it while the elasticity reached a threshold value, and the presence of a non-zero solvent viscosity has a pronounced stabilizing effect on the flow. Meulenbroek et al. [10] suggested that Poiseuille flow of viscoelastic polymer fluids exhibits a non-linear “subcritical” instability due to normal stress effects, with a threshold which decreases for increasing Weissenberg number. Nouar and Frigaard [11] studied the non-linear stability of Bingham fluid Poiseuille flows in pipes and plane channels, and showed that the critical Reynolds number for transition increases with Bingham number. Siddheshwar and Pranesh [12] investigated the Rayleigh–Benard situation in Boussinesq–Stokes suspensions using both linear and non-linear stability analyses, and delineated the effect of suspended particles on convection against the

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

background of the results of the clean fluid. Fetecau [13] presented the exact solutions corresponding to two types of unsteady flows of an Oldroyd-B fluid in a channel of rectangular cross-section, the solutions satisfy boththe associate partial differential equations and all imposed initial and boundary conditions. For r or  > 0 they tend toward similar solutions for a Maxwell or second-grade fluid. If both r and  > 0, the solutions for Navier–Stokes fluids are recovered. As regards the hydrodynamic instability in fiber suspensions, it has been studied for different type of flows. Pilipenko et al. [14] discussed the linear instability in the Couette flow of fiber suspensions between coaxial cylinders, and found that the addition of fibers to a flow can have a destabilizing effect on the flow. Macmillan [15] studied the micro-structural relaxation equation governing fiber suspension and obtained the formulae determining the linearized stability of the slow, steady, affine flow from a reduced set of constitutive parameters. Nsom [16] considered the flow of a fiber suspension in a curved channel under a pressure gradient, and computed the critical axial wavenumber and Reynolds number characterizing the occurrence of the Dean vortices in the medium against the annuluswidth, and the fiber aspect ratio in the dilute range. Galdi and Reddy [17] studied the well posedness of the equations for fiber suspensions using the second- and fourth-order orientation tensors to account for the presence of fibers, and demonstrated that the linear closure relation leads to anomalous behavior in that the rest state of the fluid is unstable for certain ranges of the fiber particle number. For the quadratic closure relation, it is shown that a unique solution exists locally in time for small data. Munganga and Reddy [18] examined the circumstances, in the energetic sense, under which the constitutive equation based on the use of orientation tensors for fiber suspensions are consistent with the second law of thermodynamics, and the conditions under which the fiber suspensions are stable. They demonstrated that the constitutive equations with use of the natural closure are consistent with the quadratic closure in the 2-D case. Azaiez [19] presented a linear stability analysis of the mixing layer in the presence of fiber additives using a formulation based on moments of the probability distribution function to determine the particle orientation, it is found that the stabilizing effects arise from the orientation diffusion due to hydrodynamic interactions, and not from the anisotropy induced by the presence of fibers in the flow, as speculated before. Azaiez [20] discussed the differences and similarities in the mechanisms of reduction of instability in free shear flows induced by polymer and fiber additives, it is found that, for a fiber suspension with hydrodynamic interactions, the shear stress disturbance induced by the misalignment of the fibers, is the main driving term behind the decrease of the flow instability. On the other hand, the normal stress disturbance acts as a destabilizing factor. Azaiez [21] considered the linear stability of the mixing layer of a fiber suspension using a model based on moments of the orientation distribution function, it is found that for a planar orientation, the trends toward decreasing the flow instability are similar in the case of the natural and hybrid closure approximations, even though the former leads to a less unstable flow. Gupta et al. [22] analyzed the linear instability of the Taylor–Couette flow of semi-dilute nonBrownian suspension. They found that the fiber additives suppress the centrifugal instability, which is attributed to the fact that the suspension develops negative first and second normal stresses in the Taylor–Couette flow. Munganga and Reddy [23] established the existence and uniqueness of solutions to the governing equations for fiber suspensions, and the existence of a unique classical solution, local in time, is proved for the cases of both linear and quadratic closure rules. They found that, by restricting consideration to the physically significant case of constant rotary diffusivity, it is possible to obtain the first results on global existence for this problem. Wan and Lin [24] performed a parametric analysis to the linear instability of the Taylor–Couette flow using the model of the anisotropic fluid,

899

and showed that the fiber additives have a stabilizing effect on the flow, and this effect is clearer for higher values of fiber volume fraction, fiber aspect ratio and lower values of two cylinder radius ratio. Lin et al. [25] demonstrated that, in the mixing layer, the existence of fibers makes the coherent structure evolve quickly. The coherent structures cause the additional stress induced by the fibers to be heterogeneous, and the degree of heterogeneity increases with the increasing of Stokes number and fiber aspect ratio. The maximum and minimum additional stresses occur on the edge of the coherent structures. You et al. [26] presented a linear instability analysis of the channel flow in the presence of fiber additives, it is found that the flow instability of fiber suspensions is governed by two parameters, i.e. a ratio of the axial stretching resistance of fibers to the inertial force of the fluid, and the orientational diffusivity coefficient which accounts for inter-fiber hydrodynamic interactions. They showed that fiber additives have a stabilizing effect on the flow, and this effect is more obvious at higher Reynolds number. Sharifi and Azaiez [27] investigated the effects of high aspect-ratio fiber additives on the temporal hydrodynamic instabilities of free shear flows, and the results showed that higher volume fraction and/or larger fiber aspect-ratio weaken the roll-up of the flow. More importantly, it is found that weak hydrodynamic interactions lead to fundamental structural changes in the flow and result in either a spatial phase shift of the vorticity or the growth of higher harmonics that compete with the fundamental and prevent the completion of the roll-up. You and Lin [28] examined the effect of different closure models of fiber orientation tensor on the flow instability of fiber suspensions in a channel flow, it is found that the difference in closure models has no effect on the behavior of hydrodynamic instability, however, the attenuation of flow instability is most distinct using 3-D hybrid model, and least apparent using its 2-D counterpart for that the fibers show a tendency toward alignment with the flow direction in this case. All these investigations mentioned above are performed based on the theory of linear instability. However, an accurate prediction and description of instability and transition is somewhat beyond the capability of a linear theory. Because of the inherent importance of non-linear effects during the evolution of perturbations, a successful predictive scheme would require some non-linear dependence on an amplitude parameter. Therefore, in spite of the evidence for effects of fiber additives on the flow instability in the channel flow based on the linear theory [24,26], the aim of the present work is to study the hydrodynamic instability of fiber suspensions in the channel flow based on the weakly non-linear theory that has been developed for the Newtonian fluids [29–31]. 2. Equations of fiber suspensions Flow of a fluid between two stationary parallel plates, when the fluid set in motion by the application of a constant pressure gradient is termed Poiseuille flow. Suppose that the fiber suspension is bounded by two parallel plates at y = −d and y = d as shown in Fig. 1. For incompressible and isothermal flow, the basic governing equations are the continuity and momentum equations: ∇ · u = 0,   ju + u · ∇u = −∇p + ∇ · s,  jt

(1) (2)

where u is the velocity of fiber suspensions,  is the fluid density, p is the isotropic pressure, and s is the stress tensor. The boundary conditions are u|y=±d = 0.

(3)

In the current work, the flow of the fiber suspensions is considered a continuum, and the apparent stress is composed of contributions from both the Newtonian suspending fluid and the anisotropic

900

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

z

y

Vorticity P

U0 (y) d o

U

θ

x

y Gradient

Fig. 1. Schematic diagram of the Poiseuille flow.

φ

particles of fibers. Hence, the stress tensor s can be divided into two terms:

s = ss + sf .

x

Flow

(4) Fig. 2. Orientation of a single fiber described by a unit vector p.

The first term, ss , corresponds to the contribution of the suspending fluid

ss = ,

(5)

where  is the dynamic viscosity,  is the rate of strain tensor:

 = ∇u + ∇uT .

(6)

The second term on the right-hand side of Eq. (4), sf , represents the contribution of the fibers. Once this stress is expressed in terms of the flow variables, the flow behavior of the whole system is then determined. 3. Equation of orientation tensors of fibers In order to determine the fiber stress f , we need to know the orientation distribution of the fibers in the flow. The fibers are regarded as rigid and axisymmetric cylinders. Assume the spatial distribution of fibers to be uniform. A unit vector p, as shown in Fig. 2, is used to denote the orientation of a single fiber. The orientation state of fibers at a position can be described by a probability distribution function, (p). This function is defined so that the probability of finding a fiber between p and (p + dp) is given by (p) dp. If we assume that fibers move with the bulk motion of the fluid then (p) may be regarded as a convected quantity. The continuity condition is then D + ∇ · (p˙ ) = 0. Dt

(7)

The second-order and fourth-order orientation tensors are defined as follows:  (a2 )ij = aij = pi pj (p) dp, (10)  (a4 )ijkl = aijkl =

pi pj pk pl (p) dp.

(11)

Combining Eq. (11) with Eqs. (7) and (8) we have the equation of change for the second-order tensor a2 : 1 Da2 = − ( · a2 − a2 · x) Dt 2 +

 2

( · a2 + a2 ·  − 2 : a4 ) + 2Dr · (I − ma2 ),

(12)

where I is the unit tensor and m is the dimension of the space. Eq. (12) contains the unknown fourth-order tensor a4 . Thus in order to obtain a closed set of equations, a closure approximation should be adopted. The second-order tensor a2 is the most concise nontrivial description of orientation. Therefore, it is conceivable that the fourth-order tensor a4 is replaced by a closure approximation based on the second-order tensor a2 . There exist a number of different closure approximations. In the recent work, the quadratic closure approximation [22,36,37] as shown in Eq. (13) is adopted. aijkl = aij akl .

(13)

By Eqs. (12) and (13), the orientation state of fibers in the flow is completely determined. It is then possible to express the fiber contribution to the total stress.

The equation of motion of a single fiber is given by [32]

 1 1 p˙ = − x · p + ( · p −  : ppp) − Dr · ∇ , 2 2 

(8)

x = ∇uT − ∇u,

(9)

4. Model of stress caused by the fibers

is the vorticity tensor,  = (H2 − 1)/(H2 + 1) with H being

where x the aspect ratio of fibers, Dr is the rotary diffusivity resulting from inter-fiber hydrodynamic interactions [33]. Folgar and Tucker [34] suggested that Dr is isotropic and can be replaced by CI ||I, where  || = tr( · )/2, and the coefficient CI depends on the fiber concentration, shape and aspect ratio. The use of (p) is computationally expensive, and a more efficient approach consists of using a set of orientation tensors, defined as dyadic products of p averaged over all possible orientations [35].

There are several models for the stress tensor sf caused by the existence of the fibers. We use the following model [38–41]:

sf = c4 : ,

(14)

where c4 is an anisotropic viscosity tensor and can be expressed as cijkl = Baijkl + B1 (aij kl + akl ij ) + B2 (aik jl + ail jk + ajl ik + ajk il ) + B3 ij kl ,

(15)

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

where ij is the Kronecker symbol, the coefficients of B and Bi are the material constants. Bi is equal to zero for the suspensions with weak or no Brownian motion, and B is B=

 H2  , 3 ln( 2/ )

(16)

where is the volume fraction of fibers. For brevity, we denote F = a4 : . Then Eq. (14) can be written as based on Eq. (13):

sf = BF,

F = a2 a2 : .

(17)

Substituting Eq. (27) into Eqs. (18), (19), (12), (17) and subtracting the base-state equations (20)–(21), (23)–(26), we obtain the equations of perturbation with non-linear terms of perturbations neglected.

ju jv + = 0, jx jy

(28)

1 ju ju dU0  jp + + v = − + jt jx dy jx Re +

5. Equation of base-state flow Define d/2 and U as the characteristic length and velocity, respectively, and then the characteristic time is d/2U and the Reynolds number is Re = Ud/2 , where is the kinematic viscosity of the fluid. Substituting Eqs. (4)–(6), (17) into Eq. (2) then non-dimensionalizing with respect to the characteristic length and velocity scale, we have ∇ · u = 0,

(18)

1 2 ju B

+ u · ∇u = −∇p + ∇ u+ ∇ · F. jt Re Re

(19)

For the Poiseuille flow as shown in Fig. 1, the base-state streamwise and transverse velocities are U0 = U0 (y) and V0 (y) = 0, respectively. Substituting the base-state velocities into Eq. (19) yield −

1 j2 U0 j P0 B jF0xy + = 0, + jx Re jy2 Re jy

(20)

−

jP0 B jF0yy + = 0, jy Re jy

(21)

where P0 is the base-state pressure, F0xy and F0yy are the components of the base-state F0 = a04 : 0 . Substituting Eq. (13) into F0 yields F0 = a02 a02 : 0 ,

(22)

where the components of the base-state orientation tensor a02 can be presented as follows based on Eq. (12): (1 + )a0xy − 2a0xx a0xy + 2CI (1 − ma0xx ) = 0,

(23)

(1 + )a0yy − (1 − )a0xx − 4a20xy − 4CI ma0xy = 0,

(24)

( − 1)a0xy − 2a0xy a0yy + 2CI (1 − ma0yy ) = 0,

(25)

and the corresponding components of Eq. (22) are given by dU F0xx = 2a0xx a0xy 0 , dy F0xy = 2(a0xy )2

dU0 . dy

B

Re

1 jv jv jp + U0 = − + jt jx jy Re +

B

Re









j2 u j2 u + jx2 jy2





  jFxy jFxx + , jx jy

j2 v j2 v + jx2 jy2

(29)





  jFxy jFyy + , jx jy

(30)

jaxx ja + U0 xx jt jx

   dU jv j u 1 −2 0 axy + 2a0xy − 2 dy jx jy     dU ju jv j u 1 +  2 0 axy + 4a0xx + 2a0xy + 2 dy jx jx jy

   jv ju  + 2CI + −2Fxx jx jy     jv ju dU − 6CI 0 axx , − 6CI a0xx + jx jy dy

=−

(31)

jaxy jaxy + U0 jt jx

   dU 1 jv ju − 0 (ayy − axx ) + (a0yy − a0xx ) − 2 dy jx jy dU0  1 (a + ayy ) + (a0xx + a0yy ) +  2 dy xx 

  jv ju dU  − 2Fxy − 6CI 0 axy × + jx jy dy     jv ju + , − 6CI a0xy jx jy

=−

(32)

jayy jayy + U0 jt jx

dU0 , dy

F0yy = 2a0xy a0yy

901

   dU 1 jv ju 2 0 axy − 2a0xy − 2 dy jx jy    dU jv j u 1 + +  2 0 axy + 2a0xy 2 dy jx jy

    jv jv j u  + 2CI +4a0yy − 2Fyy + jy jx jy    jv j u dU0  + − 6CI a , − 6CI a0yy jx jy dy yy  = 2 dU0 (a Fxx a + a0xy axx ) dy 0xx xy   j u jv + a0yy + a0xx a0xx jx jy     jv ju +a0xx a0xy + , jx jy =−

(26)

6. Equation of perturbation We examine the hydrodynamic instability of the flow by introducing perturbations into the base flow. Let u , v , p , a2 , and F  denote perturbations of u, p, a2 , and F for the Poiseuille flow, respectively. Then the velocity, pressure, a2 , and F in Eq. (19) can be represented by the base-state profile plus a perturbation: ⎧ u = U0 + u (x, y, t), ⎪ ⎪ ⎪ ⎪  ⎪ ⎪ ⎨ v = v (x, y, t), (27) p = P0 + p (x, y, t), ⎪ ⎪  ⎪ ⎪ a2 = a20 + a2 (x, y, t), ⎪ ⎪ ⎩ F = F0 + F  (x, y, t).

(33)

(34)

902

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

  dU0  ju jv  = 2 2a Fxy axy + a0xy a0xx + a0yy 0xy dy jx jy     jv ju + , +(a0xy )2 jx jy  = 2 dU0 (a Fyy a + a0yy axy ) dy 0xy yy   ju jv + a0yy a0xx + a0yy jx jy     jv ju + . +a0xy a0yy jx jy

EVF –D = (35)



(36)

jv jv jp¯ B

+ v =− + u jx

jy

jy



∗ ∗ jFxy jFxx + , jx jy



Re

(37)



∗ ∗ jFxy jFyy + , jx jy

(38)

where F ∗ = tr(a02 · 0 )a02 + tr(a02 ·  )a2 + tr(a2 · 0 )a2 + tr(a2 ·  )a02 ,

(39)

where a20 and 0 are the base-state quantities of a2 and , respectively;  is the perturbation of . When the perturbations with finite amplitude are amplifying or decaying, the mean motion also suffers modification in order to maintain the energy balance. Therefore, the ¯ t) and p¯ in Eqs. (37) U0 and P0 of base-state are replaced by u(y, and (38). 8. Energy balance relation Substituting Eq. (27) into Eqs. (18) and (19) and separating out the mean and perturbation parts of the motion, we can obtain the energy balance relation by integrating the resulting perturbation equations: ET = EP − EV –N − EVF –B − EVF –D , where

j jt

×

(40)

 

1  2 [(u ) + (v )2 ] dx dy, 2     j u 2B

(−u v ) − (a0xx − a0yy )axy EP = Re jx   

   ju jv ju + axx +2a0xy axy + jy jx jx ET =

ju¯ dx dy, jy 1 Re

EVF –B =

2B

Re

(41)

(42)

  

E V –N =

ju jv − jy jx

2 dx dy,



ju jv + jy jx

ju jv + jy jx

j u jx

 2 dx dy.

(45)

axx + ayy = 0.

(46)

The term of left-hand side, ET , in Eq. (40) denotes the rate of increase of the perturbation energy within the volume considered. On the right-hand side, the first term EP gives the rate of transfer of energy from the mean flow to the perturbation through Reynolds stress, and the product of the perturbation velocity gradient and the perturbed fiber orientation components. The other three terms EV –N , EVF –B and EVF –D , which are necessarily always positive, represent the viscous dissipation of energy due to the suspending fluid and fiber additives. Within them, EVF –B arises from the product of the base-state fiber orientations and the perturbation velocity gradient, and EVF –D denotes the effects of higher order terms associated with the perturbed fiber orientation. 9. Equation of perturbation amplitude The most obvious flaw of a linear stability analysis is that the outcome is independent of the initial perturbation amplitude. That quantity is an arbitrary constant in the eigenvalue problem. A second defect is that unsteady perturbations amplify or decay exponentially forever-behavior that is not in accord with observation. These and other defects are remedied with the weakly non-linear theory. In this section, we describe the details of explicit non-linear expansion for the instability of fiber suspension in the Poiseuille flow. Since we consider periodic modulations of the pressure and fiber orientations, so there are periodic pressure and fiber orientation modulations but their average values remain unchanged. Normally, a linear eigenmode of the equations is the form ⎛  ⎞ ⎛ ∗ ⎞ u (x, y, t) u (y) ⎜  ⎟ ⎜ ∗ ⎟ ⎝ p (x, y, t) ⎠ = ⎝ p (y) ⎠ exp[i (x − ct)], a2 (x, y, t)

(47)

a∗2 (y)

where the temporal mode is considered and the wavenumber of the streamwise is treated as real and c is a complex phase velocity c = cr + ici whose real part represents the propagating speed of the perturbation in the streamwise direction and its imaginary part is proportional to the temporal growth rate, respectively. All components of u∗ (y), p∗ (y) and a∗2 (y) in Eq. (47) can be obtained directly from the results of the linear stability analysis. In the present study, we confine the analysis to the temporal evolution of the first harmonic component of the perturbation and higher harmonic components having wavenumbers n with n being an integral, and interchanges of energy between them and both the mean flow and the first harmonic component are ignored. This is physically realistic since the high wavenumbers correspond to rapidly oscillating waves and such a function contains little energy at high wavenumbers. So the non-linear evolution of the amplitude of a single mode has a following form: u (x, y, t) = A(t)u∗ (y) exp[i (x − ct)] + c.c.,

(44)

where c.c. means their complex conjugate. Then, in an expansion in powers of A(t), whose initially exponential growth according to the linear theory will be modified as non-linearity becomes significant and probably approaches an equilibrium value. In order to determine the lowest non-linear order of the equation for A(t), we obtain the

 2 dx dy,

(axx − ayy )

(43)

  j u (a0xx − a0yy ) jx

+a0xy



a0xx + a0yy = 1,

Note that u (x, y, t), v (x, y, t) are the perturbed velocities with finite amplitude in the weakly non-linear theory, the non-linear terms of perturbations are not negligible. Substituting Eq. (27) into Eq. (19) and taking the average over the terms in Eq. (19), we have

1 j2 u¯ jp¯ B

+ + =− jx Re jy2 Re

+axy

 

To derive the previous expressions, the continuity equations and following relations are used:

7. Equation of perturbations with finite amplitude

ju¯ j u j u + u + v jt jx jy

2B

Re

(48)

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

following equation of amplitude by substituting Eq. (48) into Eqs. (40)–(45):

j A2  = 2 A2 − 3 ReA4 − 4 A2 1 jt Re

−

Table 1 Comparison with others

Orszag [43] Thomas [44] Nachtsheima Groscha Kirchner [45] Present work

+ 2B (1 A2 + 2 A4 + 3 A3 ) 2B

(1 A2 + 2 A2 + 3 A + 4 A3 + 5 A2 ), Re

903

a

c

Rec

1.02056 1.026 1.02 1.025 1.026 1.021

5772.22 5780 5767 5750 5775.99 5772.245

Reported by Kirchner [45].

(49) where the coefficients i , i and i are given in Appendix A. In particular, Eq. (49) degenerates into the amplitude equation for a Newtonian fluid derived by Stuart. 10. Numerical method In order to determine the perturbation components u∗ (y), p∗ (y) and a∗2 (y) in Eq. (47), we apply the Chebyshev spectral collocation method [42]. This method appears now to be the standard approach for solving this type of eigenvalue problem due to its simplicity, efficiency and accuracy. The collocation points in the computational direction are the Gauss–Lobatto points: yj = cos

(j − 1) N−1

,

j = 1, 2, . . . , N,

(50)

and then the perturbation variable u∗ (y), p∗ (y) and a∗2 (y) can be represented as Chebyshev series: u∗ (y) =

N−1 

uˆ j Tj (y),

(51)

vˆ j Tj (y),

(52)

j=0

v∗ (y) =

N−1  j=0

p∗ (y) =

N−1 

N−1 

11. Verification of numerical code The codes used here were validated by comparisons of the results with available literature. The critical Reynolds number Rec and the critical wavenumber c in the Newtonian flow were calculated, and the corresponding results as well as other results are shown in Table 1, in which Orszag's results [43] were given using the Chebyshev collocation, and Thomas's results [44] were obtained using fourth-order finite difference with a truncation error involving the eighth derivative because they have tabulated their calculations for a range of Reynolds numbers and wavelengths of the disturbances. From Table 1 it can be seen that our results are in very good agreement with their results. 12. Results and discussions

pˆ j Tj (y),

(53)

aˆ 2j Tj (y),

(54)

j=0

a∗2 (y) =

points. In order to obtain specific values of the amplitude at a given Reynolds number and the fiber parameters, Eq. (49) should be solved. Note that A is finite and bounded but sufficiently small. As the perturbation amplitude approaches zero, the terms with order of O(A2 ) of the non-linear solution represent the solution of the linear stability problem. Thus, Eq. (49) can be calculated using bisect method in the range of [0, M], where M is initially set a relative small value to safely get a correct amplitude.

j=0

where Tj (y) is the Chebyshev polynomials of degree N. Substituting Eqs. (47), (51)–(54) into Eqs. (28)–(36), we can determine the discrete linear system with the clamped boundary conditions imposed. The spatial derivatives of the perturbations at collocation points are accomplished by matrix operations. The derivative matrices can be computed either by directly differentiating the interpolating polynomial constructed on the set of collocation points or by means of a transform method. The governing equations in the discrete sense can be represented as a complex generalized algebraic eigenvalue problem: H · U = G · U,

(55)

where  = −i c and the vector U contains the real and imaginary parts of the coefficients uˆ j defined in Eq. (47). H and G are constant square matrices containing the derivatives of the Chebyshev polynomials, and dependent on the base-state of flow, wavenumber , complex phase velocity c, Reynolds number Re, and parameters of fiber suspensions. The number of collocation points used in the present work is sufficient to obtain nine decimal places of accuracy in the computed values of . With this procedure, the Chebyshev–Gauss–Legendre quadrature is used to compute integrals as it includes the endpoints of the finite interval among the grid

Before studying the non-linear instability on the fiber suspensions, the linear instability analysis should be done firstly. The previous computational study on the linear instability of different flows such as free mixing layer and Taylor–Couette flow shows that the fiber additives will increase the critical Reynolds number [21,22]. Fig. 3 depicts the variations of the growth rate of the perturbation r versus the wavenumber in the case of Newtonian flow and case of fiber suspensions for Re = 104 , = 0.1 × 10−4 and H = 102 . In the present work, we take the value of parameter CI , which is related to the fiber interactions, as 0.001 according to the experimental observation of the rotary diffusion between particles in the shear flow [33]. From Fig. 3 we can see that the maximum growth rate in the case of fiber suspensions is larger than that of Newtonian flow, and the region of unstable wavenumbers for fiber additives is reduced from that of the Newtonian flow, indicating that the addition of fibers in the plane Poiseuille flow leads to the reduction of the flow instability. This phenomenon, also found in the mixing layer and Taylor–Couette flows, becomes significant for large volume concentration and aspect ratio H. Corresponding to Fig. 3, Fig. 4 shows the change of the amplitude A in the equilibrium as a function of wavenumbers , in which Af , AN denote the amplitude in fiber suspension and in Newtonian flow, respectively. For the same wavenumber, Af is always associated with a larger AN but the overall behavior of the equilibrium amplitude in fiber suspension is in accordance with that of Newtonian flow except for the difference of the wavenumber achieving the peak amplitude. It is known that the equilibrium amplitude A, in the Newtonian Poiseuille shear flow, increases with the Reynolds number. The similar trends are observed in the case of the fiber suspensions as shown

904

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

0.012

0.004

A

λr

0.003

0.002

0.011

Newtonian flow Fiber suspensions

0.001

0.01

0 0.8

0.9

1

α

1.1

Fig. 3. The growth rate of the perturbation versus the wavenumber in the case of Newtonian flow and fiber suspensions (Re = 104 , = 0.1 × 10−4 and H = 102 ).

8000

10000

12000 Re

14000

16000

Fig. 5. The amplitude in the equilibrium versus the Reynolds number ( = 0.8,

= 0.1 × 10−4 and H = 102 ).

0.012 Af

1

AN

0.03

0.02

Af

0.6 u

AN

0.8 0.011

Laminar Poiseuille flow

0.4

0.01

Modified fiber flow

0.01 0.8

0.85

0.9

0.95 α

1

Modified Newtonian flow

0.2

1.05

Fig. 4. The change of the amplitude in the equilibrium as a function of wavenumbers.

in Fig. 5 where the parameters are =0.8, =0.1×10−4 and H =102 . The equilibrium amplitude A rises from primarily small value at Re = 8 × 103 to a maximum value at Re = 1.7 × 104 . In the state of equilibrium, the velocity profile of mean flow is distorted from its original laminar form. The distortion of the mean profile is crucial because it may excite the appearance of small wavelength, highly unstable waves and causes a further modification of the mean velocity profile. Integrating Eqs. (37)–(38) and applying the boundary conditions that the perturbations vanish the wall, we have u¯ = 1 − y2 + Re − 2B

 y 1

 y 1

u v dy



a − 2(axy )2 y + Re[(a0xy )2 + (axy )2 ]u v

+ (a0xx − a0yy )axy

×

j u + 2a0xy jx

ju j u jv + axy + axy axx jx jy jx

 dy,

(56)

where

a = {1 + 2B [(a0xy )2 + (axy )2 ]}−1 .

(57)

0

-1

- 0.5

0 y

0.5

1

Fig. 6. The modified velocity profiles (Re = 104 , = 1.02, = 0.1 × 10−3 , H = 102 ).

If we set fiber parameter B = 0 and neglect the Reynolds stress in Eq. (56), the parabolic velocity profile can be obtained. Since the perturbations dissipate and extract the energy from the mean motion through the perturbed Reynolds stress, the velocity profile in the perturbed fiber suspensions is different from that in the flow with original infinitesimal perturbations. Fig. 6 shows the modified mean velocity u¯ in the channel for Re = 104 , = 1.02, = 0.1 × 10−3 and H = 102 . For the sake of comparison we also plot the results for Newtonian flow and the laminar Poiseuille flow. It can be seen that the modified velocities are smaller than that in the laminar Poiseuille flow. Note that the modified velocity of fiber suspensions is larger than that in the Newtonian flow.  Fig. 7 shows the root-mean-square, (u )2 , of streamwise component of perturbation velocity u∗ , as a function of y for Re = 104 ,

= 1.02, = 0.1 × 10−4 and H = 102 . The subscript “N” denotes the values in Newtonianflow. We can find that both curves have a

similar shape, however, (u∗ )2 N has about one order of magnitudes larger than that of the fiber suspension, i.e. the root-mean-square  2  (u ) are significantly attenuated by the addition of fibers.

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

905

0.05 ED

0.03

0.0012

ED

(u*)2 (u*)2N

(u*)2

0.03

0.0008

0.02

EB

EB

0.04

0.0004

0.01

0.02 0.02

0.04

0.06

0.08

-3 0.1(x10 )

φ Fig. 9. Variations of EB and ED versus the fiber volume fraction (Re = 104 , = 0.8 and H = 102 ).

0.01

-0.5

0 y

0.5

0.0015

1 0.05

Fig. 7. Root-mean-square of streamwise component of perturbation velocity (Re=104 ,

= 1.02, = 0.1 × 10−4 , H = 102 ).

EB

0.04 ED

0.02

ED

0.001

0.03

EB

-1

0.02

0.0005

0.01

20

40

60

80

100

(v*) 2

H Fig. 10. Variations of EB and ED versus the fiber aspect ratio (Re = 104 , = 0.8 and H = 102 ).

0.01

is much smaller than ED , indicating that the basic orientation of the fibers has a smaller effect on the dissipation of the system energy. The relationships between EB , ED and fiber aspect ratio, H, are shown in Fig. 10. We can see that EB and ED also increase with the fiber aspect ratio.

(v*)2 (v*)2N

0

13. Conclusions

-1

-0.5

0 y

0.5

1 4

Fig. 8. Root-mean-square of normal component of perturbation velocity (Re = 10 ,

= 1.02, = 0.1 × 10−4 , H = 102 ).

 The root-mean-square, (v )2 , of normal component of perturbation velocity v∗ are shown in Fig. 8. It was found that the rootmean-square of perturbation velocity has a comparable magnitude with their corresponding turbulent intensity given by Lin et al. [46]. The energy balance of the perturbation can provide valuable information about the mechanisms behind the attenuation of the flow instability. We further examine the relative importance of all terms in Eq. (40) on the budget of the kinetic energy. In the equilibrium state, the time derivative of the kinetic energy is dependent on the balance between the inertial and dissipative terms. EVF –B and EVF –D in Eq. (40) are attributed to the contributions of the fibers, they are normalized with respect to the energy source Ep and expressed by EB and ED , respectively. Fig. 9 shows the variations of EB and ED versus the fiber volume fraction for Re = 104 , = 0.8 and H = 102 . It can be seen that both EB and ED increase with increasing , but EB

The evolution equation of the amplitude of a finite perturbation for non-linear instability of plane Poiseuille flow of the fiber suspensions is derived. In this process, the compact quadratic closure approximation of the fiber orientation is used. The equation is solved numerically. The results show that the maximum growth rate in the case of fiber suspensions is larger than that of the Newtonian flow, and the region of unstable wavenumbers for fiber additives is reduced compared to the Newtonian flow. In the equilibrium state, the amplitude of the perturbation in the fiber suspensions increases with the Reynolds number, and is smaller than that in the Newtonian flow. The velocity in the perturbed fiber suspensions is smaller than that in the laminar Poiseuille flow. The root-mean-squares of streamwise and normal components of perturbation velocity are significantly attenuated by the addition of fibers. The energy EB and ED , related to the fiber orientations, increase with increasing fiber volume fraction and aspect ratio. The basic orientation of the fibers has a smaller effect on the dissipation of the system energy. Based on the above results, it is concluded that the addition of fibers in the plane Poiseuille flow leads to the reduction of the flow instability. This phenomenon becomes significant for large volume concentration and aspect ratio.

906

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

Acknowledgment The authors gratefully acknowledge the support of the Major Program of the National Natural Science Foundation of China with Grant no. 10632070.

⎧  1 ⎨ 2 = 4 4 2 [(aˆ xxR )2 + (aˆ xxI )2 ](Uˆ R2 + Uˆ I2 ) −1 ⎩ " + 4(aˆ xxR aˆ xyR + aˆ xxI aˆ xyI ) 2 (Uˆ R Vˆ R + Uˆ I Vˆ I )

# dUˆ I dUˆ R ˆ ˆ − UI + UR + [(aˆ xxR )2 + (aˆ xxI )2 ] dy dy ⎧ ⎡

2

 ⎤ ⎨ ˆ ˆ 2 d U d U R I ⎦ + 2 (Vˆ 2 + Vˆ 2 ) × ⎣ + R I ⎩ dy dy

Appendix A. Coefficients of amplitude equation ˆ V, ˆ In the following expressions, the complex eigenfunctions U, ˆaxx , and aˆ xy corresponding to disturbance quantities u∗ and a∗ are all represented as real and imaginary parts such as Uˆ = Uˆ R + iUˆ I , etc.

1 =

 1 −1

[(Uˆ 2R + Uˆ 2I ) + (Vˆ 2R + Vˆ 2I )] dy,

 1 du¯ ˆ ˆ (UR VR + Uˆ I Vˆ I ) dy, 2 = −2 −1 dy

3 = 4

 1 −1

3 = (A.3)

1 = 4

+ 2 (Vˆ R2 + Vˆ I2 ) dy, ⎭

 1 du¯ ˆ a [(aˆ xyR )2 + (aˆ xyI )2 ](Uˆ R Vˆ R + Uˆ I Vˆ I ) dy, −1 dy

2 = 4Re

 1

× (Uˆ R Vˆ R + Uˆ I Vˆ I )2 dy,

3 = 2

 1 −1

(A.4)

(A.5)

ˆ a {(a0xy )2 + 2[(aˆ xyR )2 + (aˆ xyI )2 ]}

−1

(A.6)

ˆ a ˆ auv (Uˆ R Vˆ R + Uˆ I Vˆ I ) dy,

(A.7)

where

ˆ a = {1 + 2B {(a0xy )2 + 2[(aˆ xyR )2 + (aˆ xyI )2 ]}}−1 , ˆ auv = 2 (a0xx − a0yy )(Uˆ R aˆ xyI − Uˆ I aˆ xyR ) 

+ 4a0xy [(Uˆ R aˆ xxI − Uˆ I aˆ xxR ) + (Vˆ R aˆ xyI − Vˆ I aˆ xyR )]

 dUˆ R dUˆ + aˆ xyI I + aˆ xyR , dy dy

1 = 2

⎧  1 ⎨

2 (a0xx − a0yy )2 (Uˆ R2 + Uˆ I2 ) + (a0xy )2

−1 ⎩

⎧ ⎡

2

2 ⎤ ⎨ dUˆ R dUˆ I ⎦ ⎣ × + + 2 (Vˆ R2 + Vˆ I2 ) ⎩ dy dy

dUˆ I dUˆ − Vˆ I R +2 Vˆ R dy dy

⎭⎭

dy,

 1 −1

⎫ ⎬ ⎭

"

 dUˆ I dUˆ − Uˆ I R + 2a0xy (a0xx − a0yy ) Uˆ R dy dy #⎫ ⎬ dy, + 2 (Uˆ R Vˆ R + Uˆ I Vˆ I ) ⎭

(A.8)

5 = −2B

(A.9)

(A.10)

ˆ a {(a )2 {1 − 2B  0xy

ˆ auv dy, + 2[(aˆ xyR )2 + (aˆ xyI )2 ]}}(Uˆ R Vˆ R + Uˆ I Vˆ I ) ⎫ ⎬



 ⎫⎫ ⎬⎬

 1 du¯ ˆ  uv {1 − 4B ˆ [(aˆ xyR )2 + (aˆ xyI )2 ]} dy, −1 dy

4 = 2Re

⎧⎡ ⎤  1 ⎨ ˆ 2 ˆ 2 dUI ⎦ dUR ⎣ 4 = 2 + dy dy −1 ⎩ dUˆ dUˆ I +4 Vˆ I R − Vˆ R dy dy

dUˆ I dUˆ − Vˆ I R +2 Vˆ R dy dy

(A.2)

(Uˆ R Vˆ R + Uˆ I Vˆ I )2 dy,





(A.1)

 1 −1

ˆ auv )2  ˆ a dy. (

(A.11) (A.12)

References [1] K. Kannan, K.R. Rajagopal, Simulation of flow induced crystallization in fiber spinning using the radial resolution approximation, Mech. Adv. Mater. Struct. 12 (2005) 413–424. [2] K. Kannan, K.R. Rajagopal, Simulation of fiber spinning including flow-induced crystallization, J. Rheol. (2005) 683–703. [3] P.J. Schmid, D.S. Henningson, Stability and Transition in Shear Flows, Springer, New York, 2001. [4] K.R. Rajagopal, A.S. Wineman, A note on viscoelastic materials that can age, Int. J. Non-linear Mech. 10 (2004) 1547–1554. [5] J. Malek, J. Necas, K.R. Rajagopal, Global existence of solutions for flows of fluids with pressure and shear dependent viscosities, Appl. Math. Lett. 15 (2002) 961–967. [6] D. Mansutti, K.R. Rajagopal, Flow of a shear thinning fluid between intersecting planes, Int. J. Non-linear Mech. 26 (1991) 769–775. [7] J. Azaiez, G.M. Homsy, Linear stability of free shear flow of viscoelastic liquids, J. Fluid Mech. 268 (1994) 37–69. [8] I.A. Frigaard, S.D. Howison, I.J. Sobey, On the stability of Poiseuille flow of a Bingham fluid, J. Fluid Mech. 263 (1994) 133–150. [9] R. Sureshkumar, A.N. Beris, Linear stability analysis of the viscoelastic Poiseuille flow using an Arnoldi orthogonalization algorithm, J. Non-Newtonian Fluid Mech. 56 (1995) 151–182. [10] B. Meulenbroek, C. Storm, A.N. Morozov, W. van Saarloos, Weakly non-linear subcritical instability of visco-elastic Poiseuille flow, J. Non-Newtonian Fluid Mech. 116 (2004) 235–268. [11] C. Nouar, I.A. Frigaard, Nonlinear stability of Poiseuille flow of a Bingham fluid theoretical results and comparison with phenomenological criteria, J. NonNewtonian Fluid Mech. 100 (2001) 127–149. [12] P.G. Siddheshwar, S.T. Pranesh, An analytical study of linear and non-linear convection in Boussinesq–Stokes suspensions, Int. J. Non-linear Mech. 39 (2004) 165–172. [13] C. Fetecau, C. Fetecau, Unsteady flows of Oldroyd-B fluids in a channel of rectangular cross-section, Int. J. Non-linear Mech. 40 (2005) 1214–1219. [14] V.N. Pilipenko, N.M. Kalinichenko, A.S. Lemak, Stability of the flow of a fiber suspension in the gap between coaxial cylinders, Sov. Phys. Dokl. 26 (1981) 646–648. [15] E.H. Macmillan, Slow flows of anisotropic fluids, J. Rheol. 33 (1989) 1071–1105. [16] B. Nsom, Stability of fiber suspension flow in curved channel, J. Phys. II 6 (1996) 1483–1492. [17] G.P. Galdi, B.D. Reddy, Well-posedness of the problem of fiber suspension flows, J. Non-Newtonian Fluid Mech. 83 (1999) 205–230. [18] J.M.W. Munganga, B.D. Reddy, K.J. Diatezua, Aspects of the thermodynamic stability of fibre suspension flows, J. Non-Newtonian Fluid Mech. 92 (2000) 135–150. [19] J. Azaiez, Linear stability of free shear flows of fiber suspensions, J. Fluid Mech. 404 (2000) 179–209. [20] J. Azaiez, Reduction of free shear flows instability: effects of polymer versus fiber additives, J. Non-Newtonian Fluid Mech. 91 (2000) 233–254.

Z. Wan et al. / International Journal of Non-Linear Mechanics 43 (2008) 898 -- 907

[21] J. Azaiez, Stability of the mixing layer of fiber suspensions: role of the closure approximation and off-plane orientation, J. Non-Newtonian Fluid Mech. 95 (2000) 253–276. [22] V.K. Gupta, R. Sureshkumar, B. Khomami, J. Azaiez, Centrifugal instability of semidilute non-Brownian fiber suspensions, Phys. Fluids 14 (2002) 1958–1971. [23] J.M.W. Munganga, B.D. Reddy, Local and global existence of solutions to equations for flows of fiber suspensions, Math. Models Methods Appl. Sci. 12 (2002) 1177–1203. [24] Z.H. Wan, J.Z. Lin, Hydrodynamic instability of semi-concentration fiber suspensions between two rotating coaxial cylinders, J. Nonlinear Sci. Numer. Simulation 5 (3) (2004) 211–216. [25] J.Z. Lin, Z.S. Yu, X.M. Shao, Coherent structures in the mixing layers of a nonNewtonian fluid, J. Turbulence 5 (2004) 1–18. [26] Z.J. You, J.Z. Lin, Z.S. Yu, Hydrodynamic instability of fiber suspensions in channel flows, Fluid Dyn. Res. 34 (2004) 251–271. [27] F. Sharifi, J. Azaiez, Vortex dynamics of fiber-laden free shear flows, J. NonNewtonian Fluid Mech. 127 (2005) 73–87. [28] Z.J. You, J.Z. Lin, Effects of tensor closure models and 3-D orientation on the stability of fiber suspensions in a channel flow, Appl. Math. Mech. 26 (2005) 307–312. [29] J.T. Stuart, On the non-linear mechanics of wave disturbances in stable and unstable parallel flows: Part 1. The basic behaviour in plane Poiseuille flow, J. Fluid Mech. 9 (1960) 353–370. [30] J. Watson, On the non-linear mechanics of wave disturbances in stable and unstable parallel flows: Part 2. The development of a solution for plane Poiseuille flow and for plane Couette flow, J. Fluid Mech. 9 (1960) 371–389. [31] T. Herbert, On perturbation methods in nonlinear stability theory, J. Fluid Mech. 126 (1983) 167–186. [32] J.S. Cintra, C.L. Tucker, Orthotropic closure approximations for flow-induced fiber orientation, J. Rheol. 39 (1995) 1095–1122.

907

[33] D.L. Koch, A model for orientational diffusion in fiber suspensions, Phys. Fluids 7 (1995) 2086–2088. [34] F.P. Folgar, C.L. Tucker, Orientation behaviour of fibres in concentrated suspensions, J. Reinf. Plast. Comp. 3 (1984) 98–119. [35] S.G. Advani, C.L. Tucker, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 31 (1987) 751–784. [36] G. Marrucci, P.L. Maffettone, Description of the liquid-crystalline phase of rodlike polymers at high shear rates, Macromolecules 22 (1989) 4076–4082. [37] S.G. Advani, C.L. Tucker III, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 31 (1987) 751–784. [38] G.G. Lipscomb, M.M. Denn, D.U. Hur, D.V. Boger, The flow of fiber suspensions in complex geometries, J. Non-Newtonian Fluid Mech. 26 (1988) 297–325. [39] G.K. Batchelor, The stress generated in a non-dilute suspension of elongated particles by pure straining motion, J. Fluid Mech. 46 (1971) 813–829. [40] E.S.G. Shaqfeh, G.H. Fredrickson, The hydrodynamic stress in a suspension of rods, Phys. Fluids 2 (1990) 7–24. [41] S.M. Dinh, R.C. Armstrong, A rheological equation state for semi-concentrated fibre suspensions, J. Rheol. 28 (1984) 207–227. [42] L.N. Trefethen, Finite difference and spectral methods for ordinary and partial differential equations, unpublished text, 1996, available at, http://web.comlab.ox.ac.uk/oucl/work/. [43] S.A. Orszag, An accurate solution of the Orr–Sommerfeld equation, J. Fluid Mech. 50 (1971) 689–703. [44] L.H. Thomas, The stability of plane Poiseuille flow, Phys. Rev. 91 (1953) 780–783. [45] N.P. Kirchner, Computational aspects of the spectral Galerkin FEM for the Orr–Sommerfeld equation, Int. J. Numer. Methods Fluids 22 (2000) 119–137. [46] J.Z. Lin, W.F. Zhang, Z.S. Yu, Numerical research on the orientation distribution of fibers immersed in laminar and turbulent pipe flows, J. Aerosol Sci. 35 (2004) 63–82.