Homogeneous and Multiphase Analysis of Nanofluids Containing Nonspherical MWCNT and GNP Nanoparticles Considering the Influence of Interfacial Layering

Abstract

The practical implication of nanofluids is essentially dependent on their accurate modelling, particularly in comparison with the high cost of experimental investigations, yet the accuracy of different computational approaches to simulate nanofluids remains controversial to this day. Therefore, the present study is aimed at analysing the homogenous, multiphase Eulerian–Eulerian (volume of fluid, mixture, Eulerian) and Lagrangian–Eulerian approximation of nanofluids containing nonspherical nanoparticles. The heat transfer and pressure drop characteristics of the multiwalled carbon nanotubes (MWCNT)-based and multiwalled carbon nanotubes/graphene nanoplatelets (MWCNT/GNP)-based nanofluids are computed by incorporating the influence of several physical mechanisms, including interfacial nanolayering. The accuracy of tested computational approaches is evaluated by considering particle concentration and Reynolds number ranges of 0.075–0.25 wt% and 200–470, respectively. The results demonstrate that for all nanofluid combinations and operational conditions, the Lagrangian–Eulerian approximation provides the most accurate convective heat transfer coefficient values with a maximum deviation of 5.34% for 0.25 wt% of MWCNT–water nanofluid at the largest Reynolds number, while single-phase and Eulerian–Eulerian multiphase models accurately estimate the thermal fields of the diluted nanofluids at low Reynolds numbers, but overestimate the results for denser nanofluids at high Reynolds numbers.

1. Introduction

Nanofluids are nanotechnology-based novel fluids of ultrahigh thermal efficiency. Over the past several decades, these thermal fluids have gained overwhelming attention in relevance to their potential applications in microelectronics, transportation, refrigeration, solar, nuclear and space technologies [1,2]. Even at a low concentration, nanoparticles possess high effective surface area and therefore exhibit exceptional thermal efficiency, flexible thermophysical properties, suspension stability and enlarged solid–particle interface for maximum interphase heat exchange. Additionally, the reduced tendency of sedimentation, channel clogging and negligible pressure drop make nanofluids appealing for the advanced thermal engineering systems [3,4,5,6]. The hydrothermal performance of the nanofluids is characterised by nanoparticle loading, morphology, hosting fluid thermophysical properties and operational temperature [7,8,9,10]. The assessment of hydrothermal characteristics of nanofluids for practical applications is essentially dependent on their accurate modelling, particularly in comparison with high-cost experimental investigations. The accurate simulation of nanofluids requires an inclusive understanding of the underlying physical phenomena responsible for the enhanced thermal activity of these wonder fluids [11,12]. Nanoparticle Brownian motion, drag and lift forces, thermophoresis, thermal boundary layer disruption and liquid molecular layering around nanoparticles are several of the extensively accepted potential mechanisms [13,14,15,16].

Among the different computational techniques, the Navier–Stokes-based homogeneous and multiphase approaches are the most commonly adopted ones to simulate nanofluid flow and heat transfer characteristics. In the single-phase approach, nanofluids are presumed as a homogenous mixture of nanoparticles and hosting fluid with negligible slip in between. This homogenous solution is assumed to exhibit enriched thermophysical properties that are assessed using experimental observations or theoretical models. Unlike the single-phase model, the nanoparticles and the hosting fluids are modelled as distinct phases in the multiphase numerical approaches and their mutual interactive forces are computed. The multiphase approaches are classified as Eulerian–Eulerian and Lagrangian–Eulerian based on treatment of nanoparticles and liquid phases. In the Lagrangian–Eulerian approach, the nanoparticles are treated as discrete entities and the governing equations of the liquid phase are simulated in the continuum approximation. The Eulerian–Eulerian approximation is based on the assumption of continuous interpenetration of the nanoparticles and the liquid phases. The latter approach comprises three numerical models of the volume of fluid (VOF), mixture and the Eulerian models. In the VOF model, the two phases are approximated as immiscible regimes and their interphase is traced. The mixture model considers that individual phases of the solution comprise their slip velocities and volume concentration fields. However, a single set of governing equations is solved for all of the phases, whereas the transport equations of each phase are solved independently in the Eulerian model and the coupling between all phases is achieved with the pressure and interphase exchange coefficients.

Concerning the accuracy of homogeneous and multiphase approaches to compute nanofluids, so far, the literature contains immense controversies and debate about the selection of an accurate model. The homogeneous or single-phase approximation of nanofluid, credited with simplicity and cost-effectiveness, is the most extensively utilised technique [17,18,19,20,21,22,23]. Several comparative reports also emphasised the accuracy of this approach [24,25,26,27,28]. However, a group of researchers condemned homogenous approximation of highly concentrated nanofluids due to the significant influence of interphase interactive forces that are neglected in this approach [2,29,30,31,32,33,34,35,36,37]. Analysing forced convection of Al₂O₃/H₂O, TiO₂/H₂O and Cu/H₂O nanofluids, Nishat et al. [38] demonstrated that the homogeneous model accurately predicts the highly concentrated nanofluids (0.5% ≤ φ ≤ 2%) compared with the Lagrangian–Eulerian model. Similar contradictory statements are witnessed related to multiphase numerical modelling of the nanofluids. Majid et al. [39] evaluated the Eulerian and mixture models to assess the hydrothermal characteristics of the Al₂O₃/H₂O (0.1% ≤ φ ≤ 2%) nanofluid and demonstrated that the analogous results estimated by both of the models were consistent with the experimental data. However, in an assessment of homogeneous, multiphase mixture and Eulerian models for simulating hydrothermal aspects of the Al₂O₃/H₂O (0.6% ≤ φ ≤ 1.6%) nanofluid, Göktepe et al. [40] criticised the mixture model for underperforming compared with the Eulerian model in terms of accuracy. On the contrary, Lotfi et al. [41], Ehsan et al. [37] and Behroyan et al. [28] critiqued the capability of the Eulerian model to simulate the thermal behaviours of the nanofluids. Lotfi et al. [41] and Behroyan et al. [28] stated an underestimation of the heat transfer enhancement of Al₂O₃/H₂O (2% ≤ φ ≤ 7%) and Cu/H₂O (1% ≤ φ ≤ 2%) nanofluids, respectively, whereas Ehsan et al. [37] overestimated the convective heat transfer coefficient results of TiO₂/H₂O (1% ≤ φ ≤ 2.3%) nanofluid when computed utilising the Eulerian approach. In a comparative evaluation of single-phase and Eulerian–Eulerian models, Akbari et al. [25,31] stated that Eulerian, mixture and VOF models overpredict the convective heat transfer coefficient of highly concentrated Al₂O₃/H₂O (0.6% ≤ φ ≤ 1.6%) nanofluid. For the Lagrangian–Eulerian approach, several reports endorsed its superiority [32,33,34,36,42], but a few conflicting reports claimed otherwise [26]. The possible justifications of these contradictions are the insufficient nanofluid-related numerical investigations, lack of understanding of physical phenomena responsible for the enhanced thermal activity of nanofluids and inaccurate modelling. Another major influencing factor is the anomalies in the experimental data used for the validation of different numerical models. Considering the fact that thermophysical properties of nanofluids are immensely affected by the nanofluids’ stability, nanoparticle agglomeration or sedimentation in unstable nanofluid samples may mislead the evaluated results.

In view of the above controversies, a recent review emphasised conducting systematic investigations on the single- and multiphase modelling of nanofluids by incorporating all of the correlated physical mechanisms to establish fundamental results [43,44]. The report also accentuated evaluating the accuracy of different computational models to predict the hydrothermal characteristics of the nonspherical nanoparticles, given that all of the previous relevant studies are based on spherical approximation of nanoparticles. Therefore, addressing this research gap, in the present study, we investigate the accuracy of homogenous, Eulerian–Eulerian (VOF, mixture, Eulerian) and Lagrangian–Eulerian approaches for simulating the hydrothermal characteristics of the MWCNT–H₂O and MWCNT/GNP–H₂O hybrid nanofluids by incorporating several physical mechanisms, including interfacial nanolayering. The influence of nanolayering around nanoparticles is examined by introducing an effective volume fraction of nanoparticles that defines the combined effects of nanoparticle volume fraction and surfactants on the rheological properties of the solution.

2. Problem Description

The forced convection of aqueous-based nanofluids containing suspensions of nonspherical MWCNT and MWCNT/GNP hybrid nanoparticles is simulated in a horizontal minichannel subjected to a constant heat flux of 11,800 W/m² (Figure 1). The probed range of the Reynolds number is 200 ≤ Re ≤ 470 and the flow is presumed as developing and fully developed. The considered nanoparticle weight concentrations are 0.075% ≤ ψ ≤ 0.25% for the MWCNT-based nanofluids, whereas each sample contained an additional ψ = 0.035% of GNP nanoparticles for the case of MWCNT/GNP hybrid nanofluids. The MWCNT and GNP nanoparticles are assumed to be dispersed homogenously into the aqueous solution without any surfactant. The MWCNT nanoparticles are approximated as tubular-shaped with an average length and diameter of 3 um and 15 nm, respectively. The averaged diameter and thickness of disc-shaped GNP nanoparticles are 7 nm and 15 um, respectively. The minichannel hydraulic diameter and length are D_h = 0.0011 m and L = 0.27 m, respectively. The accuracy of the computed results is assessed in comparison with the experimental conditions specified by Hussein et al. [45] due to the availability of rheological properties of the nanofluid. All the specifics of the physical model and nanofluid combinations are kept in accordance with the experimental study of Hussein et al. [45] to perform a comprehensive comparison of the tested numerical models.

3. Mathematical Formulation

3.1. Single-Phase Model

The single-phase modelling of nanofluids is based on the assumption that the solution is a single homogenous fluid of enhanced thermophysical properties. The nanoparticle and dispersion phases are assumed to be thermally and hydrodynamically at equilibrium with each other and, consequently, exhibit negligible interphase momentum and energy exchange. The accuracy of this approach is predominantly dependent on the accurate prediction of the thermophysical properties of the nanofluids, which are therefore presently estimated using the same property correlations utilised in the respective experimental study [45]. The single-phase model’s governing equations are

$$\nabla \cdot \left({\rho }_{nf}\overrightarrow{V}\right)=0$$

(1)

$$\nabla \left({\rho }_{nf}\overrightarrow{V}\overrightarrow{V}\right)=-\nabla P+\nabla ·\left({\mu }_{nf}\nabla \overrightarrow{V}\right)+{\rho }_{nf}g$$

(2)

$$\nabla ·\left({\rho }_{nf}\overrightarrow{V}H\right)=\nabla ·\left(k_{nf}\nabla T\right)$$

(3)

The density (${\rho }_{nf}$) and specific heat (${(c_p)}_{nf}$) of the nanofluids are estimated by Equations (4) [14] and (5) [46], respectively, while Equations (6) [47] and (7) [48] define the temperature-dependent viscosity (${\mu }_{nf}$) and thermal conductivity ($k_{nf}$) of the nanofluids, respectively.

$${\rho }_{nf}=\phi {\rho }_s+\left(1-\phi \right){\rho }_f$$

(4)

$${(c_p\rho )}_{nf}=\phi {\left(c_p\rho \right)}_s+\left(1-\phi \right){\left(c_p\rho \right)}_f$$

(5)

$${\mu }_{nf}=\left(1+2.5\phi +6.2{\phi }^2\right){\mu }_f\left(T\right)$$

(6)

$$k_{nf}=\left(\frac{k_s+(n-1)k_f\left(T\right)+\left(n-1\right)\left(k_s-k_f\left(T\right)\right)\phi }{k_s+(n-1)k_f\left(T\right)-\left(k_s-k_f\left(T\right)\right)\phi }\right)k_f\left(T\right)$$

(7)

Here, n = 3/ϖ is the shape factor of the nanoparticle and ϖ denotes the sphericity of the particles. Sphericity is the ratio of the surface area of a sphere with the volume equivalent to the average particle to the surface area of the particle. The value of s is 1 and 0.5 for the cases of spherical and cylindrical nanoparticles.

The volume (φ) and mass fractions (ψ) of the nanoparticles in the solution are associated as

$$\phi =\frac{\psi {\rho }_f}{{\rho }_s+\psi {\rho }_f-\psi {\rho }_s}$$

(8)

The effective volume fraction φ is defined [49] as

$$\phi =\frac{{\eta }_r-1}{2.5}$$

(9)

The relative viscosity η_r is the ratio of nanofluid viscosity (μ_nf) to the viscosity of the base fluid ($k_f$).

3.2. Eulerian–Eulerian Model

The investigated Eulerian–Eulerian models include mixture, VOF and Eulerian models based on the coupling between the phases and considered interphase interactions.

3.2.1. Mixture Model

In the mixture model, the continuity, momentum and energy equations are solved for the mixture by assuming a weak interphase coupling. Each phase is presumed to exhibit its peculiar velocity field vector and volume fraction. The thermophysical properties of the mixture are defined as the weighted average of all the phases corresponding to their volume fraction. The primary phase (liquid) influences the secondary phases (nanoparticles) via drag force while the secondary phases affect the primary phase by reducing its mean momentum. The governing equations incorporated in this model are

$$\nabla \cdot \left({\rho }_{nf}{\overrightarrow{V}}_{nf}\right)=0$$

(10)

$$\nabla \left({\rho }_{nf}{\overrightarrow{V}}_{nf}{\overrightarrow{V}}_{nf}\right)=-\nabla P+\nabla ·\left({\mu }_{nf}\nabla {\overrightarrow{V}}_{nf})+\nabla ·({\phi }_f{\rho }_f{\overrightarrow{V}}_{dr,f}{\overrightarrow{V}}_{dr,f}+{\phi }_s{\rho }_s{\overrightarrow{V}}_{dr,s}{\overrightarrow{V}}_{dr,s}\right)+{\rho }_{nf}g$$

(11)

$$\nabla ·\left({\phi }_f{\overrightarrow{V}}_f{\rho }_f{C}_{p_f}{T}_f+{\phi }_s{\overrightarrow{V}}_s{\rho }_s{C}_{p_s}{T}_s\right)=\nabla ·\left(k_{nf}\nabla T_{nf}\right)$$

(12)

The summation of volume fraction of all the phases is given by ${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\phi }}_{\mathrm{i}}=1$, where n is the number of phases.

The nanofluid mixture velocity is determined using ${\overrightarrow{\mathrm{V}}}_{\mathrm{m}}=\frac{{\mathsf{\phi }}_{\mathrm{f}}{\overrightarrow{\mathrm{V}}}_{\mathrm{f}}{\mathsf{\rho }}_{\mathrm{f}}+{\mathsf{\phi }}_{\mathrm{s}}{\overrightarrow{\mathrm{V}}}_{\mathrm{s}}{\mathsf{\rho }}_{\mathrm{s}}}{{\mathsf{\rho }}_{\mathrm{n}\mathrm{f}}}$ and the drift velocity of the ith phase is approximated as ${\overrightarrow{\mathrm{V}}}_{\mathrm{d}\mathrm{r},\mathrm{i}}={\overrightarrow{\mathrm{V}}}_{\mathrm{i}}-{\overrightarrow{\mathrm{V}}}_{\mathrm{n}\mathrm{f}}$. Any of the mixture properties (J) such as thermal conductivity, density and viscosity are determined as

$$N={\displaystyle \sum }_{i=1}^n{\phi }_i{J}_i$$

(13)

$$\nabla ·\left({\phi }_s{\rho }_s{\overrightarrow{V}}_{nf}\right)=-\nabla ·\left({\phi }_s{\rho }_s{\overrightarrow{V}}_{dr,s}\right)$$

(14)

The correlation (Equation (15)) of Manninen et al. [50] is used to estimate the slip velocity, that is, the relative velocity of nanoparticles w.r.t. base fluid (${\overrightarrow{\mathrm{V}}}_{\mathrm{s}\mathrm{f}}={\overrightarrow{\mathrm{V}}}_{\mathrm{s}}-{\overrightarrow{\mathrm{V}}}_{\mathrm{f}}$):

$${\overrightarrow{V}}_{fs}=\frac{{\rho }_s{d}_s^2}{18{\mu }_f{F}_{drag}}\frac{{\rho }_f-{\rho }_{nf}}{{\rho }_f}\overrightarrow{a}$$

(15)

where

$$\overrightarrow{a}=\overrightarrow{g}-\left({\overrightarrow{V}}_{nf}·\nabla \right){\overrightarrow{V}}_{nf}$$

(16)

Drag force is estimated by the correlation of Schiller and Naumann [51] as a function of the Reynolds number:

$$F_{drag}=\{\begin{array}{c}1+0.15R{e}_s^{0.687}\\ 0.0183R{e}_s\end{array} \begin{array}{c}R{e}_s\le 1000\\ R{e}_s>1000\end{array}$$

(17)

where

$$R{e}_s=\frac{{\overrightarrow{V}}_{nf}{d}_s}{v_{nf}}$$

(18)

In Equation (18), $d_s$ denotes the nanoparticle diameter that is assumed as the diameter of a spherical particle with the equivalent surface area.

3.2.2. Volume of Fluid Model

In this model, the interface between phases is tracked by solving a continuity equation for the volume fraction of the fluid phase. A momentum equation is solved to compute the velocity components shared by all phases. Similarly, a shared temperature is estimated from the single energy equation. The governing equations of the model are

$$\nabla \cdot \left({\phi }_f{\rho }_f{\overrightarrow{V}}_f\right)=0$$

(19)

$$\nabla \left({\rho }_{nf}{\overrightarrow{V}}_{nf}{\overrightarrow{V}}_{nf}\right)=-\nabla P+\nabla ·\left({\mu }_{nf}\nabla {\overrightarrow{V}}_{nf}\right)+{\rho }_{nf}g$$

(20)

$$\nabla ·\left({\overrightarrow{V}}_{nf}{\rho }_{nf}{C}_{p_{nf}}{T}_{nf}\right)=\nabla ·\left(k_{nf}\nabla T_{nf}\right)$$

(21)

The summation of the volume fractions of all the phases is given by ${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\phi }}_{\mathrm{i}}=1$, where n is the number of phases. Similar to the mixture model, the thermophysical properties of the solution are calculated by taking a weighted average of different phases based on their volume fractions (Equation (13)).

3.2.3. Eulerian Model

In the Eulerian model, the coupling between the phases is presumed to be strong and pressure is shared by all the phases. Separate continuity, momentum and energy equations are simulated for all the primary and secondary phases. The interphase momentum transfer is given by the drag and Saffman lift force, whereas particle–particle interactions are defined in terms of collision forces. The volume of individual phases is calculated by integrating their volume fractions throughout the domain. The summation of the volume fractions of all the phases is given by ${\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathsf{\phi }}_{\mathrm{i}}=1$. The governing equations of the primary (liquid) and secondary (nanoparticle) phases are

$$\nabla \cdot \left({\phi }_f{\rho }_f{\overrightarrow{V}}_f\right)=0$$

(22)

$$\nabla \cdot \left({\phi }_s{\rho }_s{\overrightarrow{V}}_s\right)=0$$

(23)

$$\nabla \left({\phi }_f{\rho }_f{\overrightarrow{V}}_f{\overrightarrow{V}}_f\right)=-{\phi }_f\nabla P+{\phi }_f\nabla ·\left({\mu }_f\nabla {\overrightarrow{V}}_f\right)+{\phi }_f{\rho }_f\overrightarrow{g}+F_{drag,f}+{\overrightarrow{F}}_{lift,f}$$

(24)

$$\nabla \left({\phi }_s{\rho }_s{\overrightarrow{V}}_s{\overrightarrow{V}}_s\right)=-{\phi }_s\nabla P+{\phi }_s\nabla ·\left({\mu }_s\nabla {\overrightarrow{V}}_s\right)+{\phi }_s{\rho }_s\overrightarrow{g}+F_{drag,s}+{\overrightarrow{F}}_{lift,s}++{\overrightarrow{F}}_{col,s}$$

(25)

$$\nabla ·\left({\phi }_f{\overrightarrow{V}}_f{\rho }_f{C}_{p_f}{T}_f\right)=\nabla ·\left({\phi }_f{k}_f\nabla T_f\right)-h_v\left(T_s-T_f\right)$$

(26)

$$\nabla ·\left({\phi }_s{\overrightarrow{V}}_s{\rho }_s{C}_{p_s}{T}_s\right)=\nabla ·\left({\phi }_s{k}_s\nabla T_s\right)-h_v\left(T_s-T_f\right)$$

(27)

The lift ($F_{lift}$), collision ($F_{col,s}$) and the drag forces ($F_{drag}$) are estimated by Equations (28)–(31), respectively.

$$F_{lift,f}=-0.5{\rho }_s{\phi }_f\left({\overrightarrow{V}}_s-{\overrightarrow{V}}_f\right)\times \left(\nabla \times {\overrightarrow{V}}_s\right)$$

(28)

$$F_{lift,s}=-0.5{\rho }_f{\phi }_s\left({\overrightarrow{V}}_f-{\overrightarrow{V}}_s\right)\times \left(\nabla \times {\overrightarrow{V}}_f\right)$$

(29)

$$F_{col,s}=\exp \left[-600\left({\phi }_s-0.376\right)\right]\left(\nabla {\phi }_s\right)$$

(30)

$$F_{drag,f} =F_{drag,s}=\frac{{\phi }_f{\phi }_s{f}_{D}18{\mu }_f}{d_s^2}\left({\overrightarrow{V}}_f-{\overrightarrow{V}}_s\right)$$

(31)

$$f_D=\frac{C_{D}R{e}_s}{24}$$

(32)

$$C_D=\{\begin{array}{c}\frac{24\left(1+0.15R{e}_s^{0.687}\right)}{R{e}_s}\\ 0.44\end{array} \begin{array}{c}R{e}_s\le 1000\\ R{e}_s >1000\end{array}$$

(33)

$F_{lift}$ is determined using the Drew and Lahey [52] equation (Equations (28) and (29)), whereas $F_{col,s}$ is estimated by the Bouillard et al. [53] correlation. In Equation (32), the drag coefficient ($f_D$) is estimated by the Schiller and Naumann [51] correlation. The Reynolds number of the secondary phase in Equation (33) is defined as

$$R{e}_s =\frac{d_s\left|{\overrightarrow{V}}_s-{\overrightarrow{V}}_f\right|d_k{\rho }_f}{{\mu }_f}$$

(34)

The volumetric interphase heat exchange (Equations (26) and (27)) is calculated using the following equation:

$$h_v=6k_s{\phi }_s{\phi }_{f}N{u}_s/d_s^2$$

(35)

where $N{u}_s$ is estimated using the Ranz and Marshall correlation [54]:

$$N{u}_s=2+0.6R{e}_s^{0.5}P{r}_s^{0.33}$$

(36)

In all the Eulerian–Eulerian models, the diameter of the spherical particles with the equivalent surface area is defined.

3.3. Lagrangian–Eulerian Model

In the Lagrangian–Eulerian model, the trajectories of nanoparticles are traced in the Lagrangian approach while the governing equations of the primary phase are computed in the continuum (Eulerian) approximation. The momentum and energy exchange between the primary and particulate phases are simulated by incorporating source terms in the momentum ($S_m$) and energy ($S_e$) equations of the primary phase. The momentum interactions include particle Brownian motion, drag and Saffman lift forces, thermophoresis, virtual mass, gravity and pressure-gradient-induced forces. The interparticle collisions are supposed to be negligible considering the small weight concentration of nanoparticles. The MWCNT and GNP nanoparticles are assumed as tubular- and thin-disc-shaped, respectively. The governing equations of the primary phase are

$$\nabla \cdot \left({\rho }_f{\overrightarrow{V}}_f\right)=0$$

(37)

$$\nabla \left({\rho }_f{\overrightarrow{V}}_f{\overrightarrow{V}}_f\right)=-\nabla P+\nabla ·\left({\mu }_f\nabla {\overrightarrow{V}}_f\right)+S_m$$

(38)

$$\nabla ·\left({\overrightarrow{V}}_f{\rho }_f{C}_{p_f}{T}_f\right)=\nabla ·\left(k_f\nabla T_f\right)+S_e$$

(39)

where

$$S_m={{\displaystyle \sum }}^{}{n}_s\frac{m_s}{\delta \mathcal{V}} \frac{d{V}_s}{dt}$$

(40)

$$\frac{d{V}_s}{dt}=F_{drag}+F_{gravity}+F_{Brownian}+F_{thermophoresis}+F_{lift}+F_{virtual}+F_{pressure}$$

(41)

where $F_{drag}$ and $F_{lift}$ are estimated by the Stokes drag and Saffman’s lift [55] correlations, respectively.

$$F_{drag}=\frac{18{\mu }_f}{d_p^2{\rho }_p{C}_c}\left(V_f-V_s\right)$$

(42)

$$F_{gravity}=\frac{g\left({\rho }_s-{\rho }_f\right)}{{\rho }_s}$$

(43)

$$F_{Brownian}= {\zeta }_i\sqrt{\frac{\pi }{\Delta t}\frac{216{\nu }_f{K}_{B}T}{{\pi }^2{\rho }_f{d}_p^5{\left(\frac{{\rho }_s}{{\rho }_f}\right)}^2{C}_c}}$$

(44)

$$F_{thermophoresis}=\frac{7.02\pi d_s{\mu }_f^2\left(K+2.18K_n\right)}{m_s{\rho }_{f}T\left(1+3.42K_n\right)\left(1+2K+4.36K_n\right)}\times \frac{\partial T}{\partial x}$$

(45)

$$F_{lift}=\frac{5.188v_f^{0.5}{\rho }_f{d}_{ij}}{{\rho }_s{d}_s{\left(d_{lk}{d}_{kl}\right)}^{0.25}}\left(V_f-V_s\right)$$

(46)

$$F_{virtual}=0.5\frac{{\rho }_f}{{\rho }_s}\frac{d}{dt}\left(V_f-V_s\right)$$

(47)

$$F_{pressure}=\left(\frac{{\rho }_f}{{\rho }_s}\right)V_s\frac{d{V}_f}{dx}$$

(48)

In Equation (45), $K=\frac{k}{k_s}.$ The $d_{ij}$ and $K_n=\frac{2\mathsf{\lambda }}{d_s}$ are the deformation tensor and Knudsen number, respectively. ${\zeta }_i$ in Equation (44) is the zero-mean, unit-variance-independent Gaussian random numbers. The $C_c$ is the Cunningham correction factor of the spherical nanoparticles. The modified Cunningham correction factors for cylindrical (MWCNT) and thin-disc (GNP)-shaped nanoparticles are given in Equations (49) and (50) [56], respectively.

$$C_c=\frac{3d_L\left[ln\left(2\beta \right)-0.72\right]}{2}$$

(49)

$$C_c=\frac{3\pi }{8}$$

(50)

In the above equations, β denotes the aspect ratio of the particle length to diameter. The energy source term $S_e$ in Equation (39) is defined as

$$S_e={{\displaystyle \sum }}^{}{n}_s\frac{m_s}{\delta \mathcal{V}} C_p\frac{d{T}_s}{dt}$$

(51)

$$m_s{C}_s\frac{d{T}_s}{dt}=h_i{A}_s\left(T_f-T_s\right)$$

(52)

$$h_i=\frac{N{u}_s{k}_f}{d_s}$$

(53)

where $N{u}_s$ is estimated by incorporating the Ranz and Marshall model (Equation (36)) [54].

For all the tested models, the thermophysical properties of the nanoparticles are listed in

Table 1. Thermophysical properties of multiwalled carbon nanotubes (MWCNT) and graphene nanoplatelets (GNP) nanoparticles.

Material	$\mathbf{\rho }\left[\mathbf{k}\mathbf{g} {\mathbf{m}}^{-3}\right]$	$\mathit{C}\mathit{p}$$\left[\mathbf{J} \mathbf{k}{\mathbf{g}}^{-1}{\mathbf{K}}^{-1}\right]$	$\mathit{k}$$\left[\mathbf{W} {\mathbf{m}}^{-1}{\mathbf{K}}^{-1}\right]$	Ref.
MWCNT	2100	630	1500	[57]
GNP	2200	790	3000	[45]

. Base fluid (primary phase for multiphase models) temperature-dependent properties are assessed using the following expressions:

$${\rho }_f= 2446-20.674T+0.11576 T^2-3.12895\ast {10}^{-4}{T}^3-2.0546\ast {10}^{-10}{T}^5$$

(54)

$${\left(C_p\right)}_f=\exp \left(\frac{8.29041-0.012557T}{1-\left(1.52373\ast {10}^{-3}\right)T}\right)$$

(55)

$${\mu }_f=2.414\ast {10}^{-5}\ast {10}^{\left(\frac{247.8}{T-140}\right)}$$

(56)

$$k_f= 0.5706+1.756\ast {10}^{-3} T -6.46\ast {10}^{-6} T^2$$

(57)

In the single-phase model, the nanoparticle shape effect is modelled through the shape factor (n) in the effective thermal conductivity of the nanofluids (Equation (7)). The nanoparticles are assumed as spherical-shaped with the surface area equivalent to the considered nanoparticles in all the Eulerian–Eulerian models (mixture, VOF and Eulerian), whereas in the Lagrangian–Eulerian model, the nanoparticle shape is modelled through the modified Cunningham correction factor (Equations (49) and (50)) in the particle drag (Equation (42)) and Brownian motion (Equation (44)).

4. Numerical Details and Grid Independency

The numerical analysis is conducted by employing the finite-volume-method-based commercial software ANSYS Fluent [58]. The pressure-based coupled algorithm is used for pressure–velocity coupling and the second-order upwind scheme is utilised for the discretisation of the convective terms [58]. To track the secondary-phase volume fraction, the implicit formulation is used in VOF and Eulerian models. The first-order upwind scheme is applied for the discretisation of volume fraction in both mixture and Eulerian models while, in the case of the VOF model, a compressive scheme is incorporated. The nanoparticles are introduced into the computational flow domain as surface injection type in the Lagrangian–Eulerian model [10,42,59]. The trapezoidal high-order scheme and implicit low-order scheme are implemented for nanoparticle tracking.

At the inlet, outlet and wall of the minichannel, boundary conditions of velocity inlet, pressure outlet and constant heat flux are respectively defined in accordance with their respective experimental test conditions [45]. Additional escape boundary conditions at the channel inlet and outlet and reflect boundary conditions at the channel wall are defined for the nanoparticle tracing in the Lagrangian–Eulerian model. Considering the symmetry of the computational domain, half of the channel is modelled by implying symmetry boundary condition (see Figure 2). Residuals of volume fraction, continuity, momentum and energy equations are ensured to be less than ${10}^{-7}$ for the convergence criteria.

A structured hexahedral grid is generated over the whole computational domain using the ICEM CFD package as depicted in Figure 2. The mesh is refined near the channel wall to capture the thermal–viscous boundary layer. Four grid sizes (G1–G4) of different resolutions are tested for the forced convection of the MWCNT/GNP–H₂O hybrid nanofluid (ψ = 0.25%, Re = 470) by utilising all the tested single and multiphase models to ensure the grid-insensitive solution. The specifics of the tested mesh sizes are listed in

Table 2. Specifics of tested grid sizes.

Grid	$\mathbf{\Delta }\mathit{t}$	${\mathit{N}}_{\mathit{t}}$
G1	0.058	826,541
G2	0.044	1,102,676
G3	0.031	1,442,556
G4	0.023	2,302,986

. The results approximated by different grid sizes are compared in terms of average heat transfer coefficient (h_avg) and pressure drop (ΔP) as illustrated in Figure 3. The results verify that h_avg and ΔP values show nearly constant values when computed utilising grid G3 and G4 with the maximum percentage difference less than 0.5%. Thereby, grid G3 is utilised for all the subsequent computations (Y+ < 1) to assure solution accuracy at the expense of minimal computational cost.

5. Results and Discussion

In this section, the accuracies of the single-phase and multiphase Eulerian–Eulerian (VOF, mixture, Eulerian) and Lagrangian–Eulerian approaches for simulating the nonspherical MWCNT and MWCNT-/GNP-nanoparticles-based nanofluids are evaluated. In all the tested models, the influence of nanolayering around the nanoparticles is incorporated in terms of effective volume fraction that describes the combined effects of nanoparticle volume fraction and surfactants on the rheological properties of the solution. The range of Reynolds number and nanoparticle weight concentrations are 200 ≤ Re ≤ 470 and 0.075% ≤ ψ ≤ 0.25%, respectively. The accuracy of the computed hydrothermal (average and local heat transfer coefficient, pressure drop) results is assessed through comparison with the relevant experimental results [45].

5.1. Validation of the Homogenous Model

The numerically calculated local Nusselt number results (along the channel wall) of the water (primary phase) are validated with the theoretical correlation of the Shah and London [60] (Equation (58)) for laminar flow under the boundary condition of constant heat flux at Re = 470:

$$\{\begin{array}{ccc}N{u}_x=3.302x_{+}^{-1/3}& for& x_{+}\le 0.00005\\ N{u}_x=1.302x_{+}^{-1/3}-0.5& for& 0.00005\le x_{+}\le 0.0015\\ N{u}_x=4.362+8.68{\left({10}^3{x}_{+}\right)}^{-0.506}{e}^{-41x_{+}}& for& x_{+}\ge 0.0015\end{array}$$

(58)

where

$$x_{+}=\frac{X}{D_{h}Re Pr}$$

(59)

In addition, the results of friction factor in fully developed flow region of the channel are compared with those estimated with the Darcy friction factor ($f=64/Re$) correlation [61]. Figure 4 and Figure 5 validate that the numerical results are in acceptable agreement with the theoretical values with a maximum percentage deviation of 10.6% and 1.9% for the local Nusselt number and friction factor, respectively. The average percentage difference between the computed and theoretical values of the local Nusselt number is 3.1%.

5.2. Heat Transfer Characteristics

Figure 6 presents the average convective heat transfer coefficient (h_avg) results of MWCNT–H₂O and MWCNT/GNP–H₂O nanofluids predicted using single and multiphase models for the nanoparticle concentrations of 0.25%. The figures illustrate that at ψ = 0.25%, the multiphase Eulerian–Eulerian (mixture, VOF and Eulerian) and single-phase models overestimate the h_avg results of both nanofluid samples. Such overprediction of the thermal results computed by the Eulerian–Eulerian model is also highlighted in literature and was correlated with the overestimation of the volume-weighted average thermal conductivity of nanofluids [42,62]. All of the Eulerian–Eulerian models estimate quite similar values, whereas those of the single-phase model are slightly larger. The deviation of the Eulerian–Eulerian and single-phase models is more pronounced at high Reynolds numbers when the interphase interactions intensify, i.e., the percentage difference of h_avg of MWCNT–H₂O (ψ = 0.25%) nanofluid with the experimental data is 6.8% at Re = 200, which increases to 12.1% at Re = 470. Similarly, the percentage difference between the experimental and numerical values of Eulerian–Eulerian and single-phase models is more significant for the MWCNT/GNP–H₂O hybrid nanofluid than for the MWCNT/H₂O. This trend can be explained as the outcome of the stronger interphase interactions provided by the larger nanoparticle concentrations, given that the hybrid nanofluid comprises an additional 0.035 wt% of the GNP nanoparticles. Of all the single and multiphase approaches, the Lagrangian–Eulerian treatment of nanofluids results in the most accurate values with the maximum deviations of 5.34% and 4.9% for the MWCNT and MWCNT/GNP nanofluids, respectively, at Re = 470.

Figure 7 and Figure 8 reflect the h_avg results of both the nanoparticle samples for the nanoparticle concentrations of 0.125% and 0.075%, respectively. It is apparent from the figures that all the Eulerian–Eulerian models compute identical yet slightly overestimated results. The accuracy of the tested models for predicting h_avg values is in the order of Lagrangian–Eulerian, Eulerian–Eulerian and single-phase models. However, comparing the h_avg values at ψ = 0.25%, the deviation between the experimental and numerical results decreases for small nanoparticle concentrations as an outcome of insignificant interphase interactions (Figure 9). At Re = 470, the percentage differences between the experimental and numerical results computed by the single, mixture, VOF, Eulerian and Lagrangian–Eulerian models are 6.7%, 6.3%, 6.4%, 5.5% and 3.7%, respectively, for the MWCNT-based nanofluids at ψ = 0.125%. This deviation reduces to 1.9%, 2%, 3.4%, 1.9% and 3.1%, respectively, for particle concentration ψ = 0.075%. Similarly, for MWCNT/GNP-based (ψ = 0.125%, Re = 470) nanofluid, the percentage errors between the experimental and computed results are 10.73%, 9.5%, 9%, 9.4% and 3.5% for the single, mixture, VOF, Eulerian and Lagrangian–Eulerian models, respectively. This deviation becomes insignificant with 4.9%, 5.23%, 5.43%, 5% and 2.56%, respectively, at ψ = 0.075%. Irrespective of the nanofluid type, h_avg increases with the augmentation of nanoparticle concentration and Re. Such enhancement in heat transfer is consistent with the literature [2,63,64,65]. However, the rate of h_avg increment is highest for the single-phase and Eulerian–Eulerian models, followed by the Lagrangian–Eulerian model. The h_avg of the MWCNT/GNP–H₂O hybrid nanofluid is larger than that of the MWCNT–H₂O nanofluid. At optimal Re and ψ, the enhancements of h_avg of MWCNT/GNP–H₂O nanofluid compared with MWCNT–H₂O nanofluid are 14%, 12.2%, 12.9%, 13.1% and 9.1% when estimated using single, mixture, VOF, Eulerian and Lagrangian–Eulerian models, respectively.

To further elaborate the accuracy trends of the tested models in computing the thermal characteristics of the nanofluid, Figure 10 demonstrates the local heat transfer coefficient (h_local) distribution along the nondimensional streamwise wall length of the channel for MWCNT–H₂O and MWCNT/GNP-H₂O nanofluids at ψ = 0.25%. The reduction of h_local from the highest value at the entrance region of the channel to subsequently stabilised low values refers to thermally and hydrodynamically developing and fully developed flow regions. In agreement with previously stated results, a comparison of the experimental and numerical data signifies the superiority of the Lagrangian–Eulerian approach in predicting the h_local distribution of nanofluids over single and other multiphase models. The latter models predict significantly low h_local values along z/D_h > 75 while overestimating the results along the remaining tube length 75< z/D_h < 250. For the MWCNT-H₂O (ψ = 0.25%) nanofluid, the deviations of h_local results determined using the single-phase, mixture, VOF, Eulerian and Lagrangian–Eulerian models are 13.3%, 11.9%, 12.6%, 11.68% and 6.85%, respectively, while for the MWCNT/GNP-H₂O (ψ = 0.25%) nanofluid, these deviations become 17.5%, 14%, 14.4%, 13.92% and 5.47%, respectively.

The h_local for the small nanoparticle concentrations of 0.125 wt% and 0.075 wt% (Figure 11) also illustrate the accuracy trend of the tested models in the order of single, Eulerian–Eulerian and Lagrangian–Eulerian. However, for the highly diluted nanofluids, the tested models show quite insignificant deviations, that is, at ψ = 0.075%, the percentage differences between the experimental and computed h_local data of single-phase, mixture, VOF, Eulerian and Lagrangian–Eulerian models are 3%, 3.2%, 4.8%, 3% and 4.7%, respectively. Thus, it can be concluded that all the single- and multiphase numerical models estimate quite accurate thermal fields for highly diluted nanofluids when the interphase interactions also become insignificant. However, single-phase and Eulerian–Eulerian models fail to predict the thermal characteristics of dense nanofluids due to the overestimation of the effective thermal conductivity and interphase interactions.

Figure 12 depicts the temperature variation along the channel central longitudinal and cross-sectional planes (z/D_h = 0, 49, 98, 147, 196, 245) for the MWCNT (left) and MWCNT/GNP (right) nanofluids at ψ = 0.25%. As stated earlier, the thermal contours also verify that the MWCNT/GNP–H₂O hybrid nanofluid improves the heat transfer rate compared with the MWCNT–H₂O nanofluid. However, this heat transfer enhancement rate predicted by different models is in the order of single-phase, Eulerian–Eulerian and Lagrangian–Eulerian models.

5.3. Pressure Drop Characteristics

Figure 13 presents the compression of pressure drop results of MWCNT–H₂O and MWCNT/GNP–H₂O nanofluid at ψ = 0.25%. Clearly, the accuracies of the pressure drop results estimated by all the models are in acceptable range, and are in the order of single-phase, Eulerian–Eulerian and Lagrangian–Eulerian models. Of all the tested models, the Lagrangian–Eulerian approximation demonstrates the most deviated pressure drop values due to the inclusion of interphase interactions, yet the computed results are in the acceptable limit. Similarly, all the tested models determine quite similar and accurate results for the smaller nanoparticle concentration and at low Re (Figure 14 and Figure 15). As the nanoparticle concentration or the flow rate increases, the deviation of the computed results also increases. At Re = 470, the percentage errors of the single-phase, mixture, VOF, Eulerian and Lagrangian–Eulerian models are 7%, 5.7%, 6.7%, 5% and 11.1% for the MWCNT/GNP–H₂O (ψ = 0.25%) nanofluid, while decreasing to 3.4%, 2.9%, 4.6%, 2.1% and 6.5% for the nanoparticle concentration of 0.075%, respectively.

6. Conclusions

The homogenous, Eulerian–Eulerian (VOF, mixture, Eulerian) and Lagrangian–Eulerian approaches for simulating the hydrothermal characteristics of MWCNT–H₂O and MWCNT/GNP–H₂O hybrid nanofluids are assessed by incorporating the influence of nanolayering around the nanoparticles. The analysis is performed for the laminar forced convection of nanofluids in a minichannel subjected to constant heat flux by considering the nanoparticle weight concentration and Reynolds number of 0.075–0.25% and 200–470, respectively. The results are evaluated in terms of average and local convective heat transfer coefficients, thermal contours and pressure drop across the channel. The study findings can be concluded as follows.

• Single and all the multiphase numerical models estimate reasonably accurate convective heat transfer coefficient results for diluted nanofluids and at low Reynolds number when the interphase interactions are negligible.
• The single-phase and Eulerian–Eulerian models overestimate the thermal fields of the nanofluid with a more perceptible difference at high nanoparticle concentration and Reynolds number.
• The multiphase Eulerian–Eulerian approaches reveal marginal differences between average and local convective heat transfer coefficients results. The maximum deviation of the single-phase, VOF, mixture and Eulerian models are 10.73%, 9.5%, 9% and 9.4%, respectively, for the average heat transfer coefficient results of the MWCNT/GNP–H₂O (0.25 wt%) nanofluid at Reynolds number of 470.
• Of all the tested models, the Lagrangian–Eulerian approximation of nanofluid provides the most accurate convective heat transfer coefficient with a maximum deviation of 5.34% for 0.25 wt% of MWCNT–water nanofluid at Reynolds number of 470. However, the model requires high memory and computational time to compute the trajectories of individual nanoparticles.
• Despite the thermal discrepancies, all numerical models determine quite accurate pressure drop results at all the studied nanoparticle concentrations and Reynolds numbers.