Nonlinear analysis for the polygonal element

As an important method for solving boundary value problems of differential equations, the finite element method (FEM) has been widely used in the fields of engineering and academic research. For two dimensional problems, the traditional finite element method mainly adopts triangular and quadrilateral elements, but the triangular element is constant strain element, its accuracy is low, the poor adaptability of quadrilateral element with complex geometry. The polygon element is more flexible and convenient in the discrete complex geometric model. Some interpolation functions of the polygon element were introduced. And some analysis was given. The numerical calculation accuracy and related features of different interpolation function were studied. The damage analysis for the koyna dam was given by using the polygonal element polygonal element of Wachspress interpolation function. The damage result is very similar to that by using Xfem, which shows the calculation accuracy of this method is very high.


Introduction
The finite element method (FEM) is commonly utilized in engineering and academic research fields. The polygonal element is more flexible and convenient in the discrete complex geometric model and it can easily simulate the mechanical properties of a given material (see [1,2]). The polygonal element has the geometric flexibility which makes it suitable for simulating crack growth. Thus, it can be utilized to solve problems related to fracture and large-deformation, frictionless, dynamic, contact-impact of arbitrary geometries and boundaries (see [3,4]). Furthermore, the polygonal element and Voronoi FEM have favorable accuracy in solving plastic problems with small deformation (see [5]).
In this study, we exploited the afore-mentioned advantages of the polygonal element for the numerical analysis of a practical engineering problem. We compared the results yielded by the various kinds of polygonal element analysis and ultimately suggested the kinds of the shape function and numerical integrations as well as the number of integral points of polygonal element.

Polygonal mesh generation
Here, we first used the pretreatment function in ANSYS software to divide the geometric region into triangular meshes, to divide the boundary into linear meshes and to divide the corner point into point mesh. The element and node information of the above mesh were then derived according to a certain format. Next, the centroid of the elements (such as triangular elements, linear elements and point elements) with common nodes were arrayed counterclockwise and added to the point list of polygons. Finally, the polygon points (as the nodes of polygonal element) were connected to form polygonal element. Figure 1 shows the process from a triangular mesh to a polygonal mesh.

Structure of polygon element interpolation function
Wachspress interpolation function, Laplace interpolation function and Mean value coordinates interpolation function are the focus of the polygon coordinates. The general expression for the polygonal interpolation function is defined as: where ωi is the weight function of the vertex Pi of the polygonal element. i=1,2,…,n. ωi are satisfied the following nature: non negative: Ni(x,y)≥0,(i=1,2,…,n), properties of unit decomposition: ul u l linear completeness: ul c c u l. A polygon with n edges, its vertices are P1,P2,…,Pn are arranged in a counter clockwise direction. For any point P∈Ω, connecting the point P and the various vertices of the polygon, the polygon is divided into n triangles. The interior angles of a triangle are shown in Figure 2. All angles were γi=∠PPi+1Pi,γi-1=∠PPiPi-1,βi=∠PPiPi+1,βi-1=∠PPi-1Pi, αi=∠PiPPi+1, αi-1=∠Pi-1PPi, ‖P -Pi‖ is the distance between point P and Pi; A represents three points around the area of a triangle (i=1,2,…,n).
where values of ω (1) i and ω (2) i are not present at the boundary, ω (3) i to meet the nature of Kronecker strictly, it is meaningful at the node. The weight function of Laplace interpolation shape function is: The weight function of the mean value coordinates interpolation function is

Integral scheme for numerical integration
The polygon is divided into a number of triangles, respectively, numerical integration is performed on each triangle, as shown in Figure 3. Gauss numerical integral formula is expressed as where Wn,j, Wn,k is the weight function in the Gauss integral, ξn,j, ξn,k is the Gauss integral point coordinate.

Concrete simulation
We use the elastic-plastic damage evolving model in software ABAQUS for concrete. Based on the different tensile and compressive properties of concrete, the damage model is , 0 where σ and ε are stress and strain, respectively. ε pl is plastic strain, 0 el D and d are elastic stiffness without damage and damage factor, respectively, d can be defined as where dt and dc are tensile and compressive damage factor, respectively, which can be obtained according to the uniaxial tensile and compressive experiment results. st and sc are the state function of the recovery of stiffness about the direction of stress. Figure 4 shows the process of damage under uniaxial cyclic load using the elasticplastic damage evolving model in software ABAQUS. When the tensile stress reaches the peak (point A), concrete cracks. Then loading to point B, the tensile stiffness reduces E= (1-dt)E0 (in which dt is the tensile damage factor, E and E0 are defined as the elastic modulus and the initial elastic modulus). Then unloading to point C, the tensile stiffness changes to (1-dt)E0 and the loading path is changed to BC. When the reversed uniaxial compression is applied to the concrete, if wc=0, it describes that there is no stiffness recovery under uniaxial compression(E=(1-dt)E0) and the loading path is CD; if wc=1, it describes that the material fully recovers the compressive stiffness under uniaxial compression(E=E0) and the loading path is CEF. When loading to point F, the reversed uniaxial tension is applied to the concrete, if wt=0, it describes that there is no stiffness recovery under uniaxial tension (E=(1-dt) (1-dc)E0) and the loading path is GH; if wt=1, it describes that the material fully recovers the tensile stiffness under uniaxial tension(E=(1dc)E0) and the loading path is GJ. The function of non-associated flow potential is where G is the non-associated flow potential, which is adopted according to Drucker-Prager hyperbolic function. σt0 is the peak of uniaxial tensile stress. ξ is the parameter of the offset, which is the speed of the non-associated flow potential G tends to asymptote (when ξ is closed to 0, the non-associated flow potential G is closed to the straight line). ψ is the dilatancy angle in p-q plan. p is the hydrostatic pressure. q is Mises equivalent stress.  The length of the square plate is 1 m in the schematic diagram shown in Figure 5, and the uniformly distributed load q = 1 Pa．The material constant Young's modulus E = 100 MPa, Poisson ratio v = 0.3, and the square plate is in the plane stress state. The square panel is divided into polygonal meshes, as shown in Figure 6.
As shown in Figure 7, the highest accuracy of Wachspress interpolation function, and the accuracy increases with the increase of the number of integral points. The precision of Laplace interpolation function is lowest, and its precision is not improved obviously with the increase of integral point. So for convex polygonal elements, finite element analysis is performed by using Wachspress interpolation function. As shown in Figure 8, with the increase of the number of polygonal elements, the relative error tends to be stable quickly.

Damage analysis for the Koyna dam
The damage analysis for the Koyna dam is given by using the polygonal element polygonal element of Wachspress interpolation function as shown in Figure 10. The density is 2500 kg/m 3 . The Young's modulus is 34.5 GPa. The Poisson's ratio is 0.2. The material parameters for the damage model are the dilatancy angle ψ=30°, the flow potential eccentricity α=0.10, the ratio of initial equi-biaxial compressive yield stress to initial uniaxial compressive yield stress ξ=1.18, the ratio of the second stress invariant on the tensile meridian γ=0.68, the other parameters of the damage model wc=1, wt=1. The relationship between the damage factor, the tensile strength and inelastic strain for the damage model is shown in Table 1. The seismic ground accelerations are shown in Figure 9.  As shown in Figure 11, the polygonal element of Wachspress interpolation function can describe the damage distribution well. And the damage result is very similar to that by using Xfem, which shows the calculation accuracy of this method is very high.

Conclusion
In this paper, we compared the results yielded by the various kinds of polygonal element analysis and ultimately suggested the kinds of the shape function and numerical integrations as well as the number of integral points of polygonal element. The results show the highest accuracy of Wachspress interpolation function, and the accuracy increases with the increase of the number of integral points. The precision of Laplace interpolation function is lowest, and its precision is not improved obviously with the increase of integral point. So for convex polygonal elements, finite element analysis is performed by using Wachspress interpolation function. With the increase of the number of polygonal elements, the relative error tends to be stable quickly.