A new finite element method using MATLAB for Poisson’s Equation in axisymmetric domain

Article History:Received:11 January 2021; Accepted: 27 February 2021; Published online: 5 April 2021 Abstract: This paper demonstrates finite element procedure for two-dimensional axisymmetric domains. For many engineering applications like structural engineering, aerospace engineering, geo-mechanics etc., the solution domain and boundary conditions are axisymmetric. Henceforth, we can illuminate just the axisymmetric part of the solution domain that gives the data of the entire domain. This paper demonstrates the effectiveness of using MATLAB programming demonstrated by Persson et.al (2004) as the initial mesh for discretization of axisymmetric domains for higher order meshing. Further, solving some class of partial differential equations using finite element method with nodal relation given by subparametric transformations Rathod et.al (2008). In this paper a cubic order curved triangular meshing for some of the domains like ellipse and circle are demonstrated. These in turn finds its applications in the fields like stress analysis in mechanical engineering, torsion twist (shear strength) analysis in civil engineering, evaluation of stress intensity factor for quarter elliptical crack in pressure vessels in equipment industry etc,. The output data from the meshing scheme like meshing of the domain, nodal position, element connectivity and boundary edges are been used in the finite element procedures. The efficiency of the method is achieved by p-refinement scheme i.e., fixing the number of elements and increasing the polynomial order.


INTRODUCTION
The mathematical models where developed in order to analysis the physical phenomenon that depicts the system with specific assumptions and simplifications. This leads to governing mathematical expressions, which represent the behaviour of the system. The mathematical expressions are usually differential equations with boundary conditions. It is tedious to obtain the analytical solution for problems involving complex material properties and boundary conditions; hence, we resort to numerical methods that provide approximate but acceptable solution. Various numerical methods where developed to obtain the approximate solution for the differential equations. Especially, finite element method (FEM) is widely used numerical methods to solve the physical problems. The major advantages of using FEM are generalized computer program can be developed to analyse various problems and it can handle any complex geometry with the given boundary conditions. The strategy of FEM is based on meshing the complex geometry into simple elements such as triangle, quadrilateral or polygonal in 2D and tetrahedron or hexahedron in 3D. The investigation states that various endeavours were made on meshing the geometry with linear straight edged triangular, quadrilateral and polygonal elements [1,2,3]. Since, the domain of physical problem often contains curved boundaries that cause the difficulty in accurately meshing with straight edged elements. Hence, the requirement of meshing curved boundaries emerged. In addition, the use of transformations is essential to enhance the quality of meshing arbitrary domains, especially the curved boundaries. In paper [4,5] the use of isoparametric transformations by parabolic or cubic arcs to match the curved boundaries using curved triangular elements where developed. Generally, usage of isoparametric transformations encounters a complex rational integral function whose denominator and Jacobian are bivariate polynomial of higher order. Thus the difficulty arising in solving such an integral of an arbitrary curved boundaries can be overcome by the use of subparametric transformations by parabolic or cubic arcs for higher order curved triangular elements developed in the works [6,7]. Hence, the Jacobian acquired will always be bivariate polynomial of linear order. In the light of subparametric transformations, curved boundaries are meshed effectively. Numerous computer programming exist which executes meshing of two-dimensional (2D) geometries using straight edges, especially triangular elements. Persson and Gilbert [1] have presented a simple and efficient MATLAB meshing code for linear order straight edged triangular elements. Koko [8] has demonstrated a new algorithm improving the MATLAB code of Perrson et.al [1] for fast refinement upto quadratic order straight edged triangular elements. The implementation of the MATLAB code [1] for finite element procedures upto cubic order for straight edged elements is provided in the works of John Coady [9]. Though there are some excellent tools available for meshing, there is necessity for meshing in higher order curved triangular elements. Some of the two dimensional geometries for higher order curved triangular elements are meshed in the works of [10,11]. In the present work, a simple and efficient MATLAB code is proposed based on distmesh2d developed by Persson et.al [1]. The proposed meshing scheme uses the finest point transformations proposed in the works of Rathod et.al [6]. The proposed meshing scheme is implemented in discretization of some of the 2D axisymmetric geometries. As an output meshing scheme provides triangular meshing of the geometry, nodal position, element connectivity and boundary nodes, which can be used in the finite element procedures. The mathematical formulation of the finite element method is provided in section 2. The formulation of meshing cubic order curved triangular element using the subparametric transformations is briefed in section 3. Discretization of some of the axisymmetric geometries is demonstrated in section 4. The implementation of meshing scheme for the solution of Poisson's problem is provided in section 5. The numerical results and the comparison of contour plot for the stress analysis are provided in section 6.

MATHEMATICAL FORMULATION OF FINITE ELEMENT METHOD
Finite element method is a computational methodology used for numerical approximation of some class of partial differential equations that models the problems arising in applied sciences. The strategy of FEM is based on meshing the complex domains into simple elements like triangular or quadrilateral in 2D. PDEs are used to describe a good number of phenomenons in sciences and engineering. One such phenomenon in 2D boundary value problem is the stress concentration over elliptical or circular plate [12,13]. Consider a two-dimensional boundary value problem, subject to the boundary conditions, Essential boundary condition: = known value on Ω. (1a) Natural boundary condition: where Ω and Ω are the interior and boundary region of the domain respectively. (2) where refers to Lagrange shape functions, refers to the unknown nodal values and represents the order of the curved triangular elements. The underlying mathematical basis of FEM is very well explained by the theories of Galerkin weak formulation [6,7]. A key feature of this method is that the integrals of functions can be evaluated on any arbitrary domains. The integrations are carried out as the product of shape function and governing differential equation over the given domain and it is equal to zero. Now, on applying the Galerkin weak formulation to Eq. (1), we have where varies from 1 to . Any natural boundary conditions Eq. (1b) are imposed in the weak formulation. Hence, the finite element equation is as follows, where Ω is the region of each element and is the determinant of Jacobian matrix which is expressed as, The above integrals Eqs. (5) -(8) are evaluated numerically for each element by using Gauss quadrature rules over a standard triangle from the works of Rathod et.al. [14,15]. Then assembling is done to get the global matrices relating the entire domain. Any essential boundary conditions Eq. (1a) are imposed to the global matrices. The obtained system of equations are solved numerically to get the solution at each unknown nodal points.

MATHEMATICAL FORMULATION FOR MESHING OF CUBIC ORDER CURVED TRIANGULAR ELEMENT
Generally, in FEM while solving PDEs for curved boundaries encounters a complex rational integral function whose denominator and Jacobian are bivariate polynomials of higher order as discussed in [6]. Thus, the difficulty arising in solving such an integral function of an arbitrary curved boundaries can be overcome by the use of subparametric transformations with parabolic arcs developed by [6]. The Jacobian acquired will be a bivariate polynomial of linear order. In light of subparametric transformations for higher order, curved boundaries can be meshed effectively. We use [1] as the starting mesh by defining the nodal points on the interior and boundary by subparametric transformations matching parabolic arcs given in [6]. The investigation states that various endeavours where made on subparametric transformations by several authors [6,7]. As a part of this section, the way to generate curved cubic order triangular meshing is discussed. Considering, the mapping of global co-ordinates ( , ) in one-sided curved triangular element to the local co-ordinates ( , ) in the standard triangular element as evident in Fig.1. The subparametric transformations defining the nodes on the boundary and interior of the curved triangular element is represented by, where and denotes the Lagrange shape functions and nodal values for each element at ℎ node respectively. The shape functions for curved triangular element in the case of cubic order are as per the following,  The subparametric transformations which maps the global coordinates of two sided straight and one sided curved triangular element to the local coordinates of the standard triangular element for cubic order is given in Fig.1.
The above subparametric transformations is incorporated into MATLAB code in terms of and co-ordinate.

DISCRETIZATION DETAILS OF AXISYMMETRIC DOMAINS
Meshing is one of the important steps in FEM especially for some arbitrary curved domains. The proposed meshing scheme technique is based on a distmesh2d MATLAB code developed by [1] for linear triangular elements which is simple and produces high-quality meshes. The geometrical description as defined by [1] of the domain use up few terms as input such as signed distance function , edge length function ℎ, mesh size ℎ0, bounding box , fixed nodal position and additional parameters to the function and ℎ as varargin arguments. In addition, the scheme uses the nodal relationship from subparametric transformations developed by [6] for curved triangular elements. The physical domain under study in the field of FEM may possess axisymmetric geometry. Since, the solution domain and boundary conditions are axisymmetric we can illuminate just the axisymmetric part of the solution domain that gives the data of the entire domain which can be expected to give symmetric solution. The meshing has been carried out for p-refinement scheme i.e., keeping the number of elements same and increasing the order of the polynomial. Some of the axisymmetric geometries are discertized for cubic (10 -noded) ordered curved triangular elements in Figs. 2 -3. Meshing of the domain finds advantageous in the fields like stress analysis in structural industry, torsional twist (shear strength) in manufacturing industry, evaluation of stress intensity factor for quarter elliptical crack in pressure vessels in equipment industry, Prandtl stress distribution [13] etc.

NUMERICAL RESULTS AND INTERPRETATION
The present meshing scheme produces output data like nodal position, element connectivity and boundary data. These outputs for quarter circle has been extracted and given in Fig. 4 and Table 1 for cubic order curved triangular elements. Fig. 4 gives the discretization of quarter circle for cubic order with node numbers. Table 1 gives the meshing quantities like number of nodes , nodal position , element connectivity and boundary points . One of the aims of this work is to solve the PDEs using output data extracted from the meshing scheme explained in this paper by finite element procedure. The numerical experiment detailed in section 4 governed by Poisson's equation for two domains; quarter ellipse and quarter circle are studied. Both the domains are discretized with linear, quadratic and cubic order curved triangular elements. The 2 norm, which gives the overall error in the domain and CPU time, which depicts the efficiency of finite element procedure are tabulated in Tables 2 -3.  Table 2: Numerical results over quarter ellipse upto cubic order triangular element with mesh size ℎ0 = 0.9 and number of elements = 48.
Order ( )  5b. Exact solution for cubic order. 6b. Exact solution for cubic order. Figure 6: Contour plot for cubic order curved triangular elements for FEM and exact solution of circular domain.

Conclusions
The conclusions drawn on the numerical experiments are as per the following, (i) The present meshing scheme as utilized best discretization procedure for axisymmetric geometries like circle and ellipse in Figs. 2 -3.
The meshing quantities like nodal position ( ), element connectivity ( ) and boundary points ( ) are obtained as an output from the proposed meshing scheme. And these meshing quantities are given in Fig.4 and Table 1 for cubic order. (iii) The finite element procedure uses the acceptable quadrature rule and finest subparametric transformations for obtaining the solutions of Poisson's equation with the output extracted from meshing scheme. (iv) The accuracy and efficiency of this meshing scheme is evident from Tables 2 -3. We have obtained 2 norm upto three decimal accuracy for quarter ellipse and four decimal accuracy for quarter circle for cubic order curved triangular elements. The CPU time (in seconds) depicts the efficiency of meshing. (v) The exact solution matches very well with the FEM solution can be compared by contour plots Figs. 5a -5b for quarter ellipse and Figs. 6a -6b for quarter circle.
Subsequently, it is evident from the proposed scheme that it is one of the simple and computationally efficient meshing scheme for higher order curved triangular elements. These in turn finds its applications in stress analysis in mechanical engineering, torsion twist (shear strength) analysis in civil engineering, evaluation of stress intensity factor for quarter elliptical crack in pressure vessels in equipment industry, Prandtl stress distribution etc.