# Finite Element Solution of an Unsteady MHD Flow through Porous Medium between Two Parallel Flat Plates

Authors: AbdelLatif Sa’adAldin; Naji Qatanani
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

Finite element solution of unsteady magnetohydrodynamics (MHD) flow of an electrically conducting, incompressible viscous fluid past through porous medium between two parallel plates is presented in the presence of a transverse magnetic field and Hall effect. The results obtained from some test cases are then compared with previous published work using the finite difference method (FDM). Numerical examples show that the finite element method (FEM) gives more accurate results in comparison with the finite difference method (FDM).


---

## Body

## 1. Introduction

Theoretical study of magnetohydrodynamics (MHD) flow problems are frequently encountered in cooling systems of nuclear reactors, MHD generators, blood flow measurements, pumps, and accelerators.Due to coupling of the equations for electrodynamics and fluid mechanics, exact solution is possible only for some simple situations. By using several numerical techniques, such as finite element method (FEM), finite volume method (FVM), and boundary element method (BEM), approximate solution for the MHD flow problems can be obtained.Gapta and Singh [1] obtained the exact solutions for unsteady flow in some special cases. Ram and Mishra [2] investigated the unsteady flow through magnetohydrodynamic porous media. Singh and Lal [3] studied the FEM solution of time-dependent MHD flow equations. Ram and Jain [4] have discussed MHD free convective flow through a porous medium in a rotating fluid. Reddy and Bathaiah [5] have analyzed the Hall effects on MHD flow through a porous straight channel. Lee and Dulikravich [6] proposed FDM scheme for the 3-dimensional unsteady MHD flow together with temperature field. Sheu and Lin [7] presented a convection-diffusion-reaction model for solving the unsteady MHD flow using a FDM scheme. The stabilized FEM for solution of the 3-dimensional time-dependent MHD flow equations was given by Ben Salah and et al. [8]. Chauhan and Rastogi [9] have studied the Hall effects on MHD slip flow and heat transfer through a porous medium over an accelerated plate in a rotating system. Saha and Chakrabarti [10] have investigated the impact of magnetic field strength on magnetic fluid flow through a channel. Moniem and Hassanin [11] have developed a solution of MHD flow past a vertical porous plate through a porous medium under oscillatory suction. Sa’adAldin and Qatanani [12] have studied the unsteady MHD flow through two parallel porous flat plates. Sivaiah and Srinivasa-Raju [13] have discussed the finite element solution of heat and mass transfer in MHD flow of a viscous fluid past a vertical plate under oscillatory suction velocity. Yuksel and Ingram [14] have investigated the numerical analysis of a finite element method, Crank-Nicolson discretization for MHD flows at small magnetic Reynold number. Bég et al. [15] have developed a finite element and network electrical simulation of rotating magnetofluid flow in nonlinear porous media with inclined magnetic field and Hall currents. Sa’ad Aldin and Qatanani [16] have studied the analytical and finite difference methods for solving unsteady MHD flow through porous medium between two parallel flat plates.In this work, the finite element solution for the unsteady magnetohydrodynamics (MHD) flow of an electrically conducting, incompressible viscous fluid past through porous medium between two parallel plates in the presence of a transverse magnetic field and Hall effect is considered. A comparison study has been carried out between the finite difference and the finite element solutions. A case study is analyzed with both the finite element method (FEM) and the finite difference method (FDM), namely, the implicit scheme presented in [16]. It was found that the finite element method (FEM) is more accurate for solving these type of problems.


## 2. Formulation of the Problem

We consider an unsteady flow of an electrically conducting, incompressible viscous fluid past through porous medium between two parallel plates with Hall effect. Let thex-axis be taken along the plates and y-axis be normal to the plates. The fluid is subjected to a constant transverse magnetic field of strength B0 in the y direction, with the flow being considered in the x direction, as illustrated in Figure 1. The governing equations for the unsteady, viscous incompressible flow of an electrically conducting fluid for the Brinkman-extended Darcy model are as follows[2]:Figure 1
Schematic diagram of the system.Equation of continuity is(1)∇·q=0.Equation of motion is(2)∂q∂t+q·∇q=-1ρ∇p+μρ∇2q-μρqk+1ρJ×B.General Ohm’s law is(3)J+ωτB0J×B=σE+q×B+1ρene∇Pe.Gauss’s law of magnetism is(4)∇·B=0,where q is the velocity vector, ρ is the fluid density, p is the pressure, J is the current density, B is the magnetic vector, μ is the coefficient of viscosity, σ is the electrical conductivity, k is the permeability of the medium, ω is the electron frequency, τ is the electron collision time, ρe is the electric charge, ne is the number density of electron, Pe is the electron pressure, and E is the electric field.We assumeE to be negligible and the magnetic Reynold’s number is small so that magnetic induction effect is ignored. Moreover, in the absence of pressure gradient, the ion-slip effects and electron pressure gradient, we have (5)J=σq×B-mB0J×B,J=jx,jy,jz,q=u,0,0,B=0,B0,0,(6)jx=mjz,(7)jy=0,(8)jz=σB0u-mjx.Solving (6) and (8), we have(9)jx=σB0mu1+m2,jz=σB0u1+m2,1ρJ×B=-1ρjzB0.As the plates are infinite, there is nox dependence. Consequently, (2) and (3) take the following form:(10)∂u∂t=-1ρ∂p∂x+ν∂2u∂y2-σB02uρ1+m2-νuk,0=-1ρ∂p∂y,where u is the axial velocity, ν is the kinematic viscosity, and m=ωτ is the Hall parameter. The initial and boundary conditions are given by(11)u=0,t≤0,u=0,y=±h,t>0.Upon introducing the nondimensional quantities,(12)Y=yh,T=νth2,U=uV,M2=σB02h2μ,K=kh2,P=phμV,X=xh,where M is the Hartman number, K is the Darcy parameter, and V is the mean velocity of the fluid. Then, the partial differential equations (10) together with the initial and boundary conditions (11) become(13)∂U∂T=-∂P∂X+∂2U∂Y2-M21+m2+1KU,(14)0=∂P∂Y,subject to the initial and boundary conditions:(15)U=0,T≤0,U=0,Y=±1,T>0.In virtue of (14), the pressure is independent of Y; then it is a function of T only. In this case, we can take the pressure gradient as a constant quantity; that is,(16)∂P∂X=-P0,where P0>0; thus (13) becomes(17)∂U∂T=P0+∂2U∂Y2-aUsubject to the initial and boundary conditions:(18)U=0,T≤0,U=0,Y=±1,T>0,where a=M2/1+m2+1/K.


## 3. Finite Element Method

### 3.1. Variational Formulations and Galerkin Approximation

The dimensionless partial differential equation (17) subject to the initial and boundary conditions (18) is solved by weighted residual Galerkin finite element method. The standard approach to deriving a Galerkin scheme is to multiply both sides of (17) by a test function ψ∈H01[-1,1] and integrate over the domain(19)∫-11∂U∂Tψdx=∫-11P0ψdx+∫-11∂2U∂Y2ψdx-a∫-11Uψdx,where(20)H01≔ψ∈L2-1,1,∂ψ∂Y∈L2-1,1,ψ-1=ψ1=0.Integrating by parts, we obtain(21)∂U∂T,ψ=P0,ψ-∂U∂Y,∂ψ∂Y-aU,ψ,where ·,· denotes the L2-inner product and U,ψ=0 for T≤0.We shall approximate the solution of (21) by assuming that U and ψ lie in finite dimensional subspace of H01[-1,1] for each T. Let φi∈H01[-1,1] for i=1,…,W and assume that the set φ1,…,φW is linearly independent. Further, let Sh be a partition of the interval [-1,1] into subintervals -1=Y0<Y1<⋯<YW<YW+1=1. Now we define the finite dimensional space Vh spanned by φ1,…,φW as(22)Vh≔v∈H01-1,1,vis piecewise linear function onSh.To this end, the approximate solution uh is(23)uh=∑i=1WuiTφiY.Inserting (23) into (21) and selecting as trial function ψ the basis function of uh, we obtain a system of ODEs:(24)Bu˙=P-Au-Bau,where u˙=(∂u1/∂T,…,∂uW/∂T)T and u=(u1,…,uw)T. Here A is the Stiffness matrix and B is the Mas matrix defined, respectively, as(25)A=Aij=∫-11φi′Yφj′YdY,i,j=1,2,…,w,B=Bij=∫-11φiYφjYdY,i,j=1,2,…,w,P=Pi=∫-11P0φiYdY,i=1,2,…,w,with φi(Yi)=δij being the usual finite element basis corresponding to the partition Sh. Thus, to compute the entries of the Stiffness matrix A, next, we need to determine φi′(Y).Here(26)φiY=Y-Yi-1hi,Yi-1≤Y≤Yi,Yi+1-Yhi+1,Yi≤Y≤Yi+1,0,elsewhere,and then(27)φi′Y=1hi,Yi-1≤Y≤Yi,-1hi+1,Yi≤Y≤Yi+1,0,elsewhere,where hi=Yi-Yi-1.Using (26) and (27), with a uniform mesh hi=h, we get(28)A=1h2-10⋯0-1⋱⋱⋱⋮0⋱⋱⋱0⋮⋱⋱⋱-10⋯0-12,(29)B=h23160⋯016⋱⋱⋱⋮0⋱⋱⋱0⋮⋱⋱⋱160⋯01623,(30)P=hP01⋮⋮⋮1.Scheme (24) is called semidiscretization, since uh is still a continuous function of T [17].


### 3.2. Time Stepping

In this section, we consider the semidiscretization in time. We first discretize the time interval(0,T) into a uniform grid with size k=T/N. Approximating the derivative in (24) at time level Tn by the Crank-Nicolson scheme with u0=0, we have (31)Bun-un-1k=P-Aun+un-12-Baun+un-12.Then, we rewrite the Crank-Nicolson method as(32)B+k2A+Baun=B-k2A+Baun-1+kP.Thus, we have the full discretization which is simply a combination of discretization in space and time:(33)U=uhn.


## 3.1. Variational Formulations and Galerkin Approximation

The dimensionless partial differential equation (17) subject to the initial and boundary conditions (18) is solved by weighted residual Galerkin finite element method. The standard approach to deriving a Galerkin scheme is to multiply both sides of (17) by a test function ψ∈H01[-1,1] and integrate over the domain(19)∫-11∂U∂Tψdx=∫-11P0ψdx+∫-11∂2U∂Y2ψdx-a∫-11Uψdx,where(20)H01≔ψ∈L2-1,1,∂ψ∂Y∈L2-1,1,ψ-1=ψ1=0.Integrating by parts, we obtain(21)∂U∂T,ψ=P0,ψ-∂U∂Y,∂ψ∂Y-aU,ψ,where ·,· denotes the L2-inner product and U,ψ=0 for T≤0.We shall approximate the solution of (21) by assuming that U and ψ lie in finite dimensional subspace of H01[-1,1] for each T. Let φi∈H01[-1,1] for i=1,…,W and assume that the set φ1,…,φW is linearly independent. Further, let Sh be a partition of the interval [-1,1] into subintervals -1=Y0<Y1<⋯<YW<YW+1=1. Now we define the finite dimensional space Vh spanned by φ1,…,φW as(22)Vh≔v∈H01-1,1,vis piecewise linear function onSh.To this end, the approximate solution uh is(23)uh=∑i=1WuiTφiY.Inserting (23) into (21) and selecting as trial function ψ the basis function of uh, we obtain a system of ODEs:(24)Bu˙=P-Au-Bau,where u˙=(∂u1/∂T,…,∂uW/∂T)T and u=(u1,…,uw)T. Here A is the Stiffness matrix and B is the Mas matrix defined, respectively, as(25)A=Aij=∫-11φi′Yφj′YdY,i,j=1,2,…,w,B=Bij=∫-11φiYφjYdY,i,j=1,2,…,w,P=Pi=∫-11P0φiYdY,i=1,2,…,w,with φi(Yi)=δij being the usual finite element basis corresponding to the partition Sh. Thus, to compute the entries of the Stiffness matrix A, next, we need to determine φi′(Y).Here(26)φiY=Y-Yi-1hi,Yi-1≤Y≤Yi,Yi+1-Yhi+1,Yi≤Y≤Yi+1,0,elsewhere,and then(27)φi′Y=1hi,Yi-1≤Y≤Yi,-1hi+1,Yi≤Y≤Yi+1,0,elsewhere,where hi=Yi-Yi-1.Using (26) and (27), with a uniform mesh hi=h, we get(28)A=1h2-10⋯0-1⋱⋱⋱⋮0⋱⋱⋱0⋮⋱⋱⋱-10⋯0-12,(29)B=h23160⋯016⋱⋱⋱⋮0⋱⋱⋱0⋮⋱⋱⋱160⋯01623,(30)P=hP01⋮⋮⋮1.Scheme (24) is called semidiscretization, since uh is still a continuous function of T [17].


## 3.2. Time Stepping

In this section, we consider the semidiscretization in time. We first discretize the time interval(0,T) into a uniform grid with size k=T/N. Approximating the derivative in (24) at time level Tn by the Crank-Nicolson scheme with u0=0, we have (31)Bun-un-1k=P-Aun+un-12-Baun+un-12.Then, we rewrite the Crank-Nicolson method as(32)B+k2A+Baun=B-k2A+Baun-1+kP.Thus, we have the full discretization which is simply a combination of discretization in space and time:(33)U=uhn.


## 4. Numerical Results and Discussion

To show the efficiency of the FEM described in the previous parts and to draw a comparison with the FDM, we present some examples. These tests are chosen such that there exist analytical solutions for them to give an obvious overview of the methods presented in this work.Numerical Example  1.
As an application to the FEM and FDM, we consider the following test case with the valuesM=1, m=1, K=0.1, P0=1, T=0.25 fixed, and Y∈[-1,1]. Next the computed matrices A,B, and P given in (28), (29), and (30), respectively, are used in (32) to obtain the velocity U(Y,T).
Table1 compares the exact values for the velocity U(Y,T) with both the FEM and the FDM values (for more details on the exact and FDM solutions, see [16]). A further comparison between the exact, FEM, and FDM values for the velocity can be observed in Figure 2. A plot of the absolute error that resulted from the FEM can be seen in Figure 3. Figure 4 compares the absolute errors obtained from the FEM and the FDM solutions.Table 1
The exact, (FEM) and (FDM) solutions of the velocityU.
Y
Exact solutionUE
FEM solutionUFE
|
U
E
-
U
F
E
| by FEM
FDM solutionUFD
|
U
E
-
U
F
D | by FDM
-
0.8
0.0439256214739
0.0439263413482
7.1987427380237
×
1
0
-
7
0.0439284690437
0.2847569757748
×
1
0
-
5
-
0.6
0.0662966816784
0.0662977446546
1.0629762063479
×
1
0
-
6
0.0662354201183
0.6126156005674
×
1
0
-
4
-
0.4
0.0773731593814
0.0773743987214
1.2393400675991
×
1
0
-
6
0.0772434172904
0.1297420909940
×
1
0
-
3
-
0.2
0.0824053505275
0.0824066768827
1.3263552283981
×
1
0
-
6
0.0822283020686
0.1770484588176
×
1
0
-
3
0
0.0838631333124
0.0838644864314
1.3531190140581
×
1
0
-
6
0.0836692897549
0.1938435575730
×
1
0
-
3
0.2
0.0824792484052
0.0824805760964
1.3276912429471
×
1
0
-
6
0.0823013886500
0.1778597551711
×
1
0
-
3
0.4
0.0775578709141
0.0775591133358
1.2424216914302
×
1
0
-
6
0.0774266269162
0.1312439979023
×
1
0
-
3
0.6
0.0666815959272
0.0666826648556
1.0689283097276
×
1
0
-
6
0.0666185173810
0.6307854626647
×
1
0
-
4
0.8
0.0446904326399
0.0446911643550
7.3171502338459
×
1
0
-
7
0.0446921159260
0.1683286084617
×
1
0
-
5Figure 2
The exact, FEM, and FDM values for the velocity.Figure 3
The absolute error resulted from the FEM approximation.Figure 4
The absolute error resulted from the FEM and FDM approximations.Numerical Example  2.
As for another test case, we takeM=2,  m=1,  K=0.3,  P0=1,  Y=0.5 fixed, and T∈[0,1]. Figure 5 compares the exact, FEM, and FDM values for the velocity. Figure 6 presents a plot of the absolute error that resulted from the FEM. A comparison between the absolute errors obtained from the FEM and the FDM solutions can be seen in Figure 7.Figure 5
The exact, FEM, and FDM values for velocity case 2.Figure 6
The absolute error resulted from FEM approximation case 2.Figure 7
The absolute error that resulted from FEM and FDM approximations case 2.


## 5. Conclusions

MHD flow problems, which have a very important place in physics and engineering, are usually hard to solve analytically. Therefore, it is required to obtain approximate solutions using computational methods. In this work, the problem of unsteady MHD flow through porous medium in the presence of magnetic field between two parallel flat plates has been investigated and solved using the FEM.A comparison between FFM and FDM has been carried out. The exact results and the numerical results using the FEM have shown to be in closed agreement. This can clearly be seen in Figures2, 3, 5, and 6. The results of the numerical examples indicate that the FFM is more accurate than the FDM (see Figures 4 and 7). This asserts the ability and reliability of the FEM for solving these types of problems.



---

*Source: 6856470-2017-06-14.xml*