#### Abstract

The polygonal scaled boundary finite element method (PSBFEM) is a novel approach integrating the standard scaled boundary finite element method and the polygonal mesh technique. In this work, a user-defined element (UEL) for dynamic analysis based on the PSBFEM is developed in the general finite element software ABAQUS. We present the main procedures of interacting with Abaqus, updating AMATRX and RHS, defining the UEL element, and solving the stiffness and mass matrices through eigenvalue decomposition. Several benchmark problems of free and forced vibration are solved to validate the proposed implementation. The results show that the PSBFEM is more accurate than the FEM with mesh refinement. Moreover, the PSBFEM avoids the occurrence of hanging nodes by constructing a polygonal mesh. Thus, the PSBFEM can choose an appropriate mesh resolution for different structures ensuring accuracy and reducing calculation costs.

#### 1. Introduction

Structural dynamic analysis, especially the analysis of the response of a structure to an earthquake, is an important problem in engineering design. In the past, calculations were made for such problems under simplified conditions. The finite element method (FEM) has become the most powerful and routine method of numerically analyzing the time-dependent responses of structures in free and forced vibration problems. Typically, the domain of the conventional two-dimensional FEM is discretized using triangular or quadrilateral elements. The quadrilateral mesh is limited in handling complex geometries [1] whereas triangular elements are more flexible than quadrilateral elements in mesh generation while providing lower accuracy [2]. The generation of a high-quality mesh can thus be a time-consuming task.

The scaled boundary finite element method (SBFEM) was presented in the 1990s by Song and Wolf [3]. The SBFEM is a semianalytical method that attempts to fuse the advantages and characteristics of the FEM and boundary element method (BEM) into one approach. The SBFEM only discretizes at the boundary. Hence, the governing differential equations have a weak form along the circumferential direction whereas they maintain a strong form along the radial direction. The SBFEM has been applied to many physical field problems, such as those of wave propagation [4, 5], heat conduction [6, 7], fracture mechanics [8–11], acoustics [12], seepage [13, 14], elastoplastics [15], and fluids [16, 17].

The SBFEM has also been applied to the analysis of free and forced vibrations. Sepehry et al. [18] adopted the SBFEM to analyze the free and forced vibrations of piezoelectric materials. Liu et al. [19] solved the associated eigenvalue problem to determine the free vibration response of a circular cylindrical piezoelectric panel using the SBFEM. Liu et al. [20] investigated the free vibration and transient dynamic behaviors of functionally graded sandwich plates using the SBFEM. Yang et al. [21] developed a frequency-domain method for modeling general transient linear-elastic dynamic problems using the SBFEM. The SBFEM was shown to be more efficient than the conventional FEM for these problems.

The polygonal scaled boundary finite element method (PSBFEM) is a novel method integrating the standard SBFEM and the polygonal mesh or quadtree mesh technique [22]. Compared with the standard finite element, polygons are not limited in terms of their number of sides [23]. The polygon-based finite element is easily generated through Voronoi [24] and Delaunay tessellations [23]. The ability to assume any irregular shape makes polygon elements more flexible in meshing complex geometries. The discretization thus becomes more flexible for complex geometries [25], such as in modeling polycrystals [26] and biomechanics applications [27]. Furthermore, for topology optimization, Talischi et al. [28] revealed that polygons possess higher degrees of geometric isotropy and allow for greater flexibility in discretizing complex domains without suffering from numerical instabilities. Therefore, these advantages have further promoted the use of polygons as an alternative technique.

The SBFEM has not yet been included in commercial software despite the maturing of its development. It is thus difficult for the engineer to use this method in solving engineering problems. The commercial software Abaqus has powerful linear and nonlinear and static and dynamic analysis capabilities. Abaqus/Standard analysis also provides a user-defined element (UEL) to create elements with an element formulation that is not available in Abaqus/Standard. Molnar et al. [29] implemented an implicit, staggered elastoplastic version of the phase-field approach in Abaqus through the UEL. Kumbhar et al. [30] implemented the element-based smoothed finite element method (CSFEM) using UEL subroutines. Notably, the UEL can be used to extend Abaqus to solve new problems conveniently.

Recent works have focused on the implementation of the SBFEM in the commercial FEM software Abaqus. Ya et al. [31] implemented an open-source polyhedral SBFEM element for three-dimensional and nonlinear problems through the Abaqus UEL. Ye et al. [32] implemented a polygon SBFEM for two-dimensional linear elastostatic problems using the UEL subroutine. Yang et al. [33] developed a polygon SBFEM element to solve steady-state and transient heat conduction problems. These studies showed that the Abaqus UEL subroutine is suitable for implementing the SBFEM.

Yang et al. [34] conducted an explicit dynamic analysis using the VUEL subroutine. For high-speed simulations, the advantages of the explicit dynamic analysis are apparent within the desired tolerance. The cost of the explicit method is much less than that of the implicit method. However, the solutions are more unstable for low-speed simulations because high-frequency numerical noise becomes more important. The numerical damping is induced in the implicit analysis functions to remove noise and maintain accuracy. The explicit method is conditionally stable, and the stable time period is thus much shorter than that for the implicit method. Hence, an explicit method has many time increments for a long overall procedure [35]. At present, the use of the Abaqus user-defined element of the SBFEM mainly focuses on the forced vibration problem [34] whereas the free vibration problem has received little attention.

Two-dimensional models are relatively easy to set up and have reasonably short analysis times in engineering design, allowing sensitivity and optimization analyses [36]. Additionally, the models are readily available. Andersen and Jones [37] compared a two-dimensional model and a three-dimensional model for vibration analysis of two railway tunnel structures using the combined finite element method and boundary element method. They observed that the results obtained using the two-dimensional model are qualitatively consistent with those obtained using the three-dimensional model at most frequencies. It is desirable to use two-dimensional models as much as possible owing to these design advantages.

This work aims to develop a UEL of the PSBFEM for free and forced vibration analyses by relying on the user subroutine interface of general finite element software Abaqus. The remainder of the paper is organized as follows. Section 2 introduces the basic principles of the SBFEM in dynamic analysis. Section 3 describes the implementation of the PSBFEM with dynamic analysis. Section 4 compares the convergence and accuracy rates of the PSBFEM with those of the conventional FEM using several numerical benchmarks. Section 5 presents concluding remarks.

#### 2. Theory

##### 2.1. Polygonal SBFEM Formulation

The motion equation and boundary conditions are written aswhere is the differential operator, is the Cauchy stress tensor, is the body force, is the mass density, is the displacement vector, and is the acceleration vector. denotes the displacement boundary conditions, and denotes the surface traction boundary conditions.

The SBFEM has a local coordinate system as illustrated in Figure 1. The coordinates of a point along the radial line and inside the domain are expressed as [38].where is the radial coordinate and is the circumferential coordinate.

The differential operator is transformed from the Cartesian coordinate system to the scaled boundary coordinate system bywithwhere the Jacobian matrix at the boundary () is expressed as

The displacement field at any point in the SBFEM coordinates is expressed aswhere is the radial displacement function along a line connecting the scaling center and a node at the boundary and is the shape function matrix:

The strain-displacement transition matrices and are introduced:

The strain field is expressed in the scaled boundary coordinates as

The stress field is expressed aswhere is an elasticity matrix that has the formfor plane stress problems, where is Young’s modulus and is Poisson’s ratio.

According to the principle of virtual work, the radial displacement function is the solution to the SBFEM equation in displacement [39]:

In the derivation, the governing equations of linear elasticity are weakened in the circumferential direction whereas they remain strong in the radial direction. The coefficient matrices depend only on the geometry and material properties and are expressed as

When , the solution of equation (12) is a series of second-order ordinary differential equations obtained by introducing the variablewhere is the internal force vector. Equation (12) can be rewritten as a first-order ordinary differential equation with twice the number of unknowns aswhere the coefficient matrix is a Hamiltonian matrix. The solution for a bounded domain is obtained using the positive eigenvalues of . Hence, is expressed as

The solution of equation (15) can be obtained using the eigenvalue decomposition technique. The eigenvalue decomposition of is expressed as

The real parts of eigenvalues are negative whereas those of are positive. and are transformation matrices corresponding to the modal displacements and forces, respectively. The general solution of equation (15) is written aswhere and are the integration constants. For a super element with , the solution of equation (15) iswhere the integration constants are extracted from the nodal displacements on the boundary as

Eliminating the integration constants from at yields

Hence, the stiffness matrix of the S-element is written as

The equation of motion of an S-element is expressed aswhere is the mass matrix, is the damping matrix, and is the external force vector.

The mass matrix is written aswhere the coefficient matrix is introduced aswith

To perform the integration with respect to in equation (24) analytically, the abbreviationis introduced. The integration is rewritten asand equation (24) is rewritten as

Each entry of matrix is evaluated analytically, resulting in

##### 2.2. Natural Frequency and Mode Shapes

The equations of equilibrium governing the linear dynamic response are written as [40]where , , and are, respectively, the mass, damping, and stiffness matrices; is the external force vector; and , , and are, respectively, the acceleration, velocity, and displacement vectors. The relationship between the damped natural frequency and the undamped natural frequency is expressed as [41]

Equation (32) shows that the damped natural frequency is lower than the undamped natural frequency. However, in many practical cases, the damping ratio is relatively small and the difference is thus negligible [41, 42]. Equation (32) is therefore rewritten as

For the free vibrations, the external force vector is zero [43]. Hence,

For the steady-state conditions, starting from the equilibrium state, we setwhere is the vector of nodal amplitudes of vibration and (rad/s) is the circular frequency. Introducing equation (35) into equation (34) yields

We thus have a generalized eigenvalue problemwhere is an eigenvector, representing the vibrating mode, corresponding to the eigenvalue . The eigenvalue is the square of the circular frequency . The frequency in hertz is obtained as .

##### 2.3. Implicit Dynamic Analysis

Abaqus/Standard uses the Hilber–Hughes–Taylor (HHT) [44] time integration for the implicit dynamic analysis. The HHT method is an extension of the Newmark method [45].

In solving the dynamic equationthe HHT method is expressed as [44]where is time, is the time increment, and , , and are parameters that are determined to obtain integration accuracy and stability; , , and [46]. This method reduces to the Newmark method when the integration parameter is set to zero [47].

#### 3. Implementation Details

##### 3.1. Implementation of PSBFEM Using the UEL

According to the Abaqus User Subroutines Reference Manual [48], the most critical work of the UEL is to update the contribution of the element to the internal force vector RHS and the stiffness matrix AMATRX in the user subroutine interface provided by Abaqus. AMATRX is an array containing the contribution of the element of the stiffness or the other matrix of the overall system of equations whereas RHS is an array containing the contributions of the element to the right-hand-side vectors of the overall system of equations.

In free vibration analysis, the natural frequencies of an undamped system can be written as equation (37), and AMATRX is thus defined aswhere LFLAGS is an array containing flags that define the current solution procedure and requirements for element calculations.

In the case of forced vibration, Abaqus/Standard adopts implicit dynamic analysis, and AMATRX and RHS are defined as

Equation (41) can be rewritten as

Substituting equation (45) into equation (40) yields

It follows from equations (43), (45), and (46) that

Finally, the UEL subroutine updates AMATRX and RHS according to equations (47) and (44) in the forced vibration analysis.

Figure 2 shows an overall implementation of the UEL subroutine. Using the connectivity information of the input file, the UEL computes the scaling centers and transforms the global coordinates into local coordinates. Equation (13a) computes the coefficient matrices , , , and , which are used to construct the Hamilton matrix using (16). The two eigenvector matrices (, ) are constructed through eigenvalue decomposition. We finally obtain the stiffness matrix and mass matrix of the PSBFEM elements.

We need to adopt eigenvalue decomposition (see equation (17)) in solving the stiffness matrix and mass matrix. There are many mathematical libraries with which eigenvalue decomposition is performed. We use the Intel Math Kernel Library (MKL) [49] for eigenvalue decomposition in this work.

Abaqus CAE does not support the visualization of UEL results. Therefore, this paper uses the open-source program Paraview [50] to visualize the results. The specific method is to write the node information, field variable information, unit information, and unit type in the result file “odb” into the “vtu” file according to the file format provided by Paraview through a Python script and finally import the file into Paraview for visualization.

##### 3.2. Defining Elements in the UEL

Figure 3 shows the schematic diagram of the PSBFEM polygon mesh. Figure 3(a) shows the arbitrary polygon mesh, including triangle, quadrilateral, and pentagonal elements. Moreover, Figure 3(b) shows a quadtree mesh, including a traditional quadrilateral mesh and complex quadrilateral mesh. For complex quadrilateral elements, such as elements 5 and 6 in Figure 3(b), there are hanging nodes, which can be regarded as polygonal elements in the PSBFEM.

The input file of Abaqus usually contains a complete description of the numerical model, such as information on the nodes, elements, degrees of freedom of nodes, and materials. This information needs to be defined in the “*inp*” file by the user. We show a simple polygon mesh of the PSBFEM to demonstrate the defining of elements in the UEL, as shown in Figure 3(a). The mesh comprises three types of element: triangular elements (U3), quadrilateral elements (U4), and pentagon elements (U5). The pentagon element (U5) is defined as follows:(a)The keywords of basic information are shown as follows:(1)^{∗}User element, nodes = 5, type = *U*5, properties = 2, coordinates = 2(2)1,2(3)^{∗}Element, type = *U*5, elset = *E*5(4)3,2,3,4,8,7(5)^{∗}Uel property, elest = *E*5(6)2e5,0.3,2000,0,0 Here, Line 1 assigns the element type, number of nodes, number of element properties, and number of degrees of freedom at each node; Line 2 assigns the active degrees of freedom; Lines 3 and 4 define the element set E5; and Lines 5 and 6 set the element properties (i.e., Young’s modulus, Poisson’s ratio, density, and Rayleigh damping coefficient) of E5.(b)The keywords of free vibration analysis are shown as follows:(7)^{∗}Frequency, eigensolver = Lanczos(8)5, , , , , Here, Line 7 defines the analysis step of eigenvalue extraction; the eigenvalue extraction method is Lanczos. Line 8 defines the number of eigenvalues requested.(c)The keywords of forced vibration analysis are shown as follows:(9)^{∗}Dynamic, direct(10)0.02,25.

Here, Lines 9 and 10 define the implicit dynamic analysis of the fixed time increment; the total time and the time increment .

#### 4. Numerical Examples

This section validates the convergence and accuracy of the implementation by solving benchmark problems. Moreover, the results of the PSBFEM are compared with those of the FEM. The FEM analysis uses the commercial finite element software Abaqus. The comparison is evaluated on an Intel Core i7-4710MQ CPU running at 2.50 GHz with 4.0 GB of RAM. The relative error norm in the results is analyzed according towhere is the numerical result obtained using the PSBFEM and is the reference solution obtained as the analytical solution or numerical solution.

##### 4.1. Analysis of the Natural Frequency and Mode Shapes

###### 4.1.1. Free Vibration Analysis of a Horizontal Soil Foundation

In this example, we consider a horizontal soil foundation with a thickness of 5 m for free vibration analysis, as shown in Figure 4. The material properties are a dynamic shear modulus of soil = 1000 kPa, Poisson’s ratio = 0.3, and mass density = 2000 kg/m^{3}. The bottom of the foundation is fixed for all displacement components. The vertical displacement of the soil foundation is constrained. The frequency of mode of the foundation is expressed as [51]

The shear wave velocity is calculated as

The domain is discretized with quadrangle and arbitrary polygonal elements. A convergence study is conducted using both types of mesh. We refine meshes successively following the sequence = 0.5, 0.25, 0.1, and 0.05 m. Element meshes with element size of 0.5 m are shown in Figure 4(b) and 4(c).

The relative errors in the first five natural frequencies are given in Table 1. The errors decrease as the mesh is refined. The relative error of the PSBFEM is less than that of the FEM. As an example, the error of the PSBFEM is 1.057% whereas that of the FEM is 1.623% when the element size is 0.5 m. Moreover, the accuracy of the PSBFEM improves with the refining of the mesh. As an example, the error of the PSBFEM is 0.007% whereas that of the FEM is 0.017% when the element size is 0.05 m. The PSBFEM element is thus more accurate than the FEM element for the same element size.

Figure 5 shows that the relative error decreases at natural frequencies with an increasing number of degrees of freedom. The convergence rate of the PSBFEM is close to that of the FEM. Figure 6 shows the first four mode shapes of the foundation, where the displacement field is normalized. Mode shapes are virtually the same for the FEM and PSBFEM.

###### 4.1.2. Free Vibration Analysis of a Multihole Panel

A panel with four circular holes, as shown in Figure 7, is considered to demonstrate the flexibility of the PSBFEM element and quadtree algorithm. The multihole panel has a height = 1 m and length = 5 m. Each of the four circular holes has a radius of 0.3 m. The material properties are a Young’s modulus *E* = 206 GPa, mass density , and Poisson’s ratio .

The quadtree mesh is generated by setting the same number of mesh seeds on each hole, as shown in Figure 8(b), and the mesh transition between the holes of different sizes is effectively handled. The multihole panel is modeled with 1632 quadtree elements, and there are a total of 2478 nodes. Moreover, this problem is analyzed with a similar number of nodes (2565 nodes) using the Abaqus CPS4 element.

**(a)**

**(b)**

Table 2 gives the relative error norms of the FEM and PSBFEM as 0.54% and 0.80%, respectively. The two methods thus show good accuracy. However, the FEM is more accurate in this example owing to the meshes of the FEM conforming more than the quadtree meshes to the geometry. The first four mode shapes are presented in Figure 9. Modes 1, 2, and 4 are characterized by bending deformation and mode 3 by axial deformation. The present solutions are in good agreement with those obtained by the FEM.

**(a)**

**(b)**

##### 4.2. Response History Analysis

###### 4.2.1. Forced Vibration Analysis of a Cantilever Beam

To verify the validity of the proposed method for forced vibration problems, we investigate the forced vibration of a cantilever beam, as shown in Figure 10(a). The mesh is refined through *h*-refinement. The meshes of the FEM and PSBFEM are shown in Figures 10(c) and 10(d). The results are compared with the analytical solutions [52]. The cantilever beam has a length of 1.0 m and height of 0.5 m. The material properties are a Young’s modulus *E* = 100 Pa, mass density = 1.0 kg/m^{3}, and Poisson’s ratio = 0. A uniform loading is applied to the right edge with a Heaviside step loading = 1.0 , as shown in Figure 10(b).

**(a)**

**(b)**

**(c)**

**(d)**

In general, the numerical integrations using the Newmark method are accurate when is smaller than approximately 0.01 [40]. In this example, we consider that the time step *s* is used for the time integration, and the total time is 1.2 s. The convergence of the relative error in the horizontal displacement is presented in Figure 11. All methods show an optimal convergence rate. The PSBFEM element provides better accuracy than the FEM element. Figure 12 illustrates the horizontal displacement time history of the monitoring point. The PSBFEM obtains a history of acceleration that is in excellent agreement with that obtained by the FEM.

###### 4.2.2. Forced Vibration Analysis of a Multihole Plate

We again use the model of a multihole panel. The geometry and boundary conditions are presented in Figure 13. The same meshes as used in Section 4.1.2 are used in this example, as shown in Figure 8. As for the results in Section 4.1.2, the fundamental frequency is . Hence, a concentrated force is applied to the right end of the plate. The material properties are the same as those in Section 4.1.2. The response analyses of damped and undamped systems are considered. The relationship between the coefficients of the Rayleigh damping and and the damping ratio is expressed as [40]where denotes the th mode. In traditional methods, two reference vibration modes are selected, and their damping ratios and are obtained through measurement or reliable test data estimation and their frequencies and are used to solve and [53]:

When , the coefficients of Rayleigh damping are determined as

Generally, the damping ratio is 0.03–0.05 [51]. In this example, we consider that the damping ratio is 0.05. Hence, the coefficients of Rayleigh damping are and . The time step *s* is used for the time integration, and the total time is 0.4 s. The vertical displacement and velocity time history of the monitoring point calculated using the FEM and PSBFEM are compared in Figures 14 and 15. The results obtained using the two methods correspond well.

##### 4.3. Case Analysis of a Concrete-Face Rockfill Dam

We consider the vibration analysis of a simple concrete-face rockfill dam (CFRD) to show the advantage of the PSBFEM in handling hanging nodes. The geometry of the CFRD is shown in Figure 16. The height of the dam is 10 m, and the slopes upstream and downstream have a gradient of 1 : 2. The thickness of the concrete face is 0.4 m. The properties of the CFRD are given in Table 3. The bottom boundary of the dam is fixed for all displacement components in the free vibration analysis. In the forced vibration analysis, the vertical displacement of the bottom boundary of the dam is fixed. The bottom boundary of the dam is subject to a seismic load [51] in the horizontal direction.

The accuracy of modeling improves with an increasing number of degrees of freedom as shown in Sections 4.1 and 4.2. However, increasing the number of degrees of freedom requires a longer computing time. We should thus choose an appropriate mesh resolution that ensures accuracy but a low cost of calculation. Figure 17 presents the results of the mesh size sensitivity for two element types. It is clear that these elements provide more accurate results when the mesh resolution is less than 0.5 m. The size of the concrete face of the CFRD is smaller than the size of the rockfill. In this case, to well evaluate the dynamic characteristics of the face, we set the mesh resolution of the face at 0.06 m. There are three modeling cases for this example. Case 1 is FEM modeling with mesh resolutions of the face and rockfill of 0.06 m, as shown in Figure 18(a). Case 2 is FEM modeling with mesh resolutions of the face and rockfill of 0.06 and 0.5 m, respectively, as shown in Figure 18(b). We set a tie constraint between the concrete face and rockfill to compatibly couple the meshes. Case 3 is PSBFEM modeling with mesh resolutions of the face and rockfill again set at 0.06 and 0.5 m, respectively, as shown in Figure 18(c). The element type is a plane strain element.

**(a)**

**(b)**

**(c)**

Table 4 gives the first five natural frequencies of the CFRD. The relative error norms of the three cases for the first five natural frequencies are 0.03%, 0.18%, and 0.04%. The PSBFEM and the FEM with a global element size of 0.06 m is more accurate than the FEM using the tie constraint. Figure 19 shows the first three mode shapes of the dam. The mode shapes are virtually the same for the FEM and PSBFEM. In addition, the PSBFEM obtains a history of acceleration that is in good agreement with the reference solution, as shown in Figure 20. The total CPU running time is 22.5 s for Case 1, 0.6 s for Case 2, and 1.4 s for Case 3. Results show that the computational cost of the PSBFEM is slightly higher than that when using the Abaqus standard element for a similar number of elements. The PSBFEM is bound to generate relatively more computations than the standard FEM owing to the PSBFEM needing to decompose the eigenvalue and sort eigenvalues and conduct other operations [54]. However, the PSBFEM avoids the occurrence of hanging nodes by constructing a polygonal mesh. Hence, the PSBFEM can choose an appropriate mesh resolution for different structures to ensure accuracy and reduce the cost of calculation.

#### 5. Conclusions

This work implemented the PSBFEM in two-dimensional free and forced vibration analyses within an Abaqus/Standard analysis using the UEL subroutine. The work mainly focused on the main procedures of interacting with Abaqus, updating AMATRX and RHS, defining the UEL element in the input file, and solving the stiffness matrix and mass matrix through eigenvalue decomposition using the MKL mathematical libraries in the UEL implementation procedure. Moreover, meshes were generated automatically and results were visualized in Paraview using a Python script.

The implementation of the PSBFEM was validated against the FEM by solving benchmark problems. The results demonstrate that the PSBFEM is more accurate than the FEM with mesh refinement. In addition, the implementation of the PSBFEM can conveniently use arbitrary polygon elements through polygon/quadtree discretization in the commercial finite element software Abaqus. Notably, because the PSBFEM avoids hanging nodes by constructing a polygonal mesh, the PSBFEM can choose an appropriate mesh resolution for different structures to ensure accuracy and reduce calculation costs.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request. The developed UEL source code and the associated input files can be downloaded from https://github.com/hhupde/DynamicPSBFEM.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the National Natural Science Foundation of China (grant no. 51579089) and the Yunnan Postdoctoral Orientation Project. The authors thank Liwen Bianji (Edanz) (http://www.liwenbianji.cn) for editing the language of a draft of this manuscript.