Numerical Solution of Natural Convection Problems by a Meshless Method

Natural convection is a phenomenon where fluid motion is generated by density changes due to the temperature or concentration variations in a gravity field. The computational modelling of systems with natural convection (Bejan, 2004) has become a highly popular research subject due to its pronounced influence in better understanding of nature as well as in the development of the advanced technologies. Melting of the polar ice caps, the global oceans dynamics, various weather systems, water transport, soil erosion and denudation, magma transport and manufacturing of nano-materials, improving casting processes, energetic studies, exploitation of natural resources, welding, casting and advanced solidifications are two typical contemporary example groups where natural convection plays an important role. This chapter deals with the numerical approach towards solution of this type of problems by a meshless technique. The main part of the solution procedure is focused on the general transport equation treatment and the pressure velocity coupling strategy. The transport phenomena are solved by a local meshless method and explicit time stepping. The local variant of Radial Basis Function Collocation Method (LRBFCM) has been previously developed for diffusion problems (Sarler and Vertnik, 2006), convection-diffusion solid-liquid phase change problems (Vertnik and Sarler, 2006) and subsequently successfully applied in industrial process of direct chill casting (Vertnik, et al., 2006). The fluid flow, which is generally a global problem, is treated by the proposed local iterative method. Instead of solving the pressure Poisson equation or/and pressure correction Poisson equation (Divo and Kassab, 2007) a more simplified local pressure-velocity coupling (LPVC) (Kosec and Sarler, 2008a) algorithm is proposed where the pressure-correction is predicted from the local mass continuity violation similar to the SOLA algorithm (Hong, 2004). The presented solution procedure represents a variant of already developed global approach (Sarler, et al., 2004a, Sarler, 2005). In this chapter, such a local solution procedure is tested with the standard free fluid flow benchmark test (de Vahl Davis natural convection test (de Vahl Davis, 1983)). The test is especially convenient for benchmarking purposes as there are several numerical solutions published in the literature (Divo and Kassab, 2007, Hortmann, et al., 1990, Manzari, 1999, Prax, et al., 1996, Sadat and Couturier, 2000, Sarler, 2005, Wan, et al., 2001).

In addition to the basic test, the proposed local solution procedure is tested on the tall cavity, natural convection in the porous media and melting of a pure material (Gobin -Le Quéré test) driven by a natural convection tests. Numerous analyses and comparisons with the published data are performed in order to assess the characteristics of the proposed numerical approach in details. The comprehensive verification procedure shows excellent agreements with the previously published data, based on different numerical methods.

Governing equations 2.1 De Vahl Davis test
The de Vahl Davis natural convection benchmark problem is described by three coupled partial differential equations (PDEs). The PDEs are mass, momentum and energy conservation equations where all material properties are considered to be constant. The Boussinesq approximation is used for density hypothesis and the phenomena is thus described by the following system of equations where ( ) ,, , , , , , , ,a n d ( ) 1, 0 and the initial conditions as follows ( ) www.intechopen.com ( ) ,00 . 5 where Ω and Γ stand for interior and boundary nodes indexes, respectively.
where W Ω and H Ω stand for domain width and height, respectively. The problem is characterized by three dimensionless numbers; the thermal Rayleigh number ( ) Ra T , the Prandtl number ( ) Pr and the domain aspect ratio ( ) www.intechopen.com where the ratio between the Prandtl and the Rayleigh number is known as the Grashof number The de Vahl Davis benchmark is limited to the natural convection of the air in a rectangular cavity with aspect ratio R A1 = and Pr 0.71. = In this work additional tests are done for lower Prandtl number and higher aspect ratio in order to test the method in regimes similar to those in the early stages of phase change simulations of metal like materials where the oscillatory "steady-state" develops.

Porous media natural convection
A variant of the test, where instead of the free fluid, the domain is filled with porous media, is considered in the next test. Similar to the de Vahl Davis benchmark test, the porous natural convection case is also well known in the literature (Chan, et al., 1994, Jecl, et al., 2001, Ni and Beckermann, 1991, Prasad and Kulacky, 1984, Prax, et al., 1996, Raghavan and Ozkan, 1994, Šarler, et al., 2000, Šarler, et al., 2004a, Šarler, et al., 2004b and therefore a good quantitative comparison is possible. The only difference from de Vahl Davis case is in the momentum equation, and the consecutive velocity boundary conditions. Instead of the Navier-Stokes the Darcy momentum equation is used to describe the fluid flow in the porous media where K stands for permeability. The main difference in the momentum equation is in its order. The Navier-Stokes equation is of the second order while the Darcy equation is of the first order and therefore different boundary conditions for the velocity apply. Instead of the no-slip boundary condition for velocity, the slip and impermeable velocity boundary conditions are used. This is formulated as Instead of the thermal Rayleigh and Prandtl numbers, the filtration Rayleigh number defines ( )

Phase change driven by natural convection
The benchmark test is similar to the previous cases with an additional phase change phenomenon added. The solid and the liquid thermo-physical properties are assumed to be equal. In this case the energy transport is modelled through enthalpy (h) formulation. The concept is adopted in order to formulate a one domain approach. The phase change phenomenon is incorporated within the enthalpy formulation with introduction of liquid fraction (f L ). The problem is thus defined with equations (1), (2), (4) and The phase change of the pure material occurs exactly at the melting temperature which produces discontinues in the enthalpy field due to the latent heat release. The constitutive relation (24) incorporates a smoothing interval near the phase change in order to avoid numerical instabilities.
and initial state to ( ) www.intechopen.com

( )
,00 Velocity in the solid state is forced to zero by multiplying it with the liquid fraction. This approach introduces additional smoothing in the artificial "mushy" zone. This smoothing produces an error of the same magnitude as smoothing of the enthalpy jump at the phase change temperature. The problem is schematically presented in Figure 2. Additional dimensionless number to characterize the ratio between the sensible and latent heat, the Stefan number, is introduced

Solution procedure
There exist several meshless methods such as the Element free Galerkin method, the Meshless Petrov-Galerkin method, the point interpolation method, the point assembly method, the finite point method, the smoothed particle hydrodynamics method, the reproducing kernel particle method, the Kansa method (Atluri and Shen, 2002a, Atluri and Shen, 2002b, Atluri, 2004, Chen, 2002, Gu, 2005, Kansa, 1990a, Kansa, 1990b, Liu, 2003, etc. However, this chapter is focused on one of the simplest classes of meshless methods in development today, the Radial Basis Function (Buhmann, 2000) Collocation Methods (RBFCM) (Šarler, 2007). The meshless RBFCM was used for the solution of flow in Darcy porous media for the first time in (Šarler, et al., 2004a). A substantial breakthrough in the development of the RBFCM was its local formulation, LRBFCM. Lee at al. (Lee, 2003) demonstrated that the local formulation does not substantially degrade the accuracy with respect to the global one. On the other hand, it is much less sensitive to the choice of the RBF shape and node distribution. The local RBFCM has been previously developed for diffusion problems , convection-diffusion solid-liquid phase change problems ) and subsequently successfully applied in industrial process of direct chill casting . In this chapter a completely local numerical approach is used. The LRBFCM spatial discretization, combined with local pressure-correction and explicit time discretization, enables the consideration of each node separately from other parts of computational domain. Such an approach has already been successfully applied to several thermo-fluid problems (Kosec and Šarler, 2008a, Kosec and Šarler, 2008b, Kosec and Šarler, 2008c, Kosec and Šarler, 2008d, Kosec and Šarler, 2009) and it shows several advantages like ease of implementation, straightforward parallelization, simple consideration of complex physical models and CPU effectiveness. An Euler explicit time stepping scheme is used for time discretization and the spatial discretization is performed by the local meshfree method. The general idea behind the local meshless numerical approach is the use of a local influence domain for the approximation of an arbitrary field in order to evaluate the differential operators needed to solve the partial differential equations. The principle is represented in Figure 3. Each node uses its own support domain for spatial differential operations; the domain is therefore discretized with overlapping support domains. The approximation function is introduced as where ,, a n d Basis n n N θα Ψ stand for the interpolation function, the number of basis functions, the approximation coefficients and the basis functions, respectively. The basis could be selected arbitrarily, however in this chapter only Hardy's Multiquadrics (MQs) with σ C standing for the free shape parameter of the basis function, are used. By taking into account all support domain nodes and equation (32), the approximation system is obtained. In this chapter the simplest possible case is considered, where the number of support domain nodes is exactly the same as the number of basis functions. In such a case the approximation simplifies to collocation. With the constructed collocation function an arbitrary spatial differential operator (L) can be computed In this work only five node support domains are used and therefore a basis of five MQs is used as well.
Influence domain nodes The implementation of the Dirichlet boundary condition is straightforward. In order to implement Neumann and Robin boundary conditions, however, a special case of interpolation is needed. In these boundary nodes the function directional derivative instead of the function value is known and therefore the equation in the interpolation system changes to in the Neumann boundary nodes and to in the Robin boundary nodes.

www.intechopen.com
With the defined time and spatial discretization schemes, the general transport equation under the model assumptions can be written as where 0,1 0 ,, a n d Dt S θ Δ stand for the field value at current and next time step, general diffusion coefficient, time step and for source term, respectively. To couple the mass and momentum conservation equations a special treatment is required.
The equation (38) did not take in account the mass continuity. The pressure and the velocity corrections are added 1mm where , and mv P stand for pressure velocity iteration index, velocity correction and pressure correction, respectively. By combining the momentum and mass continuity equations the pressure correction Poisson equation emerges Instead of solving the global Poisson equation problem, the pressure correction is directly related to the divergence of the intermediate velocity where stands for characteristic length. The proposed assumption enables direct solving of the pressure velocity coupling iteration and thus is very fast, since there is only one step needed in each node to evaluate the new iteration pressure and the velocity correction. With the computed pressure correction the pressure and the velocity can be corrected as 11ˆ and mm m m t PPP P ζ ζ ρ where ζ stands for relaxation parameter. The iteration is performed until the criterion

Results
The results of the benchmark tests are assessed in terms of streamfunction ( ) Ψ , cavity Nusselt number ( ) Nu and mid-plane velocity components. () www.intechopen.com The Nusselt number is computed locally on five nodded influence domains, while the streamfunction is computed on one dimensional influence domains each representing an x row, where all the nodes in the row are used as an influence domain. The streamfunction is set to zero in south west corner of the domain ( ) The de Vahl Davis test represents the first benchmark test in the series and therefore some additional assessments regarding the numerical performance as well as the computational effectiveness are done. One of the tests is focused on the global mass continuity conservation, which indicates the pressure-velocity coupling algorithm effectiveness. The global mass leakage is analysed by implementing avg avg avg where avg ,a n d D N ρ ρ Δ stand for average density, number domain nodes and density change.
The pressure-velocity coupling relaxation parameter ζ is set to the same value as the dimensionless time-step in all cases. The reference values in the Boussinesq approximation are set to the initial values.

De Vahl Davis test
The classical de Vahl Davis benchmark test is defined for the natural convection of air (Pr 0.71 = ) in the square closed cavity ( R A1 = ). The only physical free parameter of the test remains the thermal Rayleigh number. In the original paper (de Vahl Davis, 1983) de Vahl Davis tested the problem up to the Rayleigh number 6 10 , however in the latter publications, the results of more intense simulations were presented with the Rayleigh number up to 8 10 . Lage and Bejan (Lage and Bejan, 1991) showed that the laminar domain of the closed cavity natural convection problem is roughly below 9 Gr<10 . It was reported (Janssen andHenkes, 1993, Nobile, 1996) that the natural convection becomes unsteady for 8 Ra 2 10 =⋅ . This section deals with the steady state solution and therefore regarding to the published analysis, a maximum 8 Ra 10 T = case is tested. A comparison of the present numerical results with the published data is stated in Table 1 where the mid (0.5,0.5) ψψ = , avg Nu , max (0.5, ) x y v p and max (, 0 . 5 ) xy vp stand for mid-point streamfunction, average Nusselt number and maximum mid-plane velocities, respectively. The results of the present work are compared to the (de Vahl Davis, 1983) (a), (Sadat and Couturier, 2000) (b), (Wan, et al., 2001) (c) and (Šarler, 2005) (d). The specifications of the simulations are stated in Table 2. The temperature contours (yellow-red continuous plot) and the streamlines are plotted in Figure 4 with the streamline contour plot step 0.05 for  Figure 5 in order to characterize the system dynamics.
Due to the completely symmetric problem formulation ( ( ) ,00 . 5 Tt Ω == p ) the cold side and the hot side average Nusselt numbers should be the same at all times and therefore the www.intechopen.com difference between the two can be understood as a numerical error of the solution procedure. A simple relative error measure is introduced as where avg max Nu and Nu stand for average and maximum Nusselt number. The Nusselt number as a function of time is presented in Figure 5. The hot-cold side errors ( ) E a r e plotted in Figure 6 and the mid-plane velocities are presented in Figure 7.

Low Prandtl number -tall cavity test
The next test is merely a generalization of de Vahl Davis with the same governing equations, initial state and boundary conditions, where the cold side of the domain is set to the initial temperature. The tall cavity with aspect ratio R A1 / 4 = is filled with metal like (Al-4.5%Cu) low Prandtl fluid Pr 0.0137 = at the Rayleigh number 5 Ra 2.81 10 =⋅ is considered. The case is especially interesting due to its oscillating »steady-state« which is a result of a balance between the buoyancy and the shear forces. This case is also relevant for fluid flow behaviour at initial stages of melting of low Prandtl materials (metals), since it has a similar geometrical arrangement. The test case has been already computed by two different numerical methods (Založnik, et al., 2005) (spectral FEM and FVM) with good mutual agreement. Respectively, these solutions have been used for assessment of the present method as well. The temperature contours and the streamlines are plotted as continuous red to yellow fill and dotted lines with a contour step 0.2, respectively, in Figure 8. To confirm the agreement of the results, the hot side Nusselt number frequency domains are compared, where the early stages of signal development are omitted. From Figure 9 one can see that the agreement with reference results is excellent. In Figure 9 the frequency domains for different node distributions are compared, as well. The Nu freq stands for Nusselt number transformation to the frequency domain and f stands for dimensionless frequency. The presented case is highly sensitive; even the smallest changes in the case setup affect the results dramatically, for example, changing the aspect ratio for less than 1 % results in completely different flow structure. Instead of two there are three major oscillating vortices. On the other hand, the presented results are computed with the completely different numerical approach in comparison with the reference solutions (meshless spatial www.intechopen.com discretization against FVM and Spectral method, explicit against implicit time discretization scheme and LPVC against SIMPLE pressure-velocity coupling) and still the comparison shows a high level of agreement. The presented comparison infers on a high level of confidence in the present novel meshless method and the local solution approach.

Porous media test
The next test is focused on the assessments of the solution procedure behaviour when working with fluid flow in the Darcy porous media. Again, a symmetrical differentially heated rectangular cavity is considered with impermeable velocity boundary condition.  Table 3. A comparison of the results and numerical parameters Three different aspect ratios are tested A R = [0.5,1,2] for filtration Rayleigh numbers up to 10 4 . The results are compared against (a) (Šarler, et al., 2000), (b) (Ni and Beckermann, 1991), and (c) (Prax, et al., 1996) with good agreement achieved (Table 3). In addition to the previously treated cases in quoted works, results for Ra F = 10 3 and Ra F = 10 4 are newly represented in this work. The pressure-velocity coupling algorithm was tested for up to 160797 uniformly distributed nodes and it behaves convergent. The temperature and the www.intechopen.com streamfunction contours are presented in Figure 11 for tall cavity, Figure 12 for low cavity and Figure 13 for square cavity with streamline steps 2, 5, 15 and 40 for Ra F = 50, Ra F = 10 2 , Ra F = 10 3 and Ra F = 10 4 , respectively. Additional comparison of the results with reference Finite Volume Method (FVM) solution, previously used in (Šarler, et al., 2000) for mid-plane velocities, hot side Nusselt number, mid-plane and top temperature profiles is done for case with A R = 1 and filtration Rayleigh number Ra F = 100 ( Figure 14). The comparison shows good agreement with the generally accepted solution.     Table 4. The reference values in the Boussinesq approximation are set to the initial values. A comparison of the results or demonstration for all four cases is done at specific dimensionless time C t , stated in Table 4, unless stated otherwise. The streamfunction and temperature contour plots are shown in Figure 15 with streamline steps: 0.1, 0.5 and 4.0 for Cases 1,2 and 3,4, respectively. The phase change front comparison is demonstrated in Figure 17.

Conclusions
In this chapter it is shown that the proposed numerical method as well as the proposed solution procedure performs well for the natural convection problems at different flow constitutive relations and regimes. A detailed analysis of the de Vahl Davis test has been performed in order to assess the method in details. The mass leakage test and the cold-hot side Nusselt number comparison both confirmed that the method is accurate when such type of problems are considered. Furthermore, a comparison with already published data shows good agreement as well. The test involving a tall cavity, where the method is compared to the two completely different approaches, gives excellent agreement in case of the oscillatory flow regimes. The time average hot side Nusselt number development shows good quantitative comparison with the benchmark data. To complete the tests, the natural convection in a Darcy momentum regime is performed. Again, good agreement with the published data is achieved. The final test is the phase change driven by a natural convection. The present results show good agreement with other approaches in terms of interphase boundary dynamics and complicated flow structures despite the simplest LRBFCM implementation. The Nusselt number oscillations in the Case 2 were already reported by Mencinger (Mencinger, 2003). The oscillations are a result of an unstable flow regime in the low Prandtl fluid, similar as in the tall cavity natural convection case, where the periodic solutions occur. The potential instabilities can occur in the natural convection in liquid metals, due to their low Prandtl number (Založnik, et al., 2005, Založnik and. Complex flow patterns and fast transients can occur already in laminar regimes at relatively low Rayleigh numbers. With another words; lowering the Prandtl number increases the nonlinearity of the natural convection. Generally, the main sources of instabilities for all presented cases stems from nonlinearities due to the complex liquid flow pattern and enthalpy behaviour at the phase change interphase. The enthalpy jump problem is partly resolved by introducing numerical smoothing of the phase change, at least enough to get stable results, but at the price of physical model accuracy. A detailed discussion on the parameter range with appearance of the flow physics based oscillations can be found in the work (Hannoun, et al., 2003). However, the results are in good agreement with the already known solutions (Gobin and Le Quéré, 2000). There is a bit higher deviance in the Case 4, still the deviations of other authors for that case is very high and it is difficult to conclude which solution is more credible. The convection and conduction heat transfer, thermal conductivity, and phase transformations are significant issues in a design of wide range of industrial processes and devices. This book includes 18 advanced and revised contributions, and it covers mainly (1) heat convection, (2) heat conduction, and (3) heat transfer analysis. The first section introduces mixed convection studies on inclined channels, double diffusive coupling, and on lid driven trapezoidal cavity, forced natural convection through a roof, convection on non-isothermal jet oscillations, unsteady pulsed flow, and hydromagnetic flow with thermal radiation. The second section covers heat conduction in capillary porous bodies and in structures made of functionally graded materials, integral transforms for heat conduction problems, non-linear radiative-conductive heat transfer, thermal conductivity of gas diffusion layers and multi-component natural systems, thermal behavior of the ink, primer and paint, heating in biothermal systems, and RBF finite difference approach in heat conduction. The third section includes heat transfer analysis of reinforced concrete beam, modeling of heat transfer and phase transformations, boundary conditions-surface heat flux and temperature, simulation of phase change materials, and finite element methods of factorial design. The advanced idea and information described here will be fruitful for the readers to find a sustainable solution in an industrialized society.