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

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 ( Pr = 6.2 ), the Rayleigh numbers varying between Ra ≤ ≤ × 3 5 10 5 10 and for different volume fractions χ varies between 0% and 10% for the nanofluid (water + Cu).


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 2 O 3 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 2 O 3 or TiO 2 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.

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: x-momentum equation: ,0 y-momentum equation: Energy equation: ( ) where ( ) 2 and .
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 ( ) The heat capacitance of the nanofluid is expressed as Abu-Nada [25] and Khanafer et al. [2] ( ) ( ) ( )( ) The effective thermal conductivity of the nanofluid is approximated by the Maxwell-Garnetts model [26] ( ) ( ) ( ) ( ) Equations ( 1)-( 4) can be converted to the dimensionless forms by definition of the following parameters as: ; and Hence, the governing equations of continuity, linear momentum and energy for unsteady laminar flow in Cartesian coordinates take the following dimensionless form: Table 1.Thermophysical properties of water and nanoparticles.

Physical properties
Pure water Cu ( ) ( ) 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 men- tioned assumptions are written in dimensionless form as: Kinematics equation: Vorticity equation: where: ( ) 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 ω .
All these terms are listed in Table 2.

Dimensionless boundary conditions for ( )
, ψ ω are: 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.

Quantities transported
For corner points, we again follow Störtkuh et al. [27] and use the following expression for calculating the vorticity values: 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: Similarly, for the corner points also shown in Figure 2(b), the vorticity is calculated as the following: ( ) 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 Nu , which are, respectively determined as follows:

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 fivepoint 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 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 ( )

2
O x ∆ code for the solution of conservation equations the mass and momentum, they can easily convert their existing code to fourthorder accuracy ( )

O x
∆ 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 , taking water as a base fluid with a Prandtl number equal to ( )

Pr =
using a very fine grid mesh to demonstrate the efficiency of this new formulation.

Principle Method of Fourth Order Compact
We will use the equations of streamlines ψ , vorticity ω and energy θ dimensionless forms are given as follows: Stream function: x y General equation of conservation: For first and second-order derivatives the following discretizations are fourth-order accurate: ( ) ( ) where x Φ and xx Φ are standard second-order central discretizations such that 1 1 If we apply the discretizations in Equations ( 24) and (25) to Equations ( 22) and ( 23), we obtain the following equation ( ) 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 x x x y And also, by taking the first and second y-derivative of the stream function Equation ( 22) we obtain Using standard second-order central discretizations given in Table 3, these equations can be written as ( ) ( ) ( ) , When we substitute Equations ( 35) and (37) into Equation ( 28) we obtain the following finite difference equation.

Derivations Discretizations
i j i j i j i j i j i j x y i j i j i j i j i j i j i j i j i j x y We note that the solution of Equation ( 38) is also a solution to stream function Equation ( 22) with fourth-order 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: 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: where: ; .
Again we note that the solution of Equation ( 43) satisfy the vorticity and energy Equation ( 23) with fourthorder accuracy.
As the final form of our FOC scheme, we prefer to write Equations ( 38) and (43) as   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 Equations ( 47) and ( 48) are the standard second-order accurate , O x y ∆ ∆ approximation of the streamfunction 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 max 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 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 de- fined as 1 , , In our calculations, for all Rayleigh numbers we considered that convergence was achieved when both 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 , 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 , 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.

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 0.1 χ = and a Rayleigh number Ra = 10 5 .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 3 and 10 6 (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 0% 10% χ ≤ ≤ and the Rayleigh number varies from 10 3 to 10 5 .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 3 ≤ Ra ≤ 10 5 , for the case of a water-Cu nanofluid and pure water.The value of solid volume fraction is set to 0.04 χ = . Figure 5 represents the same physical quantities but for a volume fraction value of 0.1 χ = . 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 3 , the deviation (relative to χ = 0) between the maximum values of vertical velocity is ( ) max 0.1 8.73% V ∆ = respectively.As far as the temperature distribution is con- cerned, 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 3 ) 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 max ψ 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 4 and Ra = 5 × 10 5 with a maximum heat transfer in the vicinity of X = 0.25 and X = 0.75.For example, for Ra = 5 × 10 4 , 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 3 , 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 4 , as χ is increased to 2%, an increase of 6.03% is observed.A heat transfer augmentation of above 27.63% is obtained for 8% χ = compared to 0% χ = 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.

Conclusions
In this study the heat transfer enhancement in a two dimensional enclosure filled with nanofluids is studied       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.

Figure 1 .
Figure 1.Physical model the coordinate system.

Figure 2 .
Figure 2. Grid points at the wall and at the corner: (a) wall points and (b) corner points.andthis 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 for different solide volume fractions χ of nanoparticles (Cu) is varied as 0% 10% χ ≤ ≤

Figure 6 .
Figure 6.Velocity profiles along the mid-plane for different Ra and different solid volume fractions (Water-Cu).

Figure 7 .
Figure 7. Local Nusselt number through the heated wall for different Ra and solid volume fractions (water-Cu).

Table 6 .
Comparison of average Nusselt number Nu for different Rayleigh number and various solid volume fractions.

Table 4 .
Grid independency results for water-Cu,

Table 5 .
Comparaison between the present work and other studies for Nu .