Unsteady Magnetohydrodynamic Heat Transfer in a Semi-Infinite Porous Medium with Thermal Radiation Flux : Analytical and Numerical Study

1 Aerospace Engineering, Department of Engineering and Mathematics, Sheaf Building, Sheffield Hallam University, Sheffield S1 1WB, UK 2 Thermal Engineering and Fluids Department, Technical University of Cartagena, Campus Muralla del Mar, 30202 Cartagena, Spain 3 Magnetohydrodynamics, Applied Mathematics Program, Department of Mathematics, Narajole Raj College, P.O. Narajole, Midnapore 721 211, WB, India 4 Institute for Advanced Studies, Tehran 14456-63543, Iran


Introduction
Transient MHD magnetohydrodynamic flows with and without heat transfer in electricallyconducting fluids have attracted substantial interest in the context of metallurgical fluid dynamics, re-entry aerothermodynamics, astronautics, geophysics, nuclear engineering, and applied mathematics.An early study was presented by Carrier and Greenspan 1 who considered unsteady hydromagnetic flows past a semi-infinite flat plate moving impulsively in its own plane.Gupta 2 considered unsteady magneto-convection under buoyancy forces.Singer 3 further assessed the unsteady free convection heat transfer with magnetohydrodynamic effects in a channel regime.Pop 4 reported on transient buoyancy-driven convective hydromagnetics from a vertical surface.Yu and Yang 5 investigated the influence of channel wall conductance on hydromagnetic convection.Rao 6 analyzed the unsteady magnetohydrodynamic convection heat transfer past an infinite plane.Further excellent studies of unsteady free convection magnetohydrodynamic flows were reported by Antimirov and Kolyshkin 7 for a vertical pipe and Rajaram and Yu for a parallel-plate channel 8 .Tokis 9 used Laplace transforms to analyze the threedimensional free-convection hydromagnetic flow near an infinite vertical plate moving in a rotating fluid when the plate temperature undergoes a thermal transient.The influence of oscillatory pressure gradient on transient rotating hydromagnetic flow was considered by Ghosh 10 .Other transient MHD studies include the papers by Sacheti et al. 11 , Attia 12 who included viscosity variation effects, Al-Nimr and Alkam 13 who considered open-ended vertical annuli, and Takhar et al. 14 who employed a numerical method to study flat-plate magnetohydrodynamic unsteady convection flow.Eswara et al. 15 examined the transient laminar magnetohydrodynamic convection in a cone due to a point sink, with the free stream velocity varying continuously with time and also for the case of an impulsive change either in the strength of the point sink or in the wall temperature.They showed numerically that magnetic field increases the skin friction but decreases heat transfer and that the transient nature of the convection flow is active for short durations with suction present and greater times with injection.Chamkha 16 has analyzed the unsteady MHD free three-dimensional convection over an inclined permeable surface with heat generation/absorption.Jha 17 presented exact solutions for transientfree convection MHD Couette channel flow with impulsive motion of one of the plates discussed.More recent communications on unsteady hydromagnetic heat transfer flows include the articles by Seddeek 18 incorporating variable viscosity effects, Zakaria 19 who considered a polar fluid, and Ghosh and Pop 20 who included Hall currents 21 .Zueco presented network simulation solutions for the transient natural convection MHD flow with viscous heating effects.Bég et al. 22 studied the free convective MHD flow from a spinning sphere with impulsive motion using the Blottner difference method.Duwairi et al. 23 analyzed the unsteady MHD natural convection for the non-Boussinesq case that is, using a nonlinear density relationship for water at low temperatures.Bozkaya and Tezer-sezgin 24 have presented boundary element numerical solutions for transient magnetohydrodynamic flow in a rectangular duct with insulating walls, showing that n increase in Hartmann number causes the formation of boundary layers for both the velocity and the induced magnetic field with the velocity becoming uniform at duct centre.With increasing magnetic field, the time for reaching steady-state solution is also reduced.In many industrial applications, hydromagnetic flows also occur at very high temperatures in which thermal radiation effects become significant.The vast majority of radiation-convection flows have utilized algebraic flux approximations to simplify the general equations of radiative transfer 25 .The most popular of these simplifications remains the Rosseland diffusion approximation which has been employed by for example, Ali et al. 26 and later by Hossain et al. 27 .Radiation magnetohydrodynamic convection flows are also important in astrophysical and geophysical regimes.Raptis and Massalas 28 considered induced magnetic field effects in their study of unsteady hydromagnetic-radiative free convection.

Mathematical Model
Consider the unsteady free convective hydromagnetic flow of an incompressible, viscous incompressible, electrically-conducting fluid along an infinite hot vertical plate moving with constant velocity, adjacent to a saturated porous regime.The physical scenario is shown in Figure 1.The x -axis is oriented along the plate from the leading edge in the vertically upward direction, and the y -axis is perpendicular to this.A uniform magnetic field, B 0 , is applied parallel to the y -axis that is, transversely to the plate.Thermal radiation acts as a unidirectional flux in the y -direction.The fluid is gray and absorbing-emitting but nonscattering, and the magnetic Reynolds number is assumed to be small enough to neglect induced magnetic field effects.Hall and ionslip current effects are also neglected.The electromotive force generated by a magnetic field is a function of the speed of the fluid and the magnetic field strength.Following Shercliff 39 , we define the electrical field intensity, E, using Maxwell's equation: where t is time.The magnetic flux density, B, is defined as follows: where H is the magnetic field strength and μ e is the magnetic permeability.The generalized Ohm's law defines the total current flow as follows: where V is the velocity vector and σ is the electrical conductivity of the fluid.The electromagnetic retarding force, F magnetic , to be incorporated into the momentum conservation equation, then takes the form: Incorporating this magnetic retarding force in the momentum boundary layer equation, the appropriate conservation equations, under the Boussinesq approximation for the flow under the above assumptions, may be expressed as The appropriate boundary conditions at the wall and in the free stream are where u is the velocity along the plate, v is velocity normal to the plate, ν is the kinematic viscosity of the conducting fluid, g denotes gravity, β is the coefficient of thermal expansion, T is fluid temperature, T ∞ is free stream temperature, k 1 is thermal conductivity of the fluid, C p is the specific heat capacity, ρ is the fluid density, q r is radiative heat flux, K is the permeability of the porous regime, B 0 the component of magnetic field in the x -direction, T w is the plate temperature isothermal , and U is the velocity of the plate.Following Isachenko et al. 40 , we employ a diffusion-type radiation heat transfer approximation, namely, where σ * and k * are respectively the Stefan-Boltzmann constant and the spectral mean absorption coefficient of the saturated medium.Assuming that the temperature differences within the saturated porous regime are sufficiently small such that T 4 may be expressed as a linear function of the temperature, a power-series expansion of T 4 about T ∞ , neglecting higher order terms leads to Implementing 2.8 and 2.9 in 2.6 , we arrive at the modified energy conservation equation:

2.10
In order to render solutions to the boundary value problem described by 2.5 and 2.10 subject to the spatial and temporal conditions specified in 2.7 , we introduce a group of nondimensional transformations, defined as follows: where u is dimensionless x -direction velocity, y is dimensionless coordinate normal to the plane, t is dimensionless time, θ is dimensionless temperature, P r is Prandtl number, G r is Grashof number, K r is the radiation-conduction parameter, K 2 p is the Darcian drag force coefficient inverse permeability parameter , and M 2 is the Hartmann magnetohydrodynamic parameter.The transformed equations are thereby reduced to the following pair of coupled, second-order partial differential equations: The corresponding transformed boundary conditions become

2.14
We note that the optically-thick radiative approximation is valid for relatively low values of the parameter, K r .The electrically nonconducting version i.e., with M 2 0 of 2.12 has recently been studied by Ghosh and Bég 41 where extensive computations were provided of the influence of thermal radiation on the flow field.In the present study, we shall consider the supplementary influence of transverse magnetic field for the case where the fluid is saturated with air for which the Prandtl number is assumed to take the value 0.7.

Analytical Solution
The Laplace transform technique is now employed to generate closed-form solutions for the coupled, linear partial differential equations 2.12 and 2.13 subject to the boundary conditions 2.14 .The solutions for the transient velocity u and transient temperature θ take the following expressions:

3.1
The spatial gradients of these functions provide expressions for the dimensionless shear stress i.e., related to the skin friction and temperature gradient i.e., related to the Nusselt number at the plate surface:

Network Numerical Solution
Numerical solutions to the two-point boundary value problem have also been obtained with the Network Simulation Method NSM .This powerful and robust computational method has been employed extensively by the authors in a wide spectrum of both linear and nonlinear steady and unsteady magnetohydrodynamic and thermal convection flows.A number of networks are connected in series to make up the whole medium.After experimenting with a few sets of mesh sizes, a region of integration of 200 cells has been selected.Boundary conditions are subsequently added by means of special electrical devices current or voltage control-sources that is, resistors, capacitors, and so forth.Once the complete network model is designed, the Pspice code is employed for the numerical simulations.This code is designated the "electric circuits simulator".Using the Fourier Law, the spatial discretization of 2.12 and 2.13 gives Δy/2 .The electrical analogy is applied to 4.1 together with Kirchhoff's law for the currents.To implement the boundary conditions at y 0 and y 1, constant voltage sources are employed for both conditions.The principal advantage of the NSM approach is that it avoids the necessity in traditional numerical difference schemes of manipulation of difference equations and the specified constraints concerning the convergence of numerical solutions.For example, the time-step used in transient problems, which is required for convergence is not a prerequisite as Pspice achieves this via sophisticated numerical algorithms largely analagous to those intrinsic to the standard difference numerical solvers, as described by Nagel 47 .Design of the model does require a comprehensive appreciation of electrical circuit theory.Momentum balance "currents" are defined systematically for each of the discretized equations and errors can be quantified in terms of the quantity of control volumes.The network model is shown in Figure 2 for the momentum equation 2.12 and Figure 3 for the energy equation 2.13 .

Results and Discussion
We have obtained a comprehensive range of solutions to the transformed conservation equations.To test the validity of our numerical NSM computations, we have compared the velocity and shear stress distributions in Tables 1 and 2 with the Laplace transform solutions.Very good correlation is apparent.In all computations the key thermophysical parameters have been prescribed as follows, unless otherwise stated: G r 10, K r 2.0, K p 1.0, M 2 5.0, P r 0.7, and t 0.25, corresponding to free convection of air in a highly porous regime with strong magnetic field and high thermal radiation flux at intermediate time.
Both Tables 1 and 2 correspond to distributions computed a short time after the initiation of motion that is, at t 0.2.In Table 1, we observe that an increase in Hartmann number M 2 from 5.0 through 8 to 10 strong magnetic flux density causes a significant decrease in the flow velocity, u with distance normal to the plate surface into the boundary layer y .This trend is consistent with many classical studies on magneto-convection showing that the hydromagnetic body force retards the flow that is, decelerates the fluid causing a thinning in the boundary layer thickness.Very high Hartmann numbers i.e., M 2 1 are usually associated with the formation of a Hartmann boundary layer 39 .In Table 2, the shear stress, ∂u/∂y| y 0 , is found to be decreased significantly with an increase in Hartmann number M 2 from 5.0 through 8.0, 10.0 to 12.0, for all values of the radiation-conduction parameter, K r .In all cases, the shear stresses are negative since the high values of Hartmann number, M 2 , retard the flow in the boundary layer to such an extent that reversal of the flow is caused.This result is significant in the design of, for example, MHD generators since a critical magnetic flux density may be applied i.e., Hartmann number to reverse the flow dynamics during operation.The change in shear stresses with an increase in radiation-conduction parameter indicates that an increase in thermal radiation has a positive effect on the flow that is, reduces the degree of flow reversal.For example, for M 2  5.0, shear stress, ∂u/∂y| y 0 increases for the NSM solutions from −2.04602 for K r 0.5 thermal conduction exceeds thermal radiation to −2.0252 for K r 1.0 for which the thermal radiation and thermal conduction mode contributions are approximately the same , to −1.9960 for K r 2.0 thermal radiation exceeds thermal conduction .Therefore, for very high-strength magnetic field operating conditions, thermal radiation mitigates to some extent flow reversal effects.In Figure 4, the dimensionless shear stress profiles in time for different Hartmann numbers M 2 are illustrated.A strong decrease is observed in shear stress from t 0 to t 0.05, after which profiles, although they continue to decrease with increasing M 2 values, tend for M 2 0, 5, 10, and 15, to the steady state.For these values of Hartmann number, profiles are always positive indicating that flow reversal does not occur.Comparing with the trends in Table 2, we can note that while the values of K r , K p , and P r are identical for the third column, the Grashof number G r is lower in Table 2 at 2.0, compared with the value of 10.0 in Figure 4.As such the flow is more strongly assisted by buoyancy forces in Figure 4 which prevents the reversal of flow for all the profiles, with the exception of M 2 20 very strong transverse magnetic field which becomes negative for t > 0.4, Thermal buoyancy force, G r θ , is directly proportional to the Grashof free convection parameter G r and therefore would appear to assist the flow, whereas magnetic field inhibits flow acceleration in the regime.In Figure 5, temperature gradient profiles in time for the influence of radiation-conduction parameter K r are presented.In all cases, profiles are a maximum initially at the isothermal plate and decay quickly from the wall with time.An increase in K r from 0.1 through 0.5, 1, 1.5 to 2.0 is seen to markedly reduce heat transfer gradient especially at shorter times 0 < t < 0.2 ; with further elapse of time all profiles converge that is, radiation effects are negligible for large times.Increasing K r implies a greater augmentation of heat transfer by thermal radiation which will serve to increase fluid temperatures in the regime; the spatial heat transfer rate ∂θ/∂y| y 0 that is, temperature gradient at the wall will therefore be reduced as greater thermal energy heat will be imparted to the fluid-saturated regime raising temperatures within the porous regime.In Figure 6, the influence of the Darcian drag force parameter, K p , on the time evolution of shear stress profiles, is depicted.As K p increases from 0.1 through 0.5, 1, 2 to 5, a very large escalation in Darcian drag force is caused, as expressed in 2.12 in the linear term, −K 2 p , which decelerates the flow and reduces the shear stress at the plate.Steady-state values are achieved faster with lower Darcian drag K p 0.1 than with higher Darcian drag K p 5 .The effect of Grashof number on velocity gradient i.e., shear stress through time is presented in Figure 7. Increasing G r for the case of very low Darcian drag i.e., highly permeable medium, K p 0.1 strongly increases shear stress values at the wall that is, accelerates the flow over time.We note that values become negative for very low G r values since the magnetic impedance force M 2 5.0 will dominate and have a greater inhibiting influence with low buoyancy that is, flow reversal accompanies lower thermal buoyancy forces for higher permeability regimes.In Figure 8, the effect of the Prandtl number P r on the temporal shear stress distribution is shown again for G r 10.0, K r 2.0, M 2 5.0 but with K p 1.0.Increasing P r strongly boosts the flow and increases shear stress profile values which remain positive for small times; however, with increasing elapse of time shear stress, values become negative indicating backflow occurs at the plate.For lower P r values 0.7, 0.1 , negative values are attained more quickly that is, backflow takes place quicker.In Figure 9, an increase in P r is observed to the enhance temperature gradient.Prandtl number controls the relative thickness of the momentum and thermal boundary layers.When P r is of low value, heat diffusion exceeds momentum diffusion.For P r < 1, the thickness of the thermal boundary layer therefore exceeds the thickness of the velocity boundary layer that is, temperatures will be greater.In Figure 10, temperatures are seen to decrease considerably with an increase in P r values for a fixed time, t 0.25 as we progress into the boundary layer regime; profiles also decay much more sharply for higher P r values since momentum diffusion exceeds energy heat diffusion for P r > 1.For the case of P r 1, the boundary layer thicknesses will be approximately of the same order of magnitude.For P r 0.1, the profile is approximately linear for a substantial distance from the plate.Spatial velocity u distributions, for two time values are illustrated in Figure 11, for the effect of Hartmann hydromagnetic parameter M 2 .This parameter represents the ratio of the hydromagnetic retarding force to the viscous hydrodynamic force in the boundary layer.The classical velocity overshoot is identified 1, 2, 4, 39 near the moving plate surface for lower values of M 2 that is, 1.0 and 5.0; with M 2 10.0 this overshoot is clearly suppressed owing to stronger resistance to the flow.We note that for t 1.0, the profiles are always greater in value than for t 0.25 that is, the flow is accelerated considerably with time, although velocities are strongly reduced with an increase in Hartmann number.All profiles decrease towards zero in the free stream, although this state is attained much faster for higher magnetic field values M 2 10 and for shorter times.In Figure 12, the combined effects of time t and radiationconduction parameter K r on spatial distribution of temperature θ through the boundary layer is shown.An increase in K r serves to supplement fluid thermal conductivity with radiation contribution and significantly heats the fluid-saturated regime that is, increases temperature values.Similar results were reported by Ali et al. 26 ,Hossain et al. 28 , and very recently by Ghosh and Bég 41 .A large difference is observed between the profiles computed at t 1.0 and t 0.25, indicating that thermal radiation effects are amplified at greater times, compared with smaller times where the flow is still developing.After greater times a greater quantity of thermal energy will be absorbed into the fluid regime via the imposed flux causing enhanced heating of the fluid.For example, for t 1.0, at y 2, for K r 5 maximum thermal radiation effect , θ reaches a value of approximately 0.65, whereas the corresponding value for t 0.25 is much lower at 0.35.Finally in Figure 13 we have plotted the spatial variation of velocity for the combined effects of radiation-conduction parameter K r and time t .Again a velocity overshoot is observed in the close vicinity of the plate; however, this overshoot is distinctly greater for the highest value of K r 2.0 and greater time values t 1.0 ; all profiles descend gradually to zero far from the wall.Thermal radiation therefore augments the flow that is, accelerates the flow in the porous regime.Velocities are minimized when thermal conduction swamps thermal radiation that is, for K r 0.1.

Conclusions
Closed form and numerical NSM solutions have been presented for the transient hydromagnetic natural convection boundary layer flow past a moving vertical plane adjacent to a Darcian porous regime with thermal radiation flux present.It has been shown that thermal radiation strongly increases fluid temperatures and accelerates the flow; conversely magnetic field as simulated via the Hartmann number serves to impede the flow and reduce velocity gradient shear stress values.The effects of both parameters are enhanced with a greater elapse of time.Darcian drag is seen to decelerate the flow, whereas increasing free convection serves to accelerate the flow owing to the assistance of thermal buoyancy forces in the regime.A rise in Prandtl number however decreases temperatures in the regime, but accelerates the flow that is, increases velocity gradient values.A velocity overshoot is observed with magnetic field effects but vanishes for very high values of the Hartmann number.

Figure 1 :
Figure 1: Physical regime and coordinate system.
Zueco 42  recently studied the periodic temperature variation effect on thermal convection in a horizontal channel.Other very recent studies employing NSM include the works of Bég et al. 43 considered unsteady rotating Couette flow in a porous medium channel, 44 analyzed the magnetohydrdoynamic rotating flow in a Darcian channel with dissipation effects and Hall/ionslip currents, and 45 examined the effects of thermal stratification and non-Darcian drag on natural convection boundary layers in a porous regime.In the NSM technique, a second-order central-difference scheme is utilized to discretize the momentum and energy conservation equations and the resulting system of finite difference equations are solved employing the Pspice program 46 .A network model is subsequently designed, with component equations which are formally equivalent to the discretized ones.The electrical analogy relates the electrical current J with the velocity flux ∂u/∂y and temperature flux ∂θ/∂y , while the electrical potential Φ is equivalent to the velocity, u and temperature θ.

Figure 2 :
Figure 2: Electronic network model for the momentum equation 2.12 .

Figure 11 :
Figure 11: Spatial velocity distribution u with K p 1, G r 10, K r 2.0, and P r 0.7 for the effect of Hartmann hydromagnetic parameter at t 0.25 and t 1.

7 Figure 12 :
Figure12: Spatial temperature distribution θ with K p 1, G r 10, M 2 10, and P r 0.7 for the effect of radiation-conduction parameter K r at t 0.25 and t 1.

Figure 13 :
Figure13: Spatial velocity distribution u with G r 10.0, M 2 5, K r 2.0, and P r 0.7 for the effect of radiation-conduction parameter K r .

Table 1 :
Electronic network model for the energy equation 2.13 .Dimensionless spatial velocity distribution u for G r 2.0K r 1.0, K p 1.0, P r 0.7, and t 0.2 for the effect of Hartmann number M 2 .

Table 2 :
Dimensionless shear stress ∂u/∂y| y 0 for G r 2.0, K p 1.0, P r 0.7, and t 0.2 for the combined effects of radiation-conduction parameters K r and Hartmann number M 2 .