- Email: [email protected]
Effect of inclined magnetic field on magneto fluid flow through a porous medium between two inclined wavy porous plates (numerical study)
Applied Mathematics and Computation 135 (2003) 85–103 www.elsevier.com/locate/amc
Effect of inclined magnetic field on magneto fluid flow through a porous medium between two inclined wavy porous plates (numerical study) E.F. Elshehawey, Elsayed M.E. Elbarbary, Nasser S. Elgazery * Department of Mathematics, Faculty of Education, Ain Shams University, Roxy, Heliopolis, Cairo, Egypt
Abstract In this paper, the problem of an incompressible viscous fluid moving through a porous medium (Brinkmain model) between two inclined wavy porous plates under the effects of a constant inclined magnetic field that makes an angle with the vertical axis and constant suction (or injection) is studied numerically by a method related to that of Takabatake–Ayukawa in 1982. The present approach is not restricted by any of the parameters appearing in the problem such as Reynolds number, Froude number, magnetic parameter, suction (or injection) parameter, permeability parameter, wave number and amplitude ratio. The effects of the above variable parameters on the velocity, stream function and pressure gradient profiles have been studied. Moreover, the effect of varying the Froude number and the inclined angle on the pressure gradient and pressure rise is studied. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Inclined magnetic field; Porous medium; Inclined wavy porous plates; Takabatake finite difference method
*
Corresponding author. E-mail address: [email protected] (N.S. Elgazery).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 3 1 4 - 9
86
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
1. Introduction The study of fluid flows through wavy plates has been the object of scientific research for a long time, since it is encountered in several phenomena, e.g. in the nature, the wavy shape is formed under the earth’s crust as a result of compressional stresses due to earthquakes, the generation of wind waves on water and the formation of sedimentary ripples in river channels and dunes in deserts, etc. Geophysicists are interested in the study of water flow between such wavy plates underground, not only underground water but also petroleum, which can assist in the research for water and oil wells [5]. In 1980 Sobey et al. [6,7] analyzed numerically the oscillatory flow in a sinusoidal wavy channel for a laminar flow. Harris and Street [2] presented a numerical model of a fully developed turbulent channel flow over a moving wavy boundary. Recently, Selvarajan et al. [5] studied numerically the shear flow in a channel whose walls are subjected to a wave-like force excitation using the pseudospectral collocation method. For some relevant work of interest, we refer to [1,3,4,8,9]. In this paper, the Navier–Stokes equations governing the flow of an incompressible viscous fluid through a porous medium between two inclined wavy porous plates that make an angle h with the vertical axis, under the effect of a constant inclined magnetic field that makes an angle H with the vertical axis, are solved numerically by using the finite difference technique related to the method of Takabatake–Ayukawa [8]. The effects of inclination angles H and h, magnetic parameter M, suction (or injection) parameter S, Froude number Fr and permeability parameter K on the velocity, stream function, pressure gradient and pressure rise have been studied without any restrictions on any of the parameters of the problem. Our results were discussed for several values of physical parameters of interest. We notice that, if the Reynolds numbers become high enough, the flow becomes turbulent and the present approach becomes invalid. 2. Governing equations and boundary conditions We consider the two-dimensional steady flow of an incompressible Newtonian fluid with constant viscosity through a porous medium between two infinite wavy porous plates separated at a mean distance 2h, under the effects of a constant magnetic field and constant suction (or injection) along the channel walls with velocity V0 ðUm V0 Þ=Um . Both the magnetic field and plates inclined at angles H and h to the vertical, respectively. The geometry of the plate surface is defined as gðxÞ ¼ h cos
2p x: k
ð1Þ
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
87
where is the amplitude of the wave, x is the same direction of the channel axis and k is the wavelength. The equations of motion are given by oU oV þ ¼ 0; ox oy U
ð2Þ
oU oU 1 oP m þV ¼ U þ mr2 U ox oy q ox K
U
b20 h cos HðU cos H V sin HÞ þ g sin h; q
ð3Þ
oV oV 1 oP m þV ¼ V þ mr2 V ox oy q oy K þ
b20 h sin HðU cos H V sin HÞ g cos h: q
ð4Þ
The equations of motion are subjected to the boundary conditions: oU ¼ 0 at y ¼ 0 ðsymmetry conditionÞ; oy
ð5aÞ
and U ¼ 0;
V ¼
V0 ðUm V0 Þ Um
at
y ¼ gðxÞ ðno-slip and impermeability conditionsÞ;
ð5bÞ
where U and V are the fluid velocities in the directions of x and y, respectively, V0 ðUm V0 Þ=Um the normal velocity at the two plates having positive value at the upper plate for suction and negative for blowing, q is the density of the fluid, P is the pressure, Um is the average fluid velocity, which is defined as Z 1 h Um ¼ U ðyÞ dy; h 0 V0 is a constant velocity, m is the kinematic viscosity, K is the permeability parameter, h is the electrical conductivity and b0 is the constant magnetic field. Moreover, since both the planes y ¼ 0 and y ¼ gðxÞ constitute the streamlines, the flow rate q is constant at all cross-sections of the channel. Thus the following equations can be obtained [8]: w¼0
on y ¼ 0;
w ¼ q ð¼ constantÞ on y ¼ gðxÞ: We notice that the source terms for the inclination of the plates in Eqs. (3) and (4) are hydrostatic terms and the flows calculated are entirely independent of
88
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
the inclination of the plates. Therefore, the inclination affects only the pressure field. The w–x forms of the governing equations after eliminating the pressure in the dimensionless form are given by o2 w o2 w þ ¼ x; ox2 oy 2
ð6Þ
2 ow ox ow ox 1 o x o2 x 1 ¼ x þ oy ox ox oy ðRe SÞ ox2 oy 2 K 2 2 ow ow o2 w 2 2 þM sin H cos H þ 2 cos H : sin H þ 2 ox2 oxoy oy ð7Þ The equations of motion are subjected to the boundary conditions: ow o2 w ¼ 0; ¼0 on y ¼ 0; ox oy 2 ow V0 ow ¼ ¼0 on y ¼ gðxÞ ¼ ð1 / cos 2paxÞ; w ¼ q; ; ox oy Um ð8Þ w ¼ 0;
where e /¼ ; h
h a¼ ; k
Re ¼
hUm ; m
M¼
b20 hh ; qUm
S¼
hV0 ; m
Fr ¼
ðUm V0 Þ2 ; hg
/ is the amplitude ratio, a is the wave number, Re is the Reynolds number, M is the magnetic parameter, S is the suction (or injection) parameter, Fr is the Froude number and K is the permeability parameter. w and x are the stream function and vorticity, respectively.
3. Numerical procedure The governing equations (6)–(8) are solved numerically in the finite region ABCD (see Fig. 1). From geometry of the present problem it is enough to consider the central part of the finite region, where the far up-stream and the far down-stream end effects can be ignored. As boundary conditions on the two end sections (BC), (DA) and the no-slip condition or the symmetry conditions on the boundaries AB, CD are considered, we obtain the boundary conditions as discussed [4]:
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
89
Fig. 1. Geometry of the problem, calculating region, oblique lattice coordinates and lattice points on the interior region.
ow o2 w ¼ 0; ¼0 on AB; ox oy 2 ow ¼0 on BC ðinflowÞ; w ¼ f ðyÞ; ox ow V0 ow ¼ ¼0 on CD; w ¼ q; ; ox oy Um ow w ¼ f ðyÞ; ¼0 on DA ðoutflowÞ; ox w ¼ 0;
ð9Þ
where q is the dimensionless flow rate and f ðyÞ is an arbitrary function [8]. We assume that f ðyÞ satisfies the boundary conditions to make sure that the solution is free of discontinuities. We consider oblique lattice coordinates as shown in Fig. 1; one wavelength region is divided by an integral number of meshes N in the x-direction and by an integral number of meshes m in the y-direction, thus the lattice points are numbered i and j. Then the mesh sizes in the x- and y-directions are given, respectively, by
k¼
1 ; N
hi ¼
gi m
i1 gi ¼ 1 / cos 2pa : N
ð10Þ
90
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
3.1. Finite difference approximation of the governing equations We now consider six lattice points ðx; yÞ; ðx þ k; y þ ij Þ; ðx; y þ hi Þ; ðx k; y i1;j Þ; ðx; y hi Þ; ðx k; y hi1 i1;j Þ numbered 0, 1, 2, 3, 4, and 5, respectively, as shown in Fig. 1, and try to determine the coefficients An so that the following difference equation is satisfied: 2 o w o2 w ¼ A0 w0 þ A1 w1 þ A2 w2 þ A3 w3 þ A4 w4 þ A5 w5 ; ð11Þ þ ox2 oy 2 0 where wn represents the value of the function w at the point numbered n. Expanding wn by Taylor series about the point numbered 0 and equating the coefficients of wn , we obtain a set of equations for the unknown A;sn . Substituting from Eq. (11) into Eq. (6) yields the following w difference equation with a second-order accuracy oðh2 ; k 2 Þ when centered differences are used, which is based upon Eq. (6), can be written as A0 wi;j þ A1 wiþ1;j þ A2 wi;jþ1 þ A3 wi1;j þ A4 wi;j1 þ A5 wi1;j1 þ xi;j ¼ 0: ð12Þ Substitute from Eq. (11) into Eq. (7), and approximate the derivatives ox=ox and ox=oy by applying the up-wind difference technique, which is superior in stability for nonlinear terms and is fast in convergence of the calculation. We set n ¼ w1 w3 , d ¼ w2 w4 , and decide the stream function from the signs of n and d. Then, to assure the dominance of the coefficient of x0 , the center point numbered 0 and the points in the up-stream are used in the difference approximations to ox=ox and ox=oy. Then we obtain ox x2 x0 ¼ ; oy 0 hi ox 1 i1;j ðx ¼ Þ ðx x Þ x ðn P 0; d P 0Þ; 0 3 2 0 ox 0 k hi ox x2 x0 ¼ ; oy 0 hi ox 1 i;j ðx1 x0 Þ ðx2 x0 Þ ¼ ox 0 k hi ox x0 x4 ¼ ; oy 0 hi ox 1 i1;j ðx ¼ Þ ðx x Þ x 0 3 0 4 ox 0 k hi
ðn P 0; d < 0Þ;
ðn < 0; d P 0Þ;
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
ox x0 x4 ¼ ; oy 0 hi ox 1 i;j ðx1 x0 Þ ðx2 x4 Þ ¼ ox 0 k hi
91
ðn < 0; d < 0Þ:
Consequently, the x difference equations are obtained as follows: B0 xi;j þ B1 xiþ1;j þ B2 xi;jþ1 þ B3 xi1;j þ B4 xi;j1 þ B5 xi1;j1 þ Hi;j ¼ 0: ð13Þ 3.2. Finite difference approximations on the boundary Now let us develop difference approximations to Eq. (6) on the boundaries. We consider the four points on and near the boundary AB, numbered 0, 1, 2, and 3, respectively, as shown in Fig. 2, and try to determine the coefficients bn : 2 o w o2 w o2 w ¼ b0 w0 þ b1 w1 þ b2 w2 þ b3 w3 þ b4 2 þ ox2 oy 2 0 oy 0 In the same manner as the determination of An , we can obtain easily a difference equation on the boundary AB. Then applying the boundary condition equation (9) to this difference equation, xi;j on AB can be expressed as xi;j ¼ 0
on AB ði ¼ 2ð1ÞN ; j ¼ 1Þ:
Similarly on and near the boundaries BC, CD, AD, we consider the five points shown in Fig. 2. Then the difference equations are approximated by using w at these five points and ow=ox or ow=oy at the point numbered 0. At the points A, B, C, D, we consider the four points shown in Fig. 2. Then the difference equations are approximated by using w at these four points and ow=ox and o2 w=oy 2 at the point numbered 0. By applying the boundary conditions (9) to these difference equations, xi;j on each boundary can be expressed as xi;j ¼ c0;i wi;j1 þ c1;i ðwi1;j1 SÞ þ c2;i q on CD ði ¼ 2ð1ÞN ; j ¼ m þ 1Þ; xi;j ¼ ðd0;j þ d3;j þ d4;j Þf ðyÞ þ d1;i wiþ1;j þ d2;i wiþ1;j on DA ði ¼ 1; j ¼ 2ð1Þm þ 1Þ; xi;j ¼ ðd0;j þ d3;j þ d4;j Þf ðyÞ þ d1;j wi1;j þ d2;j wi1;j
on BC ði ¼ N þ 1; j ¼ 2ð1Þm þ 1Þ; xi;j ¼ 0
on A ði ¼ 1; j ¼ 1Þ;
92
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Fig. 2. Lattice points at the boundary.
xi;j ¼ 0
on B ði ¼ N þ 1; j ¼ 1Þ;
xi;j ¼ e2 q þ e0 wi;j1 þ e1 wiþ1;j1 þ e3
on D ði ¼ 1; j ¼ m þ 1Þ;
xi;j ¼ e2 q þ e0 wi;j1 þ e1 wi1;j1 þ e3
on C ði ¼ 1; j ¼ m þ 1Þ:
Finally, when we consider the three points on and near the boundary AB shown in Fig. 3, the derivatives of w can be described by
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
93
Fig. 3. Lattice points at the boundary.
o2 w 1 ffi ðw 2w1 þ w2 Þ: oy 2 0 h2i 0 Then, by using the boundary condition (9), we obtain 1 wi;j ¼ wi;jþ1 2
on AB ði ¼ 2ð1ÞN ; j ¼ 2Þ:
Similarly the following formulas are obtained for other boundaries: wi;j ¼
i1;j f ðyÞ þ i;j wi1;j i;j þ i1;j
on BC ði ¼ N ; j ¼ 3ð1ÞmÞ;
1 wi;j ¼ ðwi;j1 þ 3q þ 2hi Þ on CD ði ¼ 2ð1ÞN 1; j ¼ mÞ; 4 wi;j ¼
i;j f ðyÞ þ i1;j wiþ1;j i;j þ i1;j
on DA ði ¼ 2; j ¼ 3ð1Þm 1Þ:
The constants are set out in Appendix A.
4. Analysis of the pressure field Since the flow induced by the infinite train of sinusoidal waves is periodic, it is enough to evaluate the pressure and stress fields in the central part of the domain, which covers one wavelength region.
94
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Expressions for the pressure terms in the dimensionless from Eqs. (3) and (4) can be written as oP oV oU 1 ox 1 ¼ U V þ U ox oy oy ðRe SÞ oy K 1 M cos HðU cos H V sin HÞ þ sin h; ð14aÞ Fr oP oV oU 1 ox 1 ¼ U V þ V oy ox ox ðRe SÞ ox K 1 þ M sin HðU cos H V sin HÞ cos h: ð14bÞ Fr These equations are easily approximated by the difference equations, and the values of the pressure gradients ðop=oxÞi;j and ðop=oyÞi;j at each lattice point can be calculated by using wij ; xij ; where the pressure rise per wavelength ðDPk Þ is given by Z 1 dP dx: ð15Þ DPk ¼ 0 dx The numerical computations have been done by the symbolic computation software Mathematica. The command (FindRoot) is used to solve the nonlinear system, where k ¼ 0:05, maxi ðhi Þ ¼ 0:05, the AccuracyGoal is 10 digits and Working Precision is 16 digits.
Fig. 4. Variation of the dimensionless velocity component U with the magnetic parameter M for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, H ¼ p=4 and S ¼ 1.
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
95
Fig. 5. Variation of the dimensionless velocity component U with the inclined angle H for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and S ¼ 1.
Fig. 6. Variation of the dimensionless velocity component U with the permeability parameter K for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, M ¼ 5, H ¼ p=4 and S ¼ 1.
5. Results and discussion The flow in the present problem is controlled by seven dimensionless parameters: the wave number a, the amplitude ratio /, which are determined by the geometry of the problem; the Reynolds number Re , the permeability parameter K, the suction (or injection) parameter S and the magnetic parameter M, which are defined from the governing equations and the Froude number Fr ,
96
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Fig. 7. Variation of the axial velocity component U with the permeability parameter K for a ¼ 0:01, / ¼ 0:4, M ¼ 5, H ¼ p=4 and S ¼ 1.
defined from the pressure equations. Not only these seven dimensionless parameters control our problem but also the inclined angles h and H, determined by the geometry of the problem. In this section, the velocity, stream function, pressure gradient, and pressure rise fields of the flow are calculated for various values of the above parameters.
Fig. 8. The stream function w versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=4 and S ¼ 1.
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
97
Fig. 9. Variation of the dimensionless velocity component U with the suction parameter S for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and H ¼ p=4.
5.1. Velocity field and stream function In fact, the inclined angle h of the wavy porous plates does not affect the velocity field. It affects only the pressure field. So in this section we can consider that this problem represents the flow of an incompressible viscous fluid through a porous medium between two wavy porous plates, under the effect of a constant inclined magnetic field that makes an angle H with the vertical axis.
Fig. 10. Variation of pressure gradient distribution oP =ox with Froude number Fr for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2 and S ¼ 1.
98
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
The velocity field can be calculated straight from the difference approximation of the stream function. The effect of M on the axial velocity U versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, H ¼ p=4 and S ¼ 1, is shown in Fig. 4. With increase in the values of M, there is a decrease in the values of the velocity U. The axial velocity U is plotted versus x in Fig. 5 for various values of H at a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and S ¼ 1. As shown, the velocity increases as H increases. Fig. 6 illustrates the effect of K on the axial velocity U versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, M ¼ 5, H ¼ p=4 and S ¼ 1. The velocity increases with increasing K. Fig. 7 refers to the effects of K and Re upon the axial velocity at the crest section at a ¼ 0:01, / ¼ 0:4, M ¼ 5, H ¼ p=4 and S ¼ 1. As shown, the velocity is decreasing with increasing Re . Also it is clear that the velocity increases as K increases. We noticed that the effect of K becomes smaller as K increases and the maximum value for our study on the velocity U is at x ¼ 0:5. Fig. 8 represents the stream function w versus x at a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=4 and S ¼ 1. This indicates to the reader that the periodic boundary conditions (9) have been imposed. Fig. 9 shows the effect of S on the axial velocity U versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and H ¼ p=4. The value of the velocity U is decreasing as S increases. 5.2. Pressure gradient and pressure rise We notice from Eqs. (14a) and (14b) that the pressure gradient increases with increasing sin h. Also this contribution of sin h to the pressure gradient is
Fig. 11. Variation of pressure rise DPk on the center axis with Froude number Fr for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2 and S ¼ 1.
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
99
Fig. 12. Variation of pressure gradient distribution oP =ox with Froude number Fr for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and sin h ¼ 0:5.
found decreasing at all values of sin h as Fr becomes larger. Also we notice from Eq. (15) that the pressure rise increases with increasing sin h. Also this contribution of sin h to the pressure rise is decreasing at all values of sin h as Fr becomes larger (see Figs. 10 and 11). The effect of Fr on the pressure gradient distribution for the case that a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and sin h ¼ 0:5 is illustrated in Fig. 12. As shown, the pressure gradient distribution decreases as Fr increases. The effect of sin h on the
Fig. 13. Variation of pressure gradient distribution oP =ox with sin h for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and Fr ¼ 0:1.
100
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Fig. 14. Variation of pressure gradient distribution oP =ox with suction parameter S for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, Fr ¼ 0:1 and sin h ¼ 0:5.
pressure gradient distribution is shown in Fig. 13 for the case that a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and Fr ¼ 0:1. As sin h increases, the pressure gradient increases. Fig. 14 represents the effect of S on the pressure gradient distribution with x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, Fr ¼ 0:1 and sin h ¼ 0:5. The value of pressure gradient increases with increasing S.
6. Conclusions This paper has presented a numerical study of an incompressible viscous fluid moving through a porous media (Brinkmain model) flow between two inclined wavy porous plates under the effect of a constant inclined magnetic field that makes an angle H with the vertical axis. The stream fu–vorticity formulation was chosen in conjunction with an ‘‘oblique lattice’’ coordinate system. The numerical method that was used has a second-order accuracy oðh2 ; k 2 Þ. The present approach is not restricted by any of the parameters appearing in the problem such as Reynolds, Froude and wave numbers, suction (or injection) and magnetic parameters, and amplitude ratio, but if the Reynolds numbers become high enough, the flow becomes turbulent and the present approach becomes invalid. The magnetic field and suction (or injection) have the same effect on the axial velocity component. Also the inclined angle of magnetic field and the permeability have the same effect on the axial velocity component. The effect of K becomes smaller as K increases. Physically,
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
101
this result can be achieved when the holes of the porous medium are very large so that the resistance of the medium may be neglected. The inclination angle of the wavy plates h and the Froude number Fr have the same effect on the pressure gradient and the pressure rise and they do not play any role in the fluid flow. The pressure gradient and pressure rise increase as Fr decreases and they increase with increasing sin h. Also this contribution of sin h to the pressure gradient and pressure rise is found decreasing at all values of sin h as the Froude number Fr becomes larger.
Appendix A
A1 ¼
1 ; k2
A0 ¼
2 A1 ð2 þ Ci;j Þ; h2i
A3 ¼ A1 ð1 þ Li;j Þ;
A4 ¼
A2 ¼
1 1 þ A1 ðCi;j 2Bi;j Þ; h2i 2
1 1 þ A1 ðCi;j þ 2Bi;j Þ; h2i 2
A5 ¼ A1 Li;j ;
where Bi;j ¼
i;j ; hi
Li;j ¼
i;j þ i1;j ; hi1
Si;j ¼
i;j þ i1;j ; hi
Ci;j ¼ Sij2 þ
hi1 Si;j 2B2i;j ; hi
and Ri;j ¼
i;j þ i1;j hi
1 ; K 1 B0 ¼ a0 þ Ni jnj jdj þ sgn n sgn d dRij ; 2 a 0 ¼ A0
B1 ¼ A1 þ Ni dfGðdÞ 1g; 1 B2 ¼ A2 þ Ni jnj sgn n sgn d dRij GðnÞ; 2 B3 ¼ A3 þ Ni dGðdÞ; B4 ¼ A 4 þ N i
1 jnj þ sgn n sgn d dRij ðGðnÞ 1Þ; 2
B5 ¼ A 5 ;
102
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
2 1 Hi;j ¼ Re M A0 þ 2 wi;j þ A1 wiþ1;j þ A2 þ 2 wi;jþ1 þ A3 wi1;j h hi i 1 2 2 þ A4 2 wi;j1 þ A5 wi1;j1 sin H þ ðwi1;j1 wi1;j Þ khi1 hi hi1 ðwi;jþ1 wi;j1 Þ Ei;j sin H cos H þ ½Ei;j cos2 H ; þ 2hi Re hi1 Ni ¼ ; Ei;j ¼ ð2i1;j þ hi1 Þðwiþ1;j 2wi;j þ wi;j1 Þ: 2khi 2h2i x ðx 6¼ 0Þ; 1 ðx P 0Þ GðxÞ ¼ and sgn x ¼ jxj 0 ðx < 0Þ 0 ðx ¼ 0Þ: ! 2ði;j þi1;j Þ Si1;j k c1;i ¼ ; c2;i ¼ ðc0;i þ c1;i Þ; Si;j þ 2khi1
1 S 2 i1;j 2i;j þ hi1 ðhi1 þ 2i1;j Þ c0;i ¼ 2 c1;i hi 2k A1 ð2i1;j þ 2i;j Þ þ 2 ; d0;j d2;j d4;j d0;j
2 S2;j þ B2;j ¼ 2 þ A1 þ B1;j L2;j ; R2;j hi
S2;j ; R2;j h1 B1;j 1 A1 h1 1 ¼ A1 ; d3;j ¼ 2 B1;j L2;j 1 ; hi h2 R2;j 2 h2 R2;j 1 A1 h1 1 ¼ 2 B1;j L2;j 1 þ ; hi 2 h2 R2;j 2 1 hN hN þ1 ¼ 2 þ A1 ðBN ;j 2SN ;j Þ þ BN ;j mN ;j 2 ; hN þ1 RN ;j hN þ1
¼ 2A1 d1;j
SN ;j ; RN ;j
d2;j ¼ A1
BN ;j ; RN ;j
d1;j ¼ 2A1
A1 hN hN 1 hN þ1 1 B m 1 þ ; N ;j N ;j h2N þ1 2 h2N þ1 hN RN ;j 1 A1 hN hN 1 hN þ1 1 d4;j ¼ 2 BN ;j mN ;j 2 1 ; hN þ1 2 hN þ1 hN RN ;j 2 h2 e0 ¼ 2 2A1 B1;mþ1 B1;mþ1 ; h1 h1 2 h2 h1 e2 ¼ 2 þ 2A1 B1;mþ1 B1;mþ1 þ ; h1 h1 h2
¼ d3;j
1
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
e1 ¼ 2A1
103
h1 B1;mþ1 ; h2
and e3 ¼ 2kA1 þ e0 : References [1] E.F. Elshehawey, E.M. Elbarbary, N.S. Elgazery, Numerical study of a magneto-fluid motion through a porous medium between two wavyplates, Int. J. Comput. Fluid Dyn. 15 (2001) 177– 181. [2] J.A. Harris, R.L. Street, Numerical simulation of turbulent flow over a moving wavy boundary Norris and Reynolds extended, Phys. Fluids 6 (1994) 924–943. [3] C. Pozridikis, A study of peristaltic flows, J. Fluid Mech. 180 (1987) 515. [4] B.V. Rathish Kumar, K.B. Naidu, A numerical study of peristaltic flows, Comput. Fluids 24 (2) (1995) 161–176. [5] S. Selvarajan et al., A numerical study of flow through wavy-walled channels, Int. J. Numer. Meth. Fluids 26 (1998) 519–531. [6] K.D. Setphanhoff, I.J. Sobey, B.J. Bellhouse, On flow through furrowed channels, Part 2: Observed flow patterns, J. Fluid Mech. 96 (1980) 27–32. [7] I.J. Sobey, On flow through furrowed channels, Part 1: Calculated flow patterns, J. Fluid Mech. 96 (1980) 1–26. [8] S. Takabatake, K. Ayukawa, Numerical study of two-dimensional peristaltic flows, J. Fluid Mech. 122 (1982) 439. [9] S. Takabatake, K. Ayukawa, A. Mori, Peristaltic pumpingin circular cylindrical tubes: a numerical study of fluid transport and its efficiency, J. Fluid Mech. 193 (1988) 267.
Effect of inclined magnetic field on magneto fluid flow through a porous medium between two inclined wavy porous plates (numerical study) E.F. Elshehawey, Elsayed M.E. Elbarbary, Nasser S. Elgazery * Department of Mathematics, Faculty of Education, Ain Shams University, Roxy, Heliopolis, Cairo, Egypt
Abstract In this paper, the problem of an incompressible viscous fluid moving through a porous medium (Brinkmain model) between two inclined wavy porous plates under the effects of a constant inclined magnetic field that makes an angle with the vertical axis and constant suction (or injection) is studied numerically by a method related to that of Takabatake–Ayukawa in 1982. The present approach is not restricted by any of the parameters appearing in the problem such as Reynolds number, Froude number, magnetic parameter, suction (or injection) parameter, permeability parameter, wave number and amplitude ratio. The effects of the above variable parameters on the velocity, stream function and pressure gradient profiles have been studied. Moreover, the effect of varying the Froude number and the inclined angle on the pressure gradient and pressure rise is studied. Ó 2002 Elsevier Science Inc. All rights reserved. Keywords: Inclined magnetic field; Porous medium; Inclined wavy porous plates; Takabatake finite difference method
*
Corresponding author. E-mail address: [email protected] (N.S. Elgazery).
0096-3003/02/$ - see front matter Ó 2002 Elsevier Science Inc. All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 0 1 ) 0 0 3 1 4 - 9
86
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
1. Introduction The study of fluid flows through wavy plates has been the object of scientific research for a long time, since it is encountered in several phenomena, e.g. in the nature, the wavy shape is formed under the earth’s crust as a result of compressional stresses due to earthquakes, the generation of wind waves on water and the formation of sedimentary ripples in river channels and dunes in deserts, etc. Geophysicists are interested in the study of water flow between such wavy plates underground, not only underground water but also petroleum, which can assist in the research for water and oil wells [5]. In 1980 Sobey et al. [6,7] analyzed numerically the oscillatory flow in a sinusoidal wavy channel for a laminar flow. Harris and Street [2] presented a numerical model of a fully developed turbulent channel flow over a moving wavy boundary. Recently, Selvarajan et al. [5] studied numerically the shear flow in a channel whose walls are subjected to a wave-like force excitation using the pseudospectral collocation method. For some relevant work of interest, we refer to [1,3,4,8,9]. In this paper, the Navier–Stokes equations governing the flow of an incompressible viscous fluid through a porous medium between two inclined wavy porous plates that make an angle h with the vertical axis, under the effect of a constant inclined magnetic field that makes an angle H with the vertical axis, are solved numerically by using the finite difference technique related to the method of Takabatake–Ayukawa [8]. The effects of inclination angles H and h, magnetic parameter M, suction (or injection) parameter S, Froude number Fr and permeability parameter K on the velocity, stream function, pressure gradient and pressure rise have been studied without any restrictions on any of the parameters of the problem. Our results were discussed for several values of physical parameters of interest. We notice that, if the Reynolds numbers become high enough, the flow becomes turbulent and the present approach becomes invalid. 2. Governing equations and boundary conditions We consider the two-dimensional steady flow of an incompressible Newtonian fluid with constant viscosity through a porous medium between two infinite wavy porous plates separated at a mean distance 2h, under the effects of a constant magnetic field and constant suction (or injection) along the channel walls with velocity V0 ðUm V0 Þ=Um . Both the magnetic field and plates inclined at angles H and h to the vertical, respectively. The geometry of the plate surface is defined as gðxÞ ¼ h cos
2p x: k
ð1Þ
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
87
where is the amplitude of the wave, x is the same direction of the channel axis and k is the wavelength. The equations of motion are given by oU oV þ ¼ 0; ox oy U
ð2Þ
oU oU 1 oP m þV ¼ U þ mr2 U ox oy q ox K
U
b20 h cos HðU cos H V sin HÞ þ g sin h; q
ð3Þ
oV oV 1 oP m þV ¼ V þ mr2 V ox oy q oy K þ
b20 h sin HðU cos H V sin HÞ g cos h: q
ð4Þ
The equations of motion are subjected to the boundary conditions: oU ¼ 0 at y ¼ 0 ðsymmetry conditionÞ; oy
ð5aÞ
and U ¼ 0;
V ¼
V0 ðUm V0 Þ Um
at
y ¼ gðxÞ ðno-slip and impermeability conditionsÞ;
ð5bÞ
where U and V are the fluid velocities in the directions of x and y, respectively, V0 ðUm V0 Þ=Um the normal velocity at the two plates having positive value at the upper plate for suction and negative for blowing, q is the density of the fluid, P is the pressure, Um is the average fluid velocity, which is defined as Z 1 h Um ¼ U ðyÞ dy; h 0 V0 is a constant velocity, m is the kinematic viscosity, K is the permeability parameter, h is the electrical conductivity and b0 is the constant magnetic field. Moreover, since both the planes y ¼ 0 and y ¼ gðxÞ constitute the streamlines, the flow rate q is constant at all cross-sections of the channel. Thus the following equations can be obtained [8]: w¼0
on y ¼ 0;
w ¼ q ð¼ constantÞ on y ¼ gðxÞ: We notice that the source terms for the inclination of the plates in Eqs. (3) and (4) are hydrostatic terms and the flows calculated are entirely independent of
88
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
the inclination of the plates. Therefore, the inclination affects only the pressure field. The w–x forms of the governing equations after eliminating the pressure in the dimensionless form are given by o2 w o2 w þ ¼ x; ox2 oy 2
ð6Þ
2 ow ox ow ox 1 o x o2 x 1 ¼ x þ oy ox ox oy ðRe SÞ ox2 oy 2 K 2 2 ow ow o2 w 2 2 þM sin H cos H þ 2 cos H : sin H þ 2 ox2 oxoy oy ð7Þ The equations of motion are subjected to the boundary conditions: ow o2 w ¼ 0; ¼0 on y ¼ 0; ox oy 2 ow V0 ow ¼ ¼0 on y ¼ gðxÞ ¼ ð1 / cos 2paxÞ; w ¼ q; ; ox oy Um ð8Þ w ¼ 0;
where e /¼ ; h
h a¼ ; k
Re ¼
hUm ; m
M¼
b20 hh ; qUm
S¼
hV0 ; m
Fr ¼
ðUm V0 Þ2 ; hg
/ is the amplitude ratio, a is the wave number, Re is the Reynolds number, M is the magnetic parameter, S is the suction (or injection) parameter, Fr is the Froude number and K is the permeability parameter. w and x are the stream function and vorticity, respectively.
3. Numerical procedure The governing equations (6)–(8) are solved numerically in the finite region ABCD (see Fig. 1). From geometry of the present problem it is enough to consider the central part of the finite region, where the far up-stream and the far down-stream end effects can be ignored. As boundary conditions on the two end sections (BC), (DA) and the no-slip condition or the symmetry conditions on the boundaries AB, CD are considered, we obtain the boundary conditions as discussed [4]:
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
89
Fig. 1. Geometry of the problem, calculating region, oblique lattice coordinates and lattice points on the interior region.
ow o2 w ¼ 0; ¼0 on AB; ox oy 2 ow ¼0 on BC ðinflowÞ; w ¼ f ðyÞ; ox ow V0 ow ¼ ¼0 on CD; w ¼ q; ; ox oy Um ow w ¼ f ðyÞ; ¼0 on DA ðoutflowÞ; ox w ¼ 0;
ð9Þ
where q is the dimensionless flow rate and f ðyÞ is an arbitrary function [8]. We assume that f ðyÞ satisfies the boundary conditions to make sure that the solution is free of discontinuities. We consider oblique lattice coordinates as shown in Fig. 1; one wavelength region is divided by an integral number of meshes N in the x-direction and by an integral number of meshes m in the y-direction, thus the lattice points are numbered i and j. Then the mesh sizes in the x- and y-directions are given, respectively, by
k¼
1 ; N
hi ¼
gi m
i1 gi ¼ 1 / cos 2pa : N
ð10Þ
90
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
3.1. Finite difference approximation of the governing equations We now consider six lattice points ðx; yÞ; ðx þ k; y þ ij Þ; ðx; y þ hi Þ; ðx k; y i1;j Þ; ðx; y hi Þ; ðx k; y hi1 i1;j Þ numbered 0, 1, 2, 3, 4, and 5, respectively, as shown in Fig. 1, and try to determine the coefficients An so that the following difference equation is satisfied: 2 o w o2 w ¼ A0 w0 þ A1 w1 þ A2 w2 þ A3 w3 þ A4 w4 þ A5 w5 ; ð11Þ þ ox2 oy 2 0 where wn represents the value of the function w at the point numbered n. Expanding wn by Taylor series about the point numbered 0 and equating the coefficients of wn , we obtain a set of equations for the unknown A;sn . Substituting from Eq. (11) into Eq. (6) yields the following w difference equation with a second-order accuracy oðh2 ; k 2 Þ when centered differences are used, which is based upon Eq. (6), can be written as A0 wi;j þ A1 wiþ1;j þ A2 wi;jþ1 þ A3 wi1;j þ A4 wi;j1 þ A5 wi1;j1 þ xi;j ¼ 0: ð12Þ Substitute from Eq. (11) into Eq. (7), and approximate the derivatives ox=ox and ox=oy by applying the up-wind difference technique, which is superior in stability for nonlinear terms and is fast in convergence of the calculation. We set n ¼ w1 w3 , d ¼ w2 w4 , and decide the stream function from the signs of n and d. Then, to assure the dominance of the coefficient of x0 , the center point numbered 0 and the points in the up-stream are used in the difference approximations to ox=ox and ox=oy. Then we obtain ox x2 x0 ¼ ; oy 0 hi ox 1 i1;j ðx ¼ Þ ðx x Þ x ðn P 0; d P 0Þ; 0 3 2 0 ox 0 k hi ox x2 x0 ¼ ; oy 0 hi ox 1 i;j ðx1 x0 Þ ðx2 x0 Þ ¼ ox 0 k hi ox x0 x4 ¼ ; oy 0 hi ox 1 i1;j ðx ¼ Þ ðx x Þ x 0 3 0 4 ox 0 k hi
ðn P 0; d < 0Þ;
ðn < 0; d P 0Þ;
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
ox x0 x4 ¼ ; oy 0 hi ox 1 i;j ðx1 x0 Þ ðx2 x4 Þ ¼ ox 0 k hi
91
ðn < 0; d < 0Þ:
Consequently, the x difference equations are obtained as follows: B0 xi;j þ B1 xiþ1;j þ B2 xi;jþ1 þ B3 xi1;j þ B4 xi;j1 þ B5 xi1;j1 þ Hi;j ¼ 0: ð13Þ 3.2. Finite difference approximations on the boundary Now let us develop difference approximations to Eq. (6) on the boundaries. We consider the four points on and near the boundary AB, numbered 0, 1, 2, and 3, respectively, as shown in Fig. 2, and try to determine the coefficients bn : 2 o w o2 w o2 w ¼ b0 w0 þ b1 w1 þ b2 w2 þ b3 w3 þ b4 2 þ ox2 oy 2 0 oy 0 In the same manner as the determination of An , we can obtain easily a difference equation on the boundary AB. Then applying the boundary condition equation (9) to this difference equation, xi;j on AB can be expressed as xi;j ¼ 0
on AB ði ¼ 2ð1ÞN ; j ¼ 1Þ:
Similarly on and near the boundaries BC, CD, AD, we consider the five points shown in Fig. 2. Then the difference equations are approximated by using w at these five points and ow=ox or ow=oy at the point numbered 0. At the points A, B, C, D, we consider the four points shown in Fig. 2. Then the difference equations are approximated by using w at these four points and ow=ox and o2 w=oy 2 at the point numbered 0. By applying the boundary conditions (9) to these difference equations, xi;j on each boundary can be expressed as xi;j ¼ c0;i wi;j1 þ c1;i ðwi1;j1 SÞ þ c2;i q on CD ði ¼ 2ð1ÞN ; j ¼ m þ 1Þ; xi;j ¼ ðd0;j þ d3;j þ d4;j Þf ðyÞ þ d1;i wiþ1;j þ d2;i wiþ1;j on DA ði ¼ 1; j ¼ 2ð1Þm þ 1Þ; xi;j ¼ ðd0;j þ d3;j þ d4;j Þf ðyÞ þ d1;j wi1;j þ d2;j wi1;j
on BC ði ¼ N þ 1; j ¼ 2ð1Þm þ 1Þ; xi;j ¼ 0
on A ði ¼ 1; j ¼ 1Þ;
92
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Fig. 2. Lattice points at the boundary.
xi;j ¼ 0
on B ði ¼ N þ 1; j ¼ 1Þ;
xi;j ¼ e2 q þ e0 wi;j1 þ e1 wiþ1;j1 þ e3
on D ði ¼ 1; j ¼ m þ 1Þ;
xi;j ¼ e2 q þ e0 wi;j1 þ e1 wi1;j1 þ e3
on C ði ¼ 1; j ¼ m þ 1Þ:
Finally, when we consider the three points on and near the boundary AB shown in Fig. 3, the derivatives of w can be described by
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
93
Fig. 3. Lattice points at the boundary.
o2 w 1 ffi ðw 2w1 þ w2 Þ: oy 2 0 h2i 0 Then, by using the boundary condition (9), we obtain 1 wi;j ¼ wi;jþ1 2
on AB ði ¼ 2ð1ÞN ; j ¼ 2Þ:
Similarly the following formulas are obtained for other boundaries: wi;j ¼
i1;j f ðyÞ þ i;j wi1;j i;j þ i1;j
on BC ði ¼ N ; j ¼ 3ð1ÞmÞ;
1 wi;j ¼ ðwi;j1 þ 3q þ 2hi Þ on CD ði ¼ 2ð1ÞN 1; j ¼ mÞ; 4 wi;j ¼
i;j f ðyÞ þ i1;j wiþ1;j i;j þ i1;j
on DA ði ¼ 2; j ¼ 3ð1Þm 1Þ:
The constants are set out in Appendix A.
4. Analysis of the pressure field Since the flow induced by the infinite train of sinusoidal waves is periodic, it is enough to evaluate the pressure and stress fields in the central part of the domain, which covers one wavelength region.
94
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Expressions for the pressure terms in the dimensionless from Eqs. (3) and (4) can be written as oP oV oU 1 ox 1 ¼ U V þ U ox oy oy ðRe SÞ oy K 1 M cos HðU cos H V sin HÞ þ sin h; ð14aÞ Fr oP oV oU 1 ox 1 ¼ U V þ V oy ox ox ðRe SÞ ox K 1 þ M sin HðU cos H V sin HÞ cos h: ð14bÞ Fr These equations are easily approximated by the difference equations, and the values of the pressure gradients ðop=oxÞi;j and ðop=oyÞi;j at each lattice point can be calculated by using wij ; xij ; where the pressure rise per wavelength ðDPk Þ is given by Z 1 dP dx: ð15Þ DPk ¼ 0 dx The numerical computations have been done by the symbolic computation software Mathematica. The command (FindRoot) is used to solve the nonlinear system, where k ¼ 0:05, maxi ðhi Þ ¼ 0:05, the AccuracyGoal is 10 digits and Working Precision is 16 digits.
Fig. 4. Variation of the dimensionless velocity component U with the magnetic parameter M for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, H ¼ p=4 and S ¼ 1.
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
95
Fig. 5. Variation of the dimensionless velocity component U with the inclined angle H for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and S ¼ 1.
Fig. 6. Variation of the dimensionless velocity component U with the permeability parameter K for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, M ¼ 5, H ¼ p=4 and S ¼ 1.
5. Results and discussion The flow in the present problem is controlled by seven dimensionless parameters: the wave number a, the amplitude ratio /, which are determined by the geometry of the problem; the Reynolds number Re , the permeability parameter K, the suction (or injection) parameter S and the magnetic parameter M, which are defined from the governing equations and the Froude number Fr ,
96
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Fig. 7. Variation of the axial velocity component U with the permeability parameter K for a ¼ 0:01, / ¼ 0:4, M ¼ 5, H ¼ p=4 and S ¼ 1.
defined from the pressure equations. Not only these seven dimensionless parameters control our problem but also the inclined angles h and H, determined by the geometry of the problem. In this section, the velocity, stream function, pressure gradient, and pressure rise fields of the flow are calculated for various values of the above parameters.
Fig. 8. The stream function w versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=4 and S ¼ 1.
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
97
Fig. 9. Variation of the dimensionless velocity component U with the suction parameter S for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and H ¼ p=4.
5.1. Velocity field and stream function In fact, the inclined angle h of the wavy porous plates does not affect the velocity field. It affects only the pressure field. So in this section we can consider that this problem represents the flow of an incompressible viscous fluid through a porous medium between two wavy porous plates, under the effect of a constant inclined magnetic field that makes an angle H with the vertical axis.
Fig. 10. Variation of pressure gradient distribution oP =ox with Froude number Fr for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2 and S ¼ 1.
98
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
The velocity field can be calculated straight from the difference approximation of the stream function. The effect of M on the axial velocity U versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, H ¼ p=4 and S ¼ 1, is shown in Fig. 4. With increase in the values of M, there is a decrease in the values of the velocity U. The axial velocity U is plotted versus x in Fig. 5 for various values of H at a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and S ¼ 1. As shown, the velocity increases as H increases. Fig. 6 illustrates the effect of K on the axial velocity U versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, M ¼ 5, H ¼ p=4 and S ¼ 1. The velocity increases with increasing K. Fig. 7 refers to the effects of K and Re upon the axial velocity at the crest section at a ¼ 0:01, / ¼ 0:4, M ¼ 5, H ¼ p=4 and S ¼ 1. As shown, the velocity is decreasing with increasing Re . Also it is clear that the velocity increases as K increases. We noticed that the effect of K becomes smaller as K increases and the maximum value for our study on the velocity U is at x ¼ 0:5. Fig. 8 represents the stream function w versus x at a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=4 and S ¼ 1. This indicates to the reader that the periodic boundary conditions (9) have been imposed. Fig. 9 shows the effect of S on the axial velocity U versus x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5 and H ¼ p=4. The value of the velocity U is decreasing as S increases. 5.2. Pressure gradient and pressure rise We notice from Eqs. (14a) and (14b) that the pressure gradient increases with increasing sin h. Also this contribution of sin h to the pressure gradient is
Fig. 11. Variation of pressure rise DPk on the center axis with Froude number Fr for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2 and S ¼ 1.
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
99
Fig. 12. Variation of pressure gradient distribution oP =ox with Froude number Fr for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and sin h ¼ 0:5.
found decreasing at all values of sin h as Fr becomes larger. Also we notice from Eq. (15) that the pressure rise increases with increasing sin h. Also this contribution of sin h to the pressure rise is decreasing at all values of sin h as Fr becomes larger (see Figs. 10 and 11). The effect of Fr on the pressure gradient distribution for the case that a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and sin h ¼ 0:5 is illustrated in Fig. 12. As shown, the pressure gradient distribution decreases as Fr increases. The effect of sin h on the
Fig. 13. Variation of pressure gradient distribution oP =ox with sin h for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and Fr ¼ 0:1.
100
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
Fig. 14. Variation of pressure gradient distribution oP =ox with suction parameter S for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, Fr ¼ 0:1 and sin h ¼ 0:5.
pressure gradient distribution is shown in Fig. 13 for the case that a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, S ¼ 1 and Fr ¼ 0:1. As sin h increases, the pressure gradient increases. Fig. 14 represents the effect of S on the pressure gradient distribution with x for a ¼ 0:01, / ¼ 0:4, Re ¼ 10, K ¼ 0:1, M ¼ 5, H ¼ p=2, Fr ¼ 0:1 and sin h ¼ 0:5. The value of pressure gradient increases with increasing S.
6. Conclusions This paper has presented a numerical study of an incompressible viscous fluid moving through a porous media (Brinkmain model) flow between two inclined wavy porous plates under the effect of a constant inclined magnetic field that makes an angle H with the vertical axis. The stream fu–vorticity formulation was chosen in conjunction with an ‘‘oblique lattice’’ coordinate system. The numerical method that was used has a second-order accuracy oðh2 ; k 2 Þ. The present approach is not restricted by any of the parameters appearing in the problem such as Reynolds, Froude and wave numbers, suction (or injection) and magnetic parameters, and amplitude ratio, but if the Reynolds numbers become high enough, the flow becomes turbulent and the present approach becomes invalid. The magnetic field and suction (or injection) have the same effect on the axial velocity component. Also the inclined angle of magnetic field and the permeability have the same effect on the axial velocity component. The effect of K becomes smaller as K increases. Physically,
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
101
this result can be achieved when the holes of the porous medium are very large so that the resistance of the medium may be neglected. The inclination angle of the wavy plates h and the Froude number Fr have the same effect on the pressure gradient and the pressure rise and they do not play any role in the fluid flow. The pressure gradient and pressure rise increase as Fr decreases and they increase with increasing sin h. Also this contribution of sin h to the pressure gradient and pressure rise is found decreasing at all values of sin h as the Froude number Fr becomes larger.
Appendix A
A1 ¼
1 ; k2
A0 ¼
2 A1 ð2 þ Ci;j Þ; h2i
A3 ¼ A1 ð1 þ Li;j Þ;
A4 ¼
A2 ¼
1 1 þ A1 ðCi;j 2Bi;j Þ; h2i 2
1 1 þ A1 ðCi;j þ 2Bi;j Þ; h2i 2
A5 ¼ A1 Li;j ;
where Bi;j ¼
i;j ; hi
Li;j ¼
i;j þ i1;j ; hi1
Si;j ¼
i;j þ i1;j ; hi
Ci;j ¼ Sij2 þ
hi1 Si;j 2B2i;j ; hi
and Ri;j ¼
i;j þ i1;j hi
1 ; K 1 B0 ¼ a0 þ Ni jnj jdj þ sgn n sgn d dRij ; 2 a 0 ¼ A0
B1 ¼ A1 þ Ni dfGðdÞ 1g; 1 B2 ¼ A2 þ Ni jnj sgn n sgn d dRij GðnÞ; 2 B3 ¼ A3 þ Ni dGðdÞ; B4 ¼ A 4 þ N i
1 jnj þ sgn n sgn d dRij ðGðnÞ 1Þ; 2
B5 ¼ A 5 ;
102
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
2 1 Hi;j ¼ Re M A0 þ 2 wi;j þ A1 wiþ1;j þ A2 þ 2 wi;jþ1 þ A3 wi1;j h hi i 1 2 2 þ A4 2 wi;j1 þ A5 wi1;j1 sin H þ ðwi1;j1 wi1;j Þ khi1 hi hi1 ðwi;jþ1 wi;j1 Þ Ei;j sin H cos H þ ½Ei;j cos2 H ; þ 2hi Re hi1 Ni ¼ ; Ei;j ¼ ð2i1;j þ hi1 Þðwiþ1;j 2wi;j þ wi;j1 Þ: 2khi 2h2i x ðx 6¼ 0Þ; 1 ðx P 0Þ GðxÞ ¼ and sgn x ¼ jxj 0 ðx < 0Þ 0 ðx ¼ 0Þ: ! 2ði;j þi1;j Þ Si1;j k c1;i ¼ ; c2;i ¼ ðc0;i þ c1;i Þ; Si;j þ 2khi1
1 S 2 i1;j 2i;j þ hi1 ðhi1 þ 2i1;j Þ c0;i ¼ 2 c1;i hi 2k A1 ð2i1;j þ 2i;j Þ þ 2 ; d0;j d2;j d4;j d0;j
2 S2;j þ B2;j ¼ 2 þ A1 þ B1;j L2;j ; R2;j hi
S2;j ; R2;j h1 B1;j 1 A1 h1 1 ¼ A1 ; d3;j ¼ 2 B1;j L2;j 1 ; hi h2 R2;j 2 h2 R2;j 1 A1 h1 1 ¼ 2 B1;j L2;j 1 þ ; hi 2 h2 R2;j 2 1 hN hN þ1 ¼ 2 þ A1 ðBN ;j 2SN ;j Þ þ BN ;j mN ;j 2 ; hN þ1 RN ;j hN þ1
¼ 2A1 d1;j
SN ;j ; RN ;j
d2;j ¼ A1
BN ;j ; RN ;j
d1;j ¼ 2A1
A1 hN hN 1 hN þ1 1 B m 1 þ ; N ;j N ;j h2N þ1 2 h2N þ1 hN RN ;j 1 A1 hN hN 1 hN þ1 1 d4;j ¼ 2 BN ;j mN ;j 2 1 ; hN þ1 2 hN þ1 hN RN ;j 2 h2 e0 ¼ 2 2A1 B1;mþ1 B1;mþ1 ; h1 h1 2 h2 h1 e2 ¼ 2 þ 2A1 B1;mþ1 B1;mþ1 þ ; h1 h1 h2
¼ d3;j
1
E.F. Elshehawey et al. / Appl. Math. Comput. 135 (2003) 85–103
e1 ¼ 2A1
103
h1 B1;mþ1 ; h2
and e3 ¼ 2kA1 þ e0 : References [1] E.F. Elshehawey, E.M. Elbarbary, N.S. Elgazery, Numerical study of a magneto-fluid motion through a porous medium between two wavyplates, Int. J. Comput. Fluid Dyn. 15 (2001) 177– 181. [2] J.A. Harris, R.L. Street, Numerical simulation of turbulent flow over a moving wavy boundary Norris and Reynolds extended, Phys. Fluids 6 (1994) 924–943. [3] C. Pozridikis, A study of peristaltic flows, J. Fluid Mech. 180 (1987) 515. [4] B.V. Rathish Kumar, K.B. Naidu, A numerical study of peristaltic flows, Comput. Fluids 24 (2) (1995) 161–176. [5] S. Selvarajan et al., A numerical study of flow through wavy-walled channels, Int. J. Numer. Meth. Fluids 26 (1998) 519–531. [6] K.D. Setphanhoff, I.J. Sobey, B.J. Bellhouse, On flow through furrowed channels, Part 2: Observed flow patterns, J. Fluid Mech. 96 (1980) 27–32. [7] I.J. Sobey, On flow through furrowed channels, Part 1: Calculated flow patterns, J. Fluid Mech. 96 (1980) 1–26. [8] S. Takabatake, K. Ayukawa, Numerical study of two-dimensional peristaltic flows, J. Fluid Mech. 122 (1982) 439. [9] S. Takabatake, K. Ayukawa, A. Mori, Peristaltic pumpingin circular cylindrical tubes: a numerical study of fluid transport and its efficiency, J. Fluid Mech. 193 (1988) 267.