1. Introduction

Natural convection heat transfer has taken great deal of interest in the last decade due to its various engineering applications such as thermal insulation, cooling of electronic equipments, heat exchangers and solar collectors. In order to improve hydrothermal behaviour of the conventional fluids, nanofluids are proposed in thermal systems ^[3].

A lot of numerical studies on nanofluids are reported. Khanafer et al. ^[10] used finite volume method (FVM) with alternating direction implicit procedure analyzing heat transfer performance with a model considering the solid particle dispersion. Their conclusion shows that heat transfer rate increases with the increase in solid volume fraction at any Grashof number. Tiwari et al. ^[22] performed FVM in a Cu-water nanofluid-filled cavity with moving walls in different directions. In their study, heat transfer is more reduced in case of walls moving upward. The same problem of mixed convection flow in different aspect ratios is also considered by Muthtamilselvan et al. ^[14] employing FVM with SIMPLE algorithm on a staggered grid. The linear relation between solid volume fraction and average Nusselt number is revealed. Sheremet et al. ^[21] applied the finite difference method to simulate the free convection in an inclined wavy enclosure filled with a Cu-water nanofluid under the effect of a uniform magnetic field and isothermal corner heater. Sheikholeslami ^[19] analyzes nanofluid flow and heat transfer between heated and permeable parallel plates employing Runge-Kutta integration scheme, and considers Koo-Klein-streuer-Li (KKL) model. Results show that the Nusselt number has a reverse relation with expansion ratio while it has a direct relationship with power law index. Geridonmez Pekmen ^[16] solved natural convection flow in a cavity filled with nanofluid using RBF-PS in space derivatives and differential quadrature method in time derivatives. In this study, different type of nanofluids are considered as well as using both multiquadric and inverse multiquadric RBFs.

Ferrofluids, which are magnetic nanofluids, are used in various fields such as electronic packing, mechanical engineering, thermal engineering, aerospace and bioengineering. Many researches have studied numerically and experimentally on ferrofluids in different geometries either in presence of external magnetic field or in absence of magnetic field. Aminfar et al. ^[1] utilized the control volume technique with SIMPLEC for simulation of 3D laminar ferrofluid in presence of an electric current through a wire at the bottom of the duct resulting with transverse nonuniform magnetic field. Ghasemi et al. ^[7] used FVM for simulation of a water based ferrofluid in a mini channel under the effect of constant or alternating magnetic field induced by wires carrying current. They found that heat transfer enhancement occurs with a constant magnetic field at magnetic number $1.07\times 10^8$, and Reynolds number $25$. The influence of heat dissipation and an external magnetic source on natural convection flow in a kerosene based cobalt-filled cavity is analyzed in ^[8] performing LBM. At a large Rayleigh number $Ra = 10^5$, heat transfer is decelarated by the increase in solid volume fraction of nanoparticles. LBM is also used to simulate the same problem in an inclined cavity in ^[9]. The effect of the increase in solid volume fraction on heat transfer dominates over the effect of inclination angle. In ^[18], governing equations based on ferrohydrodynamic (FHD) effect and magnetohydrodynamic (MHD) effect are solved in a semi annulus enclosure with sinusoidal hot wall with a magnetic source employing control volume based finite element method (CVFEM). The results demonstrate that the effect of Kelvin force is pronounced at a low Rayleigh number.

Nanofluids in porous medium are encountered in ^[12,20]. In ^[12], Cu-water nanofluid in a porous medium with Darcy-Brinkman-Forchheimer model considering the magnetic field is solved by FVM with SIMPLE algorithm. They reported that the heat transfer rate increases with the augmentation in Darcy number. Sheikholeslami ^[20] analyzed the nanofluid flow and heat transfer in KKL model in a stretching porous cylinder conducting the fourth order Runge Kutta integration scheme with a shooting technique. Results show that solid volume fraction and skin friction coefficient have a reverse relation, and average Nusselt number increases with the increase in the suction parameter.

Geridonmez Pekmen has also solved the natural convection in a square cavity involving magnetite nanoparticles inside of the nanofluid and a magnetic source considering an exceptional case as $\nabla \cdot H \neq 0$ in ^[17]. In that study, the determination of shape parameter in view of minimum radial distances in a different grid distribution, and the usage of perturbed-like multiquadric RBFs are focused. In the present study, natural convection flow in a porous Fe$_3$O$_4$-water-filled cavity in presence of a magnetic source is investigated considering the effects of Kelvin forces and non-Darcian in governing dimensionless equations. The Brinkman-Forchheimer-extended Darcy model is simulated utilizing the radial basis function based pseudo spectral (RBF-PS) method with multi quadric (MQ) RBFs. A way to determine the shape parameter in MQ RBF is also proposed. The effects of physical non-dimensional parameters are investigated. To the author's knowledge, this method is firstly applied to this problem.

---

2. Physical and mathematical background

The problem configuration is depicted in Figure 1. No-slip boundary conditions for velocity are imposed. The left hot wall and the cold right wall are maintained at $T_h = 1$ and $T_c = 0$, respectively. The top and bottom walls are adiabatic ($\partial T/\partial n = 0$). A magnetic source close to the left wall is placed.

Fe₃O₄-water nanofluid flow is unsteady, laminar and incompressible. The problem is two-dimensional and unsteady. The nanoparticles and the base fluid (water) are in thermal equilibrium. The porous medium is homogeneous, isotropic, electrically insulated, and the fluid and solid of porous medium are also in thermal equilibrium. Thermophysical properties of water and magnetite are given in Table 1. The constant thermophysical properties of the fluid are considered except the density variation treated by Boussinessq approximation. The radiation effect, induced currents and Hall current are neglected.

The magnetic field intensity and its components are defined as ^[23]

$\overline{H}_x = \frac{\gamma}{2\pi}\frac{y-\bar{b}}{ (\bar b - y)^2 + (\bar a - x )^2},$

(2.1a)

$\overline{H}_y = \frac{\gamma}{2\pi}\frac{\bar{a} - x}{ (\bar b - y)^2 + (\bar a - x )^2}$

(2.1b)

$\overline{H} = \sqrt{\overline{H}_x^2 + \overline{H}_y^2}$

(2.1c)

where $(\bar a, \bar b)$ is the location of the magnetic source and $\gamma$ is the strength of the magnetic field at the source.

Under the aforementioned assumptions, the physical governing equations are

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

(2.2a)

$\frac{\mu_{nf}}{\rho_{nf}\varepsilon} \nabla^2 u = \frac{1}{\varepsilon}\frac{\partial u}{\partial t} + \frac{1}{\varepsilon^2}\bigg( u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} \bigg) + \frac{1}{\rho_{nf}}\frac{\partial p}{\partial x} + \frac{\mu_0^2 }{\rho_{nf}} \overbrace{\sigma_{nf} \left(\overline{H}_y^2 u - \overline{H}_x\overline{H}_yv \right)}^\text{A} \nonumber \\ \quad \quad \quad \quad \hspace{3cm}-\frac{1}{\rho_{nf}}\underbrace{\mu_0 M \frac{\partial \overline{H}}{\partial x}}_\text{B} + \underbrace{\frac{\mu_{nf}}{\kappa\rho_{nf}} u + \frac{c_F}{\sqrt{\kappa}}|{\bf u}|u}_\text{F}$

(2.2b)

$\frac{\mu_{nf}}{\rho_{nf}\varepsilon} \nabla^2 v = \frac{1}{\varepsilon}\frac{\partial v}{\partial t} + \frac{1}{\varepsilon^2}\bigg( u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} \bigg) + \frac{1}{\rho_{nf}}\frac{\partial p}{\partial y} + \frac{\mu_0^2 }{\rho_{nf}} \overbrace{\sigma_{nf} \left(\overline{H}_x^2 v - \overline{H}_x\overline{H}_y u \right)}^\text{A'} \nonumber \\ \quad \quad \quad \hspace{1cm}-\frac{1}{\rho_{nf}} \underbrace{\mu_0 M \frac{\partial \overline{H}}{\partial y}}_\text{B'} +\underbrace{ \frac{\mu_{nf}}{\kappa\rho_{nf}} v + \frac{c_F}{\sqrt{\kappa}}|{\bf u}|v}_\text{F'} -\underbrace{g\beta_{nf}(T-T_c)}_\text{G}$

(2.2c)

$k_{nf}\nabla^2 T = (\rho c_p)_{nf}\bigg(\frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y}\bigg) - \mu_0^2 \overbrace{\sigma_{nf}(\overline{H}_x v - \overline{H}_y u)^2 }^\text{C} \\ \quad \quad + \underbrace{\mu_0 T \frac{\partial M}{\partial T} \left( u\frac{\partial \overline{H}}{\partial x} + v \frac{\partial \overline{H}}{\partial y} \right)}_\text{D} - \mu_{nf}\underbrace{\left[2\left(\frac{\partial u}{\partial x}\right)^2 + 2\left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)^2 \right]}_\text{E},$

(2.2d)

where A, A' are Lorentz force components, B, B' are Kelvin force components, C is the Joule heating, D is thermal power per unit volume resulting from the magneto-caloric effect, E is the viscous dissipation, F, F' are Brinkman term plus Forchheimer term, G is the buoyancy force term and ^[2,10,13]

$\begin{align*} & M = K' \overline{H} (T_c' - T)\ \mbox{ is the magnetization }, \\ & c_F = \frac{1.75(1-\varepsilon)}{d_p\ \varepsilon^3}\ \mbox{ is the form coefficient}, \\ & \kappa = \frac{d_p^2\ \varepsilon^3}{150(1-\varepsilon)^2}\ \mbox{ is the permeability of the porous medium}, \\ & \mu_0 = 4\pi \times 10^{-7}\ Tm/A \, \mbox{ is the magnetic permeability of vacuum}, \\ & \rho_{nf} = \rho_f(1-\phi) + \rho_s \phi\, \mbox{ is the density of } \\ & (\rho c_p)_{nf} = (\rho c_p)_f (1-\phi) + (\rho c_p)_s \phi\, \mbox{is the heat capacitance of} \\ & \mu_{nf} = \frac{\mu_f}{(1-\phi)^{2.5}}\, \mbox{ is the dynamic viscosity of} \\ & \beta_{nf} = \beta_f (1-\phi) + \beta_s \phi\, \mbox{ is the thermal expansion coefficient of} \\ & \alpha_{nf} = \frac{k_{nf}}{(\rho c_p)_{nf}}\, \mbox{ is the thermal diffusivity of} \\ & k_{nf} = k_f \left( \frac{k_s - 2\phi(k_f-k_s)+2k_f}{k_s+\phi(k_f-k_s) + 2k_f} \right)\, \mbox{ is the thermal conductivity of}\\ & \sigma_{nf} = \sigma_f \left( \frac{3(\sigma_s/\sigma_f -1)\phi}{(\sigma_s/\sigma_f+2)-(\sigma_s/\sigma_f-1)\phi} +1 \right)\, \mbox{is the electrical conductivity of} \end{align*}$

nanofluid. Also, $T_c'$ is the Curie temperature, $T$ is the fluid temperature, $K'$ is a constant, $\varepsilon$ is the porosity of the porous medium, $d_p$ is the solid particle size of the porous medium, and $\phi$ is the nanoparticle volume fraction. Subindices $f, s$ and $nf$ correspond to fluid, solid and nanofluid, respectively.

In order to obtain the dimensionless equations, the non-dimensional parameters are defined as

$\begin{align} (a', b') & = \frac{(\bar a, \bar b)}{L}, \quad (x', y') = \frac{(x, y)}{L}, \quad (H_x', H_y', H') = \frac{(\overline{H}_x, \overline{H}_y, \overline{H})}{\overline{H}_0}, \nonumber \\ p'& = \frac{pL^2}{\rho_f \alpha_f^2}, \quad t' = \frac{tL^2}{\alpha_f}, \quad u' = \frac{uL}{\alpha_f}, \quad v' = \frac{vL}{\alpha_f}, \quad T' = \frac{T - T_c}{T_h-T_c}, \end{align}$

(2.3)

where $L$ is the characteristic length, $\overline{H}_0 = \overline{H}(\bar{a}, 0) = \frac{\gamma}{2\pi |b|}$.

These parameters in Eq.(2.3) are put into the dimensional equations Eq.(2.2). Then, the prime notations are dropped and dimensionless equations in $u-v-p-T$ form are derived as

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

(2.4a)

$\begin{eqnarray*} \nabla^2 u & = & \frac{1}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f \mu_{nf}} \left( \frac{\partial u}{\partial t} + \frac{1}{\varepsilon}\bigg(u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y}\bigg)\right) + \frac{\varepsilon}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f\mu_{nf}} \frac{\partial p}{\partial x} + \varepsilon Ha^2\frac{\sigma_{nf}\mu_f}{\sigma_f\mu_{nf}}(H_y^2 u - H_xH_y v) \nonumber \\ & +& \frac{\varepsilon}{Da}u + \frac{c_f}{\sqrt{Da}}\frac{\varepsilon}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f\mu_{nf}} |{\bf u}|u - \varepsilon Mn_F \frac{\mu_f}{\mu_{nf}}(\epsilon_2-\epsilon_1-T)H\frac{\partial H}{\partial x} \end{eqnarray*}$

(2.4b)

$\begin{eqnarray*} \nabla^2 v & = & \frac{1}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f \mu_{nf}} \left( \frac{\partial v}{\partial t} + \frac{1}{\varepsilon}\bigg(u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y}\bigg)\right) + \frac{\varepsilon}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f\mu_{nf}}\frac{\partial p}{\partial y} + \varepsilon Ha^2\frac{\sigma_{nf}\mu_f}{\sigma_f\mu_{nf}}(H_x^2 v - H_xH_y u) \nonumber \\ & +& \frac{\varepsilon}{Da}v + \frac{c_f}{\sqrt{Da}}\frac{\varepsilon}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f\mu_{nf}} |{\bf u}|v -\varepsilon Ra \frac{\rho_{nf}\mu_f \beta_{nf}}{\rho_f \mu_{nf} \beta_f} T - \varepsilon Mn_F \frac{\mu_f}{\mu_{nf}}(\epsilon_2-\epsilon_1-T)H\frac{\partial H}{\partial y} \end{eqnarray*}$

(2.4c)

$\begin{eqnarray*} \nabla^2 T & = & \frac{(\rho c_p)_{nf} k_f}{(\rho c_p)_f k_{nf}} \left( \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y}\right) - Ha^2 Ec \frac{\sigma_{nf}k_f}{\sigma_f k_{nf}} (H_x v - H_y u) ^2 \nonumber \\ & -& Ec \frac{\mu_{nf}k_f}{\mu_f k_{nf}} \left[2\left(\frac{\partial u}{\partial x}\right)^2 + 2\left(\frac{\partial v}{\partial y}\right)^2 + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right)^2 \right] \label{nondimT} \nonumber \\ & -& Mn_F Ec \frac{k_f}{k_{nf}} \left(u \frac{\partial H}{\partial x} + v \frac{\partial H}{\partial y}\right) H (\epsilon_1 + T) \end{eqnarray*}$

(2.4d)

in which $|{\bf u}| = \sqrt{u^2 + v^2}, \ c_f = \frac{1.75}{\sqrt{150}\varepsilon^{1.5}}$, and the dimensionless parameters Darcy, Hartmann, Rayleigh, Eckert, Prandtl, Magnetic, temperature and Curie temperature numbers, respectively, are

$Da = \frac{\kappa}{L^2}, \ Ha = L\mu_0 \overline{H}_0\sqrt{\frac{\sigma_f}{\mu_f}}, \ Ra = \frac{g\beta_f L^3 \Delta T}{\alpha_f \nu_f}, \ Ec = \frac{\mu_f \alpha_f}{L^2 (\rho c_p)_f \Delta T},$

(2.5)

$Pr = \frac{\nu_f}{\alpha_f}, \ Mn_F = \frac{\mu_0 H_0^2 K' \Delta T L^2}{\mu_f \alpha_f}, \quad \epsilon_1 = \frac{T_1}{\Delta T}, \quad \epsilon_2 = \frac{T_c'}{\Delta T}.$

(2.6)

Pressure terms are eliminated by applying the definition of vorticity $w = \nabla \times {\bf u}$ to the momentum equations Eq.(2.4b)-Eq.(2.4c). The continuity equation is satisfied defining the velocity components in terms of stream function $\psi$ as $u = \partial \psi/ \partial y, v = -\partial \psi/\partial x$. Then, stream function, vorticity and energy equations are deduced as

$\nabla^2 \psi = -w$

(2.7a)

$\begin{eqnarray*} \nabla^2 w & = & \frac{1}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f \mu_{nf}} \left[\frac{\partial w}{\partial t} + \frac{1}{\varepsilon}\bigg( u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y}\big)\right] -\varepsilon Ra \frac{\rho_{nf}\mu_f \beta_{nf}}{\rho_f \mu_{nf} \beta_f} \frac{\partial T}{\partial x} \nonumber \\ &+&\varepsilon Ha^2 \frac{\sigma_{nf}\mu_f}{\sigma_f\mu_{nf}} \bigg[2H_x\frac{\partial H_x}{\partial x}v-\frac{\partial H_x}{\partial x}H_y u-H_x\frac{\partial H_y}{\partial x}u + H_x^2 \frac{\partial v}{\partial x}-H_xH_y \frac{\partial u}{\partial x} \nonumber\\ &-&2H_y\frac{\partial H_y}{\partial y}u - H_y^2\frac{\partial u}{\partial y} +\frac{\partial H_x}{\partial y}H_y v + H_x\frac{\partial H_y}{\partial y}v + H_xH_y\frac{\partial v}{\partial y}\bigg] \nonumber \\ & -& \varepsilon Mn_F \frac{\mu_f}{\mu_{nf}} H \left(\frac{\partial H}{\partial x}\frac{\partial T}{\partial y} - \frac{\partial H}{\partial y}\frac{\partial T}{\partial x}\right) \nonumber \\ & +& \frac{\varepsilon}{Da}w + \frac{c_f}{\sqrt{Da}}\frac{\varepsilon}{Pr}\frac{\rho_{nf}\mu_f}{\rho_f\mu_{nf}}\bigg( v \frac{\partial |{\bf u}|}{\partial x} - u \frac{\partial |{\bf u}|}{\partial y} + |{\bf u}|w\bigg). \end{eqnarray*}$

(2.7b)

$\begin{align} & {{\nabla }^{2}}T = \frac{{{(\rho {{c}_{p}})}_{nf}}{{k}_{f}}}{{{(\rho {{c}_{p}})}_{f}}{{k}_{nf}}}\left(\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y} \right)-H{{a}^{2}}Ec\frac{{{\sigma }_{nf}}{{k}_{f}}}{{{\sigma }_{f}}{{k}_{nf}}}{{({{H}_{x}}v-{{H}_{y}}u)}^{2}} \\ & \ \ \ \ \ \ \ -Ec\frac{{{\mu }_{nf}}{{k}_{f}}}{{{\mu }_{f}}{{k}_{nf}}}\left[ 2{{\left(\frac{\partial u}{\partial x} \right)}^{2}}+2{{\left(\frac{\partial v}{\partial y} \right)}^{2}}+{{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \right)}^{2}} \right] \\ & \ \ \ \ \ \ \ -M{{n}_{F}}Ec\frac{{{k}_{f}}}{{{k}_{nf}}}\left(u\frac{\partial H}{\partial x}+v\frac{\partial H}{\partial y} \right)H({{\epsilon }_{1}}+T), \\ \end{align}$

(2.7c)

where

$\frac{\partial H_x}{\partial x} = H_{xdx} = 2|b|\frac{(y-b)(a-x)}{((x-a)^2 + (y-b)^2)^2}, \quad \frac{\partial H_x}{\partial y} = H_{xdy} = |b|\frac{(x-a)^2-(y-b)^2}{((x-a)^2 + (y-b)^2)^2},$

(2.8a)

$\frac{\partial H_y}{\partial y} = H_{ydx} = 2|b|\frac{(y-b)(x-a)}{((x-a)^2 + (y-b)^2)^2}, \quad \frac{\partial H_y}{\partial x} = H_{ydy} = |b|\frac{(x-a)^2-(y-b)^2}{((x-a)^2 + (y-b)^2)^2},$

(2.8b)

$\frac{\partial H}{\partial x} = H_{dx} = |b|\frac{a-x}{((x-a)^2 + (y-b)^2)^{3/2}}, \quad \frac{\partial H}{\partial y} = H_{dy} = |b|\frac{b-y}{((x-a)^2 + (y-b)^2)^{3/2}}.$

(2.8c)

---

3. Numerical approach

Radial basis functions based pseudo spectral (RBF-PS) method generates all space derivatives in the problem by RBFs. Radial basis functions enable one to use any type of grids, and are attracted by most of the researchers in mesh-free methods. The novel books ^[4], ^[5] involve both theory and applications on RBFs.

RBFs approximate an unknown $\varphi$ ($\psi, T$ or $\omega$) in a diffusion-convection type equation, $\nabla^2 \varphi = {\bf u} \cdot \nabla \varphi$, as

$\varphi_i = \sum\limits_{j = 1}^{N_b+N_i} \alpha_j f_{ij},$

(3.1)

where $N_b$ is the number of boundary nodes, $N_i$ is the number of interior nodes, $f$'s are RBFs depending on radial distance $r = || {\bf x} -{\bf x}_j||$ (${\bf x} = (x, y)$ is the field point and ${\bf x}_j = (x_j, y_j)$ is the collocation point), and $\alpha_j$'s are initially unknown coefficients.

Eq.(3.1) can also be expressed in matrix-vector form as

$\varphi = F \alpha,$

(3.2)

where the matrix $F$ is of size $(N_b+N_i)\times (N_b+N_i)$, and the coefficient vector $\alpha$ is $\{\alpha_1, \ldots, \alpha_{N_b+N_i} \}$. Eq.(3.2) can be rewritten as

$\alpha = F^{-1} \varphi.$

(3.3)

Eq.(3.1) and Eq.(3.3) lead to the first and second order space derivatives of $\varphi$ as

$\frac{\partial \varphi}{\partial x} = \frac{\partial F}{\partial x}\alpha = \frac{\partial F}{\partial x}F^{-1}\varphi, \qquad \frac{\partial \varphi}{\partial y} = \frac{\partial F}{\partial y}\alpha = \frac{\partial F}{\partial y}F^{-1}\varphi,$

(3.4a)

$\frac{\partial^2 \varphi}{\partial x^2} = \frac{\partial}{\partial x}\left( \frac{\partial \varphi}{\partial x} \right) = \frac{\partial^2 F}{\partial x^2}F^{-1}\varphi, \qquad \frac{\partial^2 \varphi}{\partial y^2} = \frac{\partial^2 F}{\partial y^2}F^{-1}\varphi.$

(3.4b)

RBF-PS method for the space derivatives is carried out to simulate the equations (2.7a), (2.7b), (2.7c). The time derivative is managed by Backward-Euler method. Thus, an iterative system with the application of these two methods to the non-dimensional system is built as follows

$D_2 \psi^{n+1} = -w^n$

(3.5a)

$\left( D_2 - C_1\frac{\mathcal{I}}{\Delta t} - C_1 M \right)T^{n+1} = -C_1\frac{T^n}{\Delta t} - C_2 Ha^2 Ec (H_x v - H_y u)^2 \nonumber \\ \quad \quad -C_3 Ec [2(D_x u)^2 + 2(D_y v)^2 + (D_y u + D_x v)^2] \\ \quad \quad - C_7 Mn_F Ec (u H_{dx} + v H_{dy})H(\epsilon_1 + T^{n})$

(3.5b)

$\left(D_2 - \frac{C_4}{Pr}\frac{\mathcal{I}}{\Delta t} - \frac{C_4}{Pr}M - \frac{\mathcal{I}}{Da} \right)w^{n+1} = -\frac{C_4}{Pr}\frac{w^n}{\Delta t} - C_5Ra D_x T^{n+1} \nonumber \\ \quad \quad \hspace{1cm}+ C_6 Ha^2 \big[2H_x H_{xdx} v-H_{xdx} H_y u-H_x H_{ydx} u + H_x^2 (D_x v)-H_x H_y (D_x u) \nonumber \\ \quad \quad \hspace{1cm} - 2H_y H_{ydy} u - H_y^2 (D_y u)+ H_{xdy} H_y v + H_x H_{ydy} v + H_x H_y (D_y v) \big] \nonumber \\ \quad \quad \hspace{1cm} - C_8 Mn_F H \big[H_{dx}(D_y T^{n+1})-H_{dy} (D_x T^{n+1}) \big] \nonumber \\ \quad \quad \hspace{1cm} + \frac{C_9}{\sqrt{Da}Pr} \bigg( v D_x |{\bf u}| - u D_y |{\bf u}| + |{\bf{u}}| w^n \bigg)$

(3.5c)

where $n$ is the number of time iteration, $\mathcal{I}$ is the identity matrix, $C_1, C_2, C_3, C_4$, $C_5, C_6, C_7, C_8, C_9$ are constants seen in Eqs.(2.7c), (2.7a), (2.7b), and $D_2, D_x, D_y, M$ matrices are

$D_2 = \left( \frac{\partial^2 F}{\partial x^2} + \frac{\partial^2 F}{\partial y^2} \right) F^{-1}, \quad D_x = \frac{\partial F}{\partial x}F^{-1}, \quad D_y = \frac{\partial F}{\partial y}F^{-1}, \quad M = [u]_d D_x + [v]_d D_y.$

(3.6)

Initially, $\psi, T$ and $\omega$ are taken as zero except on the boundary. The velocity components are computed after Eq.(3.5a) by $u = u^{n+1} = D_y \psi^{n+1}, \ v = v^{n+1} = -D_x \psi^{n+1}$. Then, no-slip boundary conditions for velocity are imposed. Boundary conditions of vorticity are handled by using the definition of vorticity as $\omega_{bc} = D_x v^{n+1} -D_y u^{n+1}$. Once the vorticity equation is solved, a relaxation parameter $\tau (0 < \tau < 1)$ is used as $\omega^{n+1}\leftarrow \tau\omega^{n+1} + (1-\tau)\omega^n$.

The iterations continue until the criterion

$\sum\limits_{k = 1}^3 \frac{||\varphi_k^{n+1}-\varphi_k^{n}||_{\infty}}{|| \varphi_k^{n+1} ||_{\infty}} < 10^{-5},$

(3.7)

is satisfied.

The average Nusselt number through the left hot wall is computed by

$\overline{Nu} = -\frac{k_{nf}}{k_f} \int_0^1 \frac{\partial T}{\partial x}\ dy.$

(3.8)

---

4. Numerical results and discussion

The multiquadric radial basis function $f = \sqrt{r^2+c^2}$ is used. The shape parameter $c$ is determined by taking the average value of minimum of $r$ distances (different than zero) of any point. In this view, $c$ is considered depending on the node distribution. This determination is done with regard to the node distribution given in Figure 2.

The validation of the method is done computing the average Nusselt number through the heated left wall considering different well known problems. The results of Table 2 and Table 3 are in good agreement with ^[11]. This also validates that the proposed way of finding a suitable shape parameter works well.

Regarding the grid independence search in Table 4, $N_b = 120, \, N_i = 1201$ number of grid points are used in the computations unless otherwise declared.

Figure 3 illustrates the effect of Hartmann number on fluid flow and heat transfer. Primary vortex in streamlines becomes smaller and is pushed from left to right due to the intense magnetic field imposed by magnetic source. Conductive heat transfer reveals in isotherms as $Ha$ increases. The influence of the magnetic source on vorticity contours is more pronounced.

Figure 4 illustrates the effect of Darcy number. Small Darcy number corresponds to small permeability of the porous medium. This means that fluid flows very tough, so the fluid velocity decreases. As is seen in streamlines, primary vortex is smaller at $Da = 10^{-4}$ than the larger Darcy numbers. Also, the effect of magnetic source almost disappears at $Da = 10^{-4}$. Heat transfer in conduction is transferred into convection as Darcy number augments due to the increase in permeability. While vorticity contours are stagnant at the center, two new cells through the left and right walls together with the effect of magnetic source also occur. Additionally, not much alteration is noted in hydrothermal behaviour if $Da\geq 10$ as is seen from $\overline{Nu}$ values at $Da = 10$ and $Da = 100$.

Figure 5 shows the magnetic number variation at large $Mn_F$ numbers. The main vortex in streamlines for the case $Da = 10$ in Figure 4 is directed to the left wall. Isotherms exhibiting convective behaviour starts to perturb near the magnetic source forming a small cell on the left wall. Strong circulation on the middle of the left wall is also noticed in vorticity contours due to the presence of magnetic source. Additionally, $\overline{Nu}$ is $3.24$ for $Mn_F = 100$, $3.23$ for $Mn_F = 10^3$, and $3.00$ for $Mn_F = 10^4$. So, the convective heat transfer is diminished a little bit.

The influence of magnetic number in lower effect of buoyancy-driven force is examined in Figure 6. While there is no significant change at $Mn_F = 100$, secondary cells are formed in streamlines at larger magnetic numbers. Perturbation occurs in isotherms at $Mn_F = 10^4$. Vorticity contours are also divided into new cells. The increase in average Nusselt number emphasizes that the convective heat transfer rises with the augmentation of magnetic number at $Ra = 10^3$.

The change in Rayleigh number is presented in Figure 7 and Figure 8. The increase in Rayleigh number points to the increase in buoyancy-driven force. The circulation in the flow increases as is seen in streamlines. Primary vortex expands. Isotherms demonstrates the dominance of convective heat transfer forming the temperature gradient on the vertical walls. Center-stagnant vorticity contours starts to form new cells through the corners of the right and left wall. The effect of magnetic source is not noted due to the small Darcy number. In Figure 8, the inhibitive effect of smaller permeability on fluid flow and heat transfer diminishes slowly. The increase in temperature gradient in isotherms is noticed for $Ra = 10^4$ and $Ra = 10^5$. At $Ra = 10^6$, isotherms are reduced through the left top corner of the cavity. This may be due to the existence of viscous dissipation.

At a large Rayleigh number, an interesting phenomena occurs in fluid flow and heat transfer. This is depicted in Figure 9. In this case, $N_b = 160, N_i = 2001$ number of grids are performed. It is seen that as the effect of viscous dissipation increases, heat transfer starts to be changed forming almost no isotherms close to the upper left corner (which is a similar case in Figure 8 at $Ra = 10^6$). In other words, there is a change in fluid flow and heat transfer when $Ec$ number is at a critical number. Besides, average Nusselt number also becomes negative at $Ec = 10^{-5}$. A similar phenomenon is also noted in ^[6].

Table 5 indicates the increase in convective heat transfer with the rise in either porosity of the porous medium or concentration of solid nanoparticles. At a large porosity, fluid flows faster than a smaller porosity. Further, $\overline{Nu}$ increases with the augmentation of nanoparticle concentration due to the greater thermal conductivity of nanoparticle Fe$_3$O$_4$ than the base fluid.

---

5. Conclusion

The effects of Brinkman and Forchheimer terms, Kelvin force terms and magnetic source on ferrofluid flow and heat transfer in a porous cavity are numerically studied. RBF-PS method with MQ RBFs is carried out for simulation. For determining the shape parameter in MQ RBF, a way is also proposed. The procedure is validated with different studies. In numerical results, streamlines, isotherms and vorticity contours in different range of dimensionless parameters are presented. The following concluding remarks are obtained

1. Increasing Hartmann number decreases the convective heat transfer due to the retarding effect of Lorentz forces.

2. Decreasing Darcy number causes the fluid to flow slowly due to the smaller permeability. Also, conductive heat transfer is pronounced at lower values of Darcy number.

3. Larger Magnetic number has a dominant effect on heat transfer at smaller Rayleigh number.

4. Buoyancy-driven effect enhances with the increase in Rayleigh number suppressing the effect of magnetic source.

5. The augmentation of porosity of the porous medium and solid volume fraction results in rise of convective heat transfer.

6. The shape parameter $c$ decreases with the increase in number of grid points.

---

Conflict of interest

Author declares no conflict of interest in this paper.

---