# PSBFEM-Abaqus: Development of User Element Subroutine (UEL) for Polygonal Scaled Boundary Finite Element Method in Abaqus

Authors: Nan Ye; Chao Su; Yang Yang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

The polygonal scaled boundary finite element method (PSBFEM) is a novel method integrating the standard scaled boundary finite element method (SBFEM) and the polygonal mesh technique. This work discusses developing a PSBFEM framework within the commercial finite element software Abaqus. The PSBFEM is implemented by the User Element Subroutine (UEL) feature of the software. The details on the main procedures to interact with Abaqus, defining the UEL element, and solving the stiffness matrix by the eigenvalue decomposition are present. Moreover, we also develop the preprocessing module and the postprocessing module using the Python script to generate meshes automatically and visualize results. Several benchmark problems from two-dimensional linear elastostatics are solved to validate the proposed implementation. The results show that PSBFEM-UEL has significantly better than FEM convergence and accuracy rate with mesh refinement. The implementation of PSBFEM-UEL can conveniently use arbitrary polygon elements by the polygon/quadtree discretizations in the Abaqus. The developed UEL and the associated input files can be downloaded fromhttps://github.com/hhupde/PSBFEM-Abaqus.


---

## Body

## 1. Introduction

The finite element method (FEM) is a reliable computational tool to solve partial differential equations (PDE) in science and engineering [1–4]. A domain of complex geometry is partitioned into a finite number of nonoverlapping subdomains of simplex shapes by introducing the concept of discretization [5]. Typically, the shape of the conventional two-dimensional finite element method is triangles or quadrilaterals. At present, the conventional FEM also faces several problems. For example, (a) the accuracy of the solution depends on the quality of the mesh and (b) requires sophisticated discretization techniques to generate high-quality meshes and to capture topological changes [6]. These elements used in conventional FEM must conform to the domain’s boundary, which leads to difficulty in solving many complex problems.To overcome the weakness above the conventional FEM, researchers proposed other alternatives methods, such as the meshfree method [7–9], the smoothed finite element method [6, 10], the Isogeometric Analysis (IGA) [11], Deep Neural Networks (DNNs) [12], and the polygonal finite element method [10, 13]. Polygonal element with more than four sides involves more nodes in their interpolation compared with a conventional FEM. Simultaneously, they generally exhibit superior solution accuracy [5]. Moreover, it is more flexible in the discretization of complex geometry. Therefore, these advantages have further motivated polygonal elements as an alternative to conventional FEM using triangles or quadrilaterals.The scaled boundary finite element method (SBFEM) is an alternative method to construct polygonal elements. Song and Wolf developed the technique in the 1990s [14]. The scaled boundary finite element method is a semianalytical method that attempts to fuse the advantages and characteristics of FEM and the boundary element method (BEM) into one new approach. The SBFEM has been applied to many physical field problems, such as wave propagation [15, 16], heat conduction [17, 18], fracture [19, 20], acoustic [21], seepage [22, 23], elastoplastic [5], and fluid [24]. For these problems, the SBFEM presents more efficiency compared with the conventional FEM.The polygonal scaled boundary finite element method (PSBFEM) is a novel method integrating the standard SBFEM and the polygonal mesh technique. This method is flexible in meshing complex geometries, and the use of polygons to discretize the computational domain naturally complements the SBFEM. Recently, an alternative mesh technique has been widely used in geometric discretization. In computational mechanics, the quadtree algorithm is usually used in large-scale simulations typical in the modeling of earthquake and ground motions [25], flood [26], and tsunami [27]. The quadtree algorithm is fast, efficient, and capable of achieving rapid and smooth transitions of element sizes between mesh refinement regions [28]. Mesh generation and adaptive refinement of quadtree meshes are straightforward [29]. Due to hanging nodes between two adjacent elements of different sizes, it is problematic that quadtree meshes are directly used to simulate within the finite element method’s framework. However, the SBFEM only discretizes in the boundary of geometry. Hence, each cell in a quadtree mesh, regardless of hanging nodes, is treated as a generic polygon. This enables the structure of the quadtree to be exploited for efficient computations. The ability to assume any number of sides also enables the SBFEM to discretize curved boundaries better.Although the PSBFEM has become quite mature, these studies only exist as a few stand-alone codes. The PSBFEM has not yet formed a part of commercial software at present. Hence, it is not easy to use the method in a specific community or laboratory. The commercial software Abaqus has powerful linear or nonlinear, static, or dynamic analysis capabilities. Also, Abaqus/Standard analysis provides a User Element Subroutine (UEL) to define an element with a very available option to interface with the code. Liang et al. [30] developed a UEL for dynamic analysis of saturated porous media based FEM. Kumbhar et al. [6] implemented the element based smoothed finite element method (CSFEM) by the UEL subroutines. Molnar et al. [31] implemented an implicit, staggered elastoplastic version of the phase-field approach Abaqus through the UEL. Therefore, the UEL can be used to extend the Abaqus to solve new problems conveniently.Recently, Ya et al. [32] implemented an open-source polyhedral SBFEM element for three-dimensional and nonlinear problems in the commercial software Abaqus through the UEL. The code can use the polyhedral element to enhance the performances of ABAQUS for interfacial problems, and it significantly reduces the meshing burden encountered in the FEM. At the same time, it also implemented the technique of octree mesh generation, which is efficient and robust to mesh complex geometries and promise to integrate geometric models and numerical analysis in a fully automatic manner.In engineering designs, to calculate models to be used as design tools, two-dimensional (2D) models are relatively easy to set up and have reasonably short computational times, which would allow sensitivity and optimization analyses [33]. There are no available subroutines of SBFEM to solve 2D problems in the commercial software Abaqus at present. However, the implementation of two-dimension can be easily derived from three-dimensional implementation. Hence, Ya et al. [32] provided a reference for developing the UEL to solve 2D problems using the polygonal scaled boundary finite element method with the polygon/quadtree meshes.In addition, to solve the stiffness matrix, we need to perform the eigenvalue decomposition in the SBFEM. Ya et al. [32] used the Linear Algebra Package (LAPACK) [34] to perform the eigenvalue decomposition by adding the LAPACK source codes into UEL subroutines. Hence, LAPACK needs to be compiled every time. It will take extra time when calling UEL subroutines. Abaqus provides mathematical libraries (the Intel Math Kernel Library, MKL [35]) to perform the eigenvalue decomposition. The user can directly use MKL only by modifying the Abaqus environment file, which can avoid taking extra time when calling UEL subroutines.In this paper, we implement the PSBFEM technique with Abaqus and provide the details. For simplicity, we consider a two-dimensional linear elastostatics problem using the PSBFEM technique. This work is organized as follows: the basic principles of the SBFEM are introduced in Section2. Section 3 describes the implementation of PSBFEM by the Abaqus UEL subroutine. Moreover, we also develop the preprocessing and postprocessing module using the Python script. Section 4 presents a detailed convergence study by comparing the convergence and accuracy rates with the conventional FEM using several numerical benchmarks in two-dimensional linear elastostatics. Finally, the concluding remarks of this work are summarized in Section 5.


## 2. Theory

### 2.1. The SBFEM Equations in the 2D Linear Elastostatics

The fundamental difference between SBFEM and FEM is introducing a scaling centerO [36]. As shown in Figure 1, a scaling center can be directly visible from any point on the whole boundary of the S-element. By a scaling center, we can transform the scaled boundary coordinates to polar coordinates. Hence, the S-element can be described by a radial coordinate ξ and a circumferential coordinate η.Figure 1
Concept of scaled boundary finite element method.In this section, we consider a two-dimensional isotropic linear elastostatics problem. The governing equation and the boundary conditions can be expressed as(1)∇⋅σ+b=0,inΩ,u=u¯,onΓu,σ⋅n=t¯,onΓt,where ∇ is the differential operator, σ is the Cauchy stress tensor, b is the body force, Ω is the computational domain. Γu is the displacement boundary conditions, and Γt is the surface traction boundary conditions. The displacement u¯ and the surface traction t¯ are imposed on the domain boundaries Γu and Γt with the outward unit normal n.As illustrated in Figure2, the SBFEM presents a local coordinate system ξ,η. The coordinates of a point x,y along the radial line and inside the domain can be expressed as follows [36]:(2a)x=ξNηx,(2b)y=ξNηy,where ξ,η are the scaled boundary coordinates in two dimensions, ξ is a radial coordinate, and η is the circumferential coordinate.Figure 2
The SBFEM local coordinates.The differential operator can be transformed from the Cartesian coordinate systemx,y into the scaled boundary coordinates system ξ,η as follows:(3)∇=b1∂∂ξ+1ξb2∂∂η,with(4a)b1=1Jbyb,η00−xb,η−xb,ηyb,η,(4b)b2=1Jb−yb00xbxb−yb,where the Jacobian matrix at the boundary ξ=1 can be expressed as(5)Jb=xbybxb,ηyb,η=xbyb,η−ybxb,η.A comma followed by a subscript is used to denote partial differentiation to the variable in the subscript. The displacement fielduξ,η at any point in SBFEM coordinates is written as(6)uξ,η=Nuηuξ,where uξ is radial displacement functions along a line connecting the scaling center O and a node at the boundary. Nuη is the shape function matrix(7)Nuη=N1η0N2η0…0Nmη00N1η0N2η…0Nmη.The strain-displacement transition matricesB1=B1η and B2=B2η are introduced:(8a)B1=b1Nu,(8b)B2=b2Nu.η.The strain fieldε is expressed in the scaled boundary coordinates as(9)ε=B1uξ,ξ+1ξB2uξ.The stress field can be express as(10)σ=DB1uξ,ξ+1ξB2uξ,where D is an elasticity matrix and has the form(11)D=E1−v21v0v10001−v2,and for plane stress cases, E is Young’s modulus, and v is Poisson’s ratio.According to the virtual work principle, the radial displacements functionuξ is the solution of the SBFEM equation in displacement [36]:(12)E0ξ2uξ,ξξ+E0+E1T−E1ξuξ,ξ−E2uξ+Fξ=0,where uξ,ξξ is the second partial derivative with respect to the variable ξ, uξ,ξ is the first partial derivative with respect to variable ξ, Ei, i=0,1,2, are the coefficient matrices, and Fξ is a load vector.In the derivation, the governing equations of linear elasticity are weakened in the circumferential direction, while the strong form remains in the radial direction. The coefficient matricesEi, i=0,1,2, depends only on the geometry and material properties. These coefficient matrices can be given as follows:(13a)E0=∫−1+1B1TDB1Jbdη,(13b)E1=∫−1+1B2TDB1Jbdη,(13c)E2=∫−1+1B2TDB2Jbdη.WhenFξ=0, the solution of equation (12) is second-order ordinary differential equations that can be obtained by introducing the variable:(14)Xξ=uξqξ,where qξ is the internal force vector. Equation (12) can be transformed into a first-order ordinary differential equation with twice the number of an unknown as(15)ξXξ,ξ−ZpXξ=0,where the coefficient matrix Zp is a Hamiltonian matrix. The solution for a bounded domain is obtained using the positive eigenvalues of Zp. Hence, Zp can be expressed as(16)Zp=−E0−1E1TE0−1E2−E1E0−1E1TE1E0−1.The solution of the equation (15) can be obtained by the mathematical theorem and eigenvalue decomposition technique. The eigenvalue decomposition of Zp is expressed as(17)ZpΦunΦupΦqnΦqp=ΦunΦupΦqnΦqpλn00λp.The real parts of eigenvaluesλn are negative, and that of λp are positive. Φun and Φup are the transformation matrices corresponding to the modal displacements and forces, respectively. The general solution of equation (15) can be written as(18)Xξ=ΦunΦupΦqnΦqpξ−λn00ξ−λpcncp,where cn and cp are the integration constants. For a super element with 0≤ξ≤1, the solution of equation (15) is(19)uξ=Φunξ−λncn,(20)qξ=Φqnξ−λncn,where the integration constants cn can be extracted from the nodal displacements on the boundary ub=uξ=1 as(21)cn=Φun−1ub.Eliminating the integration constantscn at ξ=1 is given(22)F=qξ=1=ΦqnΦun−1ub.Hence, the stiffness matrix of the S-element can be written as(23)K=ΦqnΦun−1.The displacement fielduξ,η inside subdomains by a line element on the S-element can be given (24)uξ,η=NuηΦunξ−λncn.


### 2.2. The Element of PSBFEM

The PSBFEM is inherently appropriate for modeling polygons and has other promising capabilities. Because the PSBFEM is discretized only in the boundary, and an element of PSBFEM can assume more complex shapes than a finite element method, the PSBFEM element is much more flexible than the FEM element. With such an important advantage, the PSBFEM is more applicable to complex geometries than other methods. Figure3 shows the supported element type in the PSBFEM-UEL. In the PSBFEM-UEL, we can use the Abaqus element type, such as CPS3, CPE3, CPS4, and CPE4. Besides, the PSBFEM-UEL also provides the complex element type: polygonal elements and complex quadrilateral element (quadtree discretization). Hence, it is an effective tool to solve complex elements in numerical analysis.Figure 3
The supported element type in the PSBFEM-UEL; (a) a polygon discretization illustrating the mesh generation using PSBFEM-UEL element; (b) a quadtree discretization illustrating the mesh generation using PSBFEM-UEL element; (c) the PSBFEM-UEL element library.
(a)(b)(c)


### 2.3. Automatic Meshes Generation

This section uses a Python script to automatically generate PSBFEM elements by the Delaunay triangulation [37]. This algorithm is robust and efficient, and it has been integrated into Abaqus CAE. The algorithm can be used to triangulate any set of points on a two-dimensional plane. The triangulated mesh is then used to generate the polygon elements. The detail of the polygon generation algorithm mainly contains two steps. Firstly, we automatically generate triangular mesh by the Abaqus CAE. Secondly, considering each triangle interior node as the center of the polygon element, a polygonal element can be generated by connecting the centroids of all the triangular elements circumventing. More detail of this algorithm is presented in Figure 4. Finally, the PSBFEM-UEL can be directly used with a polygon mesh generator to analyze complex geometry problems.Figure 4
Generation of a polygon element from a triangular element: (a) arbitrary geometric model; (b) generating triangular meshes by Abaqus CAE; (c) construction of polygon meshes; and (d) final SBFEM polygon meshes.The quadtree decomposition is a tree data structure in which each parent has precisely four children [38]. The quadtree meshes are fast, efficient, and capable of achieving rapid and smooth transitions of element sizes between mesh refinement regions. The quadtree mesh can provide vital support in the preprocessing for the adaptive analysis of the SBFEM. We developed a quadtree mesh automatic generation code by a Python script. A simple example of quadtree discretization with three levels is illustrated in Figure 5. The quadtree algorithm is described in more detail in [28, 38].Figure 5
A quadtree discretization with three levels.


## 2.1. The SBFEM Equations in the 2D Linear Elastostatics

The fundamental difference between SBFEM and FEM is introducing a scaling centerO [36]. As shown in Figure 1, a scaling center can be directly visible from any point on the whole boundary of the S-element. By a scaling center, we can transform the scaled boundary coordinates to polar coordinates. Hence, the S-element can be described by a radial coordinate ξ and a circumferential coordinate η.Figure 1
Concept of scaled boundary finite element method.In this section, we consider a two-dimensional isotropic linear elastostatics problem. The governing equation and the boundary conditions can be expressed as(1)∇⋅σ+b=0,inΩ,u=u¯,onΓu,σ⋅n=t¯,onΓt,where ∇ is the differential operator, σ is the Cauchy stress tensor, b is the body force, Ω is the computational domain. Γu is the displacement boundary conditions, and Γt is the surface traction boundary conditions. The displacement u¯ and the surface traction t¯ are imposed on the domain boundaries Γu and Γt with the outward unit normal n.As illustrated in Figure2, the SBFEM presents a local coordinate system ξ,η. The coordinates of a point x,y along the radial line and inside the domain can be expressed as follows [36]:(2a)x=ξNηx,(2b)y=ξNηy,where ξ,η are the scaled boundary coordinates in two dimensions, ξ is a radial coordinate, and η is the circumferential coordinate.Figure 2
The SBFEM local coordinates.The differential operator can be transformed from the Cartesian coordinate systemx,y into the scaled boundary coordinates system ξ,η as follows:(3)∇=b1∂∂ξ+1ξb2∂∂η,with(4a)b1=1Jbyb,η00−xb,η−xb,ηyb,η,(4b)b2=1Jb−yb00xbxb−yb,where the Jacobian matrix at the boundary ξ=1 can be expressed as(5)Jb=xbybxb,ηyb,η=xbyb,η−ybxb,η.A comma followed by a subscript is used to denote partial differentiation to the variable in the subscript. The displacement fielduξ,η at any point in SBFEM coordinates is written as(6)uξ,η=Nuηuξ,where uξ is radial displacement functions along a line connecting the scaling center O and a node at the boundary. Nuη is the shape function matrix(7)Nuη=N1η0N2η0…0Nmη00N1η0N2η…0Nmη.The strain-displacement transition matricesB1=B1η and B2=B2η are introduced:(8a)B1=b1Nu,(8b)B2=b2Nu.η.The strain fieldε is expressed in the scaled boundary coordinates as(9)ε=B1uξ,ξ+1ξB2uξ.The stress field can be express as(10)σ=DB1uξ,ξ+1ξB2uξ,where D is an elasticity matrix and has the form(11)D=E1−v21v0v10001−v2,and for plane stress cases, E is Young’s modulus, and v is Poisson’s ratio.According to the virtual work principle, the radial displacements functionuξ is the solution of the SBFEM equation in displacement [36]:(12)E0ξ2uξ,ξξ+E0+E1T−E1ξuξ,ξ−E2uξ+Fξ=0,where uξ,ξξ is the second partial derivative with respect to the variable ξ, uξ,ξ is the first partial derivative with respect to variable ξ, Ei, i=0,1,2, are the coefficient matrices, and Fξ is a load vector.In the derivation, the governing equations of linear elasticity are weakened in the circumferential direction, while the strong form remains in the radial direction. The coefficient matricesEi, i=0,1,2, depends only on the geometry and material properties. These coefficient matrices can be given as follows:(13a)E0=∫−1+1B1TDB1Jbdη,(13b)E1=∫−1+1B2TDB1Jbdη,(13c)E2=∫−1+1B2TDB2Jbdη.WhenFξ=0, the solution of equation (12) is second-order ordinary differential equations that can be obtained by introducing the variable:(14)Xξ=uξqξ,where qξ is the internal force vector. Equation (12) can be transformed into a first-order ordinary differential equation with twice the number of an unknown as(15)ξXξ,ξ−ZpXξ=0,where the coefficient matrix Zp is a Hamiltonian matrix. The solution for a bounded domain is obtained using the positive eigenvalues of Zp. Hence, Zp can be expressed as(16)Zp=−E0−1E1TE0−1E2−E1E0−1E1TE1E0−1.The solution of the equation (15) can be obtained by the mathematical theorem and eigenvalue decomposition technique. The eigenvalue decomposition of Zp is expressed as(17)ZpΦunΦupΦqnΦqp=ΦunΦupΦqnΦqpλn00λp.The real parts of eigenvaluesλn are negative, and that of λp are positive. Φun and Φup are the transformation matrices corresponding to the modal displacements and forces, respectively. The general solution of equation (15) can be written as(18)Xξ=ΦunΦupΦqnΦqpξ−λn00ξ−λpcncp,where cn and cp are the integration constants. For a super element with 0≤ξ≤1, the solution of equation (15) is(19)uξ=Φunξ−λncn,(20)qξ=Φqnξ−λncn,where the integration constants cn can be extracted from the nodal displacements on the boundary ub=uξ=1 as(21)cn=Φun−1ub.Eliminating the integration constantscn at ξ=1 is given(22)F=qξ=1=ΦqnΦun−1ub.Hence, the stiffness matrix of the S-element can be written as(23)K=ΦqnΦun−1.The displacement fielduξ,η inside subdomains by a line element on the S-element can be given (24)uξ,η=NuηΦunξ−λncn.


## 2.2. The Element of PSBFEM

The PSBFEM is inherently appropriate for modeling polygons and has other promising capabilities. Because the PSBFEM is discretized only in the boundary, and an element of PSBFEM can assume more complex shapes than a finite element method, the PSBFEM element is much more flexible than the FEM element. With such an important advantage, the PSBFEM is more applicable to complex geometries than other methods. Figure3 shows the supported element type in the PSBFEM-UEL. In the PSBFEM-UEL, we can use the Abaqus element type, such as CPS3, CPE3, CPS4, and CPE4. Besides, the PSBFEM-UEL also provides the complex element type: polygonal elements and complex quadrilateral element (quadtree discretization). Hence, it is an effective tool to solve complex elements in numerical analysis.Figure 3
The supported element type in the PSBFEM-UEL; (a) a polygon discretization illustrating the mesh generation using PSBFEM-UEL element; (b) a quadtree discretization illustrating the mesh generation using PSBFEM-UEL element; (c) the PSBFEM-UEL element library.
(a)(b)(c)


## 2.3. Automatic Meshes Generation

This section uses a Python script to automatically generate PSBFEM elements by the Delaunay triangulation [37]. This algorithm is robust and efficient, and it has been integrated into Abaqus CAE. The algorithm can be used to triangulate any set of points on a two-dimensional plane. The triangulated mesh is then used to generate the polygon elements. The detail of the polygon generation algorithm mainly contains two steps. Firstly, we automatically generate triangular mesh by the Abaqus CAE. Secondly, considering each triangle interior node as the center of the polygon element, a polygonal element can be generated by connecting the centroids of all the triangular elements circumventing. More detail of this algorithm is presented in Figure 4. Finally, the PSBFEM-UEL can be directly used with a polygon mesh generator to analyze complex geometry problems.Figure 4
Generation of a polygon element from a triangular element: (a) arbitrary geometric model; (b) generating triangular meshes by Abaqus CAE; (c) construction of polygon meshes; and (d) final SBFEM polygon meshes.The quadtree decomposition is a tree data structure in which each parent has precisely four children [38]. The quadtree meshes are fast, efficient, and capable of achieving rapid and smooth transitions of element sizes between mesh refinement regions. The quadtree mesh can provide vital support in the preprocessing for the adaptive analysis of the SBFEM. We developed a quadtree mesh automatic generation code by a Python script. A simple example of quadtree discretization with three levels is illustrated in Figure 5. The quadtree algorithm is described in more detail in [28, 38].Figure 5
A quadtree discretization with three levels.


## 3. UEL Implementation of SBFEM in Abaqus

The Abaqus/Standard analysis provides a programming interface UEL to define the customized elements. In this section, we present the major implementation details of UEL for the PSBFEM in Abaqus. The Abaqus solver can be carried out with the UEL subroutine by the command:(25)abaqusjob=inputfilenameuser=UELsubroutinefile.Figure6 shows the framework sketch of implementing the PSBFEM within Abaqus. This system contains three parts: preprocessing module, PSBFEM-UEL module, and postprocessing module. A Python script is used in preprocessing to generate polygon SBFEM mesh and define the loading and boundary conditions. Finally, we can obtain an Abaqus input file by preprocessing. The Abaqus CAE does not support the visualization of the UEL element. Thus, a Python script is provided in the postprocessing to extract the VTU format results and visualize them in the software Paraview [39].Figure 6
Framework sketch of implementing the PSBFEM within Abaqus.


### 3.1. The Implementation of PSBFEM-UEL

The most critical work of UEL is to update the contribution of the element to the internal force vector RHS and the stiffness matrix AMATRX according to the information from ABAQUS/Standard analysis. In this paper, the UEL code for implementing the PSBFEM is written in FORTRAN77. Figure7 shows an overall implementation of the UEL subroutine. Based on the input file’s connectivity information, the UEL computes the scaling centers and transforms the global coordinate into the local coordinate. Equation (13) computes the coefficient matrices E0, E1, and E2, which are used to construct the Hamilton matrix Zp by equation (16). The two eigenvector matrices Φqn,Φun are constructed by eigenvalue decomposition. Finally, we can obtain the stiffness matrix K of the PSBFEM element. Moreover, to avoid calculating the stiffness matrix K at each incremental step and decrease the calculation cost, we only calculate the stiffness matrix K at the first incremental step and store it in the state variable. The stiffness matrix K is read directly at the next incremental step.Figure 7
Framework sketch of implementing the UEL.To solve the stiffness matrix, we need to employ the eigenvalue decomposition (see equation (17)). At present, many mathematical libraries to perform the eigenvalue decomposition exist. In this work, we use the Intel Math Kernel Library (MKL) [35] to decompose the eigenvalue. In the Abaqus/Standard analysis, we can directly use MKL by modifying the Abaqus environment file.


### 3.2. Defining the UEL Elements

Abaqus’s input file usually contains a model information section (such as defining nodes, elements, the active degrees of freedom, and materials). At present, this information cannot be set in Abaqus CAE and must be defined through an input file. The Abaqus provides the keyword∗USER ELEMENT to defined a new user element. The main contents of defining the user element are as follows:(1)
Assigning an element type key to a user-defined element and the number of nodes. The element type key must be of the form “Un” in Abaqus/Standard analysis, where′n′ is a positive integer that identifies the element type uniquely. In this implementation of PSBFEM, the integer is equal to the number of nodes of the element.(2)
Defining the element properties and assigning an Abaqus material to the user element.(3)
Defining the number of degrees of freedom per node and the active degrees of freedom at the nodes.Listing 1: user element definition for the arbitrary polygons c.f. Figure8.(1)
∗USER ELEMENT, NODES = 3, TYPE = U3, PROPERTIES = 2, COORDINATES = 2(2)
1, 2(3)
∗ELEMENT, TYPE = U3, ELSET = E3(4)
2,1,2,7(5)
4,5,6,4(6)
∗UEL PROPERTY, ELEST = E3(7)
1000, 0.2(8)
∗USER ELEMENT, NODES = 4, TYPE = U4, PROPERTIES = 2, COORDINATES = 2(9)
1,2(10)
∗ELEMENT, TYPE = U4, ELSET = E4(11)
1,1,7,6,5(12)
∗UEL PROPERTY, ELEST = E1(13)
1000,0.2(14)
∗USER ELEMENT, NODES = 5, TYPE = U5, PROPERTIES = 2, COORDINATES = 2(15)
1,2(16)
∗ELEMENT, TYPE = U5, ELSET = E5(17)
3,2,3,4,6,7(18)
∗UEL PROPERTY, ELEST = E1(19)
1000,0.2Figure 8
A simple example of the SBFEM polygonal mesh.In the case of PSBFEM, we present a simple polygonal mesh of PSBFEM (see Figure8). This mesh consists of three element types: triangular element (U3), quadrilateral element (U4), and Pentagon element (U5). In the input file (see listing 1), 1∼19 is the line number; the actual input file does not contain the line number. Lines 1∼7 are used to define two triangular elements (U3). Line 1 assigns the element type, the number of nodes, the number of element properties, and the number of freedom degrees per node. Line 2 sets the active degrees of freedom. Lines 3∼5 define the element sets’ E3’. Lines 6∼7 set the element properties (Young’s modulus and Poisson’s ratio) of “E3”. Similarly, other types of elements are defined using the same way.


## 3.1. The Implementation of PSBFEM-UEL

The most critical work of UEL is to update the contribution of the element to the internal force vector RHS and the stiffness matrix AMATRX according to the information from ABAQUS/Standard analysis. In this paper, the UEL code for implementing the PSBFEM is written in FORTRAN77. Figure7 shows an overall implementation of the UEL subroutine. Based on the input file’s connectivity information, the UEL computes the scaling centers and transforms the global coordinate into the local coordinate. Equation (13) computes the coefficient matrices E0, E1, and E2, which are used to construct the Hamilton matrix Zp by equation (16). The two eigenvector matrices Φqn,Φun are constructed by eigenvalue decomposition. Finally, we can obtain the stiffness matrix K of the PSBFEM element. Moreover, to avoid calculating the stiffness matrix K at each incremental step and decrease the calculation cost, we only calculate the stiffness matrix K at the first incremental step and store it in the state variable. The stiffness matrix K is read directly at the next incremental step.Figure 7
Framework sketch of implementing the UEL.To solve the stiffness matrix, we need to employ the eigenvalue decomposition (see equation (17)). At present, many mathematical libraries to perform the eigenvalue decomposition exist. In this work, we use the Intel Math Kernel Library (MKL) [35] to decompose the eigenvalue. In the Abaqus/Standard analysis, we can directly use MKL by modifying the Abaqus environment file.


## 3.2. Defining the UEL Elements

Abaqus’s input file usually contains a model information section (such as defining nodes, elements, the active degrees of freedom, and materials). At present, this information cannot be set in Abaqus CAE and must be defined through an input file. The Abaqus provides the keyword∗USER ELEMENT to defined a new user element. The main contents of defining the user element are as follows:(1)
Assigning an element type key to a user-defined element and the number of nodes. The element type key must be of the form “Un” in Abaqus/Standard analysis, where′n′ is a positive integer that identifies the element type uniquely. In this implementation of PSBFEM, the integer is equal to the number of nodes of the element.(2)
Defining the element properties and assigning an Abaqus material to the user element.(3)
Defining the number of degrees of freedom per node and the active degrees of freedom at the nodes.Listing 1: user element definition for the arbitrary polygons c.f. Figure8.(1)
∗USER ELEMENT, NODES = 3, TYPE = U3, PROPERTIES = 2, COORDINATES = 2(2)
1, 2(3)
∗ELEMENT, TYPE = U3, ELSET = E3(4)
2,1,2,7(5)
4,5,6,4(6)
∗UEL PROPERTY, ELEST = E3(7)
1000, 0.2(8)
∗USER ELEMENT, NODES = 4, TYPE = U4, PROPERTIES = 2, COORDINATES = 2(9)
1,2(10)
∗ELEMENT, TYPE = U4, ELSET = E4(11)
1,1,7,6,5(12)
∗UEL PROPERTY, ELEST = E1(13)
1000,0.2(14)
∗USER ELEMENT, NODES = 5, TYPE = U5, PROPERTIES = 2, COORDINATES = 2(15)
1,2(16)
∗ELEMENT, TYPE = U5, ELSET = E5(17)
3,2,3,4,6,7(18)
∗UEL PROPERTY, ELEST = E1(19)
1000,0.2Figure 8
A simple example of the SBFEM polygonal mesh.In the case of PSBFEM, we present a simple polygonal mesh of PSBFEM (see Figure8). This mesh consists of three element types: triangular element (U3), quadrilateral element (U4), and Pentagon element (U5). In the input file (see listing 1), 1∼19 is the line number; the actual input file does not contain the line number. Lines 1∼7 are used to define two triangular elements (U3). Line 1 assigns the element type, the number of nodes, the number of element properties, and the number of freedom degrees per node. Line 2 sets the active degrees of freedom. Lines 3∼5 define the element sets’ E3’. Lines 6∼7 set the element properties (Young’s modulus and Poisson’s ratio) of “E3”. Similarly, other types of elements are defined using the same way.


## 4. Numerical Examples

In this section, we validate the convergence and accuracy of the implementation by solving a few benchmark problems. Moreover, the results of the PSBFEM are compared with the FEM. The FEM analysis uses the commercial finite element software Abaqus. For validation, the relative errorL2 norms in the displacement are computed as follows:(26)eu=u−uL2Ωh=∫Ωu−uhTu−uhdΩ∫ΩuTudΩ,where uh is the numerical solution and u is the analytical or reference solution.The relative error in the energy as given by(27)eenery=ε−εL2Ωh=∫Ωε−εhTDε−εhdΩ∫ΩεTDεdΩ,where εh is the numerical solution and ε is the analytical or reference solution.


### 4.1. A Two-Dimensional Cantilever Beam

A two-dimensional cantilever beam of heightH=1.0 m and length L=5.0 m subjected to a uniform loading P=10kPa is considered, as shown in Figure 9(a). The material properties are given by Young’s modulus E=2GPa and Poisson’s ratio v=0. The right boundaries are constrained without displacement ΔX=0,ΔY=0. The domain is discretized with the quadrangle and arbitrary polygonal elements. A convergence study is performed by mesh refinement. The meshes are refined successively following the sequence n=2,4,8,16,32. The element size is chosen as h=H/n. A representative mesh is presented in Figures 9(b) and 9(c). In this work, the Abaqus uses the CPS4 element. The analytical solution of the stress fields can be expressed as [40–42](28a)σx=qyh4y2h2−35−6x2h2,(28b)σy=−q21−3yh+4y3h3,(28c)τxy=−3q2xh1−4y2h2.Figure 9
A cantilever beam subjected to a uniform loading: (a) geometry and boundary conditions; (b) the quadrilateral mesh; (c) the polygonal mesh; the discretization of the domain by a Python script, and visualization mesh by the software Paraview.
(a)(b)(c)The results of the relative error in the vertical displacementuy for the point O are given in Table 1. It can be observed that the errors decrease as the mesh refinement. The errors of PSBFEM are less than the FEM at the same element size. Moreover, the polygonal element shows higher accuracy than the quadrilateral element in the PSBFEM because the polygonal element has more nodes for the same element size. Figure 10 shows that the results of mesh size sensitivity for three element types. It is clear that these elements would obtain more accurate results when mesh sizeis less than 1/8 m. The convergence of the relative error in the displacement and the energy norm with mesh refinement are presented in Figures 11 and 12. It is observed that the PSBFEM converges to an exact solution with an optimal convergence rate.Table 1
The relative error in the vertical displacementuy of a cantilever beam subjected to a uniform loading with mesh refinement.
Element size (m)FEMPSBFEM-QuadPSBFEM-PolygonDOFseu%DOFseu%DOFseu%0.5661.13 × 10−1668.75 × 10−2828.20 × 10−20.252103.11 × 10−22102.37 × 10−23221.86 × 10−20.1257388.01 × 10−37386.08 × 10−312823.30 × 10−30.062527542.02 × 10−327541.53 × 10−351247.73 × 10−40.03125106265.07 × 10−4106263.83 × 10−4204801.96 × 10−4Figure 10
Sensitivity analysis of mesh sizes (Uanaly denotes the analytical solution of displacement; Unum denotes the numerical solution of displacement).Figure 11
Convergence of the relative error in the displacements of a cantilever beam is subjected to uniform loading.Figure 12
Convergence of relative error in the energy for a cantilever beam is subjected to uniform loading.Moreover, Figure13 shows the contours of the vertical displacement uy of Abaqus CPS4, PSBFEM-Quad, and PSBFEM polygon. It is clearly shown that the results are virtually the same for the FEM and PSBFEM. In addition, Figure 14 presents the computational cost comparison of the developed UEL and Abaqus standard elements. Due to the increment size setting affecting the solving time, we use the automatic incrementation type. The comparison is evaluated with an Intel Core i7-4710MQ CPU running at 2.50 GHz and 4.0 GB of RAM. The total CPU time is normalized. It is noted that computation costs for the PSBFEM-UEL and Abaqus standard elements are comparable and of the same trend. Moreover, the computational cost of PSBFEM-UEL is slightly higher than the Abaqus standard element. It is due to the semianalytical method of SBFEM being bound to generate relatively more computations in its concept than FEM [43].Figure 13
Contour plots of displacementsuy in a cantilever beam: (a) Abaqus CPS4; (b) PSBFEM Quad; (c) PSBFEM Polygon.
(a)(b)(c)Figure 14
The time consumption comparison of PSBFEM-UEL between Abaqus standard elements.


### 4.2. Infinite Plate with a Circular Hole

In this example, an infinite plate with a circular hole of radiusa under remote uniaxial tension σx is considered as shown in Figure 15(a). Considering the symmetry of geometry, we only analyze a quarter of the infinite plate, as shown in Figure 15(b). The analytical solution of stress fields in the polar coordinates r,θ is expressed as [36](29a)σx=p22−a2r23cos2θ+2−3a2r2cos4θ,(29b)σy=−pa22r2cos2θ−2−3a2r2cos4θ,(29c)τxy=−pa22r2sin2θ+2−3a2r2sin4θ,and the displacement fields can be expressed as(30a)ux=Pa8Gra1+κcosθ+2ar1+κcosθ+cos3θ−2a3r3cos3θ,(30b)uy=Pa8Graκ−3sinθ+2ar1−κsinθ+sin3θ−2a3r3sin3θ.Figure 15
Infinite plate with a circular hole: geometry and boundary conditions: (a) a circular hole; (b) a quarter of the circular hole; (c) the quadrilateral mesh; and (d) the polygonal mesh.
(a)(b)(c)(d)The shear modulusG and Kolosov’s constant κ can be expressed as(31)G=2E21+v,k=3−v1+v,where a is the radius of the hole. The circular hole radius is 0.4 m, the length of the quarter of the circle hole is 1.0 m, and remote tension σx=1 kPa. The domain is discretized using quadrilateral and polygonal meshes, respectively. The problem is also modeled by the plane stress, and the material properties are given by Young’s modulus E=1×105 Pa and Poisson’s ratio v=0.25. The vertical displacement of bottom boundaries ΔY=0 and the horizontal displacement of left boundaries ΔX=0 are constrained without displacement.The results of PSBFEM and FEM in terms of convergence behavior and accuracy are compared and demonstrated. Figures16 and 17 show the convergence of the relative error norms decrease in displacement and the energy norm with increasing elements. It is observed that the PSBFEM element is significantly more accurate than the FEM element. Besides, PSBFEM Polygon has a more slightly fast convergence rate.Figure 16
Convergence of the relative error in the displacement for an infinite plate with a circular hole.Figure 17
Convergence of the relative error in the energy for an infinite plate with a circular hole.From the results, it is clear that all the methods asymptotically converge to the analytical solution with mesh refinement. It is noted that the PSBFEM requires fewer DOFs when compared to conventional FEM. Besides, Figure18(a) provides that the contours of the vertical displacement uy of Abaqus CPS4, PSBFEM-Quad, and PSBFEM polygon in an infinite plate with a circular hole. The results show a good agreement for the FEM and PSBFEM. Moreover, Figures 18(b) and 18(c) also show that the distribution of stress component σy and strain component εyobtained from the PSBFEM and the FEM by Abaqus. Similarly, good agreement between the three sets of results is observed.Figure 18
Contour plots of results in an infinite plate with a circular hole: (a) vertical displacementuy; (b) stress component σy; and (c) strain component εy.


### 4.3. Pressurized Thick Cylinder Modeled by a Quarter-Annulus Model

In this example, we consider a benchmark problem of a quarter-thick cylinder subjected to internal pressure at the inner circular edge. The analytical solution of stress fields in the polar coordinates is [44](32a)σrr=Ra2PRa2−Rb21−Ra2r2,(32b)σθθ=Ra2PRa2−Rb21+Ra2r2,(32c)σrθ=0,and the displacement fields can be expressed as(33a)urad=Ra2PrERb2−Ra21−ν+Rbr21+ν,(33b)ux=uradcosθ,(33c)uy=uradsinθ.where Ra and Rb are the inner and the outer radius of the cylinder, respectively, with Ra=1 m and Rb=2 m. P is the pressure exerted along the inner circular edge, withP=1000Pa. The material properties are given by Young’s modulus E=1×105 Pa and Poisson’s ratio v=0.30. The geometry and boundary conditions are presented in Figure 19(a). The domain is discretized with the quadrilateral and polygonal elements, as shown in Figures 19(b) and 19(c).Figure 19
A quarter-thick cylinder subjected to internal pressure at the inner circular edge: (a) the geometry and boundary conditions for a quarter-annulus model; (b) the quadrilateral mesh; and (c) the polygonal mesh.
(a)(b)(c)Figures20 and 21 show the convergence of the relative error in both the displacement and the energy norm with mesh refinement. Results show that all the techniques asymptotically converge to analytical solutions with reducing element size, and the PSBFEM Polygon shows slightly accurate results. Moreover, a good agreement between PSFEM and FEM is observed in Figure 22.Figure 20
Convergence of the relative error in the displacement for a quarter-thick cylinder.Figure 21
Convergence of the relative error in the energy for a quarter-thick cylinder.Figure 22
Contour plots of resultant displacements in a quarter-thick cylinder: (a) Abaqus CPS4; (b) PSBFEM Quad; and (c) PSBFEM Polygon.
(a)(b)(c)


### 4.4. Square Body with Multiple Holes

To highlight the flexibility of the PSBFEM element and quadtree algorithm in handling mesh, a unit square bodyL=1.0m with four circular holes is considered, as shown in Figure 23(a). The material properties are Young’s modulus E=10 GPa and Poisson’s ratio v=0.25. The bottom edge is constrained without displacement ΔX=0,ΔY=0, and a uniform tension P=1 MPa is applied on the top edge of the square. The four circular holes locate X=0.5m,Y=0.75m,r=0.05m, X=0.25m,Y=0.5m,r=0.15m, X=0.5m,Y=0.25m,r=0.05m, and X=0.75m,Y=0.5m,r=0.15m, respectively, to form a hole cluster. For the error estimation, the relative error L2 norms in the displacement are computed as equation (26). Moreover, a simple regression is investigated to compare the results of FEM and PSBFEM.Figure 23
Schematic diagram of the square body with multiple holes under uniaxial tension: (a) the geometry and boundary conditions; (b) Abaqus mesh; and (c) Quadtree mesh.
(a)(b)(c)The quadtree mesh is generated by setting the same number of mesh seeds on every hole, as shown in Figure23(c). It is clearly present that the mesh transition between the holes of different sizes is effectively handled. The square body with multiple holes is modeled with 1632 quadtree elements, and the total nodes are 2478. Moreover, this problem is also analyzed with a similar number of nodes (2532 nodes) using the Abaqus CPS4 element. Due to the symmetry of geometry, only the results on the right side of the top edge are provided for comparison. Results in Figure 24 show the comparison between quadtree PSBFEM and Abaqus CPS4 in the vertical displacement uy. A simple regression on the Abaqus CPS4 and quadtree PSBFEM had presented R2  = 0.99, and the error in the L2 norm is 0.0097. The strain energy of FEM is 85.934 kJ, and the strain energy of PSBFEM is 86.5917 kJ. The relative error of strain energy is 0.76%. Hence, these results indicate PSBFEM accuracy and reliability for the quadtree mesh. Also, the contour plots of the vertical displacement uy obtained from the PSBFEM analysis and the FEM analysis by Abaqus are shown in Figure 25. It is noted that the contour plots present a good agreement.Figure 24
Comparison between quadtree PSBFEM and Abaqus CPS4 in the vertical displacementuy on the top edge.Figure 25
Contour plots of displacementsuy in a square body with multiple holes: (a) Abaqus CPS4 and (b) Quadtree PSBFEM.
(a)(b)


### 4.5. A Square Plate for the Complex Geometry

A square body with a rabbit shape cavity, as shown in Figure26. The size of the square plate is L=4m. The material properties are Young’s modulus E=1 MPa and Poisson’s ratio v=0.25. The left edge is constrained without displacement ΔX=0,ΔY=0, and the right edge is applied a horizontal displacement U=0.1m, as shown in Figure 26. We chose four nodes, A, B, C, D, to compare results, as shown in Figure 26. The square body with a rabbit shape cavity is modeled with 3147 quadtree elements. Moreover, this problem is also analyzed with a similar number of elements using the Abaqus CPS4 element, and the total elements are 3358, as shown in Figure 27.Figure 26
The geometry and boundary conditions of the square body with a rabbit shape cavity.Figure 27
The meshes of a square plate with a rabbit shape cavity: (a) Abaqus mesh and (b) Quadtree mesh.
(a)(b)Table2 shows the relative error in the horizontal displacement for the PSBFEM and FEM. The relative errors are less than 0.5%. Moreover, the strain energy of FEM is 4.304 kJ, and the strain energy of PSBFEM is 4.299 kJ. The relative error of strain energy is 0.12%. Figure 28 presents that the contour plots present a good agreement for FEM and PSBFEM. Therefore, the quadtree mesh of PSBFEM has sufficient accuracy and reliability.Table 2
The relative error for the PSBFEM and FEM in the horizontal displacement.
PointHorizontal displacementuxFEM (m)PSBFEM (m)eu%A0.050420.050320.20B0.012700.012640.47C0.049590.049560.06D0.086000.086160.19Figure 28
Contour plots of horizontal displacements in a square plate with a rabbit shape cavity: (a) Abaqus CPS4 and (b) Quadtree PSBFEM.
(a)(b)


## 4.1. A Two-Dimensional Cantilever Beam

A two-dimensional cantilever beam of heightH=1.0 m and length L=5.0 m subjected to a uniform loading P=10kPa is considered, as shown in Figure 9(a). The material properties are given by Young’s modulus E=2GPa and Poisson’s ratio v=0. The right boundaries are constrained without displacement ΔX=0,ΔY=0. The domain is discretized with the quadrangle and arbitrary polygonal elements. A convergence study is performed by mesh refinement. The meshes are refined successively following the sequence n=2,4,8,16,32. The element size is chosen as h=H/n. A representative mesh is presented in Figures 9(b) and 9(c). In this work, the Abaqus uses the CPS4 element. The analytical solution of the stress fields can be expressed as [40–42](28a)σx=qyh4y2h2−35−6x2h2,(28b)σy=−q21−3yh+4y3h3,(28c)τxy=−3q2xh1−4y2h2.Figure 9
A cantilever beam subjected to a uniform loading: (a) geometry and boundary conditions; (b) the quadrilateral mesh; (c) the polygonal mesh; the discretization of the domain by a Python script, and visualization mesh by the software Paraview.
(a)(b)(c)The results of the relative error in the vertical displacementuy for the point O are given in Table 1. It can be observed that the errors decrease as the mesh refinement. The errors of PSBFEM are less than the FEM at the same element size. Moreover, the polygonal element shows higher accuracy than the quadrilateral element in the PSBFEM because the polygonal element has more nodes for the same element size. Figure 10 shows that the results of mesh size sensitivity for three element types. It is clear that these elements would obtain more accurate results when mesh sizeis less than 1/8 m. The convergence of the relative error in the displacement and the energy norm with mesh refinement are presented in Figures 11 and 12. It is observed that the PSBFEM converges to an exact solution with an optimal convergence rate.Table 1
The relative error in the vertical displacementuy of a cantilever beam subjected to a uniform loading with mesh refinement.
Element size (m)FEMPSBFEM-QuadPSBFEM-PolygonDOFseu%DOFseu%DOFseu%0.5661.13 × 10−1668.75 × 10−2828.20 × 10−20.252103.11 × 10−22102.37 × 10−23221.86 × 10−20.1257388.01 × 10−37386.08 × 10−312823.30 × 10−30.062527542.02 × 10−327541.53 × 10−351247.73 × 10−40.03125106265.07 × 10−4106263.83 × 10−4204801.96 × 10−4Figure 10
Sensitivity analysis of mesh sizes (Uanaly denotes the analytical solution of displacement; Unum denotes the numerical solution of displacement).Figure 11
Convergence of the relative error in the displacements of a cantilever beam is subjected to uniform loading.Figure 12
Convergence of relative error in the energy for a cantilever beam is subjected to uniform loading.Moreover, Figure13 shows the contours of the vertical displacement uy of Abaqus CPS4, PSBFEM-Quad, and PSBFEM polygon. It is clearly shown that the results are virtually the same for the FEM and PSBFEM. In addition, Figure 14 presents the computational cost comparison of the developed UEL and Abaqus standard elements. Due to the increment size setting affecting the solving time, we use the automatic incrementation type. The comparison is evaluated with an Intel Core i7-4710MQ CPU running at 2.50 GHz and 4.0 GB of RAM. The total CPU time is normalized. It is noted that computation costs for the PSBFEM-UEL and Abaqus standard elements are comparable and of the same trend. Moreover, the computational cost of PSBFEM-UEL is slightly higher than the Abaqus standard element. It is due to the semianalytical method of SBFEM being bound to generate relatively more computations in its concept than FEM [43].Figure 13
Contour plots of displacementsuy in a cantilever beam: (a) Abaqus CPS4; (b) PSBFEM Quad; (c) PSBFEM Polygon.
(a)(b)(c)Figure 14
The time consumption comparison of PSBFEM-UEL between Abaqus standard elements.


## 4.2. Infinite Plate with a Circular Hole

In this example, an infinite plate with a circular hole of radiusa under remote uniaxial tension σx is considered as shown in Figure 15(a). Considering the symmetry of geometry, we only analyze a quarter of the infinite plate, as shown in Figure 15(b). The analytical solution of stress fields in the polar coordinates r,θ is expressed as [36](29a)σx=p22−a2r23cos2θ+2−3a2r2cos4θ,(29b)σy=−pa22r2cos2θ−2−3a2r2cos4θ,(29c)τxy=−pa22r2sin2θ+2−3a2r2sin4θ,and the displacement fields can be expressed as(30a)ux=Pa8Gra1+κcosθ+2ar1+κcosθ+cos3θ−2a3r3cos3θ,(30b)uy=Pa8Graκ−3sinθ+2ar1−κsinθ+sin3θ−2a3r3sin3θ.Figure 15
Infinite plate with a circular hole: geometry and boundary conditions: (a) a circular hole; (b) a quarter of the circular hole; (c) the quadrilateral mesh; and (d) the polygonal mesh.
(a)(b)(c)(d)The shear modulusG and Kolosov’s constant κ can be expressed as(31)G=2E21+v,k=3−v1+v,where a is the radius of the hole. The circular hole radius is 0.4 m, the length of the quarter of the circle hole is 1.0 m, and remote tension σx=1 kPa. The domain is discretized using quadrilateral and polygonal meshes, respectively. The problem is also modeled by the plane stress, and the material properties are given by Young’s modulus E=1×105 Pa and Poisson’s ratio v=0.25. The vertical displacement of bottom boundaries ΔY=0 and the horizontal displacement of left boundaries ΔX=0 are constrained without displacement.The results of PSBFEM and FEM in terms of convergence behavior and accuracy are compared and demonstrated. Figures16 and 17 show the convergence of the relative error norms decrease in displacement and the energy norm with increasing elements. It is observed that the PSBFEM element is significantly more accurate than the FEM element. Besides, PSBFEM Polygon has a more slightly fast convergence rate.Figure 16
Convergence of the relative error in the displacement for an infinite plate with a circular hole.Figure 17
Convergence of the relative error in the energy for an infinite plate with a circular hole.From the results, it is clear that all the methods asymptotically converge to the analytical solution with mesh refinement. It is noted that the PSBFEM requires fewer DOFs when compared to conventional FEM. Besides, Figure18(a) provides that the contours of the vertical displacement uy of Abaqus CPS4, PSBFEM-Quad, and PSBFEM polygon in an infinite plate with a circular hole. The results show a good agreement for the FEM and PSBFEM. Moreover, Figures 18(b) and 18(c) also show that the distribution of stress component σy and strain component εyobtained from the PSBFEM and the FEM by Abaqus. Similarly, good agreement between the three sets of results is observed.Figure 18
Contour plots of results in an infinite plate with a circular hole: (a) vertical displacementuy; (b) stress component σy; and (c) strain component εy.


## 4.3. Pressurized Thick Cylinder Modeled by a Quarter-Annulus Model

In this example, we consider a benchmark problem of a quarter-thick cylinder subjected to internal pressure at the inner circular edge. The analytical solution of stress fields in the polar coordinates is [44](32a)σrr=Ra2PRa2−Rb21−Ra2r2,(32b)σθθ=Ra2PRa2−Rb21+Ra2r2,(32c)σrθ=0,and the displacement fields can be expressed as(33a)urad=Ra2PrERb2−Ra21−ν+Rbr21+ν,(33b)ux=uradcosθ,(33c)uy=uradsinθ.where Ra and Rb are the inner and the outer radius of the cylinder, respectively, with Ra=1 m and Rb=2 m. P is the pressure exerted along the inner circular edge, withP=1000Pa. The material properties are given by Young’s modulus E=1×105 Pa and Poisson’s ratio v=0.30. The geometry and boundary conditions are presented in Figure 19(a). The domain is discretized with the quadrilateral and polygonal elements, as shown in Figures 19(b) and 19(c).Figure 19
A quarter-thick cylinder subjected to internal pressure at the inner circular edge: (a) the geometry and boundary conditions for a quarter-annulus model; (b) the quadrilateral mesh; and (c) the polygonal mesh.
(a)(b)(c)Figures20 and 21 show the convergence of the relative error in both the displacement and the energy norm with mesh refinement. Results show that all the techniques asymptotically converge to analytical solutions with reducing element size, and the PSBFEM Polygon shows slightly accurate results. Moreover, a good agreement between PSFEM and FEM is observed in Figure 22.Figure 20
Convergence of the relative error in the displacement for a quarter-thick cylinder.Figure 21
Convergence of the relative error in the energy for a quarter-thick cylinder.Figure 22
Contour plots of resultant displacements in a quarter-thick cylinder: (a) Abaqus CPS4; (b) PSBFEM Quad; and (c) PSBFEM Polygon.
(a)(b)(c)


## 4.4. Square Body with Multiple Holes

To highlight the flexibility of the PSBFEM element and quadtree algorithm in handling mesh, a unit square bodyL=1.0m with four circular holes is considered, as shown in Figure 23(a). The material properties are Young’s modulus E=10 GPa and Poisson’s ratio v=0.25. The bottom edge is constrained without displacement ΔX=0,ΔY=0, and a uniform tension P=1 MPa is applied on the top edge of the square. The four circular holes locate X=0.5m,Y=0.75m,r=0.05m, X=0.25m,Y=0.5m,r=0.15m, X=0.5m,Y=0.25m,r=0.05m, and X=0.75m,Y=0.5m,r=0.15m, respectively, to form a hole cluster. For the error estimation, the relative error L2 norms in the displacement are computed as equation (26). Moreover, a simple regression is investigated to compare the results of FEM and PSBFEM.Figure 23
Schematic diagram of the square body with multiple holes under uniaxial tension: (a) the geometry and boundary conditions; (b) Abaqus mesh; and (c) Quadtree mesh.
(a)(b)(c)The quadtree mesh is generated by setting the same number of mesh seeds on every hole, as shown in Figure23(c). It is clearly present that the mesh transition between the holes of different sizes is effectively handled. The square body with multiple holes is modeled with 1632 quadtree elements, and the total nodes are 2478. Moreover, this problem is also analyzed with a similar number of nodes (2532 nodes) using the Abaqus CPS4 element. Due to the symmetry of geometry, only the results on the right side of the top edge are provided for comparison. Results in Figure 24 show the comparison between quadtree PSBFEM and Abaqus CPS4 in the vertical displacement uy. A simple regression on the Abaqus CPS4 and quadtree PSBFEM had presented R2  = 0.99, and the error in the L2 norm is 0.0097. The strain energy of FEM is 85.934 kJ, and the strain energy of PSBFEM is 86.5917 kJ. The relative error of strain energy is 0.76%. Hence, these results indicate PSBFEM accuracy and reliability for the quadtree mesh. Also, the contour plots of the vertical displacement uy obtained from the PSBFEM analysis and the FEM analysis by Abaqus are shown in Figure 25. It is noted that the contour plots present a good agreement.Figure 24
Comparison between quadtree PSBFEM and Abaqus CPS4 in the vertical displacementuy on the top edge.Figure 25
Contour plots of displacementsuy in a square body with multiple holes: (a) Abaqus CPS4 and (b) Quadtree PSBFEM.
(a)(b)


## 4.5. A Square Plate for the Complex Geometry

A square body with a rabbit shape cavity, as shown in Figure26. The size of the square plate is L=4m. The material properties are Young’s modulus E=1 MPa and Poisson’s ratio v=0.25. The left edge is constrained without displacement ΔX=0,ΔY=0, and the right edge is applied a horizontal displacement U=0.1m, as shown in Figure 26. We chose four nodes, A, B, C, D, to compare results, as shown in Figure 26. The square body with a rabbit shape cavity is modeled with 3147 quadtree elements. Moreover, this problem is also analyzed with a similar number of elements using the Abaqus CPS4 element, and the total elements are 3358, as shown in Figure 27.Figure 26
The geometry and boundary conditions of the square body with a rabbit shape cavity.Figure 27
The meshes of a square plate with a rabbit shape cavity: (a) Abaqus mesh and (b) Quadtree mesh.
(a)(b)Table2 shows the relative error in the horizontal displacement for the PSBFEM and FEM. The relative errors are less than 0.5%. Moreover, the strain energy of FEM is 4.304 kJ, and the strain energy of PSBFEM is 4.299 kJ. The relative error of strain energy is 0.12%. Figure 28 presents that the contour plots present a good agreement for FEM and PSBFEM. Therefore, the quadtree mesh of PSBFEM has sufficient accuracy and reliability.Table 2
The relative error for the PSBFEM and FEM in the horizontal displacement.
PointHorizontal displacementuxFEM (m)PSBFEM (m)eu%A0.050420.050320.20B0.012700.012640.47C0.049590.049560.06D0.086000.086160.19Figure 28
Contour plots of horizontal displacements in a square plate with a rabbit shape cavity: (a) Abaqus CPS4 and (b) Quadtree PSBFEM.
(a)(b)


## 5. Conclusions

This paper implements the PSBFEM of two-dimensional linear elastostatic problems within the Abaqus/Standard analysis by the UEL subroutine. This work mainly focuses on the main procedures to interact with Abaqus, defining the UEL element in the input file, and solving the stiffness matrix by the eigenvalue decomposition in the UEL implementation procedure. Also, we discuss the automatic mesh generation of polygon/quadtree and the visualization of results by the Paraview.The implementation of PSBFEM is validated against the FEM by solving a few benchmark problems. The results demonstrate that PSBFEM-UEL has a significantly better than FEM convergence rate. Moreover, the polygon mesh has a higher accuracy rate than the quadrangle mesh in the PSBFEM-UEL. Notably, the implementation of PSBFEM can conveniently use arbitrary polygon elements by the polygon/quadtree discretizations in the commercial finite element software Abaqus. In the future, the scope of this approach developed here can be extended to higher-order elements. The source code of the implementation can be downloaded fromhttps://github.com/hhupde/PSBFEM-Abaqus with input files of numerical examples presented in this work.



---

*Source: 6628837-2021-09-16.xml*