Mechanical Engineering Faculty Publications Mechanical Engineering 2005 Application of Meshless Methods for hermal Analysis Darrell Pepper University of Nevada, Las Vegas, darrell.pepper@unlv.edu Bozidar Sarler Nova Gorica Polytechnic Follow this and additional works at: htp://digitalscholarship.unlv.edu/me_fac_articles Part of the Heat Transfer, Combustion Commons, and the Numerical Analysis and Computation Commons Citation Information Pepper, D., Sarler, B. (2005). Application of Meshless Methods for hermal Analysis. Journal of Mechanical Engineering, 51(7), 476-483. htp://digitalscholarship.unlv.edu/me_fac_articles/251 his Article is brought to you for free and open access by the Mechanical Engineering at Digital Scholarship@UNLV. It has been accepted for inclusion in Mechanical Engineering Faculty Publications by an authorized administrator of Digital Scholarship@UNLV. For more information, please contact digitalscholarship@unlv.edu. Strojniãki vestnik - Journal of Mechanical Engineering 5 1(2005)7-8, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005)7-8, 476-483 UDK-UDC 536.2 UDK - UDC 536.2 Izvirni znanstveni lanek - Original scientific paper (1.01) Izvirni znanstveni þlanek - Original scientific paper (1.01) Application of Meshless Methods for Thermal Analysis Darrell W. Pepper 1 and Božidar Šarler 2 1 Nevada Center for Advanced Computational Methods, University of Nevada Las Vegas, Las Vegas, NV, 89154-402 7, USA, dwpepper@nscee.edu 2 Laboratory for Multiphase Processes, Nova Gorica Polytechnic, Nova Gorica, Slovenia Abstract Many numerical and analytical schemes exist for solving heat transfer problems. The meshless method is a particularly attractive method that is receiving attention in the engineering and scientific modeling communities. The meshless method is simple, accurate, and requires no polygonalisation. In this study, we focus on the application of meshless methods using radial basis functions (RBFs) – which are simple to implement – for thermal problems. Radial basis functions are the natural generalization of univariate polynomial splines to a multivariate setting that work for arbitrary geometry with high dimensions. RBF functions depend only on the distance from some center point. Using distance functions, RBFs can be easily implemented to model heat transfer in arbitrary dimension or symmetry. Introduction For decades, finite difference, finite volume, and finite element methods (FDM/FVM/FEM) have been the dominant numerical schemes employed in most scientific computation. These methods have been used to solve numerous thermal related problems covering a wide range of applications. A common difficulty in these classic numerical methods is the considerable amount of time and effort required to discretize and index domain elements, i.e., creating a mesh. This is often the most time consuming part of the solution process and is far from being fully automated, particularly in 3D. One method for alleviating this difficulty has been to utilize the boundary element method (BEM). The major advantage of the BEM is that only boundary discretization is required rather than domain, thereby reducing the problem by one order. However, the discretization of surfaces in 3-D can still be a complex process even for simple shapes. In addition, these traditional methods are often slowly convergent, frequently requiring the solution of 100’s of thousands of equations in order to get acceptable accuracy. In recent years, a novel numerical technique called “meshless methods” (or “mesh-free methods”) has been undergoing strong development and has attracted considerable attention from both science and engineering communities. Currently, meshless 476 476 methods are now being developed in many research institutions all over the world. Various methods belonging to this family include: Diffuse Element Methods, Smooth Particle Hydrodynamics Methods, Element-Free Galerkin Methods, Partition of Unity Methods, h-p Cloud Methods, Moving Least Squares Methods, Local Petrov-Galerkin Methods, Reproducing Kernel Particle Methods, and Radial Basis Functions. A common feature of meshless methods is that neither domain nor surface polygonisation is required during the solution process. These methods are designed to handle problems with large deformation, moving boundaries, and complicated geometry. Recently, advances in the development and application of meshless techniques show they can be strong competitors to the more classical FDM/FVM/FEM approaches [1,2], and may likely become a dominant numerical method for solving st science and engineering problems in the 21 century. A recent book by Liu [3] discusses meshfree methods, implementation, algorithms, and coding issues for stress-strain problems. Liu [3] also includes Mfree2D, an adaptive stress analysis software package available for free from the web. Atluri and Shen [4] produced a research monograph that describes the meshless method in detail, including much in-depth mathematical basis. They also present comparison results with other schemes. Strojniãki vestnik - Journal of Mechanical Engineering 5 1(2005) 7-477, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483 Nomenclature c shape parameter C transport variable D diffusion coefficient f interior functional g boundary functional h convective film coefficient k thermal conductivity N number of nodes q heat flux Q heat source/sink r radial distance t time T temperature V velocity vector x horizontal distance There exists various types of meshless methods and each method has its advantages and disadvantages. Intensive research conducted in many major research institutions all over the world are now working to improve the performance of these approaches. In this study, we focus on the introduction of the basic concept of meshless methods using radial basis functions (RBFs) – which are simple to implement. Currently, there are two major approaches in this direction: (i) a domain-type meshless method that was developed by Kansa [5] in 1990; (ii) a boundarytype meshless method that has evolved from the BEM [6]. Radial basis functions are the natural generalization of univariate polynomial splines to a multivariate setting. The main advantage of this type of approximation is that it works for arbitrary geometry with high dimensions and it does not require any mesh. A RBF is a function whose value depends only on the distance from some center point. Using distance functions, RBFs can be easily implemented to reconstruct a plane or surface using scattered data in 2-D, 3-D or higher dimensional spaces. Similar to FEM techniques, meshless methods produce banded system matrices that can be handled in similar fashion. Both sets of methods can utilize either direct methods based on Gauss elimination or matrix decomposition methods or iterative methods, e.g., Gauss-Seidel or SOR techniques. When dealing with nonlinear problems, additional iterative loops are needed. Meshless methods generally require more CPU time since the creation of shape functions are more time-consuming and are performed during the computation. However, less time is spent in setting up computational nodes. Results using meshless methods are typically more accurate than y Į İ ı lateral distance thermal diffusivity (k/ȡcp ) emissivity Stefan-Boltzmann constant ij trial function; field variable ș dummy variable ȥ approximation function; streamfunction Ȧ vorticity Subscripts i,j nodal values I internal number of node points o initial value; free parameter Superscripts n, n+1 known, unknown values Ö approximate solution conventional numerical methods based on mesh discretizations -thus the ratio of accuracy to CPU is likely to be greater for meshless methods. Meshless methods hold promising alternative approaches for problems involving fluid flow and heat transfer analyses. The most attractive feature is the lack of a mesh that is required in the more conventional numerical approaches. This becomes particularly interesting in that one can begin to conduct adaptive analyses for CFD problems. The Meshless Method using RBFs In 1990, Kansa [5] extended the idea of interpolation scheme using RBFs to solving various types of engineering problems. The method is simple and direct and is becoming very popular in the engineering community. To illustrate the application of the meshless method using Kansa’s method, we first consider the elliptic problems. For simplicity, we consider the 2-D Poisson problem with Dirichlet boundary condition ’2 T f (x, y), (x, y)  :, T g(x, y), (x, y)  *. (1) Notice that the solution of Eq. (1) is in fact nothing but a surface. Techniques in surface interpolation can be applied to solve Eq. (1). To approximate T, Kansa [5] assumed the approximate solution could be obtained using a linear combination of RBFs Ö y) T(x, N ¦ T I(r ) j j (2) j 1 Application of meshless methods for thermal analysis 477 477 Strojniãki vestnik - Journal of Mechanical Engineering 5 1(2005) 478-8, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483 N where {T 1 , T2 , …,T N } are the unknown coefficients to be determined, ij(r j) is some form of RBF (trial function), and r is defined as (x  x j ) 2  (y  y j ) 2 rj (3) Since multiquadrics (MQ) are infinitely smooth functions, they are often chosen as the trial function for ij, i.e., I(rj ) rj2  c 2 (x  x j ) 2  (y  y j ) 2  c 2 (4) where c is a shape parameter provided by the user. The optimal value of c is still a subject of research. Other functions such as polyharmonic splines can also be chosen as the trial function. By direct differentiation of Eq. (6), the first and second derivatives of ij with respect to x and y can be expressed as wI x  xj wx rj  c 2 2 (y  y j )  c w I 2 2 wx , 2 rj  c 2 wI y  yj wy rj  c 2 , wy 2 2 rj  c 2 (5) 2 2 Substituting Eq. (3) into Eq. (1) and using collocation, one obtains § ¦T ¨ © N (x  x )  ( y  y )  2c 2 i 2 j i j j 1 2 2 (x  x )  ( y  y )  c i j 2 j i 2 3 / 2 j i n1 T n  D’ T 2 't 2 2 An implicit time marching scheme can be used and Eq. (7) becomes T 2 (x  x j )  c w I 2 Figure 1. Interior points and boundary points using Kansa’s method. · ¸ ¹ f (x , y ), i 1, 2, ! , N i n 1 f (x, y, T , ’T ) n n (8) where ǻW denotes the time step, superscript n+1 is the unknown (or next time step) value to be solved, and superscript n is the current known value. The approximate solution can be expressed as Ö , y , t n1 ) T(x i i N ¦T n 1 j Ij (x i , y i ) (9) j 1 Substituting Eq. (9) into Eq. (8), one obtains I N ¦T j (x i  x j )  (yi  y j )  c 2 2 2 g(x i , yi ), N ¦T j 1 i j 1 N I  1, N I  2, !, N wt 478 478  D’ 2 T f (x, y, T, ’ T) § Ij · 1 n 2 T (xi , yi )  ¨  D’ Ij ¸ (xi , yi ) 't © 't ¹ f (x i , y i , t , T (x i , y i ), ’T (x i , y i )), n (6) where NI denotes the total number of interior points and N I +1,…, N are the boundary points. Figure 1 shows two sets of interpolation points: interior and boundary points. Note that Eq. (6) is a linear system of N X N equations and can be solved by direct Gaussian elimination. Once the unknown coefficients {T1 , T2, …, TN } are found, the solution of T can be approximated at any point in the domain. For time dependent problems, we consider the following heat equation as an example: wT n 1 j N ¦T n 1 j n I (xi , yi ) n g(xi , yi , t n 1 ), i i 1, 2, !, N I NI  1,! , N j 1 (10) which produces an N X N linear system of equations n 1 for the unknown Tj . Note that the right hand side of the first equation in Eq. (10) can be updated b efore the next time step, i. e., (7) Pepper D.W. - Šarler B. Pepper D. W. – âarler B. Strojniãki vestnik - Journal of Mechanical Engineering 5 1(2005) 7-8, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483 N n T (x i , yi ) ¦ T I (x , y ), n j j i i j 1 N n Tx (x i , yi ) ¦T n j 1 N n Ty (x i , yi ) wI j wx j ¦T n j j 1 wI j wy (11) (xi , yi ), (xi , yi ) Heat Transfer Applications To illustrate the use of meshless methods, let us begin with a simple heat transfer problem. The governing equation for temperature transport can be written as wT  V<’T wt D’2T  Q (12) q  k’T  h(T  Tf )  HV(T4  Tf4 ) T(x, 0) 0 To (13) (14) where V is the vector velocity, x is vector space, T(x,t) is temperature, T’ is ambient temperature, To is initial temperature, D is thermal diffusivity (ț/ȡc p), İ is emissivity, ı is the Stefan-Boltzmann constant, h is the convective film coefficient, q is heat flux, and Q is heat source/sink. Velocities are assumed to be known and typically obtained from solution of the equations of motion (a separate program is generally used for fluid flow [7]). In this first example, a two-dimensional plate is subjected to prescribed temperatures applied along each boundary [8], as shown in Fig. 2. The temperature at the mid-point (1,0.5) is used to compare the numerical solutions with the analytical solution. The analytical solution is given as T(x, y) { T  T1 2 T2  T1 S ¦ 8 n 1 1 n 1 n 1 sin nS x sinh(n Sy / L) L sinh(n SW / L) which yields ș(1,0.5) = 0.445, or T(1,0.5) = 94.5 oC. Table 1 lists the final temperatures at the mid-point using a finite element method, a boundary element method, and a meshless method, compared with the exact solution. Table 1. Comparison of results for example 1 Method Exact FEM BEM Meshless o mid-pt ( C) 94.512 94.605 94.471 94.514 Elements 0 256 64 0 Nodes 0 289 65 325 As a second example, a two-dimensional domain is prescribed with Dirichlet and Neumann boundary conditions applied along the boundaries, as shown in Fig. 3(a,b,c). This problem, described in Huang and Usmani [2], was used to assess an h-adaptive FEM o technique for accuracy. A fixed temperature of 100 C is set along side AB; a surface convection of 0oC acts o along edge BC and DC with h = 750 W/m C and k = o 52 W/m C. The temperature at point E is used for comparative purposes. The severe discontinuity in boundary conditions at point B creates a steep temperature gradient between points B and E. Figures 3(b,c) show the initial and final FEM meshes after two adaptations using bilinear triangles. The analytical solution for the temperature at point B is T o = 18.2535 C. Table 2 lists the results for the three methods compared with the exact solution. The initial 3-noded triangular mesh began with 25 elements and 19 nodes. Table 2. Comparison of results for example 2 Method Exact FEM BEM Meshless Pt E (oC) 18.2535 18.1141 18.2335 18.253 1 Elements 0 256 32 0 Nodes 0 155 32 83 A simple irregular domain is used for the third example and results compared with the three methods. Results from a fine mesh FEM technique (without adaptation) are used as a reference benchmark [7]. The discretized domain and accompanying boundary conditions set along each Application of meshless methods for thermal analysis Application of meshless methods for thermal analysis 479 479 Strojniãki vestnik - Journal of Mechanical Engineering 5 1(2005) 7-480, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483 surface are shown in Fig. 4. The FEM results are displayed as contour intervals. Figure 5(a,b) shows meshless results (using FEM fine mesh nodes for contouring) versus FEM solutions using adapted quadrilateral elements. Heat conduction occurs as a result of constant temperatures set on the top and bottom surfaces, adiabatic faces in the upper right cutout and lower cutout portions, and convective heating along the right and left vertical walls. Adaptive meshing occurs in the corners as a result of steep temperature gradients; this is not evident when using meshless methods. (c) Figure 3. Problem (a) geometry - boundary conditions, (b) initial FEM mesh, and (c) final FEM adapted mesh (from [2]) Figure 4: Problem specification for heat transfer in a user-defined domain. The FEM, BEM, and meshless mid-point values at (0.5,0.5) are listed in Table 3. Table 3. Comparison of results for example 3 Method FEM BEM Meshless mid-pt (oC) 75.899 75.885 75.893 Elements 138 36 0 Nodes 178 37 96 All three techniques provide accurate results for the three example cases. The meshless method was clearly the fastest, simplest, and least storage demanding method to employ. Advances being made in meshless methods will eventually enable the scheme to compete with the FEM and BEM on a much broader range of problems [3,4]. Dr. Y. C. 1 Hon is a leading expert in the application of Kansa’s method. Much work in engineering 480 480 Pepper D.W. - Šarler B. Pepper D. W. – âarler B. Strojniãkivestnik vestnik--Journal JournalofofMechanical MechanicalEngineering Engineering51(2005) 5 1(2005)7-8, 7-8,476-483 476-483 Strojniški modeling using Kansa’s method has been done by his research group. Natural Convection Test Case Natural convection within a 2-D rectangular enclosure is a well-known problem commonly used to test the ability of a numerical algorithm to solve for both fluid flow and heat transfer. The equations are strongly coupled through the buoyancy term in the momentum equations and the temperature. There are various ways to nondimensionalize the equations, and numerous references can be found in the literature and on the web regarding these various forms. The solution to the problem generally splits between solving either the primitive equations for velocity or the vorticity equation, coupled with the transport equation for temperature. The issue in this early development of the meshless approach is not to dwell on various schemes to deal with pressure (e.g., projection methods or the SIMPLE scheme both of which are well known). Hence, most researchers that have developed meshless approaches use the streamfunction-vorticity and temperature equations [3]. These equations are the well-known set generally formulated as follows: Figure 5: FEM solutions (a) meshless (on FEM fine mesh) and (b) adapted mesh. Department of Mechanical Engineering, Hong Kong University, Hong Kong, China 1 wZ  V<’Z Pr ’2 Z  Pr <Ra < ’T wt (20) wT  V<’T ’ 2T wt (21) ’2 \ Z (22) where Ȧ is vorticity and ȥ is streamfunction, with the conventional definitions for velocity in terms of the streamfunction gradients. Pr is the Prandtl number (Ȟ/Į) and Ra is the Rayleigh number. Figure 6 shows the physical and computational domain with accompanying boundary conditions. Two types of nodal configurations are shown in Fig. 7 (a,b) utilizing 256 nodes. Results are in excellent agreement with well-known results in the literature for 103” Ra ” 105 [3]. Figure 8 (a,b) shows streamlines and isotherms for the differentially heated enclosure for Ra = 10 5. Convergence rates showing the difference in rates between a conventional FDM and applications of two meshless techniques is discussed in Liu [3]. The two meshless methods converged more rapidly than the finite difference scheme. Application of meshless methods for thermal analysis Application of meshless methods for thermal analysis 481 481 Strojniãki vestnik - Journal of Mechanical Engineering 51(2005) 482-8, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483 w (UC)  ’ <(U VC) ’ <( D’C)  S wt (23) where ȡ, ij, V, t, D, and S denote density, transport variable, velocity, time, diffusion matrix, and source. The transport variable C consisted of enthalpy C(h(ij = T)), velocity C(ij = u,v), and pressure C(ij = p), with a pressure correction Poisson equation used to resolve the pressure. The nonlinear equations solved with the meshless technique were of the form ’ 2I Figure 6. Boundary conditions for natural convection within a rectangular enclosure (from [3]). Figure 7. Nodal configurations for a) uniform distribution and b) arbitrary distribution for 256 nodes (after [3]). T  ’ <Ĭ (24) T ªw º (U 'C  S) » / D ¬ wt ¼ (25) Ĭ >U 'V C  D ' ’I@ / D (26) where ȡƍ denotes density, C(ij) the transport variable, t is time, V is velocity, and D is the diffusion matrix with D' being the nonlinear anisotropic part. The variable C(ij) denotes the relation between the transported and the diffused variable. The solution requires the use of an iterative technique. The final form of the transformed Poisson equation is ’ 2I T  TI ( I  I)  ’< Ĭ  ’ <Ĭ I (I  I ) (27) where the bar denotes values from the previous iteration. Time discretization utilizes the relation T| ¬ ª UC(I )  U C(I o ) º  S» / D 't ¬ ¼ (28) with the unknown field ij approximated by the N global approximation functions ȥn (p) and their coefficients Ȣn , Figure 8. Natural convection results showing a) i.e., streamlines and b) isotherms for Ra = 10 5 using the I(p) | \ n (p)9 n , n 1, 2,! , N * (29) MLPG method (after [3]). Sarler et al [9] simulated natural convection within a rectangular enclosure using the RBF approach of Kansa [5]. Solving a nonlinear Poisson re-formulation of the general transport equation representing mass, energy, and momentum, the problem was solved by dividing the physical domain into two parts consisting of an internal array of nodes and a set of boundary nodes for the Dirichlet and Neumann conditions. The governing equation for the transport variable is of the form (with C(ij)) 482 482 The global radial basis function approximation was based on multiquadrics with the free parameter ro: \n rn  ro 2 2 1/ 2 (30) The coefficients were calculated from the N collocation equations of which Nī were equally distributed over boundary ī and Nȍ over the domain ȍ. Separate relations were established for the boundary condition indicators. Pepper D.W. - Šarler B. Pepper D. W. – âarler B. Strojniãki vestnik - Journal of Mechanical Engineering 51(2005) 483-8, 476-483 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 476-483 The computational domain was discretized into 80 boundary nodes and 361 domain nodes. The multiquadrics constant ro was set to 0.2. Steady state 3 results were achieved after 34 iterations for Ra = 10 , 4 187 iterations for Ra = 10 , and 293 iterations for Ra = 10 5. The calculated values for temperature and velocity were in excellent agreement with results obtained using a fine grid FDM [10]. Acknowledgements We wish to thank Professor C. S. Chen and Professor Jichun Li from the Department of Mathematics at UNLV for their helpful insight and assistance. References 1. R. W. Lewis, K. Morgan, H. R. Thomas, and K. N. Seetharamu, The Finite Element Method in Heat Meshless methods are a unique and novel numerical Transfer Analysis, J. Wiley & Sons: Chichester, UK technique now making inroads into various fields. (1996). Their advantages in solving problems associated with crack propagation and stress/strain including 2. H-C. Huang and A. S. Usmani, Finite Element Analysis for Heat Transfer, Springer-Verlag: London, deformation over more conventional numerical UK (1994). schemes have been demonstrated repeatedly in the literature. The application of meshless methods for 3. G. R. Liu, Mesh Free Methods: Moving Beyond the heat transfer is equally advantageous; such methods Finite Element Method, CRC Press, Boca Raton, FL have become very competitive with both finite volume (2002). and finite element methods for problems involving 4. S. N. Atluri and S. Shen, The Meshless Local Petrov irregular geometries. The requirements for creating Galerkin (MLPG) Method, Tech Science Press, grids as well as the detailed input necessary for Encino, CA (2002). establishing volume or element properties is greatly reduced or eliminated. However, much has yet to be 5. E. J. Kansa, Multiquatric – A Scattered Data Approximation Scheme with Applications to done before meshless methods can handle a wide Computational Fluid Dynamics II, Computers Math. range of fluid flow problems and produce results with Appl., 19, 8/9, (1990), 147-161. confidence and surety. Only a small portion of incompressible flow problems has been addressed; 6. S. N. Atluri, H. K. Kim, and J. Y. Cho, J.Y., A Critical advances are just now being made in the area of Assessment of the Truly Meshless Local Petrovporous media flows. Application to compressible Galerkin (MLPG), and Local Boundary Integral flows has yet to be addressed. Equation (LBIE) Methods, Computational Mechanics, While meshless methods may be more accurate than 24, (1999), 348-372. FDM/FVM/FEM techniques, they can be much slower with regards to computational time to achieve 7. D. W. Pepper, D. B. Carrington, and L. Gewali, A Web-based, Adaptive Finite Element Scheme for Heat convergence. This is due in part to some of the effort Transfer and Fluid Flow. ISHMT/ASME 4th Conf. on needed for numerical integration and subsequent use Heat and Mass Transfer, Jan. 12-14, Pune, India of a direct matrix solver. However, meshless methods (2000). do not need any prior knowledge of their nodal arrangement, as in conventional numerical schemes. 8. F. P. Incropera and D. P. DeWitt, Fundamentals of This makes the method particularly attractive for Heat and Mass Transfer, 5th Ed., J. Wiley & Sons: New developing adaptive capabilities. Since much of a York (2002). modeler’s efforts are generally spent on developing a 9. B. Sarler, J. Perko, C. S. Chen, and G.Kuhn, A good mesh that will lead to a converged solution, the Meshless Approach to Natural Convection, S. N. overall time for obtaining problem solutions using Atluri and D. W. Pepper (Eds.), Proceedings of the meshless methods can be significantly less. 22nd ICES Conference, July 31-Aug. 2, 2002, Reno, Additional information regarding FDM, FVM, and NV (2002). FEM algorithms and some of the meshless techniques can be obtained from the web site 10.D. W. Pepper and C. S. Chen, A Meshless http://www.unlv.edu/NCACM. Method for Modeling Heat Transfer, S. N. Atluri nd and D. W. Pepper (Eds.), Proceedings of the 22 ICES Conference, July 31-Aug. 2, 2002, Reno, NV (2002). Conclusions Application of meshless methods for thermal analysis Application of meshless methods for thermal analysis 483 483