New Study of the Stability of Fluid Flow in a Porous Channel under Effects of Magnetic Field and Radiation

Iman Hashim AL-Obeide, Alaa Abdul-Raheem Hammodat

College of Education for Pure Science, University of Mosul, Department of Mathematic, Mosul, Iraq.

DOI: 10.4236/oalib.1106282                

Abstract

In this research, we presented a study of the stability of flow systems and heat transfer in horizontal channels under the influence of the vertical magnetic field as well as thermal radiation where we built a mathematical model and put the partial differential equations that control the model and we analyze these equations after transformation into non-dimensional equations into two parts: the first unsteady state equations and second steady state equations, hence analyze the stability on the unsteady state equations. It’s noticed that the increase of Reynolds number Re causes the increase to unstable probability of the system. But the increase of Schmidt number has an effect on the system towards unstable, as well as the increase of grash of number makes the system stable. Finally, the increase of wave number k has a positive effect towards stable.

Keywords

Horizontal Channels, Unsteady State Equations, Reynolds Number and Stability

1. Introduction

The stability became very important and a study axis for many researchers in past last years, because of the industrial and technological development, for example, the existing solar system now depended on time condition which the planets move around the sun in a regular shape.

As known, if an additional small planet was entered into this system, the original state doesn’t shiver to an important degree, so, it can say, that the original state is stabilized. In the following, reference review for some last works in this field.

In 2000, Mahdi, F. analyzed the stability case for a sample of transfer the heat by convection and by connecting in porous area, and from the analysis, it turns out that the supposed hassle on the sample decays with passing the time, and the sample under the probability is always stable [1] .

In 2004, Mahdi, F. and Muthana, A. studied the stability for a system as a fluid between two parallel infinite slabs which are heated from down, considering that the heat transfer by the connection, convection and radiation, so it’s clear that the stability of the system depends on the ratio of the heat between the two slabs, also at the thermal expansion factor [2] .

In 2005, Hamsa, D. studied the stability of the Navier-stocks equation after disturbing this equation, and we find the regions where the flow is stable or unstable [3] .

In 2010, Osama, T. and Ahmed, M. did a study the natural convection inside a glass cavity, they found that the suitable ankle to be higher portability to isolate the outer medium from the inner medium [4] .

In 2012, Ala’a, A. and Ahmed, M. discussed the fluid flow matter in a horizontal channel under the effect of a vertical magnetic field on the level of the channel, when the slop ankle of the channel is: 0, 30, 60, 90 degrees, so he noted that the increasing and decreasing in Brickman’s values Fs, the wave number k, a wave number w effect in the stability of the system [5] .

In 2014, Ala’a, A. and Taghread, H. discussed the stability analysis in the glass cavity and this analysis was done by finding the self-values for the system which we are could find the growth of the disturbance or not, that after making the linear equations [6] .

In 2016, Mahantesh, M. and Shilpa, J. reported the stagnation point flow of Non-Newtonian fluid and heat transfer over a stretching/shrinking sheet in a porous medium and we discussed the numerical values of skin friction coefficient and local Nusselt [7] . In this paper, the stability of unsteady state solution of horizontal channel with the presence of magnetic field and radiation have been investigated and analyzed, it’s found that the parameters Re, Sc, Gr as well as k have a significant effect on the stability of the system.

2. Mathematical Formulation

Consider a fully-developed, steady laminar flow of the fluid in the horizontal channel, the distance between the walls of the channel is h apart. Choosing the coordinate system such that the x-axis in the direction of the flow, y-axis is measured perpendicular to the plane of the channel, whilst the z-axis is in the direction mutually orthogonal to the other two axes.

In the model under consideration, the magnetic field has a component $B_x$ induced along the channel in the direction of the flow, $B_z$ is zero and the component parallel to y-axis denoted by $B_0$ , the velocity u and v are zero at the edge, and $h,g,B_0,T_0,T_1$ are The distance between the two walls, gravitational acceleration, Boltzmann Number, Lower wall temperature, Upper wall temperature respectively as shown in Figure 1. Under these assumptions, the geometry and governing equations of the problem are:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (1)

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\frac{-1}{\rho }\frac{\partial p}{\partial x}+\upsilon \left[\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right]-\frac{\upsilon }{K}u$ (2)

$\begin{array}{l}\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\\ =\frac{-1}{\rho }\frac{\partial p}{\partial y}+\upsilon \left[\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right]+\frac{\sigma }{\rho }{B}_0^{2}v-g\beta \left(T-T_1\right)-g{\beta }^{*}\left(C-C_1\right)\end{array}$

$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{k^{*}}{\rho C_p}\left[\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right]-\frac{1}{\rho C_P}\left[\frac{\partial q_r}{\partial y}\right]+\frac{\mu }{\rho C_p}\left[{\left(\frac{\partial u}{\partial y}\right)}^2\right]$ (4)

$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=D\left[\frac{{\partial }^{2}C}{\partial x^2}+\frac{{\partial }^{2}C}{\partial y^2}\right]$ (5)

where u, v are the velocity components t, t is the time and $T,g,\beta ,{\beta }^{*},k^{*},K,\mu ,\rho ,C_p,D,B_0,\upsilon ,C,q_r,\sigma$ , are Temperature, gravitational acceleration, Thermal expansion Coefficient, Concentration Expansion, Absorption coefficient, Permeability of medium, Fluid Viscosity, Density, Specific heat, The mass diffusion coefficient, Boltzmann Number, Kinematic Fluid Viscosity, Concentration, Radiation flux, Electric Conductivity respectively.

With boundary conditions:

$\begin{array}{l}u=v=0,T=T_0,\frac{\partial T}{\partial y}=0,C=C_0,\frac{\partial C}{\partial y}\text{at}y=0\\ u=v=0,T=T_1,\frac{\partial T}{\partial y}=0,C=C_1,\frac{\partial C}{\partial y}\text{at}y=h\end{array}\}$ (6)

Figure 1. Mathematical model.

By using the Roseland approximations consider the radiative heat flux for optically thick fluid is given by [8] .

$q_r=\left(\frac{-4{\sigma }^{*}}{3k_1}\right)\frac{\partial T^4}{\partial y}$ (7)

where ${\sigma }^{*}$ is the Stefan-Boltzmann constant and $k_1$ is the mean absorption coefficient. Assume that the difference in temperature within the flow is sufficiently small such as that $T^4$ can be expressed as a linear function of the temperature, we expand $T^4$ in a Taylor’s series about $T_1$ and neglected higher order terms, thus [9] .

$T^4\cong 4T_1^{3}T-3T_1^4$ (8)

Hence the equation of energy Equation (4) becomes:

$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\frac{k^{*}}{\rho C_p}\frac{{\partial }^{2}T}{\partial x^2}+\left(\frac{16{\sigma }^{*}{T}_1^3}{3k_1\rho C_P}+\frac{k^{*}}{\rho C_p}\right)\frac{{\partial }^{2}T}{\partial y^2}+\frac{\mu }{\rho C_p}\left[{\left(\frac{\partial u}{\partial y}\right)}^2\right]$ (9)

Let us introduce the following similarity transformation [10]

$\begin{array}{l}U=\frac{uh}{\upsilon \sqrt{Gr}},V=\frac{vh}{\upsilon \sqrt{Gr}},\theta =\frac{T-T_1}{T_1-T_0},\varphi =\frac{C-C_1}{C_1-C_0},\\ p=P\rho u^2,\tau =\frac{t\upsilon \sqrt{Gr}}{h^2},X=\frac{x}{h},Y=\frac{y}{h}\end{array}$ (10)

And non-dimensional parameters:

$\begin{array}{l}M={\left(\frac{\sigma }{\mu }\right)}^{\frac{1}{2}}{B}_{0}h,\alpha =\frac{k^{*}}{\rho C_p},R^{*}=\frac{3k^{*}{k}_1+16{\sigma }^{*}{T}_1^3}{3\mu C_{\rho }{k}_1},\\ Da=\frac{K}{h^2},Re=\frac{hu}{\upsilon },Pr=\frac{\upsilon }{\alpha },Sc=\frac{\upsilon }{D},N=\frac{\upsilon }{h^2\Delta T},\\ Gr=\frac{g\beta h^3\left(T_1-T_0\right)}{{\upsilon }^2},G{r}^{*}=\frac{g{\beta }^{*}{h}^3\left(C_1-C_0\right)}{{\upsilon }^2},\epsilon =\frac{\upsilon }{C_p}\end{array}\}$ (11)

where M is Hartmann number, $\alpha$ is Thermal diffusivity, $R^{*}$ is Radiation coefficient, N is new physical quantity, Sc is Schmidt number, Da is Darcy number, Re is Reynolds number and Pr is Prandtle number, Gr is Gratshof number for heat transfer, $G{r}^{\text{*}}$ is Gratshof number for mass transfer, $\epsilon$ is dispersion parameter.

The above Equations (10), (11) reduce the Equations (1), (2), (3), (5), (9) into the following system of non-dimensional equations:

$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$ (12)

$\frac{\partial U}{\partial \tau }+U\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial Y}=\frac{-{\left(Re\right)}^2}{Gr}\frac{\partial P}{\partial X}+\frac{1}{\sqrt{Gr}}\left[\frac{{\partial }^{2}U}{\partial X^2}+\frac{{\partial }^{2}U}{\partial Y^2}\right]-\frac{1}{Da\sqrt{Gr}}U$ (13)

$\begin{array}{l}\frac{\partial V}{\partial \tau }+U\frac{\partial V}{\partial X}+V\frac{\partial V}{\partial Y}\\ =\frac{-{\left(Re\right)}^2}{Gr}\frac{\partial P}{\partial Y}+\frac{1}{\sqrt{Gr}}\left[\frac{{\partial }^{2}V}{\partial X^2}+\frac{{\partial }^{2}V}{\partial Y^2}\right]+\left(\frac{M^2}{\sqrt{Gr}}\right)V-\theta -\left(\frac{G{r}^{*}}{Gr}\right)\varphi \end{array}$ (14)

$\frac{\partial \theta }{\partial \tau }+U\frac{\partial \theta }{\partial X}+V \frac{\partial \theta }{\partial Y}=\left(\frac{1}{Pr\sqrt{Gr}}\right)\frac{{\partial }^2\theta }{\partial X^2}+\left(\frac{R^{*}}{\sqrt{Gr}}\right)\frac{{\partial }^2\theta }{\partial Y^2}+ \epsilon N\sqrt{Gr}\left[{\left(\frac{\partial U}{\partial Y}\right)}^2\right]$ (15)

$\frac{\partial \varphi }{\partial \tau }+U\frac{\partial \varphi }{\partial X}+V\frac{\partial \varphi }{\partial Y}=\frac{1}{Sc\sqrt{Gr}}\left[\frac{{\partial }^2\varphi }{\partial X^2}+\frac{{\partial }^2\varphi }{\partial Y^2}\right]$ (16)

3. Fourier Mode Stability Analysis

Assume that the solution of Equations (12), (13), (14), (15) and (16) can be written in the form [11] .

$\begin{array}{l}U=U_1\left(x,y\right)+U_2\left(x,y,t\right)\\ V=V_1\left(x,y\right)+V_2\left(x,y,t\right)\\ \theta ={\theta }_1\left(x,y\right)+{\theta }_2\left(x,y,t\right)\\ \varphi ={\varphi }_1\left(x,y\right)+{\varphi }_2\left(x,y,t\right)\\ P=P_1\left(x,y\right)+P_2\left(x,y,t\right)\end{array}\}$ (17)

where $U_1,V_1,{\theta }_1,{\varphi }_1,P_1$ are the steady state solution and $U_1,V_1,{\theta }_1,{\varphi }_1,P_1$ are the disturbance.

Substituting Equation (17) into Equations (12), (13), (14), (15) and (16), with its boundary conditions, we get the following equations:

$\frac{\partial U_2}{\partial X}+\frac{\partial V_2}{\partial Y}=0$ (18)

$\frac{\partial U_2}{\partial \tau }=\frac{-{\left(Re\right)}^2}{Gr}\frac{\partial P_2}{\partial X}+\frac{1}{\sqrt{Gr}}\left[\frac{{\partial }^2{U}_2}{\partial X^2}+\frac{{\partial }^2{U}_2}{\partial Y^2}\right]-\frac{1}{Da\sqrt{Gr}}{U}_2$ (19)

$\frac{\partial V_2}{\partial \tau }=\frac{-{\left(Re\right)}^2}{Gr}\frac{\partial P_2}{\partial Y}+\frac{1}{\sqrt{Gr}}\left[\frac{{\partial }^2{V}_2}{\partial X^2}+\frac{{\partial }^2{V}_2}{\partial Y^2}\right]+\left(\frac{M^2}{\sqrt{Gr}}\right)V_2-{\theta }_2-\left(\frac{G{r}^{*}}{Gr}\right){\varphi }_2$ (20)

$\frac{\partial {\theta }_2}{\partial \tau }=\left(\frac{1}{Pr\sqrt{Gr}}\right)\frac{{\partial }^2{\theta }_2}{\partial X^2}+\left(\frac{R^{*}}{\sqrt{Gr}}\right)\frac{{\partial }^2{\theta }_2}{\partial Y^2}+\epsilon N\sqrt{Gr}\left[2\left(\frac{\partial U_1}{\partial Y}\right)\left(\frac{\partial U_2}{\partial Y}\right)\right]$ (21)

$\frac{\partial {\varphi }_2}{\partial \tau }=\frac{1}{Sc\sqrt{Gr}}\left[\frac{{\partial }^2{\varphi }_2}{\partial X^2}+\frac{{\partial }^2{\varphi }_2}{\partial Y^2}\right]$ (22)

With the boundary conditions:

$\begin{array}{l}{U}_2=V_2=0,{\theta }_2=0,\frac{\partial {\theta }_2}{\partial y}=0\text{at}y=0,1\\ U_2=V_2=0,{\varphi }_2=0,\frac{\partial {\varphi }_2}{\partial y}=0\text{at}y=0,1\end{array}\}$ (23)

4. Stability Analysis in the Case of the Variable Amplitude

To solve the linearized system (or to analyze the stability) and because the coefficient in the differential equations is independent of the attempt to find the solution of the form [5] :

$\begin{array}{l}{U}_2=U\left(y\right){\text{e}}^{ikx}{\text{e}}^{at}\\ V_2=V\left(y\right){\text{e}}^{ikx}{\text{e}}^{at}\\ P_2=P\left(y\right){\text{e}}^{ikx}{\text{e}}^{at}\\ {\theta }_2=\theta \left(y\right){\text{e}}^{ikx}{\text{e}}^{at}\\ {\varphi }_2=\varphi \left(y\right){\text{e}}^{ikx}{\text{e}}^{at}\end{array}\}$ (24)

where $\theta \left(y\right),\varphi \left(y\right),P\left(y\right),V\left(y\right),U\left(y\right)$ are the amplitude functions, k is wave number in the direction of x, and a is the complex number which has the form $a=a_1+i{a}_2,a_1,a_2\in R$ is speed number, when $a_1>0$ the system is unstable while $a_1<0$ , the system is stable.

From Equations (18), (19), (20), (21), (22) and (24), we get:

$\left(ikU\left(y\right)+V^{\prime }\left(y\right)\right){\text{e}}^{ikx}{\text{e}}^{at}=0$ (25)

$\begin{array}{l}(\left(a+\frac{k^2}{\sqrt{Gr}}+\frac{1}{Da\sqrt{Gr}}\right)U\left(y\right)\\ +\left(\frac{ik{\left(Re\right)}^2}{Gr}\right)P\left(y\right)-\frac{1}{\sqrt{Gr}}{U}^{″}\left(y\right)){\text{e}}^{ikx}{\text{e}}^{at}=0\end{array}$ (26)

$\begin{array}{l}(\left(a+\frac{k^2}{\sqrt{Gr}}-\frac{M^2}{\sqrt{Gr}}\right)V\left(y\right)+\frac{{\left(Re\right)}^2}{Gr}{P}^{\prime }\left(y\right)\\ -\frac{1}{\sqrt{Gr}}{V}^{″}\left(y\right)+\theta \left(y\right)+\left(\frac{G{r}^{*}}{Gr}\right)\varphi \left(y\right))){\text{e}}^{ikx}{\text{e}}^{at}=0\end{array}$ (27)

$\left(\left(a+\frac{k^2}{Pr\sqrt{Gr}}\right)\theta \left(y\right)-\left(\frac{R^{*}}{\sqrt{Gr}}\right){\theta }^{″}\left(y\right)\right){\text{e}}^{ikx}{\text{e}}^{at}=0$ (28)

$\left(\left(a+\frac{k^2}{Sc\sqrt{Gr}}\right)\varphi \left(y\right)-\frac{1}{Sc\sqrt{Gr}}{\varphi }^{″}\left(y\right)\right){\text{e}}^{ikx}{\text{e}}^{at}=0$ (29)

Since ${\text{e}}^{ikx}{\text{e}}^{at}\ne 0$ , then:

$\left(ikU\left(y\right)+V^{\prime }\left(y\right)\right)=0$ (30)

$\left(\left(a+\frac{k^2}{\sqrt{Gr}}+\frac{1}{Da\sqrt{Gr}}\right)U\left(y\right)+\left(\frac{ik{\left(Re\right)}^2}{Gr}\right)P\left(y\right)-\frac{1}{\sqrt{Gr}}{U}^{″}\left(y\right)\right)=0$ (31)

$\begin{array}{l}(\left(a+\frac{k^2}{\sqrt{Gr}}-\frac{M^2}{\sqrt{Gr}}\right)V\left(y\right)+\frac{{\left(Re\right)}^2}{Gr}{P}^{\prime }\left(y\right)\\ -\frac{1}{\sqrt{Gr}}{V}^{″}\left(y\right)+\theta \left(y\right)+\left(\frac{G{r}^{*}}{Gr}\right)\varphi \left(y\right))\end{array}$ (32)

$\left(\left(a+\frac{k^2}{Pr\sqrt{Gr}}\right)\theta \left(y\right)-\left(\frac{R^{*}}{\sqrt{Gr}}\right){\theta }^{″}\left(y\right)\right)=0$ (33)

$\left(\left(a+\frac{k^2}{Sc\sqrt{Gr}}\right)\varphi \left(y\right)-\left(\frac{1}{Sc\sqrt{Gr}}\right){\varphi }^{″}\left(y\right)\right)=0$ (34)

Hence, we get the following system:

$\begin{array}{l}{U}^{\prime }\left(y\right)=H\left(y\right)\\ H^{\prime }\left(y\right)=\left(a\sqrt{Gr}+k^2+\frac{1}{Da}\right)U\left(y\right)+\left(\frac{ik{\left(Re\right)}^2}{\sqrt{Gr}}\right)P\left(y\right)\\ V^{\prime }\left(y\right)=-ikU\left(y\right)\\ P^{\prime }\left(y\right)=-\left(\frac{aGr+k^2\sqrt{Gr}-M^2\sqrt{Gr}}{{\left(Re\right)}^2}\right)V\left(y\right)-\left(\frac{ik\sqrt{Gr}}{{\left(Re\right)}^2}\right)H\left(y\right)\\ -\left(\frac{Gr}{{\left(Re\right)}^2}\right)\theta \left(y\right)-\left(\frac{G{r}^{*}}{{\left(Re\right)}^2}\right)\varphi \left(y\right)\\ {\theta }^{\prime }\left(y\right)=Q\left(y\right)\\ Q^{\prime }\left(y\right)=\frac{\left(aPr\sqrt{Gr}+k^2\right)}{Pr{R}^{*}}\theta \left(y\right)\\ {\varphi }^{\prime }\left(y\right)=S\left(y\right)\\ S^{″}\left(y\right)=\left(aSc\sqrt{Gr}+k^2\right)\varphi \left(y\right)\end{array}\}$

And $\Omega$ system matrix of transaction such that

$A=\left(a\sqrt{Gr}+k^2+\frac{1}{Da}\right),B=\frac{ik{\left(Re\right)}^2}{\sqrt{Gr}},C=-ik,D=-\left(\frac{ik\sqrt{Gr}}{{\left(Re\right)}^2}\right)$

$E=-\left(\frac{aGr+k^2\sqrt{Gr}-M^2\sqrt{Gr}}{{\left(Re\right)}^2}\right),F=-\left(\frac{Gr}{{\left(Re\right)}^2}\right),G=-\left(\frac{G{r}^{*}}{{\left(Re\right)}^2}\right)$

$H=-\left(\frac{aPr\sqrt{Gr}+k^2}{Pr{R}^{*}}\right),W=-\left(aSc\sqrt{Gr}+k^2\right)$

By using $\left|\Omega -\delta I\right|=0$ we get the following equation:

$\begin{array}{c}f\left(\lambda \right)={\delta }^8-\left(DB+A+W+H\right){\delta }^6\\ +\left(WBD+WA-CBE+HW+BDH+AH\right){\delta }^4+\cdots \\ +\left(BCEH+WBCE-AWH-WHBD\right){\delta }^2-WHBCE\\ =0\end{array}$ (36)

where

Now, we solve Equation (36) numerically using (Maple 11) [12] , to find the roots of these equations as shown in Figures 2-5.

Figure 2. Effect of Reynolds number (Re = 250, 1000, 4000).

Figure 3. Effect of Schmidt number (Sc = 0.1, 0.3, 0.5).

Figure 4. Effect of Grash of number (Gr = 0.1, 0.5, 1).

Figure 5. Effect of wave number (k = 2, 4, 8).

Acknowledgements

The authors is grateful to Assist. Prof. Dr. Ala’a Abdul-Raheem Ahmed Hammodat for his valuable remarks.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Mosa, F. (2000) On the Stability of Heat Transfer through MHD Porous Medium. AL-Rafidain Journal of Computer Sciences and Mathematics, 11, 90-94.
[2]	Mosa, M.F. and Ali, A.M. (2004) Stability Analysis for Fluid Flow between Two Infinite Parallel Plates I. AL-Rafidain Journal of Computer Sciences and Mathematics, 1, 8-19. https://doi.org/10.33899/csmj.2004.164093
[3]	Saleem, H.D. (2005) Stability of Liquid Film with Negligible Viscosity. AL-Rafidain Journal of Computer Sciences and Mathematics, 2, 11-19. https://doi.org/10.33899/csmj.2005.164063
[4]	Osama, T. and Ahmed, M. (2010) Heat Transfer in Glazing Cavities. M.Sc. Thesis, University of Mosul, Mosul.
[5]	Ala’a, A. and Ahmed, M. (2012) A Theoretical Study of Stability of Flowing and Heat Transfer Systems in Vessels and Channels. Ph.D. Thesis, University of Mosul, Mosul.
[6]	Hammodat, A.A. and Shuker, T.H. (2014) Stability Analysis of a Fluid in Horizontal and Oblique Glass Cavitis. AL-Rafidain Journal of Computer Sciences and Mathematics, 11, 89-97. https://doi.org/10.33899/csmj.2014.163758
[7]	Namdeppanavar, M.M. and Shilpa, J.M. (2016) Stagnation Point Flow of Non- Newtonian Fluidand Heat Transfer over Stretching/Shrinking Sheet. Chemical and Process Engineering Research, 46, 27-34.
[8]	Brewster, M. (1972) Thermal Radiative Transfer Properties. John Wiley and Sons, Hoboken.
[9]	Bahattcharyya, K. (2013) MHD Stagnation-Point Flow of Casson Fluid and HeatTransfer over a Stretching Sheet with Thermal Radiation. Journal of Thermodynamics, 2013, Article ID: 169674. https://doi.org/10.1155/2013/169674
[10]	Cookey, C. and Omubo-Pepple, V. (2010) On Steady Hydro Genetic Flow of a Radiating Viscous Fluid through a Horizontal Channal in a Porous Media. AMERICAN Journal of Science and Industrial Research, 1, 203-208.
[11]	Logan, J. (1987) Applied Mathematics. John Wiley and Sons, Hoboken.
[12]	Martha, L. and James, P. (2005) Maple by Example. 3rd Edition, Elsevier Academic Press, Cambridge.