# Free and Forced Vibration Analysis in Abaqus Based on the Polygonal Scaled Boundary Finite Element Method

Authors: Nan Ye; Chao Su; Yang Yang
Journal: Advances in Civil Engineering
 (2021)

---

## 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.


---

## Body

## 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. Section2 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 as(1)∇σ+b=ρu..inΩ,u=u−inΓu,σ⋅n=t−inΓt,where ∇ is the differential operator, σ is the Cauchy stress tensor, b is the body force, ρ is the mass density, u is the displacement vector, and u.. is the acceleration vector. Γu denotes the displacement boundary conditions, and Γt denotes the surface traction boundary conditions.The SBFEM has a local coordinate systemξ,η as illustrated in Figure 1. The coordinates of a point x,y along the radial line and inside the domain are expressed as [38].(2a)x=ξNηx,(2b)y=ξNηy,where ξ is the radial coordinate and η is the circumferential coordinate.Figure 1
Scaled boundary finite element coordinate system on a generic polygon.The differential operator is transformed from the Cartesian coordinate system to the scaled boundary coordinate system by(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) is expressed as(5)Jb=xbybxb,ηyb,η=xbyb,η−ybxb,η.The displacement fielduξ,η at any point in the SBFEM coordinates is expressed as(6)uξ,η=Nuηuξ.where uξ is the radial displacement function along a line connecting the scaling center O and a node at the boundary and 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σ is expressed as(10)σ=DB1uξ,ξ+1ξB2uξ,where D is an elasticity matrix that has the form(11)D=E1−v21v0v10001−v2,for plane stress problems, where E is Young’s modulus and v is Poisson’s ratio.According to the principle of virtual work, the radial displacement functionuξ is the solution to the SBFEM equation in displacement [39]:(12)E0ξ2uξ,ξξ+E0+E1T−E1ξuξ,ξ−E2uξ+Fξ=0.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 matricesEi,i=0,1,2 depend only on the geometry and material properties and are expressed as(13a)E0=∫−1+1B1TDB1Jbdη,(13b)E1=∫−1+1B2TDB1Jbdη,(13c)E2=∫−1+1B2TDB2Jbdη,(13d)M0=∫−1+1NuTρNuJbdη.WhenFξ=0, the solution of equation (12) is a series of second-order ordinary differential equations obtained by introducing the variable(14)Xξ=uξqξ,where qξ is the internal force vector. Equation (12) can be rewritten as a first-order ordinary differential equation with twice the number of unknowns 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 is expressed as(16)Zp=−E0−1E1TE0−1E2−E1E0−1E1TE1E0−1.The solution of equation (15) can be obtained using the 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 whereas those of λp are positive. Φun and Φup are transformation matrices corresponding to the modal displacements and forces, respectively. The general solution of equation (15) is 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,qξ=Φqnξ−λncn,where the integration constants cn are extracted from the nodal displacements on the boundary ub=uξ=1 as(20)cn=Φun−1ub.Eliminating the integration constantscn from at ξ=1 yields(21)F=qξ=1=ΦqnΦun−1ub.Hence, the stiffness matrix of the S-element is written as(22)K=ΦqnΦun−1.The equation of motion of an S-element is expressed as(23)Mu¨+Cu˙+Ku=P,where M is the mass matrix, C is the damping matrix, and P is the external force vector.The mass matrix is written as(24)M=Φbu−T∫01ξλbΦbuTM0ΦbuξλbξdξΦbu−1,where the coefficient matrix M0 is introduced as(25)M0=∑eM0e,with(26)M0e=∫−1+1NuTρNuJbdη.To perform the integration with respect toξ in equation (24) analytically, the abbreviation(27)m0=ΦbuTM0Φbu,is introduced. The integration is rewritten as(28)m=∫01ξλbm0ξλbξdξ,and equation (24) is rewritten as(29)M=Φbu−TmΦbu−1.Each entry of matrixm is evaluated analytically, resulting in(30)mij=∫01ξλbim0ijξλbjξdξ=m0ijλbi+λbj+2.


### 2.2. Natural Frequency and Mode Shapes

The equations of equilibrium governing the linear dynamic response are written as [40](31)Mu¨+Cu˙+Ku=P,where M, C, and K are, respectively, the mass, damping, and stiffness matrices; P is the external force vector; and u¨, u˙, and u are, respectively, the acceleration, velocity, and displacement vectors. The relationship between the damped natural frequency fd and the undamped natural frequency fud is expressed as [41](32)fd=1−ζ2fud.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(33)Mu¨+Ku=P.For the free vibrations, the external force vectorP is zero [43]. Hence,(34)Mu¨+Ku=0.For the steady-state conditions, starting from the equilibrium state, we set(35)u=Usinωt,where U is the vector of nodal amplitudes of vibration and ω (rad/s) is the circular frequency. Introducing equation (35) into equation (34) yields(36)KU=ω2MU.We thus have a generalized eigenvalue problem(37)KU=λMU,where U is an eigenvector, representing the vibrating mode, corresponding to the eigenvalue λ. The eigenvalue λ is the square of the circular frequency ω. The frequency f in hertz is obtained as f=ω/2π.


### 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 equation(38)Mu¨+Cu˙+Ku=P,the HHT method is expressed as [44](39)Mu¨t+Δt+1+αCu˙t+Δt+Kut+Δt−αCu˙t+Kut=Pt+Δt,(40)u˙t+Δt=u˙t+1−γu¨t+γu¨t+ΔtΔt,(41)ut+Δt=ut+u˙tΔt+12−βu¨t+βu¨t+ΔtΔt,where t is time, Δt is the time increment, and α, β, and γ are parameters that are determined to obtain integration accuracy and stability; −1/3≤α<0, β≥1/4γ+1/22, and γ≥1/2 [46]. This method reduces to the Newmark method when the integration parameter α is set to zero [47].


## 2.1. Polygonal SBFEM Formulation

The motion equation and boundary conditions are written as(1)∇σ+b=ρu..inΩ,u=u−inΓu,σ⋅n=t−inΓt,where ∇ is the differential operator, σ is the Cauchy stress tensor, b is the body force, ρ is the mass density, u is the displacement vector, and u.. is the acceleration vector. Γu denotes the displacement boundary conditions, and Γt denotes the surface traction boundary conditions.The SBFEM has a local coordinate systemξ,η as illustrated in Figure 1. The coordinates of a point x,y along the radial line and inside the domain are expressed as [38].(2a)x=ξNηx,(2b)y=ξNηy,where ξ is the radial coordinate and η is the circumferential coordinate.Figure 1
Scaled boundary finite element coordinate system on a generic polygon.The differential operator is transformed from the Cartesian coordinate system to the scaled boundary coordinate system by(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) is expressed as(5)Jb=xbybxb,ηyb,η=xbyb,η−ybxb,η.The displacement fielduξ,η at any point in the SBFEM coordinates is expressed as(6)uξ,η=Nuηuξ.where uξ is the radial displacement function along a line connecting the scaling center O and a node at the boundary and 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σ is expressed as(10)σ=DB1uξ,ξ+1ξB2uξ,where D is an elasticity matrix that has the form(11)D=E1−v21v0v10001−v2,for plane stress problems, where E is Young’s modulus and v is Poisson’s ratio.According to the principle of virtual work, the radial displacement functionuξ is the solution to the SBFEM equation in displacement [39]:(12)E0ξ2uξ,ξξ+E0+E1T−E1ξuξ,ξ−E2uξ+Fξ=0.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 matricesEi,i=0,1,2 depend only on the geometry and material properties and are expressed as(13a)E0=∫−1+1B1TDB1Jbdη,(13b)E1=∫−1+1B2TDB1Jbdη,(13c)E2=∫−1+1B2TDB2Jbdη,(13d)M0=∫−1+1NuTρNuJbdη.WhenFξ=0, the solution of equation (12) is a series of second-order ordinary differential equations obtained by introducing the variable(14)Xξ=uξqξ,where qξ is the internal force vector. Equation (12) can be rewritten as a first-order ordinary differential equation with twice the number of unknowns 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 is expressed as(16)Zp=−E0−1E1TE0−1E2−E1E0−1E1TE1E0−1.The solution of equation (15) can be obtained using the 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 whereas those of λp are positive. Φun and Φup are transformation matrices corresponding to the modal displacements and forces, respectively. The general solution of equation (15) is 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,qξ=Φqnξ−λncn,where the integration constants cn are extracted from the nodal displacements on the boundary ub=uξ=1 as(20)cn=Φun−1ub.Eliminating the integration constantscn from at ξ=1 yields(21)F=qξ=1=ΦqnΦun−1ub.Hence, the stiffness matrix of the S-element is written as(22)K=ΦqnΦun−1.The equation of motion of an S-element is expressed as(23)Mu¨+Cu˙+Ku=P,where M is the mass matrix, C is the damping matrix, and P is the external force vector.The mass matrix is written as(24)M=Φbu−T∫01ξλbΦbuTM0ΦbuξλbξdξΦbu−1,where the coefficient matrix M0 is introduced as(25)M0=∑eM0e,with(26)M0e=∫−1+1NuTρNuJbdη.To perform the integration with respect toξ in equation (24) analytically, the abbreviation(27)m0=ΦbuTM0Φbu,is introduced. The integration is rewritten as(28)m=∫01ξλbm0ξλbξdξ,and equation (24) is rewritten as(29)M=Φbu−TmΦbu−1.Each entry of matrixm is evaluated analytically, resulting in(30)mij=∫01ξλbim0ijξλbjξdξ=m0ijλbi+λbj+2.


## 2.2. Natural Frequency and Mode Shapes

The equations of equilibrium governing the linear dynamic response are written as [40](31)Mu¨+Cu˙+Ku=P,where M, C, and K are, respectively, the mass, damping, and stiffness matrices; P is the external force vector; and u¨, u˙, and u are, respectively, the acceleration, velocity, and displacement vectors. The relationship between the damped natural frequency fd and the undamped natural frequency fud is expressed as [41](32)fd=1−ζ2fud.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(33)Mu¨+Ku=P.For the free vibrations, the external force vectorP is zero [43]. Hence,(34)Mu¨+Ku=0.For the steady-state conditions, starting from the equilibrium state, we set(35)u=Usinωt,where U is the vector of nodal amplitudes of vibration and ω (rad/s) is the circular frequency. Introducing equation (35) into equation (34) yields(36)KU=ω2MU.We thus have a generalized eigenvalue problem(37)KU=λMU,where U is an eigenvector, representing the vibrating mode, corresponding to the eigenvalue λ. The eigenvalue λ is the square of the circular frequency ω. The frequency f in hertz is obtained as f=ω/2π.


## 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 equation(38)Mu¨+Cu˙+Ku=P,the HHT method is expressed as [44](39)Mu¨t+Δt+1+αCu˙t+Δt+Kut+Δt−αCu˙t+Kut=Pt+Δt,(40)u˙t+Δt=u˙t+1−γu¨t+γu¨t+ΔtΔt,(41)ut+Δt=ut+u˙tΔt+12−βu¨t+βu¨t+ΔtΔt,where t is time, Δt is the time increment, and α, β, and γ are parameters that are determined to obtain integration accuracy and stability; −1/3≤α<0, β≥1/4γ+1/22, and γ≥1/2 [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 as(42)AMATRX=K,M,when LFLAGS3=2,when LFLAGS3=4,where 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(43)AMATRX=Mdu¨du+1+αCdu˙du+1+αK,(44)RHS=−Pt+Δt=−Mu¨t+Δt−1+αCu˙t+Δt+Kut+Δt+αCu˙t+Kut,Equation (41) can be rewritten as(45)u¨t+Δt=1βΔt2ut+Δt−1βΔt2ut−1βΔtu˙t−12β−1u¨t.Substituting equation (45) into equation (40) yields(46)u˙t+Δt=γβΔtut+Δt−γβΔtut−γβ−1u˙t−Δtγ2β−1u¨t.It follows from equations (43), (45), and (46) that(47)AMATRX=1βΔt2M+γβΔt1+αC+1+αK.Finally, the UEL subroutine updates AMATRX and RHS according to equations (47) and (44) in the forced vibration analysis.Figure2 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 E0, E1, E2, and M0, which are used to construct the Hamilton matrix Zp using (16). The two eigenvector matrices (Φqn, Φun) are constructed through eigenvalue decomposition. We finally obtain the stiffness matrix K and mass matrix M of the PSBFEM elements.Figure 2
Framework of implementing the UEL in free vibration and forced vibration analyses.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

Figure3 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.Figure 3
Schematic diagram of polygon meshes: (a) random polygon meshes and (b) quadtree meshes.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 = U5, properties = 2, coordinates = 2(2)
1,2(3)
∗Element, type = U5, elset = E5(4)
3,2,3,4,8,7(5)
∗Uel property, elest = E5(6)
2e5,0.3,2000,0,0Here, 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 timet=25 and the time increment Δt=0.02.


## 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 as(42)AMATRX=K,M,when LFLAGS3=2,when LFLAGS3=4,where 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(43)AMATRX=Mdu¨du+1+αCdu˙du+1+αK,(44)RHS=−Pt+Δt=−Mu¨t+Δt−1+αCu˙t+Δt+Kut+Δt+αCu˙t+Kut,Equation (41) can be rewritten as(45)u¨t+Δt=1βΔt2ut+Δt−1βΔt2ut−1βΔtu˙t−12β−1u¨t.Substituting equation (45) into equation (40) yields(46)u˙t+Δt=γβΔtut+Δt−γβΔtut−γβ−1u˙t−Δtγ2β−1u¨t.It follows from equations (43), (45), and (46) that(47)AMATRX=1βΔt2M+γβΔt1+αC+1+αK.Finally, the UEL subroutine updates AMATRX and RHS according to equations (47) and (44) in the forced vibration analysis.Figure2 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 E0, E1, E2, and M0, which are used to construct the Hamilton matrix Zp using (16). The two eigenvector matrices (Φqn, Φun) are constructed through eigenvalue decomposition. We finally obtain the stiffness matrix K and mass matrix M of the PSBFEM elements.Figure 2
Framework of implementing the UEL in free vibration and forced vibration analyses.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

Figure3 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.Figure 3
Schematic diagram of polygon meshes: (a) random polygon meshes and (b) quadtree meshes.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 = U5, properties = 2, coordinates = 2(2)
1,2(3)
∗Element, type = U5, elset = E5(4)
3,2,3,4,8,7(5)
∗Uel property, elest = E5(6)
2e5,0.3,2000,0,0Here, 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 timet=25 and the time increment Δt=0.02.


## 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 to(48)e=Unum−UrefL2UrefL2,where Unum is the numerical result obtained using the PSBFEM and Uref 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 Figure4. The material properties are a dynamic shear modulus of soil G  = 1000 kPa, Poisson’s ratio ν  = 0.3, and mass density ρ  = 2000 kg/m3. The bottom of the foundation is fixed for all displacement components. The vertical displacement of the soil foundation is constrained. The frequency of mode j of the foundation is expressed as [51](49)ωj=2j−1π2Hvs.Figure 4
Geometry and meshes of a horizontal soil foundation: (a) geometry of a horizontal soil foundation, (b) Abaqus CPE4 meshes, and (c) polygonal meshes.The shear wave velocity is calculated as(50)vs=Gρ.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 sequenceh  = 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 Table1. 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.Table 1
The first five natural frequencies (Hz) of a horizontal soil foundation.
MethodElement size (m)Degrees of freedomNatural frequencies (Hz)Relative error (%)Mode 1Mode 2Mode 3Mode 4Mode 5Analytical solution——0.34671.04011.73342.42683.1202—FEM0.51260.34661.03771.72232.39643.05561.6230.254100.34671.03951.73072.41923.10400.4070.114580.34671.04001.73302.42563.11760.0650.0584420.34671.04001.73332.42653.11950.017PSBFEM0.51640.34681.04161.74062.44663.16231.0570.256440.34671.04041.73492.43083.12860.2120.139980.34671.04011.73362.42733.12130.0270.05159980.34671.04011.73342.42683.12050.007Figure5 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.Figure 5
Convergence behavior of the soil foundation in terms of the relative error norm at the first five natural frequencies.Figure 6
First four mode shapes of the horizontal soil foundation. (a) FEM. (b) PSBFEM.


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

A panel with four circular holes, as shown in Figure7, is considered to demonstrate the flexibility of the PSBFEM element and quadtree algorithm. The multihole panel has a height H  = 1 m and length L  = 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 ρ=7800kg/m3, and Poisson’s ratio v=0.3.Figure 7
Geometry and boundary conditions of a multihole panel.The quadtree mesh is generated by setting the same number of mesh seeds on each hole, as shown in Figure8(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.Figure 8
Elements of a multihole panel: (a) Abaqus CPS4 elements and (b) quadtree elements.
(a)(b)Table2 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.Table 2
The first four natural frequencies (Hz) of a multihole panel.
MethodNatural frequenciesf (Hz)Relative error norm (%)Mode 1Mode 2Mode 3Mode 4Reference solutions31.811150.19209.31330.64—FEM31.916151.15210.03333.400.54PSBFEM31.992151.66210.07334.930.80Figure 9
The first four modes of a multihole panel: (a) modes of the FEM and (b) modes of the PSBFEM.
(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 Figure10(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/m3, and Poisson’s ratio ν  = 0. A uniform loading is applied to the right edge with a Heaviside step loading Pt  = 1.0 N/m2, as shown in Figure 10(b).Figure 10
Heaviside step loading applied to a cantilever beam: (a) geometry and boundary condition, (b) Heaviside step loading, (c) Abaqus CPS4 element, and (d) PSBFEM element.
(a)(b)(c)(d)In general, the numerical integrations using the Newmark method are accurate whenΔt/T is smaller than approximately 0.01 [40]. In this example, we consider that the time step Δt=0.001s 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.Figure 11
Convergence of the relative error in the horizontal displacement at the monitoring point.Figure 12
Time history of the horizontal displacement of the monitoring point.


#### 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 Figure13. 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 ω=2πf=199.77. Hence, a concentrated force pt=10000sin199.77πt 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 αR and βR and the damping ratio ςn is expressed as [40](51)ςn=αR2ωn+βRωn2,where n denotes the nth mode. In traditional methods, two reference vibration modes are selected, and their damping ratios ςi and ςj are obtained through measurement or reliable test data estimation and their frequencies ωi and ωj are used to solve αR and βR [53]:(52)αRβR=2ωiωjωj2−ωi2ωj−ωi−1ωj1ωiςiςj.Figure 13
Geometry and boundary conditions of a multihole panel subject to a concentrated dynamic loading.Whenςi=ςj=ς, the coefficients of Rayleigh damping are determined as(53)αRβR=2ςωi+ωjωiωj1.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 αR=16.494 and βR=0.0001. The time step Δt=0.001s 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.Figure 14
Time history of the vertical displacement.Figure 15
Time history of the vertical velocity.


### 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 Figure16. 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.Figure 16
Geometry model of a simple CFRD.Table 3
Properties of the CFRD.
MaterialYoung’s modulus (MPa)Poisson’s ratioDensity (kg/m3)Concrete face3.0 × 1040.1672500Rockfill5200.32000The accuracy of modeling improves with an increasing number of degrees of freedom as shown in Sections4.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.Figure 17
Sensitivity analysis of meshes.Figure 18
Mesh model of the CFRD: (a) Case 1, FEM model with a global element size of 0.06 m; (b) Case 2, FEM model with a tie constraint; and (c) Case 3: PSBFEM model.
(a)(b)(c)Table4 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.Table 4
The first five natural frequencies (Hz) of a simple CFRD.
CaseNumber of elementsNatural frequencies (Hz)Relative error norm (%)Time consumption (s)Mode 1Mode 2Mode 3Mode 4Mode 5Reference—12.41316.67320.49222.10625.580——Case 16932912.40316.67020.48222.10325.5770.0322.5Case 2221212.39516.65820.44722.07125.5110.180.6Case 3278812.40616.67620.49822.11225.6020.041.4Figure 19
The first three modes of the dam in the three cases.Figure 20
History of acceleration at the top of the dam.


## 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 Figure4. The material properties are a dynamic shear modulus of soil G  = 1000 kPa, Poisson’s ratio ν  = 0.3, and mass density ρ  = 2000 kg/m3. The bottom of the foundation is fixed for all displacement components. The vertical displacement of the soil foundation is constrained. The frequency of mode j of the foundation is expressed as [51](49)ωj=2j−1π2Hvs.Figure 4
Geometry and meshes of a horizontal soil foundation: (a) geometry of a horizontal soil foundation, (b) Abaqus CPE4 meshes, and (c) polygonal meshes.The shear wave velocity is calculated as(50)vs=Gρ.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 sequenceh  = 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 Table1. 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.Table 1
The first five natural frequencies (Hz) of a horizontal soil foundation.
MethodElement size (m)Degrees of freedomNatural frequencies (Hz)Relative error (%)Mode 1Mode 2Mode 3Mode 4Mode 5Analytical solution——0.34671.04011.73342.42683.1202—FEM0.51260.34661.03771.72232.39643.05561.6230.254100.34671.03951.73072.41923.10400.4070.114580.34671.04001.73302.42563.11760.0650.0584420.34671.04001.73332.42653.11950.017PSBFEM0.51640.34681.04161.74062.44663.16231.0570.256440.34671.04041.73492.43083.12860.2120.139980.34671.04011.73362.42733.12130.0270.05159980.34671.04011.73342.42683.12050.007Figure5 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.Figure 5
Convergence behavior of the soil foundation in terms of the relative error norm at the first five natural frequencies.Figure 6
First four mode shapes of the horizontal soil foundation. (a) FEM. (b) PSBFEM.


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

A panel with four circular holes, as shown in Figure7, is considered to demonstrate the flexibility of the PSBFEM element and quadtree algorithm. The multihole panel has a height H  = 1 m and length L  = 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 ρ=7800kg/m3, and Poisson’s ratio v=0.3.Figure 7
Geometry and boundary conditions of a multihole panel.The quadtree mesh is generated by setting the same number of mesh seeds on each hole, as shown in Figure8(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.Figure 8
Elements of a multihole panel: (a) Abaqus CPS4 elements and (b) quadtree elements.
(a)(b)Table2 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.Table 2
The first four natural frequencies (Hz) of a multihole panel.
MethodNatural frequenciesf (Hz)Relative error norm (%)Mode 1Mode 2Mode 3Mode 4Reference solutions31.811150.19209.31330.64—FEM31.916151.15210.03333.400.54PSBFEM31.992151.66210.07334.930.80Figure 9
The first four modes of a multihole panel: (a) modes of the FEM and (b) modes of the PSBFEM.
(a)(b)


## 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 Figure4. The material properties are a dynamic shear modulus of soil G  = 1000 kPa, Poisson’s ratio ν  = 0.3, and mass density ρ  = 2000 kg/m3. The bottom of the foundation is fixed for all displacement components. The vertical displacement of the soil foundation is constrained. The frequency of mode j of the foundation is expressed as [51](49)ωj=2j−1π2Hvs.Figure 4
Geometry and meshes of a horizontal soil foundation: (a) geometry of a horizontal soil foundation, (b) Abaqus CPE4 meshes, and (c) polygonal meshes.The shear wave velocity is calculated as(50)vs=Gρ.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 sequenceh  = 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 Table1. 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.Table 1
The first five natural frequencies (Hz) of a horizontal soil foundation.
MethodElement size (m)Degrees of freedomNatural frequencies (Hz)Relative error (%)Mode 1Mode 2Mode 3Mode 4Mode 5Analytical solution——0.34671.04011.73342.42683.1202—FEM0.51260.34661.03771.72232.39643.05561.6230.254100.34671.03951.73072.41923.10400.4070.114580.34671.04001.73302.42563.11760.0650.0584420.34671.04001.73332.42653.11950.017PSBFEM0.51640.34681.04161.74062.44663.16231.0570.256440.34671.04041.73492.43083.12860.2120.139980.34671.04011.73362.42733.12130.0270.05159980.34671.04011.73342.42683.12050.007Figure5 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.Figure 5
Convergence behavior of the soil foundation in terms of the relative error norm at the first five natural frequencies.Figure 6
First four mode shapes of the horizontal soil foundation. (a) FEM. (b) PSBFEM.


## 4.1.2. Free Vibration Analysis of a Multihole Panel

A panel with four circular holes, as shown in Figure7, is considered to demonstrate the flexibility of the PSBFEM element and quadtree algorithm. The multihole panel has a height H  = 1 m and length L  = 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 ρ=7800kg/m3, and Poisson’s ratio v=0.3.Figure 7
Geometry and boundary conditions of a multihole panel.The quadtree mesh is generated by setting the same number of mesh seeds on each hole, as shown in Figure8(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.Figure 8
Elements of a multihole panel: (a) Abaqus CPS4 elements and (b) quadtree elements.
(a)(b)Table2 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.Table 2
The first four natural frequencies (Hz) of a multihole panel.
MethodNatural frequenciesf (Hz)Relative error norm (%)Mode 1Mode 2Mode 3Mode 4Reference solutions31.811150.19209.31330.64—FEM31.916151.15210.03333.400.54PSBFEM31.992151.66210.07334.930.80Figure 9
The first four modes of a multihole panel: (a) modes of the FEM and (b) modes of the PSBFEM.
(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 Figure10(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/m3, and Poisson’s ratio ν  = 0. A uniform loading is applied to the right edge with a Heaviside step loading Pt  = 1.0 N/m2, as shown in Figure 10(b).Figure 10
Heaviside step loading applied to a cantilever beam: (a) geometry and boundary condition, (b) Heaviside step loading, (c) Abaqus CPS4 element, and (d) PSBFEM element.
(a)(b)(c)(d)In general, the numerical integrations using the Newmark method are accurate whenΔt/T is smaller than approximately 0.01 [40]. In this example, we consider that the time step Δt=0.001s 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.Figure 11
Convergence of the relative error in the horizontal displacement at the monitoring point.Figure 12
Time history of the horizontal displacement of the monitoring point.


### 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 Figure13. 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 ω=2πf=199.77. Hence, a concentrated force pt=10000sin199.77πt 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 αR and βR and the damping ratio ςn is expressed as [40](51)ςn=αR2ωn+βRωn2,where n denotes the nth mode. In traditional methods, two reference vibration modes are selected, and their damping ratios ςi and ςj are obtained through measurement or reliable test data estimation and their frequencies ωi and ωj are used to solve αR and βR [53]:(52)αRβR=2ωiωjωj2−ωi2ωj−ωi−1ωj1ωiςiςj.Figure 13
Geometry and boundary conditions of a multihole panel subject to a concentrated dynamic loading.Whenςi=ςj=ς, the coefficients of Rayleigh damping are determined as(53)αRβR=2ςωi+ωjωiωj1.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 αR=16.494 and βR=0.0001. The time step Δt=0.001s 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.Figure 14
Time history of the vertical displacement.Figure 15
Time history of the vertical velocity.


## 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 Figure10(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/m3, and Poisson’s ratio ν  = 0. A uniform loading is applied to the right edge with a Heaviside step loading Pt  = 1.0 N/m2, as shown in Figure 10(b).Figure 10
Heaviside step loading applied to a cantilever beam: (a) geometry and boundary condition, (b) Heaviside step loading, (c) Abaqus CPS4 element, and (d) PSBFEM element.
(a)(b)(c)(d)In general, the numerical integrations using the Newmark method are accurate whenΔt/T is smaller than approximately 0.01 [40]. In this example, we consider that the time step Δt=0.001s 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.Figure 11
Convergence of the relative error in the horizontal displacement at the monitoring point.Figure 12
Time history of the horizontal displacement of the monitoring point.


## 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 Figure13. 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 ω=2πf=199.77. Hence, a concentrated force pt=10000sin199.77πt 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 αR and βR and the damping ratio ςn is expressed as [40](51)ςn=αR2ωn+βRωn2,where n denotes the nth mode. In traditional methods, two reference vibration modes are selected, and their damping ratios ςi and ςj are obtained through measurement or reliable test data estimation and their frequencies ωi and ωj are used to solve αR and βR [53]:(52)αRβR=2ωiωjωj2−ωi2ωj−ωi−1ωj1ωiςiςj.Figure 13
Geometry and boundary conditions of a multihole panel subject to a concentrated dynamic loading.Whenςi=ςj=ς, the coefficients of Rayleigh damping are determined as(53)αRβR=2ςωi+ωjωiωj1.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 αR=16.494 and βR=0.0001. The time step Δt=0.001s 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.Figure 14
Time history of the vertical displacement.Figure 15
Time history of the vertical velocity.


## 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 Figure16. 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.Figure 16
Geometry model of a simple CFRD.Table 3
Properties of the CFRD.
MaterialYoung’s modulus (MPa)Poisson’s ratioDensity (kg/m3)Concrete face3.0 × 1040.1672500Rockfill5200.32000The accuracy of modeling improves with an increasing number of degrees of freedom as shown in Sections4.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.Figure 17
Sensitivity analysis of meshes.Figure 18
Mesh model of the CFRD: (a) Case 1, FEM model with a global element size of 0.06 m; (b) Case 2, FEM model with a tie constraint; and (c) Case 3: PSBFEM model.
(a)(b)(c)Table4 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.Table 4
The first five natural frequencies (Hz) of a simple CFRD.
CaseNumber of elementsNatural frequencies (Hz)Relative error norm (%)Time consumption (s)Mode 1Mode 2Mode 3Mode 4Mode 5Reference—12.41316.67320.49222.10625.580——Case 16932912.40316.67020.48222.10325.5770.0322.5Case 2221212.39516.65820.44722.07125.5110.180.6Case 3278812.40616.67620.49822.11225.6020.041.4Figure 19
The first three modes of the dam in the three cases.Figure 20
History of acceleration at the top of the dam.


## 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.



---

*Source: 7664870-2021-12-31.xml*