Modeling the Natural Convection Flow in a Square Porous Enclosure Filled with a Micropolar Nanoﬂuid under Magnetohydrodynamic Conditions

: The laminar, natural convective ﬂow of a micropolar nanoﬂuid in the presence of a magnetic ﬁeld in a square porous enclosure was studied. The micropolar nanoﬂuid is considered to be an electrically conductive ﬂuid. The governing equations of the ﬂow problem are the conservation of mass, energy, and linear momentum, as well as the angular momentum and the induction equations. In the proposed model, the Darcy–Brinkman momentum equations with buoyancy and advective inertia are used. Experimentally obtained forms of the dynamic viscosity, the thermal conductivity, and the electric conductivity are employed. A meshless point collocation method has been applied to numerically solve the ﬂow and transport equations in their vorticity-stream function formulation. The e ﬀ ects of characteristic dimensionless parameters, such as the Rayleigh and Hartmann numbers, for a range of porosity and solid volume fraction of Al 2 O 3 particles in a water-based micropolar nanoﬂuid on the ﬂow and heat transfer in the cavity are investigated. The results indicate that the intensity of the magnetic ﬁeld signiﬁcantly a ﬀ ects both the ﬂow and the temperature distributions. Moreover, the addition of nanoparticles deteriorates the heat-transfer e ﬃ ciency under speciﬁc conditions. increases, and this slope declines as the Rayleigh number increases. As the Hartmann number increases, the Nusselt number recedes further. K.M., and V.N.B.; funding acquisition, N.P.K., K.M., and V.N.B.; methodology, N.P.K., G.C.B., E.D.S., and V.C.L.; project administration, K.M. and V.N.B.; software, N.P.K. and G.C.B.; supervision, K.M. and V.N.B.; validation, N.P.K., G.C.B., E.D.S., and V.C.L.; visualization, N.P.K.; writing—original draft, N.P.K. and G.C.B.; writing—review and editing, E.D.S., V.C.L., K.M.,


Introduction
Magnetohydrodynamics (MHD) focuses on electrically conducting fluids, which move under the influence of a magnetic field. The induced electrical current interacts with the magnetic field to produce a body force acting on the fluid. The interaction between the flow field (fluid velocity) and the externally applied magnetic field gives rise to a magnetic (induced) field inside the fluid. The MHD-oriented research field has become an important branch of fluid dynamics. The study of MHD flow applies to many industrial and engineering systems, such as nuclear reactors, heat exchangers [1], home ventilation systems, cooling of electronic equipment, solar energy collectors, nuclear reactors, chemical processing equipment, geothermal reservoirs, magnetic behavior of plasmas in fusion reactors, petroleum industries [2], boundary layer control in aerodynamics, crystal growth, and ship propulsion [3], to name a few [4,5].
Furthermore, MHD has been studied by many researchers, and considerable effort has been directed toward MHD flow combined with heat transfer. In fact, many authors have studied the effects of magnetic field on mixed, natural, and forced convective heat-transfer problems. Chandra and Gosh [6] studied the effect of a magnetic field on an electrically conducted viscoelastic fluid, reporting on the flow decrease through the increase in the magnetic field strength. Jha [7] studied the effects of natural convection and applied magnetic field on the Couette flow of an electrically conducting fluid between two parallel plates. Duwairi and Duwairi [8] reported on the thermal radiation heat transfer effects on the MHD-Rayleigh flow. Di Piazza and Ciofalo [9] studied the magneto-convection in a lead-lithium mixture, showing that a magnetic field, directed perpendicular to the temperature gradient, produced a complex three-dimensional flow pattern. Kenjeress and Hanjalic [10] studied the Rayleigh-Benard convection, subjected to various magnetic fields, over a range of Rayleigh and Hartmann numbers, and reported on the impact of the orientation and distribution of externally applied magnetic fields on the heat transfer and the flow regime (convective structures) in thermal convection of electrically conductive fluids. They studied the flow laminarization over specific range of Hartmann numbers. Through laminarization, it is possible to control thermal convection and associated heat transfer.
Over the last few decades, there has been a growing interest in the heat-transfer properties of nanofluids. Nanofluids are engineered colloids, i.e., a mixture of a base fluid (usually water, oil, or glycerin) and nanosized solid particles [11][12][13][14][15]. Nanofluids exhibit novel properties that make them potentially useful in many applications in heat transfer, including fuel cells, pharmaceutical processes, engine cooling/vehicle heating, and heat exchangers. They possess increased thermal conductivity and heat-transfer coefficients, compared to the base fluid. However, a critical parameter for their suitability in many applications is their rheological behavior.
In numerous computational studies, nanofluids are modeled as single-phase fluids, having altered physical properties compared to the ones of the base fluid [16][17][18]. Mahmoudi and Abu-Nada [19] studied the combined effect of the magnetic field and the nanofluids' behavior on natural convection within a square cavity. They reported that the presence of the nanoparticles has an adverse effect on the heat transfer at a high-volume fraction of the nanoparticles. Makinde [20] studied the effects of viscous dissipation. Their results show that the increase of the volume fraction of nanoparticles enhances the heat-transfer rate. However, intensification of viscous dissipation has a negative effect in the heat-transfer rate. Jani [21] studied the free convection heat transfer in rectangular enclosure. The obtained results showed that the rate of heat transfer increased only at low Rayleigh numbers. This contradicting behavior has not been studied with adequate consideration and is one of the objectives of the present work.
Additionally, intensive research is conducted in order to comprehend the underlying physics of nanofluids, particularly the phenomena related to their increased thermal conductivity, compared to the thermal conductivity of the base fluid. Seven main mechanisms have been proposed toward a more comprehensive understanding of their physical properties, and, among them, nanoparticle microrotation seems to play a crucial role [22][23][24]. A mathematical model that takes into consideration the microrotation has been used in the study of Bourantas and Loukopoulos [25]. Therein, authors applied the micropolar flow theory to model the nanofluid. They used an experimentally obtained formula for the thermal conductivity and the dynamic viscosity of the nanofluids, and compared the numerical results with experimental data, with the accuracy being excellent.
Considerable effort has been directed toward heat-transfer studies for natural convection in porous media [26][27][28][29][30]. Oztop [31] studied the heat transfer with natural convection of an inclined partially cooled rectangular enclosure filled with a porous medium. Natural convection heat transfer in a cavity saturated with porous media in the presence of magnetic fields is a new branch of thermo-fluid mechanics. Shehadeh and Duwairi [32] studied MHD natural convective heat transfer in a rectangular enclosure filled with a porous medium, considering Joule and viscous heating effects.
In the present study, the micropolar flow theory is utilized to study the natural convection of a micropolar nanofluid in a square cavity filled with a porous medium and subjected to a magnetic field. Given its practical interest in the engineering and medical fields, including blood flow and purification [24], exportation of geothermal energy, purifications of metals, and solar heating systems [33], the topic needs to be further explored. The work on natural convection in a square cavity filled by a micropolar nanofluid [25] is extended here to account for the applied magnetic field in the presence of a porous medium. We developed a robust, efficient, and accurate meshless strong form collocation solver to numerically solve the governing flow equations. The latter are a set of nonlinear, fully coupled partial differential equations (PDEs), which express the conservation of mass, energy, linear momentum, and angular momentum, in addition to the magnetic induction. In the present study, the velocity-vorticity formulation is used, in order to avoid the problem of determining the pressure boundary conditions explicitly. To this end, a sophisticated five-step algorithm is implemented numerically. The criterion for terminating the algorithm is the normalized root-mean-square error of all field variables becoming less that a predefined value. An advanced Discretization Corrected Particle Strength Exchange (DC PSE) interpolation method is reformulated to calculate the spatial derivatives up to the second order in an Eulerian framework.
The physical aim of this study is to investigate the effect of the magnetic field magnitude and the intensity of the flow field, as well as the volume fraction of the nanoparticles and the permeability of the medium, on the convective heat transfer of an Al 2 O 3 /water nanofluid. The streamlines, temperature, microrotation, and vorticity contours are presented for characteristic ranges of solid volume fractions (0 ≤ ϕ ≤ 0.03), Hartmann numbers (0 ≤ Ha ≤ 1000), Rayleigh numbers (10 6 ≤ Ra ≤ 5 × 10 8 ), and porosity of the medium (0.5 ≤ ε ≤ 0.7). In the proposed model, the Darcy-Brinkman momentum equations with buoyancy and advective inertia terms are used. The microrotation vector is affected by the presence of the particles and the pore walls. For this reason, a single-fluid equation in a porous setting is implemented, incorporating the effective microrotational properties of the nanofluid, using volume-averaging in representative elementary volumes (REVs), avoiding the necessity of considering the effect of the detailed pore structure on the form of the equation itself. Similar manipulation has been applied for the presence of the magnetic field. In this way, a magnetic field term is not required in the transport equation for the microrotation vector.
The rest of the paper is as follows. In Section 2, the governing equations of the proposed micropolar nanofluid model are presented, while the validation of the numerical method used for the solution of the flow equations and of the proposed micropolar nanofluid model, as well, are depicted in Section 3. In Section 4, the heat-transfer performance of a micropolar nanofluid enclosed in a rectangular cavity is studied. For all simulations, pure water is considered as the base fluid, with Pr = 6.2. Finally, in Section 5, the conclusions complete the paper.

Problem Formulation
In the present study, we consider the viscous, laminar, incompressible, natural convective flow in a square porous enclosure of length (L) filled with a micropolar Al 2 O 3 /water nanofluid in the presence of a magnetic field. The geometry and the configuration of the system are schematically shown in Figure 1.
The top and bottom walls are insulated, and the fluid is isothermally heated and cooled at the left and right vertical walls at uniform temperatures, T h and T c , respectively. In the case of natural convection, the buoyancy, which is a result of local fluid density variations due to temperature differences, is the dominant driving force, and the flow is permeated by an applied uniform external magnetic field, B, in the −y direction of strength, B 0. The fluid is in thermal equilibrium with the porous medium; the physical properties of the fluid are assumed constant, except of the density in the buoyancy force, which is calculated by using the Boussinesq approximation. We consider a granular porous medium of glass spheres. The thermophysical properties for the nanofluid and the glass spheres are listed in Table 1. The fluid is hydrodynamically, thermally, and electrically isotropic. No-slip and impermeable boundary conditions are imposed on all boundaries. Joule heating, viscous dissipation, pressure work, and Hall effects are neglected, as they are assumed negligible. Table 1. Thermophysical properties of water, alumina nanoparticles, and glass spheres, namely, the density ( ), the heat capacity ( ), the thermal conductivity ( ), the thermal expansion coefficient ( ), and the electric conductivity ( ) [34][35][36]. From Ohm's law, the current density, J, relates to the electric field, E, as follows: where is the electrical conductivity, and = ( ( , ), ( , )) is the superficial velocity vector, with u and v being the velocity components along x and y axes, respectively. Using Ampere's law, × = , we get the following: omitting the inclusion of electromagnetic waves propagation and neglecting displacement currents, is the magnetic permeability (the nanofluid is assumed to have the same magnetic permeability with that of the base fluid). Using Faraday's law and after algebraic manipulation, the magnetic induction equation is written as follows [37]:  Table 1. Thermophysical properties of water, alumina nanoparticles, and glass spheres, namely, the density (ρ), the heat capacity (C p ), the thermal conductivity (k), the thermal expansion coefficient (β), and the electric conductivity (σ) [34][35][36]. The fluid is hydrodynamically, thermally, and electrically isotropic. No-slip and impermeable boundary conditions are imposed on all boundaries. Joule heating, viscous dissipation, pressure work, and Hall effects are neglected, as they are assumed negligible.
From Ohm's law, the current density, J, relates to the electric field, E, as follows: where σ is the electrical conductivity, and u = (u(x, y), v(x, y)) is the superficial velocity vector, with u and v being the velocity components along x and y axes, respectively. Using Ampere's law, ∇ × B = µ m J, we get the following: omitting the inclusion of electromagnetic waves propagation and neglecting displacement currents, µ m is the magnetic permeability (the nanofluid is assumed to have the same magnetic permeability with that of the base fluid). Using Faraday's law and after algebraic manipulation, the magnetic induction equation is written as follows [37]: Micropolar fluid theory was introduced by Eringen [38] as a generalization of the Navier-Stokes equations. Micropolar theory examines the microstructure of the fluid, along with the inertial characteristics of the particles, which can rotate under the influence of local forces. Microrotation, which is a new kinematic variable introduced in Eringen's model, is related the rotation of the particles. The micropolar vector describes the rotation of the particles contained in a volume element [22,23]. A detailed analysis of the micropolar fluid flow theory is presented in [24].
In the proposed model, the Darcy-Brinkman momentum equations with buoyancy and advective inertia terms are used (see Equation (5)). The quantities that appear in the equations are the classical averaged quantities for flows in porous media. On the other hand, Equation (6), which is the transport equation for the microrotation vector, n, represents quantities which are not averaged directly. It is applied within the pores and, therefore, the correct equation to be used would be an averaged form of these equations. The microrotation vector is affected by the presence of particles and pore walls; however, a single-phase fluid equation in a porous setup was used, incorporating the effective microrotational properties of the particulate phase of the nanofluid, assuming volume-averaged values in porous representative elementary volumes (REVs) and avoiding the incorporation of the details of the pore structure in the expression of the angular momentum equation itself. The same concept applies for the presence of the magnetic field, which should affect the properties of the micropolar nanofluid due to their intrinsic nature; however, the magnetic field is not required to be implemented as a term in the transport equation for n, since its influence on the flow can be regarded solely through the particular term in the linear momentum equation. Thus, the governing equations of the problem, i.e., the conservation of mass, energy, and linear and angular momentum, as well as the induction equation, are written as Equations (4)-(8) [39][40][41].
ρ n f j ∂n ∂t ρC p e f f ∂T ∂t The boundary conditions for the natural convection problem studied are set to the following: p is the pressure, T = T(x, y) is the fluid temperature, n is the superficial microrotation vector, is the induced magnetic field, and g is the acceleration due to gravity. ρ n f is the density, µ n f is the dynamic viscosity, k is the vortex viscosity, γ n f is the spin-gradient viscosity, and j is the micro-inertia density. For the present simulation study, it is further assumed that γ n f has the following form, as proposed in [42,43].
The effective density for the nanofluid is given by an Effective Medium Theory (EFM) approach [30]: with ϕ being the volume fraction of the nanoparticles. The effective dynamic viscosity for the nanofluid is given by the following: which was measured experimentally by Pak and Cho [44]. The heat capacity of the nanofluid, along with the effective heat capacity of the medium, are given as follows [34]: ρC where ε is the porosity of the medium. The effective thermal diffusivity for the nanofluid is as follows: while the thermal expansion coefficient of the nanofluid can be determined as follows: The thermal conductivity of the nanofluid was calculated experimentally [44] and is given by the following: The electrical conductivity of the alumina nanofluid, for the range of volume fractions used (0.005 ≤ ϕ ≤ 0.03), was calculated experimentally [35]: The effective thermal conductivity, k e f f , and the effective electrical conductivity, σ e f f , of the porous medium can be calculated by the Maxwell correlation [39]: The Carman-Kozeny relationship was chosen for the permeability of the medium [36]: where D p is the average diameter of the glass spheres, considered here as D p = 1.6 mm. The total length of the cavity is L = 100 D p = 160 mm.
Since pressure boundary conditions are inherently difficult to be defined with the present meshless method, the velocity-vorticity formulation was used. Taking the curl of the vorticity definition, ω = ∇ × u, and using the continuity equation (Equation (4)), a vector Poisson equation can be used to correlate the vorticity and velocity fields (Equation (22)). The curl operator applies to the linear momentum conservation equation, Equation (5), and with the application of ∇·ω = 0, ∇·n = 0, following the vorticity and microrotation definitions, and the mass conservation ∇·u = 0, one gets Equation (23). Additionally, the magnetic potential, A, is used, given as B = ∇ × A, which satisfies the solenoidal nature of the magnetic field B, (∇·B = 0). Thus, B x = ∂A z /∂y and B y = −∂A z /∂x, where it is assumed that A = A zẑ , withẑ being the z-unit vector, as expected in two-dimensional (x, y) fields. Therefore, vector Equation (7) can be replaced by the scalar Equation (24). Finally, in two-dimensional flows, taking into account all previous assumptions, the governing equations in the velocity-vorticity formulation can be written as follows: The following non-dimensional variables are introduced: for the case of two-dimensional flows. Moreover, considering all previous assumptions, the final form of the equations is written as follows: Appl. Sci. 2020, 10, 1633 8 of 23 The dimensionless boundary conditions are: is the Hartmann number, and Da = K D L 2 is the Darcy number. The Nusselt number can be expressed as follows: while the average Nusselt number is defined as follows:

Numerical Method
The steady-state flow equations were numerically solved by using a meshless point collocation method [45,46]. The intermediate velocity, U * , computed from the Poisson equations for the velocity components (Equations (28) and (29)), did not satisfy the continuity equation. Therefore, the velocity field was updated, using the velocity correction method [47,48]. Assuming that the velocity correction (∂U) is irrotational, a correction potential, ψ n+1 , can be defined as ∇ψ n+1 = ∂U n+1 . The correction potential should satisfy a Poisson-type equation to satisfy the continuity equation.
The updated velocity, U n+1 , is calculated as U n+1 = U * + ∂U n+1 . Following this, we briefly describe the solution procedure applied, placing emphasis on the calculation of the convective terms and on the imposition of Neumann boundary conditions.

Solution Procedure
The governing equations for the magnetic potential, the temperature, the microrotation and the vorticity (Equations (28)-(33)) can be expressed compactly as follows: where L is a differential operator; q accounts for A, θ, N, and Ω; and f is the right-hand side of the equation. All the field functions considered in Equation (38) are written in bold, since they are column vectors of dimensions N p × 1, and N p is the total number of nodes (interior and boundary) used to represent the spatial domain. The notation, n, defines the iteration step. For easier convergence, the convective term in the n + 1 step is expressed as U(∇·q) = U n (∇·q) n+1 + U n − U n+1 (∇·q n ). For the variables under consideration, the L operator and the right-hand side vector are given by the following equations: The solution procedure commences with an initial guess for all the field variables, U (0) , V (0) , Ω (0) , N (0) , A (0) , and Θ (0) . The algorithmic steps in the solution procedure may be ordered as follows:

•
Step 1: Compute the intermediate velocity components, U * and V * , through Equations (28) and (29), using vorticity values at the previous step (for the first iteration, the prescribed initial conditions are used).

•
Step 2: Compute the correction potential and update the velocity field (U n+1 and V n+1 ).

•
Step 4: Solve the coupled equations for the microrotation and vorticity, Equations (32) and (33), using the updated magnetic field and temperature gradient.
The Discretization Corrected Particle Strength Exchange (DC PSE) method was implemented for the calculation of the derivatives [49] of the unknown field functions up to the second order. DC PSE was originally introduced as a Lagrangian particle-based solution method. The method has been reformulated [47] to numerically solve partial differential equations in Eulerian frameworks. The DC PSE interpolation method relies on Particle Strength Exchange (PSE) operators, i.e., kernel functions that approximate differential operators which conserve particle strength in particle-particle interactions. The PSE method was first introduced by Degond and Mas-Gallic [50] for diffusion and convection-diffusion problems, while Eldredge et al. [51] extended the method to compute approximations of arbitrary derivatives. Details on the DC PSE method can be found in [47].

Vorticity Boundary Conditions
Several formulae, based on finite difference methods, for the imposition of vorticity boundary values [52], usually referred to as local, were developed. Vorticity values on the boundaries were computed, using vorticity values of interior nodes, avoiding the involvement of boundary nodes. The most commonly used formula was Thom's formula. Unfortunately, the application of these formulas in physical applications had limited success.
Global vorticity boundary conditions are more accurate, since both boundary and interior nodes are involved in their computation. In the present study, the DC PSE differential operators were applied for the computation of the vorticity boundary conditions [45,53]. In the DC-PSE-meshless collocation context, computing and imposing vorticity boundary conditions were considered as an interpolation problem. Using the updated velocity components from Step 2, we computed U n+1 and V n+1 , the updated vorticity values, ω n+1 (in the entire spatial domain, including the boundaries): The vorticity values computed are used as boundary condition for Step 4, Equation (33).

Zero Flux Temperature Boundary Conditions
Frequently, for natural convection flow problems, Neumann boundary conditions (no-flux boundary conditions) are encountered, given as follows: where n is the unit normal vector of the boundary, and q is the unknown scalar field function at the boundary (the temperature T in the present case). We imposed such Neumann boundary conditions, using interior nodes, only in the support domain of the nodes with Neumann boundary conditions, shown in Figure 2, when computing the derivatives in Equation (44).
Appl. Sci. 2020, 10, 1633 10 of 22 Figure 2. Outline of the support domain used for imposing no-flux boundary conditions. A boundary node is denoted with a red circle, and the corresponding nodes in its support domain with green squares.

Results and Discussion
For the flow case studies considered here, the microrotation number, , is set to = 2, since it provides results that are closer to the experimental findings [25]. The Prandtl number is set to = 6.2, and the magnetic Prandtl to = 10 −5 .

Code Verification
In the present simulations, uniform Cartesian grids with increasing nodal resolution are used, ranging from 51 × 51 up to 201 × 201. A grid-independent numerical solution was obtained by a grid of 101 × 101 resolution, as seen in Figure 3. The average Nusselt number ( ) on the left vertical wall (hot wall) was computed for validation purposes. The average Nusselt number for the natural convection in the square porous cavity is presented for = 100, = 0.7, = 0.01, and two different values of = 5 10 7 and = 5 10 8 . For higher accuracy, a denser grid consisting of 151 × 151 nodes was used. The accuracy of the proposed scheme was tested against verified Boundary nodes having only interior nodes in their support domain allow the direct computation of the values on the boundary nodes based on known values of interior nodes and, thus, decoupling the Neumann boundary conditions. For a two-dimensional system with unit normal vector n = (n x , n y ), the derivatives that appear in the Neumann boundary conditions are estimated by using the DC PSE method, with the following expression at each boundary node, j: with w x i and w y i being the weights of the first-order spatial derivatives. The Neumann boundary condition for node j becomes the following: allowing the computation of the value that has to be imposed at the boundary node as follows:

Results and Discussion
For the flow case studies considered here, the microrotation number, K, is set to K = 2, since it provides results that are closer to the experimental findings [25]. The Prandtl number is set to Pr = 6.2, and the magnetic Prandtl to Pr m = 10 −5 .

Code Verification
In the present simulations, uniform Cartesian grids with increasing nodal resolution are used, ranging from 51 × 51 up to 201 × 201. A grid-independent numerical solution was obtained by a grid of 101 × 101 resolution, as seen in Figure 3. The average Nusselt number (Nu ave ) on the left vertical wall (hot wall) was computed for validation purposes. The average Nusselt number for the natural convection in the square porous cavity is presented for Ha = 100, ε = 0.7, ϕ = 0.01, and two different values of Ra = 5 × 10 7 and Ra = 5 × 10 8 . For higher accuracy, a denser grid consisting of 151 × 151 nodes was used. The accuracy of the proposed scheme was tested against verified numerical results obtained for MHD natural convection in a square cavity [54] (Table 2). Table 2. Average Nusselt numbers at the hot wall for characteristic Grashof and Hartmann numbers obtained by the present method (underlined) and [54].

Effect of Porosity (ε)
Initially, the flow regime in the case of varying porosity of the medium was examined. The porosity is strongly related to the permeability of the system (Equation (21)). By increasing porosity, the permeability and the Darcy number of the porous medium increases. In the present simulations, the volume fraction of nanoparticles was set to = 0.01, the Rayleigh number to = 5 10 8 , and the Hartmann number to = 0. In all porosity values considered, a central clockwise-rotating flow pattern was formed along the diagonal of the square enclosure, having two anticlockwise circulating regions in both the upperleft and lower-right parts of the diagonal, as shown in Figure 4. Three flow regions were formed, a hot region bounded between the left and top walls, a cold region between the right and bottom walls, and a third one along the diagonal area where temperature varies. By increasing porosity, the flow   [54] and the results of the proposed method. In these cases, the medium consisted only of the fluid (ε = 1, Da = ∞, ϕ = 0), the Hartmann number ranged from Ha = 0 to Ha = 100, and the Grashof number varied in the range 10 3 ≤ Gr ≤ 10 5 . The Grashof number is inherently related to the Rayleigh and Prandtl numbers through the relationship Gr = Ra/Pr. The comparative results are in good agreement.

Effect of Porosity (ε)
Initially, the flow regime in the case of varying porosity of the medium was examined. The porosity is strongly related to the permeability of the system (Equation (21)). By increasing porosity, the permeability and the Darcy number of the porous medium increases. In the present simulations, the volume fraction of nanoparticles was set to ϕ = 0.01, the Rayleigh number to Ra = 5 × 10 8 , and the Hartmann number to Ha = 0.
In all porosity values considered, a central clockwise-rotating flow pattern was formed along the diagonal of the square enclosure, having two anticlockwise circulating regions in both the upper-left and lower-right parts of the diagonal, as shown in Figure 4. Three flow regions were formed, a hot region bounded between the left and top walls, a cold region between the right and bottom walls, and a third one along the diagonal area where temperature varies. By increasing porosity, the flow patterns intensify (higher stream function values), and the anticlockwise rotating cells tend to dominate the flow. Isotherms expand gradually, as stronger circulation favors concentration of isotherms lines. The microrotation isolines form a central vortex diagonally, increasing in length as the porosity increases. Vorticity variations are observed at the upper-right and lower-left edges of the cavity, whereas their thickness decreases and the intensity of the vorticity increases with increasing porosity.
patterns intensify (higher stream function values), and the anticlockwise rotating cells tend to dominate the flow. Isotherms expand gradually, as stronger circulation favors concentration of isotherms lines. The microrotation isolines form a central vortex diagonally, increasing in length as the porosity increases. Vorticity variations are observed at the upper-right and lower-left edges of the cavity, whereas their thickness decreases and the intensity of the vorticity increases with increasing porosity.

Effect of the Hartmann Number (Ha)
Next, the flow regime in the case of different Hartmann numbers was investigated, related with the intensity of the applied magnetic field. The volume fraction was set to φ = 0.03, the Rayleigh number Ra = 5 10 8 , and the porosity to ε = 0.7.
After the onset of the magnetic field, notable variations in the stream function, the microrotation, and the temperature contours were observed. The anticlockwise rotating vortex formations gradually disappeared, and two secondary rotating cells (having low stream function values) appeared initially along the diagonal, shifting horizontally and disappearing as the Hartmann number increased ( Figure 6). Boundary layers (Hartmann layers) were established on the vertical walls of the cavity. Microrotation contours formed a vortex that moved toward the center of the cavity as the Hartmann number increased. Initially, the isotherms swirled according to the flow field, and as the Hartmann number increased, they became parallel to the vertical walls. The vorticity formed two vortices at the vertical walls for high values of the Hartmann number.
A decrease in the heat transfer was observed, along with an increase in the intensity of the magnetic field, through the variation of the Nuave along the hot wall (Figure 7). The average Nusselt number decreased with the increase in the Hartmann number. The decrease of the average Nusselt number was more evident for elevated values of the Rayleigh number.

Effect of the Hartmann Number (Ha)
Next, the flow regime in the case of different Hartmann numbers was investigated, related with the intensity of the applied magnetic field. The volume fraction was set to ϕ = 0.03, the Rayleigh number Ra = 5 × 10 8 , and the porosity to ε = 0.7.
After the onset of the magnetic field, notable variations in the stream function, the microrotation, and the temperature contours were observed. The anticlockwise rotating vortex formations gradually disappeared, and two secondary rotating cells (having low stream function values) appeared initially along the diagonal, shifting horizontally and disappearing as the Hartmann number increased ( Figure 6). Boundary layers (Hartmann layers) were established on the vertical walls of the cavity. Microrotation contours formed a vortex that moved toward the center of the cavity as the Hartmann number increased. Initially, the isotherms swirled according to the flow field, and as the Hartmann number increased, they became parallel to the vertical walls. The vorticity formed two vortices at the vertical walls for high values of the Hartmann number. Appl. Sci. 2020, 10, 1633 14 of 22

Effect of the Rayleigh Number (Ra)
Cases with varying Rayleigh number were also under investigation. For these simulations, the volume fraction was kept at = 0.03, while the porosity was kept at = 0.7. The Hartmann number was set to = 100. For the low Rayleigh values considered, a large clockwise rotating vortex was formed, as shown in Figure 8. The fluid rose from the left vertical wall, the location of the heat source, and flowed downward along the right vertical wall. By increasing the Rayleigh number, isotherms expanded gradually, as they emerged from the heat source and reached the rest of the walls. Furthermore, the flow regime became stronger, and the streamlines were intensified. For higher values of the Rayleigh number ( = 5 10 8 ), the two anticlockwise rotating formations appeared in the stream function contour, and the distribution of the isothermal lines changed from conduction-controlled to convection-dominated. For small values of the Rayleigh number, a minute variation on the vorticity was observed near the walls of the cavity. Microrotation contours had a central vortex at the middle of the cavity that elongated diagonally for elevated values of the Rayleigh number.

Effect of the Rayleigh Number (Ra)
Cases with varying Rayleigh number were also under investigation. For these simulations, the volume fraction was kept at ϕ = 0.03, while the porosity was kept at ε = 0.7. The Hartmann number was set to Ha = 100.
For the low Rayleigh values considered, a large clockwise rotating vortex was formed, as shown in Figure 8. The fluid rose from the left vertical wall, the location of the heat source, and flowed downward along the right vertical wall. By increasing the Rayleigh number, isotherms expanded gradually, as they emerged from the heat source and reached the rest of the walls. Furthermore, the flow regime became stronger, and the streamlines were intensified. For higher values of the Rayleigh number (Ra = 5 × 10 8 ), the two anticlockwise rotating formations appeared in the stream function contour, and the distribution of the isothermal lines changed from conduction-controlled to convection-dominated. For small values of the Rayleigh number, a minute variation on the vorticity was observed near the walls of the cavity. Microrotation contours had a central vortex at the middle of the cavity that elongated diagonally for elevated values of the Rayleigh number. , on is shown in Figure 9 for two values, and two porosities of the medium, while the volume fraction was set to = 0.03. The increase of resulted in an increase of . For = 100, a smaller increase of the Nusselt number was observed, as the presence of the magnetic field had a negative effect on the strength of the flow field. The dependence of the average Nusselt number, Nu ave , on Ra is shown in Figure 9 for two Ha values, and two porosities of the medium, while the volume fraction was set to ϕ = 0.03. The increase of Ra resulted in an increase of Nu ave . For Ha = 100, a smaller increase of the Nusselt number was observed, as the presence of the magnetic field had a negative effect on the strength of the flow field.

Effect of the Volume Fraction (φ)
The influence of the volume fraction of the Al2O3 nanoparticles on the flow regime and the heattransfer properties of the nanofluid was under examination. The Hartmann number was set equal to = 100 , the Rayleigh number to = 5 × 10 8 , and the porosity to = 0.7 . In the present simulations, the nanoparticles volume fraction, , ranges from = 0 to = 0.03.
The volume fraction of the nanoparticles added in the base fluid alters the heat-transfer properties of the nanofluid. When nanoparticles were added, the thermal and electrical conductivity of the base fluid (water) increased. Even with low (1%-5% volumetric) concentration of the nanoparticles, the thermal conductivity increased by 10%-20%.
As the nanoparticles volume fraction increased, the intensity of the flow regime became weaker, and that of the flow field decreased, Figure 10. This is attributed to the increase of the electrical conductivity, as the effect of the applied magnetic field on the nanoparticles becomes stronger. The temperature contours follow the alteration of the flow field, and the microrotation isolines form a vortex at the center of the cavity.
The dependence of Nuave on the solid volume fraction is shown in Figure 11. At low-flow conditions, ≤ 5 10 6 , conduction is the dominant mechanism for heat transfer, and the Nusselt number increases as the volume fraction of the particles increases, as expected by the increased thermal conductivity offered by the nanoparticles with respect to the base fluid thermal properties. As the flow conditions intensify, convection is at play, and increases, as well, provided the volumetric concentration of nanoparticles is constant. However, at elevated flow conditions, = 5 10 8 , even in the absence of a magnetic field, the Nusselt number is reduced as the volume fraction of the nanoparticles increases, and this slope declines as the Rayleigh number increases. As the Hartmann number increases, the Nusselt number recedes further.
The decrease in the heat transfer with the addition of nanoparticles under elevated flow conditions was not an anticipated result, given the intended purpose of nanofluids to increase the efficiency of thermal processes. However, similar behavior has been observed and reported in previous experimental [36] and numerical [21] investigations of the thermal behavior of buoyancy-

Effect of the Volume Fraction (ϕ)
The influence of the volume fraction of the Al 2 O 3 nanoparticles on the flow regime and the heat-transfer properties of the nanofluid was under examination. The Hartmann number was set equal to Ha = 100, the Rayleigh number to Ra = 5 × 10 8 , and the porosity to ε = 0.7. In the present simulations, the nanoparticles volume fraction, ϕ, ranges from ϕ = 0 to ϕ = 0.03.
The volume fraction of the nanoparticles added in the base fluid alters the heat-transfer properties of the nanofluid. When nanoparticles were added, the thermal and electrical conductivity of the base fluid (water) increased. Even with low (1%-5% volumetric) concentration of the nanoparticles, the thermal conductivity increased by 10%-20%.
As the nanoparticles volume fraction increased, the intensity of the flow regime became weaker, and that of the flow field decreased, Figure 10. This is attributed to the increase of the electrical conductivity, as the effect of the applied magnetic field on the nanoparticles becomes stronger. The temperature contours follow the alteration of the flow field, and the microrotation isolines form a vortex at the center of the cavity.
presence of the porous medium. The nanofluid moved with additional hindrance through the pores, due to the quadratic increase of the effective fluid viscosity (Equation (11)), and an increase of the volume fraction reduced the fluid flow and, thus, the convection. Therefore, even though the thermal conductivity increased, the convective heat transfer was reduced, and the Nusselt number decreased. Moreover, in the case of an applied magnetic field, the decrease was steeper, since the magnetic field had a negative effect on the intensity of the flow field.  The dependence of Nu ave on the solid volume fraction is shown in Figure 11. At low-flow conditions, Ra ≤ 5 × 10 6 , conduction is the dominant mechanism for heat transfer, and the Nusselt number increases as the volume fraction of the particles increases, as expected by the increased thermal conductivity offered by the nanoparticles with respect to the base fluid thermal properties. As the flow conditions intensify, convection is at play, and Nu ave increases, as well, provided the volumetric concentration of nanoparticles is constant. However, at elevated flow conditions, Ra = 5 × 10 8 , even in the absence of a magnetic field, the Nusselt number is reduced as the volume fraction of the nanoparticles increases, and this slope declines as the Rayleigh number increases. As the Hartmann number increases, the Nusselt number recedes further. Figure 11. Average Nusselt number as a function of the volume fraction of nanoparticles, for characteristic Hartmann and Rayleigh numbers, and = 0.6.

Conclusions
In the present work, the laminar, natural convective flow of an Al2O3-water micropolar nanofluid in a rectangular porous enclosure under the presence of a magnetic field was numerically investigated. The micropolar flow theory was applied and extended to study nanofluidic natural convection phenomena subjected to magnetic fields in the interior of porous media. The nanofluid was assumed hydrodynamically, thermally, and electrically isotropic. Moreover, it was assumed to be in constant thermal equilibrium with the porous medium. The governing equations were solved in the stream function-vorticity description. A robust, efficient, and accurate meshless strong-form solver was presented, to numerically solve the governing set of coupled partial differential equations. The steady-state equations were solved, using a meshless point collocation method for the discretization, and a five-step algorithm was developed. A reformulated DC PSE method was employed to calculate the derivatives of the unknown field functions up to the second order, in a Eulerian framework, with increased accuracy.
The flow characteristics and, thus, the convective heat transfer inside the enclosure were determined by the intensity of the magnetic field, the volume fraction of the nanoparticles, the porosity, and the strength of the flow field. Circulation and convection phenomena intensified with increasing Rayleigh numbers, though they were suppressed considerably by the presence of a strong magnetic field. Formation of rotating vortexes greatly influenced the temperature distribution. The average Nusselt number increased substantially with the Rayleigh number, since the circulation grew stronger. The magnetic field reduced the local Nusselt number considerably by the suppression of the convection currents. The porosity of the medium also affected the convective heat transfer. The average Nusselt number increased with increasing porosity, as the flow intensity increased. The addition of nanoparticles had either a negative or positive effect on the heat transfer, depending on the intensity of the flow field. For small Rayleigh numbers, the heat transfer increased due to a higher thermal conductivity. For large Rayleigh numbers, the decrease in the heat transfer was caused by the elevated viscosity of the nanofluid, which reduces the fluid flow and the heat convection The decrease in the heat transfer with the addition of nanoparticles under elevated flow conditions was not an anticipated result, given the intended purpose of nanofluids to increase the efficiency of thermal processes. However, similar behavior has been observed and reported in previous experimental [36] and numerical [21] investigations of the thermal behavior of buoyancy-driven nanofluids in rectangular enclosures. In [36], the Nusselt number of the base fluid in a porous medium was found to be higher than that of the nanofluid for elevated Rayleigh numbers. In [21], the Nusselt number decreased (~10-30%) with the addition of nanoparticles for high Rayleigh numbers. In the present case, the phenomenon was more intense (~50%) due to the additional presence of the porous medium. The nanofluid moved with additional hindrance through the pores, due to the quadratic increase of the effective fluid viscosity (Equation (11)), and an increase of the volume fraction reduced the fluid flow and, thus, the convection. Therefore, even though the thermal conductivity increased, the convective heat transfer was reduced, and the Nusselt number decreased. Moreover, in the case of an applied magnetic field, the decrease was steeper, since the magnetic field had a negative effect on the intensity of the flow field.

Conclusions
In the present work, the laminar, natural convective flow of an Al 2 O 3 -water micropolar nanofluid in a rectangular porous enclosure under the presence of a magnetic field was numerically investigated. The micropolar flow theory was applied and extended to study nanofluidic natural convection phenomena subjected to magnetic fields in the interior of porous media. The nanofluid was assumed hydrodynamically, thermally, and electrically isotropic. Moreover, it was assumed to be in constant thermal equilibrium with the porous medium. The governing equations were solved in the stream function-vorticity description. A robust, efficient, and accurate meshless strong-form solver was presented, to numerically solve the governing set of coupled partial differential equations. The steady-state equations were solved, using a meshless point collocation method for the discretization, and a five-step algorithm was developed. A reformulated DC PSE method was employed to calculate the derivatives of the unknown field functions up to the second order, in a Eulerian framework, with increased accuracy.
The flow characteristics and, thus, the convective heat transfer inside the enclosure were determined by the intensity of the magnetic field, the volume fraction of the nanoparticles, the porosity, and the strength of the flow field. Circulation and convection phenomena intensified with increasing Rayleigh numbers, though they were suppressed considerably by the presence of a strong magnetic field. Formation of rotating vortexes greatly influenced the temperature distribution. The average Nusselt number increased substantially with the Rayleigh number, since the circulation grew stronger. The magnetic field reduced the local Nusselt number considerably by the suppression of the convection currents. The porosity of the medium also affected the convective heat transfer. The average Nusselt number increased with increasing porosity, as the flow intensity increased. The addition of nanoparticles had either a negative or positive effect on the heat transfer, depending on the intensity of the flow field. For small Rayleigh numbers, the heat transfer increased due to a higher thermal conductivity. For large Rayleigh numbers, the decrease in the heat transfer was caused by the elevated viscosity of the nanofluid, which reduces the fluid flow and the heat convection mechanisms. This decrease in the heat transfer of the nanofluids under elevated flow conditions was not an anticipated result, given the intended purpose of nanofluids to increase the efficiency of thermal processes; however, it was corroborated by previous experimental [36] and numerical [21] observations. Furthermore, the very presence of the porous medium affected the heat transfer noticeably, as an increase in the effective viscosity of the fluid further impeded the nanofluid motion through the pores. The presence of the magnetic field intensified the phenomenon.