Abstract
This paper presents a numerical study of natural convection within a wavy enclosure heated via corner heating. The considered enclosure is a square enclosure with left wavy side wall. The vertical wavy wall of the enclosure and both of the corner heaters are maintained at constant temperature, and , respectively, with while the remaining horizontal, bottom, top and side walls are insulated. A penalty element-free Galerkin approach with reduced gauss integration scheme for penalty terms is used to solve momentum and energy equations over the complex domain with wide range of parameters, namely, Rayleigh number (Ra), Prandtl number (Pr), and range of heaters in the - and -direction. Numerical results are represented in terms of isotherms, streamlines, and Nusselt number. It is observed that the rate of heat transfer depends to a great extent on the Rayleigh number, Prandtl number, length of the corner heaters and the shape of the heat transfer surface. The consistent performance of the adopted numerical procedure is verified by comparison of the results obtained through the present meshless technique with those existing in the literature.
1. Introduction
Natural convection within closed cavities with different fluids has attracted considerable attention of many researchers [1–6] due to its immense applications in engineering such as the cooling of electronic devices, refrigerators, room ventilating, heat exchangers, and solar collectors. A good amount of literature is available for convection patterns in differentially heated rectangular enclosures. The basic benchmark study related to this subject was reported by De Vahl Davis [2] and De Vahl Davis and Jones [3]. Later, many authors [4–6] studied the impact of different boundary conditions on convection phenomenon within an enclosure. In these studies, different shapes of heat transfer surfaces were not considered.
However, the analysis of natural convection heat transfer with different shapes of heat transfer surfaces (e.g., wavy surfaces or rough surfaces) is extremely important in many applications where natural convection is the only feasible mode of cooling. The shape of the surface influences the performances of the heat transfer equipment. Hence, there has been considerable interest in the intentional roughening of the surface to enhance the heat transfer. This small scale roughness may be represented by periodic functions (sine and cosines). A literature review of free convection heat transfer from vertical surfaces with surface roughness elements in continuum fluid is provided by Bhavnani and Bergles [7]. Varol and Oztop [8] performed a comparative numerical study on natural convection in wavy and flat plate solar collectors and it was observed that heat transfer rate is higher for wavy collector as compared to the flat plate collector. Dalal and Das [9] discussed laminar natural convection in an inclined complicated wavy cavity with spatially variable wall temperature. Saha [10] analyzed magnetoconvection inside a sinusoidal corrugated enclosure where bottom surface is heated with discrete isoflux using finite element method.
Besides, the linear distribution of temperature along the heated surface behavior of many engineering systems such as building construction elements with passive solar heating, solar pond, and liquid fuel storage tank may be characterized with concentrated heat source on one side or having partial heaters. Chu et al. [11] investigated the natural convection in an enclosure with a partial heater located at the left vertical wall and the cooling condition is applied on right vertical wall, both experimentally and numerically. Zhao et al. [12] analyzed double-diffusive natural convection in a porous enclosure which is partially heated from one side. It was observed that both the location and length of the heater have considerable impact on temperature and flow field. More work on this can be found in the literature [13–17].
In the present study, our aim is to investigate the effect of corner heaters on flow field and temperature distribution within an enclosure which has complex geometry due to wavy side wall. An advanced meshfree numerical technique, element free Galerkin method (EFGM), is used to solve the nonlinear mathematical model of the problem. Although finite element method and finite difference methods are quite efficient and general techniques, they behave poorly for the simulation of problems where complex geometry is involved, due to tedious and expensive meshing, remeshing procedure. To overcome these difficulties, a number of meshfree methods such as element free Galerkin method, meshless local Petrov-Galerkin method, and smooth particle hydrodynamics method have been proposed.
Nikfar and Mahmoodi [18] applied meshless local Petrov-Galerkin method for the analysis of natural convection of nanofluids inside a cavity having wavy side walls. The element free Galerkin method has also been successfully used to solve various problems in different areas such as heat transfer [19] and fracture mechanics [20]. Recently, Singh and Bhargava [21] have applied EFGM for the simulation of an unsteady micropolar squeeze film flow.
Hence, in the present study EFGM has been used as a tool to solve the coupled, nonlinear partial differential equations governing the flow inside a wavy enclosure. Obtained results are compared with benchmark results and excellent agreement has been observed between them.
2. Mathematical Model and Definition of the Problem
A schematic view of the two-dimensional enclosure with one wavy side wall is depicted in Figure 1. The left wavy wall of the enclosure is described by The wavy wall of the enclosure is maintained at a constant temperature and its temperature is lower than that of the corner heaters, maintained at constant temperature . The remaining top, bottom, and side walls are insulated. All the physical properties of the fluid are assumed to be constant except the density which is determined according to Boussinesq approximation, primarily used in the field of buoyancy-driven flow. It is assumed that viscous dissipation and radiation effects are negligible and gravity acts in the negative -direction. For a two-dimensional flow of an incompressible Newtonian fluid in steady state regime, the governing continuity, momentum, and energy equations can be obtained with the following dimensionless variables [4, 5]: The governing equations are written as The physical boundary conditions as shown in geometric model (Figure 1) can be defined as follows.
On the wavy wall:
On the adiabatic walls:
On the corner heaters:
2.1. Nusselt Number Evaluation
The heat transfer rate evaluation is the most important for many engineering applications. Heat transfer rates on the surface of isothermal horizontal, vertical heaters, and cold wavy surface are computed by local Nusselt number. Mean Nusselt number is calculated by averaging local nusselt number along the cold wavy side of the enclosure.
The local nusselt number along the cold wavy surface is given as The heat transfer rate in terms of local Nusselt number along the horizontal (HS1) or vertical (HS2) heaters is expressed as The mean nusselt number is computed along the cold wavy surface and is given as follows,
3. Solution Methodology and Postprocessing
The momentum and energy conservation equations are solved using a meshfree numerical technique known as element free Galerkin method. The continuity equation is used as a constraint due to mass conservation and pressure distribution [4, 5]. The penalty approach [24, 25] is employed to impose this constraint. In order to solve momentum equations (4)-(5), pressure is eliminated by a penalty parameter and the incompressibility condition given in (3) results For large values of penalty parameter , the continuity equation (3) is automatically satisfied. Using (13), the momentum conservation equations (4)-(5) result in Now, the system of simultaneous partial differential equations (6), (14), and (15) is solved numerically using element free Galerkin method.
3.1. Meshfree Method
Numerical simulation of fluid flow problems within complex geometries is a computational challenge. Although finite element method (FEM) is one of the advanced and general numerical technique, meshing, remeshing of complex geometries is quite difficult and expensive in FEM. To overcome these problems, a meshfree numerical technique element free Galerkin method is used for the numerical simulation. The essence of meshfree techniques is that, instead of a predefined mesh, they use only a set of nodes scattered in the whole problem domain without any fixed connectivity which is quite easy process as compared to the tedious meshing procedure in FEM.
3.2. Element-Free Galerkin Method
The element free Galerkin method (EFGM) requires moving least square (MLS) interpolation functions to approximate an unknown function. The MLS approximant requires only a set of nodes for its construction and is made up of these components: a compact support weight function associated with each node, a polynomial basis function, and a set of coefficients that depends on node position. The weight function is nonzero over a small neighborhood of a particular quadrature point or evaluation point, in which small neighborhood area is called support domain of the quadrature point. A view for selecting support domain of a quadrature point is shown in Figure 2.
Using MLS approximation, the unknown field variable is approximated over the two-dimensional domain as (details can be seen in [26]) where , is the number of terms in the basis, is the monomial basis function, is the nonconstant coefficient functions, and .
The coefficients are determined by minimizing the functional given by where is a weight function which is nonzero over a small domain, called support domain, and is the number of nodes that are included in the support domain of a point . The minimization of with respect to leads to the following set of equations: where and are given as Substituting (19) in (16), the MLS approximant is obtained as where the shape function is defined by
3.2.1. Weight Function Description
The choice of weight function affects the resulting approximation in EFGM and other meshless methods. In EFGM, the continuity of MLS approximants is governed by the continuity of weight function. Singh et al. [25] have studied these weight functions and reported that cubic spline weight function gives more accurate results as compared to others. Therefore, in the present work, cubicspline weight function has been used.
3.2.2. Cubicspline Weight Function
Consider the following: where .
The shape of the support domain can be circular or rectangular or both, but rectangular support domain is more general, and therefore is used in the present simulation. In rectangular support domain, the weight function at any given point can be calculated as where , , , and , , are the sizes of support domain in the - and -direction. is a scaling parameter, and , are the distances to the nearest neighbors in the - and -direction, respectively. The size of the support domain at a particular node is only controlled by scaling parameter since the distance between nearest neighbors for an evaluation point (or quadrature point) remains unchanged for a given nodal data distribution. The minimum value of should be greater than 1 so that , and the maximum value of should be such that it preserves the local character of MLS approximation. It has been shown by Singh [27] that is the optimum range of scaling parameter for heat transfer problems. In the present simulation has been fixed as 1.2.
3.3. Variational Formulation
The weighted integral form of (14), (15), and (6) over the entire problem domain can be written as where , , are arbitrary test functions and may be viewed as the variation in , , , respectively.
3.4. Element Free Galerkin Model and Imposition of Boundary Conditions
The element free Galerkin model of (24) is obtained by substituting MLS approximation for the unknown field variables , , using (20)-(21) and can be obtained as Since MLS shape functions do not satisfy the kronecker delta property, so we cannot directly impose the essential boundary conditions. To remove this problem, different numerical techniques have been proposed to enforce the essential boundary conditions in EFG method such as Lagrange multiplier method and the penalty method. In the present simulation penalty method is applied.
3.5. Penalty Method
Using penalty method to enforce the essential boundary conditions, the variational form is written as where , , are to be replaced by the MLS shape functions for obtaining stiffness matrix expressions. is the penalty parameter imposed to enforce the essential boundary conditions. The values of , , along the boundaries are prescribed in Figure 1 and (7)–(9).
Using the EFGM model given by (25), into (26), the system of equations can be defined as follows: where each of these , , are given as follows: where and denotes the total number of nodes in the whole problem domain.
3.6. Background Integration
The whole domain is discretized with , that is, 1681 nodes in such a manner that nodes are denser near the cold wavy surface as compared to the hot flat surface. For numerical integration purpose, a background mesh of size is utilized and four-point Gauss quadrature formula has been used to evaluate the viscous terms in , , , , whereas for the evaluation of penalty terms in , , , the reduced Gauss integration scheme that is one-point Gauss quadrature formula has been used. Large values of penalty parameter , can satisfy the constraint condition or continuity condition absolutely, but due to high penalty parameter , the contribution from viscous terms in , , , (given in (28)) becomes negligible as compare to the penalty terms. Hence, it implies that, for infinitely large values of penalty parameter , the governing equations are left only with constraint condition or continuity condition and the contribution from momentum and energy conservation equations is completely lost. In order to avoid such a situation, reduced Gauss integration scheme [24] has been applied for the evaluation of penalty terms in , , , .
At each node, three functions , , are to be evaluated; hence after assembly, we obtain a nonlinear system of equations of order , as given in (27). Owing to the nonlinearity of the system, an iterative scheme has been used to solve it with an initial guess. The system of equations is linearized by incorporating known functions , , as given in (29), which is solved using Gauss elimination method. This gives a new set of values of unknowns , , and the process continues till the required accuracy (0.0005) is achieved.
3.7. Validation of the Results
We have validated our results for a canonical problem available in the literature reported by De Vahl Davis [2], Barakos et al. [22], Fusegi et al. [23], and Wan et al. [28] in Figures 3 and 4. The spatial structure of local nusselt number and mean nusselt number obtained with present EFGM technique is compared with benchmark results in Figure 4 and Table 1, respectively, and a good agreement has been observed between them.
![(b) Results reported by Barakos et al. [22] for](https://static.hindawi.com/articles/tswj/volume-2014/842401/figures/842401.fig.003b.jpg)
![(d) Results reported by Barakos et al. [22] for](https://static.hindawi.com/articles/tswj/volume-2014/842401/figures/842401.fig.003d.jpg)
![(f) Results reported by Barakos et al. [22] for](https://static.hindawi.com/articles/tswj/volume-2014/842401/figures/842401.fig.003f.jpg)
4. Results and Discussions
In the present work, numerical calculations are carried out over a wide range of parameters, Rayleigh number, , Prandtl number, , and dimensionless length of the heaters in the - and -direction to investigate their impact on thermal and flow fields. The flow and temperature fields are graphically presented in terms of stream lines and isotherm contours, respectively. Heat transfer characteristics are examined in terms of mean and local nusselt numbers at the cold wavy surface and corner heaters.
4.1. Effect of Rayleigh Number
The impact of Rayleigh number on flow and temperature field is represented in terms of streamlines and isotherms, respectively, in Figure 5. From the contours of the isotherms, it is observed that the dominant heat transfer mechanism changes from conduction to convection with increase in Ra, due to increased buoyancy. For low Ra-values the heat is transferred by conduction between hot wall and the cold wavy wall because these almost vertical isotherms lines are more concentrated in the vicinity of hot corner. For high Ra-values, the heat transfer mechanism changes from conduction to convection and the isotherm lines depart from the hot corner towards the cold wavy surface. For higher values of Rayleigh number, convection strength is increased. Comparing the results shown in Figure 3 for higher Ra-values, with the results shown in Figure 5, it could be noticed that isotherms in Figure 3 at the centre of the cavity are becoming almost horizontal showing dominant convection heat transfer mechanism which is not so in isotherm contours of Figure 5, because heating is being done on a part of the vertical and horizontal side of cavity. Therefore, for controlled heat transfer partial heaters could be used in many physical applications.
| (a) Maximum value of stream function |
| (b) Maximum value of stream function |
| (c) Maximum value of stream function |
From the stream function plots in Figure 5 for , it can be seen that fluid is unicellular and rotates in counterclockwise direction. For low Ra-values, a single circulating eddy appears as the dominant characteristic of the flow. The magnitude of the stream functions increases significantly with increase in from . Hence, higher flow strength is observed for higher values of since increases the buoyancy effect. Due to asymmetric heating, higher values of the velocities are observed in the region of corner heaters as compared to the velocities obtained in the upper half (insulated part) of the enclosure and it gives more intended streamline towards the heated corner.
4.2. Effect of Prandtl Number
Figure 6 depicts the impact of Prandtl number on isotherms and streamlines (for ). For the lowest Prandtl number , the stream function values are very low and isotherm lines are more concentrated towards the hot corner. However, for higher Prandtl number, isotherms are denser near the cold wavy side wall and higher values of stream functions are observed. In comparison with, considering the contour plots of with , we observe that they look qualitatively similar which shows that the nature of streamlines and isotherms does not change considerably with higher Prandtl number. However, the maximum stream function values are obtained with .
| (a) Maximum value of stream function |
| (b) Maximum value of stream function |
| (c) Maximum value of stream function |
4.3. Effect of Corner Heating
The impact of heater length on flow and temperature field is illustrated in Figure 7 for , . In Figures 7(a)–7(c), the horizontal length of the heater, that is, , is kept fixed at 0.6 and, results are obtained with different values of vertical heater length, that is, . From these contours, higher temperature gradient near the heater is observed with increase in heater length in vertical direction due to stronger convection and therefore heat transfer rate increases with increase in heater length. In Figure 7(d), higher heater length is chosen in horizontal direction ( and ). It shows steep temperature gradient near the heater due to stronger convection and shape of the isotherms shows plume-like distribution. Comparing Figure 7(c) ( and ) with (d) ( and ), higher values of stream function are obtained for and . It indicates that given the same total heater length, the flow strength is increased with an increase of the length in horizontal direction. At the centre of the enclosure, the highest values of the stream functions are obtained.
| (a) Maximum value of stream function |
| (b) Maximum value of stream function |
| (c) Maximum value of stream function |
| (d) Maximum value of stream function |
4.4. Effect of Wavy Surface Amplitude and Number of Undulations
In Figure 8 contour plots of stream lines and isotherms are shown for , , with different wavy amplitudes (). In Figure 9, contours of stream lines and isotherms are drawn for , and wavy surface amplitude with different numbers of undulation (). From Figure 8, it is observed that with increase in wavy surface amplitude, the inward protrusion of the isotherm lines increases. This narrows down the straight vertical passage for the fluid to move under the influence of buoyancy force and the isotherms show the steep temperature gradient near the wavy surface. Due to which higher heat transfer rate along the cold wavy surface is observed with increased wavy surface amplitude. Therefore, rough surfaces can be used for higher heat transfer rate while designing physical devices. It appears from Figures 9(a)–9(c) that varying the number of undulations does not change the global flow and isotherms pattern except in the vicinity of the vertical wavy wall. In the vicinity of the vertical wavy wall, the isotherm lines and stream lines adopt the profile of wavy wall. Increment in the number of undulations along the wavy surface has a negative impact on the heat transfer rate. Heat transfer rate is decreased with the increase of number of undulations but not significantly.
| (a) Maximum value of stream function |
| (b) Maximum value of stream function |
| (c) Maximum value of stream function |
| (a) Maximum value of stream function |
| (b) Maximum value of stream function |
| (c) Maximum value of stream function |
4.5. Heat Transfer Characteristics
The implications of various parameters, Rayleigh number, Prandtl number, heater length, wavy surface amplitude and number of undulations on energy transport process, are reported in terms of mean and local nusselt number along the cold wavy wall and corner heaters. The average nusselt number obtained along the cold wavy surface is given in Table 2 for all the parameters and spatial variation of the nusselt number for different parameters is depicted in Figures 10–14.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Due to wavy geometry, local nusselt number profiles along the cold wavy surface are also wavy in nature. Figure 10 indicates that heat transfer rate increases rapidly for higher Rayleigh number due to strong convection. The smallest value of the mean nusselt number is obtained with . Figure 11 shows the variation of heat transfer rate with prandtl number. We observe from Figure 11 that prandtl number is more significant on heat transfer rate for and that the difference between heat transfer rate obtained with and is much bigger as compared to the heat transfer rate calculated with and . Heater length both in vertical and horizontal directions has considerable impact on heat transfer rate, as shown in Figures 12 and 13. At the intersection point of the heaters, the highest value of nusselt number is obtained. Higher heat transfer rate is observed with increased heater length in vertical or horizontal direction. Figure 14 depicts the impact of wavy surface amplitude and number of undulations on heat transfer rate. Both average nusselt number along the cold way wall and local nusselt number are higher for higher wavy amplitude. It shows that heat transfer rate is higher for wavy surface as compared to the smooth surface while the impact of the number of undulations is not so significant on average nusselt number. A slight decrease in the heat transfer rate is observed with more undulations along the wavy surface.
5. Conclusions
The problem of natural convection heat transfer and fluid flow within a wavy enclosure with corner heating effect has been studied. The meshfree element free Galerkin model is demonstrated as an alternative approach to eliminate the well known drawbacks (meshing and remeshing) of grid based methods such as FDM and FEM. Obtained results are verified with available benchmark results and a good agreement has been observed. It is observed that Rayleigh number has significant impact on heat transfer rate and heat transfer rate is an increasing function of Rayleigh number. Another significant parameter which influences the heat transfer is the heater’s length in vertical and horizontal directions. However, given the same total heater length, the flow strength is increased with an increase of the length in horizontal direction. The impact of increasing prandtl number is remarkable on flow, temperature, and heat transfer for . Heat transfer rate increases with increase of prandtl number. The shape of the heat transfer surface also plays an important role. Heat transfer rate is higher for wavy surface as compared to the smooth surface.
Nomenclature
| : | gravitational acceleration m/s2 |
| : | normal direction on the plate |
| Ra: | Rayleigh number |
| Pr: | Prandtl Number |
| Nu: | mean nusselt number |
| : | wavy surface amplitude |
| : | number of undulations |
| : | temperature, K |
| : | dimensionless heater length in the -direction |
| : | dimensionless heater length in the -direction |
| : | velocities in the - and -direction, m/s |
| : | nondimensional velocities in the - and -direction |
| : | nondimensional coordinates |
| : | side of the cavity, m |
| HS1: | heated surface 1 (horizontal) |
| HS2: | heated surface 2 (vertical) |
| CS: | cold surface. |
| : | kinematic viscosity, m2 s−1 |
| : | thermal diffusivity, m2 s−1 |
| : | nondimensional temperature |
| : | thermal expansion coefficient, K−1 |
| : | nondimensional stream function. |
| : | cold |
| : | hot. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
One of the authors (Sonam Singh) would like to thank the Council of Scientific and Industrial Research (CSIR) and Government of India for their financial support through the award of a research grant.