Numerical study of normal pressure distribution in entrance flow between parallel plates, I. Finite difference calculations.

Abstract. This paper deals with the computation of flow between two parallel plates, including the pressure distribution in the entrance region. There are few implicit solutions available for the pressure distribution in the normal or y-direction of such flows. The pressure distribution in the y-direction is thus computed for the first time for flow between parallel plates. The minimum critical Reynolds number for laminar-turbulent transition is known to be in the range from 1300 to 1400, and we have thus focused our finite difference computations of the pressure gradient in the y-direction at Reynolds numbers (Re) between 100 and 5000. Our results have enabled us to conclude that a large difference in pressure between the wall and the centerline exists near the inlet for a low Re and decreased as Re increased. The pressure at the wall is lower than that in the central core for Re [less than or equal to] 5000, indicating that the pressure distribution is contrary to Bernoulli's law across parallel plates, although the law does not apply to viscous flow.

Key words. computational fluid dynamics, numerical analysis.

AMS subject classifications. 15A15, 15A09, 15A23

1. Introduction. In this section we briefly describe the history of the problem considered in this paper, related literature, and content of our paper.

1.1. Background and objectives. Numerous investigations of laminar incompressible fluid flow have been carried out both experimentally and theoretically in the entrance region between parallel plates as well as for pipe flow. The flow between parallel plates is called "plane Poiseuille flow" or "channel flow." Shah and London  have presented an excellent overall review of previous research studies on such problems. Generally, thus far, three major variables have been studied : (i) the velocity distribution at any section, (ii) the entrance length (Le), and (iii) the pressure difference between any two sections. The results of previous research studies on the velocity distribution, entrance length, and pressure difference in dimensionless X coordinates are approximately the same at Re > 500, independent of the Reynolds number, i.e., these quantities are not functions of Re for Re > 500; see .

Regarding the critical Reynolds number (Rc) for the transition between laminar and turbulent flows between parallel plates, Davies and White  observed a minimum critical Rc of 1440 when using a channel of rectangular cross section with the aspect ratio of 37 through 165. The measurements of Davies and White were made rather close to the entry to the channel so that the flow did not have sufficient length to become fully developed and turbulent. Moreover, Patel and Head  carried out experiments in a rectangular parallel-sided channel 1/4 in. high and 12 in. wide in the fully developed region. The aspect ratio of 48 was considered large enough for the flow to be assumed two-dimensional. Patel and Head state that the approximate value of 1300 may be accepted as representing the lower or minimum critical Rc for channel flow. Reynolds  found in his pipe flow experiments that the transition occurs in the entrance region under natural calm conditions. For channel flow, the transition should occur in the entrance region as Davies and White observed.

Accordingly, the objective of our investigation is to find variables that vary or decrease in the entrance region as Re increases, particularly in the range of 100 [less than or equal to] Re [less than or equal to] 5000. As a result, we found that the pressure drop in the y-direction near the inlet decreases as Re increases. We expect that this pressure drop may lead to calculation of the minimum Rc (Kanda, ). Therefore, the primary objective of this investigation is to accurately study the pressure distribution in the y-direction for Re values ranging from 100 to 5000. To this end, we aim to answer the following question: Which is larger, pressure near the wall (Pw) or the pressure at the centerline (Pc)? The secondary objective of our paper is to develop a more accurate algorithm for the calculation of the pressure distribution without any assumptions made for pressure distribution, particularly at the wall, since Peyret and Taylor  state that the two most troublesome boundary conditions to prescribe and satisfy are (i) downstream flow conditions and (ii) pressure conditions at a solid surface. Hence, we present the boundary conditions and good numerical solutions to Poisson's equation for pressure distribution. In a follow up paper, we will perform these calculations via the use of Sinc methods.

1.2. Literature review. In the entrance region, the velocity of fluid particles near the wall decreases due to the effect of internal friction, and the velocity of fluid particles near the centerline increases until finally an equilibrium between pressure drop and friction resistance adjusts itself . That is, the pressure gradient in the entrance region differs from that of fully developed flow. The initial velocity profile, which is constant across the channel inlet, develops into a parabolic profile in a fully developed downstream region.

Since an analytical solution near the channel inlet is not possible, various approximate solutions, mostly involving Prandtl's boundary-layer approximation, have been developed. The single most important part of the boundary-layer assumption is that the static pressure may be taken to be constant across the boundary layer. Thus, in the boundary-layer assumption, although streamwise variables are retained, one critical assumption is made: the static pressure across each section is uniform, i.e., the normal pressure gradient [partial derivative]P/[partial derivative]y = 0. In most previous investigations, it was assumed that the pressure P is constant across the channel.

Unfortunately, we found little information on the normal pressure gradient [partial derivative]P/[partial derivative]y. Wang and Longwell  have investigated a numerical solution that does not use the boundary-layer assumptions for a Reynolds number of 300. They calculated two cases: Case I in which upstream effects are not considered and Case II in which upstream effects are considered. They reported that [partial derivative]P/[partial derivative]y is significant near the channel inlet. For the pressure gradients at x = 0.147 in Case II, -[partial derivative]P/[partial derivative]x was approximately 0.53 at y = 0.8 and 0.15 at y = 0.0, and -[partial derivative]P/[partial derivative]y was approximately 0.32 at y = 0.9 and -0.01 at y = 0.3. Furthermore, the excess pressure drops at x = 0.147 and y = 0.9, 0.5, and 0.1 were 0.1579, 0.0660, and 0.0449 in Case I and 0.1134, 0.0715, and 0.0778 in Case II, respectively. Accordingly, note that the pressure drop near the wall is higher than that near the centerline in both cases. This indicates that, near the channel inlet, the pressure near the wall is lower than that near the centerline. Morihara and Cheng  have investigated the numerical solution of viscous flow in the entrance region of parallel plates for Reynolds numbers between 0 and 2000. They reported that the streamwise pressure gradient is very large near the inlet and the normal pressure gradient [partial derivative]P/[partial derivative]y is also far from zero near the wall for a small x. For example, for Re = 20 and x/Re = 0.01, -([partial derivative]P/[partial derivative]x)Re/6 values were more than 4 at y = 0.75, approximately 1 at y = 0.5, 0.4 at y = 0.25, and 0.3 at y = 0. Morihara and Cheng also reported the existence of a delta-shaped adverse pressure gradient zone that extends into the entrance region at a small Reynolds number, and it is practically nonexistent at Re = 2000.

Accordingly, the pressure distribution in the y-direction has not been well determined to date for Reynolds numbers between 100 and 5000.

1.3. Nomenclature. We end this introduction by presenting the notation we use in the paper.

H = channel height

h = one-half of channel height

i = grid point in x-direction

IO = maximum number of grid points in x-direction

j = grid point in y-direction

JO = maximum number of grid points in y-direction

Le = dimensionless entrance length= [x'.sub.e]/(HRe)

Lp = pressure agreement length = x' (Lp)/(HRe)

P = dimensionless pressure= p/(1/2 * [rho][U.sup.2.sub.0])

Pc = pressure at centerline

Pw = pressure at wall

p = pressure

Rc = critical Reynolds number for transition

Re = Reynolds number = [U.sub.0]H/v

t = dimensionless time = ([U.sub.0]H)t'

t' = time

[U.sub.0] = average velocity in x-direction

u = dimensionless x component of velocity = u'/[U.sub.0]

u' = x component of velocity

v = dimensionless y component of velocity = v'/[U.sub.0]

v' = y component of velocity

x = dimensionless coordinate along channel = x'/H

x' = coordinate along channel

[x'.sub.e] = entrance length

x' (Lp) = actual pressure agreement length

X = dimensionless x-coordinate = x'/(HRe) = x/Re

y = dimensionless coordinate across channel = y'/H

y' = coordinate across channel

[psi] = dimensionless stream function= [psi]'/([U.sub.0]H)

[psi] = stream function

[omega] = dimensionless vorticity = (H/[U.sub.0])[omega]'

[omega]' = vorticity

2. Governing equations. Fig. 2.1 shows the entrance region between parallel plates at y = [+ or -] h in two dimensions. We have assumed that at the inlet x = 0, the fluid enters the channel with a flat axial velocity profile [U.sub.0] across parallel plates, and that there is no velocity component in the y-direction.

First, consider dimensionless variables. All lengths and velocities in the problem are normalized by the channel height H(= 2h) and the mean velocity [U.sub.0], respectively. The pressure is normalized by (1/2)[rho] [U.sup.2.sub.0], not [rho] [U.sup.2.sub.0] . The Reynolds number is based on the channel height H and the mean velocity [U.sub.0]. Note that x is used for calculation and X(=x/Re) for presentation in figures and tables.

2.1. Governing equations. The equations that govern the incompressible laminar flow are the vorticity transport equation,

[partial derivative][omega] / [partial derivative]t + [partial derivative][psi] [partial derivative][omega] / [partial derivative]y [partial derivative]x - [partial derivative][psi] [partial derivative]w / [partial derivative]x [partial derivative]y = 1 / Re [[nabla].sup.2][omega], (2.1)

[FIGURE 2.1 OMITTED]

and Poisson's equation for the stream function,

[[nabla].sup.2] [psi] = -[omega]. (2.2)

Hence, it is possible to avoid some assumptions for pressure distribution by introducing the vorticity and stream function as dependent variables.

The relationships between stream function and velocity are defined as

u = [partial derivative][psi] / [partial derivative]y, v = [partial derivative][psi] / [partial derivative]x (2.3)

In a two-dimensional flow field, only the z component of vorticity, [[omega].sub.z], is effective; thus, [omega] denotes [[omega].sub.z] in this study,

[omega] = [[omega].sub.z] = [[[nabla] x V].sub.z] = [partial derivative]v / [partial derivative]x - [partial derivative]u / [partial derivative]y (2.4)

The [psi] - [omega] solution does not give any information regarding the pressure field. The pressure can be calculated using the Navier-Stokes equations in a steady state : the pressure distribution for the x derivative is

[partial derivative]P / [partial derivative]x = -2 (u [partial derivative]u / [partial derivative]x + v [partial derivative]u / [partial derivative]y) + 2 / Re [[nabla].sup.2] u, (2.5)

and that for the y derivative is

[partial derivative]P / [partial derivative]y = -2 (u [partial derivative]v / [partial derivative]x + v [partial derivative]v / [partial derivative]y) + 2 / Re [[nabla].sup.2] u, (2.6)

Since u and v are known at every point from (2.3), the derivatives on the right-hand sides of (2.5) and (2.6) can be obtained. Hence, note that the result of (2.5) must satisfy the result of (2.6). Furthermore, a smooth pressure distribution which satisfies both (2.5) and (2.6) is calculated using Poisson's equation ,

[[nabla].sup.2]P = -4[([partial derivative]v / [partial derivative]x)([partial derivative]u /[partial derivative]y) - ([partial derivative]u /[partial derivative]x) ([partial derivative]v / [partial derivative]y)]. (2.7)

In the numerical method, the vorticity transport equation is first solved and then the pressure distribution equation is solved without any assumptions made for pressure distribution. For the calculation of pressure distribution, it is important to make no assumptions. In this study, initial values are obtained using (2.5) and then (2.7) is used to obtain better solutions.

2.2. Axial pressure drop at centerline. For fully developed flow where [partial derivative]P/[partial derivative]y = 0, the pressure gradient at the centerline  is given by

-dP/dx = 24/Re.

The total pressure drop from the channel inlet is expressed as the sum of the pressure drop that would occur if the flow were fully developed, plus the excess pressure drop K(X) to account for the developing region,

P(0) - P(X) = 24X + K(X). (2.8)

2.3. Normal pressure gradient at wall. Here, we consider the normal pressure gradient [partial derivative]P/[partial derivative]y. The dimensionless N-S equation in vector form  is written as

[partial derivative]V / [partial derivative]t = -grad (P / 2 + [V.sup.2] / 2) [nabla] x W. (2.9)

Since the velocity vector V = 0 at the wall, that is, the normal component of (2.9) at the wall, reduces to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.10)

the normal pressure gradient is derived from the negative normal component of the curl of vorticity at the wall. This normal pressure gradient is also presented by

[partial derivative]P / [partial derivative]n = - 2 / Re [partial derivative]W / [partial derivative]s (2.11)

where (n, s) are the normal and tangents to the wall , . Since n = -y and s = x at the wall, (2.10) and (2.11) are the same. (2.10) is, however, clearer than Eq. (2.11) when we consider a physical force mechanism in vector form. The normal component of the curl of vorticity at the wall hereafter is called normal wall strength (NWS). Since only the z component of the vorticity [[omega].sub.z] is effective in a two-dimensional flow field, NWS is expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2.12)

The following characteristics of NWS are considered.

(i) NWS is effective near the channel inlet where the vorticity gradient in the x-direction is large and decreases inversely with Re. In the fully developed region, NWS does not exist since the curl of vorticity disappears.

(ii) It is clear from (2.12) that NWS causes a pressure gradient in the y-direction, that is, the pressure gradient results from the curl of vorticity. NWS and the normal pressure gradient [partial derivative]P/[partial derivative]y have the same magnitude at the wall, but have opposite directions. When [partial derivative]P/[partial derivative]y < 0, the direction of NWS is from the wall to the centerline, as shown in Fig. 2.2. NWS causes the fluid particles near the wall to move towards the centerline in the normal direction.

(iii) When using the boundary-layer assumptions, NWS vanishes since [partial derivative]P/[partial derivative] is always neglected in the assumptions.

[FIGURE 2.2 OMITTED]

3. Numerical methods. The rectangular mesh system used is schematically shown in Fig. 2.2, where I0 and J0 are the maximum numbers of mesh points in the x- and y-directions, respectively, and I1 = I0 -1 and I2 = I0 - 2. In this paper, I0 = 1001, J0 = 51 and 101, and the dimensionless X grid space is constant: [DELTA]X = [DELTA]x/Re = 0.0001. Accordingly, since [DELTA]x = 0.0001 * Re, [DELTA]x and the maximum x-distance, x = (I0 - 1)[DELTA]x, are proportional to Re.

3.1. Vorticity transport equation. This computational scheme involves the forward-time, centered-space (FTCS) method. For unsteady problems, (2.1) in finite difference form can be solved efficiently in time using an explicit or implicit Gauss-Seidel iteration method (this study).

The implicit form for vorticity is written as

[[omega].sup.n+1] - [[omega].sup.n] / [DELTA]t + [partial derivative][[psi].sup.n] / [partial derivative]y [partial derivative] [[omega].sup.n+1] / [partial derivative]y = 1 / Re ([[partial derivative].sup.2] [[omega].sup.n+1] / [partial derivative][x.sup.2] + [[partial derivative].sup.2] [[omega].sup.n+1] / [partial derivative][y.sup.2])

Here n is the time step.

The initial condition for the stream function is given by

[psi](i,j) = (j - 1) [DELTA]y, 1[less than or equal to] i [less than or equal to] I0, 1 [less than or equal to] j [less than or equal to] J0.

Within the boundaries the initial vorticity is obtained by solving (2.2). The velocities u and v are set using (2.3) whenever the stream function is newly calculated.

The following are the boundary conditions.

(i) At the centerline: [[psi].sub.i,1] = 0, [[omega].sub.i,1] = 0, 1 [less than or equal to] i [less than or equal to] I1.

(ii) At the inlet: [[psi].sub.1,j] = (j - 1) [DELTA]y, [[omega].sub.1,j] = 0, 2 [less than or equal to] j [less than or equal to] J1.

(iii) At the wall: [[psi].sub.i,J0] = (J0 - 1) [DELTA]y, 1 [less than or equal to] i [less than or equal to] I1.

The vorticity boundary condition at no-slip walls is derived from (2.4):

[omega] = -[partial derivative]u/[partial derivative]y.

A three-point, one-sided approximation for derivatives is used to maintain second-order accuracy:

[[omega].sub.i,J0] = - 3[u.sub.i,J0] - 4[u.sub.i,J0-1] + [u.sub.i,J0-2] / 2 [DELTA]y = 4[u.sub.i,J0-1] - [u.sub.i,J0-2] / 2 [DELTA]y. (3.1)

(iv) At the outlet, the linear extrapolation method is used:

[[psi].sub.I0,j] = 2[[psi].sub.I1,j] - [[psi].sub.I2,j], [[omega].sub.I0,j] = 2[[omega].sub.I1,j] - [[omega].sub.I2,j].

3.2. Pressure distribution. Before starting iterative calculations, first, we need to set initial values for pressure distribution using (2.5). Next, we calculate a smooth pressure distribution using Poisson's equation (2.7) and the Gauss-Seidel iteration method. The following are the boundary conditions for pressure.

(i) For the pressure at the centerline, we use the three-point finite difference form; since [partial derivative]P/[partial derivative]y = 0 at y = 0,

[P.sub.i,1] = (4[P.sub.i,2] - [P.sub.i,3]) / 3, 1 [less than or equal to] i [less than or equal to] I0.

(ii) The pressure at the channel inlet is given as zero without the leading edge:

[P.sub.1,j] = 0, 1 [less than or equal to] j [less than or equal to] J1.

(iii)The pressure at the wall is derived from (2.10). For the leading edge with i = 1 and j = J0, using the three-point approximation for w, the pressure gradient is expressed as

3[P.sub.1,J0] - 4[P.sub.1,J1] + [P.sub.1,J2] / 2[DELTA]y = 2 / Re -[[omega].sub.3,J0] + 4[[omega].sub.2,J0] - 3[[omega].sub.1,J0] / 2[DELTA]x

For the wall with 2 [less than or equal to] i [less than or equal to] I1, and j = J0,

3[P.sub.i,J0] - 4[P.sub.i,J1] + [P.sub.i,J2] / 2[DELTA]y = 2 / Re -[[omega].sub.i+1,J0] - [[omega].sub.i-1,J0]/ 2[DELTA]x

(iv) For the outflow boundary conditions, the linear extrapolation method is used:

[P.sub.I0,j] = 2[P.sub.I1,j] - [P.sub.I2,j], 1 [less than or equal to] j [less than or equal to] J0.

The numerical calculations are carried out on an NEC SX-7/232H32 supercomputer that has a peak performance of 8.83 G-FLOPS/processor. At I0 = 1001, J0 = 101, and Re = 1000, the CPU time was 5306 sec.

4. Results and discussion. To evaluate the accuracy of the calculations, the calculated velocity distribution, entrance length, and excess pressure drop were compared with those obtained by the previous researchers listed in Tables 4.1 and 4.2. The calculated results for J0 = 51 and 101 are approximately the same. Since the J0 = 101 mesh system has twice refinement of the J0 = 51 mesh system, we use the calculated results at J0 = 101 when mentioning specific values hereafter. The accuracy of the calculations in this study was verified as described in the following subsections.

4.1. Velocity distribution. The calculated profiles of the axial velocity component at Re = 1000 are shown in Fig. 4.1. The circles and squares show the velocity profiles given by Schlichting (, ) and Wang and Longwell . The values on the curves are listed in Table 4.1. Schlichting solved the problem with the aid of boundary-layer theory. Wang obtained a numerical solution in a steady state. It will be noted in Fig. 4.1 and Table 4.1 that the velocity distribution is concave in the central portion for X [less than or equal to] 0.0005 at Re = 1000 as Wang found. In the case of Wang, the concave distribution can be seen for X [less than or equal to] 0.002 at Re = 300 (in the present study, X [less than or equal to] 0.0015 at Re = 300). This difference may be caused by the mesh refinement and Re; [DELTA]X = 0.0001 for our mesh system and 0.000146 for Wang's mesh system. For the velocity distribution, after X = 0.002, Schlichting's results and our results are approximately the same. After X = 0.01, Schlichting's, Wang's, and our results are approximately the same.

4.2. Entrance length (Le). The entrance length (Le) is defined as the distance from the entrance to the point where the centerline velocity reaches 98 or 99% of the fully developed value. Table 4.2 gives the Le values for 98%, 99%, and 99.9% of the fully developed value computed by other researchers, most of which are obtained regardless of Re. Chen  provides (4.1) for an entrance length for channel flow,

Le = X / Re = 0.63 / Re(0.035Re + 1) + 0.044. (4.1)

At Re [greater than or equal to] 500, Le has a constant value of 0.044, as determined using (4.1).

Our calculated Le values are also shown in Table 4.2 for 98%, 99%, and 99.9% of the fully developed value. The following are our main conclusions. The calculated results are approximately the same at Re [greater than or equal to] 300, where the Le values are approximately 0.0333 for 98%, 0.0424 for 99%, and 0.0732 for 99.9%. The Le value of 0.0424 for 99% is slightly smaller than that calculated using (4.1) and approximately equal to the value 0.0429 of Morihara and Cheng . At Re [less than or equal to] 300, Le increases as Re increases. For 99% of the fully developed value, the Le values are 0.0407 at Re = 100 and 0.0421 at Re = 300. These values are smaller than those calculated using (4.1) (see Table 4.2). According to Eq. (18), Le decreases as Re increases for Re [less than or equal to] 300. In contrast, our results show that Le increases as Re increases for Re [less than or equal to] 300. Thus, according to our results, it seems that Le increases as Re increases for Re [less than or equal to] 300 because of a viscous effect, and that Le becomes constant at Re [greater than or equal to] 300.

[FIGURE 4.1 OMITTED]

4.3. Excess pressure drop (K([infinity])). The pressure drop can be conveniently represented by (2.8). If X is larger than 0.1, K(X) is assumed to be K([infinity]) for the fully developed region. Table 4.2 also gives the K([infinity]) values computed by other researchers, most of which are obtained regardless of Re. Chen  provides the following, (4.2) for K([infinity]):

K ([infinity]) = 0.64 + 19 / Re (4.2)

In our study, the excess pressure drop K (X) at the centerline is shown in Fig. 4.2 for Re = 1000, in which the curve of the pressure drop approaches a straight line with a gradient of 24 as X increases. The calculated K([infinity]) values for 100 [less than or equal to] Re [less than or equal to] 5000 are listed in Table 4.2. The following are our main conclusions. The calculated results are approximately the same at Re [greater than or equal to] 300, where the K([infinity]) value is approximately 0.643. At Re [less than or equal to] 300, K([infinity]) decreases as Re increases. The K([infinity]) values are 0.721 at Re = 100 and 0.643 for Re = 300. These values are slightly smaller than those calculated using Eq. (19) (see Table 4.2). This difference is attributed to the fact that Chen invoked the momentum integral method and boundary-layer assumptions. However, the present behavior that K([infinity]) decreases to a constant value as Re increases is the same as the result of Chen. For example, K([infinity]) decreases to 0.659 at Re = 1000 in Chen's case and K([infinity]) decreases to 0.643 at Re = 300 in the present study.

Generally, the results of this study for Le and K([infinity]) are in agreement with the previous results listed in Table 4.2.

[FIGURE 4.2 OMITTED]

4.4. Pressure distribution in the y-direction. Let us consider the question: "Which is higher, the pressure at the wall (Pw), or the pressure at the centerline (Pc) in the y-direction?" At the wall, the derivative of v with respect to x is 0, so that, from (2.4) and (3.1) the vorticity is reduced to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.3)

Near the wall, the x component of velocity, u, can be linearly approximated as follows:

[u.sub.i,J0-1] [approximately equal to] ([u.sub.i,J0] + [u.sub.i,J0-2]) / 2 = 1 / 2 [u.sub.i,J0-2] (4.4)

From (4.3) and (4.4), the vorticity at the wall is simply approximated as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4.5)

Substituting (4.5) into (2.10) gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thus, since [u.sub.i+1,J0-1] < [u.sub.i-1,J0-1] in the entrance region, the normal pressure gradient at the wall becomes negative. On the other hand, in the fully developed region, since [u.sub.i+1,J0-1] = [u.sub.i-1,J0-1], the normal pressure gradient at the wall becomes 0, so that the pressure distribution results in uniform in the y-direction.

As for the value of the vorticity at the wall, the velocity distribution in the fully developed region is given by

u(y) = 3/2 (1- [y.sub.2] / [h.sub.2]). (4.6)

Differentiating (4.6) with respect to y gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the dimensionless value of h is 0.5. Thus, the value of [omega] decreases monotonously from a large positive value at the leading edge to 6 in the fully developed region (see Table 4.3). Hence, the positiveness of the vorticity in (4.5) is also verified. In Table 4.3, there exists a singular point at the leading edge (i = 1, j = J0): [omega](1, J0) is 151.5 and 301.5 for J0 = 51 and 101, respectively. In this calculation method, however, after i = 2, the vorticity is smoothed regardless the number of normal grid points.

Next, let us verify the above deductions using the calculated results. The axial pressure drops [DELTA]P and the pressure distribution in the y-direction, P, at X [less than or equal to] 0.002 are shown in Figs. 4.3 through 4.9. Since the initial and boundary conditions for pressure at the inlet X = 0 are assumed to always be zero, the pressure drop from the inlet [DELTA]P is

[DELTA]P(X) = 0 - P(X) = -P(X),

where the pressure P(X) is negative and the pressure drop [DELTA]P is positive. The absolute value of pressure is equal to the value of the pressure drop. In Figs. 4.3 through 4.9 (a), the squares and circles show the pressure drops at the wall [DELTA]Pw and at the centerline [DELTA]Pc, respectively.

We consider the pressure drop. For example, at Re = 1000, it is clear in Fig. 4.6 (a) that (i) there is a large difference between [DELTA]Pw and [DELTA]Pc across the width of the channel at X < 0.001. Note that [DELTA]Pw is larger than [DELTA]Pc. This indicates that (ii) the pressure at the wall Pw is lower than the pressure at the centerline Pc, i.e., Pw < Pc. This is in contradiction to the results obtained using the boundary layer theory, in which there is a flat distribution in the central region, and also to Bernoulli's law, although the law does not apply to viscous flow. Moreover, it is seen from Figs. 4.3 through 4.9 (a) that (iii) the difference, Pc - Pw, decreases as Re increases. The values of Pc - Pw are listed in Table 4.4. It is clear from Table 4.4 that at Re [less than or equal to] 3000 and X = 0.0005, there exists a significant pressure gradient in the y-direction. The distance where [DELTA]Pw approximately equals [DELTA]Pc is hereafter called the pressure agreement length Lp. In Table 4.4, it is seen that (iv) Lp decreases as Re increases. At Re [greater than or equal to] 300, Le and K([infinity]) are approximately constant. However, even if Re [greater than or equal to] 300, the pressure difference Pc - Pw and Lp decrease as Re increases, and disappear at Re [greater than or equal to] 5000.

Next, we consider the pressure distribution in the y-direction as shown in Figs. 4.3 through 4.9 (b). At Re [less than or equal to] 1000 and X [less than or equal to] 0.0005, the difference in pressure between the wall and the centerline exists and the absolute value of normal pressure gradient [partial derivative]P/[partial derivative]y near the wall is larger than that in the central core. This normal pressure gradient decreases as Re and X increase. At Re [greater than or equal to] 5000 and X [greater than or equal to] 0.0005, the pressure difference in the y-direction disappears as shown in Fig. 4.9.

Consider the relationship between the pressure distribution and Re. In Fig. 4.10, the straight and the dotted curves denote the pressure drops at the wall and at the centerline at Re = 500, 1000, and 3000. Fig. 4.10 shows that the difference Pc - Pw depends strongly on Re for Re [less than or equal to] 3000 near the inlet. For Re [greater than or equal to] 1000, the pressure drops at the wall and at the centerline becomes the same after X [greater than or equal to] 0.00015 regardless of Re. Strictly, as X increases, the pressure drop at Re = 500 becomes approximately the same as that at Re [greater than or equal to] 1000: At X = 0.02, [DELTA]P = 1.0242 at Re = 500 and [DELTA]P = 1.0225 at Re = 1000. Similarly, as X increases far downstream, the pressure drop at Re = 300 becomes approximately the same as that at Re [greater than or equal to] 1000. However, the pressure drop at Re = 100 is higher than that at Re = 1000: At X = 0.03, [DELTA]P = 1.3996 at Re = 100, [DELTA]P = 1.3200 at Re = 300, and [DELTA]P = 1.3113 at Re = 1000. Thus, K([infinity]) becomes approximately the same at Re [greater than or equal to] 300 as shown in Table 4.2.

In summary, there exists a large difference in pressure in the y-direction near the inlet, which depends on Re. The difference disappears downstream and the pressure drop at Re [greater than or equal to] 300 becomes the same regardless of Re.

[FIGURE 4.3 OMITTED]

[FIGURE 4.4 OMITTED]

5. Conclusions. An analysis of flow development at Reynolds numbers from 100 to 5000 in the entrance region of parallel plates is presented. The flow field was calculated with time using the vorticity transport equation and Poisson's equation for the stream function independently of pressure. After the computation reached the steady state, the pressure was calculated in two steps: (i) initial pressure values were obtained using the Navier-Stokes equation in a steady state (2.5), and then, (ii) smooth solutions were calculated using Poisson's equation (2.7). The Navier-Stokes equation can be expressed in vector form using (2.9). At the wall, the viscous term is expressed by the curl of vorticity so that the pressure gradient in the y-direction is given by the vorticity gradient in the x-direction (2.10). Hence, the calculation procedure for pressure distribution was carried out without any preliminary assumptions.

[FIGURE 4.5 OMITTED]

[FIGURE 4.6 OMITTED]

As a result, the pressure distribution in the y-direction was obtained for the first time for the above Reynolds numbers. For the fully developed region, the calculated velocity distribution, entrance length, and excess pressure drop were in good agreement with those reported previously by the researchers listed in Table 4.2.

The conclusions obtained can be summarized as follows.

1. There is a large difference between Pw and Pc across the width of the channel for a small X, where Pw is smaller than Pc. This is in contradiction to the results obtained using the boundary layer theory, in which there is a flat distribution in the central region, and also to Bernoulli's law, although the law does not apply to viscous flow. The difference between Pw and Pc disappears at Re [greater than or equal to] 5000. This indicates that the boundary-layer assumptions do not hold for Re [less than or equal to] 5000. Note that NWS causes the difference between Pw and Pc and forces the fluid particles to move towards the centerline.

2. The calculated Le and K([infinity]) values are approximately constant at Re [greater than or equal to] 300. Since the minimum critical Reynolds number is in the neighborhood of 1300, it is important to find variable that vary at Re [greater than or equal to] 300. Hence, we found that the pressure difference in the y-direction exists even if Re [greater than or equal to] 300 and varies inversely as Re increases, and disappear at Re > 5000. The relation between NWS and Rc may thus be a next research item.

[FIGURE 4.7 OMITTED]

[FIGURE 4.8 OMITTED]

Acknowledgments. We wish to thank Professor Y Nemoto and the staff of the Information Synergy Center, Tohoku University, Japan, for their outstanding professional services and computational environment. The valuable advice and comments of Professor Frank Stenger of the University of Utah are greatly appreciated.

REFERENCES

 J. R. BODOIA AND J. F. OSTEREE, Finite difference analysis of plane Poiseuille and Couette flow developments, Appl. Sci. Res., 10 (1961), pp. 265-276.

 R-Y. CHEN, Flow in the entrance region at low Reynolds numbers, Trans. ASME J. Fluids Engrg., 95 (1973), pp. 153-158.

 M. COLLINS AND W. R. SCHOWALTER, Laminar flow in the inlet region of a straight channel, Phys. Fluids, 5 (1962), pp. 1122-1124.

 S. J. DAVIES AND C. M. WHITE, An Experimental Study of the flow of water in pipes of rectangular section, Proc. Roy. Soc., A119 (1928), pp. 92-107.

[FIGURE 4.9 OMITTED]

[FIGURE 4.10 OMITTED]

 S. GOLDSTEIN, Modern Developments in Fluid Dynamics, Vol. I, Dover, 97, 299, 1965.

 L. S. HAN, Hydrodynamic entrance lengths for incompressible laminar flow in rectangular ducts, Trans. ASME J. Appl. Mech., 27 (1960), pp. 403-409.

 H. KANDA, Difference in critical Reynolds number between Hagen-Poiseuille and plane Poiseuille flows, Proc. of ASME Fluids Eng. Division - 2001, U. S. Rohaygi, ed., ASME FED-Vol. 256, American Society of Mechanical Engineers, New York, 2001, pp. 189-196.

 H. KANDA, Verification of a model for occurrence of turbulent flow between parallel plates, Proc. of Mechanical Engineering Congress, 2002 Japan (MECJ-02), Vol. VII, M. Tanaka and Y Matsumoto, eds., Japan Society of Mechanical Engineers, Tokyo, 2002, pp. 65-66.

 M. KIYA, S. FUKUSAKO AND M. ARIE, Effect of non-uniform inlet velocity profile on the development of a laminar flow between parallel plates, Bulletin of JSME, 15, No. 81, (1972), pp. 324-336.

 H. MORIHARA AND R. T. CHENG, Numerical solution of the viscous flow in the entrance region of parallel plates, J. Comput. Phys.,11 (1973), pp. 550-572.

 B. S. NARANG AND G. KRIS HNAMOORTHY, Laminar flow in the entrance region of parallel plates, Trans. ASME J. Appl. Mech., 43 (1976), pp. 186-188.

 T. V. NGUYEN AND I. L. MACLAINE-CROSS, Incremental pressure drop number in parallel-plate heat exchangers, Trans. ASME J. Fluids Engrg., 110 (1988), pp. 93-96.

 R. L. PANTON, Incompressible Flow, Second ed., Wiley-Interscience, New York, 1996, pp. 376-378.

 V. C. PATEL AND M. R. HEAD, Some observations on skin friction and velocity profiles in fully developed pipe and channel flows, J. Fluid Mech., 38 (1969), pp. 181-201.

 C. E. PEARSON, A computational method for viscous flow problems, J. Fluid Mech., 21 (1965), pp. 611-622.

 R. PEYRET AND T. D. TAYLOR, Computational Methods for Fluid Flow, Springer-Verlag, Berlin, 236, 1983.

 L. PRANDTL AND O. G. TIETJENS, Applied Hydro- and Aeromechanics, Dover, 22, 1957.

 O. REYNOLDS, 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, Trans. Roy. Soc. Lond., 174 (1883), pp. 935-982.

 P. J. ROACHE, Fundamentals of Computational Fluid Dynamics, Hermosa, 1998, pp. 196-200.

 R. M. SADRI AND J. M. FLORYAN,Accurate evaluation of the loss coefficient and the entrance length of the inlet region of a channel, Trans. ASME J. Fluids Engrg., 124 (2002), pp. 685-693.

 H. SCHLICHTING, Laminar channel entrance flow, ZAMM,14(1934), pp. 368-373.

 H. SCHLICHTING, Boundary-Layer Theory, Seventh ed., McGraw-Hill, New York, 1979, p. 186.

 R. K. SHAH AND A. L. LONDON, Laminar Flow Forced Convection In Ducts, Academic Press, 1978, pp. 160-169.

 E. M. SPARROW AND S. H. LIN, Flow development in the hydrodynamic entrance region of tubes and ducts, Phys. Fluids, 7 (1964), pp. 338-347.

 F. STENGER, Sinc-Pack, A package of Sine methods for every type of computation, including solution of differential and integral equations.

 Y. L. WANG AND P. A. LONGWELL, Laminar flow in the inlet section of parallel plates, AIChE Journal, 10 (3), (1964), pp. 323-329.

 F. M. WHITE, Fluid Mechanics, Fourth ed., McGraw-Hill, 1999, p. 260.

* Received December 24, 2004. Accepted for publication January 17, 2006. Recommended by F. Stenger.

KENSHU SHIMOMUKAI, ([dagger]) AND HIDESADA KANDA ([double dagger])

([dagger]) SGI Japan, Ltd., Yebisu Garden Place Tower 31F, 4-20-1 Ebisu Shibuya-ku, Tokyo 150-6031 Japan, (shimomukai@sgi.co.ip).

([double dagger]) Department of Computer Science, University of Aizu, Aizu-Wakamatsu, Fukushima 965-8580, Japan, (kanda@u-aizu.ac.jp).
TABLE 4.1
Velocity distributions

x = 0.0002

Present Work Schlichting Wang
y (Re = 1000) (Re = 300)

0.5 0.0 0.0 0.0
0.4 1.078 1.042 1.044
0.3 u 1.076 1.043 1.029
0.2 1.061 1.043 1.021
0.1 1.053 1.043 1.018
0.0 1.050 1.043 1.017
x = 0.001
0.5 0.0 0.0 0.0
0.4 0.x98 1.004 0.924
0.3 u 1.124 1.108 1.053
0.2 1.126 1.109 1.071
0.1 1.126 1.109 1.069
0.0 1.126 1.109 1.068
x=0.005
0.5 0.0 0.0 0.0
0.4 0.670 0.697 0.723
0.3 u 1.099 1.087 1.142
0.2 1.221 1.254 1.235
0.1 1.234 1.261 1.248
0.0 1.235 1.262 1.249
x=0.05
0.5 0.00 0.00 -
0.4 0.544 0.540 -
0.3 u 0.928 0.960 -
0.2 1.236 1.260 -
0.1 1.422 1.440 -
0.0 1.492 1.500 -

x = 0.0005

Present Work Schlichting Wang
y (Re = 1000) (Re = 300)

0.5 0.0 0.0 0.0
0.4 0.989 1.057 1.044
0.3 u 1.107 1.075 1.068
0.2 1.101 1.075 1.051
0.1 1.097 1.075 1.044
0.0 1.096 1.075 1.042
x = 0.002
0.5 0.0 0.0 0.0
0.4 0.x02 0.x10 0.x47
0.3 u 1.136 1.140 1.107
0.2 1.160 1.160 1.142
0.1 1.161 1.160 1.138
0.0 1.161 1.160 1.135
x=0.01
0.5 0.0 0.0 0.0
0.4 0.589 0.626 0.655
0.3 u 1.033 1.048 1.035
0.2 1.246 1.264 1.233
0.1 1.309 1.333 1.306
0.0 1.319 1.335 1.322
x=0.1
0.5 0.00 0.00 -
0.4 0.539 0.540 -
0.3 u 0.923 0.960 -
0.2 1.235 1.260 -
0.1 1.427 1.440 -
0.0 1.499 1.500 -

TABLE 4.2
Entrance lgth Le & excess prev. obtained press drop K ([infinity])

Entrance Length Le

Author Year 98% 99% 99.9%

Schlichting 1934 0.034 0.04 -
Han 1960 0.0297 0.0396 -
Bodoia & Osterle 1961 0.034 0.044 0.076
Collins & Schowalter 1962 0.034 - -
Sparrow & Lin 1964 - - -
Wang & Longwell (Case 1) 1964 0.034 - -
Kiya, Fukusako, & Arie 1972 0.0348 0.0445 -
Morihara & Cheng 1973 - 0.0429 -
Chen 1973 - 0.0454 -
- 0.0442 -
- 0.044 -
- 0.044 -
Narang & Krishnamoorthy 1976 - 0.0455 -
- 0.0422 -
Nguyen & Maclaine-cross 1988 - 0.0551 -
Sadri & Floryan 2002 - 0.0399 -

Present Work (JO = 51) 2006 0.0315 0.0406 0.0731
0.0328 0.0419 0.0745
0.0330 0.0421 0.0746
0.0331 0.0422 0.0747
0.0331 0.0422 0.0747
0.0331 0.0421 0.0747
0.0330 0.0421 0.0746
0.0330 0.0420 0.0746
Present Work (JO = 101) 2006 0.0317 0.0407 0.0716
0.0331 0.0421 0.0729
0.0333 0.0423 0.0731
0.0333 0.0424 0.0732
0.0334 0.0424 0.0732
0.0333 0.0424 0.0732
0.0333 0.0423 0.0731
0.0332 0.0423 0.0731

Entrance Length Le

Author K([infinity]) Re

Schlichting 0.601 -
Han 0.871 -
Bodoia & Osterle 0.676 -
Collins & Schowalter 0.676 -
Sparrow & Lin 0.686 -
Wang & Longwell (Case 1) 0.676 300
Kiya, Fukusako, & Arie 0.666 -
Morihara & Cheng - 2000
Chen 0.830 100
0.703 300
0.659 1000
0.650 2000
Narang & Krishnamoorthy - 1000
- 2000
Nguyen & Maclaine-cross 0.682 1000
Sadri & Floryan 0.686 1000

Present Work (JO = 51) 0.744 100
0.671 300
0.663 500
0.664 700
0.664 1000
0.664 2000
0.663 3000
0.662 5000
Present Work (JO = 101) 0.721 100
0.643 300
0.641 500
0.641 700
0.643 1000
0.643 2000
0.643 3000
0.641 5000

TABLE 4.3
Vorticity distributions at wall (j = JO)

x
Re 0.00 0.0001 0.0002

Present Work 100 151.50 105.67 72.16
(JO = 51) 300 151.50 62.24 30.96
500 151.50 47.97 23.28
1000 151.50 36.70 18.92
2000 151.50 30.79 17.82
3000 151.50 28.87 18.07
5000 151.50 27.57 18.72
Present Work 100 301.50 143.61 74.25
(JO= 101) 300 301.50 65.40 28.49
500 301.50 48.67 22.17
1000 301.50 36.79 18.54
2000 301.50 30.77 17.68
3000 301.50 28.81 17.99
5000 301.50 27.50 18.69

x
Re 0.0005 0.001 0.0015

Present Work 100 29.16 15.13 11.54
(JO = 51) 300 14.42 11.14 10.34
500 13.18 11.55 11.02
1000 13.56 12.54 11.64
2000 14.66 12.90 11.73
3000 15.04 12.93 11.71
5000 15.25 12.93 11.69
Present Work 100 26.47 14.66 11.30
(JO= 101) 300 14.00 10.96 10.23
500 13.01 11.48 10.98
1000 13.52 12.53 11.63
2000 14.65 12.90 11.72
3000 15.04 12.92 11.70
5000 15.25 12.92 11.68

x
Re 0.002 0.01 0.02

Present Work 100 9.96 7.18 6.58
(JO = 51) 300 9.95 7.46 6.65
500 10.56 7.50 6.66
1000 10.90 7.51 6.66
2000 10.92 7.51 6.66
3000 10.90 7.50 6.65
5000 10.88 7.49 6.65
Present Work 100 9.80 7.15 6.55
(JO= 101) 300 9.88 7.44 6.62
500 10.53 7.48 6.63
1000 10.89 7.49 6.64
2000 10.92 7.49 6.64
3000 10.89 7.48 6.63
5000 10.87 7.47 6.63

x
Re 0.05 0.07 0.10

Present Work 100 6.11 6.07 6.06
(JO = 51) 300 6.12 6.07 6.06
500 6.12 6.07 6.06
1000 6.12 6.07 6.06
2000 6.12 6.07 6.06
3000 6.12 6.07 6.06
5000 6.12 6.07 6.06
Present Work 100 6.08 6.04 6.03
(JO= 101) 300 6.09 6.04 6.03
500 6.09 6.04 6.03
1000 6.09 6.04 6.03
2000 6.09 6.04 6.03
3000 6.09 6.04 6.03
5000 6.09 6.04 6.03

TABLE 4.4
Pressure difference [P.sub.c]-[P.sub.w] and pressure agreement length
[L.sub.p]

X

0.0002 0.0005 0.001
Re [P.sub.c]-[P.sub.w]

Present Work 100 1.286 1.206 0.882
(JO = 51) 300 0.436 0.393 0.227
500 0.290 0.224 0.081
700 0.227 0.139 0.030
1000 0.174 0.071 0.008
2000 0.093 0.009 0.00
3000 0.058 0.002 0.00
5000 0.027 0.00 0.00
Present Work 100 1.288 1.205 0.882
(JO = 101) 300 0.437 0.395 0.227
500 0.293 0.226 0.083
700 0.226 0.138 0.036
1000 0.178 0.071 0.007
2000 0.092 0.008 0.00
3000 0.057 0.002 0.00
5000 0.025 0.00 0.00

X

0.0015 0.002

Re [P.sub.c]-[P.sub.w] [L.sub.p]

Present Work 100 0.704 0.575 0.0213
(JO = 51) 300 0.122 0.064 0.0051
500 0.027 0.011 0.0028
700 0.008 0.00 0.0020
1000 0.00 0.00 0.0015
2000 0.00 0.00 0.0007
3000 0.00 0.00 0.0006
5000 0.00 0.00 0.0005
Present Work 100 0.703 0.575 0.0212
(JO = 101) 300 0.123 0.062 0.0051
500 0.026 0.010 0.0027
700 0.008 0.00 0.0020
1000 0.00 0.00 0.0014
2000 0.00 0.00 0.0006
3000 0.00 0.00 0.0005
5000 0.00 0.00 0.0004
COPYRIGHT 2006 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.