Generalization of FEM Using Node-Based Shape Functions

: In standard FEM, the stiffness of an element is exclusively influenced by nodes associated with the element via its element-based shape functions. In this paper, the authors present a method that can be viewed as a generalization of FEM for which the influence of a node is not limited by a hat function around the node. Shape functions over an element can be interpolated over a predefined set of nodes around the element. These node-based shape functions employ Kriging Interpolations commonly found in geostatistical technique. In this study, a set of influencing nodes are covered by surrounding layers of elements defined as its domain of influence (DOI). Thus, the element stiffness is influenced by not only the element nodes, but also satellite nodes outside the element. In a special case with zero satellite nodes, the method is specialized to the conventional FEM. This method is referred to as Node-Based Kriging FEM or K-FEM. The K-FEM has been tested on 2D elastostatic, Reissner-Mindlin’s plate and shell problems. In all cases, exceptionally accurate displacement and stress fields can be achieved with relatively coarse meshes. In addition, the same set of Kringing shape functions can be used to interpolate the mesh geometry. This property is very useful for representing the curved geometry of shells. The distinctive advantage of the K-FEM is its inheritance of the computational procedure of FEM. Any existing FE code can be easily extended to K-FEM; thus, it has a higher chance to be accepted in practice.


Introduction
Simulation of physical phenomena is very useful and important both in academic researches and in industrial product designs. The underlying mathematical models of the simulation are usually so complex that it is very difficult or even impossible to obtain analytical solutions. Thus, numerical methods have become indispensable in simulations. Among various numerical techniques, the finite element method (FEM) has been widely used in industries. Its versatility and robustness have been tested by several decades of real engineering practices.
Motivated by the desire to minimize efforts in preparing finite element meshes, various mesh-free methods have been proposed. Their common advantages are as follows: (1) No element mesh is required for the construction of approximate functions; (2) High-order continuity of the approximate functions can be achieved; (3) Superior performance can normally be expected over the standard FEM. A detailed review is presented in [1][2][3].
Among countless mesh-free methods, the authors were interested in the methods of which formulation basis is the same as that of the FEM, i.e., those employing a global Galerkin weak form. 1  One earliest mesh-free method in this category is the element-free Galerkin methods (EFGM) presented by T. Belytschko et al. in 1994 [4], which is an improved version of the diffuse element method proposed by B. Nayroles et al. two years earlier [5]. The mesh-free character of the EFGM is made possible by the use of moving least-squares (MLS) approximant for the test and trial functions in the Galerkin weak form. The MLS approximation is essentially a least-squares regression with a local weighting function. Therefore, it is generally not passing through the data nodes. In other words, MLS shape functions do not possess the Kronecker delta property. Because of this, the enforcement of essential boundary conditions has been a major issue in the EFGM; a special constraint technique must be utilized to impose essential boundary conditions.
In 2003, L. Gu [6] proposed an EFGM with moving Kriging (MK) interpolation to replace the MLS because of its two key properties: the Kronecker delta property and the consistency (polynomial reproducing) property. Following this work, P. Tongsuk and W. Kanok-Nukulchai [7] in 2004 found that, with the same number of nodes, the EFGM with MK interpolation consistently outperformed the original EFGM in terms of accuracy. A further application of the method to shell problems was presented by V. Sayakoummane and W. Kanok-Nukulchai [8] in 2007.
Even though EFGM is claimed to be "element free", a mesh of background cells, a term used to differentiate from "elements", is still needed for numerical integration. In problems dealing with material and geometric discontinuities, the need for a mesh to outline these discontinuities is practically unavoidable. Another disadvantage of the EFGM and its variants is the difficulty in their implementation based on existing general purpose FEM codes. Due to these inconveniences, their acceptance in real engineering practices seems to be unsatisfactory.
In 2005, K. Plengkhom and W. Kanok-Nukulchai [9] proposed a more convenient implementation of the EFGM with Kriging interpolation (KI). In their method, the field variables (trial and test functions) are approximated by "element-by-element" piecewise KI. For each element, KI is constructed from a set of nodes in its domain of influence (DOI) defined over surrounding layers of elements, as illustrated in Figure 1. Like FEM, elements are also used as subdomains for numerical integration. The method is named Kriging-based FEM (K-FEM). This variant of EFGM can be viewed as a generalization of FEM for which the influence of a node is not limited only to hat functions. In standard FEM, the element stiffness is exclusively influenced by its element nodes, whereas the element stiffness in K-FEM can also be influenced by satellite nodes not directly connected to the element. If we limit the DOI to only one element layer with no satellite nodes, K-FEM is then identical to the conventional FEM.

Kriging Interpolation
Named after Danie G. Krige, a South African mining engineer, Kriging is a well-known geostatistical technique for spatial data interpolation in geology and mining. Using this interpolation, unknown at any point can be interpolated from known values at scattered points in its specified neighborhood. The basic concepts of the KI in the context of K-FEM are presented in the following. A detail explanation and derivation of Kriging can be found in the geostatistics literatures (e.g. [10,11]).
Consider a two-dimensional domain modeled by a mesh of triangular elements (Figure 1). Suppose there is a single field variable over the domain, u(x). For each element, the KI is constructed over a set of nodes in a sub-domain E    encompassing a predetermined number of layers of elements. The KI over sub-domain E  can be expressed in the usual is the 1 n  matrix of Kriging shape functions and d is the 1 n matrix of field values at the nodes. In contrast to the FEM, here n is not necessarily only the number of nodes associated with the element, but also includes all its satellite nodes.
In Kriging formulation, the field variable u(x), which is a deterministic function, is viewed as the realization of a random function U(x). The shape function matrix can be expressed as From the above formulation, constructing Kriging shape functions requires a polynomial basis function and a correlation function. For the basis function, besides complete polynomial bases, it is also possible to use incomplete polynomial bases such as bi-linear, bi-quadratic and bi-cubic bases. A widely used correlation function in the area of computational mechanics is the Gaussian correlation function [6][7][8][9]. This function contains an important parameter affecting the quality of KI, known as the correlation parameter θ. In order to obtain reasonable results in K-FEM, K. Plengkhom and W. Kanok-Nukulchai [9] suggested a criterion for choosing a stable range of θ.
Recently the authors introduced a new correlation function [12] in the form of a quartic spline (QS) correlation function. Our studies indicate a superior performance of QS to the Gaussian correlation function, as the resulting Kriging shape functions are less sensitive to the change of θ.
In Figure 1, to illustrate the concept of elementlayered DOI, suppose that the element of interest in a square domain is Element 1, the choices of DOI, comprising one up to four element layers, are shown in the Figure 1. It is noted that the DOI does not have to be convex. If one uses quadratic basis function (m=6) and choose to use three-layered DOI to construct KI over Element 1, the DOI will encompass 30 (n=30) nodes. The plot of Kriging shape function associated with node I, based on QS correlation function, is shown in the right-hand side of Figure 1.

Key Advantages of K-FEM The Stress Field can be Obtained with Remarkable Accuracy and Global Smoothness
Using the same mesh size, K-FEM yields a stress field with higher accuracy and better smoothness than that of the standard FEM. This is because one can freely adopt a higher-order basis function and a larger DOI for any fixed mesh. To show this, a cantilever plane-stress beam, Figure 2, under end parabolic shear is modeled with a crude mesh of 6x10 triangular elements. In the same figure, the quality of stress output obtained by K-FEM using cubic basis and three-layered DOI is demonstrated by the stress contours generated directly from nodal values with no post-processing manipulation. Like FEM, there is no guarantee for stress field to be perfectly continuous across the inter-element boundaries; however, the degree of discontinuity is found to be rather insignificant.

Solution Refinements can be Achieved with no Re-meshing
In K-FEM, quality improvement of solutions can be achieved by: (a) increasing the order of the basis function or p-refinement, or (b) enlarging the element-layered DOI or l-refinement. For illustration, the cantilever plane-stress beam is modeled with 3 mesh sizes, i.e., with 6x10, 12x20 and 24x40 triangular elements. Each mesh is tested with linear (P1), quadratic (P2) and cubic (P3) polynomial basis functions. For P1, it is possible to use 1, 2, or 3 element layers for the DOI. However, at least 2 layers must be used for P2 and at least 3 layers for P3, following the general rule that the number of nodes covered in the DOI must not be fewer than the number of terms in the polynomial basis. Results of the end deflection, normalized by the exact solution, for all cases are presented in Table 1 together with the corresponding computational times

Normal Stress Shear Stress
Accuracy performance and computational times over the matrix of the h-refinement and the l-refinement, all using linear basis function, are presented in Figure 3. For relatively crude meshes, the accuracy can be enhanced by adopting a larger DOI with more layers of elements. Almost the same accuracy can be achieved by h-refinement from 6x10 to 24x40 mesh sizes, or by l-refinement from 1 to 3 element layers. The latter requires about 20% more computational time. However for the case of h-refinement, we do not consider engineer's time for the remesh. A more detailed comparison of beam displacement profile between h-refinement and l-refinement is illustrated in Figure 4. Accuracy performance and computational times over the matrix for h-refinement and p-refinement, all using DOI of three element layers, are presented in Figure 5. From the figure, higher accuracy can be achieved for a fixed mesh by simply adopting a higher order basis function without significantly increasing the computing time.

Geometry of Curved Domain can be Represented More Accurately by KI Isoparametric Mapping
The same set of Kriging shape functions for field variable can be used to interpolate the geometric field. This is very useful for curved shell problems. To demonstrate this advantage, a cantilever quarter cylinder shell under pure bending is modeled by triangular elements as shown in Figure 6. K-FEM is used to solve the shell problem with quartic basis functions and a DOI of 4 element layers. This shell problem will be tested for two different situations, one with and the other without isoparametric mapping. In the first case, the geometry of individual shell elements shall be interpolated by Kriging shape functions. In the latter case, the geometry of individual shell element is basically a flat facet. The results clearly confirm the advantage of the Kriging interpolated shell geometry.

Implementation of K-FEM can be Easily Incorporated into Existing FEM Codes.
As K-FEM inherits the computational procedure of FEM, existing general-purpose FE programs can be easily modified for this new concept. Figure 7 shows the flow diagram of a typical FEM code extended for K-FEM. After the modification, the standard FEM becomes in fact a subclass of K-FEM. With this convenience, K-FEM has a high chance to be widely accepted in practice.

Conclusions
The basic concept and the advantages of K-FEM have been described. The present method is as simple as the conventional FEM in terms of its implementation; yet it retains much of the advantages of mesh-free methods.
K.Y. Dai et al. [13] pointed out that the method using standard Galerkin weak form with KI is noncomforming and so is K-FEM. This means the elemental piecewise KI is not fully compatible across the interelement boundaries. Its effect on the convergence was studied in the context of 2D elastostatic problems [14,15]. It was found that K-FEM with appropriate choice of correlation function passes the weak patch test and therefore the convergence can be guaranteed.
One possible drawback of K-FEM is its excessive demand of the computational time, as Kriging shape functions are constructed element by element during the computation. Moreover, a larger DOI means a longer time for stiffness formation and for solving a system with larger average bandwidth. However under the current trend, the cost of running a FEM project is heavily weighted on the engineer's time for preparing meshes, rather than on the computational time.
Several investigations have been carried out successfully on different applications of K-FEM. Aside from plane elasticity problems [9,16,17], so far Krigingbased finite elements have been developed for degenerated solid beams, plates and shells [12,18,19]. The results confirmed that K-FEM is indeed a viable alternative to the conventional FEM and has great potential in engineering applications. Future research may be directed at (1) applications of K-FEM to nonlinear problems and (2) improvement of its computational efficiency