A combined method for simulating MHD convection in square cavities through localized heating by method of line and penalty-artificial compressibility

A new combined technique is developed for studying free convection of magnetohydrodynamic unsteady incompressible flow in a square cavity that partially heated from lower wall and regular cooling from lateral walls. Convection has been studied for several lengths of heat source, Hartman and the Rayleigh numbers. The local and average Nusselt numbers are calculated and presented graphically along the heat source portion. The proposed technique is a combination of the method of lines (MOL) and a developed penalty-artificial compressibility technique (PACT). Unsteady penalty compressibility governing equations are used instead of solving steady-state balances. The penalty and artificial compressibility factors are estimated using a deduced stability restriction. The validation of the MOL-PACT is verified by comparing the current data with other numerical and experimental ones existing in the literature. The proposed technique provides a new choice through the investigation of MHD mass and heat transfer.


Introduction
Free convection, induced by internal heat generation, in square cavities has received significant attention because of its many applications in geophysics and energy related to engineering problems. Indeed, natural convection in closed enclosures with various heated boundaries is a model of various engineering applications like the safety and operation of a nuclear reactor and convective heat transfer related to electronic cooling devices. In some practical applications such as metal casting, the crystal growth in a fluid, the extractions of geothermal energy and the fusion reactors, the free convection has occurred within the magnetic field effect.
Analysis of the literature reveals comparatively huge attention in the convective phenomenon study in a square and rectangular cavity heated from the lower wall and subjected to dissimilar isothermal conditions. As example, and not as a limitation, Sharma et al. [1] investigated the turbulent free convection in a square cavity containing a localized heater at a lower wall and identical cooling at vertical sidewalls. Their research models the situation of the unintended heat generation resulted from a fire in an isolated body of an electronic components cabin or a nuclear reactor. In [2], the nat-ural convection in a rectangle cavity which was partly heated from the lower wall and equivalently cooled at the vertical sidewalls has been numerically investigated. Calcagni et al. [3] experimentally and numerically investigated natural convection and heat transfer in a square cavity that happened by a heater positioned on the bottom wall and cooled from the sidewalls. In this research, they examined how the transferring of heat grows in the cavity due to the increase of the heat source length. A systematic numerical study of 2D steady laminar free convection in a rectangle cavity due to localized heaters of various walls thermal conditions is considered by Deng et al. [4]. Yigit et al. [5] numerically studied the steady laminar free convection of power-law fluids in a square cavity containing a partial heater at the bottom wall and identical cooling at sidewalls by ANSYS-FLUENT solver. Recently, many researchers have investigated many problems related to MHD free convection. In [6], Colak et al. have examined magnetohydrodynamic mixed convection in a chamfered lid-driven enclosure containing a fractional heater. Sajjadi et al. [7] have done a 3D mesoscopic simulation of magnetohydrodynamic free convection in a cubic enclosure using the Lattice-Boltzmann method and with two multi relaxation time models. In [8], Usman et al. have investigated the magnetohydrodynamic heat and mass transfer flow inside a square enclosure heated from the top wall and equipped with cold and heated square obstacles. Sivaraj and Sheremet [9] investigated the free convection in a tending porous enclosure containing heated compact obstacle located at the centre in presence of an applied magnetic field of diverse orientations. Zhang et al. [10] studied the MHD free convection in square and cube cavities in presence of the effects of thermal radiation using the spectral collocation-artificial compressibility method (ACM). In [11], Shahid and Tunç have obtained exact expressions of temperature and velocity of an elastoviscous magnetohydrodynamic fluid slipping past an unbounded upright channel with a vacillating temperature.
Method of lines [12][13][14][15][16] as modified (G /G)-expansion method [17,18] is one of the efficient and accurate numerical methods that have been handled successfully to solve many partial differential equations (PDEs) with high accuracy. The idea of the MOL is summarized in reducing the given partial differential equations to a system of ordinary differential equations (ODEs) in the time by discretization of space variables and spatial derivatives using finite difference schemes. Then the resultant ODEs system can be solved using a suitable ODE solver like the 4th order Runge-Kutta method (RK4). As a global approach, the MOL is used in this work for solving PDEs of Navier-Stokes that describes the flow of incompressible fluids. In order to utilize this explicit method in solving such PDEs, the ellipticity in Navier-Stokes equations must be removed. The leading method to achieve this is using the penalty method [19], which was received significant attraction in the last years, especially in the finite element and free convection community [20][21][22][23]. The idea of this technique is based on eliminating the pressure variable (p) in the momentum conservation equations with the aid of a mass conservation equation (continuity equation) that will be used as a restriction. In the penalty method, the errors in the divergence of velocity are transferred to the boundaries [24]. An alternative related technique is the artificial compressibility technique (ACT) [25][26][27][28][29], in which the deviations from incompressibility determine ∂p/∂t rather than p itself. Also, the ACT diffuses errors in velocity divergence. The penalty method and ACT converge to the Navier-Stokes equations in the limit as their parameters approach to infinity [30].
Therefore, we try to utilize a combination between MOL, penalty method and ACT for the purpose of spatial discretizing of PDEs of unsteady penalty-artificial compressibility and present the MOL-PACT for solving PDEs governing natural convection of MHD laminar incompressible fluid in a square cavity with a partial heater at the bottom wall and cooled from sidewalls.
Optimal values of penalty and artificial compressibility parameters are calculated using a deduced stability condition. The validation of the developed technique is tested by comparing the obtained results with other experimental and numerical results existing in the literature as well. As we will see in the model validation section, an excellent agreement between the obtained results and the results that exciting in the literature [1,2], when neglecting the magnetic field effect, is achieved although the simplicity of the method. This can be considered one of the advantages and benefits of the proposed method.

Mathematical and physical model
Suppose the motion of viscous fluid in a square encloser of length L · as shown in Figure 1. The heat source of length l is centrally positioned on the lower wall and varied from 0.2 to 0.8 of L · . The l L · is termed as ε (dimensionless heated length piece). The upper wall and the unheated piece of the lower wall are adiabatic; however, the vertical sidewalls are maintained at a fixed temperature T c and the heated piece of the lower wall is at a temperature T h . The change in the fluid density with temperature is combined using the Boussinesq approximation. The encloser is subjected to a horizontal uniform magnetic field with a magnitude B 0 , [31] while the brought magnetic field is neglected being small as compared to B 0 . Enclosure walls are supposed to be electrically deposed with no Hall impacts. Under these hypotheses, the transient form of continuity, energy and momentum equations can be obtained from [21] by adding the term which represents the magnetic field effect or from [32] by removing the term which represents the medium permeability effect, wherever u and v are the velocities in x and y directions sequentially, p is the pressure, ρ is the density, β is the coefficient of thermic expansion, α is the thermic diffusivity, ν is the viscosity, σ is the conductivity and g is the gravity acceleration. The boundary conditions (BCs) for our case are: • vertical sidewalls walls: Substituting the nondimensional variables: into Equations (1)-(4) yields dimensionless equations ∂ϑ ∂τ +Ù ∂ϑ ∂χ And hence the BCs are reduced to: • vertical sidewalls walls: • upper wall: HereÙ, V, P and ϑ are the nondimensional velocities components, pressure and temperature, respectively. The parameters P r , R a and H a represent the Prandtl, Rayleigh and Hartman numbers, respectively. The rate of the heat transfer through the heated portion of the lower wall can be quantified using the averaged Nusselt symbol Nu aνg and local Nusselt symbol (Nu(χ)) that given by, respectively, The motion of the fluid can be presented by the stream function ψ that can be perceived from the veloc-itiesÙ and V as: Equation (11) leads to the following Poisson's equation, for obtaining the stream function ψ.
Because of a nonslip condition, ψ = 0 is considered at all boundaries of the cavity.

Penalty-artificial compressibility technique (PACT)
Based on the penalty method and the artificial compressibility technique, instead of continuity equation (5), we present the following penalty-artificial compressibility equation where the constant c represents the artificial compressibility parameter and b represents the penalty parameter. When c = 0, Equation (13) 8) and (13), firstly we construct the MOL finite difference scheme and then we apply a suitable stability technique as the conventional linearized Fourier stability analysis.

MOL and stability analysis
As indicated by the MOL, the coordinate χ and Y in the system of Equations (6)- (8) and (13) is discretized with N × N regularly 2D spaced grid points where . . , N and χ = Y = 1/N. Central finite difference scheme with a second-order accuracy used to approximate the first and second derivatives corresponding to the spatial variables χ and Y at each grid point (χ i , Y j ) for i, j = 2, 3, . . . , N − 1.Then the 4th Runge-Kutta technique is utilized as a time integrator to solve the ODEs system upon the pre-described initial conditions with an appropriate step time τ . The central first-and second-ordered finite difference operators regarding χ and Y are defined as: In order to drive a stability condition needed to estimate the constants b and c, we write the finite difference scheme of Equations (6)-(8) and (13) in a forward Euler explicit representation with respect to time τ , as where k is the time level.
Ignoring the convection terms in the present finite difference scheme and applying a linearized von Neumann stability analysis [33] yield the following stability condition [34,35]: , and defining a mixture ratio M = c 2 τ 2(b+P r ) , the stability restriction (18) can be written as The values M = 0 and M = ∞ correspond to the pure penalty technique and ACT, respectively, while values near M = 1 correspond to a roughly equal mixture of the two techniques. In order to estimate a value of the parameter b form inequality (19), we introduce a safety factor f where 0 < f < 1. Thus, we can obtain b and c 2 that are safely consistent with stability condition as From the relations (20) and (21), that satisfying the stability condition, one can estimate the values of b and c 2 for a given grid by setting arbitrary values of Prandtl number P r , mixture ratio M, safety factor f and selecting a small value of τ .
All results that will be discussed in Section 5 are computed using a grid of size 40 × 40 points for P r = 0.71 (air), M = 1, f = 0.8 and τ = 10 −5 . These arbitrary values give that b = 5.54 and c 2 = 1.25 × 10 6 . The optimal value of the mixture ratio M is selected upon analysing maximum absolute deviations (relative error) in the continuity equation for several values of M ranging from M = 10 −4 to M = 10 4 . This analysis is not shown here. A detailed analysis of choosing the optimal value of the mixture ratio M and the safety factor f can be found in [34,35].
The grid of size 40 × 40 is decided upon the next grid independence study. RK4 technique is used as an integrator of the resulting ODEs system, instead of Euler technique, for the purpose of obtaining more accurate results. The stream function ψ is obtained by applying the successive over-relaxation technique, with suitable accuracy, to the system of linear algebraic equations resulted from the following equation The steady-state solution is checked by employing the following condition where stands forÙ, V, ϑ, P or ψ.

Study of grid independency
Calculations have been gotten on three meshes to test the grid independency of the numerical scheme: G1 : 20 × 20, G2 : 40 × 40 and G3 : 80 × 80. The grid independence study is done for the average Nusselt number along with the heated piece of the lower wall when R a = 10 6 , H a = 0, ε = 0.2 and ε = 0.8. The results of the grid study for the three grids are presented in Table 1.
As shown in Table 1, the variance in Nu avg. between the grids G1 and G2 is 3.59% for ε = 0.2 and 3.25% for ε = 0.8, while between G2 and G3 is only 0.73% for ε = 0.2 and 0.39% for ε = 0.8. So, a grid pattern of G2 is considered in all study computations.

Validation of the model
Through the objective of consequences validation, the results got from the present work are discussed with the experimental and numerical ones of Sharma  [3] with identical boundary conditions and R a values in absence of the magnetic field effect, i.e. H a = 0. The average Nusselt number at the heated portion of the lower wall of dimensionless length ε = 0.8 are discussed in Table 2, with the obtained results. It can be noted that the agreement is good. For further results validation, some of the obtained results are validated against ones obtained in [2,3] when H a = 0 as shown in Figure 2. The upper row of Figure 2 contains plots of temperature contours when considering ε = 0.8 at R a = 10 5 , while the lower row contains plots of the Nusselt number variations along the heated piece of the lower wall with the Rayleigh symbol for all values of ε. From Figure 2, one can state that, a good agreement is detected between present consequences and numerical experimental results existing in [2,3]. This agreement is enough to ensure the validation and accuracy of the present computational scheme.

Results and discussions
Here, we discuss the results of the MOL-PACT numerical solution for simulating 2D steady laminar flow under the effects of MHD inside a square enclosure with a partially heating from the lower wall and symmetrical cooling from vertical sidewalls with an isolated upper side. The results are investigated by obtaining isotherms, streamlines, local and average Nusselt numbers for several  values of R a (10 3 to 10 6 ), H a (0 to 100) and ε(0.2 to 0.8) while P r is remained fixed at 0.71. The simulation is done to study the impact of Hartman number, Rayleigh number, and dimensionless heated length part on fluid stream function, temperature distribution (isotherms) and rate of heat transfer.  Moreover, temperature distribution will be stratified and the stratification degree grows with the growing of R a . The larger R a means further fluid circulation, and so, further heat is added to the fluid to upsurge the convection of the fluid. The Hartmann number presents a measure of the magnetic field effect on thermal free convection. The amount of the heat transfer dramatically reduces with the growth of H a due to the fact, that the magnetic field tends to resist the fluid motion inside the cavity as shown in Figures 3 and 5.
From Figures 4 and 6 when H a = 0, one can detect that the growing of the Rayleigh number will grow the strength of the circulation inside the cavity and make the cores of the vortex's cells move upward as concluded in [2]. On other hand as H a increases, the main recirculation cell moves downward because of the magnetic field strength effect. For a fixed H a , the streamlines remain almost similar but the boundary layer near the midpoint of the bottom wall becomes slightly compressed.
To show the impacts of R a , H a and ε on the heat transfer rate inside the cavity, it is important to visualize the local and average Nusselt number. Figure 7 represents Nu along heated part of the lower wall for various values of H a and ε at R a = 10 6 . For a fixed ε, increasing H a resists the fluid convection. Moreover, growing ε for a fixed value of H a produces a raise in heat transfer. These outcomes can be obviously illustrated under the temperature distribution demonstrated in Figures 3 and 5. Figure 8 illustrates the variations of Nu avg. along heated part of the lower wall versus R a for all values of ε, H a = 0 and H a = 100. This figure summarizes the effects of R a , H a and ε on the heat transfer rate inside the fluid. It can be shown that the rate of heat transfer enhances by increasing Rayleigh number and dimensionless heated length, while it reduces by decreasing the magnetic field parameter H a . When H a = 0, the shown results in Figure 8, are consistent with the results obtained in [2,3]. On other hand, when H a = 100, the values of Nu avg. remain approximately constant till R a = 10 5 then begin to increase when R a = 10 6 .

Conclusion
In this work, the mixed technique MOL-PACT is developed for modelling MHD laminar flows and heat transfer problems. Although the developed technique is proposed to generate accurate transient solutions, it is also applicable to calculate steady-state solutions. The developed technique is combining hyperbolic and parabolic nature to get improved convergence behaviour than either one alone may offer. The penalty and artificial compressibility parameters that consistent with the deduced stability restriction are estimated.
As an example, to verify the efficiency of the MOL-PACT, the technique is used to model the MHD free convection in a square cavity with a localized heating. The simulation is done to investigate the influence of magnetic field parameter, Rayleigh number and heated source length on the flow, isotherms and rate of the heat transfer. The main conclusions of the simulation are as follows. The flow and temperature distribution are symmetrical around mid-length of the cavity because of the symmetry nature of considered conditions. By increasing the Rayleigh number, the flow recirculation growths, and the cores of the vortex's cells move upward. The heat transfer rate enhances by increasing Rayleigh number and heated source length, while it reduces by decreasing the magnetic field effect. Some of the obtained results are consistent and in a good agreement with ones found in the literature.
Finally, the novelty of the paper can be summarized in that the proposed method provides a new efficient choice for the investigating of MHD laminar flow and heat transfer related problem. The generalizations of MOL-PACT to be applicable for handling compressible flows of low Mach number and nanofluids are of possible attention and are being pursued.