Numerical Study of Natural Convection Heat Transfer in a Porous Annulus Filled with a Cu-Nanofluid

Natural convection heat transfer in a porous annulus filled with a Cu nanofluid has been investigated numerically. The Darcy–Brinkman and the energy transport equations are employed to describe the nanofluid motion and the heat transfer in the porous medium. Numerical results including the isotherms, streamlines, and heat transfer rate are obtained under the following parameters: Brownian motion, Rayleigh number (103–105), Darcy number (10−4–10−2), nanoparticle volume fraction (0.01–0.09), nanoparticle diameter (10–90 nm), porosity (0.1–0.9), and radius ratio (1.1–10). Results show that Brownian motion should be considered. The nanoparticle volume fraction has a positive effect on the heat transfer rate, especially with high Rayleigh number and Darcy number, while the nanoparticle diameter has an inverse influence. The heat transfer rate is enhanced with the increase of porosity. The radius ratio has a significant influence on the isotherms, streamlines, and heat transfer rate, and the rate is greatly enhanced with the increase of radius ratio.


Introduction
Fluid flow and heat transfer due to natural convection in porous annulus is one of the most considerable research issues due to its wide applications in science and engineering [1][2][3], such as thermal insulators, chemical catalytic convectors, thermal storage systems, geothermal energy utilization, electronic cooling, and nuclear reactor systems. Therefore, natural convection in porous annulus has been extensively investigated during the past decades. Caltagirone [4] was the first to study natural convection in a saturated porous medium bounded by two concentric, horizontal, isothermal cylinders experimentally and numerically. The author used the Christiansen effect in order to visualize a fluctuating three-dimensional thermal field for Rayleigh number exceeding some critical value. Rao et al. [5,6] performed steady and transient investigations on natural convection in a horizontal porous annulus heated from the inner face using Galerkin method. The effects of Rayleigh number and Darcy number on heat transfer characteristics were studied. In addition, the bifurcation point was obtained numerically, which compared very well with that from experimental observation.
Himasekhar [7] examined the two-dimensional bifurcation phenomena in thermal convection in horizontal, concentric annuli containing saturated porous media. The fluid motion is described by the Darcy-Oberbeck-Boussinesq equations, which were solved using regular perturbation expansion. The flow structure was obtained under different parameters, such as the radius ratio and Rayleigh-Darcy number. A parametric study was performed by Leong et al. [8] to investigate the effects of Rayleigh number, Darcy number, porous sleeve thickness, and relative thermal conductivity on heat transfer characteristics. Braga et al. [9] presented numerical computations for laminar and turbulent natural convection within a horizontal cylindrical annulus filled with a fluid saturated porous medium. Computations covered the range 25 < Ra m < 500 and 3.2 × 10 −4 > Da > 3.2 × 10 −6 and made use of the finite volume method. Khanafer et al. [10] carried out a numerical simulation in order to examine the parametric effects of Rayleigh number and radius ratio on the role played by natural convection heat transfer in the porous annuli. The model was governed by Darcy-Oberbeck-Boussinesq equations and solved using the Galerkin method. In order to investigate the buoyancy-induced flow as affected by the presence of the porous layer, Alloui and Vasseur [11] studied natural convection in a horizontal annular porous layer filled with a binary fluid under the influence of the Soret effect using the Darcy model with the Boussinesq approximation. Numerical solutions of the full governing equations are obtained for a wide range of the governing parameters, such as Rayleigh number, Lewis number, buoyancy ratio, radius ratio of the cavity, and normalized porosity.
Belabid and Cheddadi [12] solved natural convection heat transfer within a twodimensional horizontal annulus filled with a saturated porous medium using ADI (Alternating Direction Implicit) finite difference method. This work placed emphasis on the mesh effect on the determination of the bifurcation point between monocellular and bicellular flows for different values of the aspect ratio. In a recent work, Belabid and Allali [13,14] studied the effects of a periodic gravitational and temperature modulation on the convective instability in a horizontal porous annulus. Results showed that the convective instability is influenced by the amplitude and the frequency of the modulation. Rostami et al. [15] provided a review of recent natural convection studies, including experimental and numerical studies. The effects of the parameters, such as nanoparticle addition, magnetic fields, and porous medium on the natural convection were examined.
The traditional fluids in engineering, such as water and mineral oils, have a primary limitation in the enhancement of heat transfer due to a rather low thermal conductivity. The term nanofluids, which was first put forward by Choi [16], is used to describe the mixture of nanoparticles and base fluid. Due to the relatively higher thermal conductivities, nanofluids are considered as an effective approach to meet some challenges associated with the traditional fluids [17]. In the last few decades, studies on the natural convection of nanofluids in enclosure were conducted by a number of researchers [18]. For instance, Jou and Tzeng [19] investigated the numerically natural convection heat transfer enhancement utilizing nanofluids in a two-dimensional enclosure. Results showed that increasing the buoyancy parameter and volume fraction of nanofluids caused an increase in the average heat transfer coefficient. Ghasemi et al. [20] simulated the natural convection heat transfer in an inclined enclosure filled with a CuO-water nanofluid. The effects of pertinent parameters such as Rayleigh number, inclination angle, and solid volume fraction on the heat transfer characteristics were studied. The results indicated that the heat transfer rate is maximized at a specific inclination angle depending on Rayleigh number and solid volume fraction.
In a related work, Abu-Nada and Oztop [21] found that the effect of nanoparticles concentration on Nusselt number was more pronounced at low volume fraction than at high volume fraction and the inclination angle could be a control parameter for nanofluid filled enclosure. Soleimani et al. [22] studied the natural convection heat transfer in a semiannulus enclosure filled with a Cu-water nanofluid using the Control Volume based Finite Element Method. The numerical investigation was carried out for different governing parameters, such as Rayleigh number, nanoparticle volume fraction, and angle of turn for the enclosure. The results revealed that there was an optimal angle of turn in which the average Nusselt number was maximum for each Rayleigh number. Seyyedi et al. [23] simulated the natural convection heat transfer of Cu-water nanofluid in an annulus enclosure using the Control Volume-based Finite Element Method. The Maxwell-Garnetts and Brinkman models were employed to estimate the effect of thermal conductivity and viscosity of nanofluid. The results showed the effects of the governing parameters on the local Nusselt number, average Nusselt number, streamlines, and isotherms. Boualit et al. [24] and Liao [25] studied respectively natural convection heat transfer of Cu-and Nanomaterials 2021, 11, 990 3 of 24 Al 2 O 3 -water nanofluids in a square enclosure under the horizontal temperature gradient. Both the flow structure and the corresponding heat transfer characteristics at different Rayleigh numbers and nanoparticle volume fractions were obtained. Wang et al. [26] investigated numerically the natural convection in a partially heated enclosure filled with Al 2 O 3 nanofluids. The results indicated that at low Rayleigh numbers, the heat transfer performance increased with nanoparticle volume fraction, while at high Rayleigh numbers, there existed an optimal volume fraction at which the heat transfer performance had a peak. In a recent work, Mi et al. [27] examined the effects of graphene nano-sheets (GNs) nanoparticles by comparing the thermal conductivity of graphene nano-sheets (GNs)/ethylene glycol (EG) nanofluid with EG thermal conductivity. Results showed that the presence of nanoparticles improved the thermal conductivity, and with increasing temperature, the effect of adding GNs was strengthened.
Although several numerical and experimental studies on the natural convection heat transfer were published, most of them concentrated on traditional fluids in cavities, and only a few of them consider a nanofluid in a porous annulus [28][29][30]. In the present work, steady natural convection heat transfer in a porous annulus filled with a Cu-nanofluid has been investigated, and the governing equations, including the Darcy-Brinkman equation, were solved using the Galerkin method. This paper presented a systematical examination on the effects of Brownian motion, solid volume fraction, nanoparticle diameter, Rayleigh number, Darcy number, porosity on the flow pattern, temperature distribution, and heat transfer characteristics. To the best of our knowledge, no study on this problem has been considered before, and accordingly, the current paper will address this topic.

Physical Description
We consider porous annulus filled with a Cu-water nanofluid between a horizontal inner and outer cylinder of radius r i and r o , respectively, as shown in Figure 1. The inner and outer cylinders are kept at uniform high temperature T i and low temperature T o , respectively. It is taken into consideration that the flow is two-dimensional, steady, and laminar due to the low velocity. The porous medium is considered as isotropic, homogeneous, and filled with a nanofluid, which is thermal equilibrium with the solid matrix, and the Darcy-Brinkman equation without inertia item is adopted. For the nanofluid, the effect of Brownian motion [31][32][33] is considered. The thermophysical properties of the base water, copper nanoparticles, and solid structure of the porous medium used in this study are presented in Tables 1 and 2 [34,35].  [24] and Liao [25] studied respectively natural convection heat transfer of Cu-and Al2O3water nanofluids in a square enclosure under the horizontal temperature gradient. Both the flow structure and the corresponding heat transfer characteristics at different Rayleigh numbers and nanoparticle volume fractions were obtained. Wang et al. [26] investigated numerically the natural convection in a partially heated enclosure filled with Al2O3 nanofluids. The results indicated that at low Rayleigh numbers, the heat transfer performance increased with nanoparticle volume fraction, while at high Rayleigh numbers, there existed an optimal volume fraction at which the heat transfer performance had a peak. In a recent work, Mi et al. [27] examined the effects of graphene nano-sheets (GNs) nanoparticles by comparing the thermal conductivity of graphene nano-sheets (GNs)/ethylene glycol (EG) nanofluid with EG thermal conductivity. Results showed that the presence of nanoparticles improved the thermal conductivity, and with increasing temperature, the effect of adding GNs was strengthened.
Although several numerical and experimental studies on the natural convection heat transfer were published, most of them concentrated on traditional fluids in cavities, and only a few of them consider a nanofluid in a porous annulus [28][29][30]. In the present work, steady natural convection heat transfer in a porous annulus filled with a Cu-nanofluid has been investigated, and the governing equations, including the Darcy-Brinkman equation, were solved using the Galerkin method. This paper presented a systematical examination on the effects of Brownian motion, solid volume fraction, nanoparticle diameter, Rayleigh number, Darcy number, porosity on the flow pattern, temperature distribution, and heat transfer characteristics. To the best of our knowledge, no study on this problem has been considered before, and accordingly, the current paper will address this topic.

Physical Description
We consider porous annulus filled with a Cu-water nanofluid between a horizontal inner and outer cylinder of radius ri and ro, respectively, as shown in Figure 1. The inner and outer cylinders are kept at uniform high temperature Ti and low temperature To, respectively. It is taken into consideration that the flow is two-dimensional, steady, and laminar due to the low velocity. The porous medium is considered as isotropic, homogeneous, and filled with a nanofluid, which is thermal equilibrium with the solid matrix, and the Darcy-Brinkman equation without inertia item is adopted. For the nanofluid, the effect of Brownian motion [31][32][33] is considered. The thermophysical properties of the base water, copper nanoparticles, and solid structure of the porous medium used in this study are presented in Tables 1 and 2 [34,35].

Governing Equations and Boundary Conditions
By introducing the Boussinesq approximation, the governing equations for the heat transfer and fluid flow can be written as follows: ∂u ∂x where (x, y) are the Cartesian coordinates of the geometry, (u, v) are the velocity components, T and p are the temperature and pressure, respectively, ρ, (ρc p ), ν, β and α denote the density, heat capacitance, kinematic viscosity, thermal expansion coefficient, and thermal diffusion, respectively, k is the thermal conductivity, ε is the porosity, and K is the medium permeability. The subscripts nf and mnf designate nanofluid and porous medium filled with nanofluid. The effective thermal conductivity of the porous medium filled with nanofluid can be as modeled as: Using the following dimensional variables, The governing Equations (1)-(4) reduce to a dimensionless form:

∂U ∂X
where (X, Y) are the dimensionless Cartesian coordinates of the geometry, (U, V) are the dimensionless velocity components, and θ and P are the temperature and pressure, respectively. Ra, Pr, and Da denote respectively the Raleigh number, Prandtl number, and Darcy number. The boundary conditions for this problem are as follows: On the outer cylinder surface: On the inner cylinder surface: On the outer and inner cylinder surfaces: The local and overall heat transfer rate along the inner cylinder surface are estimated using the local Nusselt number and average Nusselt number, respectively: Nu loc (14) where S is the non-dimensional length along the inner cylinder surface.

Numerical Procedure
The dimensionless governing Equations (6)-(9) together with the boundary conditions (10) and (11) have been solved numerically using the commercial software tool, which is known as COMSOL Multiphysics (Version 5.5, COMSOL Inc., Stockholm, Sweden). The software employs the Galerkin finite element method, which enforces the orthogonality of residuals to all basis functions in a basis. In Galerkin formulation, weighting functions are chosen to become identical to basis functions [36]. In this paper, we have employed a segregated and parallel direct (Pardiso) solver to solve those equations. As convergence criteria, 10 −6 has been chosen for all dependent variables.

Grid Generation and Independence Test
In the finite element method, grid generation is the technique to discretize the computational domain into subdomains. In this study, we have adopted unstructured triangular elements in the interior and structured quadrilateral elements on the boundary. Figure 2 is the grid generation of the structure with a legend of quality measure. A quality of 1 represents the best possible grid quality, and a quality of 0 represents the worst possible grid quality. For purpose of ensuring the grid independence of the numerical solution, different grid levels in the Comsol Multiphysics are examined. As shown in Table 3, the average Nusselt number of on the inner cylinder at different grids is presented for Cu nanofluid when the nanoparticle volume fraction (φ) is 0.5, nanoparticle diameter (d sp ) is 50 nm, Rayleigh number (Ra) is 10 5 , porosity (ε) is 0.5, and Darcy number (Da) is 10 −2 . The difference between normal and extremely fine is within 1.35%. By comprehensive considering the calculation accuracy and cost, the extra fine level is chosen in this study. different grid levels in the Comsol Multiphysics are examined. As shown in Table 3, the average Nusselt number of on the inner cylinder at different grids is presented for Cu nanofluid when the nanoparticle volume fraction (ϕ) is 0.5, nanoparticle diameter (dsp) is 50 nm, Rayleigh number (Ra) is 10 5 , porosity (ε) is 0.5, and Darcy number (Da) is 10 −2 . The difference between normal and extremely fine is within 1.35%. By comprehensive considering the calculation accuracy and cost, the extra fine level is chosen in this study.

Code Validation
Due to the lack of experimental data for conjugate heat transfer of nanofluid in an annulus filled with porous medium, we have compared our results with the numerical results of Abhishek Kumar Singh and Tanmay Basak et al. [37] for a square cavity filled with base fluid (ϕ = 0). The comparisons are presented in Table 4, when ε = 0.4, and the deviations are 0%, 0%, and 0.84% with Ra = 10 3 , 10 4 , and 10 5 , respectively. When ε = 0.9, the deviations are

Code Validation
Due to the lack of experimental data for conjugate heat transfer of nanofluid in an annulus filled with porous medium, we have compared our results with the numerical results of Abhishek Kumar Singh and Tanmay Basak et al. [37] for a square cavity filled with base fluid (φ = 0). The comparisons are presented in Table 4, when ε = 0.4, and the deviations are 0%, 0%, and 0.84% with Ra = 10 3 , 10 4 , and 10 5 , respectively. When ε = 0.9, the deviations are 0.5%, 1.2%, and 1.3% with Ra = 10 3 , 10 4 , and 10 5 , respectively. It is clear that the current results are in good agreement with the earlier work, and the maximum deviation is 1.3%. In addition, it can be seen that the deviation is increased with the increase of Rayleigh number and porosity. The validation work has enhanced the confidence in the numerical solution of the current study.

Effects of Brownian Motion
In this section, we have investigated the effect of Brownian motion on the heat transfer characteristics using two applied models in Table 2. One ignores the Brownian motion of the nanoparticles. The other takes Brownian motion into consideration and the modified effective thermal conductivity and effective dynamic viscosity are adopted. Figure 3 presents the effect of Brownian motion on the average Nusselt number along the inner wall under different parameters. From Figure 3a, the average Nusselt number increases generally as Brownian motion is considered. With the increase of nanoparticle volume fraction, the influence of Brownian motion becomes more noticeable. In addition, with the decrease of nanoparticle diameter, the influence of Brownian motion becomes more remarkable. For instance, when Ra = 5 × 10 3 , Da = 10 −2 , d sp = 90 nm, Nu avg with Brownian motion increased by 0.21% compared with that without Brownian motion at φ = 0.01, while the growth rate is 2.7% at φ = 0.09. When Ra = 5 × 10 3 , Da = 10 −2 , d sp = 10 nm, Nu avg with Brownian motion increased by 1.98% compared with that without Brownian motion at φ = 0.01, while the growth rate is 23.76% at φ = 0.09. Comparing Figure 3a,b and Figure 3b with Figure 3c, the effect of Brownian motion becomes more remarkable with the increase of Rayleigh number and the decrease of Darcy number. For example, when Da = 10 −2 , φ = 0.09, d sp = 10 nm, Nu avg with Brownian motion increased by 23.76% compared with that without Brownian motion at Ra = 5 × 10 3 , while the growth rate is 35.92% at Ra = 5 × 10 4 . When Ra = 5 × 10 4 , φ = 0.09, d sp = 10 nm, Nu avg with Brownian motion increased by 35.92% compared with that without Brownian motion at Da = 10 −2 , while the growth rate is 46.29% at Da = 10 −3 . Therefore, the effect of Brownian motion on the natural convection heat transfer of the nanofluid should be considered. Furthermore, the positive effect of Brownian motion on the overall heat transfer rate is different at different parameters. Figure 4 shows the isotherms and streamlines for different nanoparticle volume fraction and Rayleigh number at Da = 10 −2 , d sp = 50 nm, ε = 0.5, and RR = 2. The color scales on the left represent the dimensionless temperature and those on the right represent the dimensionless velocity, as in other sections of the article. From Figure 4, for Ra = 10 3 , the isotherms have a uniform distribution; this is because the buoyancy force is weak compared with the viscous force, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The effect of volume fraction on the isotherms is weak. For Ra = 10 4 , a slight thermal disturbance appeared, which indicates that the flow is enhanced and the transition from conduction to natural convection takes place. In this case, the effect of volume fraction becomes more important. For Ra = 10 5 , the isotherms are almost horizontally distributed, especially when φ = 0.9, which means that the natural convection heat transfer turns out to be more significant and the effect of volume fraction is more pronounced. With the increase of Rayleigh number, the streamlines become denser near the walls and the cell becomes bigger and has a tendency to move upward due to the enhanced buoyance force. In addition, for Ra = 10 3 and Ra = 10 4 , the effect of volume fraction on the streamlines is very weak, while for Ra = 10 5 , the effect is more pronounced. Figure 5 displays the effect of volume fraction on the overall heat transfer rate along the inner wall at different Rayleigh numbers. The figure shows that an increase in volume fraction leads to heat transfer enhancement for all considered Rayleigh numbers, and the effect of volume fraction is more pronounced when the Rayleigh number is high.   Figure 4 shows the isotherms and streamlines for different nanoparticle volume fraction and Rayleigh number at Da = 10 −2 , dsp = 50 nm, ε = 0.5, and RR = 2. The color scales on the left represent the dimensionless temperature and those on the right represent the dimensionless velocity, as in other sections of the article. From Figure 4, for Ra = 10 3 , the isotherms have a uniform distribution; this is because the buoyancy force is weak compared with the viscous force, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The effect of volume fraction on the isotherms is weak. For Ra = 10 4 , a slight thermal disturbance appeared, which indicates that the flow is enhanced and the transition from conduction to natural convection takes place. In this case, the effect of volume fraction becomes more important. For Ra = 10 5 , the isotherms are almost horizontally distributed, especially when ϕ = 0.9, which means that the natural convection heat transfer turns out to be more significant and the effect of volume fraction is more pronounced. With the increase of Rayleigh number, the streamlines become denser near the walls and the cell becomes bigger and has a tendency to move upward due to the enhanced buoyance force. In addition, for Ra = 10 3 and Ra = 10 4 , the effect of volume fraction on the streamlines is very weak, while for Ra = 10 5 , the effect is more pronounced.  Figure 5 displays the effect of volume fraction on the overall heat transfer rate along the inner wall at different Rayleigh numbers. The figure shows that an increase in volume fraction leads to heat transfer enhancement for all considered Rayleigh numbers, and the effect of volume fraction is more pronounced when the Rayleigh number is high.   Figure 6 illustrates the isotherms and streamlines for different nanoparticle volume fraction and Darcy number at dsp = 50 nm, Ra = 10 5 , ε = 0.5, and RR = 2. From Figure 6, for Da = 10 −4 , the isotherms have a uniform distribution due to the low permeability, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The effect of volume fraction on the isotherms is slight. For Da = 10 −3 , due to the enhanced permeability, the flow in the annulus is strengthened, and the transition from conduction to natural convection takes place. In this case, the effect of volume fraction becomes more important. For Da = 10 −2 , the isotherms are almost horizontally distributed, especially when ϕ = 0.9, which means the natural convection heat transfer plays a significant role and the effect of volume fraction is more pronounced. With the increase of Darcy number, the streamlines become denser near the walls and the cells become bigger and have a tendency to move upward due to the enhanced flow. In addition, for Da = 10 −4 and Da = 10 −3 , the effect of volume fraction on the streamlines is very weak, while for Da = 10 −2 , the effect is more pronounced. Figure 7 displays the effect of volume fraction on the overall heat transfer rate along the inner wall at different Darcy numbers. The figure shows that an increase in volume fraction leads to heat transfer enhancement for all considered Darcy numbers and the effect of volume fraction is increased with the increase of Darcy number. Figure 8 presents the evolution of the local Nusselt number along the inner wall for different nanoparticle volume fractions at Ra = 10 5 , Da = 10 −2 , dsp = 50 nm, ε = 0.5, and RR = 2. The increase of Nuloc in the whole region means that the local heat transfer rate is enhanced. It can be found that the maximum Nuloc occurred at γ = 180°, which means the natural convection heat transfer is more intense in the bottom half of the inner wall. In addition, the heat transfer is enhanced with the increase of volume fraction.   Figure 6, for Da = 10 −4 , the isotherms have a uniform distribution due to the low permeability, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The effect of volume fraction on the isotherms is slight. For Da = 10 −3 , due to the enhanced permeability, the flow in the annulus is strengthened, and the transition from conduction to natural convection takes place. In this case, the effect of volume fraction becomes more important. For Da = 10 −2 , the isotherms are almost horizontally distributed, especially when φ = 0.9, which means the natural convection heat transfer plays a significant role and the effect of volume fraction is more pronounced. With the increase of Darcy number, the streamlines become denser near the walls and the cells become bigger and have a tendency to move upward due to the enhanced flow. In addition, for Da = 10 −4 and Da = 10 −3 , the effect of volume fraction on the streamlines is very weak, while for Da = 10 −2 , the effect is more pronounced. Figure 7 displays the effect of volume fraction on the overall heat transfer rate along the inner wall at different Darcy numbers. The figure shows that an increase in volume fraction leads to heat transfer enhancement for all considered Darcy numbers and the effect of volume fraction is increased with the increase of Darcy number. Figure 8 presents the evolution of the local Nusselt number along the inner wall for different nanoparticle volume fractions at Ra = 10 5 , Da = 10 −2 , d sp = 50 nm, ε = 0.5, and RR = 2. The increase of Nu loc in the whole region means that the local heat transfer rate is enhanced. It can be found that the maximum Nu loc occurred at γ = 180 • , which means the natural convection heat transfer is more intense in the bottom half of the inner wall. In addition, the heat transfer is enhanced with the increase of volume fraction. Nanomaterials 2021, 11, x FOR PEER REVIEW 11 of 24 ϕ = 0.01 ϕ = 0.09    Figure 9 shows the isotherms and streamlines for different nanoparticle diameters and Rayleigh numbers at Da = 10 −2 , ϕ = 0.05, ε = 0.5 and RR = 2. From Figure 9, for Ra = 10 3 , the isotherms have a uniform distribution due to the weak buoyancy force, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The nanoparticle diameter has a weak effect on the isotherms. For Ra = 10 4 , the natural convection heat transfer is strengthened due to the enhanced buoyancy force. The isotherms near the top half of the inner wall are disturbed, and the change of the isotherms distribution is more remarkable at dsp = 10 compared with dsp = 90. For Ra = 5 × 10 4 , the isotherms are almost horizontally distributed, especially when dsp = 10, which means that the natural convection dominates the heat transfer and the effect of the nanoparticle diameter is more pronounced. In addition, for Ra = 10 3 and Ra = 10 4 , the effect of nanoparticle diameter on the streamlines is very weak, while for Ra = 10 5 , the effect is more pronounced, and the cell becomes bigger and has a tendency to move upward. Figure 10 displays the effect of nanoparticle diameter on the overall heat transfer rate along the inner wall at different Rayleigh numbers. The figure shows that an increase in the nanoparticle diameter leads to reduced heat transfer for all considered Rayleigh numbers. For Ra = 10 3 , the effect of nanoparticle diameter is less pronounced. For Ra = 5 × 10 4 , the effect of nanoparticle diameter is remarkable, especially when it is at a low value. For instance, Nuavg decreased by 11.79% from dsp =   Figure 9 shows the isotherms and streamlines for different nanoparticle diameters and Rayleigh numbers at Da = 10 −2 , ϕ = 0.05, ε = 0.5 and RR = 2. From Figure 9, for Ra = 10 3 , the isotherms have a uniform distribution due to the weak buoyancy force, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The nanoparticle diameter has a weak effect on the isotherms. For Ra = 10 4 , the natural convection heat transfer is strengthened due to the enhanced buoyancy force. The isotherms near the top half of the inner wall are disturbed, and the change of the isotherms distribution is more remarkable at dsp = 10 compared with dsp = 90. For Ra = 5 × 10 4 , the isotherms are almost horizontally distributed, especially when dsp = 10, which means that the natural convection dominates the heat transfer and the effect of the nanoparticle diameter is more pronounced. In addition, for Ra = 10 3 and Ra = 10 4 , the effect of nanoparticle diameter on the streamlines is very weak, while for Ra = 10 5 , the effect is more pronounced, and the cell becomes bigger and has a tendency to move upward. Figure 10 displays the effect of nanoparticle diameter on the overall heat transfer rate along the inner wall at different Rayleigh numbers. The figure shows that an increase in the nanoparticle diameter leads to reduced heat transfer for all considered Rayleigh numbers. For Ra = 10 3 , the effect of nanoparticle diameter is less pronounced. For Ra = 5 × 10 4 , the effect of nanoparticle diameter is remarkable, especially when it is at a low value. For instance, Nuavg decreased by 11.79% from dsp =  Figure 9 shows the isotherms and streamlines for different nanoparticle diameters and Rayleigh numbers at Da = 10 −2 , φ = 0.05, ε = 0.5 and RR = 2. From Figure 9, for Ra = 10 3 , the isotherms have a uniform distribution due to the weak buoyancy force, and it indicates that the heat transfer in the annulus is dominated by thermal conduction. The nanoparticle diameter has a weak effect on the isotherms. For Ra = 10 4 , the natural convection heat transfer is strengthened due to the enhanced buoyancy force. The isotherms near the top half of the inner wall are disturbed, and the change of the isotherms distribution is more remarkable at d sp = 10 compared with d sp = 90. For Ra = 5 × 10 4 , the isotherms are almost horizontally distributed, especially when d sp = 10, which means that the natural convection dominates the heat transfer and the effect of the nanoparticle diameter is more pronounced. In addition, for Ra = 10 3 and Ra = 10 4 , the effect of nanoparticle diameter on the streamlines is very weak, while for Ra = 10 5 , the effect is more pronounced, and the cell becomes bigger and has a tendency to move upward. Figure 10 Figure 11 illustrates the isotherms and streamlines for different nanoparticle diameter and Darcy numbers at Ra = 5 × 10 4 , ϕ = 0.05, ε = 0.5, and RR = 2. From Figure 11, for Da = 10 −4 , the isotherms have a uniform distribution due to the low permeability, which indicates that the flow is weak and thermal conduction plays a leading role. The effect of nanoparticle diameter on the heat transfer is weak. For Da = 10 −3 , due to the enhanced permeability, the natural convection heat transfer is strengthened. For Da = 10 −2 , the isotherms are almost horizontally distributed, especially when dsp = 10, which means that the natural convection heat transfer plays a significant role and the effect of nanoparticle diameter is more pronounced. In addition, for Da = 10 −4 and Da = 10 −3 , the effect of nanoparticle diameter on the streamlines is very weak, while for Da = 10 −2 , the effect is more pronounced; the cell becomes bigger and has a tendency to move upward. Figure 12 displays the effect of the nanoparticle diameter on the overall heat transfer rate along the inner wall at different Darcy numbers. The figure shows that an increase in the nanoparticle diameter leads to heat transfer enhancement for all considered Darcy numbers, and it can be found that the effect of nanoparticle diameter is increased with the increase of the Darcy number. Figure 13 presents the evolution of the local Nusselt number along the inner wall for different nanoparticle diameters at Ra = 10 5 , Da = 10 −2 , ϕ = 0.05, ε = 0.5, and RR = 2. It can be found that the local Nusselt number has a symmetrical distribution due to the symmetrical geometry and boundary conditions, and the maximum Nuloc occurred at γ= 180°, which means that the natural convection heat transfer is more intense in the bottom half of the inner wall. Furthermore, the local Nusselt number is increased with the decrease of nanoparticle diameter, which indicates that the nanoparticle diameter at high value will weaken the natural convection heat transfer.  Figure 11 illustrates the isotherms and streamlines for different nanoparticle diameter and Darcy numbers at Ra = 5 × 10 4 , φ = 0.05, ε = 0.5, and RR = 2. From Figure 11, for Da = 10 −4 , the isotherms have a uniform distribution due to the low permeability, which indicates that the flow is weak and thermal conduction plays a leading role. The effect of nanoparticle diameter on the heat transfer is weak. For Da = 10 −3 , due to the enhanced permeability, the natural convection heat transfer is strengthened. For Da = 10 −2 , the isotherms are almost horizontally distributed, especially when d sp = 10, which means that the natural convection heat transfer plays a significant role and the effect of nanoparticle diameter is more pronounced. In addition, for Da = 10 −4 and Da = 10 −3 , the effect of nanoparticle diameter on the streamlines is very weak, while for Da = 10 −2 , the effect is more pronounced; the cell becomes bigger and has a tendency to move upward. Figure 12 displays the effect of the nanoparticle diameter on the overall heat transfer rate along the inner wall at different Darcy numbers. The figure shows that an increase in the nanoparticle diameter leads to heat transfer enhancement for all considered Darcy numbers, and it can be found that the effect of nanoparticle diameter is increased with the increase of the Darcy number. Figure 13 presents the evolution of the local Nusselt number along the inner wall for different nanoparticle diameters at Ra = 10 5 , Da = 10 −2 , φ = 0.05, ε = 0.5, and RR = 2. It can be found that the local Nusselt number has a symmetrical distribution due to the symmetrical geometry and boundary conditions, and the maximum Nu loc occurred at γ= 180 • , which means that the natural convection heat transfer is more intense in the bottom half of the inner wall. Furthermore, the local Nusselt number is increased with the decrease of nanoparticle diameter, which indicates that the nanoparticle diameter at high value will weaken the natural convection heat transfer. Nanomaterials 2021, 11, x FOR PEER REVIEW 15 of 24     Figure 14 shows the isotherms and streamlines for different porosity and Rayleigh numbers at dsp = 50 nm, Da = 10 −2 , and RR = 2. From Figure 14, for Ra = 10 3 , the isotherms have a uniform distribution, and the isotherms are almost unchanged when the porosity increased from ε = 0.1 to ε = 0.9. This is because at low Rayleigh numbers, the buoyancy force is very weak, and the heat transfer in the annulus is dominated by thermal conduction. The porosity has a slight effect on the isotherms. For Ra = 10 4 , the isotherms have a similar trend with that for Ra = 10 3 at ε = 0.1, while the isotherms have an obvious change at ε = 0.9. For Ra = 5 × 10 4 , the natural convection heat transfer plays a leading role due to the enhanced flow, and the heat transfer is more intense at ε = 0.9 compared with ε = 0.1. The streamlines have a similar distribution, but some details are different. For Ra = 10 3 , the streamlines have little change at ε = 0.1 and ε = 0.9. For Ra = 10 4 and Ra = 5 × 10 4 , it can be found the cells become bigger and have a tendency to move upward when the porosity increases. Figure 15 displays the effect of porosity on the overall heat transfer rate along the inner wall at different Rayleigh numbers. It can be found that a continuous increase of the overall heat transfer rate occurs with the increase of porosity for all considered Rayleigh numbers and the effect of porosity is more pronounced at high Rayleigh numbers.    Figure 14 shows the isotherms and streamlines for different porosity and Rayleigh numbers at dsp = 50 nm, Da = 10 −2 , and RR = 2. From Figure 14, for Ra = 10 3 , the isotherms have a uniform distribution, and the isotherms are almost unchanged when the porosity increased from ε = 0.1 to ε = 0.9. This is because at low Rayleigh numbers, the buoyancy force is very weak, and the heat transfer in the annulus is dominated by thermal conduction. The porosity has a slight effect on the isotherms. For Ra = 10 4 , the isotherms have a similar trend with that for Ra = 10 3 at ε = 0.1, while the isotherms have an obvious change at ε = 0.9. For Ra = 5 × 10 4 , the natural convection heat transfer plays a leading role due to the enhanced flow, and the heat transfer is more intense at ε = 0.9 compared with ε = 0.1. The streamlines have a similar distribution, but some details are different. For Ra = 10 3 , the streamlines have little change at ε = 0.1 and ε = 0.9. For Ra = 10 4 and Ra = 5 × 10 4 , it can be found the cells become bigger and have a tendency to move upward when the porosity increases. Figure 15 displays the effect of porosity on the overall heat transfer rate along the inner wall at different Rayleigh numbers. It can be found that a continuous increase of the overall heat transfer rate occurs with the increase of porosity for all considered Rayleigh numbers and the effect of porosity is more pronounced at high Rayleigh numbers.  Figure 14 shows the isotherms and streamlines for different porosity and Rayleigh numbers at d sp = 50 nm, Da = 10 −2 , and RR = 2. From Figure 14, for Ra = 10 3 , the isotherms have a uniform distribution, and the isotherms are almost unchanged when the porosity increased from ε = 0.1 to ε = 0.9. This is because at low Rayleigh numbers, the buoyancy force is very weak, and the heat transfer in the annulus is dominated by thermal conduction. The porosity has a slight effect on the isotherms. For Ra = 10 4 , the isotherms have a similar trend with that for Ra = 10 3 at ε = 0.1, while the isotherms have an obvious change at ε = 0.9. For Ra = 5 × 10 4 , the natural convection heat transfer plays a leading role due to the enhanced flow, and the heat transfer is more intense at ε = 0.9 compared with ε = 0.1. The streamlines have a similar distribution, but some details are different. For Ra = 10 3 , the streamlines have little change at ε = 0.1 and ε = 0.9. For Ra = 10 4 and Ra = 5 × 10 4 , it can be found the cells become bigger and have a tendency to move upward when the porosity increases. Figure 15 displays the effect of porosity on the overall heat transfer rate along the inner wall at different Rayleigh numbers. It can be found that a continuous increase of the overall heat transfer rate occurs with the increase of porosity for all considered Rayleigh numbers and the effect of porosity is more pronounced at high Rayleigh numbers.    Figure 16, for Da = 10 −4 , the isotherms have a uniform distribution at ε = 0.1 and ε = 0.9; the reason is that the flow is constrained by the low permeability. In this case, porosity has little effect on the isotherms. For Da = 10 −3 , due to the enhanced permeability, the flow in the annulus is strengthened, and a disturbance occurs in the isotherms. For Da = 10 −2 , the isotherms at ε = 0.1 and ε = 0.9 have a remarkable difference. The isotherms are almost horizontally distributed at ε = 0.9, which means that the overall heat transfer is dominated by natural convection. For Da = 10 −4 , the streamlines are almost the same at ε = 0.1 and ε = 0.9. For Da = 10 −3 and Da = 10 −2 , the cell becomes bigger and has a tendency to move upward when porosity increases from ε = 0.1 to ε = 0.9. Figure 17 displays the effect of porosity on the overall heat transfer rate along the inner wall at different Darcy number. The figure shows that an increase in porosity leads to heat transfer enhancement for all considered Darcy numbers, and the effect of porosity is more pronounced at high Darcy numbers. Figure 18 presents the evolution of the local Nusselt number along the inner wall for different porosity at Ra = 5 × 10 4 , Da = 10 −2 , ϕ = 0.05, ε = 0.5, dsp = 50 nm, and RR = 2. The local Nusselt number has a symmetrical distribution due to the symmetrical geometry and boundary conditions, and the maximum Nuloc occurred at γ = 180°. Furthermore, the local Nusselt number is increased with the increase of porosity, so it can be concluded that the porosity has a positive effect on the overall heat transfer rate.   Figure 16, for Da = 10 −4 , the isotherms have a uniform distribution at ε = 0.1 and ε = 0.9; the reason is that the flow is constrained by the low permeability. In this case, porosity has little effect on the isotherms. For Da = 10 −3 , due to the enhanced permeability, the flow in the annulus is strengthened, and a disturbance occurs in the isotherms. For Da = 10 −2 , the isotherms at ε = 0.1 and ε = 0.9 have a remarkable difference. The isotherms are almost horizontally distributed at ε = 0.9, which means that the overall heat transfer is dominated by natural convection. For Da = 10 −4 , the streamlines are almost the same at ε = 0.1 and ε = 0.9. For Da = 10 −3 and Da = 10 −2 , the cell becomes bigger and has a tendency to move upward when porosity increases from ε = 0.1 to ε = 0.9. Figure 17 displays the effect of porosity on the overall heat transfer rate along the inner wall at different Darcy number. The figure shows that an increase in porosity leads to heat transfer enhancement for all considered Darcy numbers, and the effect of porosity is more pronounced at high Darcy numbers. Figure 18 presents the evolution of the local Nusselt number along the inner wall for different porosity at Ra = 5 × 10 4 , Da = 10 −2 , φ = 0.05, ε = 0.5, d sp = 50 nm, and RR = 2. The local Nusselt number has a symmetrical distribution due to the symmetrical geometry and boundary conditions, and the maximum Nu loc occurred at γ = 180 • . Furthermore, the local Nusselt number is increased with the increase of porosity, so it can be concluded that the porosity has a positive effect on the overall heat transfer rate.    Figure 19 illustrates the isotherms and streamlines for different radius ratio at Ra = 5 × 10 4 , Da = 10 −2 , ϕ = 0.05, dsp = 50 nm, and ε = 0.5. The radius ratio is an important parameter that affects the fluid flows and natural convection heat transfer according the previous research [10][11][12]. From Figure 19, when the radius ratio increased from RR = 1.2 to RR= 6, the flow is strengthened, a main cell is formed in the middle of the annulus, and the cell becomes bigger and has a tendency to move upward. It is worth noting that the color scales on the right have a downward trend with the increase of RR, which is inverse to the above conclusion. In fact, the reason is that the change of AR leads to a change of the characteristic length (ro-ri). The isotherms distribution is changed from vertical to horizontal, which means that the heat transfer is dominated by natural convection. When RR = 1.2, multicellular flow structures are formed. Figure 20 displays the effect of radius ratio on the overall heat transfer rate along the inner wall. It can be found that heat transfer in the annulus is enhanced with the increase of radius ratio from RR = 1.3 to RR = 6, but a fall occurred in the range of RR = 1.2 to RR = 1.3. A conclusion can be drawn that a bifurcation point exists when the radius ratios is in the range of RR = 1.2 to RR = 1.3 for the considered parameters.     Figure 19 illustrates the isotherms and streamlines for different radius ratio at Ra = 5 × 10 4 , Da = 10 −2 , ϕ = 0.05, dsp = 50 nm, and ε = 0.5. The radius ratio is an important parameter that affects the fluid flows and natural convection heat transfer according the previous research [10][11][12]. From Figure 19, when the radius ratio increased from RR = 1.2 to RR= 6, the flow is strengthened, a main cell is formed in the middle of the annulus, and the cell becomes bigger and has a tendency to move upward. It is worth noting that the color scales on the right have a downward trend with the increase of RR, which is inverse to the above conclusion. In fact, the reason is that the change of AR leads to a change of the characteristic length (ro-ri). The isotherms distribution is changed from vertical to horizontal, which means that the heat transfer is dominated by natural convection. When RR = 1.2, multicellular flow structures are formed. Figure 20 displays the effect of radius ratio on the overall heat transfer rate along the inner wall. It can be found that heat transfer in the annulus is enhanced with the increase of radius ratio from RR = 1.3 to RR = 6, but a fall occurred in the range of RR = 1.2 to RR = 1.3. A conclusion can be drawn that a bifurcation point exists when the radius ratios is in the range of RR = 1.2 to RR = 1.3 for the considered parameters.  Figure 19 illustrates the isotherms and streamlines for different radius ratio at Ra = 5 × 10 4 , Da = 10 −2 , φ = 0.05, d sp = 50 nm, and ε = 0.5. The radius ratio is an important parameter that affects the fluid flows and natural convection heat transfer according the previous research [10][11][12]. From Figure 19, when the radius ratio increased from RR = 1.2 to RR= 6, the flow is strengthened, a main cell is formed in the middle of the annulus, and the cell becomes bigger and has a tendency to move upward. It is worth noting that the color scales on the right have a downward trend with the increase of RR, which is inverse to the above conclusion. In fact, the reason is that the change of RR leads to a change of the characteristic length (r o -r i ). The isotherms distribution is changed from vertical to horizontal, which means that the heat transfer is dominated by natural convection. When RR = 1.2, multicellular flow structures are formed. Figure 20 displays the effect of radius ratio on the overall heat transfer rate along the inner wall. It can be found that heat transfer in the annulus is enhanced with the increase of radius ratio from RR = 1.3 to RR = 6, but a fall occurred in the range of RR = 1.2 to RR = 1.3. A conclusion can be drawn that a bifurcation point exists when the radius ratios is in the range of RR = 1.2 to RR = 1.3 for the considered parameters.

Conclusions
Natural convection heat transfer in a porous annulus filled with a Cu-Nanofluid has been investigated numerically. The effects of Brownian motion, solid volume fraction, nanoparticle diameter, Rayleigh number, Darcy number, porosity on the flow pattern, temperature distribution, and heat transfer characteristics are discussed in detail. The following conclusions could be drawn: (1) Brownian motion should be considered in the natural convection heat transfer of a nanofluid.

Conclusions
Natural convection heat transfer in a porous annulus filled with a Cu-Nanofluid has been investigated numerically. The effects of Brownian motion, solid volume fraction, nanoparticle diameter, Rayleigh number, Darcy number, porosity on the flow pattern, tem-perature distribution, and heat transfer characteristics are discussed in detail. The following conclusions could be drawn: (1) Brownian motion should be considered in the natural convection heat transfer of a nanofluid.
The effect of Brownian motion becomes more remarkable with the increase of Rayleigh number and nanoparticle volume fraction, while it is less pronounced with the increase of Darcy number and nanoparticle diameter.
(2) The increase of nanoparticle volume fraction results in an improvement of the overall heat transfer rate, and the effect is more remarkable when the Rayleigh and Darcy numbers are at high value.
(3) Increasing the nanoparticle diameter has a negative effect on the overall heat transfer rate and the effect has a limit when the nanoparticle diameter reaches a high value.
(4) The porosity affects the flow pattern, temperature distribution, and heat transfer rate. The flow motion is limited when the porosity is too low. The heat transfer rate is strengthened with the increase of porosity, especially with high Rayleigh and Darcy numbers.
(5) The radius ratio has a significant influence on the isotherms, streamlines, and heat transfer rate. The rate is greatly enhanced with the increase of radius ratio. Additionally, when the radius ratio is too low, multicellular flow structures are formed, and a bifurcation point exists.