Nanofluid Heat Transfer in Wavy-Wall Channels with Different Geometries: A Finite-Volume Lattice Boltzmann Study

In this work, we perform an extensive numerical investigation of the heat transfer behavior of nanofluid laminar flows, in wavy-wall channels. The adopted computational approach is based on a finite-volume formulation of the lattice Boltzmann method constructed on a fully-unstructured mesh. We show the validity and effectiveness of this numerical approach to deal with realistic problems involving nanofluid flows, and we employ it to analyze the effects of the wavy-wall channel geometry on the rate of heat transfer, thus providing useful information to the design of efficient heat transfer devices. Results show that an increasing of the wavy surface amplitude has a positive effect on the heat transfer rate, while a phase shift between the wavy walls leads to a decreasing of the mean Nusselt number along the channel. The addition of solid nanoparticles within a base liquid significantly contributes to increase the rate of heat transfer, especially when a relatively high value of nanoparticles volume fraction is employed. The present analysis then suggests that the use of nanofluids within an axis-symmetric configuration of the wavy-wall channel, with high wavy surface amplitude, may represent an optimal solution to enhance the thermal performances of heat transfer devices.


Introduction
The low thermal conductivity of common heat transfer fluids such as water, oil and ethylene glycol represents a major limitation to the development of high-efficient heat exchangers, that are required in many industrial applications. Enhancing the thermal performance of heat transfer devices is instrumental to reduce both the energy consumptions and the associated operating costs. Therefore, a significant effort has been making to design and develop advanced heat transfer fluids that exhibit improved thermal behavior. In this sense, an interesting solution is represented by nanofluids [1], that are colloidal suspensions of nanometer sized particles. Since their first introduction by Choi [2] few decades ago, nanofluids have continued to attract the research interest of the scientific community due to their remarkable thermophysical properties. Commonly used nanofluids present solid particles of three kinds, mainly: metal oxides, pure metal and carbon-based materials, such as graphene. The dispersion of these solid nanoparticles in a base liquid leads indeed to superior characteristics compared to conventional heat transfer fluids, as demonstrated by several experimental [3,4] and numerical studies [5,6]. A review on the application of nanofluids in heat exchangers can be found in [7]. The application of nanoparticles to a liquid matrix has a great potential for the energy management and control in a variety of application cases. Among these, the use of nanofluids in thermal energy storage systems operating with phase change materials [8][9][10] or involving the use of open-cell metal foams structures [11][12][13][14][15][16] represents one of the most interesting topic in the field of engineering research [17,18].
The numerical methods adopted to investigate the heat transfer mechanisms in nanofluids can be divided into two main classes: single phase (homogeneous) and multiphase (suspended particles) [19,20]. In the first one, the liquid carrier and the solid nanoparticles are assumed to be in thermal equilibrium to each other and, therefore, the nanofluid is treated as an homogeneous fluid with thermophysical properties determined on the basis of empirical or theoretical models [21,22]. Differently, in the multiphase approach, both the base fluid and the solid nanoparticles are modeled individually, thus allowing a detailed analysis of their interaction [23,24]. Most of the studies available in literature make use of the single phase method, which is often accurate enough for engineering problem with a higher simplicity and a lower computational cost with respect to the multiphase approach [1]. In this paper we employ the single-phase approach with the Brinkman model [25] and the Maxwell-Garnetts model [26] to estimate the thermophysical properties of the considered nanofluids, as will be detailed later.
Relatively recent advances in computational fluid dynamics have lead to the development of high efficient numerical methods based on the kinetic theory of gases. Among these, the lattice Boltzmann methods (LBMs) [27][28][29][30][31][32][33] have been successfully used in a broad range of engineering applications, such as, for instance, free-surface water entry problems [34,35], turbulent flows [36][37][38], reactive flows [39] and multiphase flows [40][41][42]. Different LBM realizations have been also applied to study the thermal behavior of nanofluids, either by using the single phase or the multiphase approach [43][44][45][46][47][48]. The majority of these works adopt the original formulation of LBMs, which is based on a uniform space discretization. This feature represents, however, a strict limitation for those problems involving complex domain geometries. This is the case, for instance, of the heat exchangers design. To cope with this issue, several extensions of the LBM have been proposed in literature, which are based on body-fitted mesh approaches. These includes finite-differences [49], finite-element [50], finite-volume [51][52][53][54], semi-lagrangian [55] and hybrid [56][57][58] methods. In particular, in [59] and in [60] finite-volume formulations of the lattice Boltzmann equation on non-uniform structured meshes are proposed and applied to model the thermo-fluid behavior of nanofluids.
Despite the numerous studies on off-lattice Boltzmann methods, there is still a lack of research in the field of nanofluids to deal with realistic problems. Therefore, in this paper, we extend a fully-unstructured method based on a finite-volume lattice Boltzmann formulation, originally developed by the the authors in [51], to account for the thermo-fluid dynamics phenomena occurring within nanofluid flows. We then apply the validated method to study a nanofluid flow through a wavy-wall channel. We perform an extensive parametric analysis by considering two types of water-based nanofluids, with aluminium oxide (Al 2 O 3 ) and copper oxide (CuO) nanoparticles, respectively. For each type of nanofluid we consider three concentration levels of nanoparticles within the base fluid, namely 2%, 5% and 10% and we analyze the heat transfer behavior of the nanofluids in different geometrical configurations of the channel, which are obtained by varying both the phase shift between the walls and the amplitude of the wavy walls themselves. The adopted numerical approach is meant to represent a viable option, within the LBM realm, to deal with the geometrical complexity of practical engineering problems involving heat transfer in nanofluids. The fully-unstructured mesh allows indeed for a more efficient geometrical representation of complex boundaries, with respect to the original Cartesian formulation of LBM, while keeping the computational burden at a relatively low level. In this sense, the inherent geometrical complexity of the wavy-wall channel does represent an interesting benchmark problem to evaluate the capability of the proposed approach of capturing the relevant physics behind heat transfer in nanofluids. The results show that our method has significant potential for the accurate calculation of nanofluid flows in complex geometries domain.
The remaining part of the paper is organized as follows. In Sect. 2 we present the mathematical model and the numerical formulation of the thermal lattice Boltzmann finite-volume approach used in this work. In Sect. 3 we provide the thermophysical properties of the considered nanofluids. In Sect. 4 we describe the problem under investigation and we provide details on the adopted numerical setup. In Sect. 5 the results of this work are presented and discussed, while in Sect. 6 we summarize the key points and the main findings of the study.

Mathematical Model
We model the nanofluid through a single phase approach, that is the base fluid and the nanoparticles are assumed to be in thermal equilibrium between each other and to move with the same local velocity. Therefore, we consider the nanofluid flow to be incompressible, Newtonian and to have homogeneous and constant thermophysical properties. Thus, the problem is described by the following governing equations, that are the Navier-Stokes momentum equations and the advection-diffusion heat transfer equation, respectively: where ρ is the fluid density, u is the velocity vector, p is the pressure, c p is the specific heat, k is the thermal conductivity and μ is the dynamic viscosity. Instead of numerically solving the macroscopic equations (1) and (2), the lattice Boltzmann method deals with a specific discretization of the Boltzmann kinetic equation. Thus, in its differential form, the lattice Boltzmann equation that recovers the physics described by (1) reads as follows: In Eq. 3, x represents the position of a particle in the discretized space at time t, while f i (x, t) and c i are the density distribution functions and the set of discrete speeds, respectively, along the i-th allowed lattice directions. The collision operator is characterized here by a single parameter, that is the relaxation time τ (Bhatnagar-Gross-Krook approximation). The equilibrium density distribution function f eq i (x, t) is given by: with being w i a set of weight parameters normalized to unity, while c s indicates the lattice speed of sound. Heat transfer phenomena in the fluid flow are here modelled through the thermal lattice Boltzmann method based on the double distribution function model [61], where the temperature is considered as a passive scalar. Therefore, denoting by g i (x, t) the energy distribution function, we introduce the following equation: where τ g is the thermal relaxation time. The energy equilibrium distribution function takes the same form of Eq. 4, therefore: where T is the temperature. The macroscopic variables, that are ρ, u and T , in lattice units, are calculated as follows: where Q denotes the number of discrete velocities. In this work, we use a two dimensional lattice with nine allowed lattice directions (D2Q9 scheme).

Numerical Formulation
Equations (3) and (5) are numerically solved through a finite-volume scheme of the cell-vertex type based on a triangular tessellation, originally proposed in [51]. A schematic representation of the space discretization is reported in Fig. 1. The control volume for a generic node is obtained by joining the centroids of the elements around the pivotal node P with the midpoints of the edges having P as a common node. The space-time discretization of Eq. (3) leads to the following expression for each node P of the mesh: where the summation runs over the K triangular elements which share the node P as a common vertex. The terms ik and ik represent streaming and collisional fluxes of the i-th population related to the k-th volume, respectively: where V P is the total area (in 2D) of the control volume related to the node P. The resulting finite-volume equation takes the following general form, which applies to each node P of the unstructured mesh [51]: (10) where S ik and C ik are the streaming matrix and the matrix associated with the collisional operator, respectively. Further details of the present numerical approach can be found in [56].
The same finite-volume scheme so far described for the evolution equation of the density distribution function is then applied to the temperature transport equation 5, which results as follows: In general, the theoretical kinematic viscosity of the lattice Boltzmann equation 3 consists of two separate contributions, namely ν = c 2 s (τ + τ n ), where c 2 s τ is the physical contribution due to collisional relaxation, while c 2 s τ n is the term related to the numerical diffusion. The latter depends on the details of the numerical scheme. For standard LBM, the numerical diffusion amounts to a constant value, given by τ n = −0.5. This is not the case for the unstructured lattice Boltzmann formulation here employed, which instead show a very low degree of numerical diffusivity, so that it results ν = c 2 s τ [51]. The same applies for the relation between the thermal diffusivity and the thermal relaxation time, that is α = c 2 s τ g .

Thermophysical Properties of Nanofluids
In this study we consider two types of water-based nanofluids, that are: Al 2 O 3 -water and CuO-water nanofluid. We report in Table 1 the thermophysical properties of water and the selected solid nanoparticles used in the present analysis, which we consider to be constant within the range of temperature characterizing the problem under investigation.
In addition, we set the dynamic viscosity of water equal to 0.001 Pa·s such that the water Prandtl number is equal to 6.9. To obtain the density and the specific heat of nanofluids we apply the following formulas: where subscripts b f , np and n f refer to base fluid (water), nanoparticles and nanofluid, respectively, while φ is the volume fraction of nanoparticles. The effective dynamic viscosity of the nanofluid is calculated through the Brinkman model [25], which reads as follows: As far as the effective thermal conductivity of nanofluid is concerned, here we adopt the Maxwell-Garnetts formulation [26]: We note that the effective thermal conductivity of nanofluid in the Maxwell-Garnetts model depends only on the thermal conductivities of solid nanoparticles and base liquid, and on their relative volume fraction. Therefore, we do not need to consider the solid nanoparticle size.

Problem Statement and Numerical Setup
We study the forced convective heat transfer of nanofluids in a two-dimensional wavy-wall channel [62]. The computational domain, for a generic configuration of the channel, is shown in Fig. 2.
With reference to the characteristic dimensions reported in Fig. 2, the shapes of top (S t ) and bottom (S b ) wavy-wall profiles are described by the following mathematical functions: where a is the amplitude of the wavy surfaces, H is the half separation distance between the walls, ϕ and ϕ are two phase shift parameters, X in T and X in B are the lengths of the flat walls at the entrance region of the domain. Specifically, in accordance with [62], we fix X in B = 3H and X w T = X w B = X w = 12H (both walls are made of six complete sinusoidal waves) in the whole set of numerical simulations performed in this study. Moreover, the overall extension of the computational domain along the horizontal direction is L = 18H , while the amplitude a and the phase shift parameters ϕ and ϕ are varied parametrically in this study. We define the following dimensionless quantities: where a is the non-dimensional amplitude of the wavy surface and x is the non-dimensional coordinate along the axis of the channel. The present numerical analysis is performed by considering three different values of a, 0.1, 0.2 and 0.3, two values of ϕ , 0 and π, and three values of ϕ , 0, π/4 and π/2. In the whole set of numerical simulations concerning the wavy-wall channel flow problem we set Re = 100, being the Reynolds number defined as in [62]: where u in is the mean velocity at inlet section and H is the half separation distance between the walls (Fig. 2). Also, in all of the analyzed scenarios, we select the half channel height H to be equal to 120 lattice units. The computational domain is discretized by means of a fully-unstructured mesh of triangular elements. The mesh, for each channel configuration, is constructed to minimize the computational cost, while retaining high resolution in proximity of the wavy-walls, where the highest temperature gradients occur. In particular, the mesh resolution is chosen by preliminary performing a grid sensitivity analysis for the case of water flow through a wavy-wall channel in its base configuration (ϕ = ϕ = 0) with non-dimensional amplitude of the wavy surface equal to 0.2. A detail of the final mesh used for the base configuration of the channel is drawn in Fig. 3, as an example. Similar meshes are used for all the other channel configurations.
The adopted meshes are highly stretched, with the ratio between the largest and the smallest cell size being approximately equal to 5. This leads to an overall amount of roughly 30000 nodes, for each mesh used in the present analysis.
The computational domain is divided into three parts along the horizontal direction: an entrance region, which is considered from the inlet section to the start of the bottom wavysurface (x < X in B ); a middle region, which encloses the whole bottom wavy-surface (X in B ≤ x ≤ X in B + X w ); and an outlet region, between the end of the bottom wavy-surface and the outlet section. We assume that the temperature of both walls in the middle region is held at a constant value T w , while we impose the adiabatic condition to the walls in entrance and outlet regions. Therefore, in the base configuration of the channel (ϕ = ϕ = 0) both wavy-surfaces are characterized by a fixed value of the temperature, while all of the flat walls are considered adiabatic, as in [62]. At the inlet section, the fluid temperature T 0 < T w is assumed to be constant along the channel height. At the outlet, we assume a developed flow, such that the temperature gradient along the horizontal direction is zero. Moreover, a developed velocity profile is set at inlet, while the no-slip condition is applied at both walls. These boundary conditions are summarized in Table 2. In particular, we set T w = 350 K and T 0 = 300 K. The described boundary conditions hold at each computational domain configuration of the wavy-wall channel investigated in this work.
The boundary conditions are applied by means of a simple, yet effective, approach, developed by the authors in [51]. According to this approach, at each boundary, the computational domain is augmented with one buffer of regular, straight triangles, in order to ensure that the last-but-one line of boundary nodes faces a corresponding neighbor along the x or y direction (depending on the type of boundary). Therefore, at the fictitious constructed nodes, we impose (i) the same field (temperature, pressure or velocity) of the adjacent line of nodes, such that a zero-gradient boundary condition automatically results, for the case of a Neumann boundary condition, or (ii) a fixed value for the specific macroscopic variable (temperature, pressure or velocity), for the case of a Dirichlet boundary condition. Thus, we apply a local equilibrium to reconstruct the unknown distribution functions at fictitious nodes, which are required for computing fluxes at the actual boundary nodes. This method was found to yield acceptable results in this work, as well as in previous works where the same finite-volume scheme was employed.
Conversion factors from physical to lattice units The lattice Boltzmann method is based on non-dimensional quantities. Therefore, conversion factors must be defined for space, time, mass and temperature in order to lead the physical problem into the lattice Boltzmann framework. To this aim, we consider pure-water as the reference fluid and then we set the numerical properties, in lattice units (l.u.), of any nanofluid by adopting the defined conversion factors. To set the numerical properties of water, we initially perform a sensitivity analysis to select a reference Mach number (Ma) such that compressibility effects can be neglected in the numerical simulations. Therefore, with Re = 100, H = 120 l.u. and a fixed density of 1 l.u., we carry out a set of simulations for the case of water flow in a wavy-channel in its base configuration by varying Ma between 0.24 and 0.06 and by computing the relative error on the mean Nusselt number (defined in next Section). As a result, for Ma = 0.12 we obtain a percentage error related to the case of Ma = 0.06 roughly equal to 1%. This indicates that compressibility effects may effectively considered not relevant already at Ma = 0.12, which is therefore chosen as the reference value for water flow. Since the speed of sound for the present lattice Boltzmann model is given by c s = 1/ √ 3, the resulting mean velocity at inlet is about 0.07, while the values for kinematic viscosity and relaxation time τ are 0.08 and 0.24 in lattice units, respectively.
The transport of temperature is instead governed by the thermal relaxation time τ g , which depends on the thermal diffusivity, defined as: Therefore, by indicating with subscript LU quantities expressed in lattice units, it results: where k, ρ and c p correspond to quantities in physical units. Equation 21 is obtained by defining the following conversion factors, for space, time, mass and temperature, respectively: and by expressing the thermal conductivity and the specific heat, in lattice units, as follows: Therefore, for water, we obtain α LU = 0.01158 lattice units, which leads to τ g = 0.0347. To obtain the corresponding properties of nanfluids in lattice units we then make use of the conversion factors defined in Eq. (22).

Results and Discussion
The heat transfer behavior of nanofluids is analyzed in this study by evaluating several fundamental parameters, such as local Nusselt number, mean Nusselt number, non-dimensional pressure drop and thermal enhancement factor. Specifically, we define the local Nusselt number at each solid wall of the channel as referred to the thermal conductivity of the base fluid [45,63]. Therefore: where h c is the convective heat transfer coefficient, the temperature gradient at wall is computed by: and T avg (x) indicates the bulk temperature (i.e. average temperature along a vertical section of the channel), which it is obtained by the following expression: We note that the ratio k n f /k b f in Eq. (24) reduces to 1 for water. The integrals in Eq. 26 are computed by selecting 40 equally spaced control points along each considered section of the channel, and then by performing a barycentric interpolation to obtain the values of u and T at each control point.
To compute a reliable value of the local Nusselt number, at each cross section of the channel, we define the following local average Nusselt number: where subscript B and T indicate the bottom and top walls, respectively. We also define a mean Nusselt number along the channel, which characterizes the global heat transfer behavior of the nanofluid, as: where s indicates the length of the channel centerline. Moreover, in order to evaluate the effect of nanoparticles on the heat transfer, we define the enhancement factor as [46]: where Nu m (φ) is the mean Nusselt number related to a nanofluid with a given volume fraction φ and Nu m (φ = 0) indicates the corresponding mean Nusselt number for the base fluid (i.e. water). Finally, we compute the pressure drop along the channel to understand the effects of nanoparticles volume fraction and channel shape on the pumping requirements. Therefore, we define the following non-dimensional pressure drop: where p is the pressure difference between inlet and outlet of the channel, ρ b f is the reference density of the base fluid and u in is the mean inlet velocity.

Preliminary Validation
The numerical model is initially validated by simulating the fluid flow and heat transfer convection of water between two parallel plates. In this case, the computational domain consists of a rectangular geometry with aspect ratio equal to 15. The channel is subjected to constant wall temperature T w , while the temperature of the inflow fluid is fixed at T 0 . Boundary conditions are the same as those presented in Table 2 for the wavy-wall problem. For the present case only, we set the Reynolds number equal to 20, the relaxation time is set to 0.48, while the thermal relaxation time results equal to 0.0695.
As expected, the local Nusselt number, in the fully developed region of the flow, approaches the asymptotic value of 7.54, which is the theoretical value reported in literature [64]. The result is shown in Fig. 4, where the local Nusselt number is plotted against the Graetz number. This is defined as: with being D h the hydraulic diameter, given by twice the spacing of the plates.

Forced Convection of Water in Wavy-Wall Channels
We further validate the method by analyzing the case of forced convection in a wavy-wall channel for pure water, and by comparing the obtained local Nusselt number distribution with  results provided in [62]. The numerical setup corresponds to the one described in previous Section and also adopted for the analysis of nanofluid flows presented later. Validation is performed by analyzing the base configuration of the channel, that is ϕ = ϕ = 0, with the non-dimensional amplitude of the wavy surface ranging from 0.1 to 0.3. Later on, we shall extend the analysis by considering several values of the phase shift between the wavy surfaces. These results will provide a base reference for evaluation of the nanofluid flows thermal performance. In Fig. 5 we show the obtained local Nusselt number distributions along the middle region of the channel for different values of the wavy surface amplitude, with a comparison between our numerical solution and the one reported in [62] for the case a = 0.2.  The local Nusselt number is expected to be higher in the converging section of each wave than in the diverging section. This can be inferred by considering that the converging sections are characterized by higher velocity gradients and, therefore, by higher temperature gradients, which in turn leads to enhanced heat transfer rate. This behavior is correctly captured by our model and, despite a slight overestimation of the local maximum values of Nu toward the end of the middle region of the channel, our numerical results are in good agreement with those from literature. Also, we observe a significant increase of the local maximum value of Nu, at each converging section, when the wavy surface amplitude is varied from 0.1 to 0.3. Conversely, the local minimum values of Nu, occurring at the diverging sections of the waves, decrease when a goes from 0.1 to 0.2, but they remain fairly the same when a is varied further, namely up to 0.3. These findings are in line with those provided in [62], for the same configurations, obtained for different values of Re. In Table 3 we report the computed mean Nusselt number for the three cases considered.
The mean Nusselt number values obtained via the adopted method match well with those reported by Wang and Chen [62]. Further, we show in Fig. 6 the computed local Nusselt number distributions along the middle region of the channel when a phase shift between top and bottom wavy walls is set. In this case we have kept the non-dimensional amplitude of both wavy surfaces equal to 0.2.
We note that a phase shift between the walls leads to a corresponding shift of the local Nusselt number profile along the channel. Moreover, the amplitude of the Nusselt number decreases when the phase shift ϕ is varied from 0 to π/2. In particular, higher values of ϕ lead to lower values of the local maximum Nusselt number, at the converging sections, while the local minimum values of Nu, at diverging sections, remain almost unchanged. The case of ϕ = π and ϕ = 0, which corresponds to antiphase wavy surfaces, reveals instead a different behavior. In such a configuration the cross section of the channel is constant, therefore the flow is subjected to low velocity gradient variations. This results in a lower Nusselt number along all the middle region of the channel then in the previous cases. Table  4 summarizes the mean Nusselt numbers as obtained for the analyzed configurations.
The results show a significant impact of the phase shift parameters on the heat transfer performance of the water flow, with a drop of the mean Nusselt number when the two wavy surfaces are in the antiphase configuration.

Forced Convection of Nanofluids in Wavy-Wall Channels
In this Section, we analyze the flow and heat transfer of two metal oxide nanofluids, namely water-alumina (Al 2 O 3 ) and water-copper oxide (CuO), by considering variable concentration of solid nanoparticles into the base fluid, and different configurations of wavy-wall channel. In particular, we first investigate on the effects of the wavy surface amplitude and, subsequently, we analyze the effects due to a phase shift between top and bottom wavy surfaces.

Effect of Wavy Surface Amplitude
We vary the non-dimensional amplitude of both wavy surfaces from 0.1 to 0.3, without any phase shift. In Fig. 7 we show the obtained distributions for the local Nusselt number along the middle region of the channel, as for the case of Al 2 O 3 -water nanofluid. In particular, we consider several nanoparticles volume fractions.
The local Nusselt number profiles for the considered nanofluid follow those related to water (φ = 0), as expected. However, we observe a significant enhancement in heat transfer with the increasing of the nanoparticles volume fraction, for all of the three cases analyzed. In fact, the local reduction of the temperature gradients occurring within the nanofluid flow is compensated by an increase of the thermal conductivity, thus leading to an overall increasing of the heat transfer rate. We obtain similar trends to those depicted in Fig. 7 also for the CuO-water nanofluid. Nevertheless, we emphasize that the effect due to the presence of nanoparticles within the base fluid is more pronounced in Al 2 O 3 -water nanofluid, which exhibits the most enhanced thermal behavior among the two considered nanofluids. This result can be justified by considering that the Al 2 O 3 -water nanofluid is characterized by higher thermal conductivity an higher thermal diffusivity than the CuO-water nanofluid [65,66].
Next, Fig. 8 shows a direct comparison between the local Nusselt distribution obtained for water and for Al 2 O 3 -water nanofluid with an elevated concentration of nanoparticles, that is φ = 0.10, in three different channel configurations, namely those with a equal to 0.1, 0.2 and 0.3.
Our results indicate that for wavy surface amplitudes equal to 0.2 and 0.3 the local maximum values of Nu, occurring at the converging section of each wave, are significantly higher for Al 2 O 3 -water nanofluid than those for pure water, while the local miminum values of Nu, at the diverging sections, are nearly the same. This leads to an overall enhancement of the heat transfer performance of the fluid flow. Further, Fig. 9 reports the variation of the local Nusselt number along the middle region of a wavy channel characterized by a = 0.2, for both the considered nanofluids, with φ = 0.10, in comparison with water.
As mentioned before, in the present flow regime the Al 2 O 3 -water nanofluid is the one exhibiting the highest values of Nu, along the channel. Figure 10 illustrates the temperature fields for the CuO-water nanofluid flow, at steady-state conditions, as a representative example of the flow through the wavy-wall channel, in the three different geometrical configurations so far investigated. From Fig. 10 we observe that an increasing of the wavy surface amplitude causes a significant difference in the temperature distribution inside the channel. In particular, in the a = 0.1 configuration we note that the thickness of the thermal boundary layer is quite regular along the channel and the isothermal lines are contracted near the walls. On the other hand, at higher values of the wavy surface amplitude, a more intense mixing of the fluid flow occurs at the diverging sections of the channel.
In fact, results show that a flow reversal occurs in the regions near the wall surfaces enclosed by each wave. This phenomenon becomes evident at high values of the wavy surface amplitude (a = 0.3), while it is almost negligible at lower values (a = 0.1): as the wavy  amplitude is increased, the strength of the flow reversal also increases. In order to depict this behavior, in Fig. 11, where only a portion of the domain is represented, we report the streamlines contours for the case of CuO-water nanofluid flow with φ = 0.10 in a channel configuration characterized by wavy surface amplitude a = 0.3.
Such a recirculating flow, characterized by low velocity gradients, causes the mixing of the fluid flow, which decreases the local convective heat transfer coefficient. This is the reason why the local minimum values of Nu, which occur in the proximity of the diverging sections of the channel, are higher in the case of a = 0.1 than those for the channel with a = 0.2 and a = 0.3. The presence of nanoparticles within the flow however, leads to an enhancement of   From Fig. 12 we note that the nanoparticles volume fraction plays a major role in the characterization of the thermal performance of the fluid flow. Specifically, increasing the volume fraction of the nanoparticles beyond 0.02 results in a significant enhancement in the heat transfer rate. This stands for all the geometrical configurations analyzed and it is especially true for the Al 2 O 3 -water nanofluid. For this specific nanofluid, the computed enhancement factor (Eq. 29) goes from a minimum value of 2.6%, obtained for the case of a = 0.3 with φ = 0.02, to a maximum value of 19.9%, related to the case of a = 0.2 with φ = 0.10. The enhancement factor obtained for the CuO-water nanofluid flow is generally lower, however, it is still substantial, with a maximum value of around 13%. Finally, Fig.  13 displays the variation of the non-dimensional pressure drop with nanoparticles volume fractions for different wavy surface amplitudes, for the case of Al 2 O 3 -water nanofluid. In Fig. 13, the values reported for the case of a = 0 correspond to those from the analytical solution of Navier-Stokes equations for the plane Poiseuille flow. As expected, for each channel configuration, higher values of the nanoparticles volume fraction correspond to higher pressure drops. This is due to the increase of density and dynamic viscosity of nanofluid with respect to the water. Also, we note that the channel with wavy surface amplitudes a = 0.3 has the highest pressure drop over all the range of nanoparticles volume fraction considered, followed by the channels with a = 0.2 and a = 0.1. This indicates that, despite high values of the wavy surface amplitude (up to a = 0.3) lead to an enhancement of the heat transfer rate, they cause a significant increment in pumping power of the channel.

Effect of Phase Shift
This section reports the effects of the phase shift between the two wavy walls for Al 2 O 3 -water and CuO-water nanofluid flows in four different channel configurations. These are defined through the two phase shift parameters previously introduced, namely ϕ and ϕ in Eqs. (16) and (17). The non-dimensional amplitude of both wavy surfaces is set to 0.2. Figures 14 and  15 report the distribution of the local Nusselt number along the middle region of the channel  The phase shift between the wavy walls leads to a shifting of the local Nusselt number profile as well. Moreover, the amplitude of Nu decreases with the increasing of the asymmetry of the channel, and it drops when the two wavy surfaces are in antiphase configuration. Both the Al 2 O 3 -water and CuO-water nanofluids outperform the pure water flow in terms of heat transfer. In fact, we observe that, in both cases, the presence of nanoparticles leads to an enhancement of the local maximum values of Nu along the channel. Figure 16 illustrates the temperature fields for the CuO-water nanofluid flow with φ = 0.10 for several channel configurations, as an illustrative example. As depicted in Fig. 16, the phase shift between the wavy walls gives rise to an asymmetric temperature distribution with respect to the central axis of the channel. Such a phenomenon is responsible for a significant decreasing of the heat transfer with respect to the case of symmetric channel configuration. Finally, Fig. 17 shows the effect of nanoparticles volume fraction on the mean Nusselt number for different configurations of the wavy-wall channel and for both the considered nanofluids.
The mean Nusselt number varies significantly across the different channel configurations, with the axis-symmetric geometry being the one that exhibit the best thermal performance. Also, the results indicate that the heat transfer may be significantly enhanced with addition of nanoparticles. Specifically, the enhancement factor reaches 20% for Al 2 O 3 -water nanofluid with φ = 0.10 in the antiphase wavy surfaces configuration, thus suggesting that the employment of such a fluid may represent an effective solution in heat transfer devices. As far as the pressure drop it concerns, we obtain quite similar results among all the considered geometries, for each nanoparticles volume fraction, with the anthiphase configuration being the one exhibiting slightly lower values of p with respect to the other cases. This result indicates that a phase shift between the wavy walls does not bring any significant advantage in terms of pumping power.
The presented investigation provides a quantitative analysis of the nanofluid heat transfer in wavy-wall channels of different geometries. The wavy-wall channel configuration represents indeed an interesting option to generate flow mixing and enhance the heat transfer rate in several engineering applications. Among these, for instance, microscale cooling devices, such as microchannel heat sinks, are of particular relevance for heat removal. In these cases, the flow generally occurs at relatively low Reynolds numbers, therefore an optimal design of the channel, in conjunction with the use of a nanofluid, is instrumental to achieve improved cooling performance. In this context, our study may provide useful information on the thermal behavior of nanofluid flows through corrugated walls, by proving that the irregularity of the channel surfaces may lead to an enhancement or a decreasing of the thermal performance depending on the specific geometrical arrangement. The method used in this work finds therefore applicability in the development and optimization of advanced heat transfer devices involving laminar flows with improved thermal properties.

Conclusions
A numerical investigation on the forced convection of nanofluid flows through wavy-wall channels has been conducted in this work. The employed computational approach is based on a finite-volume formulation of the lattice Boltmann method, which is applied to fully unstructured meshes. This model allows to deal with complex geometries, thus overcoming a severe limitation of traditional LBMs (i.e. LBMs on uniform Cartesian mesh), and extending its range of applicability toward the study of realistic problems of engineering interest.
In this paper, we consider two types of water-based nanofluids, that are: aluminum oxide and copper oxide nanofluids. These are modeled by following a single phase homogeneous approach. We then study the effects due to geometrical parameters such as the wavy surface amplitude and the phase shift between the wavy surfaces.
As a result, we observe that increasing the wavy surface amplitude (up to a nondimensional value of 0.3) leads to an improvement of the thermal performances of the fluid flow, in terms of heat transfer rate. In contrast, a phase shift between the wavy surfaces causes a decrease of the local Nusselt number, which results in a lower heat transfer rate. These findings suggest that the axis-symmetric configuration of the channel, with a high value of the wavy surface amplitude, represents the most suitable solution to enhance the heat transfer rate. In addition, we find that the employment of nanofluids leads to a significant increase of the local Nusselt number, especially when the nanoparticles volume fraction exceed the 5%. Among the different nanofluids considered, the Al 2 O 3 -water nanofluid provides the most effective way to improve thermal performances, by leading to an enhancement factor up to 20%.