- Email: [email protected]
Computational dynamics of unsteady flow of a variable viscosity reactive fluid in a porous pipe
Mechanics Research Communications 37 (2010) 347–353
Contents lists available at ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
Computational dynamics of unsteady flow of a variable viscosity reactive fluid in a porous pipe T. Chinyoka a,*, O.D. Makinde b a b
Center for Research in Computational and Applied Mechanics, University of Cape Town, Rondebosch 7701, South Africa Faculty of Engineering, Cape Peninsula University of Technology, P.O. Box 1906, Bellville 7535, South Africa
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 12 July 2009 Received in revised form 15 February 2010 Available online 24 February 2010 Keywords: Unsteady flow Porous pipe Variable viscosity Arrhenius reaction Finite difference method
The coupled nonlinear equations in cylindrical Cartesian coordinates governing pressure-driven unsteady flow of a reactive variable viscosity fluid and heat transfer in a circular pipe whose walls are porous, are derived and solved numerically using a semi-implicit finite difference scheme under axisymmetric conditions. The boundary conditions along the centerline of the pipe are rebuilt via an assumption on the continuity of derivatives at each stage of the computation and results are validated against the results obtained using well documented boundary conditions for flow with no suction/injection. The chemical kinetics is assumed to follow Arrhenius rate law while the fluid viscosity is an exponentially decreasing function of temperature. Both numerical and graphical results are presented and discussed quantitatively with respect to various parameters embedded in the problem. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The motion of reactive variable viscosity fluid through porous media at low Reynolds numbers has long been an important subject in the field of chemical, biomedical, and environmental engineering and science. This phenomenon is fundamental in nature and is of great practical importance in many diverse applications, including the production of oil and gas from geological structures, the gasification of coal, the retorting of shale oil, filtration, groundwater movement, regenerative heat exchange, surface catalysis of chemical reactions, adsorption, coalescence, drying, ion exchange and chromatography. Some of the applications mentioned above involve two or even three fluids, and multidimensional and unsteady flows. In some of applications the details of the local velocity field are of concern. Moreover, the flow of fluid through porous media and in the tube of circular or rectangular cross-section with the porous wall has long been investigated in many engineering applications. The following cited works are in no way exhaustive but are a glimpse of the representative works. Most of the investigations have been performed under steady and isothermal conditions and solved via analytical techniques. For such investigations in which the flow is confined to porous channels, we note the works of Wang et al. (2001), Makinde (1995), Oxarango et al. (2004), Deng and Martinez (2005). An
* Corresponding author. E-mail addresses: [email protected] (O.D. Makinde).
(T.
Chinyoka),
[email protected]
0093-6413/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2010.02.007
examples of steady isothermal investigations of flow in a porous pipes is Makinde (1996). Non-isothermal and/or unsteady flow investigations with porous boundaries include the works of Fang (2004), Makinde (2006, 2007) and Makinde and Ogulu (2008). Due to the increased complexities, some of these works necessarily made a recourse to numerical solution techniques. All these works however were limited to Newtonian fluids. An extension to nonNewtonian fluids was investigated in Ariel (2002). The mechanics of nonlinear flows presents a special challenge to engineers, physics and mathematicians since the nonlinearity can manifest itself in a variety of ways as is the case in the analysis of a reactive variable viscosity fluid flow in a slit with wall suction/ injection. In our case, the nonlinearity arises due to the temperature dependent viscosity and Arrhenius kinetic. To date, a pressure-driven laminar flow of a reactive variable viscosity fluid confined in a circular pipe with porous walls has not been analyzed to the authors’ knowledge under unsteady conditions. The goal here is to develop a model with which to determine the transient velocity profiles of a reactive and temperature dependent viscosity fluid whose flow is primarily driven by an adverse pressure gradient and also slightly disturbed by the transverse mass injection/ suction in the upwards/downwards direction. The investigation is organized as follows: in Section 2 describes the physical and mathematical formulation of the model problem. In Section 3, we employ a semi-implicit finite difference scheme similar to that described in Chinyoka (2008) to compute the transient solutions of the problem. The results are presented graphically and discussed quantitatively in Section 4.
348
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353
2. Mathematical model Consider pressure-driven flow of a variable viscosity, reactive and incompressible fluid in a porous pipe. In the present study, it is assumed that amount of the mass injection into the pipe equals that of the mass suction at the opposite end of the pipe (see Fig. 1). In other words, velocities of the mass injection and suction at the plates are equal to one another. It is also assumed that the injection/suction velocity is not strong enough to significantly alter the predominantly unidirectional flow (parallel to the z direction) as well as the assumption of axisymmetric flow. Under these ¼u ðr; tÞ, assumptions, the longitudinal velocity component is u the transverse velocity component (in the radial direction, r ) is maintained at a constant strength V and the angular velocity is neglected. In short, the velocity vector corresponding to the cylindri ; 0). Since the suction and cal coordinate system ðr ; h; zÞ is (V; u injection are maintained at equal strengths around the pipe, the fluid maximum velocity is expected to migrate away from the centerline and into the section of the pipe closest to the suction. The axisymmetry assumption will thus in general not hold around the pipe for large suction/injection velocities. We invoke the assumption mainly to simplify the mathematical modelling. We assume an Arrhenius rate law for the chemical kinetics and hence the rate of chemical reaction increases exponentially with the absolute temperature. Generally, it well known in the literature that the viscosity of a liquid decreases with increasing temperature. The model employed in our case for variation of dynamic viscosity as a function of temperature is in agreement with this notion. Since the chemical kinetics is modeled by an Arrhenius law, we choose a viscosity functions that decreases exponentially (as opposed to linearly, quadratically, etc.) with temperature. We introduce the following dimensionless variables and quantities,
l0t l E ; l¼ exp ; RT w l0 qR20 EðT T w Þ R2 p RT w l cp ; p¼ 0 ; b¼ ; Pr ¼ 0 ; T¼ 2 L l E j RT w 0V qVR0 Rew ¼ ; l0 QC 0 AER20 E l EV 2 E exp k¼ ; Br ¼ 0 2 exp ; 2 RT w RT w jRT w jRT w r¼
r ; R0
z z¼ ; L
u¼
u ; V
t¼
ð1Þ
where l0, q, j are the dynamic viscosity, fluid density and thermal conductivity respectively at channel wall temperature Tw, R0 is the pipe radius, L is the pipe characteristic length, E is the activation energy, R is the gas constant, A is the pre-exponential factor, Q is the heat of reaction, C0 is the reacting species concentration, cp is the specific heat at constant pressure and V is the fluid uniform suc-
tion/injection velocity (V > 0 implies fluid injection into bottom section and suction out of top section of pipe, the reverse is the case when V < 0). The dimensionless equations that describe the physical situation are similar to those presented in Makinde (2006, 2007), Schlicting and Gersten (2000), Hofffmann and Chiang (1995), Lamb (1945) and in particular are given by:
@u @u @p 1 @ @u þ Rew ¼ þ rlðTÞ ; @t @r @z r @r @r 2 @T @T 1 @ @T @u ; Pr þ Rew Pr ¼ r þ kHðTÞ þ lðTÞBr @t @r r @r @r @r
ð2Þ ð3Þ
where u is the longitudinal velocity, T is the fluid temperature, p is the pressure, t is the time, Br is the Brinkman number, Rew is the Reynolds number, @p/@z = G is the constant axial pressure gradient, k is the Frank-Kamenetski (or reaction rate) parameter, Pr is the Prandtl number and H(T) represents the chemical kinetics. The initial and boundary conditions for the flow are,
uðr; tÞ ¼ Tðr; tÞ ¼ 0;
for t 6 0;
uð1; tÞ ¼ Tð1; tÞ ¼ 0;
for t > 0:
ð4Þ ð5Þ
Boundary condition (5) is the no slip condition for velocity. Due to the suction/injection, a slip velocity is possible at the wall but this will however be neglected for the low suction/injection velocities of interest to this study. For an impermeable pipe, the boundary conditions at the pipe centerline are
@u @T ð0; tÞ ¼ ð0; tÞ ¼ 0; @r @r
for t > 0:
ð6Þ
For a porous pipe, however, cognisance is taken of the fact that the velocity and temperature fields would shift (however slightly) in such a way that their maximum values are (slightly) displaced away from the centerline. In light of this, we employ a slightly modified version of the centerline boundary conditions:
@u þ @u ð0 ; tÞ ¼ ð0 ; tÞ; @r @r @T þ @T ð0 ; tÞ ¼ ð0 ; tÞ; @r @r
for t > 0; ð7Þ for t > 0;
where the (+) and () represent forward and backward derivatives respectively. In particular, we assume that the variation of flow quantities within the single step size vicinity of the pipe centerline has constant rate. On a relatively fine mesh, this should be a reasonable assumption. The chemical kinetics H(T) is modelled via an Arrhenius rate law i.e. exponentially increases with temperature (Makinde, 2006; Schlicting and Gersten, 2000; Batchelor, 1967). As noted earlier, the model used for liquid viscosity l(T) decreases (exponentially) with an increase in temperature. This is physically expected, since as the
Fig. 1. Schematic diagram of the problem showing bottom injection and top suction.
349
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353
chemical reaction rate increases, more heat is generated in the flow system which invariably causes a corresponding decrease in the liquid viscosity. For this reason, the constants in both reaction kinetic and viscosity model are taken to be the same (Makinde, 2006; Schlicting and Gersten, 2000; Batchelor, 1967),
lðTÞ ¼ exp
T ; 1 þ bT
HðTÞ ¼ exp
T : 1 þ bT
ð8Þ
Here b is the activation energy parameter. In the following section, Eqs. (2)–(8) are solved numerically using a semi-implicit finite difference scheme similar to that in Chinyoka (2008). 3. Numerical solution Our numerical algorithm is based on the semi-implicit finite difference scheme given in Chinyoka (2008) for the non-isothermal viscoelastic case. The discretization of the governing equations is based on a linear Cartesian mesh and uniform grid on which finite-differences are taken. We approximate both the second and first spatial derivatives with second-order central differences. The equations corresponding to the first grid point (corresponding to the pipe centerline) are modified to incorporate the boundary conditions. The semi-implicit scheme for the velocity component reads:
ðnÞ ðnþaÞ uðnþ1Þ uðnÞ @ l @u @u l @u @2u þ lðnÞ 2 ; Rew þ ¼Gþ @r r @r @r Dt @r @r
ð9Þ
where 0 6 a 6 1. Here, as in the Crank–Nicolson scheme, terms given at (n + a) are taken as the averages of the corresponding terms at (n + 1) and n, i.e. a#(n+1) + (1 a)#(n). The scheme in Chinyoka (2008) corresponds to a = 1/2 as in the Crank–Nicolson method. The equation for u(n+1) then becomes: ðnþ1Þ
ðnþ1Þ
aDlðnÞ uj1 þ ð1 þ 2alðnÞ DÞuj
ðnþ1Þ
aDlðnÞ ujþ1
¼ e:t:;
ð10Þ
where e.t. = explicit terms and D = Dt/(Dr)2. The solution procedure for u(n+1) thus reduces to inversion of tri-diagonal matrices which is an advantage over a full implicit scheme. The semi-implicit integration scheme for the temperature equation is similar to that for the velocity component. Unmixed second partial derivatives of the temperature are treated implicitly:
T ðnþ1Þ T ðnÞ @ 2 T Pr ¼ 2 @r Dt
ðnþaÞ
1 @T ðnÞ @T ðnÞ þ Rew Pr þ kHðTÞðnÞ r @r @r !2 @uðnÞ ðnÞ þ l Br : @r
ð11Þ
The equation for T(n+1) thus becomes: ðnþ1Þ DT j1
a
DÞT ðnþ1Þ j
þ ðPr þ 2a
a
ðnþ1Þ Djþ1
¼ explicit terms;
4. Results and discussion Unless otherwise stated we employ the parameter values: k = 0.1, b = 0.1, Rew = 0, Pr = 1, Br = 1, a = 0.5, Dt = 0.001, t = 30 and Dy = 0.02. The alternative choice, say a = 1, gives similar results. We also remark that at around t = 8, all solutions (under the parameter values given herein) were noted to have converged to steady state! 4.1. Code validation, V = 0 If the channel walls are impermeable and hence there is no suction/injection, the velocity and temperature profiles are expectedly parabolic with maximum values along the centerline, see Fig. 2. We also utilize this impermeable wall case to validate our choice of boundary conditions. In this case, since as mentioned above, the maximum flow quantities (velocity and temperature) occur along the pipe centerline, the resultant problems admits the boundary conditions (6);
@u ðr; tÞ ¼ 0 and @r
@T ðr; tÞ ¼ 0; @r
for r ¼ 0 and all t:
Solving Eqs. (2) and (3) with boundary conditions (6) on an even finer mesh (Dy = 0.01 with Dt = 0.0001) leads to the same profiles as those shown in Fig. 2 at convergence (t = 30). In fact, under these conditions, the solutions converge to a maximum temperature of 0.0419 and maximum velocity of 0.2556 which are practically indistinguishable from the maximum temperature of 0.0419 and maximum velocity of 0.2557 obtained under our choice of boundary conditions (7) and either this fine mesh or the rougher one used in Fig. 2. 4.2. Uniform injection through bottom, V > 0 When there is fluid injection in the bottom to top direction, the maximum values of the flow quantities notably increase from their impermeable wall values as shown in Fig. 3. The reason for this is the fact that the fluid injection causes an increase in the heat transfer across the upper wall as will be shown later. 4.3. Uniform injection through top, V0 < 0 An opposite scenario to that described in Fig. 3 is the case when fluid injection is now from top wall to bottom. As shown in Fig. 4, the flow profiles are now shifted downwards (into the lower half of the pipe, all our graphs only show the upper half) leading to expected decreases in the maximum values of flow quantities in the upper half of the pipe. 4.4. Effects of the Reynolds and Prandtl numbers
ð12Þ
and the solution procedure reduces to inversion of a constant tridiagonal matrix. The schemes (10) and (12) were checked for consistency and it was shown that (10) and (12) are second-order accurate in both space and time. Thus the algorithm applied to the full coupled problem is also second-order accurate in both space and time. The algorithm was also tested for both spatial and temporal convergence and shown to be independent of both mesh size and time step size. In particular, if say the values Dr = 0.025, 0.02, 0.01 are employed, our algorithm converges to the same solutions both qualitatively and quantitatively. Similarly using say 20,000 time steps and Dt = 0.001 leads to the same quantitative/qualitative results as when using 200,000 time steps and Dt = 0.0001.
Figs. 5 and 6 illustrate in particular the decrease in fluid temperature with increasing Prandtl number, when all other parameter values are fixed. This has an obvious explanations since, an increase in Prandtl number can be thought of as a corresponding decrease in fluid thermal conductivity. For the current case of uniform injection through the lower plate, increasing suction Reynolds number causes a decrease in the fluid temperature due to an increase in the rate of heat transfer across the upper plate. 4.5. Effects of b in convergence, injection through bottom We also naturally expect that any increase in the strength of the activation energy and variable viscosity parameter (b) leads to a corresponding decrease in the fluid temperature and velocity.
350
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353 1
1
t t t t t
0.9
0.02 0.08 0.14 0.2 30
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
0.05
0.1
0.15
0.2
0.25
0.3
t t t t t
0.9
r
r
0.8
= = = = =
0 0
0.35
0.005
0.01
0.015
0.02
u
0.025
0.03
= = = = =
0.035
0.02 0.08 0.14 0.2 30
0.04
0.045
T
Fig. 2. Evolution of velocity and temperature with no suction/injection.
1
1
t t t t t
0.9
0.2 1 3 6 30
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.05
0.1
0.15
0.2
0.25
0.3
t t t t t
0.9
r
r
0.8
= = = = =
0
0.35
0
0.02
0.04
0.06
0.08
= = = = =
0.2 1 3 6 30
0.1
0.12
T
u Fig. 3. Evolution of velocity and temperature when Rew = 1 and Pr = 3.
1
1
t t t t t
0.9
0.2 1 3 6 30
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.05
0.1
0.15
0.2
t t t t t
0.9
r
r
0.8
= = = = =
0.25
u
0
0
0.005
0.01
T
0.015
0.02
= = = = =
0.2 1 3 6 30
0.025
Fig. 4. Evolution of velocity and temperature when Rew = 1.
The earlier is due to decreased strength of source terms through the increase in the quantity 1 + bT and the later is a combined effect of the decreased temperature and increased viscosity. Fig. 7 illustrates these conclusions.
4.6. Effects of k and Br in convergence, injection through bottom We naturally expect that any increase in the strength of the reaction parameter leads to a corresponding increase in the fluid
351
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353 1
1
Pr Pr Pr Pr
0.9
1 2 4 5
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.05
0.1
0.15
0.2
0.25
0.3
Pr Pr Pr Pr
0.9
r
r
0.8
= = = =
0
0.35
0
0.05
0.1
0.15
0.2
T
0.25
0.3
0.35
= = = =
1 2 4 5
0.4
0.45
u Fig. 5. Uniform injection through bottom, Rew = 1.
1
1
Re Re Re Re
0.9
0 1 1.5 2
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.02
0.04
0.06
0.08
0.1
0.12
Re Re Re Re
0.9
r
r
0.8
= = = =
0
0.14
0
0.05
0.1
0.15
0.2
T
0.25
0.3
0.35
= = = =
0.4
0 1 1.5 2
0.45
0.5
u Fig. 6. Uniform injection through bottom, Pr = 1.
1
β β β β
0.9 0.8
0 1 5 10
0.8 0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.02
0.04
0.06
T
0.08
0.1
β β β β
0.9
r
r
0.7
1
= = = =
0.12
0
0
0.05
0.1
0.15
u
0.2
0.25
= = = =
0.3
0 1 5 10
0.35
Fig. 7. Dependence on the activation energy and viscosity parameters, Pr = 3.
temperature. Since the Brinkman number is associated with the viscous dissipation term in the energy equation, higher values of the Brinkman parameter should necessarily lead to an increase in the amount of mechanical heat being generated by the shear forces in the fluid leading to increases in the fluid temperature. These observations are confirmed in Fig. 8.
4.7. Blow up of solutions As is well illustrated in the literature, see for example Chinyoka (2008), excessive increases of the reaction parameter beyond certain critical values leads to finite time blow up of solutions. From the parameters in Figs. 5 and 6, the maximum temperatures say
352
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353 1
λ λ λ λ
0.9 0.8
0.1 0.2 0.3 0.4
0.8
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
= = = =
0.5 1 2 3
0.7
0.6
0
Br Br Br Br
0.9
r
r
0.7
1
= = = =
0
0.4
0
0.05
0.1
T
0.15
0.2
0.25
T Fig. 8. Temperature dependence on the k and Br when Pr = 3.
14 Re = 1, Pr = 1 Re = 1, Pr = 3 Re = 1, Pr = 5
12
T
max
10 8 6 4 2 0
0
0.5
1
1.5
λ Fig. 9. Blow up of solutions.
2.5
0.7 0.6
2
|dT/dy|wall
|dT/dy|
wall
0.5 1.5
1
0.4 0.3 0.2
0.5 0.1 0
0
0.5
λ
1
1.5
0
0
0.5
λ
1
1.5
Fig. 10. Evolution of heat transfer rates.
for fixed Reynolds/Prandtl number increase with increasing Prandtl/Reynolds number. We expect thus that if we then progressively increase the reaction parameter, blow up of solutions should occur first in the case where the product RewPr is largest. Similarly the case where the product RewPr is smallest should have the largest critical value of k, or in other words should be the last case to experience any blow up. Fig. 9 clearly illustrates this phenomena. In Fig. 9, for each value of k we run the solution to conver-
gence (in this case we used t = 10) and then record the maximum temperature. This procedure is repeated until just before the critical value of k is reached. 4.8. Wall heat transfer In this case we compute the solution at t = 20. For fluid injection through bottom we use, Rew = 1 and in case fluid injection is
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353
through the top, we use, Rew = 1. The rate of heat transfer at the upper wall qw is directly proportional to the jump in temperature at the wall:
@T qw ¼ j : @y w Fig. 10 shows that as the reaction parameter increases, the rate of heat transfer at the upper wall increases in both cases and more appreciably in the case when fluid is injected through the bottom. Explanations of these observations easily follow from all the earlier results which recorded higher fluid temperatures in the upper half of the pipe (and thus in the vicinity of the upper wall) when fluid is injected through the bottom.
5. Conclusion We investigate the dynamics of a temperature dependent viscosity fluid subjected to axisymmetric one-dimensional pressuredriven flow and Arrhenius kinetics in a circular pipe. To this end, we demonstrate (at low values of the Reynolds and Prandtl numbers) both: (i) a decrease in attainable temperatures and (ii) an increase in the critical values of the reaction parameter. Both of these observations are largely attributable to the low thermal conductivities related to low Prandtl numbers. We have also checked our finite difference codes for both temporal and spatial convergence.
353
References Ariel, P.D., 2002. On exact solutions of flow problems of a second grade fluid through two parallel porous walls. Int. J. Eng. Sci. 40, 913–941. Batchelor, G.K., 1967. An Introduction to Fluid Dynamics. Cambridge University Press. pp. 302–325. Chinyoka, T., 2008. Computational dynamics of a thermally decomposable viscoelastic lubricant under shear. Trans. ASME, J. Fluids Eng. 130 (12), 121201 (7pages). Deng, C., Martinez, D.M., 2005. Viscous flow in a channel partially filled with a porous medium and with wall suction. Chem. Eng. Sci. 60, 329–336. Fang, T., 2004. A note on the incompressible Couette flow with porous walls. Int. Commun. Heat Mass Transfer 31 (1), 31–41. Hofffmann, K.A., Chiang, S.T., 1995. Computational fluid dynamics for engineering. Engineering Education System, 1. Wichita, KS, USA. Lamb, H., 1945. Hydrodynamics. Dover. pp. 684–696. Makinde, O.D., 1995. Laminar flow in a channel of varying width with permeable boundaries. Romanian J. Phys. 40 (4–5), 403–417. Makinde, O.D., 1996. Computer extension and bifurcation study by analytic continuation of porous tube flow. J. Math. Phys. Sci. 30, 1–24. Makinde, O.D., 2006. Thermal ignition in a reactive viscous flow through a channel filled with a porous medium. Trans. ASME, J. Heat Transfer 128, 601–604. Makinde, O.D., 2007. Thermal stability of a reactive viscous flow through a poroussaturated pipe. Int. J. Numer. Methods Heat Fluid Flow 17 (8), 836–844. Makinde, O.D., Ogulu, A., 2008. The effect of thermal radiation on the heat and mass transfer flow of a variable viscosity fluid past a vertical porous plate permeated by a transverse magnetic field. Chem. Eng. Commun. 195 (12), 1575–1584. Oxarango, L., Schmitz, F., Quintard, M., 2004. Laminar flow in channels with wall suction or injection: a new model to study multi-channel filtration systems. Chem. Eng. Sci. 59, 1039–1051. Schlicting, H., Gersten, K., 2000. Boundary Layer Theory, 8th Revised and Enlarged Edition (English). Springer-Verlag, New York. Wang, J., Gao, Z., Gan, G., Wu, D., 2001. Analytical solution of flow coefficients for a uniformly distributed porous channel. Chem. Eng. Sci. 84, 1–6.
Contents lists available at ScienceDirect
Mechanics Research Communications journal homepage: www.elsevier.com/locate/mechrescom
Computational dynamics of unsteady flow of a variable viscosity reactive fluid in a porous pipe T. Chinyoka a,*, O.D. Makinde b a b
Center for Research in Computational and Applied Mechanics, University of Cape Town, Rondebosch 7701, South Africa Faculty of Engineering, Cape Peninsula University of Technology, P.O. Box 1906, Bellville 7535, South Africa
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 12 July 2009 Received in revised form 15 February 2010 Available online 24 February 2010 Keywords: Unsteady flow Porous pipe Variable viscosity Arrhenius reaction Finite difference method
The coupled nonlinear equations in cylindrical Cartesian coordinates governing pressure-driven unsteady flow of a reactive variable viscosity fluid and heat transfer in a circular pipe whose walls are porous, are derived and solved numerically using a semi-implicit finite difference scheme under axisymmetric conditions. The boundary conditions along the centerline of the pipe are rebuilt via an assumption on the continuity of derivatives at each stage of the computation and results are validated against the results obtained using well documented boundary conditions for flow with no suction/injection. The chemical kinetics is assumed to follow Arrhenius rate law while the fluid viscosity is an exponentially decreasing function of temperature. Both numerical and graphical results are presented and discussed quantitatively with respect to various parameters embedded in the problem. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction The motion of reactive variable viscosity fluid through porous media at low Reynolds numbers has long been an important subject in the field of chemical, biomedical, and environmental engineering and science. This phenomenon is fundamental in nature and is of great practical importance in many diverse applications, including the production of oil and gas from geological structures, the gasification of coal, the retorting of shale oil, filtration, groundwater movement, regenerative heat exchange, surface catalysis of chemical reactions, adsorption, coalescence, drying, ion exchange and chromatography. Some of the applications mentioned above involve two or even three fluids, and multidimensional and unsteady flows. In some of applications the details of the local velocity field are of concern. Moreover, the flow of fluid through porous media and in the tube of circular or rectangular cross-section with the porous wall has long been investigated in many engineering applications. The following cited works are in no way exhaustive but are a glimpse of the representative works. Most of the investigations have been performed under steady and isothermal conditions and solved via analytical techniques. For such investigations in which the flow is confined to porous channels, we note the works of Wang et al. (2001), Makinde (1995), Oxarango et al. (2004), Deng and Martinez (2005). An
* Corresponding author. E-mail addresses: [email protected] (O.D. Makinde).
(T.
Chinyoka),
[email protected]
0093-6413/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechrescom.2010.02.007
examples of steady isothermal investigations of flow in a porous pipes is Makinde (1996). Non-isothermal and/or unsteady flow investigations with porous boundaries include the works of Fang (2004), Makinde (2006, 2007) and Makinde and Ogulu (2008). Due to the increased complexities, some of these works necessarily made a recourse to numerical solution techniques. All these works however were limited to Newtonian fluids. An extension to nonNewtonian fluids was investigated in Ariel (2002). The mechanics of nonlinear flows presents a special challenge to engineers, physics and mathematicians since the nonlinearity can manifest itself in a variety of ways as is the case in the analysis of a reactive variable viscosity fluid flow in a slit with wall suction/ injection. In our case, the nonlinearity arises due to the temperature dependent viscosity and Arrhenius kinetic. To date, a pressure-driven laminar flow of a reactive variable viscosity fluid confined in a circular pipe with porous walls has not been analyzed to the authors’ knowledge under unsteady conditions. The goal here is to develop a model with which to determine the transient velocity profiles of a reactive and temperature dependent viscosity fluid whose flow is primarily driven by an adverse pressure gradient and also slightly disturbed by the transverse mass injection/ suction in the upwards/downwards direction. The investigation is organized as follows: in Section 2 describes the physical and mathematical formulation of the model problem. In Section 3, we employ a semi-implicit finite difference scheme similar to that described in Chinyoka (2008) to compute the transient solutions of the problem. The results are presented graphically and discussed quantitatively in Section 4.
348
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353
2. Mathematical model Consider pressure-driven flow of a variable viscosity, reactive and incompressible fluid in a porous pipe. In the present study, it is assumed that amount of the mass injection into the pipe equals that of the mass suction at the opposite end of the pipe (see Fig. 1). In other words, velocities of the mass injection and suction at the plates are equal to one another. It is also assumed that the injection/suction velocity is not strong enough to significantly alter the predominantly unidirectional flow (parallel to the z direction) as well as the assumption of axisymmetric flow. Under these ¼u ðr; tÞ, assumptions, the longitudinal velocity component is u the transverse velocity component (in the radial direction, r ) is maintained at a constant strength V and the angular velocity is neglected. In short, the velocity vector corresponding to the cylindri ; 0). Since the suction and cal coordinate system ðr ; h; zÞ is (V; u injection are maintained at equal strengths around the pipe, the fluid maximum velocity is expected to migrate away from the centerline and into the section of the pipe closest to the suction. The axisymmetry assumption will thus in general not hold around the pipe for large suction/injection velocities. We invoke the assumption mainly to simplify the mathematical modelling. We assume an Arrhenius rate law for the chemical kinetics and hence the rate of chemical reaction increases exponentially with the absolute temperature. Generally, it well known in the literature that the viscosity of a liquid decreases with increasing temperature. The model employed in our case for variation of dynamic viscosity as a function of temperature is in agreement with this notion. Since the chemical kinetics is modeled by an Arrhenius law, we choose a viscosity functions that decreases exponentially (as opposed to linearly, quadratically, etc.) with temperature. We introduce the following dimensionless variables and quantities,
l0t l E ; l¼ exp ; RT w l0 qR20 EðT T w Þ R2 p RT w l cp ; p¼ 0 ; b¼ ; Pr ¼ 0 ; T¼ 2 L l E j RT w 0V qVR0 Rew ¼ ; l0 QC 0 AER20 E l EV 2 E exp k¼ ; Br ¼ 0 2 exp ; 2 RT w RT w jRT w jRT w r¼
r ; R0
z z¼ ; L
u¼
u ; V
t¼
ð1Þ
where l0, q, j are the dynamic viscosity, fluid density and thermal conductivity respectively at channel wall temperature Tw, R0 is the pipe radius, L is the pipe characteristic length, E is the activation energy, R is the gas constant, A is the pre-exponential factor, Q is the heat of reaction, C0 is the reacting species concentration, cp is the specific heat at constant pressure and V is the fluid uniform suc-
tion/injection velocity (V > 0 implies fluid injection into bottom section and suction out of top section of pipe, the reverse is the case when V < 0). The dimensionless equations that describe the physical situation are similar to those presented in Makinde (2006, 2007), Schlicting and Gersten (2000), Hofffmann and Chiang (1995), Lamb (1945) and in particular are given by:
@u @u @p 1 @ @u þ Rew ¼ þ rlðTÞ ; @t @r @z r @r @r 2 @T @T 1 @ @T @u ; Pr þ Rew Pr ¼ r þ kHðTÞ þ lðTÞBr @t @r r @r @r @r
ð2Þ ð3Þ
where u is the longitudinal velocity, T is the fluid temperature, p is the pressure, t is the time, Br is the Brinkman number, Rew is the Reynolds number, @p/@z = G is the constant axial pressure gradient, k is the Frank-Kamenetski (or reaction rate) parameter, Pr is the Prandtl number and H(T) represents the chemical kinetics. The initial and boundary conditions for the flow are,
uðr; tÞ ¼ Tðr; tÞ ¼ 0;
for t 6 0;
uð1; tÞ ¼ Tð1; tÞ ¼ 0;
for t > 0:
ð4Þ ð5Þ
Boundary condition (5) is the no slip condition for velocity. Due to the suction/injection, a slip velocity is possible at the wall but this will however be neglected for the low suction/injection velocities of interest to this study. For an impermeable pipe, the boundary conditions at the pipe centerline are
@u @T ð0; tÞ ¼ ð0; tÞ ¼ 0; @r @r
for t > 0:
ð6Þ
For a porous pipe, however, cognisance is taken of the fact that the velocity and temperature fields would shift (however slightly) in such a way that their maximum values are (slightly) displaced away from the centerline. In light of this, we employ a slightly modified version of the centerline boundary conditions:
@u þ @u ð0 ; tÞ ¼ ð0 ; tÞ; @r @r @T þ @T ð0 ; tÞ ¼ ð0 ; tÞ; @r @r
for t > 0; ð7Þ for t > 0;
where the (+) and () represent forward and backward derivatives respectively. In particular, we assume that the variation of flow quantities within the single step size vicinity of the pipe centerline has constant rate. On a relatively fine mesh, this should be a reasonable assumption. The chemical kinetics H(T) is modelled via an Arrhenius rate law i.e. exponentially increases with temperature (Makinde, 2006; Schlicting and Gersten, 2000; Batchelor, 1967). As noted earlier, the model used for liquid viscosity l(T) decreases (exponentially) with an increase in temperature. This is physically expected, since as the
Fig. 1. Schematic diagram of the problem showing bottom injection and top suction.
349
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353
chemical reaction rate increases, more heat is generated in the flow system which invariably causes a corresponding decrease in the liquid viscosity. For this reason, the constants in both reaction kinetic and viscosity model are taken to be the same (Makinde, 2006; Schlicting and Gersten, 2000; Batchelor, 1967),
lðTÞ ¼ exp
T ; 1 þ bT
HðTÞ ¼ exp
T : 1 þ bT
ð8Þ
Here b is the activation energy parameter. In the following section, Eqs. (2)–(8) are solved numerically using a semi-implicit finite difference scheme similar to that in Chinyoka (2008). 3. Numerical solution Our numerical algorithm is based on the semi-implicit finite difference scheme given in Chinyoka (2008) for the non-isothermal viscoelastic case. The discretization of the governing equations is based on a linear Cartesian mesh and uniform grid on which finite-differences are taken. We approximate both the second and first spatial derivatives with second-order central differences. The equations corresponding to the first grid point (corresponding to the pipe centerline) are modified to incorporate the boundary conditions. The semi-implicit scheme for the velocity component reads:
ðnÞ ðnþaÞ uðnþ1Þ uðnÞ @ l @u @u l @u @2u þ lðnÞ 2 ; Rew þ ¼Gþ @r r @r @r Dt @r @r
ð9Þ
where 0 6 a 6 1. Here, as in the Crank–Nicolson scheme, terms given at (n + a) are taken as the averages of the corresponding terms at (n + 1) and n, i.e. a#(n+1) + (1 a)#(n). The scheme in Chinyoka (2008) corresponds to a = 1/2 as in the Crank–Nicolson method. The equation for u(n+1) then becomes: ðnþ1Þ
ðnþ1Þ
aDlðnÞ uj1 þ ð1 þ 2alðnÞ DÞuj
ðnþ1Þ
aDlðnÞ ujþ1
¼ e:t:;
ð10Þ
where e.t. = explicit terms and D = Dt/(Dr)2. The solution procedure for u(n+1) thus reduces to inversion of tri-diagonal matrices which is an advantage over a full implicit scheme. The semi-implicit integration scheme for the temperature equation is similar to that for the velocity component. Unmixed second partial derivatives of the temperature are treated implicitly:
T ðnþ1Þ T ðnÞ @ 2 T Pr ¼ 2 @r Dt
ðnþaÞ
1 @T ðnÞ @T ðnÞ þ Rew Pr þ kHðTÞðnÞ r @r @r !2 @uðnÞ ðnÞ þ l Br : @r
ð11Þ
The equation for T(n+1) thus becomes: ðnþ1Þ DT j1
a
DÞT ðnþ1Þ j
þ ðPr þ 2a
a
ðnþ1Þ Djþ1
¼ explicit terms;
4. Results and discussion Unless otherwise stated we employ the parameter values: k = 0.1, b = 0.1, Rew = 0, Pr = 1, Br = 1, a = 0.5, Dt = 0.001, t = 30 and Dy = 0.02. The alternative choice, say a = 1, gives similar results. We also remark that at around t = 8, all solutions (under the parameter values given herein) were noted to have converged to steady state! 4.1. Code validation, V = 0 If the channel walls are impermeable and hence there is no suction/injection, the velocity and temperature profiles are expectedly parabolic with maximum values along the centerline, see Fig. 2. We also utilize this impermeable wall case to validate our choice of boundary conditions. In this case, since as mentioned above, the maximum flow quantities (velocity and temperature) occur along the pipe centerline, the resultant problems admits the boundary conditions (6);
@u ðr; tÞ ¼ 0 and @r
@T ðr; tÞ ¼ 0; @r
for r ¼ 0 and all t:
Solving Eqs. (2) and (3) with boundary conditions (6) on an even finer mesh (Dy = 0.01 with Dt = 0.0001) leads to the same profiles as those shown in Fig. 2 at convergence (t = 30). In fact, under these conditions, the solutions converge to a maximum temperature of 0.0419 and maximum velocity of 0.2556 which are practically indistinguishable from the maximum temperature of 0.0419 and maximum velocity of 0.2557 obtained under our choice of boundary conditions (7) and either this fine mesh or the rougher one used in Fig. 2. 4.2. Uniform injection through bottom, V > 0 When there is fluid injection in the bottom to top direction, the maximum values of the flow quantities notably increase from their impermeable wall values as shown in Fig. 3. The reason for this is the fact that the fluid injection causes an increase in the heat transfer across the upper wall as will be shown later. 4.3. Uniform injection through top, V0 < 0 An opposite scenario to that described in Fig. 3 is the case when fluid injection is now from top wall to bottom. As shown in Fig. 4, the flow profiles are now shifted downwards (into the lower half of the pipe, all our graphs only show the upper half) leading to expected decreases in the maximum values of flow quantities in the upper half of the pipe. 4.4. Effects of the Reynolds and Prandtl numbers
ð12Þ
and the solution procedure reduces to inversion of a constant tridiagonal matrix. The schemes (10) and (12) were checked for consistency and it was shown that (10) and (12) are second-order accurate in both space and time. Thus the algorithm applied to the full coupled problem is also second-order accurate in both space and time. The algorithm was also tested for both spatial and temporal convergence and shown to be independent of both mesh size and time step size. In particular, if say the values Dr = 0.025, 0.02, 0.01 are employed, our algorithm converges to the same solutions both qualitatively and quantitatively. Similarly using say 20,000 time steps and Dt = 0.001 leads to the same quantitative/qualitative results as when using 200,000 time steps and Dt = 0.0001.
Figs. 5 and 6 illustrate in particular the decrease in fluid temperature with increasing Prandtl number, when all other parameter values are fixed. This has an obvious explanations since, an increase in Prandtl number can be thought of as a corresponding decrease in fluid thermal conductivity. For the current case of uniform injection through the lower plate, increasing suction Reynolds number causes a decrease in the fluid temperature due to an increase in the rate of heat transfer across the upper plate. 4.5. Effects of b in convergence, injection through bottom We also naturally expect that any increase in the strength of the activation energy and variable viscosity parameter (b) leads to a corresponding decrease in the fluid temperature and velocity.
350
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353 1
1
t t t t t
0.9
0.02 0.08 0.14 0.2 30
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0
0.05
0.1
0.15
0.2
0.25
0.3
t t t t t
0.9
r
r
0.8
= = = = =
0 0
0.35
0.005
0.01
0.015
0.02
u
0.025
0.03
= = = = =
0.035
0.02 0.08 0.14 0.2 30
0.04
0.045
T
Fig. 2. Evolution of velocity and temperature with no suction/injection.
1
1
t t t t t
0.9
0.2 1 3 6 30
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.05
0.1
0.15
0.2
0.25
0.3
t t t t t
0.9
r
r
0.8
= = = = =
0
0.35
0
0.02
0.04
0.06
0.08
= = = = =
0.2 1 3 6 30
0.1
0.12
T
u Fig. 3. Evolution of velocity and temperature when Rew = 1 and Pr = 3.
1
1
t t t t t
0.9
0.2 1 3 6 30
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.05
0.1
0.15
0.2
t t t t t
0.9
r
r
0.8
= = = = =
0.25
u
0
0
0.005
0.01
T
0.015
0.02
= = = = =
0.2 1 3 6 30
0.025
Fig. 4. Evolution of velocity and temperature when Rew = 1.
The earlier is due to decreased strength of source terms through the increase in the quantity 1 + bT and the later is a combined effect of the decreased temperature and increased viscosity. Fig. 7 illustrates these conclusions.
4.6. Effects of k and Br in convergence, injection through bottom We naturally expect that any increase in the strength of the reaction parameter leads to a corresponding increase in the fluid
351
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353 1
1
Pr Pr Pr Pr
0.9
1 2 4 5
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.05
0.1
0.15
0.2
0.25
0.3
Pr Pr Pr Pr
0.9
r
r
0.8
= = = =
0
0.35
0
0.05
0.1
0.15
0.2
T
0.25
0.3
0.35
= = = =
1 2 4 5
0.4
0.45
u Fig. 5. Uniform injection through bottom, Rew = 1.
1
1
Re Re Re Re
0.9
0 1 1.5 2
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.02
0.04
0.06
0.08
0.1
0.12
Re Re Re Re
0.9
r
r
0.8
= = = =
0
0.14
0
0.05
0.1
0.15
0.2
T
0.25
0.3
0.35
= = = =
0.4
0 1 1.5 2
0.45
0.5
u Fig. 6. Uniform injection through bottom, Pr = 1.
1
β β β β
0.9 0.8
0 1 5 10
0.8 0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0
0
0.02
0.04
0.06
T
0.08
0.1
β β β β
0.9
r
r
0.7
1
= = = =
0.12
0
0
0.05
0.1
0.15
u
0.2
0.25
= = = =
0.3
0 1 5 10
0.35
Fig. 7. Dependence on the activation energy and viscosity parameters, Pr = 3.
temperature. Since the Brinkman number is associated with the viscous dissipation term in the energy equation, higher values of the Brinkman parameter should necessarily lead to an increase in the amount of mechanical heat being generated by the shear forces in the fluid leading to increases in the fluid temperature. These observations are confirmed in Fig. 8.
4.7. Blow up of solutions As is well illustrated in the literature, see for example Chinyoka (2008), excessive increases of the reaction parameter beyond certain critical values leads to finite time blow up of solutions. From the parameters in Figs. 5 and 6, the maximum temperatures say
352
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353 1
λ λ λ λ
0.9 0.8
0.1 0.2 0.3 0.4
0.8
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
= = = =
0.5 1 2 3
0.7
0.6
0
Br Br Br Br
0.9
r
r
0.7
1
= = = =
0
0.4
0
0.05
0.1
T
0.15
0.2
0.25
T Fig. 8. Temperature dependence on the k and Br when Pr = 3.
14 Re = 1, Pr = 1 Re = 1, Pr = 3 Re = 1, Pr = 5
12
T
max
10 8 6 4 2 0
0
0.5
1
1.5
λ Fig. 9. Blow up of solutions.
2.5
0.7 0.6
2
|dT/dy|wall
|dT/dy|
wall
0.5 1.5
1
0.4 0.3 0.2
0.5 0.1 0
0
0.5
λ
1
1.5
0
0
0.5
λ
1
1.5
Fig. 10. Evolution of heat transfer rates.
for fixed Reynolds/Prandtl number increase with increasing Prandtl/Reynolds number. We expect thus that if we then progressively increase the reaction parameter, blow up of solutions should occur first in the case where the product RewPr is largest. Similarly the case where the product RewPr is smallest should have the largest critical value of k, or in other words should be the last case to experience any blow up. Fig. 9 clearly illustrates this phenomena. In Fig. 9, for each value of k we run the solution to conver-
gence (in this case we used t = 10) and then record the maximum temperature. This procedure is repeated until just before the critical value of k is reached. 4.8. Wall heat transfer In this case we compute the solution at t = 20. For fluid injection through bottom we use, Rew = 1 and in case fluid injection is
T. Chinyoka, O.D. Makinde / Mechanics Research Communications 37 (2010) 347–353
through the top, we use, Rew = 1. The rate of heat transfer at the upper wall qw is directly proportional to the jump in temperature at the wall:
@T qw ¼ j : @y w Fig. 10 shows that as the reaction parameter increases, the rate of heat transfer at the upper wall increases in both cases and more appreciably in the case when fluid is injected through the bottom. Explanations of these observations easily follow from all the earlier results which recorded higher fluid temperatures in the upper half of the pipe (and thus in the vicinity of the upper wall) when fluid is injected through the bottom.
5. Conclusion We investigate the dynamics of a temperature dependent viscosity fluid subjected to axisymmetric one-dimensional pressuredriven flow and Arrhenius kinetics in a circular pipe. To this end, we demonstrate (at low values of the Reynolds and Prandtl numbers) both: (i) a decrease in attainable temperatures and (ii) an increase in the critical values of the reaction parameter. Both of these observations are largely attributable to the low thermal conductivities related to low Prandtl numbers. We have also checked our finite difference codes for both temporal and spatial convergence.
353
References Ariel, P.D., 2002. On exact solutions of flow problems of a second grade fluid through two parallel porous walls. Int. J. Eng. Sci. 40, 913–941. Batchelor, G.K., 1967. An Introduction to Fluid Dynamics. Cambridge University Press. pp. 302–325. Chinyoka, T., 2008. Computational dynamics of a thermally decomposable viscoelastic lubricant under shear. Trans. ASME, J. Fluids Eng. 130 (12), 121201 (7pages). Deng, C., Martinez, D.M., 2005. Viscous flow in a channel partially filled with a porous medium and with wall suction. Chem. Eng. Sci. 60, 329–336. Fang, T., 2004. A note on the incompressible Couette flow with porous walls. Int. Commun. Heat Mass Transfer 31 (1), 31–41. Hofffmann, K.A., Chiang, S.T., 1995. Computational fluid dynamics for engineering. Engineering Education System, 1. Wichita, KS, USA. Lamb, H., 1945. Hydrodynamics. Dover. pp. 684–696. Makinde, O.D., 1995. Laminar flow in a channel of varying width with permeable boundaries. Romanian J. Phys. 40 (4–5), 403–417. Makinde, O.D., 1996. Computer extension and bifurcation study by analytic continuation of porous tube flow. J. Math. Phys. Sci. 30, 1–24. Makinde, O.D., 2006. Thermal ignition in a reactive viscous flow through a channel filled with a porous medium. Trans. ASME, J. Heat Transfer 128, 601–604. Makinde, O.D., 2007. Thermal stability of a reactive viscous flow through a poroussaturated pipe. Int. J. Numer. Methods Heat Fluid Flow 17 (8), 836–844. Makinde, O.D., Ogulu, A., 2008. The effect of thermal radiation on the heat and mass transfer flow of a variable viscosity fluid past a vertical porous plate permeated by a transverse magnetic field. Chem. Eng. Commun. 195 (12), 1575–1584. Oxarango, L., Schmitz, F., Quintard, M., 2004. Laminar flow in channels with wall suction or injection: a new model to study multi-channel filtration systems. Chem. Eng. Sci. 59, 1039–1051. Schlicting, H., Gersten, K., 2000. Boundary Layer Theory, 8th Revised and Enlarged Edition (English). Springer-Verlag, New York. Wang, J., Gao, Z., Gan, G., Wu, D., 2001. Analytical solution of flow coefficients for a uniformly distributed porous channel. Chem. Eng. Sci. 84, 1–6.