Background {#Sec1}
==========

Nanofluid is a suspension containing a certain quantity of nanoscaled solid particles in a conventional cooling liquid, such as water and ethylene glycol \[[@CR1]\]. Nanofluid shows considerably better heat transfer performance than single-phase mediums due to particle's Brownian motion and interaction \[[@CR2],[@CR3]\]. Furthermore, because of the ultra-small particle size, nanofluid is also remarkably better than normal multi-phase fluid to eliminate erosion and clogging problems in micro channels \[[@CR4],[@CR5]\]. Recently, nanofluid is increasingly used in natural convection applications for wide areas \[[@CR6]\], such as electronic cooling, heat exchangers, boilers, nuclear reactor systems and energy storage devices \[[@CR7]\].

To gain a better understanding of nanofluid natural heat convection, many studies have been carried out in both experimental and numerical ways during the past decade \[[@CR8]-[@CR10]\]. However, some apparently different conclusions can be found in those experimental and numerical investigations \[[@CR11]\]. Briefly about nanofluid natural heat convection, deterioration was usually illustrated by experimental studies, while enhancement was always reported by numerical studies.

By experimental study, Putra et al. \[[@CR12]\] found heat transfer deterioration in Al~2~O~3~/water and water-based copper oxide (CuO/water) nanofluids (with volume fraction 1% and 4%). They ascribed the possible reasons to particle-fluid slip and nanoparticle sedimentation. A similar observation was also reported by Wen and Ding \[[@CR13]\]. In their experiments, 0.19% \~ 0.57% water-based titanium oxide (TiO~2~/water) nanofluids had a lower natural convective heat transfer coefficient than pure water, and the deterioration increased with volume fraction. They supposed the convection induced by particle concentration difference to be a possible reason. Li and Peterson \[[@CR14]\] reported a natural heat convection deterioration in 0.5% \~ 6% Al~2~O~3~/water nanofluids. They inferred the reason could be the nanoparticle's Brownian motion smoothing the temperature gradient leading to the delay of natural convection. Also, higher viscosity of nanofluids could also induce such an effect. Ni et al. \[[@CR15]\] reported deteriorated natural convection after their experiments for 1.08% Al~2~O~3~/water nanofluid in a Rayleigh-Bénard configuration. They suggested that the significant decrease might be caused by the mass diffusion of nanoparticles. In Nnanna's experiment \[[@CR16]\], it was found that the presence of Al~2~O~3~ nanoparticles did not impede the water-free heat convection when the volume fraction was in the range of 0.2% \~ 2%. However, the heat convection declined due to increase of kinematic viscosity since the volume fraction was larger than 2%. Ho et al. \[[@CR17]\] also reported up to 18% natural convective heat transfer enhancement in 0.1% Al~2~O~3~/water nanofluid, but degradation was found when the volume fraction was larger than 2%.

In numerical studies, some excellent molecular dynamic works have been carried out recently to analyse the possible heat transfer mechanisms between nanoparticle and fluid \[[@CR18],[@CR19]\], such as Chiavazzo's research on nanofins \[[@CR20],[@CR21]\]. However, today, normal computational fluid dynamics (CFD) investigation is still playing a dominant role in this area. Khanafer et al. \[[@CR22]\] found the water heat transfer rate in a two-dimensional enclosure could be substantially increased by adding more copper nanoparticles (*ϕ* was up to 20%). Oztop and Abu-Nada \[[@CR23]\] found natural heat transfer enhancement by using Al~2~O~3~/water, TiO~2~/water and Cu/water nanofluids (*ϕ* was up to 20%) in two-dimensional rectangular enclosures with different aspect ratios. Furthermore, the enhancement was more pronounced at a low aspect ratio than at a high aspect ratio. Aminossadati and Ghasemi \[[@CR24]\] reported that adding copper (Cu), silver (Ag), Al~2~O~3~ and TiO~2~ nanoparticles (*ϕ* was up to 20%) could improve cooling performance of pure water in a bottom-heated two-dimensional enclosure, especially when the Rayleigh number was low. Ghasemi and Aminossadati \[[@CR25]\] reported a larger CuO/water nanofluid volume fraction (*ϕ* = 1% \~ 4%) led to Nusselt number enhancement in a two-dimensional triangular enclosure. Oueslati et al. \[[@CR26]\] found nanofluid natural heat convection enhancement in a two-dimensional cavity when the volume fraction of Al~2~O~3~, TiO~2~ and Cu nanoparticles was lower than 5%. Ternik et al. \[[@CR27]\] examined the heat transfer enhancement of water-based gold (Au), Al~2~O~3~, Cu and TiO~2~ nanofluids (*ϕ* was up to 10%) in a two-dimensional cavity. They indicated that the average Nusselt number was an increasing function of nanofluid volume fraction.

Actually, besides the above controversial conclusions from experimental and numerical investigations, dimension issue is another problem that should be noticed. In the past years, nearly all those experiments were performed in cylinders or tubes. They were of course three-dimensional investigations. However, previous numerical simulations were always performed for two-dimensional cases \[[@CR22]-[@CR27]\]. In fact, due to the gravity-induced force acting perpendicularly to the horizontal cylinder wall, the fluid movement in the cylinder could be apparently three-dimensional \[[@CR12]\]. This could be an important factor influencing the internal flow and natural heat transfer behaviour in a horizontal cylinder. In previous publications, however, it is rare to see three-dimensional CFD investigations for correlative cases.

The aim of this study is to investigate Al~2~O~3~/water nanofluid natural heat convection in a horizontal cylinder by a three-dimensional CFD approach. For this, the numerical CFD package OpenFOAM \[[@CR28]\] will be employed for case simulation in the present work. Furthermore, a new OpenFOAM solver will be developed to count in the possible impacts from temperature-dependent fluid property variation and make the numerical simulation more realistic.

Methods {#Sec2}
=======

Problem configuration {#Sec3}
---------------------

In this investigation, a part of the experiment reported by Putra et al. \[[@CR12]\] in 2003 is selected as the modelling prototype and will be repeated numerically. The working session of this experimental device (Figure [1](#Fig1){ref-type="fig"}A) is simplified to be an insulated cylinder with heating and cooling walls at the two ends (Figure [1](#Fig1){ref-type="fig"}B). The cylinder length and diameter are given as *L* = 0.04 m and *D* = 0.04 m, respectively. The ratio of cylinder length to diameter (*L*/*D*) is 1.0. A constant temperature *T*~C~ is given at the cooling end, and the input power is controlled at the heating end with temperature *T*~H~ to obtain the cases with different Rayleigh numbers. Non-permeable and non-slip boundaries are assumed for velocity. Blender (an open-source animation suite) \[[@CR29]\] is used to create the three-dimensional geometry model. In this horizontal cylinder, natural heat convection of water and 1% and 4% Al~2~O~3~/water nanofluids will be investigated.Figure 1**Modelling prototype used in the present study.** Working session of Putra's experimental device **(A)** \[[@CR12]\]. Details: On 19th August 2014, Prof. Nandy Putra at the University of Indonesia authorized us to use his original figure which was published in \[[@CR12]\]. Schematic experimental model, in which *D* is the cylinder diameter, *L* is the cylinder length, *T* ~C~ is the temperature of the cooling wall and *T* ~H~ is the temperature of the heating wall **(B)**.

Recently, some researchers began to try different numerical methods and treated nanofluid as a multi-phase mixture, e.g. Eulerian \[[@CR30],[@CR31]\] and Lagrangian approaches \[[@CR32],[@CR33]\]. Even in some single-phase CFD simulations, the effect of nanoparticle-fluid slip began to be included, such as in Aminfar's study \[[@CR34]\]. However, there are still some arguments on the issue as to whether the multi-phase approach is better or not for nanofluids \[[@CR35],[@CR36]\]. As the mechanisms of nanofluid thermal conductivity enhancement are still not very clear \[[@CR37],[@CR38]\], the single-phase method is still used in the present work and Al~2~O~3~/water nanofluid is assumed to be a stable and homogeneous mixture as in traditional numerical investigations.

Flow model {#Sec4}
----------

Since a stable purely laminar flow state could be eventually obtained in Putra's experiment \[[@CR12]\], original OpenFOAM solver 'buoyantBoussinesqSimpleFoam' is selected for this steady-state, buoyant, incompressible fluid case. In this solver, there are three coupled partial differential equations to describe flow and heat transfer problems; they are mass equation ([1](#Equ1){ref-type=""}), momentum equation ([2](#Equ2){ref-type=""}) and energy equation ([3](#Equ3){ref-type=""}) \[[@CR39]\]:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ \mathit{\nabla}\cdot u=0 $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ \rho \mathit{\nabla}\cdot (uu)=-\mathit{\nabla}\rho +\mathit{\nabla}\cdot \left(\mu \mathit{\nabla}u\right)+{g}_k $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ \rho \mathit{\nabla}\cdot \left({c}_pTu\right)=\mathit{\nabla}\cdot \left(k\mathit{\nabla}T\right) $$
                        \end{document}$$where *u*, *p*, *T*, *ρ*, *c*~*p*~, *k*, *μ* and *g*~*k*~ indicate the velocity, pressure, temperature, density, specific heat capacity, thermal conductivity, dynamic viscosity and body force, respectively. '∇' and '∇' indicate divergence and gradient operations, respectively. Body force *g*~*k*~ is evaluated by Boussinesq approximation ([4](#Equ4){ref-type=""}) \[[@CR40]\]:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {g}_k=\rho \left[1.0-\beta \left(T-{T}_{\mathrm{Ref}}\right)\right]g $$
                        \end{document}$$where *β*, *T*~Ref~ and *g* are the thermal expansion coefficient, reference temperature and gravity, respectively.

Al~2~O~3~/water nanofluid density *ρ*~nf~ is evaluated by ([5](#Equ5){ref-type=""}) \[[@CR22],[@CR23]\]:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {\rho}_{\mathrm{nf}}=\left(1-\varnothing \right){\rho}_{\mathrm{f}}+\varnothing {\rho}_{\mathrm{s}} $$
                        \end{document}$$

Al~2~O~3~/water nanofluid heat capacity *c*~*p*nf~ is evaluated by ([6](#Equ6){ref-type=""}) \[[@CR22],[@CR23]\]:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {\rho}_{\mathrm{nf}}{c}_{p\mathrm{n}\mathrm{f}}=\left(1-\varnothing \right){\left(\rho {c}_p\right)}_{\mathrm{f}}+\varnothing {\left(\rho {c}_p\right)}_{\mathrm{s}} $$
                        \end{document}$$where subscripts f, s and nf indicate fluid, solid and nanofluid, respectively.

Water properties are collected from a heat transfer textbook \[[@CR41]\]. Density and specific heat capacity of Al~2~O~3~ nanoparticles are used as 3,970 kg m^−3^ and 765 J kg^−1^ K, respectively \[[@CR23]\]. To ensure the simulation is more reliable, thermal conductivity *k* and kinetic viscosity *v* of 1% and 4% Al~2~O~3~/water nanofluids are collected from Das's study \[[@CR12],[@CR42]\] instead of prediction models. Based on regression analysis, thermal conductivity enhancement variation will be given by ([7](#Equ7){ref-type=""}) \~ ([9](#Equ9){ref-type=""}), while kinetic viscosity variation will be given by ([10](#Equ10){ref-type=""}) \~ ([12](#Equ12){ref-type=""}). For original solver 'buoyantBoussinesqSimpleFoam' , fluid properties are considered as constant for each case and evaluated at the mean of the heating and cooling wall temperatures (in Celsius degree).$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {k}_{\mathrm{water}}=0.0027\left(T-273\right)+0.9716 $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {k}_{1\%}=0.0028\left(T-273\right)+0.9672 $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {k}_{4\%}=0.0053\left(T-273\right)+0.977 $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {v}_{\mathrm{water}}=0.0002{\left(T-273\right)}^2-0.0305\left(T-273\right)+1.52 $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {v}_{1\%}=0.0003{\left(T-273\right)}^2-0.045\left(T-273\right)+1.96 $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {v}_{4\%}=0.0003{\left(T-273\right)}^2-0.0448\left(T-273\right)+2.16 $$
                        \end{document}$$

Based on 'buoyantBoussinesqSimpleFoam' , a new solver, 'buoyantBoussinesqSimpleTDFoam' , is developed to include possible impacts from fluid property variation due to temperature change. This temperature-dependent solver assumes that the fluid properties in a computational region are not uniform but decided by each volume cell's temperature. During numerical simulation, strongly temperature-dependent fluid properties (e.g. thermal conductivity and viscosity) will be updated for each volume cell at each iteration after the energy equation has been solved.

Dimensionless parameter Nusselt number *Nu* is used to describe the fluid natural heat convection performance in those cases with different Rayleigh numbers *Ra*. They are defined in ([13](#Equ13){ref-type=""}) and ([14](#Equ14){ref-type=""}), respectively \[[@CR22],[@CR23]\]:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ Nu=\frac{h{L}_{\mathrm{c}}}{k} $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ Ra=\frac{\beta g\left({T}_{\mathrm{H}}-{T}_{\mathrm{C}}\right){L_{\mathrm{c}}}^3}{v\alpha } $$
                        \end{document}$$where *h*, *L*~c~, *v* and *α* indicate the heat transfer coefficient, characteristic length, kinetic viscosity and thermal diffusivity, respectively. Convection heat transfer coefficient *h* is defined in ([15](#Equ15){ref-type=""}) \[[@CR22],[@CR23]\]:$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ h=\frac{4Q}{\pi {D}^2\left({T}_{\mathrm{H}}-{T}_{\mathrm{C}}\right)} $$
                        \end{document}$$where *Q* is the changeable input power, by which different Rayleigh numbers can be obtained for different cases.

Numerical implementation and grid independence check {#Sec5}
----------------------------------------------------

The set of three-dimensional coupled non-linear differential equations are discretized by control volume technique \[[@CR43]\]. The semi-implicit method for pressure-linked equations (SIMPLE) algorithm \[[@CR44]\] is employed in both 'buoyantBoussinesqSimpleFoam' and 'buoyantBoussinesqSimpleTDFoam' to solve Navier-Stokes equations. The spatial schemes for gradient, Laplacian and divergence are Gauss linear, Gauss linear corrected and Gauss linear schemes, respectively. During the iterative process, the absolute residuals of *u*, *p* and *T* are carefully monitored and convergence criteria for every parameter is restricted below 10^−6^.

An extensive testing procedure is carried out to guarantee a grid-independent solution. By enGrid (an open-source mesh generation package) \[[@CR45]\], unstructured tetras are used to fill the three-dimensional cylinder. Four mesh strategies are generated for the grid independence check; their maximal cell edge lengths Δ*d* on heating and cooling walls are 4 mm (Figure [2](#Fig2){ref-type="fig"}A), 3 mm (Figure [2](#Fig2){ref-type="fig"}B), 2 mm (Figure [2](#Fig2){ref-type="fig"}C) and 1 mm (Figure [2](#Fig2){ref-type="fig"}D). To ensure simulation accuracy for near-wall regions, a non-uniform strategy is used to refine the mesh near heating and cooling walls. Original OpenFOAM solver 'buoyantBoussinesqSimpleFoam' is employed to test the four mesh strategies for a water case *Ra* = 10^8^. Since this solver has been validated in Corzo et al.'s work \[[@CR46]\], a similar validation will not be repeated and the solver will be used directly.Figure 2**Four mesh strategies for grid independence check:** Δ***d***  **= 4 mm (A), 3 mm (B), 2 mm (C) and 1 mm (D).**

Dimensionless temperature *T*\* and location *X*\* are defined in ([16](#Equ16){ref-type=""}) and ([17](#Equ17){ref-type=""}) for the grid independence check, respectively. Results of *T*\* and *X*\* on the cylinder longitudinal central line are compared to find the most appropriate mesh strategy. By the comparison in Figure [3](#Fig3){ref-type="fig"}, mesh strategies Δ*d* = 2 mm and Δ*d* = 1 mm are found to predict nearly exactly the same results. This indicates the mesh strategy Δ*d* = 2 mm is good enough for the present study. However, the strategy Δ*d* = 1 mm is eventually selected to capture even more detailed velocity and temperature features in near-wall regions. Compared to normal two-dimensional simulations, the cell amount is considerably increased in this work. Although non-uniform grid strategies are employed to reduce the total cell number, there are still about 0.4 million cells used in the present work (Figure [2](#Fig2){ref-type="fig"}D).Figure 3**Grid independence check.**$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {T}^{*}=\frac{T-{T}_{\mathrm{C}}}{T_{\mathrm{H}}-{T}_{\mathrm{C}}} $$
                        \end{document}$$$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                            $$ {X}^{*}=\frac{X}{L_{\mathrm{c}}} $$
                        \end{document}$$

Results and discussion {#Sec6}
======================

In this study, natural heat convection of water and 1% and 4% Al~2~O~3~/water nanofluids is investigated numerically by both 'buoyantBoussinesqSimpleFoam' and 'buoyantBoussinesqSimpleTDFoam'. For each fluid, five simulations are performed in the Rayleigh number range *Ra* = 10^7^ \~ 0.8 × 10^8^. Figures [4](#Fig4){ref-type="fig"}, [5](#Fig5){ref-type="fig"}, [6](#Fig6){ref-type="fig"}, [7](#Fig7){ref-type="fig"} and [8](#Fig8){ref-type="fig"} show the results of the average Nusselt number against the Rayleigh number in different cases. Basically, by both the present numerical study and previous experimental study, the natural Nusselt number of water and Al~2~O~3~/water nanofluids is found to increase with the Rayleigh number. However, some more information also can be found by further comprehensive comparisons in the present work.Figure 4**Results of Nusselt number predicted by the solver 'buoyantBoussinesqSimpleFoam'.**Figure 5**Results of Nusselt number predicted by the solver 'buoyantBoussinesqSimpleTDFoam'.**Figure 6**Comparison of Nusselt number in water cases.**Figure 7**Comparison of Nusselt number in 1% Al** ~**2**~ **O** ~**3**~ **/water nanofluid cases.**Figure 8**Comparison of Nusselt number in 4% Al** ~**2**~ **O** ~**3**~ **/water nanofluid cases.**

Figure [4](#Fig4){ref-type="fig"} shows the Nusselt number results predicted by the original solver 'buoyantBoussinesqSimpleFoam' for the cases of water and 1% and 4% Al~2~O~3~/water nanofluids. It can be found that 1% Al~2~O~3~/water nanofluid has no apparently different natural Nusselt numbers with water, but 4% Al~2~O~3~/water nanofluid has slightly lower Nusselt numbers. The possible reason is that although a nanofluid with a larger volume fraction has larger thermal conductivity, a larger volume fraction also induces larger fluid viscosity. The increased viscosity probably plays a predominant role to impede nanofluid natural heat convection.

Figure [5](#Fig5){ref-type="fig"} shows the Nusselt number results predicted by the temperature-dependent solver 'buoyantBoussinesqSimpleTDFoam' for the cases of water and 1% and 4% Al~2~O~3~/water nanofluids. It can be found that 1% Al~2~O~3~/water nanofluid has no apparently different natural Nusselt numbers with water. But 4% Al~2~O~3~/water nanofluid has lower natural Nusselt numbers. Compared to Figure [4](#Fig4){ref-type="fig"}, the two solvers give a very similar conclusion: the nanofluid does not have a better natural convective heat transfer performance than the basefluid. This actually conforms to those experimental results \[[@CR12],[@CR13]\], but in contradiction to some previous numerical conclusions \[[@CR24],[@CR25]\]. For the difference between the present simulations and other numerical studies, the reason inferred is the possible impacts from the cylinder wall (in previous two-dimensional investigations, wall effect from the third direction could not be included to influence fluid flow and heat transfer).

Figure [6](#Fig6){ref-type="fig"} shows the Nusselt number predictions from the original solver, temperature-dependent solver and experimental study for water cases. It can be found that the original solver 'buoyantBoussinesqSimpleFoam' predicts an apparently lower Nusselt number than the experimental study, but the results from the temperature-dependent solver 'buoyantBoussinesqSimpleTDFoam' have good agreement to experimental data, particularly when the Rayleigh number is larger than 2.0 × 10^7^. Figure [7](#Fig7){ref-type="fig"} shows a similar comparison for 1% Al~2~O~3~/water nanofluid cases. It can be found that the original solver 'buoyantBoussinesqSimpleFoam' does not predict good values in the whole Rayleigh number range, but the results from the new solver 'buoyantBoussinesqSimpleTDFoam' begin to have good agreement with experimental data since the Rayleigh number is larger than 2.0 × 10^7^. By the comparisons in Figures [6](#Fig6){ref-type="fig"} and [7](#Fig7){ref-type="fig"}, it can be concluded that the temperature-dependent solver 'buoyantBoussinesqSimpleTDFoam' is better for water and Al~2~O~3~/water nanofluids with a low volume fraction.

Figure [8](#Fig8){ref-type="fig"} shows the Nusselt number predictions from the original solver, new solver and experimental study for 4% Al~2~O~3~/water nanofluid cases. It can be found that the predictions from the original solver 'buoyantBoussinesqSimpleFoam' begin to have good agreement with experimental data since the Rayleigh number is larger than 2.0 × 10^7^. However, the newly developed solver 'buoyantBoussinesqSimpleTDFoam' predicts a larger Nusselt number than the experimental study. For this phenomenon, nanoparticle sedimentation is thought to be a possible reason. This is also reported by some recent nanoparticle sedimentation observations \[[@CR37],[@CR47]\]. Due to an improper nanoparticle dispersion method after fabrication, the nanofluid with a higher volume fraction will show a settlement layer quickly at the vessel bottom, leading to considerable deterioration of heat transfer performance. However, this mechanism is not included in the present single-phase CFD approach.

More details of cylinder internal flow features are disclosed in this work. With the help of ParaView (an open-source visualization application) \[[@CR48]\], Figures [9](#Fig9){ref-type="fig"}, [10](#Fig10){ref-type="fig"} and [11](#Fig11){ref-type="fig"} show the velocity vector features on three typical cylinder cross section positions - *X* = 0.01 m, *X* = 0.02 m and *X* = 0.03 m - for a 4% Al~2~O~3~/water nanofluid case with *Ra* = 4 × 10^7^. In the three figures, considerably horizontal and asymmetric velocity components can be observed, which indicate that the temperature-driven flow in the cylinder is actually different than that in broad rectangular ducts \[[@CR49],[@CR50]\]. Therefore, due to the strong three-dimensional flow features, applying two-dimensional simplification and neglecting the possible impacts from the cylinder wall may be not so realistic. Briefly, for natural heat convection problems in a horizontal cylinder, a three-dimensional CFD simulation is assumed to be more reliable than two-dimensional simplifications \[[@CR51]\].Figure 9**Velocity features on cylinder cross section** ***X***  **= 0.01 m.**Figure 10**Velocity features on cylinder cross section** ***X***  **= 0.02 m.**Figure 11**Velocity features on cylinder cross section** ***X***  **= 0.03 m.**

Conclusions {#Sec7}
===========

In this study, natural heat convection of Al~2~O~3~/water nanofluids in a horizontal cylinder (*L*/*D* = 1) is numerically investigated. The whole three-dimensional CFD process is performed in a completely open-source way. Blender, enGrid, OpenFOAM and ParaView are employed for geometry model creation, mesh generation, case simulation and post process, respectively. Natural heat convection of water and 1% and 4% Al~2~O~3~/water nanofluids is investigated by both the original OpenFOAM solver 'buoyantBoussinesqSimpleFoam' and the newly developed temperature-dependent solver 'buoyantBoussinesqSimpleTDFoam'. Based on the obtained results in this investigation, several conclusions can be drawn as follows:The natural convective Nusselt number of both water and Al~2~O~3~/water nanofluids in the horizontal cylinder (*L*/*D* = 1) increases with the Rayleigh number in the range of *Ra* = 0.7 × 10^7^ \~ 5 × 10^7^. With a given Rayleigh number, a larger nanofluid volume fraction induces a lower Nusselt number.Numerical solvers in the OpenFOAM frame begin to give better agreements to experimental investigations since the Rayleigh number is larger than 2 × 10^7^. The temperature-dependent solver is better for water and 1% Al~2~O~3~/water nanofluid cases, but the original solver is better for 4% Al~2~O~3~/water nanofluid cases.Due to strong three-dimensional flow features being observed, the three-dimensional CFD simulation is recommended for natural heat convection problems in a horizontal cylinder. This would be the best way to account for the possible impacts from the cylinder wall.

In a future study, a new OpenFOAM solver will be designed based on the present temperature-dependent one and treat nanofluid as a multi-phase mixture. Furthermore, it will be a transient solver instead of the currently steady-state one in this work. By applying proper mixture models, the interaction between nanofluid natural heat transfer and nanoparticle sedimentation will be observed (Additional file [1](#MOESM1){ref-type="media"}). This is necessary to figure out whether the nanoparticle sedimentation and nanoparticle-fluid slip have a remarkable impact on nanofluid natural heat convection in a horizontal cylinder.

Nomenclature {#Sec8}
------------

*c*~*p*~, Specific heat capacity (J kg^−1^ K^−1^)Δ*d*, Maximum cell edge length (m)*D*, Cylinder diameter (m)*g*, Gravity (m s^−2^)*g*~*k*~, Boussinesq body force (kg m^−2^ s^−2^)*h*, Heat transfer coefficient (W m^−2^ K^−1^)*k*, Thermal conductivity (W m^−1^ K^−1^)*L*, Cylinder length (m)*L*~c~, Characteristic length (m)$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                                        $$ Nu=\frac{h{L}_{\mathrm{c}}}{k} $$
                                    \end{document}$, Nusselt number*p*, Pressure (kg m^−1^ s^−2^)*Q*, Input power (W)*u*, Velocity (m s^−1^)*T*, Temperature (K)*T*~C~, Temperature of cooling wall (K)*T*~H~, Temperature of heating wall (K)*T*~Ref~, Reference temperature (K)*T*\*, Dimensionless temperature*X*\*, Dimensionless position$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}
                                        $$ Ra=\frac{\beta g\varDelta T{L_{\mathrm{c}}}^3}{v\alpha } $$
                                    \end{document}$, Rayleigh number

Greek symbols*α*, Thermal diffusivity (m^2^ s^−1^)*β*, Thermal expansion coefficient (K^−1^)*μ*, Dynamic viscosity (NS m^−2^)*v*, Kinetic viscosity (m^2^ s^−1^)*ρ*, Density (kg m^−3^)*ϕ*, Nanofluid volume fraction

Subscriptsf, Fluidnf, Nanofluids, Solid

Additional file {#Sec9}
===============

Additional file 1:**Cover letter to referee 3 and referee 4's concerns.**

**Competing interests**

The authors declare that they have no competing interests.

**Authors' contributions**

XM carried out the literature review, code design and numerical simulations. YL participated in code optimization and data processing. Both authors discussed the results and approved the final manuscript.

**Authors' information**

Dr. XM is a teaching fellow at the School of Marine Science and Technology, Newcastle University, UK. His research focuses mainly on multi-phase fluid flow, heat transfer and flow-induced vibration. Dr. YL is a postdoc researcher at Dalian Institute of Chemical Physics, Chinese Academy of Sciences. His research focuses mainly on molecular modelling, design and optimization.

This research is supported by DORA funding provided by the School of Marine Science and Technology, Newcastle University. We also thank Prof. Nandy Putra at the University of Indonesia for his kindly permission to use his published experimental data for our numerical solvers' validation. Dr. Xiangyin Meng would like to thank Prof. Qingwei Ma at the City University London for his continuous farther-like support in the past three years, as well as the kindly help from Mr. Neil Forbes at Newcastle University.
