Error estimation of bilinear Galerkin finite element method for 2D thermal problems

This study demonstrates a two-dimensional steady state heat conduction Laplace partial differential equation solution using the bilinear Galerkin finite element method. Heat transfer analysis is of vital importance in many engineering applications and devising computationally inexpensive numerical methods while maintaining accuracy is one of the primary concerns. The method uses structured mesh grid over a two-dimensional rectangular domain and solved using a stiffness matrix for the bilinear elements, calculated using the proposed modified numerical scheme. Several numerical experiments are conducted by controlling the number of nodes and changing element sizes of the presented scheme, and comparison made between analytical solution and software generated solution.


INTRODUCTION
Finite element Analysis is a robust numerical approach for solving partial differential equations (PDEs) used in various branches of science such as solid mechanics, fluid mechanics, electromagnetics, thermodynamics etc. (Kreyszig, 2011;Deb, Babuška, & Oden, 2001). One of the most powerful techniques for solving PDEs with weak formulation is using a weighed residual method called the Galerkin finite element method (GFEM). The formulation requires generating a basis function (Ainsworth & Oden, 1997) based on the elemental boundary conditions. This trial function is substituted in the partial differential equation and the first derivative of the trial function is taken for each nodal variable (Burden & Faires, 2001) to construct the residual function. The weighed form of the residual function for the whole domain is integrated by setting it equal to zero. Green's theorem can be applied over the boundary if necessary (Afzal, Sulaeman, & Okhunov, 2016). The corresponding numerical model is set up by discretizing the rectangular domain into smaller elements where each element consists of nodal coordinates and nodal variables, which are used to perform the Galerkin approximation of the PDEs. This is followed by the generation of an element matrix and vector matrix of the boundary by integrating the total number of elements. The set of linear equations represented by the matrices are consecutively solved using the Galerkin approach. For comparison purposes, an exact solution already available for the non-linear PDE (time independent and no heat source) used in 2D heat conduction rectangular domain with both Dirichlet and Neumann boundary condition, is presented. Finally, a stiffness matrix applicable for homogenous rectangular domain consisting of structured mesh grid elements is presented, the solution scheme of which significantly reduces the CPU performance cost. or Dirichlet boundary conditions as shown in Figure 1 (Gerald & Wheatly, 2003;Gockenbach, 2002;Burden & Faires, 2001;Kreyszig, 2011). Two dimensional time independent heat-conduction problem can be represented by the following partial differential equation.

MATHEMATICAL CONSTRUCTION
(1) where, and The boundary conditions are given by: where u b and f b are Dirichlet and Neumann boundary conditions respectively.
The weighed residual w is applied on Eqs (1) and (2) to generate its strong formulation as given: (3) Consequently, the weak formulation is generated using integration by part over eq. (3) as given below: The stiffness matrix is generated by discretising the solution domain into smaller elements and performing integration for each element.

BASIS FUNCTION AND STIFFNESS MATRIX
The 2D computational domain is represented by bilinear rectangular elements where each element is constructed with a node at each of its four edges as shown in Figure  The above-mentioned basis function can be represented in matrix form: Considering each bilinear element composed of four nodal points and nodal variables, the corresponding basis function in matrix form is represented as: The corresponding shape function is formulated as: The individual shape function is based on the four geometric bilinear coordinates for each element located in a structured mesh grid: Where the rectangle Area, A is given by: As it can be observed from Eq. (5) to (11) In order to determine the elemental stiffness matrix, the first integral of the weak formulation in Eq. (4) is carried out sequentially given below: The matrix generated is a 4x4 matrix as given by: The diagonal values as follows: Since: k 22 =k 11 , k 23 =k 14 , k 24 =k 13 , k 33 =k 11 , k 34 =k 12 and k 44 =k 11 the matrix is symmetric about its diagonal.  Solving the above mentioned problem required calculating the stiffness matrix of the bilinear elements using the Galerkin finite element approach. The accuracy of the approach largely depends on the element size and number of elements considered. It has been observed increasing the number of elements and decreasing the element size increases the accuracy of the solution by achieving higher convergence rates. Comparison is made with the exact solution of the Eq. (23) below and Figure 4 shows the graphical analysis.
(23) The Figure 5 represents the solution domain discretized into a numbers of smaller bilinear elements each consisting of nodes at the corners. The errors analysis of the Galerkin approach solution and analytical solution are given by the following equations:

CONCLUSION
It is shown that GFEM is an inexpensive and powerful technique to accurately model and has demonstrated solving a time dependent and external heat source independent 2D heat conduction Laplace partial differential equation, governing a rectangular planar domain, using bilinear elements.
Increasing the number of elements while decreasing the element size representing the solution domain shows optimal convergence towards the exact solution, thereby validating the accuracy of the scheme.

ACKNOWLEDGEMENT
The work has been financially supported by the ministry of higher education, Malaysia, under the grant FRGS 19-039-0647 and is gratefully acknowledged.