Development length of sinusoidally pulsating laminar pipe flows in moderate and high Reynolds number regimes

International Journal of Heat and Fluid Flow 37 (2012) 167–176

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Development length of sinusoidally pulsating laminar pipe flows in moderate and high Reynolds number regimes Subhashis Ray a, Bülent Ünsal b,⇑, Franz Durst c a

Institute of Thermal Engineering, Technische Universität Bergakademie Freiberg, Gustav-Zeuner Straße 7, 09596 Freiberg, Germany National Metrology Institute, TUBITAK UME, Pk.54 Gebze, Kocaeli, Turkey c FMP Technology GmbH, Am Weichselgarten 34, 91058 Erlangen, Germany b

a r t i c l e

i n f o

Article history: Received 17 November 2011 Received in revised form 30 January 2012 Accepted 4 June 2012 Available online 25 June 2012 Keywords: Development length Unsteady Pipe flow Pulsating Developing

a b s t r a c t Recent studies on laminar, fully-developed, sinusoidally pulsating pipe flows have revealed the existence of a unique signature map that can be used for the measurement of the arbitrary time-varying, instantaneous mass flow rate from the recorded axial pressure gradient data and vice versa. This measuring technique is, however, valid for the hydrodynamically fully-developed flow of an incompressible fluid. The present study, therefore, deals with the numerical evaluation of the required development length as functions of the mean Reynolds number, the amplitude of mass flow rate pulsation and the pulsation frequency in the moderate and high Reynolds number regimes. The investigation shows that in the low-frequency, quasi-steady regime, the instantaneous variations of L/D can be predicted by the steady-state results for corresponding instantaneous Reynolds numbers. On the other hand, at higher pulsation frequencies, considerable deviation from the pure sinusoidal signal occurs for the development length and its amplitude decreases with increase in pulsation frequency. Finally, using the results of the present simulations, a simple correlation is proposed that can be used in order to predict the maximum _ A and F. development during a cycle as functions of ReM ; m  2012 Elsevier Inc. All rights reserved.

1. Introduction and aim of work Pulsating flows are encountered in many engineering applications and also in Nature and biological systems. Owing to their importance in diverse fields, such flows have been dealt with in considerable detail over the past several decades. For fundamental investigations of such pulsations in time, it is also a common practice to introduce sinusoidal time variations in the driving pressure gradient and assume the flow to be hydrodynamically fully-developed. These type of flows may be identified as ‘‘pressure gradientdriven’’ flows. To the best of our knowledge, (Sexl, 1930) was the first to provide an analytical solution for such problems. Similar analytical studies were later presented by Szymanski (1932), Lambossy (1952) and Womersley (1955, 1957) in order to study various aspects of the flow. More complex time variations in the pressure gradient pulsation were subsequently considered by Uchida (1956), who showed that if the prescribed periodic signal can be expressed in terms of a Fourier series, in which each of the terms represents a harmonic pressure variation, using the linear properties of the governing conservation equation and the boundary conditions, the final solution can be constructed by a simple ⇑ Corresponding author. E-mail addresses: [email protected] (S. Ray), [email protected]. gov.tr (B. Ünsal), [email protected] (F. Durst). 0142-727X/$ - see front matter  2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.ijheatfluidflow.2012.06.001

linear superposition. This leads to the very important conclusion that each of these terms in the Fourier series can then be studied on their own. It may, therefore, be further concluded that the solution of arbitrary ‘‘pressure gradient-driven’’, laminar, fully-developed, time-periodic flows is completely known and the results are available in the literature. Several years later, (Ray et al., 2005) presented the solution for sinusoidally pulsating, laminar pipe flows in suitable dimensionless form and recognized the existence of a unique signature map that relates the dependence of the ratio of dimensionless amplitudes of the mass flow rate and the pressure gradient signals and the corresponding phase lag as functions of dimensionless pulsation frequency. They also showed that using the signature map, complex ‘‘mass flow rate-controlled’’ flows, involving higher order harmonics, can also be solved and the same signature map can be utilized for reconstruction of the instantaneous mass flow rate from the measured (and hence known) pressure gradient time signals and vice versa. In order to verify the findings of Ray et al. (2005), Ünsal et al. (2005) performed detailed experimental measurements on the velocity distribution for the ‘‘mass flow rate-driven’’, sinusoidally pulsating, pipe flows and observed excellent agreement with the theory. Encouraged by these findings, (Ünsal et al., 2006; Durst et al., 2007) adopted the unique signature map for the measurement of the mass flow rate as a function of time from the instantaneously

168

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

Nomenclature CL D f, F i J0, J1 L Lz ðL=DÞA _ m _ m; p, p⁄ ^ p r, r⁄ R ReA Recrit ReR ReX StR t, t⁄ t p

asymptotic value of ðL=DÞA at F ? 1 pipe diameter (m) dimensional (Hz) and dimensionless frequency, R2f/m pffiffiffiffiffiffiffi imaginary unit 1 zeroth and first order Bessel functions of first kind development length (m) dimensionless length of the computational domain normalized amplitude of L/D, (L/D)A/(L/D)M dimensional (kg/s) and dimensionless mass flow rate, _ m _M m= dimensional (Pa) and dimensionless effective pressure, p=qu2M modified pressure gradient dimensional (m) and dimensionless radial coordinate, r/ R pipe radius (m) normalized amplitude of Reynolds number, ReA/ReM critical Reynolds number Reynolds number based on pipe radius, UMR/m Reynolds number based on pipe diameter, UXD/m, X = A, M Strouhal number, fR/uM dimensional (s) and dimensionless time, tuM/R dimensionless time period, 1/StR

recorded pressure gradient data in the hydrodynamically fullydeveloped section of pipes. While (Ünsal et al., 2006) used Diesel oil and water as working fluids, (Durst et al., 2007) conducted their study with a ‘‘mass flow rate-controlled’’ air flow. Earlier, (Womersley, 1955) also suggested a similar measurement technique for the flow of blood through arteries; however, he did not recognize the existence of the signature map in dimensionless form. Although the success of these studies established the use of the signature map of Ray et al. (2005) beyond doubt, drawbacks of the method, however, lie in the assumption of incompressibility1 and the existence of fully-developed flow. As the first problem cannot be avoided in real-life physical situations and hence the method has to be used with restrictions for compressible fluids, the major concern, therefore, boils down to the evaluation of the axial length required for the hydrodynamically fully-developed flow to be established. In the present study, an attempt has been made to address this issue for moderate and high Reynolds number laminar flow through pipes under sinusoidal mass flow rate variations in time. Concerning the steady flow through pipes and channels, many researchers have dealt with the hydrodynamic development length of laminar flows in the past. An excellent overview of these studies has been presented by Durst et al. (2005). Most of the studies on the requirement of the development length for steady flows were conducted either theoretically or experimentally. As pointed out by Durst et al. (2005), theoretical studies were unable to provide correct solutions owing to the presence of non-linear convective terms in the momentum equations. These terms were either neglected or linearized using grossly simplifying assumptions in most of the theoretical investigations. Experimental studies, on the other hand, were also unable to produce satisfactory results with sufficient accuracy due to limitations in the employed test set-up. In addition, most of the measurement methods yielded enormous inaccuracies when identifying the development length, owing to the asymptotic nature of the physical phenomenon. It may, therefore, be concluded 1 Which fails for very high-frequency pulsations of compressible fluids such as air, where the pressure waves distort the original pressure gradient signals.

u, u⁄ z, z⁄

dimensional (m/s) and dimensionless velocity, u/uM dimensional (m) and dimensionless axial coordinate, z/R

Greek symbols Dh phase difference (rad) m kinematic viscosity (m2/s) q density (kg/m3) s dimensionless time, mt/R2 X angular frequency, 2pF w complex function Subscripts A amplitude in at inlet m pertaining to mass M cycle averaged mean max maximum p period r in radial direction u pertaining to velocity z in axial direction Superscript ⁄ dimensionless quantity

that only numerical simulations can provide sufficiently accurate results for developing flows through pipes and channels that eventually determine the axial length required for the establishment of fully-developed flows. For steady pipe and channel flows, such numerical investigations have already been carried out and the results have been reported by Durst et al. (2005). Particularly for pipe flows, their correlation for development length is

L ¼ ½ð0:619Þ1:6 þ ð0:0567ReÞ1:6 1=1:6 D

ð1Þ

where L is the required development length, D is the pipe diameter and Re is the Reynolds number based on the diameter of the pipe. Compared with the theoretical and experimental investigations on the development length of steady pipe flows, there have been very few studies on similar aspects of pulsating flows. Several decades ago, (Atabek and Chang, 1961) developed an analytical solution by linearizing the axial momentum equation, with an assumption that the convective velocities, appearing in the inertia terms, remain the same at all instants in time for the entire flow and are equal to the velocity at the inlet. They further assumed a plug-flow profile for the velocity. Based on these grossly simplifying assumptions, they were able to solve only for sinusoidal pulsations without any flow reversal. The second analytical solution of the problem was proposed by Caro et al. (1978). They derived their final solution by neglecting the interaction between the oscillatory and the steady parts of the flow and, in addition, by using a boundary layer-type approximation. It is worthwhile and also surprising to note that in spite of obtaining a simplified analytical solution, (Caro et al., 1978) did not provide any information regarding the development length of pulsating flow. Experimental investigations on both the developing and pulsating flows were also performed by Atabek et al. (1964), Florio and Mueller (1968) and Gerrard and Hughes (1971). Their investigations, however, were unable to provide any conclusive information on the requirement of the development length since the asymptotic nature of the flow towards its hydrodynamic fully-developed state always creates experimental difficulties.

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

To the best of authors’ knowledge, there has been no other numerical study on developing pulsating flow apart from the investigation of Krijger et al. (1991), who numerically solved the governing equations for sinusoidally pulsating channel flow. They performed computations up to 16 Hz by varying the dimensionless pulsation amplitude from 0.25 to 0.75. They concluded that for frequencies lower than 1 Hz, the flow is quasi-steady and hence the maximum required inlet length to ensure established fully-developed flow can be calculated based on the maximum Reynolds number during the cycle and the assumption of steady flow. On the other hand, for higher pulsation frequencies f > 10 Hz, they argued that since the flow is inertia dominated where the development length does not vary strongly with the pulsation frequency, it can be estimated according to the development length for the cycle-averaged, mean Reynolds number2 under steady flow conditions. The brief literature review presented here clearly shows that detailed investigations of the development length of sinusoidally pulsating, laminar pipe flows are not available in the literature, although such information is of the utmost importance in order to use the unique signature map, proposed by Ray et al. (2005), for the measurement of the instantaneous mass flow rates through tubes from the recorded pressure gradient time signals in the fullydeveloped regime. This justifies the present research effort to carry out detailed numerical investigations on this topic. In the present study, therefore, extensive numerical simulations were carried out in order to determine the development length of sinusoidally pulsating, laminar pipe flows as functions of pulsation frequency, amplitude of the mass flow rate pulsation and the mean Reynolds number in the moderate and high Reynolds number regimes. The present paper is organized in the following manner. Section 2 describes the mathematical formulation, where the dimensionless governing equations and the boundary and the initial conditions are presented in Section 2.1; the analytical solution for the laminar, fully-developed, sinusoidally pulsating pipe flow is briefly reported in Section 2.2 and in Section 2.3, the numerical method and the calculation procedure for the required development length are described. The results of the numerical simulations are presented and discussed in Section 3. The correlation, derived from the simulated results, for the evaluation of the amplitude of the development length (or its maximum value during a cycle) is also presented in Section 3. Finally, the conclusions from the present study are reported in Section 4.

2. Mathematical formulation 2.1. Governing equations, boundary and initial conditions The geometry of the present problem is fairly simple as it consists of a circular pipe of radius R (diameter D = 2R) and length Lz. At the inlet of the pipe, an incompressible fluid, with a density q, enters with a plug-flow velocity profile that has been assumed to vary in time as follows;

uz;in ðtÞ ¼ uM þ uA sinð2pftÞ ¼ uM ½1 þ uA sinð2pF sÞ

ð2Þ

where uM is the cycle-averaged and also cross-section-averaged mean velocity through the pipe, uA is the amplitude of velocity pulsation and f is the pulsation frequency. Further, uA ¼ uA =uM is the dimensionless value of uA, F = R2f/m is the dimensionless pulsation frequency and s = mt/R2 is the dimensionless time, defined in a manner same as that introduced by Ray et al. (2005). 2 Henceforth this will be referred to as the ‘‘mean’’ Reynolds number throughout the text.

169

Assuming that all the thermophysical properties remain constant in space and time (which is also a requirement for the existence of the fully-developed flow), and the flow to be axisymmetric (i.e., there is no circumferential velocity), the following dimensionless variables have been defined in order to represent the governing mass and momentum conservation equations for laminar flow in non-dimensional form:

z r ; r ¼ R R uz ur ; ur ¼ uz ¼ uM uM p p ¼ 2 quM tuM t ¼ R

z ¼

ð3aÞ ð3bÞ ð3cÞ ð3dÞ

where z and r are the axial and the radial coordinates, respectively, uz and ur are the axial and the radial velocity components, respectively, and p is the effective pressure that also includes the hydrostatic pressure variation. With these definitions, the dimensionless conservation equations have been written as follows (Bird et al., 2002): Continuity equation:

1 @     @uz r ur þ  ¼ 0 r @r  @z

ð4Þ

Axial momentum equation:

" #      @uz @p 1 1 @ @ 2 uz  @uz  @uz  @uz þ þ u þ u ¼  þ r r z @t @r  @z @z ReR r  @r @r  @z2

ð5Þ

Radial momentum equation:  @ur @u @p  @ur þ uz r ¼    þ ur  @t @r @z @r

" #   1 @ 1 @    @ 2 ur þ r ur þ 2 ReR @r  r  @r @z

ð6Þ

where ReR is the Reynolds number based on the cycle-averaged and also cross-section-averaged mean velocity and the radius of the pipe, defined as

ReR ¼

uM R

m

¼

ReM 2

ð7Þ

where ReM is the conventional mean Reynolds number, defined on the basis of the pipe diameter. It is evident from the governing Eqs. (4)–(6) that the solution for the velocity distribution depends on the cycle-averaged mean Reynolds number. It may be mentioned here that the present analysis has been carried out under the assumption of laminar flow and the existence of axi-symmetric condition, which require appropriate justification. As far as the steady pipe flows are concerned, the linear stability analysis shows that pipe flows are unconditionally stable and they can remain laminar irrespective of the Reynolds number (see for example, Lessen et al., 1968; Davey and Drazin, 1969; Salwen and Grosch, 1972; Salwen et al., 1980; Schmid and Henningson, 1994). Under the experimental conditions, however, Osborne (Reynolds, 1883) showed that laminar flow is generally expected up to Re  2000. Later more sophisticated and detailed measurements revealed that pipe flows are globally stable up to Re  1800 (Peixinho and Mullin, 2006). In order to resolve the controversy and to address the discrepancy between the theoretical analysis and experimental observations, several attempts have been made since then (Asen, 2007). Notable among them deal with the stability analysis of flows inside corrugated tubes, where the corrugations represent the surface roughness (see for example, Cotrell et al., 2008; Loh and Blackburn,

170

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

2010). These studies clearly showed that the critical Reynolds number (Recrit), below which the pipe flow remains laminar, is strongly dependent on the height of the roughness and Recrit ? 1 for smooth pipes. This observation leads to a conclusion that if the test set-up is prepared with precision (i.e., with a pipe as smooth as possible) and if the experiments are conducted carefully (i.e., by avoiding external disturbances), laminar flow can be sustained up to a very high Reynolds number (Re  105 was reported by Pfenninger (1961)). For pulsating flows, however, the stability analysis is more complex owing to the presence of two additional parameters in the form of pulsation frequency and amplitude. (Yang and Yih, 1977) showed that single harmonic oscillatory flows are linearly stable for all parameters. The analysis was later extended by Nebauer and Blackburn (2009) for pulsating flows and a similar observation was reported. These conclusions, however, contradict the experimental observation on transition of pulsatile pipe flows (Stettler and Hussain, 1986). Several investigations were carried out and are still being conducted in order to find out the reason for apparent contradiction (see for example, Tozzi and von Kerczek, 1986; Straatman et al., 2002; Fedele et al., 2005; Nebauer and Blackburn, 2010; Smith and Blackburn, 2010). Although these detailed studies included non-linear stability analysis and also considered the effect of non-axi-symmetric disturbances, exact quantification of Recrit still remains largely inclusive. Nevertheless, the careful experimental study, conducted by Ünsal et al. (2005) for sinusoidally pulsating pipe flows, showed that laminar flow could be maintained at least up to a maximum Reynolds number of 4,650. As will be shortly apparent in Section 3, the maximum cycle-averaged Reynolds number in the present study has been taken as 2,000 and the maximum dimensionless pulsation amplitude for _ A , defined later in Eq. (10a)) has been considered mass flow rate (m to be 0.8. Hence the maximum Reynolds number during the cycle is only 3,600. Therefore, it can be authentically claimed that the present assumption of laminar flow is fully supported not only by the results of linear stability analysis, but also by the detailed experiments, conducted by the present authors. Another pertinent question that may be raised is regarding the validity of the assumption of axi-symmetric flow. This assumption, however, cannot be justified by any theoretical analysis. Nevertheless, the experimental measurements of the instantaneous velocity profiles, conducted by Ünsal et al. (2005), showed that the flow, under all tested conditions, remained indeed axi-symmetric. This substantiates the present assumption. The boundary conditions used for the present simulations are as follows: ⁄

uz;in ¼ 1 þ uA sinð2pSt R t  Þ uM

ð8Þ

where StR is the Strouhal number, defined as

StR ¼

fR F ¼ uM ReR

_ m _ A sinð2pF sÞ ¼1þm _M m

ReðsÞ ¼ ReM þ ReA sinð2pF sÞ ¼ ReM 1 þ ReA ð2pF sÞ _  ðsÞ ¼ m

ð10aÞ ð10bÞ

_ M ¼ quM pR2 is the cycle-averaged mass flow rate through where m _ A ¼ m _ A =m _ M ¼ uA is the dimensionless amplitude of the the pipe, m mass flow rate pulsation, ReM is given in Eq. (7) and ReA ¼ ReA =ReM is the normalized Reynolds number amplitude, _ A . which is also equal to m   2. On the solid wall (r⁄ = 1, "z⁄), no slip uz ¼ 0 and impermeable conditions ur ¼ 0 have been imposed. 3. On the line of symmetry ðr ¼ 0; 8z Þ; @uz =@r  ¼ 0 and ur ¼ 0 (implying no mass flux through the line of symmetry) have been set.   4. At the exit of the computational domain z ¼ Lz ; 8r  , zero diffusion condition for both the velocity components, i.e., @ 2 /=@z2 ¼ 0 for / ¼ uz and ur , has been assumed. It may be noted that this condition is expected to be less restrictive than the fully-developed condition in the presence of back flow that is generally expected at higher pulsation frequencies. As far as the initial condition is concerned, the velocity solution for the steady developing flow through a pipe of exactly the same geometric dimensions, with a uniform inlet velocity condition uz;in ¼ 1 and for the same mean Reynolds number has been considered to be the solution at t⁄ = 0. Such solutions for a given ReM have been obtained numerically using exactly the same grid distribution in the computational domain as employed later for the transient problem. 2.2. Analytical solution for fully-developed flow

⁄

1. At the inlet (z = 0, "r ), a periodic plug-flow profile in time has been assumed for the axial velocity which is given by Eq. (2). In the dimensionless form, it can be rewritten as

uz;in ðt Þ ¼

icant variations in the experimental measurement of steady development length for laminar pipe flows (Durst et al., 2005). It is, however, known that a partially developed inlet velocity profile requires shorter development length as compared to that for an uniform inlet condition. Owing mainly to this reason, in the present study, the inlet velocity profile has been assumed in the form given by Eq. (8) in order to obtain the most conservative estimate on the required development length. This boundary condition clearly implies that the velocity solution also depends on the dimensionless pulsation frequency F = StRReR (according to Eq. (9), which is a combination of the Strouhal number and the mean Reynolds number) and the amplitude of pulsation for the velocity uA . Using the definitions of the mass flow rate through the tube and the instantaneous Reynolds number, one may also rewrite Eq. (8) as follows:

ð9Þ

It is obvious that according to the present definitions, StRt⁄ = Fs. Further, for a given set of F and ReR, the dimensionless time period can be calculated as t p ¼ 1=StR ¼ ReR =F. In addition, the radial velocity component at the inlet has been set to zero. The inlet boundary condition, used here, may appear to be somewhat artificial since under experimental conditions it would be very difficult to maintain a plug-flow velocity profile at the inlet of a pipe. Since the extent to which the flow already develops before it enters the test section could not be quantified, there were signif-

The detailed analytical solution for laminar, fully-developed, sinusoidally pulsating pipe flow has been already reported in the dimensionless form by Ray et al. (2005) and Ünsal et al. (2005). In this section, the same solution is briefly presented not only for the sake of completeness, but also, as it will shortly be apparent, because this solution has been invoked for comparison with the numerical solution at different axial locations in order to determine the required development length as a function of time for a _ A and ReM. prescribed set of F; m The axial pressure gradient required to drive a mass flow rate through a tube that is given by Eq. (10a) may be written as



 _ 1 @p m ^M 1 þ A sinð2pF s þ Dhm Þ ¼p q @z jwm j

ð11Þ

where, from the steady, fully-developed flow consideration, ^M ¼ 8lm _ M =pq2 R4 can be shown as the modified (in the sense that p it is divided by the density) pressure gradient required to drive a _ M through the pipe. mass flow rate m

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

The same equation may also be written in the dimensionless form as follows:



 _ A @p 16 m ¼ 1 þ p F s þ D h Þ sinð2 m @z ReM jwm j

ð12Þ

Through analytical considerations, Ray et al., 2005 showed that for sinusoidally pulsating, laminar pipe flows, it is possible to define a complex function for the mass flow rate in the following form that depends solely on the dimensionless pulsation frequency:

wm ðFÞ ¼

8

X

"

3=2

1

2J 1 ðX1=2 i

Þ

#

X1=2 i3=2 J0 ðX1=2 i3=2 Þ

¼ Rðwm Þ þ iIðwm Þ

ð13Þ

where X = 2pF is the dimensionless rotational frequency, R and I are the real and the imaginary parts, respectively, i is the square root of unity and J0 and J1 are the zeroth and the first-order Bessel functions of the first kind, respectively. This function forms the basis of the unique signature map proposed by Ray et al. (2005). In fact, jwmj and Dhm in Eqs. (11) and (12) are the ratio of the dimensionless amplitudes and the phase difference of the mass flow rate and the pressure gradient time signals, which can be derived from the function wm(F) as

jwm j ¼ ½fRðwm Þg2 þ fIðwm Þg2 1=2  Rðwm Þ Dhm ¼ tan1 Iðwm Þ

ð14aÞ ð14bÞ

Once the required pressure gradient along with the amplitude ratio and the phase lag are known, the solution for the axial velocity in the fully-developed section may be obtained as

uz ¼ 2ð1  r2 Þ þ

_ A m jw j sinð2pF s þ Dhm  Dhu Þ jwm j u

ð15Þ

where, similar to the mass flow rate, a complex function for velocity wu(r⁄, F) may be defined in the following manner, which is a function of both r⁄ and F:

3 3=2 J 0 X1=2 i r  5 ¼ Rðwu Þ þ iIðwu Þ wu ðr  ; FÞ ¼ 1 3=2 X J 0 ðX1=2 i Þ 2 84

ð16Þ

On the basis of Eq. (16), jwuj and Dhu, appearing in Eq. (15) are calculated as

jwu j ¼ ½fRðwu Þg2 þ fIðwu Þg2 1=2  Rðwu Þ Dhu ¼ tan1 Iðwu Þ

ð17aÞ ð17bÞ

171

Both the convective and diffusive terms in the momentum equations Eqs. (5) and (6) were discretized using the second-order accurate central differencing scheme. Use of such a scheme for the convective terms is, however, prone to numerical instabilities, leading to divergence of the solution during the iterative calculation procedure. In order to overcome this problem, the coefficients of the discretized equations were formed using the more stable and the first-order accurate upwind differencing scheme, while the second-order accuracy of the central differencing scheme was restored by adding an appropriate source (correction) term, following the deferred correction approach proposed by Khosla and Rubin (1974). The transient terms in Eqs. (5) and (6) were discretized using the second-order accurate, three-time-level, implicit scheme (see, for example, Ferziger and Peric´, 1999 for details) owing to the requirement for a smaller number of time steps per cycle for attaining the time-step independent results. Finally, the discretized algebraic equations, resulting from a given conservation equation, were solved using the Strongly Implicit Procedure (SIP), proposed by Stone (1968). The entire computational domain, described in Section 2.1, was divided into Nr  Nz number of non-overlapping control volumes. For all the computations in the present study, a uniform grid layout was chosen in the radial direction, whereas the numerical grid in the axial direction was allowed to expand with a common ratio of 1.05, in order to accommodate the rapid changes in the velocity components that are inevitable close to the inlet of the pipe. Before carrying out the final numerical simulations, a careful grid- and time step-independent study was carried out. The study revealed that 80  400 control volumes are required for the spatial resolution (which is the same as that used by Durst et al. (2005) for their study), whereas 180 time steps per time cycle are necessary for the temporal discretization. The number of time cycles required to obtain repetitive results was found to increase considerably with increase in pulsation frequency. From the converged solution after each time step, the development length was calculated by comparing the axial velocity solutions at different axial locations with the analytical solution for the fully-developed pulsating flow, given by Eq. (15). It was observed that the center-line velocity requires the longest distance to attain the periodic fully-developed state and hence, at a given time, the axial distance from the inlet, where the center-line velocity attains 99% of the analytical solution in Eq. (15), was considered to be the development length that ensures the establishment of the fully-developed velocity profile in the entire cross-section of the pipe. This observation is similar to that mentioned also by Durst et al. (2005). As expected, the development length was found to be a function of time, that also depends upon the pulsation frequency, the amplitude of mass flow rate and the mean Reynolds number.

2.3. Numerical simulations and calculation of development length The Computational Fluid Dynamics (CFD) code, used for the present study, is the same as that used by Durst et al. (2005) for the calculation of the development length for pipe and channel flows under steady-state operating conditions, a study which eventually led to the correlation given in Eq. (1). The dimensionless governing equations Eqs. (4)–(6), were discretized using a control volume-based finite difference method on a collocated grid layout. In the present approach, the cell-face velocities were calculated using the momentum interpolation method, suggested by Rhie and Chow (1983), in order to avoid pressure velocity decoupling. These cell-face velocities and also the nodal velocities were corrected using the nodal pressure corrections. The equations for the nodal pressure corrections were derived from the mass conservation equation Eq. (4) following the SIMPLE algorithm, proposed originally by Patankar and Spalding (1972) and later well documented by Patankar (1980).

3. Results and discussion In the present study, numerical simulations were carried out in order to determine the development length for sinusoidally pulsating, laminar pipe flows in the moderate and high Reynolds number regimes, with ReM ranging from 100 to 2000. For each of these cases, _ A ¼ 0:2, three different dimensionless mass flow rate amplitudes, m 0.4 and 0.8, were considered and the dimensionless pulsation frequency F was varied over a wide range, from 0.1 to 20, covering the quasi-steady to the inertia-dominated regimes in pulsating flows. In this section, the development of the axial velocity profiles is first presented, followed by the determination of the required development length and its quantification through correlation. The axial velocity profiles at different axial locations 1 6 z/ D 6 100 are presented in Fig. 1 for ReM = 1000, F = 10 and _ A ¼ 0:8, where parts (a), (b), (c) and (d) correspond to the different m

172

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

_ A ¼ 0:8. Lines represent numerical simulations at different z/D Fig. 1. Development of axial velocity profiles at selected phases during the cycle for ReM = 1000, F = 10 and m and symbols represent analytical solutions in the fully-developed section, according to Eq. (15).

selected phases / = 0, p/2, p and 3p/2, respectively, during the cycle. In these figures, the analytical solutions for the fully-developed flow at corresponding phases are also presented with different symbols for comparison. It is evident that at least for the given flow conditions, the axial velocity profiles for z/D = 100 at different phases are already in the fully-developed regime. The profiles for z/D = 1 (i.e., very close to the inlet) in these plots, on the other hand, are similar to those observed by Durst et al. (2005) for the steady, developing flow through pipes, with maxima appearing near the wall region. Further downstream, however, under the effects of accelerating and decelerating pressure fields required to maintain the prescribed mass flow rate through the pipe, the development of the axial velocity profiles is observed to be significantly different to the steady flow situation. In Fig. 2, the instantaneous variations of the development length over two representative time periods, as functions of the phase during the cycle, are presented for two different Reynolds numbers, ReM = 100 and 2000, and for two different amplitudes _ A ¼ 0:2 and 0.8. Variations in the developof the mass flow rate, m ment length, calculated according to Eq. (1) for the steady flow, corresponding to Re (s) given by Eq. (10b) for different phases, are also presented in Fig. 2 with different symbols for comparison. As Fig. 2 indicates, L/D increases with increase in ReM for a given set _ A and F, because the mean development length increases with of m increase in the mean Reynolds number. Also, as expected, the maximum development length increases, whereas its minimum value decreases during the cycle with increase in the mass flow rate amplitude for a given combination of ReM and F.

It is, however, interesting to note that for specified values of ReM _ A , the instantaneous variation of L/D with the phase and m (or dimensionless time) decreases considerably with increase in pulsation frequency. For very low frequency pulsations, e.g., for F = 0.1, L/D nearly follows the prediction according to Eq. (1) for the steady flow with instantaneous Re(s), which is expected in the quasi-steady regime. With subsequent increase in the pulsation frequency, deviations are apparent and are fairly prominent for very high-frequency pulsations (e.g., for F = 10 or 20). By observing the variations in L/D in Fig. 2, one may now define the amplitude of L/D as follows:

ðL=DÞA ¼ ðL=DÞmax  ðL=DÞM

ð18Þ

where (L/D)max is the maximum development length during a cycle and (L/D)M is the mean development length calculated from Eq. (1) for a given ReM. It is important to note that (L/D)A in Eq. (18) has been defined on the basis of the maximum L/D during a cycle since from an engineering point of view, a designer would be interested in the knowledge of the maximum development length for pulsating flow rather than its minimum value. Further, one may also define a phase lag (say, D/), which is given by the difference between the phases where the maximum development lengths occur according to Eq. (1) for instantaneous Re(s) (i.e., for truly quasi-steady flow or for extremely low-frequency pulsations) and for a given set of _ A and F (i.e., for a given physical flow situation). ReM ; m In view of these definitions, it is now apparent from Fig. 2 that the amplitude of the development length decreases, whereas the corresponding phase lag increases, with increase in the pulsation

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

173

Fig. 2. Variations of instantaneous L/D with phase during the cycle. Lines represent numerical simulations for different F and symbols represent L/D with Re(s), according to Eq. (1).

frequency. For very high-frequency pulsations, D/ asymptotically approaches a value of p, whereas (L/D)A tends to become insensitive to the increase in F. Other than these, it is also very important to note that with the increase in the pulsation frequency, particu_ A , the instantaneous variations in L/D larly for higher values of m deviate considerably from the expected signals for the pure sinusoidal pulsation. This can be readily observed from the plots in Fig. 2 for frequencies varying between 2 and 20, irrespective of the value of ReM. For this reason, it may be recognized that the conventional amplitude and phase difference analysis with respect to the pulsation frequency cannot be performed in a straightforward manner, although the imposed mass flow rate through the tube always remains sinusoidal in time. Therefore, in the subsequent part of this study, variations of (L/D)A were analyzed disregarding the changes in D/, since (L/D)A provides the measure of the maximum required development length during a cycle for the pulsating flow, which is useful for design purposes. One may now have a closer look at Fig. 2 and compare the results for the same dimensionless amplitudes of mass flow rate. This reveals that although the instantaneous variations in L/D differ quantitatively by significant magnitudes, due to difference in ReM, they indeed appear to be qualitatively similar. It would therefore be worthwhile to study the behavior of L/D, normalized with respect to (L/D)M, i.e., the variations of (L/D)⁄. Fig. 3 shows the instantaneous variations of (L/D)⁄ for three different mass flow rate amplitudes and for different pulsation frequencies. The results for ReM = 2000 are shown with different types of lines, whereas the variations for ReM = 100 are exhibited with symbols. It is very interesting that with the applied normalization, the instantaneous

variations of (L/D)⁄ appear to be independent of ReM and depend _ A and F. only upon m In view of this important observation and the great simplification, it is now possible to examine the variations of the normalized amplitude of development length that is defined as follows:

ðL=DÞA ¼

ðL=DÞA ¼ ðL=DÞmax  1 ðL=DÞM

ð19Þ

_ A as functions of F The variations of ðL=DÞA for different values of m are presented in Fig. 4, where the symbols represent the results of the numerical simulations and the lines were obtained from the predictive correlation, given by Eq. (20). Some additional results, obtained for ReM = 2000, are also presented. These additional simu_ A ¼ 0:2; 0.4 and 0.8 with F = 3, 4 and 7 lations were carried out for m _ A ¼ 0:1, 0.2, 0.4, (shown with solid symbols), and for F = 50 with m _ A ¼ 0:1, 0.5 and 1 are shown with empty 0.5, 0.8 and 1 (results for m circles). Fig. 4 clearly shows the existence of two asymptotes for ðL=DÞA – one for the low-frequency quasi-steady regime, where ðL=DÞA tends _ A , and the other for the high-frequency inertia-dominated reto m gime, where the asymptotic value depends upon the dimensionless amplitude of mass flow rate. It is further evident from the figure that as far as the variations of ðL=DÞA are concerned, the quasi-steady regime exists for F 6 0.5, whereas beyond F = 20, the inertiadominated regime is observed. Finally, from the present study, the correlation for ðL=DÞA is obtained in the following form:

ðL=DÞA ¼ C L þ

_ A  C L m 1 þ ðF=3:25Þ3

ð20Þ

174

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

Fig. 3. Instantaneous variations of normalized L/D. Lines represent cases for ReM = 2000 and symbols represent cases for ReM = 100.

can be calculated using Eqs. (20) and (21). Further, from the relationship given in Eq. (19), the maximum normalized development length can be obtained. Finally, by calculating (L/D)M from Eq. (1) using the known mean Reynolds number, the maximum value of L/D during the cycle can be evaluated from Eq. (19). In other words, using Eqs. (1), (19), (20) and (21), the maximum development length can be calculated from the following relationship; " # h i1=1:6 _ A  C L m 1:6 1:6 ðL=DÞmax ¼ ð0:619Þ þ ð0:0567ReM Þ 1 þ CL þ 1 þ ðF=3:25Þ3

ð22Þ

Fig. 4. Variations of normalized amplitude of development length with pulsation frequency. Symbols are for computed data and the lines are obtained from Eq. (20).

_ A . Eq. (20) clearly shows that as where CL is a function of m   _ F ! 0; ðL=DÞA tends to mA , whereas for F ! 1; ðL=DÞA assumes a value equal to CL. Hence, CL is the asymptotic value of ðL=DÞA in the inertia-dominated regime, which is found to be a nearly linear _ A and is correlated as follows; function of m

  1:1881 _A C L ¼ 0:1671 m

ð21Þ

From the knowledge of the pulsation frequency and the amplitude of mass flow rate, the normalized amplitude of development length

As explained previously, owing to the considerable deviations of the instantaneous development length from the true sinusoidal signal, particularly for higher amplitudes of mass flow rate and pulsation frequencies, it is impossible to study these variations using the conventional amplitude and phase difference analysis. Nevertheless, the most relevant question regarding the maximum length of the pipe required in order to ensure the establishment of the hydrodynamically fully-developed flow condition can be readily answered from Eq. (22). 4. Conclusions and final remarks In the present study, extensive numerical simulations were carried out for the laminar, developing flow through pipes under sinusoidally pulsating mass flow rate variations in time. All the computations were performed in the moderate and high Reynolds number regimes (100 6 ReM 6 2000). In general, three different

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

175

amplitudes of mass flow rate were considered and for each of these cases, the pulsation frequency has been varied from the quasi-steady regime (F = 0.1) to the inertia-dominated regime (F = 20). As a result of the present investigation, the following conclusions can be drawn;

cro-channel flows. This study would result in a working correlation for the normalized amplitude of the development length during a cycle, which is expected also to be a function of the mean Reynolds number, in addition to the amplitude of mass flow rate and the pulsation frequency.

1. The development of the center-line velocity dictates the requirement for an established fully-developed flow condition over the entire cross-section of the pipe since the center-line velocity requires the longest distance from the entrance in order to attain the periodic fully-developed state. The axial velocity profiles close to the inlet of the pipe appear to be similar to those observed for the steady, developing flow with maxima occurring near the pipe wall. Far downstream, however, the velocity profiles are significantly different to the steady flow situation. 2. In the low-frequency, quasi-steady regime, the variations in the instantaneous development length are truly sinusoidal and follow the predictions based on the steady flow condition with the instantaneous Reynolds number. This observation is similar to that reported by Krijger et al. (1991). With increase in pulsation frequency, the amplitude of the development length decreases and the corresponding phase lag increases. Therefore, the maximum development length, calculated on the basis of the maximum Reynolds number during the cycle and the steady flow assumption, provides the most conservative estimate. 3. At higher pulsation frequencies, particularly for higher amplitudes of mass flow rate, the instantaneous development length differs significantly from the pure sinusoidal signal. As a consequence, the variations in L/D with time cannot be represented using the conventional amplitude and phase lag analysis. 4. Although the variations of the development length for different mean Reynolds number differ considerably, its normalized variations are independent of ReM and depend only on the dimensionless amplitudes of mass flow rate and the pulsation frequency. 5. The correlation for the normalized amplitude of the development length was obtained from the present study, which can be used for the prediction of the maximum development length during the cycle. The asymptotic value of the normalized amplitude for the quasi-steady regime equals the dimensionless amplitude of the mass flow rate. For the inertia-dominated regime, on the other hand, the asymptotic value depends upon the mass flow rate amplitudes, hence the mean development length, calculated based on the mean Reynolds number, cannot be used in order to estimate the maximum development length during the cycle. This observation quite significant since it differs considerably from that reported by Krijger et al. (1991).

References

As a final remark, it may be mentioned that in a real-life situation, pulsating flows are generally more complex in nature and hence do not follow sinusoidal mass flow rate variations. More complex pulsations, e.g., triangular, trapezoidal, saw-tooth or square-wave type pulsations, contain higher order harmonics which are multiples of the base frequency. Although the problem of developing flow through pipes is non-linear, owing to the presence of the convective terms in the momentum equations, since the maximum development length decreases with increase in the pulsation frequency, the higher order harmonics are expected to develop faster than the base frequency. Therefore, the present predictions and the correlation for sinusoidal pulsations can be recommended as a conservative estimate for complex pulse shapes. Finally, it may be mentioned that the present study covers only the moderate and high Reynolds number regimes of laminar flow. It would, therefore, be worthwhile to extend the investigation to the low Reynolds number regime, applicable particularly for mi-

Asen, P.-O., 2007. Stability of plane Couette Flow and pipe Poiseuille Flow. PhD thesis, Kungliga Tekniska hgskolan (KTH) Stockholm, Sweden. Atabek, H.B., Chang, C.C., 1961. Oscillatory flow near the entry of a circular tube. Z. Angew. Math. Phys. 12, 185–201. Atabek, H.B., Chang, C.C., Fingerson, L.M., 1964. Measurement of laminar oscillatory flow in the inlet length of a circular tube. Phys. Med. Biol. 9, 219–228. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2002. Transport Phenomena. John Wiley & Sons, New York. Caro, C.G., Pedley, T.J., Schroter, R.C., Seed, W.A., 1978. The Mechanics of the Circulation. Oxford University Press, Oxford. Cotrell, D., McFadden, G., Alder, B., 2008. Instability in pipe flow. Proc. Nat. Acad. Sci. USA 105, 428–430. Davey, A., Drazin, P., 1969. The stability of poiseuille flow in a pipe. J. Fluid Mech. 36, 209–218. Durst, F., Ünsal, B., Ray, S., Saleh, O., 2005. The development lengths of laminar pipe and channel flows. J. Fluids Eng. 127, 1154–1160. Durst, F., Ünsal, B., Ray, S., Trimis, D., 2007. Method for defined mass flow variations in time and its application to test a mass flow rate meter for pulsating flows. Meas. Sci. Technol. 18, 790–802. Fedele, F., Hitt, D., Prabhu, R., 2005. Revisiting the stability of pulsatile pipe flow. Eur. J. Mech. B/Fluids 24, 237–254. Ferziger, J.H., Peric´, M., 1999. Computational Methods for Fluid Dynamics. Springer, Berlin. Florio, P.J., Mueller, W.K., 1968. Development of a periodic flow in a rigid tube. Trans. ASME J. Basic Eng. 90, 395–403. Gerrard, J.H., Hughes, M.D., 1971. The flow due to an oscillating piston in a cylindrical tube: a comparison between experiment and a simple entrance flow theory. J. Fluid Mech. 50, 97–106. Khosla, P.K., Rubin, S.G., 1974. A diagonally dominant second order accurate implicit scheme. Comput. Fluids 2, 207–209. Krijger, J.K.B., Hillen, B., Hoogstraten, H.W., 1991. Pulsating entry-flow in a plane channel. Z. Angew. Math. Phys. 42, 139–153. Lambossy, P., 1952. Oscillations forcées d’un liquide incompressibile et visqueux dans un tube rigide et horizontal calcul de la force frottement. Helv. Phys. Acta 25, 371–386. Lessen, M., Sadler, S., Liu, T., 1968. Stability of pipe poiseuille flow. Phys. Fluids 11, 1404–1409. Loh, S., Blackburn, H., 2010. Stability of steady flow through a corrugated pipe. In: Proceedings of the 17th Australasian Fluid Mechanics Conference. Curran Associates, Inc., pp. 329–332. Nebauer, J., Blackburn, H., 2009. Stability of oscillatory and pulsatile pipe flow. In: Proceedings of the 7th International Conference on CFD in the Minerals and Process Industries. CSIRO, Melbourne, Australia, pp. 051ZHU-1–6. Nebauer, J., Blackburn, H., 2010. On the stability of timeperiodic pipe flow. In: Proceedings of the 17th Australasian Fluid Mechanics Conference. Curran Associates, Inc., pp. 325–328. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere, Washington, DC. Patankar, S.V., Spalding, D.B., 1972. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. IJHMT 15 (10), 1787–1806. Peixinho, J., Mullin, T., 2006. Decay of turbulence in pipe flow. Phys. Rev. Lett. 96, 094501. Pfenninger, W., 1961. Boundary layer suction experiments with laminar flow at high reynolds numbers in the inlet length of a tube by various suction methods. In: Lachmann, G. (Ed.), Boundary Layer and Flow Control, vol. 2. Pergamon, pp. 961–980. Ray, S., Ünsal, B., Durst, F., Ertunç, O., Saleh, O., 2005. Mass flow rate controlled fully developed laminar pulsating pipe flows. J. Fluids Eng. 127 (3), 405–418. Reynolds, O., 1883. An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels. Phil. Trans. Royal Soc. Lond. 174, 935– 982. Rhie, C.M., Chow, W.L., 1983. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA J. 21, 1525–1532. Salwen, H., Cotton, F., Grosch, C., 1980. Linear stability of poiseuille flow in a circular pipe. J. Fluid Mech. 98, 273–284. Salwen, H., Grosch, C., 1972. The stability of poiseuille flow in a pipe of circular cross-section. J. Fluid Mech. 54, 113–127. Schmid, P., Henningson, D., 1994. Optimal energy density growth in hagenpoiseuille flow. J. Fluid Mech. 227, 197–225. Sexl, T., 1930. Über den von e.g. richardson entdeckten annulareffekt. Zeitschrift für Physik 61, 349–362. Smith, D., Blackburn, H., 2010. Transient growth analysis for axisymmetric pulsatile pipe flows in a rigid straight circular pipe. In: Proceedings of the 17th Australasian Fluid Mechanics Conference. Curran Associates, Inc., pp. 321–324.

176

S. Ray et al. / International Journal of Heat and Fluid Flow 37 (2012) 167–176

Stettler, J., Hussain, A., 1986. On transition of the pulsatile pipe flow. J. Fluid Mech. 170, 169–197. Stone, H.L., 1968. Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J. Num. Anal. 5, 530–541. Straatman, A., Khayat, R., Haj-Qasem, E., Steinman, D., 2002. On the hydrodynamic stability of pulsatile flow in a plane channel. Phys. Fluids 14, 1938–1944. Szymanski, P., 1932. Quelques solutions exactes des équations de l’hydrodynamique de fluide visqueux dans le cas d’un tube cylindrique. J. Math. Pure Appl. Ser. 9 (11), 67–107. Tozzi, J., von Kerczek, C., 1986. The stability of oscillatory hagenpoiseuille flow. ASME J. Appl. Mech. 53, 187–192. Uchida, S., 1956. The pulsating viscous flow superposed on the steady laminar motion of incompressible fluid in a circular pipe. J. Appl. Math. Phys. (ZAMP) 7, 403–422.

Ünsal, B., Ray, S., Durst, F., 2005. Pulsating laminar pipe flows with sinusoidal mass flux variations. Fluid Dyn. Res. 37, 317–333. Ünsal, B., Trimis, D., Durst, F., 2006. On-line instantaneous mass flow rate measurements through fuel injection nozzles. Int. J. Engine Res. 7, 371–380. Womersley, J.R., 1955. Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J. Physiol. 127, 553–563. Womersley, J.R., 1957. An elastic tube theory of pulse transmission and oscillatory flow in mammalian arteries. Technical Report TR56-614, WADG Technical Report. Yang, W., Yih, C.-S., 1977. Stability of time-periodic flows in a circular pipe. J. Fluid Mech. 82, 497–505.