1. Introduction

Natural convective heat transfer has seen a great deal of interest in the last decades due to many engineering applications such as insulation of buildings, solar energy collectors, cooling systems for electronic devices, etc. In order to improve the heat transfer characteristics of traditional liquids, nanosized metallic particles are inserted into the liquid. This increases the thermal conductivity of the fluid.

There are a considerable amount of contributions modelling the heat transfer of the nanofluid by using appropriate numerical approaches. The most common numerical approach is the finite volume method ^{[6,8,10,12,15,17,19]}. Jou et al. ^[9] implemented a finite difference solution with the line it-erative method. Fattahi et al. ^[5] and Ashorynejad et al. ^[1] employed the Lattice Boltzmann method to simulate the heat transfer comparing different nanofluids, and considering the presence of magnetic field, respectively. Gumgum and Tezer-Sezgin showed the effeciency of dual reciprocity boundary element method on natural convection flow of nanofluids. A novel approach is seen in Serna et al. ^[16] with the usage of the network simulation method to show the influence of a pulsating flow in a heated lid-driven cavity filled with a nanofluid.

To the best of the author’s knowledge, this is the first application of RBF-PS and DQM to simulate natural convection in an enclosure filled with a nanofluid. The main effect of this process is to be able to use a small number of grid points in the space domain, and to find the solution both in the space and the time domain at once. The average Nusselt number variation in different nanofluids with different solid volume fractions, and in different aspect ratios is presented. Streamlines, isotherms and vorticity contours are also illustrated.

---

2. Problem Setup

The problem geometry is described in Figure 1. No-slip boundary conditions on the walls ($u=v=\psi=0$) are imposed while the left wall is the hot wall and right wall is the cold wall. The top and bottom walls are adiabatic. Aspect ratio of the cavity is $W:H$.

The nanofluid in the enclosure is laminar, incompressible and Newtonian. The nanoparticles are assumed to be in spherical shape, and the nanoparticles and the fluid are assumed to be in thermal equilibrium. Thermophysical properties of the nanofluid are also assumed constant except the density variation in the buoyancy force term treated by the Boussinesq approximation. The radiation effect and viscous dissipation are negligible.

For some nanoparticles copper ($Cu$), silver ($Ag$), aluminium oxide ($Al_2O_3$), silicon oxide ($SiO_2$), titanium oxide ($TiO_2$), copper oxide ($CuO$), thermophysical properties are given in the Table 1.

---

3. Governing Equations

The continuity, momentum and energy equations are ^[2,13,19]

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

(3.1a)

$\frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial y}}=-\frac{1}{{{\rho _{nf}}}}\frac{{\partial P}}{{\partial x}} + \frac{{{\mu _{nf}}}}{{{\rho _{nf}}}}\left({\frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}}} \right)$

(3.1b)

$\frac{{\partial v}}{{\partial t}} + u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial y}}=-\frac{1}{{{\rho _{nf}}}}\frac{{\partial P}}{{\partial y}} + \frac{{{\mu _{nf}}}}{{{\rho _{nf}}}}\left({\frac{{{\partial ^2}v}}{{\partial {x^2}}} + \frac{{{\partial ^2}v}}{{\partial {y^2}}}} \right) + \frac{{{{(\rho \beta)}_{nf}}}}{{{\rho _{nf}}}}g(T-{T_c})$

(3.1c)

$\frac{{\partial T}}{{\partial t}} + u\frac{{\partial T}}{{\partial x}} + v\frac{{\partial T}}{{\partial y}}={\alpha _{nf}}\left({\frac{{{\partial ^2}T}}{{\partial {x^2}}} + \frac{{{\partial ^2}T}}{{\partial {y^2}}}} \right),$

(3.1d)

where

the effective density of the nanofluid is $\rho_{nf}=(1-\chi)\rho_f + \chi \rho_s\, ,$

the effective dynamic viscosity is $\mu_{nf}=\frac{\mu_f}{(1-\chi)^{2.5}} \, ,$

the thermal expansion coefficient of the nanofluid is $\beta_{nf},$

the heat capacity of the nanofluid is $(C_p)_{nf},$

the thermal diffusivity of the nanofluid is $\alpha_{nf}=\frac{k_{nf}}{(\rho C_p)_{nf}},$

the effective thermal conductivity of the nanofluid is $k_{nf}=k_f \frac{k_s + 2k_f-2\chi(k_f-k_s)}{k_s+2k_f+\chi(k_f-k_s)},$

and $\chi$ is the solid volume fraction, subindices $f$ and $s$ refer to fluid and nanosized solid particle, respectively. Also, the nanofluid has the following properties

${(\rho \beta)_{nf}}=(1-\chi){(\rho \beta)_f} + \chi {(\rho \beta)_s}{\mkern 1mu},$

(3.2)

$(\rho C_p)_{nf}=(1-\chi)(\rho C_p)_f + \chi(\rho C_p)_s\, .$

(3.3)

In order to get the non-dimensional form of the governing equations, the non-dimensional parameters are defined as

$x'=\frac{x}{\ell}, \ y'=\frac{y}{\ell}, \ u'=\frac{u\ell}{\alpha_f}, \ v'=\frac{v\ell}{\alpha_f}, \ P=\frac{p\ell^2}{\rho_{nf}\alpha_f^2}, \ T'=\frac{T-T_c}{T_h-T_c},$

(3.4)

in which $\ell$ is the characteristic length.

Substituting these parameters in Eq.(3.4) into the governing equations Eq.(3.1), and then dropping the prime notation, the following dimensionless equations are obtained

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

(3.5a)

$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y}=-\frac{\rho_f}{\rho_{nf}}\frac{\partial P}{\partial x} + Pr\frac{\mu_{nf}\rho_f}{\mu_f\rho_{nf}}\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right)$

(3.5b)

$\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y}=-\frac{\rho_f}{\rho_{nf}}\frac{\partial P}{\partial y} + Pr\frac{\mu_{nf}\rho_f}{\mu_f \rho_{nf}}\left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) + \frac{(\rho \beta)_{nf}}{\rho_{nf}\beta_f} RaPr T \label{momev}$

(3.5c)

$\frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y}=\frac{\alpha_{nf}}{\alpha_f}\left(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2}\right).$

(3.5d)

The definitions of velocity components in terms of stream function $\psi$, $u=\partial \psi/\partial y, \ v=-\partial \psi/\partial x$ (which satisfy the continuity equation) and of vorticity result in $\omega=\nabla \times {\bf u}=\dfrac{\partial v}{\partial x}-\dfrac{\partial u}{\partial y}=-\nabla^2 \psi$. Pressure terms in Eqs.(3.5b)-(3.5c) are eliminated by applying the definition of vorticity to these equations. By this way, vorticity equation is derived. Thus, the dimensionless governing equations in terms of stream function $\psi$, temperature $T$ and vorticity $\omega$ are

$\nabla^2 \psi=-\omega$

(3.6a)

$\nabla^2 T=\frac{\alpha_f}{\alpha_{nf}} \left(\frac{\partial T}{\partial t} + u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} \right)$

(3.6b)

$\nabla^2 \omega=\frac{\rho_{nf}\mu_f}{\mu_{nf}\rho_f Pr} \left(\frac{\partial \omega}{\partial t} + u\frac{\partial \omega}{\partial x} + v\frac{\partial \omega}{\partial y} \right)-\frac{\mu_f (\rho \beta)_{nf}}{\mu_{nf}\rho_f \beta_f} Ra \frac{\partial T}{\partial x},$

(3.6c)

where $Pr=\dfrac{\nu_f}{\alpha_f}$ is the Prandtl number and $Ra=\dfrac{g\beta_f (T_h-T_c)\ell^3}{\alpha_f \nu_f}$ is the Rayleigh number.

---

4. RBF-PS in space and DQM in time

The diffusion-convection type partial differential equation

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

(4.1)

may be approximated by RBFs of the form

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

(4.2)

where $\varphi$ is an unknown ($\psi, T$ or $\omega$), $\alpha_j$'s are initially unknown coefficients, $f$'s are approximating functions formed by RBFs depending on radial distance $r=||{\bf x}-{\bf x}_j||$ in which ${\bf x}=(x, y)$ is the field point and ${\bf x}_j=(x_j, y_j)$ is the collocation point, $N_i$ is the number of interior nodes, and $N_b$ is the number of boundary nodes.

In matrix-vector form, Eq.(4.2) may also be written as

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

(4.3)

The matrix $F$ of size $(N_i+N_b)\times (N_i+N_b)$ is the matrix formed by $f_j$'s columnwise, and $\alpha=\{\alpha_1, \alpha_2, \ldots, \alpha_{N_i+N_b}\}$ is the coefficient vector.

The first and second order space derivatives of $\varphi$ are derived by using $F$ and Eq.(4.3) 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,$

(4.4)

$\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.$

(4.5)

Using Eqs.(4.4)-(4.5), Eq.(4.1) is expressed in matrix-vector form as

$D_2 \varphi=\left([u]_d D_x + [v]_d D_y \right) \varphi$

(4.6)

$\mbox{or equivalently, }\, (D_2-M)\varphi=0 \label{mvf},$

(4.7)

where the matrices are $\displaystyle D_2=\left(\frac{\partial^2 F}{\partial^2 x} + \frac{\partial^2 F}{\partial y^2} \right)F^{-1}, \ D_x=\frac{\partial F}{\partial x}F^{-1}, \ D_y=\frac{\partial F}{\partial y}F^{-1}$ and $M=\left([u]_d D_x + [v]_d D_y \right)$, and the subscript $d$ refers to diagonal.

The Dirichlet type boundary conditions ($\{ d_{bc} \}$) are inserted to the system matrix of Eq.(4.7) as

$\begin{bmatrix} [D_2-M]_{N_i} & [D_2-M]_{N_b} \\ 0 & I \end{bmatrix} \begin{bmatrix} \varphi_i \\ \varphi_b \end{bmatrix}=\begin{bmatrix} 0 \\ d_{bc} \end{bmatrix},$

(4.8)

with the identity matrix $I$ of size $N_b\times N_b$. In case of Neumann type boundary conditions ($n_{bc}$), the given boundary conditions are added to the system matrix of Eq.(4.7) with the help of Eq.(4.4) as

$\begin{bmatrix} [D_2-M]_{N_i} \\ D_n \end{bmatrix} \begin{bmatrix} \varphi_i \\ \varphi_b \end{bmatrix}=\begin{bmatrix} 0 \\ n_{bc} \end{bmatrix},$

(4.9)

in which $D_n$ is either $D_x$ or $D_y$ with $+$ or $-$ sign with respect to normal vector on the boundary.

In order to handle time derivatives, differential quadrature method (DQM) is used. DQM approximates the derivatives of a function at a grid point by a linear summation of all functional values in the whole problem domain.

The time interval is firstly considered as $[0, t_{\max}]$ with a maximum value of time, $t_{\max}$. System matrices of the governing equations are formed by space and time domains as a box, and these matrices are not too large for small values of $t_{\max}$. However, for large values of $t_{\max}$, any system matrix become too large to compute. To overcome this difficulty for a large system, time is divided into equal time subintervals, and each subinterval is also divided into $L$ number of nonuniform grid points.

DQM manages the first order time derivative as

$\frac{\partial \varphi}{\partial t}=\sum_{k=1}^L a_{lk}\varphi_{ijk},$

(4.10)

where $a_{lk}$ is the first order weighting coefficient given explicitly in ^[18] as

$a_{lk}=\displaystyle \frac{M^{(1)}(t_l)}{(t_l-t_k)M^{(1)}(t_k)}, \quad l\neq k, \quad a_{ll}=\displaystyle-\sum_{k=1, k\neq l}^{L} a_{lk},$

(4.11)

where $M^{(1)}(t_l)=\displaystyle \prod_{k=1, k\neq l}^{L}(t_l-t_k),$ and $l=1, 2, \ldots, L$.

The system in each time subinterval is regarded as a block consisting of space and time domain. Aniterative system at each block on the dimensionless governing equations are built as follows

$D_2 \psi^{n+1}_{ijl}=-\omega^n_{ijl},$

(4.12a)

$u^{n+1}_{ijl}=D_y \psi^{n+1} _{ijl}, \quad v^{n+1}_{ijl}=-D_x \psi^{n+1} _{ijl},$

(4.12b)

$\left(D_2-\dfrac{\alpha_{f}}{\alpha_{nf}}(At + M) \right) T^{n+1}_{ijl}=0,$

(4.12c)

$\left(D_2-\dfrac{\rho_{nf}\mu_f}{Pr\rho_f\mu_{nf}}(At+M) \right) w^{n+1}_{ijl}=-\dfrac{\mu_f (\rho \beta)_{nf}}{\mu_{nf}\rho_f\beta_f} Ra D_x T^{n+1}_{ijl},$

(4.12d)

where $M=[u]^{m+1}_d D_x + [v]^{m+1}_d D_y$, $At$ is the matrix of size $L \times L$ formed by weighting coefficient in Eq.(4.11), and $n$ shows the iteration level.

In the first block, $\psi, T$ and $\omega$ are attained as zero at $n=0$ except the known boundary at $n=0$. Once this iteration is completed at a block, the next block starts with the initial values taken as the results of the previous block. This enables one to reach aimed value of $t_{\max}$.

The unknown vorticity boundary conditions are handled by using the definition of vorticity as

$\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}=D_x v^{n+1}-D_y u^{n+1}.$

(4.13)

The resulting systems of equations in the form $Ax=b$ are solved by Gaussian elimination with partial pivoting, and QR factorization which is for overdetermined system of temperature equation because of adiabatic boundary conditions.

The average Nusselt number through the heated left wall is defined by $\displaystyle \overline{Nu}=\int_0^1-\dfrac{\partial T}{\partial x} dy,$ and computed by Clenshaw Curtis quadrature due to the usage of Chebyshev non-uniform grid distribution.

Multiquadric (MQ) $f=\sqrt{r^2 + c^2}$ and inverse multiquadric (IMQ) RBFs $f=1/\sqrt{r^2 + c^2}$ are employed in this study. The shape parameter $c$ controls the shape of the basis functions. As $c$ gets larger, the shape becomes flat and the matrix becomes more ill-conditioned. MQ collocation matrices are conditionally positive definite ^[14], and the exponential convergence of the error of MQ approximation have been demonstrated by Madych et al. ^[11]. IMQs are strictly positive definite.

In order to determine a suitable shape parameter depending on the problem parameters, an initial interval for the shape parameter where the iterative system converges (or catches the expected behaviour) is determined. This interval is divided into equal $c$ values.

The shape parameter value $c$ providing $\psi_{\max}$ zero or closest to zero is chosen. If more than one $c$ values give $\psi_{max}=0$, then the first $c$ value giving $0$ is taken.

---

5. Numerical Results and Discussion

The implementations and computations are done in Matlab. In implementation of space derivatives Eqs.(4.4)-(4.5), right back slash operator, which makes use of Gaussian elimination with partial pivoting, is used instead of taking the inverse of $F$ directly.

The base fluid is taken as pure water with Prandtl number 6.2.

A relaxation parameter $0 < \gamma < 1$ is used as $\varphi^{n+1} \leftarrow \gamma \varphi^{n+1} + (1-\gamma)\varphi^n$ once eqs.(4.12a), (4.12c) or (4.12d) are performed. In particular, $\gamma=0.1$ is employed for $Ra=10^3, \, Ra=10^4$ and $Ra=10^5$ once eqs.(4.12a), (4.12c) and (4.12d) are solved, and $\gamma=0.01$ is carried out for $Ra=10^6$ once eq.(4.12d) is solved.

$t_{\max}$ is fixed at $20$ for $AR=1$, at $21$ for $AR\neq 1$, and each block has an up time level $L=5$ for $AR=1$, $L=3$ for $AR\neq 1$.

Table 2 describes the well agreement of the proposed scheme with the benchmark problem in ^[3] in which $41\times 41$ number of grid points are used. Both MQ RBF and IMQ RBF are carried out, and $c$ values are larger in IMQ than in MQ. In both cases, $c$ values decrease as $Ra$ increases.

When there is no nanoparticle inside of the water, the average Nusselt number is $\overline{Nu}=4.7062$ in $Ra=10^5$ utilizing MQ RBF. As is seen from Table 3, the presence of any type of nanoparticle enhances heat transfer in the cavity. Also, the increase in solid volume fraction is resulted with the increase in $\overline{Nu}$. The improved heat transfer is pronounced better by the insertion of $Cu$ nanoparticles than the other type of nanoparticles.

In Figure 2, the variation in $Ra$ is shown in $Cu$-water nanofluid with $\chi=0.2$. Streamlines in the center of the cavity expands with the increase in $Ra$. Also, the strong temperature gradient in isotherms through the left and right walls is noticed due to the increase in buoyancy force. The centered vorticity contour is divided into two cells, and these two cells are shrunk through the left and right walls. At a large Rayleigh number ($Ra=10^6$), the circulation in the main centered cell of streamlines in case of $\chi=0$ (Figure 3) is not much different than the case of $\chi=0.2$ (Figure 2). This indicates that laminar flow regime is preserved.

In Table 4, fluid velocity and $\overline{Nu}$ values are presented in different concentration of $Cu$-water nanofluid. The values $N_b, N_i$ and $c$ are also the same for $\chi=0.02$ and $\chi=0.2$ with $\chi=0$. In each concentration, the increase in fluid velocity as $Ra$ increases is noted. Fluid velocity does not change from $\chi=0$ to $\chi=0.02$. However, the remarkable change occurs in $\chi=0.2$. This also points out that fluid flows faster in highly concentrated $Cu$-water nanofluid than pure water. For small $Ra$ values ($Ra=10^3, \, 10^4$), the change in $\chi$ has no much effect on fluid velocity. For $Ra=10^5$ and $Ra=10^6$, the increase in $|\psi |_{\max}$ and $u_{\max}$ are seen in presence of nanoparticle in solid volume fraction $\chi=0.2$. In each $Ra$, the existence of nanoparticle inside water causes $\overline{Nu}$ to increase, and so the convective heat transfer is pronounced.

The increase in fluid velocity with the increase in $\chi$ is also illustrated in Figure 4 with $u$ and $v$ velocity profiles at mid sections of the unit square cavity. The fluid moves at higher velocities close to the boundaries than the center of the cavity as can be seen in Figure 4(a)-(b). In each cases, absolute value of velocities is the largest at the highest concentration of nanofluids $\chi=0.2$.

In Figures 5 and 6, different aspect ratios of the Cu-water-filled cavity are taken into account performing IMQ RBF. Similar behaviour as in $Ra=10^5$ in Figure 2 occurs. Table 5 indicates the average Nusselt number values through the heated wall. On one hand, convective heat transfer decreases inside the cavity when the cavity sizes increase. On the other hand, $\overline{Nu}$ values in case of $0 < AR < 1$ are greater than the case of $AR>1$ which may be a result of the effect of longer heated wall on the fluid flow and heat transfer inside the cavity.

---

6. Conclusion

In this study, a numerical investigation on natural convection in a cavity filled with a nanofluid is presented. The space derivatives in the dimensionless equations are discretised by RBF-PS and the time derivatives are handled by DQM. The nature of this approach provides one to obtain the results with small number of grid points. Instead of solving the system of time and space as an entire large system, time is divided into subintervals and the system is observed block by block. Also, the solution at each block is obtained at once. Convective heat transfer increase with the insertion of nanoparticles into the water. Further, the decrease in the convective heat transfer is noticed with the increase in Cu-water-filled cavity size.

---

Conflict of Interest

I have no financial or other relationship with other people or organizations in my study.

---