Received 11 April 2016; accepted 18 June 2016; published 21 June 2016

ABSTRACT

In the present work, we numerically study the laminar natural convection of a nanofluid confined in a square cavity. The vertical walls are assumed to be insulated, non-conducting, and impermeable to mass transfer. The horizontal walls are differentially heated, and the low is maintained at hot condition (sinusoidal) when the high one is cold. The objective of this work is to develop a new height accurate method for solving heat transfer equations. The new method is a Fourth Order Compact (F.O.C). This work aims to show the interest of the method and understand the effect of the presence of nanofluids in closed square systems on the natural convection mechanism. The numerical simulations are performed for Prandtl number (), the Rayleigh numbers varying between and for different volume fractions varies between 0% and 10% for the nanofluid (water + Cu).

Keywords:

Nanofluid, Heat Transfer, Natural Convection, Fourth-Order Compact (F.O.C) Formulation, Numerical Performance, Sinusoidal Boundary Thermal Condition

1. Introduction

In a better description, nanofluids are engineered colloidal suspensions of nanoparticles (1 - 100 nm) in a base fluid. Common base fluids include water, oil, and Ethylene Glycol while nanoparticles are typically made of chemically stable metals, metal oxides or carbon in various forms. The use of particles of nanometer dimension was first continuously studied by a research group at the Argonne National Laboratory a decade ago. S. Choi [1] in 1995 was probably the first one who called the fluids with particles of nanometer dimensions “nanofluids”. He showed substantial augmentation of heat transported in suspensions of copper or aluminum nanoparticles in water and other liquids. Compared with suspended particles of millimeter or micrometer dimensions, nanofluids show better stability and rheological properties, dramatically higher thermal conductivities and no penalty in pressure drop. Several published literatures have mainly focused on the prediction and measurement techniques in order to evaluate the thermal conductivity of nanofluid. It is noticeable that only a few papers have discussed the convective heat transfer of nanofluids, including the experimental and theoretical investigation.

A numerical study of natural convection of copper-water nanofluide in a two dimensional enclosure was conducted by Khanafer et al. [2] . The nanofluid in the enclosure was assumed to be in single phase. It was found in any given Grashof number, heat transfer in the enclosure increased with the volumetric fraction of the copper nanoparticles in water. Lee et al. [3] measured the thermal conductivity of Al₂O₃ water and Cu-water nanofluids and indicated that the thermal conductivity of nanofluids increases with solid volume fraction. He concluded that any new models of nanofluide thermal conductivity should contain the effect of surface area and structure dependent behavior as well as the size effect. Xie et al. [4] added spherical and cylindrical shaped nano sized SiC particles to water and Ethylene Glycol, separately and found that cylindrical nanoparticles increased thermal conductivity more than spherical ones. The dependence of thermal conductivity of nanoparticles-fluid mixture was estimated by Xie et al. [5] . Some of the theoretical and experimental studies have been reported on convective heat transfer coefficient [6] - [9] .

Sandeep Naramgari and C. Sulochana [10] analyzed the momentum and heat transfer behavior of MHD nanofluid embedded with conducting dust particles past a stretching surface in the presence of volume fraction of dust particles. They solved equations numerically using Runge-Kutta based shooting technique and showed that the increase in the interaction between the fluid and particle phase enhanced the heat transfer rate and reduced the friction factor. Nader Ben-Cheikh et al. [11] studied natural convection in a square enclosure filled with a water based nanofluid (water with Ag, Cu, Al₂O₃ or TiO₂ nanoparticles) with non-uniform (sinusoidal) temperature distribution maintained at the bottom wall. An accurate finite volume scheme along with a multi-grid technique is devised for the purpose of solution of the governing equations. Tiwari et al. [12] numerically investigated the behavior of copper-water nanofluid in a two sided lid-driven differentially heated cavity. They considered different cases characterized by the direction of movement of walls and found that both the Richardson number and direction of moving walls influenced the fluid flow and thermal behavior. Yadil et al. [13] study the Cu/Water nanofluids filled baffled square cavity.

The effects of Rayleigh number, volume fraction and partitions location on the average Nusselt number are studied. I.El Bouihi and R. Sehaqui [14] she simulate the flow features of nanofluids for a range of solid volume fraction χ and a sinusoidal thermal boundary condition, and we obtained correlations of heat transfer in enclosures for two different thermal boundary conditions on the left wall. Ami-Nossadati and Ghasemi [15] studied natural convection cooling of a localized heat source at the bottom of a nanofluid filled enclosure. Ogut [16] investigated natural convection of water-based nanofluids in an inclined enclosure with a heat source using the expression for calculating the effective thermal conductivity of solid-liquid mixtures proposed by Yu and Choi [17] . Ghasemi and Aminossadati [18] considered periodic natural convection in a nanofluid filled enclosure with oscillating heat flux. Non-uniform heating of surfaces in buoyancy driven flow in a cavity has significant effect of the flow and heat transfer characteristics and finds applications in various areas such as crystal growth in liquids, energy storage, geophysics, solar distillers and others. In a relatively recent study, Sarris et al. [19] reported that the sinusoidal wall temperature variation produced uniform melting of materials such as glass in their detailed study on the effect of sinusoidal top wall temperature variations in a natural convection within a square enclosure where the other walls are insulated. Corcione [20] studied natural convection in an air filled rectangular enclosure heated from below and cooled from above for a variety of thermal boundary conditions at the side walls. Roy and Basak [21] studied numerically natural convection flows in a square cavity with non-uniformly (sinusoidal) heated wall(s) using the finite element method. The bottom wall and one vertical wall were heated (uniformly and non-uniformly) and the top wall was insulated while the other vertical wall was cooled by means of a constant temperature bath. Sathiyamoorthy et al. [22] investigated steady natural convection flows in a square cavity with linearly heated side wall(s).

In order to optimize and improve heat transfer by natural convection in closed square cavity. Although extensive research has been given to cases of rectangular cavity filled nanofluid, few studies have focused on the study, theoretical or numerical discretizations high order (≥2). The numerical study of systems of equations in the heat transfer area is usually treated by various methods, sometimes numérical classics like Finite Elements (FE), Volume Finite (VF) and Finite Differences (DF). or using some software adapted as “FLUENT”. To solve fluid mechanics problems such as conductive heat transfer, convective or mixed into regular geometries.

More specific order four schemes have been used to solve the Navier-Stokes in enclosures without considering the energy equation by Ecran Erturk et C. Gokcol [2] . The objective of this work is to develop a new method for solving heat transfer equations in convection. The new method is a Fourth Order Compact. This work aims to show the interest of the method and understand the effect of the presence of nanofluids in closed square systems on the natural convection mechanism.

2. Mathematical Formulation

Consider a square cavity filled with a nanofluid. The vertical walls are assumed to be insulated, non-conducting, and impermeable to mass transfer. The horizontal walls are differentially heated, the low is maintained at hot condition (sinusoidal) when the high one is cold (Figure 1). The nanofluid in the enclosure is Newtonian, incompressible, and laminar. The nanoparticles are assumed to have a uniform shape and size. Moreover, it is assumed that both the fluid phase and nanoparticles are in thermal equilibrium state and they flow at the same velocity. The thermophysical properties of the nanofluid are assumed to be constant except for the density variation in the buoyancy force, which is based on the Boussinesq approximation.

We have considered the continuity, momentum and energy equations for a Newtonian, Fourier constant property fluid governing an unsteady, two-dimensional flow. It is further assumed that radiation heat transfer among sides is negligible with respect to other modes of heat transfer. Under the assumption of constant thermal properties, the Navier-Stokes equations for an unsteady, incompressible, two-dimensional flow are:

Continuity equation:

(1)

x-momentum equation:

(2)

y-momentum equation:

(3)

Energy equation:

(4)

where

The viscosity of the nanofluid can be estimated with the existing relations for the two-phase mixture. The equation given by Brinkman [23] has been used as the relation for effective viscosity in this problem, as given by Xuan and Li [24] have experimentally measured the apparent viscosity of the transformer oil-water nanofluid and of the water-copper nanofluid in the temperature range of 20˚C - 50˚C. The experimental results reveal relatively good agreement with Brinkman’s theory. The thermophysical properties of fluid and the solid phases are shown in Table 1.

The effective density of the nanofluid at reference temperature is

(5)

The heat capacitance of the nanofluid is expressed as Abu-Nada [25] and Khanafer et al. [2]

(6)

The effective thermal conductivity of the nanofluid is approximated by the Maxwell-Garnetts model [26]

(7)

Equations (1)-(4) can be converted to the dimensionless forms by definition of the following parameters as:

Hence, the governing equations of continuity, linear momentum and energy for unsteady laminar flow in Cartesian coordinates take the following dimensionless form:

(8)

(9)

(10)

(11)

The enclosure boundary conditions consist of no-slip and no penetration walls, U = V = 0 on all four walls. The thermal boundary conditions on the bottom wall is such that. The left and right vertical walls are assumed to be insulated, non-conducting, and impermeable to mass transfer

and the bottom wall are at the cold temperature

The governing equations for the present study in formulation taking into the account the above mentioned assumptions are written in dimensionless form as:

Kinematics equation:

(12)

Vorticity equation:

(13)

where:

Energy equation:

(14)

Before turning to the application of the method of fourth order in the equations governing our problem we will combine Equations (12) and (13) in a condensed form by introducing a dummy variable, which replace either the temperature is the vorticity.

(15)

All these terms are listed in Table 2.

Dimensionless boundary conditions for are:

For: and,

For: and,.

For vorcicite Störtkuh et al. [27] have presented an analytical asymptotic solution near the corners of cavity and using finite element bilinear shape functions they also have presented a singularity-removed boundary condition for vorticity at the corner points as well as at the wall points. For the boundary conditions, in both of the numerical methods described above we follow Störtkuh et al. [27] and use the following expression for calculating vorticity values at the wall.

(16)

For corner points, we again follow Störtkuh et al. [27] and use the following expression for calculating the vorticity values:

(17)

where V is the speed of the wall in our case which is equal to 0 for the four stationary walls.

In explicit notation, for the wall points shown in Figure 2(a), the vorticity is calculated as the following:

(18)

Similarly, for the corner points also shown in Figure 2(b), the vorticity is calculated as the following:

(19)

The reader is referred to Störtkuh et al. [27] for details on the boundary conditions.

The local and averaged heat transfer rates at the bottom hot wall of the cavity are presented by means of the local and averaged Nusselt numbers, Nu and, which are, respectively determined as follows:

(20)

(21)

3. Method of the Fourth Order Compact

3.1. Introduction

High-Order Compact (HOC) formulations are becoming more popular in computational fluid dynamics (CFD) field of study. Compact formulations provide more accurate solutions in a compact stencil. In finite differences, a standard three-point discretization provides second-order spatial accuracy and this type of discretization is very widely used. When a high-order spatial discretization is desired, i.e. fourth-order accuracy, then a five- point discretization has to be used. However, in a five point discretization there is a complexity in handling the points near the boundaries. High-order compact schemes provide fourth-order spatial accuracy in a 3 × 3 stencil

and this type of compact formulations does not have the complexity near the boundaries that a standard wide (five- point) fourth-order formulation would have. Dennis and Hudson [28] , MacKinnon and Johnson [29] , Gupta et al. [30] , Spotz and Carey [31] , Li et al. [32] , have demonstrated the efficiency of the HOC schemes on the stream function and vorticity formulation of two dimensions. In the literature, it is possible to find numerous different types of iterative numerical methods for the momentum equations. These numerical methods, however, could not be easily used in HOC schemes because of the final form of the HOC formulations used in References [28] - [32] . This fact might be counted as a disadvantage of HOC formulations that the coding stage is rather complex due to the resulting stencil used in these studies. It would be very useful if any numerical method for the solution of momentum equations described in books and papers could be easily applied to HOC formulations. E. Erturk and C. Gokcol [33] present a new Fourth-Order Compact Formulation. The difference of this formulation with References [28] - [32] is not in the way that the Fourth-Order Compact scheme is obtained. The main difference, however, is in the way that the final forms of the equations are written. The main advantage of this formulation is that, any iterative numerical method used for Navier-Stokes equations, can be easily applied to this new FOC formulation, since the final form of the presented FOC formulation is in the same form with the Navier?Stokes equations. Moreover, if someone already has a second-order accurate code for the solution of conservation equations the mass and momentum, they can easily convert their existing code to fourth- order accuracy by just adding some coefficients into their existing code. In this study, we will applied Using this new compact formulation, we have solved the conservation equations of mass, momentum and energy in square cavity at the Rayleigh number varies in the range and for different solide volume fractions of nanoparticles (Cu) is varied as, taking water as a base fluid with a Prandtl number equal to using a very fine grid mesh to demonstrate the efficiency of this new formulation.

3.2. Principle Method of Fourth Order Compact

We will use the equations of streamlines, vorticity and energy dimensionless forms are given as follows:

Stream function:

(22)

General equation of conservation:

(23)

For first and second-order derivatives the following discretizations are fourth-order accurate:

(24)

(25)

where and are standard second-order central discretizations such that

(26)

(27)

If we apply the discretizations in Equations (24) and (25) to Equations (22) and (23), we obtain the following equation

(28)

(29)

In these equations we have third and fourth derivatives of stream function and general equation of conservation bringing together the equations of vorticity and energy. In order to find an expression for these derivatives we use Equations (22) and (23).

For example, when we take the first and second x-derivative of the stream function Equation (22) we obtain

(30)

(31)

And also, by taking the first and second y-derivative of the stream function Equation (22) we obtain

(32)

(33)

Using standard second-order central discretizations given in Table 3, these equations can be written as

(34)

(35)

(36)

(37)

When we substitute Equations (35) and (37) into Equation (28) we obtain the following finite difference equation.

(38)

We note that the solution of Equation (38) is also a solution to stream function Equation (22) with fourth-or- der spatial accuracy. Therefore, if we numerically solve Equation (38), the solution we obtain will satisfy the stream function equation up to fourth order accuracy.

In order to obtain a fourth-order approximation for the vorticity equation and energy (23), we follow the same procedure. When we take the first and second derivatives of the general equation of conservation (23) with respect to x- and y-coordinates we obtain:

(39)

(40)

(41)

(42)

If we substitute Equations (39) and (41) for the third derivatives of the general equation of conservation and into Equations (29), (40) and (42) and also if we substitute Equations (34) and (36), for the third derivatives of stream function into Equations (29), (40) and (42) and finally, if we substitute Equations (40) and (42) for the fourth derivatives of the general equation of conservation into Equation (29), then we obtain the following:

(43)

where:

and

Again we note that the solution of Equation (43) satisfy the vorticity and energy Equation (23) with fourth- order accuracy.

As the final form of our FOC scheme, we prefer to write Equations (38) and (43) as

(44)

(45)

where

(46)

We note that the finite difference Equations (44) and (45) are fourth-order accurate approximation of the stream function, vorticity and energy Equations (22) and (23). In Equations (44) and (45), however, if A, B, C, D, E and F are chosen to be equal to 0 then the finite difference Equations (44) and (45) simply become

(47)

(48)

Equations (47) and (48) are the standard second-order accurate approximation of the stream-

function and the general equation of conservation (22) and (23). When we use Equations (44) and (45) for the numerical solution of the stream function and general equation of conservation, we can easily switch between second and fourth-order accuracy just by using homogeneous values for the coefficients A, B, C, D, E and F or by using the expressions defined in Equation (46) in the code.

We note that the numerical solutions of Equations (44) and (45), strictly provided that second-order discretizations in Table 3 are used and also strictly provided that a uniform grid mesh with and is used, are fourth-order accurate to streamfunction and the general equation of conservation (21) and (22). The only difference between Equations (44) and (45) and Equations (22) and (23) are the coefficients A, B, C, D, E and F. So these equations are of the same form, therefore, all the iterative numerical methods (such as SOR, ADI, factorization schemes, pseudo time iterations, etc.) used to solve stream-function, vorticity and energy Equations (22) and (23) can also be easily applied to fourth-order Equations (44) and (45). In our work we apply the ADI method on the equations of 4th order.

As a measure of convergence to the steady state, during the iterations we monitored three residual parameters. The first residual parameter, RES1, is defined as the maximum absolute residual of the finite difference equations of steady stream function and general Equations (44) and (45). These are, respectively, given as

(49)

(50)

The magnitude of RES1 is an indication of the degree to which the solution has converged to steady state. In the limit RES1 would be zero. The second residual parameter, RES2, is defined as the maximum absolute difference between two iteration steps in the stream function, vorticity and energy variables. These are, respectively, given as

(51)

(52)

RES2 gives an indication of the significant digit on which the code is iterating. The third residual parameter, RES3, is similar to RES2, except that it is normalized by the representative value at the previous time step. This then provides an indication of the maximum percent change in and in each iteration step. RES3 is defined as

(53)

(54)

In our calculations, for all Rayleigh numbers we considered that convergence was achieved when both and were achieved. Such a low value was chosen to ensure the accuracy of the solution. At these convergence levels, the second residual parameters were in the order of and, which means that the stream function, vorticity and energy variables are accurate to the 10^th and 9^th digit accuracy, respectively, at a grid point and even more accurate at the rest of the grids. In addition, at these convergence levels the third residual parameters were in the order of and, which means that the stream function, vorticity and energy variables are changing with 10% - 9% and 10% - 8% of their values, respectively, in an iteration step at a grid point and even with less percentage at the rest of the grids. These very low residuals ensure that our solutions are indeed very accurate.

4. Results and Discussion

In the present grid independence test, the Prandtl number is set to Pr = 6.2 (pure water). The nanoparticles are chosen to be copper (Cu) with a solid volume fraction and a Rayleigh number Ra = 10⁵. Numerical computations have been carried on five different grid sizes 32 × 32, 42 × 42, 62 × 62, 82 × 82 and 102 × 102 grid sizes. Table 4 regroups the values of the averaged Nusselt number through the hot wall and the maximum value of the stream function. Uniform grid has been used for all the computations. The distribution of the u-velocity in the vertical mid-plane and v-velocity in the horizontal mid-plane are shown in Figure 3. It is observed that the curves overlap with each other for 82 × 82 and 102 × 102. So a grid number of 82 × 82 is chosen for further computation.

Our code has been tested for natural convection fluid flows in differentially heated cavities and in Rayleigh- Bénard configuration for Rayleigh numbers between 10³ and 10⁶ (Table 5) and gave excellent results (see ref. [12] [14] [34] - [36] ).

In this section, the nanofluid-filled enclosure is studied for a range of solid volume fraction and the Rayleigh number varies from 10³ to 10⁵. For all simulations the considered base fluid is water (Pr = 6.2). In Figure 4, we present the streamlines (top) and isotherms (bottom) for 10³ ≤ Ra ≤ 10⁵, for the case of a water- Cu nanofluid and pure water. The value of solid volume fraction is set to. Figure 5 represents the same physical quantities but for a volume fraction value of. Due to the temperature distribution imposed at the bottom wall and to the boundary conditions on vertical walls, we observe a symmetry behavior in both the streamlines and in the contour maps of the isotherms. We can see that whatever the Rayleigh number and value of solid volume fraction, the flow is mainly composed of two counter-rotating circulating cells.

Figure 6 presents the velocity profiles V(X) and U(Y) along the mid-section of the enclosure X = 0.5 and Y = 0.5 for different values of χ and is in good concordance with the fact that the nanofluid moves slower than a pure water. Indeed, for Ra = 10³, the deviation (relative to χ = 0) between the maximum values of vertical velocity is for χ = 4% and for χ = 10%. By increasing the Rayleigh number, these deviations decrease. For Ra = 5 × 10⁴ and Ra = 10⁵ these deviations are, , and respectively. As far as the temperature distribution is concerned, clear differences are observed in the isotherm contour plots compared to the case χ = 0. These differences are accentuated as the solid volume fraction increases. These differences mean that the presence of nanoparticles affect especially the heat transfer rate through the enclosure.

The heat transfer distribution through the hot wall is displayed in Figure 7 through the plotted lines of the local Nusselt number for different values of Ra. One can see that for all combinations of Ra and χ, the local Nusselt number behavior is symmetric with respect to the plane X = 0.5. For low Rayleigh numbers, (Ra = 5 × 10³) and χ = 0, the transfer of heat through the hot wall is relatively low with a slight curvature at X = 0.5. This curvature is due to relatively higher intensity of the counter rotating cells represented by the highest value of when χ = 0. When χ increases to χ = 0.1, the curvature at the center disappears because the fluid velocity decreases. The heat transfer in this case is maximum at X = 0.5 and is higher due to the presence of nanoparticles whose thermal conductivity is much greater than that of water. The same phenomena are observed almost on the curves related to Ra = 5 × 10⁴ and Ra = 5 × 10⁵ with a maximum heat transfer in the vicinity of X = 0.25 and X = 0.75. For example, for Ra = 5 × 10⁴, the maximum Nusselt number value is Nu_max = 7.374 and is situated at both locations X = 0.317 and X = 0.682 for χ = 0. For χ = 0.1, Nu_max = 10.394 and is located at X = 0.317 and X = 0.695.

The variations of average Nusselt number (Nu) with Ra and χ are shown in Table 6. For Ra = 10³, there is a substantial increase in Nu as χ is increased above 2%. In general, Nu increases with χ. When χ is 2%, the increase is approximately 9.75%. When χ is 4%, the increase is approximately 19.51%. When χ is 8%, the increase is above 42.07%. For Ra = 10⁴, as χ is increased to 2%, an increase of 6.03% is observed. A heat transfer augmentation of above 27.63% is obtained for compared to or more is observed for different Rayleigh number Ra. Thus, one can conclude that the Nusselt number increases with the increase of the volume fraction χ and Rayleigh number Ra.

5. Conclusions

In this study the heat transfer enhancement in a two dimensional enclosure filled with nanofluids is studied

Figure 4. Streamlines and isotherms for the enclosure filled with water-Cu nanofluid (- - -), and pure water (―) at different Rayleigh numbers.

Figure 5. Streamlines and isotherms for the enclosure filled with water-Cu nanofluid (- - -), and pure water (―) at different Rayleigh numbers.

numerically. This study presented a new fourth-order compact formulation and investigated the effect of a sinusoidal thermal boundary condition, for different Rayleigh number Ra and volume fractions of nanoparticles. The flow and temperature fields are symmetric near the middle plane of the enclosure due to the imposed symmetry condition on the bottom wall boundary. From the results of this work, the following main conclusions may be drawn:

・ The fourth-order accurate compact formulation was developed and was in agree with previous studies.

・ Our numerical code has been validated for different Rayleigh number.

・ A comparative study illustrates that the suspended nanoparticles substantially increase the heat transfer rate with an increase in the nanoparticles volume fraction for different Rayleigh number Ra Rayleigh number. Moreover, the nanofluid flows as well as the cooper nanoparticles increase.

In the near future, this study will be extended for different geometry studies and other types of base fluids and nanoparticles.

Cite this paper

Mostafa Zaydan,Naoufal Yadil,Zoubair Boulahia,Abderrahim Wakif,Rachid Sehaqui, (2016) Fourth-Order Compact Formulation for the Resolution of Heat Transfer in Natural Convection of Water-Cu Nanofluid in a Square Cavity with a Sinusoidal Boundary Thermal Condition. World Journal of Nano Science and Engineering,06,70-89. doi: 10.4236/wjnse.2016.62009

References

Nomenclature

i x-direction grid point

j y-direction grid point

C_p Specific heat capacity (J・K⁻¹)

g Gravitational cceleration (m・s⁻²)

h ocal heat transfer coefficient (m⁻²・K⁻¹)

H Height of cavity (m)

q_w Heat flux (W・m⁻²)

t Dimensional time (s)

τ Non-dimensional time

T Temperature (K)

p pressure (pa)

(x, y) Dimensionless Cartesian coordinates (m)

u, v velocity components in x, y directions (m・s⁻¹)

Greek-Symbols

α Fluid thermal diffusivity (m²・s⁻¹)

β Thermal expansion coefficient (K⁻¹)

ν Kinematic viscosity (m²・s⁻¹)

ρ Density (kg・m⁻³)

μ Dynamic viscosity (N・s・m⁻²)

κ Thermal conductivity (W・m⁻¹・K⁻¹)

ψ dimensional stream function

θ non-dimensional temperature

ω dimensionalvorticity

χ nanoparticle volume fraction

Subscripts

c cold wall

eff effective

h hot wall

s solid

f pure fluid

Submit your manuscript at: http://papersubmission.scirp.org/