 Engineering, 2009, 1, 151-160 doi:10.4236/eng.2009.13018 Published Online November 2009 (http://www.scirp.org/journal/eng). Copyright © 2009 SciRes. ENGINEERING Heat Distribution in Rectangular Fins Using Efficient Finite Element and Differential Quadrature Methods ShahNor BASRI, M. M. FAKIR, F. MUSTAPHA, D. L. A. MAJID, A. A. JAAFAR Department of Aerospace Engineering, University Putra Malaysia, Putra, Malaysia E-mail: shahnor@eng.upm.edu Received September 2, 2009; revised September 24, 2009; accepted September 28, 2009 Abstract Finite element method (FEM) and differential quadrature method (DQM) are among important numerical techniques used in engineering analyses. Usually elements are sub-divided uniformly in FEM (conventional FEM, CFEM) to obtain temperature distribution behavior in a fin or plate. Hence, extra computational com-plexity is needed to obtain a fair solution with required accuracy. In this paper, non-uniform sub-elements are considered for FEM (efficient FEM, EFEM) solution to reduce the computational complexity. Then this EFEM is applied for the solution of one-dimensional heat transfer problem in a rectangular thin fin. The ob-tained results are compared with CFEM and efficient DQM (EDQM), with non-uniform mesh generation). It is found that the EFEM exhibit more accurate results than CFEM and EDQM showing its potentiality. Keywords: Efficient Finite Element Method, Efficient Differential Quadrature Method, Heat Transfer Problem 1. Introduction Presently there are many numerical solution techniques known to the computational mechanics community. FEM is one of those numerical solution techniques to solve structural, mechanical, heat transfer, and fluid dynamics which arise in problems of engineering and physical sci-ences [1–5]. Here, conventional FEM (CFEM) means the used elements are of same size and uniformly distributed. In its application to the solution of engineering problems, the finite element discretization has been implemented almost to the spatial problems. For dynamic or time de-pendent problems whose solutions as functions of time are of interest, a step by step procedure of finite differ-ence is usually employed with computational complexity. For heat transfer problems, rapid changes of heat/temperature distributions take place near the ele-ment boundary (and at the boundary). It is very impor-tant to know these temperature change behavior of an element prior to its use. Hence, to get an actual picture using FEM, the element is usually subdivided into very small sub-elements uniformly (conventional FEM, CFEM), which leads to huge amount of complexity, memory consumption and computational time . Oth-erwise, error flow occurs with unreliable results [1,2,6]. On the other hand, to get a clear picture about the tem-perature changes near (and at) the element boundary, better to subdivide the elements into very small sub- elements at the boundary only, followed by relatively bigger elements gradually towards the mid-point of the element non-uniformly (efficient FEM, EFEM). This may serve the intended purpose without any additional burden and this is highlighted in this paper with im-proved accuracy (approximately 65%) compared to CFEM. Hence, here, focus is given to develop and apply efficient (non-uniform mesh density) nodal points distri-bution algorithm for automatic mesh (elements) genera-tion to optimize CFEM solution. DQM is another numerical solution technique to solve above mentioned problems efficiently [7–13]. The es-sence of the DQM is that the partial derivative of a func-tion is approximated by a weighted linear sum of the function values at given discrete points. Bellman and Casti [7,8] developed this numerical solution technique in the early 1970s and since then, the technique has been successfully employed in a variety of problems in engi-neering and physical sciences. To make the DQM more efficient with less computational complexity, efficient DQM (EDQM) was proposed in [11–13] with non-uni- formly distributed mesh points. Hence, in this paper, one-dimensional (1-D) heat con-duction problems in a thin rectangular fin are solved us-ing EFEM by means of the accurate discretization and solver (code) and then the results are compared with S. N. BASRI ET AL. Copyright © 2009 SciRes. ENGINEERING 152 CFEM and EDQM to verify EFEM efficiency. The paper is organized as follows. Section II presents the governing equation with efficient FEM rules, fol-lowed by simulation set-up and assumptions, results and discussions, and finally conclusion of the paper. 2. The One-Dimensional Efficient Finite Element Method Here, the considered one dimensional (1-D) heat conduc-tion problem is [2,3,14–18] 0ddTkQdx dx (1) with the boundary conditions 00TT x and (LxLqhTT) as shown in Figure 1. Here, heat flux dtqkdx . Figure 1 shows the 1-D element discretiza-tion in the x-direction. The temperature T at various nodal points are the unknowns except at node 1, where, with initial temperature. Within a typical element ‘01 TT 0Teoreiix12elx’ the local node numbers are iand with coordinates and and element length, . For example, whose local node num-bers are 1.and 2 with coordinates and , and ele-ment length respectively. 1i1eiilxix1x1ix1x1e2xAn one-dimensional thin rectangular fin as shown in Figure 2 is considered here. Heat is transmitted along its length by conduction and dissipated from its lateral sur-faces to the surroundings by convection. The governing equation for the temperature in the fin is given in Equa-tion (1). The parameter, Mis given by2ChpMkA, where, p is the fin perimeter (meter) and Ac is the cross sectional area of the fin [meter2]. Fin length, width and thickness are , and t respectively. LwIn this case, dTqhTT kdx , 2pwt, CAwt and 22CwtpAwt t. The convection heat loss in the fin is equivalent to negative heat source and can be expressed as follows:  ()CCpdxh TTphQTAdx A T Now Equation (1) becomes 0CddTphkTTdx dxA (2) Figure 1. Boundary conditions for 1-D heat conduction. Figure 2. Thin rectangular fin. To calculate the approximate solution T, the mathe-matical formulation using Galerkin’s approach [2,3] is 00LCddTphkTTdxdxdx A (3) where  is a test function constructed from the same basis functions as those of T, with . 00Integrating by parts Equation (3) becomes, 0000LLLCdTd dTphkkdx TTdxdxdx dxA  (4) The 1st term of Equation (4) is,    000 0LdT dTdTkLkLLkdx dxdx  (5) Since 00 and  LdTqkL LhTTdx  , we get,  0LLdTkLhTdx T Equation (4) becomes  000LLLCddT phLhTTkdxT Tdxdx dxA  (6) A global virtual temperature vector is defined as then within each element, th 123,,,..., TL  S. N. BASRI ET AL. Copyright © 2009 SciRes. ENGINEERING 153test function becomes According to Reference ,eTdTdx BT , we have TddxB () iiiN (7) Here, is the element shape function and N1LNat the element boundary  (Figure 1). Therefore Equation (7) gives For matrix multiplication validity, we have TTT eiTddTdx dx TBψBTLLLN (8)  and 2211 111111 1111111() ()TTii iiiieixx xxxxl 1111  TBB  The element conductivity matrix is The element heat rate vector due to the heat source is 1111112ei eieiTTTeikl kkdlTBB  (9) 111122eieiei eiQQl Qld TRr N (10) where,  varies from to 11 and 1122()1iiiNow, Equation (2) can be transformed into xiixwith ddxxx xx. 1111022eei eiei eilL TTei eikl QlhT Tdd TTTTψBB TψN (11) or TLL LhT hTTTψKT ψR (12) where, global matrices KT , R, and T are respec-tively, Tψ1112 13 12122 23 213132 33 31122.........2.......................................................LLei eiTTT LeiLLLLKKKKKKKKkl dKK K KKKKKTKBB212223221 02231323333331 041424344441 0123 10.................................................. ....LLLLLLL LLLLKKK KKTTRKKKKT RKTKKKKT RKTTRhTKKKK hKT  (16) L (13) This equation needs to be solved to obtain the 1-D FEM numerical temperature distribution in the consid-ered rectangular fin. 11231...2ei eiLeiQl dRRRRTTRN (14) Using Equations (11-16) and the efficient FEM (EFEM) algorithm, the approximate solution T has been obtained. The 1-D EFEM algorithm (rule) is depicted in terms of self-explanatory flow chart in Figure 3. The non-uniform and uniform mesh distribution scenarios are shown in Figures 4 and 5 respectively. 23000...00 000...0000...00010... 0000...00001... 00................... ............. ................ ..................... ............ ................. ......000... 01000...0LTψ 2.1. Example Problem 1: 1-D Insulated Tip Thin Rectangular Fin (15) and with constant. 123... LTTT TTT10TT When the base of the fin is held at constant temperature, T0 and the tip of the fin is insulated, the boundary condi-tions are then given by Finally, the global matrix is formed and the Equation (12) becomes S. N. BASRI ET AL. 154 Start and Initialization Input: Fin Length: L, no. of element: N, error threshold: eh Figure 3. Efficient discretization and solution rule for 1-D FEMNo Yes No No Yes No. nodal point: Z=N+1 Calculation of mesh distribution i = 1 to Z 1()1 cos21Lixi ZElement length. i = 1 to N le(i) = x(i+1) – x(i) Non-uniform ? No nodal points: Z=N+1.Element length: le = LMesh distribution calcula-tion i = 1 to z, x(i) = (i – 1) le Numerical solution and error calculation Max. |Tn – Texact| ≤eh? Set N=N+2 Discretization and Stiff-ness matrix calculation using Galerkin approach Set N=N+2 END Copyright © 2009 SciRes. ENGINEERING S. N. BASRI ET AL. 1550TT at 0x0q atxL, where L is the length of the fin. In this case, the final form of the global matrix in Equation (16) becomes 2223221 022323333331 023 1............... ...............LLLL LLLLLxx = Lx = 0 x x = L x = 0 0AAAATR TAAATRATAAATRA       T0 (17) Example of non-uniform and uniform mesh distribu-tions and element lengths are depicted in Figures 4 and 5 respectively. 2.2. Example Problem 2: 1-D Convection Tip Thin Rectangular Fin When the base of the fin is held at a constant temperature, T0 and the tip of the fin is a convection surface, then the boundary conditions are T=T0 at x=0 LqhT T at x=L And the final global matrix shown in Equation (16) becomes 2223221 022323333331 023 1............... ...............LLLL LLLLLAA AATTRAAA TRATAA AhTRhTAT (18) 3. Simulation Set-up and Assumptions Table 1 shows the considered parameters and their cor-responding values used to obtain simulation results using FORTRAN 90 software. We used these values to obtain the temperature distribution for EFEM, CFEM, EDQM and exact methods. Figure 4: Example 1-D efficient mesh distribution and ele-ment lengths. Figure 5. Example 1-D conventional mesh distribution and element lengths. We considered, 21hPMkA and the associated assump-tions (in Table 1) to compare the obtained FEM results with DQM  and exact solution . Here to mention that, to obtain 1-D DQM solutions, element material properties, fin-width and fin-thickness are not required (which is the shortcoming of the method). The errors in FEM and DQM solutions are computed compared to exact solution . 4. Results and Discussions 4.1. Results and Discussions of 1-D Insulated Tip Thin Rectangular Fin The results of the present problem, shown in Figure 6, contain the maximum absolute percentage errors in the FEM and DQM solutions obtained with uniformly (con-ventional) and non-uniformly (efficient) distributed no- dal (mesh) points. It is essential to know, how many mesh points (elements) are required to obtain a conver-gent FEM solution in the solution domain. Hence, the comparison of convergence of fin-tem- perature in terms of maximum % error versus number of nodal (mesh) points for CFEM, EFEM and EDQM solu-tions is shown in Figure 6. Initially, all the solutions in terms of maximum % errors show a monotonic conver-gence with the increasing number of mesh points (shown 11to 104Z). It is apparent that EFEM results show bit less accuracy for 30Zand similar accuracy for compared to EDQM, but yields result with higher accuracy, of one order of magnitude or more with in-creasing 30ZZ compared to CFEM. EDQM converge up to 100Zand then saturated, whereas the EFEM solutions converge smoothly for all N within the solution domain, showing best converging result at . On the other hand, uniform FEM (CFEM) results converge slowly throughout the solution domain and then diverge without showing the best results like EFEM. It happens due to the mesh point distribution strategy of equally spaced and unequally spaced nodal points in the compu-tational domain and the inherited complexity to compute the stiffness matrix for equally spaced nodal points. Hence, the efficiency of EFEM results is apparent. 100and 101ZFigure 7 shows the convergent numerical and exact solutions (fin temperature) and the corresponding per-centage errors for 100N elements (FEM case) which is equivalent to 101Z mesh points (both FEM and DQM cases). These results are obtained at an interval of 0.1x along the fin length,01x, using cubicCopyright © 2009 SciRes. ENGINEERING S. N. BASRI ET AL. 156 Table 1. Input parameters and assumptions for 1-d rectangular fin. spline interpolation. It is seen that all the solutions are very close to exact solutions throughout the length of the fin with temperature variations at 001TC0xm to at00.648LTC1.0xm. Figure 8 shows the percentage errors at the base of the fin () are 0 for all solutions due to initial tempera-ture (Figure 7). The percentages errors remain the same with EFEM except little bit increase (with maximum error) at the middle of the fin due to nodal point distribution with maximum spacing there. Whereas, with CFEM, it increases gradually along the length of the fin with the maximum percentage error at the fin-tip (x = 1). In other case, the oscilla-tions (instability) of DQM results appear clearly com-pared to FEM results. The average percentage error in 0x01T410 0C62.44 101.9CFEM, EDQM  and EFEM are41.2 10, and respectively, which shows approximately 99% and 49% improvements in EFEM results demonstrating its superiority over CFEM and EDQM. 62.24 1061.12 104.2. Results and Discussion of 1-D Convection Tip Thin Rectangular Fin Here the results exhibit the same nature like insulated-tip fin but yield results with higher accuracy, of two order of magnitude or more with increasing Z due to different material properties and fin-thickness (as the FEM solu-tion and its accuracy depend on fin dimension, materials used and associated boundary conditions). In Figure 9, the comparison of convergence versus number of mesh points of exact, FEM and DQM solu-tions for convection-tip fin with uniform and non-uni-form mesh distributions is shown. It is apparent that for all cases, the solutions converge smoothly for all Zwithin the solution domain. The comparison shows similar results as in Figure 6 except EFEM yields result with higher accuracy, of one order of magnitude or more with increasing Z (for ) compared to that with CFEM. Here, EFEM results converge from 20Z80Z showing best result at90to 101Z, EDQM  shows similar results with some oscillations, whereas CFEM does not exhibit any best convergence. Input Parameters Assumed value for Insulated-Tip Fin Assumed value for Convec-tion-Tip Fin Boundary and other values: Initial temperature (T0) Ambient temperature (T∞) Heat flux (q) % Error threshold (eh) 1 OC 0 OC 0 at x = 1 0 - 0.1 1 OC 0 OC Variable 0 - 0.1 Element Type (NNODE): Linear (for 1-D) 2 2 Element material properties: Thermal conductivity (ke = k) Convective heat transfer coefficient (h) Heat source (Q) Variable to make M = 1 9 W/m2 0C 0 W/m3 0C 7.03125 W/(m 0C) 9 W/m2 0C 0 W/m3 0C Element (Fin) dimension: length (L) along x-axis width (w) thickness (t) Number of elements (N) 1 m Variable to make M = 1 0.001 m 11 - 104 1 m Variable to make M = 1 Variable to make M = 1 11 - 104 Copyright © 2009 SciRes. ENGINEERING S. N. BASRI ET AL. 157 Figure 6. Comparison of convergence of insulated-tip fin-temperature in terms of maximum % error for CFEM, EFEM and EDQM solutions (). 11to 104Z Figure 7. Insulated-tip fin-temperature distribution for exact, EFEM, CFEM and EDQM along with its respective % errors (). Z=101Figure 10 depict the comparison of CFEM, EFEM, EDQM  numerical and exact convection-tip fin tem-peratures and the corresponding percentage errors for conventional (uniform) and efficient (non-uniform) mesh point distribution respectively for 100 elements (i.e., ). Same as Insulated-tip fin, the results are ob-tained at an interval of  along the fin length, 101Z00.1x1x, using cubic spline interpolation. Figure 10 shows that, all numerical solutions are very close to ex-act solutions throughout the length of the fin with tem-perature variations 0at base of the fin to at the tip of the fin. Here the reduction of fin temperature is more compared to insu-lated-tip fin (Figure 7) as expected. 01T00.32 CCC00.328LTFEM versus DQM maximum % error comparison for convection-tip fin-temperature are shown in Figure 11. The comparison of CFEM, EFEM and EDQM  per- Copyright © 2009 SciRes. ENGINEERING S. N. BASRI ET AL. 158 Figure 8. Percentage error comparison of EFEM, EDQM and CFEM for along the fin-length. Z=101 Figure 9. Comparison of convergence of convection-tip fin-temperature in terms of maximum % error for CFEM, EFEM and EDQM solutions (). Z=11to 104 centage errors for convection-tip fin is shown in Figure 11. There is no error at the base of the fin and it almost remain the same with EFEM and EDQM except negligi-ble increase at the middle of the fin, whereas, with CFEM, it increases gradually along the length of the fin with the maximum percentage error at the tip (x = 1). In this case the EDQM converges with oscilla-tions throughout the solution domain. The average % error in CFEM, EDQM  and EFEM are 63.31 101.69 610, and respectively. This shows nearly 100% and 99% improvements in EFEM results 93.08 10112.24 10compared to CFEM and EDQM respectively demonstrat- ing its superiority. 5. Conclusions Here, the solutions of the temperature distribution in in-sulated-tip and convection-tip 1-D rectangular fin are computed numerically using FEM and the results are found to agree very well with the exact solution and show the efficiency of the method. Investigating the various mesh points distribution for equal and unequal spacing, it is found that, for FEM solution, unequally spaced mesh points distribution give better and more Copyright © 2009 SciRes. ENGINEERING S. N. BASRI ET AL. 159 Figure 10. Convection-tip fin-temperature distribution for exact, EFEM, CFEM and EDQM along with its respective % er-rors (Z). =101 Figure 11. Convection-tip fin error comparison of EFEM, EDQM and CFEM for along the fin-length. Z=101 accurate results than equally spaced and the solution converges smoothly as the number of nodal points (or elements) is increased. In general, this study has im-proved the stability and accuracy of EFEM results for practical consideration and implementation. Finally, the results of EFEM shows remarkable en-hancement compared to CFEM and agree very well with EDQM with very small errors or difference showing its potentiality. Hence EFEM is suitable to test the tem-perature distribution scenario in any thin metal fin prior to its design and practical implementation. 7. References  G. Strang and G. J. Fix, “An analysis of the finite element method,” Prentice-Hall, Inc., 1997.  R. Tirupathi and P. E. Chandrupatla, Introduction to finite elements in engineering, Prentice-Hall International, 1997.  “Mathematics of finite element method,” 2004, Available at: http://math.nist.gov/mcsd/savg/tutorial/ansys/FEM/ (Accessed on October 19, 2006).  B. Li, “Numerical method for a parabolic stochastic partial differential equation,” Chalmers Finite Element Center, 2004. Copyright © 2009 SciRes. ENGINEERING S. N. BASRI ET AL. Copyright © 2009 SciRes. ENGINEERING 160  F. Fairag, “Numerical computations of viscous, incom-pressible flow problems using a two-level finite element method,” arXiv: Math.NA/0109109, Vol. 1, No. 17, Sep-tember 2001.  S. Park, “Development and applications of finite elements in time domain,” PhD Thesis, Faculty of the Virginia Polytechnic Institute and State University, Virginia, De-cember, 1996, Available in the Library, Virginia Poly-technic Institute and State University.  R. Bellman and J. Casti, “Differential quadrature and long-term integration,” J Math and Appl., Vol. 34, pp. 235–238, 1971.  R. Bellman and J. Casti, “Differential quadrature: A tech-nique for the rapid solution of nonlinear partial differen-tial equations,” J. Comput Phys, Vol. 10, pp. 40–52, 1972.  C. W. Bert, S. K. Jang, and A. G. Striz, “Nonlinear bend-ing analysis of orthotropic rectangular plates by the method of differential quadrature,” J. Comput Mech, Vol. 5, pp. 217–226, 1989.  C. W. Bert and M. Malik, “Fast computing technique for the transient response of gas-lubricated journal bearings,” Proceedings of the U.S. National Congress of Applied Mechanics, Seattle, WA, pp. 298, June 26-July 1, 1994.  C. Shu, W. Chen, H. Xue, and H. Du, “Numerical analy-sis of grid distribution effect on the accuracy of differen-tial quadrature analysis of beams and plates by error es-timation of derivative approximation,” Int. Journal of Numer, Methods Engineering, Vol. 51, No. 2, pp. 159–179, 2001.  M. M. Fakir, M. K. Mawlood, W. Asrar, S. Basri, and A. A. Omar, “Triangular fin temperature distribution by the method of differential quadrature,” Journal Mekanikal, Malaysia, Vol. 15, pp. 20–31, 2003.  M. M. Fakir, M. K. Mawlood, W. Asrar, S. Basri, and A. A. Omar, “Rectangular fin temperature distribution by the method of differential quadrature,” The Journal of the In-stitute of Engineers, Malaysia (IEM), Vol. 63, No. 4, pp. 41–47, 2002.  E. Hinton and D. R. J. Owen, “An introduction to finite element computations,” Pineridge Press Limited, Swan-sea, U.K., 1985.  S. Tiwari, G. Biswas, P. L. N. Prasad, and S. Basu, “Nu-merical prediction of flow and heat transfer in a rectan-gular channel with a built-in circular tube,” ASME Jour-nal of Heat Transfer, Vol. 125, No. 3, pp. 413–421, 2003.  B. L. Wang, and Z. H. Tian, “Application of finite ele-ment-finite difference method to the determination of transient temperature field in functionally graded materi-als,” Journal of Finite Elements in Analysis and Design, Vol. 41, pp. 335–349, 2005.  S. H. Lo and W. X. Wang, “Finite element mesh genera-tion over intersecting curved surfaces by tracing of neighbours,” Journal of Finite Elements in Analysis and Design, Vol. 41, pp. 351–370, 2005.  M. N. Ozisik, “Heat transfer: A basic approach,” McGraw-Hill, 1985.