Abstract
The method presented in this work is a 3-dimensional polyhedral finite element (3D PFEM) based on virtual node method. Novel virtual node polyhedral elements (termed as VPHE) are developed here, particularly virtual node hexahedral element (termed as VHE). Stiffness matrices of these polyhedral elements consist of simple polynomials. Thus, a new algorithm is introduced in this paper, which enables exact integration of monomials without a need for high number of integration points and weights. The number of nodes for VHE elements is not restricted, as opposed to the conventional hexahedral elements. This feature enables formulation of transition elements (termed as T-VHE) which are useful to adaptive computation. Performances of the new VHE elements in solid mechanics and conductive heat transfer phenomena are examined through numerical simulations. The new T-VHE elements are utilized in octree mesh. The VHE elements are found to produce good results and T-VHE elements help to reduce number of global nodes for the analysis.
1. Introduction
Tetrahedral and hexahedral elements are the most commonly utilized in three-dimensional finite element method (FEM), due to their simplicity and well establishment. Several comparison studies between these elements have been conducted. The comparison studies are briefly highlighted below.
Performances of tetrahedral and hexahedral elements in analysis of elastic and elastoplastic solid continuum are reported in [1, 2]. Simulation results indicate that linear hexahedral elements are superior to linear and even quadratic tetrahedral elements when it comes to analysis of problems in which the shear stress is dominant, such as in linear torsional problems. Comparison study in structural problems [3] showed that both elements performed equally in terms of accuracy and computational time. The authors in [4] recommended tetrahedral elements for applications in human body meshing for impact analysis and replace the commonly utilized hexahedral elements, due to meshing requirements for highly complex geometries of human organs. The authors showed that the tetrahedral elements can be superior to hexahedral elements, by utilizing reduced integration scheme. Another similar application is reported in [5], where the authors utilized tetrahedral and hexahedral elements to mesh human body part, which is the femur. Simulation results showed that both elements yield similar results for simulation of realistic proximal femur.
Application of tetrahedral and hexahedral elements in eddy current analysis [6] showed that both elements yield accurate results, but hexahedral elements require less computational time and memory storage compared to tetrahedral elements. In some cases, conventional tetrahedral elements are modified in order to yield more reliable and accurate results compared to hexahedral elements for certain applications (in order to facilitate meshing of complex geometries). Such example can be seen in forging simulations. In [7], the authors have reported that the specific element known as tetrahedral MINI element could perform as well as hexahedral elements despite having a touch of numerical uncertainty. The accuracy for the tetrahedral MINI element is attained with the help of intelligent remeshing technique. These three-dimensional elements (tetrahedral and hexahedral elements) have also been used in analysis of foot and footwear biomechanics [8]. The simulation results showed that hybrid hexahedral elements perform better compared to the hybrid tetrahedral elements. Some summaries based on the brief review above are given below.
Drawbacks of tetrahedral elements are the high computational cost (since large number of elements needed to mesh irregular geometries), high memory storage requirement, capability of being too stiff for certain application, and having high tendency to volumetric and shear locking. Main advantage of tetrahedral elements is that they are highly automated in mesh generation schemes. Therefore, tetrahedral elements are mostly used in mesh generation and adaptive mesh refinement schemes. Hexahedral elements on the other hand are efficient in terms of stability, solution accuracy, convergence rate, computational cost, and memory storage requirement. Main drawback of the hexahedral elements is that they require user intervention for mesh generation. Therefore, applications of hexahedral elements are limited due to the difficulties in meshing highly complex geometries. Initiatives have been taken to improve accuracy of tetrahedral elements to match that of hexahedral elements, due to the tetrahedron’s capability to generate automated mesh. At the same time, many new meshing algorithms have been proposed to generate hexahedral mesh with good quality [9].
Later, arbitrary 2-dimensional polygonal elements were developed [10–12], due to their advantages over the conventional elements, such as immunity to mesh distortion, flexibility during mesh generation, good candidate to represent polycrystalline materials, and being used as transition element during meshing. Tremendous amount of works in development of 2-dimensional polygonal finite element method (2D PFEM) can be seen in the literature, particularly by the author in [11]. Other developments on PFEM can be seen within smoothed finite element method (SFEM) [13] and virtual node method [14]. Subsequently, the 2D PFEM has been extended to 3D PFEM [15, 16]. 3D PFEM have been developed and implemented in various methods and fields, such as variable-element topology finite element method (VETFEM) [17], computer graphics [18–20], nonlinear continuum mechanics [21], electromagnetic simulations [22], boundary-element-based finite element methods (BEM-based FEM) [23, 24], flow simulations, material design, and fracture [25].
The 3D PFEM consist of nonpolynomial functions, which require special integration techniques. New methods have been proposed to perform the integration over polyhedral elements by utilizing quadrature points [26, 27]. However, numerical integration of nonpolynomials over polyhedral elements is still immature, since the polyhedral elements need to be partitioned into numerous numbers of subelements in which the numerical integration will be carried out or a lot of quadrature points and weights are needed to achieve good accuracy in integration over the polyhedral elements. The error in the numerical integration leads to nonconverging solutions [28, 29]. Thus, recently the 2D and 3D PFEM have been coupled with virtual element method (VEM) to avoid the errors in the numerical integration [30]. VEM is closely related to mimetic finite difference and many recent works on VEM coupled with PFEM can be seen in literature [31–37].
In this work, novel virtual node polyhedral elements (termed as VPHE) are developed, by utilizing virtual node [38]. Later, new 14-node virtual node hexahedral elements (termed as VHE) are developed and validated through numerical simulations using MATHEMATICA and ANSYS software. The new VHE elements are applied in two applications, which are the solid mechanics and conductive heat transfer phenomena. Geometries with linear and curved boundaries are selected for the analyses. Comparisons between the new VHE elements and the conventional elements are also provided for each case.
Stiffness matrices of the VPHE and VHE elements consist of simple polynomials (with variables , , and in terms of global coordinate system and , , and in terms of local coordinate system). Thus, a new technique is introduced in this paper, which enables exact integration of monomials of the form or . This new exact integration technique can be utilized to integrate the stiffness matrices of the VPHE and VHE elements within the element geometries, without a need for high number of integration points and weights. Finally, a more efficient meshing technique known as octree meshing is presented for the 14-node VHE elements. New conforming transition VHE elements (T-VHE) are developed to overcome the hanging nodes issue. The octree meshing technique helps to reduce number of global nodes for the analysis.
This paper is organized as follows. Formulation and properties of the new VPHE and VHE elements are described in Section 2. Finite element equations for solid mechanics and conductive heat transfer phenomena are provided in Section 3. New exact integration technique and octree meshing are presented in Section 4. Patch tests and application of the new elements in solid mechanics and conductive heat transfer phenomena are carried out in Section 5. Simulation results are also presented in this section. The paper is finally concluded in Section 6.
2. Element Formulation
This section describes formulation of the new VPHE and VHE elements by utilizing virtual node method [38]. This section covers the element geometry formulation, relevant formulas for the field variables and shape functions, and properties of the new elements.
Shape functions for a VPHE are formulated by coupling nodal shape functions of tetrahedral element with least-squares approximation. For that reason, the arbitrary VPHE is partitioned into number of tetrahedral elements. Each plane of the VPHE is partitioned into triangles which share a common centre, . The number of triangles on a face is equal to total number of nodes on the particular face. The fourth node for each triangle (vertex, ) which forms tetrahedral is located at the centre of the VPHE. An example of partitioning of an arbitrary VPHE is shown in Figure 1.
(a)
(b)
(c)
Linear tetrahedral element is utilized to partition the VPHE, since this element has linear sides and the shape function for a node is always independent from the other nodes, for any geometry configuration. The vertex, , acts as an imaginary or virtual node, which is utilized to combine all the tetrahedral elements to form the VPHE. The node is utilized only to develop shape functions for the VPHE nodes and does not contribute to computational cost.
The new 14-node VHE element consists of 14 nodes (8 corner nodes and 6 midnodes) and one virtual node (centre of the hexahedral element), as shown in Figure 2. The midnodes are placed at the centre of each face. The virtual node is utilized only to develop shape functions for the 14 nodes and does not contribute to any additional computational burden.
Each face of the hexahedral element is partitioned by utilizing 4 linear tetrahedral elements, as shown in Figure 3 for face 1-3-5-7. Similar partitioning is done for the other 5 faces and thus generates total of linear tetrahedral elements. These linear tetrahedral elements share a common node at the centre of the hexahedral element, which is the virtual node (node 15).
Field variable within each tetrahedral element within the polyhedron/hexahedron can be calculated by using the formula:where represents a specific tetrahedron within the polyhedron (i.e., ), , , , and represent shape functions of a specific tetrahedron for nodes , , , and , respectively, as described in Section 2.1, and , , , and represent field variables for the specific tetrahedral nodes , , , and , respectively. Nodes , , and are located on the outer surface of the polyhedron, while node is located at inside of the polyhedron. Node , which is located at inside of the polyhedron, is the common point shared by all the tetrahedrons (which form the polyhedron/hexahedron). Location of the common point (, , and ) is simply the centre of the polyhedron/hexahedron (node or 15 in Figures 1 and 3) which can be found usingwhere represents total number of nodes for the polyhedron/hexahedron, which is 14 for VHE element. Calculation of field variable for a polyhedron/hexahedron should involve the outer surface nodes (nodes , , and ), of all the tetrahedrons which form the polyhedron/hexahedron. However, it is seen that calculation of field variable by (1) introduces one additional variable, which is the field variable at node 4 .
The additional node (virtual node) can be eliminated from (1) by introducing the term , which represents field variable within the polyhedron/hexahedron. The field variable for node () in (1) can then be replaced by this new term. Calculation of the field variable within the polyhedron/hexahedron, , involves shape functions and field variables of all the external nodes of the polyhedron/hexahedron, as shown in Section 2.2. This eliminates the need for definition of the field variable at the virtual node. The other three nodal field variables , , and can be represented by a common term (in order to combine all the tetrahedral elements to form the polyhedron/hexahedron). Equation (1) can then be rewritten aswhere represents field variable within a tetrahedron, as described in Section 2.1. Equation (3) can then be extended asVirtual node in (4) () introduces additional degree of freedom and it can be eliminated by utilizing (10):Finally, (5) can be rewritten to get expressions for shape functions of each polyhedron/hexahedron node with respect to the tetrahedron as number of shape functions ( is total number of nodes of the polyhedron/hexahedron, which is 14 for VHE element) would be generated by using (6), for each tetrahedron. Thus, there would be total of number of shape functions, where represents total number of tetrahedrons within the polyhedron/hexahedron. Each set of shape functions which are generated based on a particular tetrahedron would then be integrated within the respective tetrahedron.
2.1. Shape Functions (, , , and ) and Field Variable () Corresponding to a Tetrahedron
Shape functions for a tetrahedron are expressed in terms of the natural coordinate system, in order to ease integration of the stiffness matrices, as described in Section 2.3(d). The shape functions areThus, each tetrahedron within the polyhedron/hexahedron (with coordinates ()) is mapped to the parent element geometry in the natural coordinate system (with coordinates ()) according to the equation below:An example is shown in Figure 4, for a specific tetrahedron 7-3-12-15 within the hexahedron (Figure 3).
Field variable within a tetrahedron with respect to the tetrahedron, , is found by
2.2. Field Variable within the Polyhedron/Hexahedron,
Field variable within the polyhedron/hexahedron is found by utilizing the least-square method [14, 39, 40], given bywhere is the matrix of nodal shape functions, is the vector of field variables, and is a set of basis functions with number of polynomials (selected based on the Pascal Triangle). Four terms are selected for the hexahedron : and are basis matrices given byConstrained mean least-square (CO-MLS) is later utilized in order to prevent calculation of inverse matrix in (10) [41]. The calculation of inverse matrix causes high computational time and does not yield solution for certain matrices.
2.3. Properties of the VPHE and VHE Shape Functions
Shape functions for the VPHE and VHE elements exhibit the following properties.
(a) The field variables are continuous within the domain and linear on the boundaries. This is because each shape function is developed based on the linear tetrahedral elements that form the hexahedron. Each side of these arbitrary tetrahedrons (6 sides, -, -, -, -, -, and -) is mapped to their respective sides on the parent element (Figure 4), and thus the junctions of adjacent tetrahedrons share similar linear functions as well (continuous at the junctions). Contour plot of the geometry mapping for an arbitrary hexahedron based on (13) is shown in Figure 5. It can be seen that the interpolation functions yield linear mapping. - and -coordinates are dependent on all the three variables , , and , while -coordinates are independent of : where is hexahedral shape functions corresponding to tetrahedron .
(a)
(b)
(c)
(b) The shape functions fulfil Kronecker delta property:Vector plots for three shape functions of an arbitrary hexahedron corresponding to the base of the tetrahedron are shown in Figure 6. The three shape functions correspond to nodes 1, 2, and 10, respectively (, , and ).
(a)
(b)
(c)
From Figure 6, it can be seen that converges towards the node (, , and ) in the parent element (Figure 4(b)), converges towards the node ( and ) in the parent element, and finally converges towards the node ( and ) in the parent element.
(c) The shape functions fulfil linear completeness:where , , , or , , , and represent nodal coordinates at the parent element corresponding to node (Figure 4(b)), and represents -coordinates of the element in global coordinate system corresponding to tetrahedron (Figure 4(a)).
(d) The shape functions consist of simple polynomials of the form . Therefore, it simplifies numerical integration of the stiffness matrices. Conventional Symmetric Quadrature rules for the parent tetrahedral element (Figure 4(b)) can be used to perform the numerical integrations. Generation of high number of integration points and weights for higher order elements can be prevented by utilizing the new exact integration technique presented in Section 4.
3. Finite Element Equations
The finite element equations for solid mechanics are given aswhere represents stiffness matrix corresponding to tetrahedron with dimensions (3 degrees of freedom per node), represents partial derivatives of the interpolation functions with dimensions () corresponding to tetrahedron , represents elastic material property matrix, represents modulus of elasticity, and represents Poisson’s ratio. Once stiffness matrices corresponding to each of the tetrahedrons within a polyhedron/hexahedron () are found, the stiffness matrix for a single polyhedron/hexahedron, , is then calculated by merging the stiffness matrices: Displacements and/or forces can then be computed for the global system using the following equations:where represents global stiffness matrix, represents global displacement matrix, and represents global force matrix. Postprocessing can be done by using the formulas:where represents displacement matrix for a polyhedron/hexahedron, represents strain components for a polyhedron/hexahedron corresponding to tetrahedron , and represents stress components for a polyhedron/hexahedron corresponding to tetrahedron .
The finite element equations for conduction heat transfer arewhere is the global stiffness matrix, is the global nodal temperatures matrix, is the global forcing matrix, is the thermal conductivity matrix, and , , are the thermal conductivities in -, -, and -directions, respectively. Stiffness matrices for individual elements, , are found using similar formula in (16).
4. A New Exact Integration Technique and Octree Meshing
Integration of the stiffness matrices (16) of the 14-node VHE elements within the element geometries can be done by utilizing Conventional Symmetric Quadrature rules for the parent tetrahedral element (Figure 4(b)). However, higher number of integration points and weights would be required to integrate higher order polynomials. A new technique is introduced here, which enables exact integration of monomials of any order, without requiring any integration points or weights. Later, an octree meshing technique is described with the utilization of transition elements. Two examples of transition elements are shown for the new VHE elements.
4.1. Generalized Equations for Exact Integration
Integration limits of a function to be integrated need to be converted as required by quadrature rule of choice with range . General equation to convert integration limits of a single integration can be represented numerically by [42]where is upper limit, is lower limit, , and .
It can be seen that the conversion of the integration limits above is done by fully numerical techniques (symbolic manipulations are avoided), in order to reduce computational time. Equation (21) can then be repeated three times, in order to get the generalized equations for 3-dimensional domains:Equation (22) is an example for integration order . Sequences of the differentiation in (22) can be altered accordingly for other integration orders. General equations for exact integration of trivariate monomials over a three-dimensional domain can be derived when the integration limits , , , , , and in (22) are . Thus, the generalized equation (22) can be utilized to convert the integration limits by letting and then the exact integration of trivariate monomials can be represented numerically using the relationIt can be seen that the method shown in (23) requires handling of symbols to identify the variables , , and their respective powers, , , and .
Effectiveness of the exact integration technique above is now tested and verified. Two different polynomial functions are integrated over a domain which is enclosed by a linear plane and a curved plane (parabolic curve). The results are later compared with the analytical solution. Percentage error is calculated based onThe simulations are run on a computer with 2.93 GHz Dual Core CPU, 32-bit operating system, and 2 GB of memory. Results are shown in Table 1.
From the simulation results, it can be seen that the exact integration technique yields accurate solutions at lower computational time compared to the analytical solutions. The new exact integration technique is suitable for FEM based on stiffness matrices which consist of simple polynomials, such as virtual node method and strain based elements. This technique provides an alternative for existing numerical integration techniques, especially Gaussian quadrature which requires high number of integration points and weights to integrate higher order trivariate monomials.
4.2. Octree Meshing Technique
Fine regular meshes lead to more accurate solutions. However, computational cost will be increased when the entire domain is analysed by using fine mesh. One of the methods to overcome this problem is by utilizing the octree meshing technique, which utilizes cubic elements for adaptive mesh generation. Octree mesh enables fine regular cubic elements to be generated at the regions where the rate of change of field variable is the highest. Minimum number of cubic elements can be used to mesh regions in which the rate of change of field variable is the smallest. These two different regions can then be coupled by utilization of transition elements. The transition elements need to be compatible, in order to avoid hanging nodes. There have been many works done to generate compatible transition elements. There are two ways to generate the transition elements. One of the ways is to split the hexahedral element into tetrahedral or pyramid elements [43]. However, this technique requires additional effort to find the elements with hanging nodes and partition them into tetrahedral or pyramid elements and, hence, increases computational cost. The other method is to develop specific transition elements [44] but usually limited to 1-irregular mesh. One-irregular mesh means that a single transition element is coupled with 4 regular elements at one of the faces of the transition element.
In this work, new transition VHE (termed as T-VHE) elements are developed to be utilized in octree meshing technique. These new T-VHE elements are compatible and general, since they are not limited to 1-irregular mesh [44]. Two specific conforming transition elements are developed as examples, which are the 17-node and 22-node T-VHE elements. The first transition element is used to couple two regular cubic elements onto one of the faces of the transition element, while the latter is used in 1-irregular mesh.
Geometries of the transition elements are as shown in Figures 7(a) and 7(b). Partitioning is shown only for the transition faces in these figures.
(a)
(b)
Each face of the T-VHE element is partitioned into 4 tetrahedral elements (similar to the 14-node VHE elements), except for face 1-3-5-7. The face 1-3-5-7 would be used to couple with the regular hexahedral elements from different meshing region (fine meshing region). As can be seen in Figures 7(a) and 7(b), the face 1-3-5-7 of the 17-node T-VHE element is partitioned using 8 tetrahedral elements, while similar face of the 22-node T-VHE element is partitioned using 16 tetrahedral elements. Therefore, these T-VHE elements can be used to couple with two and four regular 14-node VHE elements from fine meshing regions, respectively. Simulation result for the transition element is shown in Section 5.2.1.
5. Patch Tests and Numerical Simulations
The new VPHE, 14-node VHE, and the T-VHE elements are validated through patch tests. Performances of the 14-node VHE and T-VHE elements are later examined through numerical simulations in solid mechanics and heat transfer phenomena. Material properties used in the patch tests as well as in all the numerical simulations are modulus of elasticity, = 106 Pa; Poisson ratio, = 0.3; and thermal conductivities, = = = 25 W/m°C.
5.1. Patch Tests
The new VPHE elements are tested in patch consisting of 6 VPHE elements with combination of 6, 7, and 8 faces as shown in Figure 8. Four different patches (with same domain) are used for the new 14-node VHE elements, as shown in Figure 9.
(a)
(b)
(c)
(d)
The external nodes for each case were assigned with fixed displacement/temperature values according to the linear functions:The linear displacements/temperatures on the external nodes represent constant strain/heat flux within the interior of the elements. The displacements/temperatures of the inner nodes were calculated using finite element equations given in Section 3 and the difference between the calculated and exact values was found using error norm in displacement/temperature usingThe patch test results for the VPHE elements are shown in Table 2, while the patch test results for 14-node VHE elements are shown in Table 3.
It is seen that the range for the error norms of the VPHE elements is similar to those obtained using Wachspress and mean value coordinate for 2 dimensions [45] and 3D PFEM [30] for 3 dimensions, but the integration of the stiffness matrices is much simpler. It is also observed that error norm for VPHE in general polyhedral form is higher compared to VPHE in the form of hexahedrons. It is also observed that the error norm in heat transfer phenomena is lower. These are caused by higher degree of freedom and nodal distortions in general VPHE.
5.2. Applications in Solid Mechanics and Conductive Heat Transfer Phenomena
A beam with linear boundaries and a hollow cylinder with curved boundaries are selected to carry out numerical simulations. Geometries are shown in Figures 10(a) and 10(b). Dimensions of the beam and hollow cylinder are length = 5 m, width = height = 1 m, inner radius = 1 m, outer radius = 4 m, and thickness of the cylinder = 1 m. Due to the symmetry, only quarter of the hollow cylinder is analysed, as shown in Figure 10(c).
(a)
(b)
(c)
Total of four cases are examined (for solid mechanics and conductive heat transfer phenomena). Material properties are set similar to what is described in Section 5. Method, geometry, and boundary conditions for each case are given below.
Case 1 (solid mechanics; a beam under bending load).
Method. A beam is subjected to bending loads and the maximum tip displacement due to the applied loads is determined.
Geometry. See Figure 10(a).
Boundary Conditions. Four concentrated vertical loads, N, are applied onto nodes , , , and . Surface of the beam is fixed and constrained from moving in all three directions.
Case 2 (solid mechanics; a hollow cylinder subjected to internal pressure).
Method. A hollow cylinder is subjected to uniform internal pressure. Displacements and stresses at inner and outer nodes are determined.
Geometry. See Figure 10(c).
Boundary Conditions. The internal surface of the cylinder 1-2-5-6 is subjected to pressure of 100 Pascals. The surface 1-3-5-7 is constrained in x-direction, surface 2-4-8-6 is constrained in y-direction, and surface 5-6-7-8 is constrained in z-direction.
Case 3 (conductive heat transfer; a hollow cylinder under conductive heat transfer).
Method. A hollow cylinder is subjected to fixed thermal loads at the inner and outer surfaces. Steady-state temperatures within the hollow cylinder are then calculated numerically.
Geometry. See Figure 10(c).
Boundary Conditions. The temperature of the inner surface 1-2-5-6 is fixed at 500°C and the outer surface 3-4-7-8 is fixed at 100°C. Other surfaces are considered adiabatic.
Case 4 (Cases 1 and 3 with octree mesh). Cases and are repeated here but by utilizing octree meshing technique.
Four types of partitioning are considered, which are classified as Mesh 1, Mesh 2, Mesh 3, and Mesh 4. Mesh 1 is generated by utilizing the 14-node VHE element. The beam is partitioned by using undistorted elements (cube) and the hollow cylinder is partitioned by using distorted elements. Three levels of partitioning are considered for each geometry, which are 1 × 1 × 5 mesh, 2 × 2 × 10 mesh, and 3 × 3 × 15 mesh for the beam (Figure 11) and 2 × 2 × 1 mesh, 3 × 3 × 1 mesh, and 4 × 4 × 1 mesh for the hollow cylinder (Figure 12).
(a)
(b)
(c)
(a)
(b)
(c)
Mesh 2 is generated by utilizing conventional tetrahedral elements. For this case, all the 24 tetrahedrons that form the 14-node VHE element are considered as individual linear tetrahedral elements. Thus, each 14-node hexahedral element in Mesh 1 is converted to 24 independent linear tetrahedral elements, as shown in Figures 13 and 14.
(a)
(b)
(c)
(a)
(b)
(c)
Mesh 3 is classified as random mesh generation using ANSYS software. Mesh 3 is also classified into three levels, in which each level contains similar/close number of global nodes used for , , , and cases in previous Mesh 1 and Mesh 2.
Mesh 4 is classified as octree mesh generation. The two types of octree mesh generation utilized in the simulations (Case 4) are shown in Figures 15 and 16.
Numerical simulations for Mesh 1, Mesh 2, and Mesh 4 are performed by utilizing MATHEMATICA software, while solutions for Mesh 3 are obtained from ANSYS software.
5.2.1. Simulation Results
Simulation results and analysis for the 4 cases above are presented in this section.
Case 1 (solid mechanics; a beam under bending load). The maximum tip displacement due to the applied loads is determined numerically by running simulations in MATHEMATICA (for Mesh 1 and Mesh 2) and ANSYS (for Mesh 3) software and the results are compared with the reference solution. Solution which is obtained by using high number of elements in ANSYS software is taken as the reference solution. The reference solution for this case is obtained by using 1717 global nodes and 5563 linear tetrahedral elements. Simulation results are shown in Figure 17.
It can be seen that the magnitude of maximum tip deflection converges towards the reference solution, when number of global nodes are increased, except for Mesh 2. It can also be seen that the convergence of the new 14-node VHE elements (Mesh 1) is faster compared to the conventional linear tetrahedral elements (Mesh 3).
Case 2 (solid mechanics; a hollow cylinder subjected to internal pressure). Stresses and displacements at nodes 1 to 8 as shown in Figure 10(c) are determined by performing numerical simulations in MATHEMATICA and ANSYS software. Five types of elements were used: linear tetrahedral element (Linear Tetra), quadratic tetrahedral element (Quad Tetra), linear hexahedral element (Linear Hexa), quadratic hexahedral element (Quad Hexa), and the new 14-node VHE element (VHE). Mesh 3 was used for all the conventional elements and Mesh 1 is used for the VHE elements. The values obtained from the numerical simulations are then compared with analytical solutions. The analytical solutions are given by the formulas below [46]:where represents displacement, represents inner radius, represents outer radius, represents arbitrary radius, represents the inner pressure, represents modulus of elasticity, and represents Poisson’s ratio. The conventional elements provide accurate solutions for the displacements but failed to provide correct value for surface stresses, as shown in Figure 18.
Another simulation was run by utilizing Mesh 2 (Linear Tetra) and the results are as shown in Figures 19 and 20.
The new 14-node elements are able to provide converging solutions for inner surface stress, as shown in Figure 18. From the results of Figure 19 for displacements, it is seen that the new elements perform better compared to the conventional linear tetrahedral elements. The 14-node VHE elements showed symmetrical results for the displacements at inner and outer surfaces of the hollow cylinder. The conventional linear tetrahedral elements failed to converge towards the analytical solutions. The VHE elements showed similar performances for the stresses (Figure 20), except for the stresses at the outer surface of the cylinder. This is because the FEM minimizes the weighted-integral of the error in the differential equation over each element, by adjusting the nodal values. The FEM does not minimize the error in the nodal values [47]. However, the VHE element showed convergence towards the analytical solution.
(a)
(b)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Case 3 (conductive heat transfer; a hollow cylinder under conductive heat transfer). Steady-state temperatures within the hollow cylinder are determined numerically by running simulations in MATHEMATICA (for Mesh 1) and ANSYS (for Mesh 3) software. Five types of elements were used, similar to Case 2. The results are later compared with the analytical solution which is given by [48]where represents arbitrary radial distance, is temperature at the inner surface, is temperature at the outer surface, is radius of the inner surface, and is radius of the outer surface. Percentage error is found using (24).
Simulation results are shown in Figure 21.
From the results, it can be seen that the 14-node VHE elements yield smaller percentage of error compared to the conventional elements.
Case 4 (Cases 1 and 3 with octree mesh). Two of the previous cases, which are Cases 1 and 3, are repeated here but by using the octree mesh as shown in Figures 15 and 16. The 22-node T-VHE element is utilized to generate the octree mesh (element 3). For the beam, fine meshing is used at the end which is subjected to load variation. The octree mesh consists of total of 19 elements with 140 global nodes. The maximum tip deflection obtained is 0.1374 m. From the results in Figure 17, about 300 global nodes are needed to achieve similar accuracy. For the hollow cylinder, fine meshing is used near the surface subjected to the temperature load. The octree mesh consists of total of 30 elements with 205 global nodes. The temperature at radius of 1.5 m is obtained as 383.769°C, with error of 0.2%. It is seen that the octree meshing technique can be used to obtain good results at the region of interest, without the need to mesh entire regions by using fine mesh.
6. Conclusions
Novel VPHE and 14-node VHE elements have been successfully developed and validated. The elements were successfully developed by utilizing 3-dimensional polyhedral finite element (3D PFEM) based on virtual node method. These elements are found to produce good results. Shape functions for these elements are simple polynomials; therefore, exact integration techniques can be utilized. The new exact integration technique presented in this work yields accurate solutions at lower computational time compared to the analytical solutions. This technique provides an alternative for existing numerical integration techniques, especially Gaussian quadrature which requires high number of integration points and weights to integrate higher order trivariate monomials. More efficient meshing technique, which is octree mesh, is introduced for the new 14-node VHE element. New transition elements (T-VHE) are developed and applied in octree mesh. Simulation results showed that the octree mesh helps to reduce number of global nodes.
Competing Interests
The author declares that there are no competing interests related to this paper.
Acknowledgments
The author would like to thank Research Management Centre (RMC) of Multimedia University, Malaysia, for providing financial support through Mini Funds with Grant nos. MMUI/130070 and MMUI/160047, which enabled purchase of required software and equipment for this work. The author would also like to express sincere appreciation to Dr. Xu Hai Tang and Dr. Yongtao Yang, for the valuable discussions.