Numerical investigation of nanoparticles slip mechanisms impact on the natural convection heat transfer characteristics of nanofluids in an enclosure

This study intends to give qualitative results toward the understanding of different slip mechanisms impact on the natural heat transfer performance of nanofluids. The slip mechanisms considered in this study are Brownian diffusion, thermophoretic diffusion, and sedimentation. This study compares three different Eulerian nanofluid models; Single-phase, two-phase, and a third model that consists of incorporating the three slip mechanisms in a two-phase drift-flux. These slip mechanisms are found to have different impacts depending on the nanoparticle concentration, where this effect ranges from negligible to dominant. It has been reported experimentally in the literature that, with high nanoparticle volume fraction the heat transfer deteriorates. Admittingly, classical nanofluid models are known to underpredict this impairment. To address this discrepancy, this study focuses on the effect of thermophoretic diffusion and sedimentation outcome as these two mechanisms turn out to be influencing players in the resulting heat transfer rate using the two-phase model. In particular, the necessity to account for the sedimentation contribution toward qualitative modeling of the heat transfer is highlighted. To this end, correlations relating the thermophoretic and sedimentation coefficients to the nanofluid concentration and Rayleigh number are proposed in this study. Numerical experiments are presented to show the effectiveness of the proposed two-phase model in approaching the experimental data, for the full range of Rayleigh number in the laminar flow regime and for nanoparticles concentration of (0% to 3%), with great satisfaction.

the work of Keblinski et al. 6 and Cui et al. 7 amongst others) and nanoscale level models using Boltzmann transport equation (e.g., the study conducted by Sheikholeslami et al. 8 ) to methods based on continuum mechanics. The latter can be also subdivided into; single-phase approach, two-phase approach based on the Eulerian-Lagrangian method (see for example the work of Sharaf et al. 9 ), or two-phase approach based on the Eulerian-Eulerian method, which is the focus of the present study. Considering this, a brief review of the evolution of Eulerian-Eulerian based models used for the prediction of nanofluid heat transfer features in the square cavity is presented here.
At the nascent stage of numerical methods development for nanofluid heat transfer, the single-phase modeling approach took the center stage mainly due to its simplicity, even though the nanofluid heat transfer process is a multiphase problem. Quite many researchers have performed numerical investigations of the natural convection heat transfer of nanofluid using a single-phase model. A non-exhaustive list of studies reported in the open literature comprises; square cavities with straight or wavy walls, cylindrical cavities, cavities of different shapes with a heat source, and so forth. However, for consistency with the scope of the present study, the review here is limited to the differentially heated quadrate enclosure problems. Several characteristics of natural convection in a quadrate enclosure such as; effects of Rayleigh number (or Grashof number), nanoparticles concentration, the enclosure inclination angle, heating methods, cavity size, and varying thermophysical properties have been investigated by many authors using the single-phase model. In this model, the nanofluid is treated as a single-phase fluid but with bulk nanofluid thermophysical properties. Thus, the basic assumption beneath the single-phase approach is that nanoparticles are homogeneously distributed in the base fluid and therefore can be assumed as a single fluid. Khanafer et al. 10 used the single-phase model to study nanofluid heat transfer enhancement in a two-dimensional enclosure considering the solid nanoparticles dispersion. In that pioneering work, the model was used to analyze the impact of the nanoparticles on the fluid flow and heat transfer processes within the enclosure. The paper had a few flaws including the consideration of very high nanoparticles concentration (up to 25%) without accounting neither for agglomeration effects nor for the nanofluid stability (from a practical point of view), model validation using only commercial code data for pure fluid and so on. Nevertheless, this might be understandable by admitting the sparsity of experimental data and the little information available to the research society to work with then. Certainly, the authors were able to identify the fact that the thermophysical properties expressions, used for the nanofluid, significantly affected the results. In the same vein, Abu-Nada and Oztop 11 adopted a single-phase model to study the effects of inclination angle on natural convection in enclosures filled with Cu-water nanofluid. The angle of inclination served as the control parameter for fluid flow and heat transfer processes. Also, Aminossadati and Ghasemi 12 investigated the influences of; Rayleigh number, location and geometry of the heat source, and volume fraction of nanoparticles on the natural convection cooling of a heat source placed at the bottom wall of an enclosure filled with nanofluid. With the same single-phase modeling approach, Oztop and Abu-Nada 13 studied the buoyancy-driven nanofluid heat transfer in a partially filled enclosure using different kinds of nanoparticles where it was found that the location of the heater has a significant impact on the nanofluid heat transfer and the associated fluid flow. More recently, Ibrahim et al. 14 performed a comprehensive study to demonstrate how the loading of nanoparticles could reduce the positive effect on thermal conductivity.
Considering the above, it could be deduced that some of the authors of previous studies were able to qualitatively predict the observed deterioration of the nanofluid heat transfer but the accuracies of their predictions were not verified with experimental data as pointed out by Chen et al. 15 even though the investigated characteristics (Rayleigh number, inclination angle, and so forth) might still be valid in terms of trends. Therefore, it might be premature to conclude that the single-phase model is adequate for the prediction of nanofluid thermal behavior based on these studies. To buttress the doubt concerning the accuracy of the single-phase model's prediction of nanofluid thermal behavior; an experimental study performed by Ho et al. 4 has indicated that the observed variable thermos-physical properties are insufficient to explain the heat transfer behavior of nanofluid at high concentrations.
To transcend the limitations, stated above, of the single-phase model, some researchers have explored alternative models that can capture the actual physics of the nanofluid thermal behavior with better explanations for the deterioration of heat transfer previously observed in experimental studies at high nanoparticle concentrations. For instance, Aminfar and Haghgoo 16 have recently provided an interpretation of this heat transfer degradation by suggesting that the slip motion occurring between the nanoparticles and base fluid could lead to the formation of stagnant thin layers of settled nanoparticles at the bottom adiabatic wall and pure base fluid at the top adiabatic wall. This was claimed to cause a significant deterioration of the heat transfer across the enclosure. Thus, a mechanistic way to capture the nanofluid heat transfer deterioration behavior is to consider the relative motion of nanoparticles to the base fluid. This could result in non-uniform distribution of nanoparticles which in its turn, could have a significant impact on the energy and momentum transfer within the nanofluid. To actualize this, the use of a two-phase model then becomes imperative since the nanoparticles' slip motions cannot be captured with the classical single-phase model. Different researchers have used different two-phase modeling approaches to capture various slip mechanisms they deemed dominant. For instance, a two-phase mixture model that considers nanoparticle slip due to only resistance to flow by particles (drag force) has been used by Chen et al. 15 in which a reasonable agreement with experimental data was obtained. Moreover, a two-phase mixture model that considers sedimentation (i.e., gravitydriven movement) of nanoparticles as the sole slip mechanism has been used by Meng et al. 17 . Unfortunately, the accuracy of the model was not reported as its' prediction was not compared with experimental data. Moreover, based on the strong model put forth by Buongiorno 18 concerning the effect of nanoparticles Brownian diffusion and thermophoretic diffusion, a two-component model accounting for the influence of Brownian diffusion and thermophoretic diffusion has been used by Haddad et al. 19 (although without validation against experimental data) and by Corcione et al. 20 who reported a validation to experimental data within ± 10% error margin where the effect of slip mechanisms in the momentum and energy equations were not considered in their model.  Giddings et al. 22 to substantially underestimate the contribution of thermophoretic diffusion to the nanoparticle migration phenomenon. Therefore, to address the shortcoming of the McNab-Meisen relation, Corcione et al. 20 have attempted to leverage experimental data to develop an empirical correlation for the prediction of the thermophoretic diffusion coefficient. However, their correlation only depends on the nanoparticle concentration and does not take into account the Rayleigh number effects. More so, a detailed review of the control parameters of natural convection in differently shaped cavities with nanofluid recently performed by Rostami et al. 23 has revealed that previous two-phase analysis of nanofluid heat transfer could still not predict accurate data consistent with experimental observations. Besides, a clear-cut impact of the nanoparticle slip mechanisms has not been presented in contrast with the experimental data. Hence, further studies of nanoparticle deposition in terms of model development for the various nanoparticle transport phenomena have been encouraged.
Considering this, the three key slip mechanisms (sedimentation, Brownian diffusion, and thermophoretic diffusion) are incorporated into the two-phase drift flux model. The aim here is to perform an elaborate investigation of the impact of these slip mechanisms on the thermal behavior of nanofluid in a quadrate enclosure. Prior to assessing these nanoparticle slip motions; an evaluation of existing single-phase and two-phase models is performed to ascertain the findings of the literature review. Furthermore, robust empirical correlations for the thermophoretic and sedimentation coefficients are proposed to improve the prediction of the nanoparticle migration due to both; thermophoretic diffusion and sedimentation. The correlations proposed here account, not only, for the nanoparticle concentration but also take into account the Rayleigh number variation effects.

Thermophysical properties of nanofluid
The thermal conductivity of the base fluid gets enhanced by the suspended conductive nanoparticles such as Al 2 O 3 with a diameter of less than 100 nm but this could also cause an undesirable increase in fluid viscosity as explained by Das et al. 24 . The thermophysical properties of the base fluid (water) and alumina nanoparticles which are extracted from the article of Ho et al. 4 are summarized in Table 1. Over the years, several correlations have been suggested to generalize the estimation of the thermophysical properties such as thermal conductivity, specific heat capacity, viscosity, and density. However, some of these correlations were developed from specific experimental data which implies that they may be limited by the corresponding experimental conditions. For instance, Ghanbarpour et al. 25 found that their experimental data of nanofluid thermal conductivity and viscosity were under-predicted by previous correlations of Maxwell 26 and Einstein 27 and this was also corroborated by Albojamal and Vafai 28 . Therefore, due to lack of generalization, care must be taken concerning an acceptable range of validity when selecting appropriate correlations for calculating the nanofluid thermal conductivity and viscosity.
The thermophysical properties of the resulting nanofluid can be derived from the thermophysical properties of the individual components (Table 1). From the literature, two methods are commonly used to calculate these thermophysical properties. The first one assumes that the properties are only dependent on the volume fraction ( ϕ ) of the nanoparticles. The second method is more general in the sense that the physical properties are said to also depend on the temperature of the nanofluid as reported by Nguyen et al. 29 . Typical correlations for calculating temperature-dependent viscosity and thermal conductivity are reported in the work of Abu-Nada 30 , Abu-Nada and Chamkha 31 , Khanafer and Vafai 32 , Bianco et al. 33 , and Palm et al. 34 . It is worth mentioning, however, that dependence of the nanofluids thermophysical properties on the base fluid pH, the surfactant type, and surfactant concentration has been also reported in the literature (see Yoo et al. 35 and Das et al. 36 ). However, to avoid the uncertainty that could result from such temperature and/or stabilizer dependencies, the relations provided in the paper of Ho et al. 4 (reference data in this study) are used and these correlations are summarized in Table 2. The exclusion of the physical properties dependence on temperature is further justified by the fact that a low-temperature difference of 1 K is used in the differential heated enclosure considered herein.

Mathematical model
The natural convection flow considered in this study is presumed to be laminar as the Rayleigh number is less than 7 × 10 6 and three modeling approaches are employed in the simulations. These modeling approaches are; the single-phase model, and two variants of the two-phase modeling approach (the mixture model and the  The Boussinesq gravity term g k in Eq. (2) is computed using the Boussinesq approximation given by Eq. (4) where it is assumed that density variations are small to extent that they have no effects on the flow field except that they give rise to buoyancy force. Finally, the nanofluid is assumed to be incompressible.
Here u , t , T , P , and g represent the velocity vector, time, temperature pressure, and gravity vector respectively. Additionally, the physical properties of the nanofluid: dynamic viscosity, density, specific heat capacity, thermal conductivity, and volumetric expansion coefficient are represented by µ m , ρ m , C pm , k m , and β m and are defined in Table 2.
Mixture model governing equations. Arguably, the non-uniformity of nanoparticle concentration in the quadrate enclosure disqualifies the use of a single-phase model for the prediction of the nanofluid heat transfer. This is since the basic single-phase model assumption (homogenous distribution of nanoparticles) no longer holds. Thus, one way to investigate the potential effect of nanoparticles' non-uniformity is to employ a mixture modeling approach. In mixture modeling approaches, aside from non-uniform nanoparticle distribution, other assumptions of the single-phase model still hold. In addition, the slip motion of nanoparticles caused by the gravity force is considered while thermophysical properties of the nanofluid are evaluated based on the local concentration of the nanoparticles.
The gravity force is considered in the mixture model based on the previous experimental observation by Chang et al. 3 where they reported that the sedimentation of nanoparticles could be responsible for the deterioration of nanofluid heat transfer as concentration increases in the enclosure. As a developmental basis, the existing drift flux solver "driftFluxFoam" in OpenFOAM 6 is used to capture the described mixture modeling concept. In addition to the variable thermo-physical properties, only the energy equation is added to the existing drift flux solver to capture the thermal characteristics of the nanofluid. In this solver, the mixture flow is described by the following governing equations of mass, momentum, nanoparticle volume fraction, and energy in Eqs. (5), (6), (7), and (8), respectively.
Specific heat capacity where u m , u pm , u bfm , ρ p , ρ bf , T m represent the mixture velocity vector, the relative velocity vector of nanoparticles, the relative velocity vector of the base fluid, the density of nanoparticles, the density of the base fluid, and the mixture temperature, respectively. The previous definitions of other symbols still hold. As stated previously, the sedimentation or settling of nanoparticles causes the slip motion of nanoparticles. Therefore, the relative velocity of nanoparticles and base fluid resulting from the settling of nanoparticles is assumed to follow the Vesilind 37 sedimentation model given by Eqs. (9) and (10). A schematic sketch representing the migration of nanoparticles in the gravity direction to settle at the bottom of the enclosure is illustrated in Fig. 1.
where u s and A represent the reference settling velocity vector and settling coefficient and can be determined from experimental data. In absence of dedicated experimental data, the reference nanoparticle settling velocity is calculated from the balance of the buoyancy and viscous force as shown in Eq. (11).

Two-component model governing equations.
The two-component model is a variant of the twophase modeling approach that was proposed by Buongiorno 18 . This model is predicated on the assumption that nanoparticle slip velocities resulting from various slip mechanisms are responsible for the convective heat transfer enhancement in the nanofluid. The two dominant slip mechanisms, identified by Buongiorno 18 are; the Brownian diffusion and the thermophoresis diffusion. The random motion of nanoparticles dispersed in a base fluid is termed Brownian diffusion resulting from the continuous collision between nanoparticles and the base fluid molecules while thermophoresis refers to the migration of nanoparticles under the influence of temperature gradient. Schematic sketches illustrating the thermophoresis and the Brownian diffusions of nanoparticles are shown in Fig. 2a,b, respectively. A detailed description of these two slip mechanisms can be found in the textbook of Michaelides 38 and the review paper of Kleinstreuer and Xu 39 . Accounting for these two dominant slip mechanisms, a two-component mixture model (base fluid and nanoparticles) has been formulated by Buongiorno as shown in Eqs. (12)(13)(14)(15) for mass, momentum, energy, and nanoparticle volume fraction conservation equations, respectively. Starting from the buoyancy-driven convection  where J p and C pp represent total nanoparticle mass flux and nanoparticle specific heat capacity, respectively. The total nanoparticle mass flux accounting for the two slip mechanisms, of Brownian diffusion and thermophoresis, is defined by Eq. (16) below.
where D B = k B T m 3πµ bf d p and D T = β µ m ρ m ϕ are the Brownian diffusion coefficient and thermophoretic diffusion coefficient, respectively. Therefore, the slip velocity due to Brownian diffusion is obtained as u B = −D B ∇ϕ while the slip velocity due to the thermophoretic diffusion is obtained as Additionally, the thermophoretic coefficient ( S T ) is determined using the expression given by McNab and Meisen 21 shown in Eq. (17). It is worth mentioning here that the McNab and Meisen expression is independent of the particle size and was originally developed using data of large size particles. Thus, it cannot be scaled down to nanoscale in its present form and might be inappropriate for the prediction of thermophoretic coefficient in a nanofluid.
here k p is the nanoparticle thermal conductivity. To address the inappropriateness of the McNab and Meisen 21 model for nanofluid, a correlation was derived by Corcione et al. 20 (shown in Eq. (18)) that described the generated dataset of S T which minimized deviation of the errors between simulation data and experimental data of the nanofluid heat transfer. Although the correlation proposed by Corcione et al. 20 did not account for the effect of the Rayleigh number, it is tested in this study to assess its accuracy in the prediction of the nanoparticle thermophoretic diffusion.  18). Unfortunately, none of these values is strong enough to significantly influence the nanofluid heat transfer as it will be shown further down. However, Eq. (17) was developed more mechanistically for large particles. Therefore, it would be more logical to extend this thermophoretic coefficient equation (Eq. (17)) to nanoparticle scale.

Geometric configuration, boundary conditions, and numerical solution scheme
In this study, a two-dimensional (2D) differentially heated enclosure representing the experimental configuration of Ho et al. 4 is considered. The quadrate enclosure has equal width and height of 40 mm. A constant temperature, T H and T C boundary conditions are imposed on the right hot side, and left cold side, respectively. Additionally, adiabatic boundary conditions are imposed on the top and bottom sides while non-slip velocity boundary condition is imposed on all four sides. For calculation involving the volume fraction equation, zero particle flux boundary conditions (Neumann boundary condition) in place of constant value boundary conditions (Dirichlet boundary condition) are imposed on all four sides. These boundary conditions are depicted in Fig. 3. As can be seen in Fig. 3, a non-uniform grid is generated with fine mesh near the wall to capture the near-wall behavior of nanofluid and relatively coarse mesh in the central region, where the velocity and temperature gradients are minimal. The nanofluid is made of alumina nanoparticles (Al 2 O 3 ) of average particle size of 33 nm dispersed in ultra-pure water (serving as the base fluid) at different volume fractions ranging from 0.1 to 3%.
The two-dimensional governing equations of the modeling approaches are discretized using the finite volume method. Then a PIMPLE algorithm in OpenFOAM is used, which is a combination of the pressure implicit with the splitting of operator (PISO) algorithm and semi-implicit method for pressure linked equations (SIMPLE) algorithm. More detailed descriptions of these algorithms can be found in the work of Issa 40  The time step is therefore directly proportional to the mesh spacing ( x i ). If the spacing is reduced by half, then the time step must be also reduced by half. To assess the solver stability criterion stated above (Eq. (19)), www.nature.com/scientificreports/ the maximum Courant number was changed to 2, 5, and 10 to ensure that the presently developed solver is sufficiently robust and that the time discretization does not affect the final steady-state results. As it can be seen in Fig. 5a, only the transient part (up to the physical time of ~ 1000 s) of the simulation is affected when increasing the Courant number maximum value limit to 2, 5, and 10. All the Courant numbers produce a converged temperature profile at the mid-point of the computational domain to an asymptotic value. Moreover, the disparity observed using different Courant numbers results from setting a maximum number for the PIMPLE iterations equal to 3. This sensitivity to time-step is eliminated by setting instead a convergence criterion of 10 -6 for each   Fig. 5b, irrespective of the maximum Courant number, the temperature profiles at the midpoint of the domain converge throughout the simulation physical time. This way, the PIMPLE iteration continues until these convergence criteria are achieved. The time-step discretization sensitivity shown in Fig. 5 is for the mixture model. However, this behavior is also true for both the single-phase model and two-component model. Based on the above findings, a maximum Courant Number is kept to 0.5 (as this had no dramatic impact on the simulation time) and the absolute residuals of velocity, temperature, pressure, and nanoparticle volume fraction are restricted to below 10 -6 as a convergence criterion during the iteration process in the present study. At the attainment of the asymptotic solution, the heat flux at the hot boundary surface is equal to the one at the cold boundary surface which is calculated using Eq. (20).
where n denotes the normal to the surface. Subsequently, the Rayleigh number, the heat transfer coefficient, and the Nusselt number are computed using Eqs.  www.nature.com/scientificreports/ where q " is the spatial average of the heat flux on the hot boundary surface (which is the same for the cold boundary surface), T is the temperature difference between the cold surface ( T C ) and hot surface ( T H ).
Using the single-phase and the two variants of the two-phase model described earlier in "Mathematical model", a grid independence test is performed on six grids 60 × 60, 80 × 80, 100 × 100, 120 × 120, 130 × 130, and 140 × 140 as shown in Table 3. The predicted average Nusselt number at the hot side for Ra = 6 × 10 6 (the maximum Rayleigh number considered in this study) at a nanoparticles volume fraction of 3% shows that variation with grid becomes infinitesimal after the grid size of 100 × 100 as there is less than 0.1% difference in the Nusselt number between the 100 × 100 and 120 × 120 grid size. Therefore, subsequent computations are performed using 100 × 100 grid size. Table 3 further shows a grid-independence study for the three models.

Validation of the models and discussion of results
In this section, the three modeling approaches are first validated against the experimental data of Ho et al. 4 for the ultra-pure water cases. Hence, simulations have been conducted for laminar natural convection in the square cavity (see Fig. 3). This elementary yet compulsory step is required to evaluate the accuracy and reliability of the models' implementation in the open-source code. As can be seen in Fig. 6, all three modeling approaches are able to reproduce the experimentally obtained average Nusselt numbers, reported by Ho et al., within the uncertainty range. Furthermore, as to be expected, the predictions from the three approaches converge to a single curve since the working fluid does not contain nanoparticles. That is, in the absence of suspended nanoparticles, the solution www.nature.com/scientificreports/ reverse back to the pure water predictions. This implies that all terms containing the void fraction variable (φ), including the ones in the thermophysical properties correlations, are correctly eliminated for the solver to reduce to its basic form; valid for pure water. Added in the figure are the experimental data of Holland et al. 42 , and Churchill and Chu 43 which are represented by the correlations given in Eqs. (24) and (25), respectively. At first sight, the discrepancy observed between the values obtained using these two correlations and the experimental measurements of Ho et al. might be of concern. However, by recognizing the fact that these correlations were generated from many data, spanning over a wide range of Rayleigh numbers, and extending from the laminar to the turbulent flow regimes then, such differences can be justified.
Also, for the cases with nanoparticles concentration of 1%, all three models can reasonably predict the experimental measurements within the ± 7% error band (see Fig. 7). It must be stated that the Brownian diffusion and Thermophoretic diffusion (computed using the correlation proposed by Corcione et al. 20 for the thermophoretic parameter) has no significant impact on the nanofluid heat transfer rate, as the resulting deviation of the twocomponent model from the single-phase model is no more than ~ 0.2%. A second interesting deduction to be made from these plots is that there is no apparent added value by using the more complex two-component model instead of the basic one for such a low concentration. As for the third nanoparticle slip mechanism (sedimentation), the deviation of the mixture model from the single-phase model is ~ 3%. Although this model is also able to predict the experimental data within the uncertainty limit as can be seen in Fig. 7, the model returns somewhat underpredicted values for the Nusselt number in comparison to the other two models and the experimental data. This is a twofold observation, firstly accounting for sedimentation does result in heat transfer impairment and secondly, these underpredicted values for the Nusselt number give a hint to the fact that the nanoparticle sedimentation term might be requiring some adjustment to better reflect the reference experimental data. These last two points are further investigated and discussed in "Proposed combination of the nanoparticle's thermophoresis diffusion, Brownian diffusion, and sedimentation in a two-phase model" below.
As shown in Fig. 8, once the volume fraction of the nanoparticles is increased to 3%, clear deviations are observed between the predicted data of the single-phase model and the mixture model. Comparing the predictions obtained by these two models with the experimental data, very intriguing observations can be made. Firstly, in the low Rayleigh number range (up to 2.3 × 10 6 approximately), the mixture model predictions are much closer to the experimental data than the ones obtained with the single-phase model. Then, for the cases with higher Rayleigh numbers, the single-phase model predictions become superior as the mixture model is observed to significantly underpredict the experimental data. This is suggesting that the current sedimentation term, in the mixture model, is providing the right magnitude for the slip mechanism at the low Rayleigh number range, but the same cannot be said for Rayleigh number values larger than 2.3 × 10 6 . Hence, by reference to the Eqs. (6), (9), (10), and (11), the sedimentation term should also account for the Rayleigh number effects in addition to www.nature.com/scientificreports/ the nanoparticle concentration ones. This finding is in agreement with the common understanding that at slow flow motion (i.e., low Rayleigh number) the nanosuspensions are more prone to sedimentation, but by increasing the flow motion this effect becomes less significant and the nanoparticles get carried away by the flow stream. Also, by analyzing the two-component model predictions, it seems that the thermophoretic diffusion and Brownian diffusion terms are not strong enough to influence the heat transfer rate. Actually, the differences in the predictions by the two-component model from the ones by the single-phase model are hardly visible. Here again, the upgrade of the modeling approach, by accounting for these two slip mechanisms and solving an additional equation for the nanoparticles volume faction, does not seem to provide any apparent benefits. However, the same cannot be said for the field distribution of the nanoparticle volume fractions illustrated in Fig. 9. In this figure is not easy to distinguish between the velocity and temperature fields predicted by the three models but the difference in the nanoparticle volume fraction is obvious. As can be seen in Fig. 9c, both the two-component model and mixture model are able to capture the nanoparticle non-uniformity in the concentration. The apparent drift of the nanoparticles from the hot to the cold regions in the two-component model results is due to the inclusion of the thermophoretic diffusion term in this model, while the impact of the Brownian diffusion term is less noticeable. As for the other two models; the nanoparticles are uniformly distributed in a single-phase model (see Fig. 9a) conforming with its basic assumption, while the nanoparticles sedimentation is captured by the mixture model as can be seen in Fig. 9b. As intended, the mixture model is predicting the nanoparticles settling at the bottom of the enclosure due to the sedimentation term effects. It is interesting to note that these high/low nanoparticle concentrations are observed to be taking place locally at the vicinity of the top/bottom isothermal walls, yet the resulting effect is propagated to vertical hot and cold walls. Furthermore, by comparing the temperature and velocity fields in Fig. 9a,b, these effects do not seem to have a significant impact on the core of the natural convection rotational motion.
To closely assess the three modeling approaches' ability to predict the heat transfer impairment due to the nanoparticles increased concentration in the nanofluid, the predictions of the models for the Nusselt number (normalized by the corresponding pure water case) as a function of nanoparticle concentration is plotted in Fig. 10. Interestingly, in conformity with the experimental data, the three modeling approaches predict the deterioration of nanofluid heat transfer as the nanoparticle concentration increases. Quantitatively, for the case of ϕ = 3% , the average values of the predicted deterioration by the mixture, single-phase, and two-component models are 19.9%, 9.7%, and 9.5%, respectively. Comparing these values with the corresponding deterioration in Ho et al. 4 experimental data which is 18.3%, it can be said that the mixture model captures better the heat transfer deterioration, while the other two models are of equal performance.
In light of the above analysis, the single-phase model is confirmed to be inadequate for a quantitative prediction of nanofluid natural convection heat transfer at high nanoparticles concentration (e.g., at ϕ = 3% ), especially at low Rayleigh numbers, where sedimentation effects can no longer be neglected. The two-phase mixture model, on the other hand, is over predicting the heat transfer deterioration for nanofluid with high nanoparticles concentration, but only for the high Rayleigh number cases. This model seems to be also missing physical flow phenomena associated with the thermophoretic diffusion and, to some extent, the Brownian motion slip mechanisms. Lastly, the two-component model, on its current version, does not seem to provide significant benefits to justify the additional modeling complexity. From this model results, the dominant slip mechanism is the thermophoretic diffusion, while from the mixture model results the sedimentation slip mechanism seems www.nature.com/scientificreports/ to be also playing an important role. Hence, to further assess the impact that these two terms (sedimentation and thermophoretic diffusion) might have on the heat transfer rate, it is possible to adjust them through their controlling parameters (settling coefficient (A) and thermophoretic coefficient ( S T )) as detailed in "Proposed combination of the nanoparticle's thermophoresis diffusion, Brownian diffusion, and sedimentation in a twophase model" below.

Proposed combination of the nanoparticle's thermophoresis diffusion, Brownian diffusion, and sedimentation in a two-phase model
Observably, in the two-component model (see Fig. 11b), the slip velocity due to the thermophoretic diffusion is comparable (in terms of magnitude) with the slip velocity (shown in Fig. 11a) resulting from nanoparticle sedimentation in the mixture model. In agreement with the observation made above, the slip velocity resulting from the thermophoretic diffusion is far higher than that of Brownian diffusion as can be seen in Fig. 11b. As nanoparticles migrate from the hot region to the cold region due to the thermophoresis diffusion causing nonuniformity in the particle distribution, the Brownian diffusion is supposed to move the nanoparticles in the opposite direction of the particle concentration gradient to create a homogeneous distribution of the nanoparticles. However, by comparing the magnitude of these slip velocities, it appears that the Brownian diffusion does not have any significant contribution to the nanoparticle migration in the nanofluid. Therefore, it can be argued that the two dominant slip mechanisms are the thermophoretic diffusion and sedimentation of nanoparticles. Thus, there is a need for a two-phase formulation able to account for both; the thermophoretic diffusion and sedimentation of nanoparticles. However, as noted above, adjustment for these two terms, through their respective coefficients A and S T parameters, is required. Previously, Esfandiary et al. 44 , and Pakravan and Yaghoubi 45 have attempted to account for nanoparticle slip effect due to Brownian and thermophoretic diffusions (without sedimentation) using a different formulation from Buongiorno 18   www.nature.com/scientificreports/ previous definitions of all parameters still hold. In this formulation, the nanoparticle slip velocities represented by u pm , u B and u T correspond to sedimentation, Brownian diffusion, and thermophoretic diffusions, respectively. These parameters have been previously defined in "Mathematical model" of this article. Using the OpenFOAM mixture model solver, "driftFluxFoam" as a developmental basis, a new solver is developed to implement the proposed model.
Using the default settling coefficient A = 285.84 and thermophoretic coefficient S T = 0.01 (between the value predicted by Eqs. (17) and (18)) as the baseline parameters, Fig. 12 shows that the prediction by the proposed model is now clearly different from the single-phase model but almost the same as the mixture model since the thermophoretic diffusion is weak at this value of the thermophoretic coefficient.
With this solution framework, the two coefficients ( A and S T ) can be adjusted to increase the strength of their respective slip mechanisms in nanofluids to investigate their impact on the nanofluid natural convection heat transfer. At first, keeping the settling coefficient constant ( A = 285.84 ), different computations are performed with a different thermophoretic coefficient ( S T ) values of 0.01, 2, 5, and the impact of these values on the distribution of the nanoparticle concentration are shown in Fig. 13a,b, Fig. 13d,e,f, respectively. As the settling coefficient increases, more nanoparticles are settled at the bottom of the enclosure while creating a thin layer of purer fluid at the top. Additionally, the nanofluid velocity profile and temperature profile are also affected by these slip mechanisms as can be seen in Fig. 14a,b, respectively. Perhaps the most important observation to note here is   www.nature.com/scientificreports/ Figure 15 shows the nanoparticle volume fraction distribution and temperature distribution along horizontal and vertical mid-plane. These results indicate that nanoparticle volume fraction is lesser at the hot wall where the temperature is higher due to the thermophoretic diffusion of nanoparticles from the hot wall to the cold wall as shown in Fig. 15a. Likewise, since the top region is hotter than the bottom region as can be seen in Fig. 15b, thermophoretic diffusion also contributes to the downward migration of the nanoparticles from the top to the bottom of the cavity which subsequently results in nanoparticle deposition on the bottom wall. Coincidentally, the gravity acts in the downward direction as well, thereby causing nanoparticle sedimentation on the bottom in addition to the thermophoretic diffusion of nanoparticles. This implies that the migration of nanoparticles from the hot wall to the cold wall through thermophoretic diffusion favor and enhance heat transfer across the enclosure while the deposition of nanoparticles at the bottom wall causes the formation of a stagnant thin layer at the top and bottom walls which could lead to the deterioration of the nanofluid heat transfer.
To gain further insight into the nanoparticle transport in nanofluid flow, the terms in the nanoparticle volume fraction equation corresponding to the contribution of the three nanoparticle transport phenomena are examined as shown in Fig. 16. As earlier stated, the contribution of the Brownian diffusion to the migration of the nanoparticle is relatively insignificant as can be seen in Fig. 16. Moreover, the thermophoretic diffusion is stronger in the horizontal direction near the hot and cold walls than near the top and bottom walls as shown in Fig. 16a,b. This is because the temperature gradients near the hot and cold walls are steeper than the temperature gradient near the top and bottom walls (see Fig. 15). Additionally, only thermophoretic diffusion contributes to the horizontal migration of nanoparticles as shown in Fig. 16a whereas both thermophoretic diffusion and sedimentation by gravity are responsible for the movement of nanoparticles to the bottom and top of the cavity as shown by the inset in Fig. 16b. www.nature.com/scientificreports/ The impact of the non-uniformity of nanoparticle distribution caused by the various transport mechanisms on the local nanofluid heat transfer represented by the local Nusselt number is shown in Fig. 17. As the thermophoretic diffusion increases, the local Nusselt number curve of the hot and cold walls is shifted upward indicating enhancement of the nanofluid heat transfer as can be seen in Fig. 17a. On the contrary, the increase in sedimentation of nanoparticles due to gravity causes the local Nusselt number to decrease near the top and bottom wall as shown in Fig. 17b. This is because the effective conductivity of the nanofluid is reduced at the top while the effective viscosity increased at the bottom. This is due to the sharp reduction of nanoparticle concentration at the top (almost pure base fluid) and a large increase of nanoparticle concentration at the bottom. The corresponding effect of this variation of nanoparticle concentration at the top and bottom could be seen from the empirical correlations of the effective thermal conductivity and viscosity earlier shown in Table 2.
After ascertaining that the two dominant slip mechanisms have been properly captured by the proposed twophase model, the impact of these slip motions on the nanofluid natural convection heat transfer is investigated next. As shown in Fig. 18a, as the thermophoretic diffusion increases while sedimentation of nanoparticles is suppressed, the nanofluid heat transfer is enhanced. This is so because the energy transport by the migrating nanoparticles from hot regions to cold regions generally favors the essence of natural convection heat transfer. Although thermophoresis also causes nanoparticle deposition on the bottom wall, the migration of nanoparticles from hot wall to cold wall is far stronger as previously shown in Fig. 16. This can also be corroborated by the work of Buongiorno 18 in which thermophoretic diffusion alongside Brownian diffusion has been used to explain the abnormal enhancement of nanofluid heat transfer that cannot be captured by dynamic thermophysical properties. Conversely, as the sedimentation of nanoparticles increases, while the thermophoretic diffusion is suppressed, the nanofluid heat transfer deteriorates as can be seen in Fig. 18b. The reason for this has been previously discussed and illustrated in Fig. 17. Moreover, the sedimentation of nanoparticles is more strongly dependent on the size of the nanoparticles than the settling coefficient as shown in Eq. (11). Therefore, it can be said that the thermophoretic diffusion and sedimentation of nanoparticles have opposing effects on the nanofluid heat  Having observed the impact of both thermophoretic diffusion and sedimentation of particles, the two parameters ( A and S T ) appear to depend on nanoparticle concentration and Rayleigh number. Therefore, a generalized expression for the two parameters in terms of ϕ and Ra are formulated and are written as shown in Eqs. (30) and (31).
where a 1 , a 2 , a 3 , b 1 , b 2 , and b 3 are constants that can be determined through multiple regressions of data for A and S T obtained using experimental measurements of nanofluid heat transfer as reference data. Using Ho et al. 4 experimental data set for nanofluid concentrations of 1% and 3%, and Rayleigh number ranging from 1 × 10 6 to 1 × 10 7 , the values of the unknown parameters A and S T that minimize the standard deviation of errors between the predicted heat transfer data and experimental heat transfer data are determined. The generated values of these parameters are then used to determine the constants in Eqs. (30) and (31) through multiple regression analysis. Thus, the values of a 1 , a 2 , a 3 , b 1 , b 2 , and b 3 computed are shown in Table 4.
Deploying Eqs. (30) and (31) in the proposed two-phase model yields a much better prediction of experimental data as shown in Fig. 19. However, further calibration of these two important parameters using a wider range of experimental data is required for a more accurate prediction of nanofluid heat transfer at high nanoparticle (30) A = a 1 + a 2 ln ϕ avg + a 3 ln(Ra)   Fig. 20, the associated computational cost of using this new two-phase model over the single-phase model is acceptable especially, considering the derivable benefits in terms of its accuracy.

Conclusion
This work proposes an enhanced two-phase drift-flux nanofluid model, by accounting for the slip nanoparticles mechanisms: Brownian and thermophoresis diffusions and more importantly the nanoparticle sedimentation. Numerical experiments of buoyancy-driven nanofluid in a differentially heated square cavity showed that the proposed model can predict the heat transfer rates for a variety of nanoparticle void fractions (from 0 to 3%) and a wide range of Rayleigh numbers. It is worth noticing that assessments of existing single-phase and twophase models against experimental data performed in this study have revealed the inadequacies of these models in obtaining the correct heat transfer rate magnitude. Their respective outcomes are seriously affected by the nanoparticle void fraction and Rayleigh number. Above and beyond, as the thermophoresis diffusion enhances the heat transfer, the sedimentation of nanoparticles deteriorates it. The Brownian motion term, on the other hand, is found to have a negligibly small www.nature.com/scientificreports/ contribution. Considering this, an empirical correlation is proposed for the controlling parameters appearing in the thermophoresis diffusion and sedimentation terms (A, S T ) and demonstrated better predictions of nanofluid heat transfer. Furthermore, the associated additional cost required to achieve this numerical prediction improvement remains within the acceptable range for industrial applications. Future work, planned by the team, will consist in further calibration of these empirical correlations using a wider range of experimental data and flow conditions.