Abstract
The node moving and multistage node enrichment adaptive refinement procedures are extended in mixed discrete least squares meshless (MDLSM) method for efficient analysis of elasticity problems. In the formulation of MDLSM method, mixed formulation is accepted to avoid second-order differentiation of shape functions and to obtain displacements and stresses simultaneously. In the refinement procedures, a robust error estimator based on the value of the least square residuals functional of the governing differential equations and its boundaries at nodal points is used which is inherently available from the MDLSM formulation and can efficiently identify the zones with higher numerical errors. The results are compared with the refinement procedures in the irreducible formulation of discrete least squares meshless (DLSM) method and show the accuracy and efficiency of the proposed procedures. Also, the comparison of the error norms and convergence rate show the fidelity of the proposed adaptive refinement procedures in the MDLSM method.
1. Introduction
Adaptivity needs appropriate numerical solutions of problems in order to describe high gradient regions and an anisotropic behavior of the solution. In an adaptive procedure, a good error estimator plays a very important role. The error estimation in numerical methods is obviously as old as the numerical computations themselves. The earliest paper by Richardson [1] proposed an error estimation procedure to use in finite difference method. This was followed in finite element method (FEM) [2]. A particular strength of the FEM is the well-developed theories of error estimation and adaptivity. Three h-refinement procedures, namely, mesh movement, mesh enrichment, and remeshing have been proposed for adaptivity [3]. In the mesh movement, the total number of nodes remains constant, but the location of the nodes can change in order to achieve a better overall distribution of the error. In the mesh enrichment, the original nodes hold fix and hierarchical nodes or simply more nodes add to the problem domain based on error distribution. In the remeshing, a completely new nodes is constructed based on the information acquired from the previous computation, and hence, it is required to implement a suitable node generator. on one hand, the mesh movement is more suitable than mesh enrichment because the problem scale remains constant, and on the other hand, its interpolations become too distort in the mesh-based methods [4].
In order to avoid these problems, an alternative approach, known as meshless methods (MMs), has been developed in recent decades to discretize a continuum body only by a finite number of nodes. In MMs the unknowns are interpolated from the nodal values that constitute the problem degrees-of-freedom. The main advantage of MMs is the fact that the interpolation accuracy is much less affected by the nodal distribution. Many meshless methods have been introduced since Gingold and Monaghan [5] proposed smoothed particle hydrodynamics (SPH) method. Nayroles et al. [6] implemented the diffuse element method (DEM). Belytschko et al. [7] presented the Element-Free Galerkin (EFG) method. Liu et al. [8] suggested the reproducing kernel particle method (RKPM). The other meshless methods that have been developed in recent years are the Finite Point (FP) method [9], the HP clouds method [10], the meshless local Petrov-Galerkin (MLPG) method [11], the local boundary integral equation (LBIE) method [12], the finite cloud (FC) method [13] and the discrete least squares meshless (DLSM) method [14].
Researchers used the advantages of MMs for developing efficient error estimate and adaptivity procedures. Rabczuk and Belytschko [15, 16] proposed an adaptive continuum/discrete crack approach for meshfree particle methods and also an adaptivity procedure for structured meshfree particle methods in 2D and 3D problems. Yoon et al. [17] worked on enriched meshfree collocation method with diffuse derivatives for elastic fracture. Zi et al. [18] investigated extended meshfree methods without branch enrichment for cohesive cracks. Bordas et al. [19] proposed enriched meshfree methods without asymptotic enrichment for 3D nonlinear fracture mechanics. Rabczuk and Samaniego [20] worked on discontinuous modelling of shear bands using adaptive meshfree methods. Zhuang et al. [21–23] investigated error control in the EFG method and adaptivity for structured meshfree particle methods in 2D and 3D problems. The DLSM method was extendted for error estimate and adaptivity in solid [24, 25] and fluid [26] problems.
Two different formulations, namely irreducible and mixed formulations have been introduced and used for the solution of engineering problems. With the mixed formulation, the continuity requirement decreases by one order compared to the irreducible formulation [27]. Use of mixed formula may result in an improved approximation, in particular, for the gradient variables, which in turn could result in higher accuracy than possible with the irreducible formulation [27]. In the standard mixed FEM, in order to obtain a coefficient matrix which leads to the system of equations with a unique and stable solution, the polynomial functions chosen for approximation of stresses and displacements must satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB or inf-sup) condition [28, 29]. The stability of mixed discretization does not allow FEM to choose independently the approximation spaces, so these spaces are restricted in the stability condition which is known as the LBB condition. However, the least squares approximation has the advantage that it does not require satisfying the LBB condition [30–32]. Hence, this advantage was used by Amani et al. [33] to implement a mixed meshless method named mixed discrete least squares Meshless (MDLSM) method which is formulated based on the least squares residuals functional of the governing partial differential equations of planar elasticity problem and its boundary conditions at the nodal points, and hence, it is stable and is not required to satisfy the LBB condition between the displacements and stresses approximations. Hence, the approximation spaces of the displacements and stresses can be choosen independently while they are obtained simuletanously.
In this paper, the MDLSM method is extended for the residual based error estimation and for the two types of adaptive refinement procedures. The node moving adaptive refinement procedure based on the spring analogy [24] and the node enrichment adaptive refinement procedure [25] are formulated and used in the MDLSM method for efficient analysis of the elasticity problems.
The present paper is organized as follows. Formulation of the mixed discrete least squares meshless method for solving the planar elasticity problems is given in Section 2. In Section 3, an error estimator based on the least square functional residuals is formulated for the MDLSM method to use in the node moving and node enrichment adaptive refinement procedures. In Section 4, we present some numerical benchmark examples which illustrate the proposed adaptive refinement process as well as the efficiency of the error estimator. Finally, some concluding remarks are addressed in Section 5.
2. Formulation of Mixed Discrete Least Squares Meshless Method for Elasticity
Consider the following two-dimensional linear elasticity problem with displacement and traction boundary conditions as follow: where is a bounded domain representing the region occupied by an elastic body, and , are the Lame constants which are defined as where is the Poisson ratio, is the Young modulus, and , are the displacement and traction boundaries, respectively. , , , and prescribed respectively the displacements and tractions in the and directions and , are direction cosines of the normal vector to the boundary.
By using the following definition of stresses in terms of the displacement components: we can rewrite (1) in term of stresses as
The compact form of (1) can be written by substituting (4)-(5) into the second-order problem of (1) in the form of where is a first-order differential operator defined as and is the vector of unknowns defined as and vector contains the forcing terms which has the form In (7), , , and are defined by the following matrices: The displacement and traction boundary conditions (2) can be written in terms of the unknown vector as where and are defined as follow: The plane elasticity problem is now defined as solving the first-order differential equation subjected only to the Dirichlet type boundary condition
The application of the proposed MDLSM method for solving problem of (13) starts with the definition of residuals as follows: where and are domain and boundary residuals, respectively.
Now the penalty approach is used to form the least square residuals functional which is defined as where is the total number of sampling (or collocation) points, is total number of domain sampling points, is total number of boundary sampling points, and the penalty coefficient is a positive scalar constant that must be large enough in order to impose the essential boundary condition with the desired accuracy. A note should be made here regarding the value of the penalty parameters. To impose the boundary conditions exactly, the penalty factor must be infinite, which is not possible in practical numerical analysis. Therefore, the boundary conditions could not be satisfied exactly but only approximately. In general, the use of a larger penalty factor will lead to better enforcement of the constraint. The proper value of the penalty parameter is determined prior to the main calculation via a trial and error process and it is problem dependent.
Minimizing the functional in (15) with respect to the nodal unknown vector leads to the following system of equation: where and is unknown matrix that contains displacements and stresses of all nodes. is the right hand side vector and the stiffness matrix in (16) is square matrix where is the number of unknowns per each node and is the moving least squares (MLS) shape functions. The proposed MDLSM method has 2.5 times unknowns compared to the irreducible DLSM method with only displacements as unknowns. This should drastically reduce the computational efficiency by (where depends on the type of linear solver used) times. But, since matrix is symmetric and positive definite, therefore, the final system of equation can be solved directly via efficient solvers. The MDLSM formulation for the plane elasticity problem has the following advantages that increase its efficiency [33].(1)The order of shape function derivatives is reduced by one order, thus complex and costly second-order derivative calculations of the MLS shape function in the irreducible DLSM method are avoided.(2)The stresses are obtained directly, while calculation of stresses in the irreducible DLSM method requires some postprocessing.(3)Only a linear complete polynomial basis is needed to construct the MLS shape functions, while in the irreducible DLSM formulation, the second order polynomial basis is required to achieve second-order consistency due to the presence of second-order derivatives in the irreducible formulation. This leads to lower computational effort for mixed method in construction of MLS shape functions and its derivatives compared to irreducible one. Furthermore, this in turn leads to the smaller number of points to be included in the support domain for the construction of the MLS shape functions in the mixed method and hence, less computational effort.(4)Both the displacement and stress boundary conditions are of the Dirichlet type which requires the specification of a single penalty parameter if a penalty method is used to enforce them as used in this work. In the irreducible DLSM method both Dirichlet and Neumann type boundary conditions are required for the determination of two types of penalty coefficient.(5)The mixed formulation, when used with the standard weighted residual methods both mesh-based and meshless forms, requires the LBB condition because the resulting problem is a saddle point problem. The least squares method, however, is a minimization method and therefore is not subject to the LBB condition.
3. Error Estimator and Adaptive Refinement
3.1. Error Estimator
In numerical methods, a problem is solved by discretization of the problem domain into the subdomains, hence, the governing equations only apply into these subdomains, so numerical methods always come with discretization error. Discretization error is one of the most important challenges in the numerical methods. Discretization error is theoretically decreased by refining the discretization domain but perfunctory refinement imposes the heavy computational cost without supplying the expected accuracy. Adaptive refinement methods mean balances between refinement procedure and its computational cost. These methods only refine locally the regions of the domain which has higher error. Adaptive procedure has two main parts: error estimation and adaptive refinement. Any success adaptive refinement needs a reliable error estimation procedure. Real error distribution can not be practicably used because the exact solution is not available for any practical problems. Several methods are used for error estimation with different numerical methods and these methods are categorized into two classes, namely, the residual based method [34] and recovery based method [35]. In residual based method, the residuals of differential equation and its boundaries are used as a criterion of error. The gradient of the solutions is used in recovery based method as the error criterion.
In this paper, the relative least square residuals functional for each node is defined as follows: where is the least square residuals functional in (15) and is unknowns obtained from the main solution. It is noticed that most of the computations of the least square residuals functional can be obtained from the main solution of the MDLSM method.
3.2. Adaptive Refinement Procedures
3.2.1. Node Moving
Mesh movement strategy can be easily and efficiently used with meshless methods since no element distortion is associated with the method. It should be noted that the mesh movement technique can be used in conjunction with the MDLSM method to adaptively adjust nodal points to improve the quality of the solution obtained with a prespecified number of nodal points. Here, a nodal refinement procedure is used that is called node moving adaptive refinement approach. When a node refinement is required, springs of prescribed stiffness are placed between each pair of nodes belonging to the same subdomain and the nodes are then moved until the spring system is in equilibrium.
In the node moving procedure, first, all nodes are connected with springs in which the neighbor nodes are defined using Voronoi diagram [36] (see Figure 1). Voronoi diagram is defined as where is Euclidean distance between and . The above equation means that neighbor nodes to node are the closest nodes to the node rather than other nodes.
Spring forces are defined as where is stiffness of spring between , and , are coordinates of , in equilibrium, respectively. The free body diagram is shown in Figure 2.
Spring stiffness is defined as a function of errors between two points , as follows: where and are the values of the error estimators obtained from (18) at nodes and , respectively, and is distance between these two connected nodes. In matrix form we have where and are the components of the force exerted at node in and directions, respectively, and and are the coordinates of node in and directions, respectively. The spring systems work as a two-dimensional truss such that nodes and springs are points and elements of truss, respectively. We mention that the system of algebraic equations (22) can be assembled in its standard finite element concept to yield global force vector for whole system of springs as follows: where is the stiffness matrix of the system calculated by assembling the stiffness matrices of all the springs defined in the system and represents the vector of nodal forces. In the equilibrium condition, the vector of assembled spring forces should be equal to zero. This requirement leads to the following system of algebraic equation which should be solved for the unknown vector of nodal position , that is, solve It is obvious that the equation system defined in (24) is singular before any boundary conditions are considered. The boundary conditions used here for solving this system of equations are defined by the requirement that the boundary nodes should not be allowed to move perpendicularly to the boundaries. In other words, boundary nodes only can be displaced along the boundaries which they have been placed on. Mathematical representation of these boundary conditions can be defined as where and are the initial coordinates of boundary node ; and represent the displaced final position of node and is outward unit vector normal to the boundary at node . This condition guarantees that the nodes initially located at the intersection of two boundary lines must remain on its initial position. Boundary conditions for a simple net of springs are shown in Figure 3.
Upon solving the system of equations with appropriate boundary conditions, the refined position of nodes is obtained leading to substantial reduction of the local and global error of the numerical solution in the subsequent analysis. The efficiency and effectiveness of the proposed adaptive refinement technique is verified in the next section by its application to benchmark test examples in plane elasticity.
3.2.2. Node Enrichment
With the meshless methods, the enrichment strategy only requires that the locations of new nodes to be added are determined without requiring to define the connectivity of the resulting configuration. Different methods can be thought to be defining the location of the nodes to be added to the current nodal configuration. Here the new nodes are added in the neighborhood of existing nodal points defined by a Voronoi diagram.
Once the Voronoi cells are defined, the vertices of Voronoi cells corresponding to the nodes with higher error than the average error over the domain are considered as the new nodes to be added to the current nodal configuration. The method is schematically illustrated in Figure 4. The value of average error over the domain is obtained by The above node enrichment adaptive refinement procedure based on error estimation is schematically shown in Figure 5.
4. Numerical Experiments
In this section, we are solving the benchmark examples by using the proposed node moving and node enrichment refinement procedures for the MDLSM method and comparing the results with the refinement procedure in the irreducible DLSM method and the exact analytical solutions or finite element results with very fine mesh.
Example 1 (an infinite plate with a circular hole). In the first example, consider the case of an infinite plate with a circular hole subjected to a uniaxial traction at infinity, as shown in Figure 6. Due to symmetry, only the upper right square quadrant of the plate is modeled. The edge length of the square is , where is the radius of the circular hole. This example is chosen because the exact analytical solution is available from Timoshenko and Goodier [37]. The solutions for the displacements and the stresses under a unit uniaxial stress along the axis are given as follows:
which is the shear modulus and where is the Poissons ratio. In this example, the constant values are , , and and .
For the node moving procedure, initial and refined nodal configurations with 122, 229, and 305 nodes are shown in Figure 7. The nodal points are refined with respect to the proposed error estimator based on the least square residuals functional. In Figure 8(a) numerical results of normal stress on the left edge is compared with the exact analytical solution. Figures 8(b) and 8(c) compare the convergence curve and the convergence rate of the MDLSM method for the node moving adaptive refinement strategy and plotted using the error norms in Table 1.
(a) 122 initial nodes
(b) 122 refined nodes
(c) 229 initial nodes
(d) 229 refined nodes
(e) 305 initial nodes
(f) 305 refined nodes
(a)
(b)
(c)
Also, the node enrichment adaptive refinement strategy is applied to the MDLSM method for obtaining accurate results. As shown in Figure 9, first the problem domain is discretized by using 122 initial nodal distributions, and then the nodal points are enriched in three steps (182, 288, and 521 nodes) on the region determined by the error indicator. Figure 10(a) compares the numerical results of normal stress on the left edge with exact analytical solution. The error norms in Table 2 are used for plotting and comparison of the convergence curve and the convergence rate of the node enrichment procedure for the MDLSM method (see Figures 12(b) and 12(c)). The results clearly show that the node moving and multistage node enrichment adaptive refinement strategies in the MDLSM method are more efficient compared to the refinement procedures in the DLSM method. Figures 8(b) and 8(c), 10(b), and 10(c) indicate that by using the node moving and the node enrichment adaptive refinement procedures, the convergence rate of the MDLSM method is increased.
(a) Initial nodes = 122
(b) First step enrichment nodes = 182
(c) Second step enrichment nodes = 288
(d) Third step enrichment nodes = 521
(a)
(b)
(c)
Example 2 (a cylinder subjected to an internal pressure). As a second elastostatic benchmark example a cylinder subjected to an internal pressure is considered. Due to the symmetry, only a quarter of the cylinder is modeled; see Figure 11. The boundary conditions are illustrated in Figure 11. The exact analytical solution of this problem is where the constant values are , , , , and . As shown in Figure 12, four types of nodal distributions with 106, 214, 359, and 543 nodes are distributed to solve and refine the nodes in the problem domain.
(a) 106 initial nodes
(b) 106 refined nodes
(c) 214 initial nodes
(d) 214 refined nodes
(e) 359 initial nodes
(f) 359 refined nodes
(g) 543 initial nodes
(h) 543 refined nodes
Figure 13(a) compares the normal stress at for initial and refined nodal configurations with 543 nodes. It is clear that the result of refined nodal configuration is more similar to the exact analytical solution than initial nodal configuration. Table 3 compares the error norms of the node moving procedure for the initial and refined nodal configurations and Figures 13(b) and 13(c) compare the convergence curve and the convergence rate, respectively.
(a)
(b)
(c)
In Figure 14, the initial and refined nodal distributions for the node enrichment strategy are shown. Figure 15(a) compares the normal stress at for initial and last steps refined with 543 nodal distributions. In Table 4 the error norm of the node enrichment refinement procedure is shown. Figures 15(b) and 15(c) compare the convergence curve and the convergence rate of the node enrichment adaptive refinement strategy.
(a) Initial nodes = 106
(b) First step enrichment nodes = 164
(c) Second step enrichment nodes = 276
(d) Third step enrichment nodes = 491
(a)
(b)
(c)
Example 3 (a reservoir fully filled with water). In this example, consider that the wall of a reservoir fully filled with water is investigated. The geometry of the wall is irregular as given in Figure 16. The material properties of the wall are given as Young’s modulus and Poissons ratio . The bottom of the wall is fixed and the curved edge of the wall is subjected to a hydrostatic pressure . Since the analytical solution of this problem is not available, a very fine mesh (with 59,400 linear triangular elements) FEM solution will be considered as our reference solution.
As shown in Figure 17, four types of nodal distribution with 84, 138, 218, and 299 nodes are distributed to solve and refine the points in the problem domain. Figure 18(a) compares the displacement in -direction along the vertical edge for initial and refined configurations. Table 5 compares the error norms of the node moving procedures based on MDLSM method and Figures 18(a) and 18(b) compare the convergence curve and convergence rate, respectively. For node enrichment adaptive refinement procedure, the problem domain is discretized with initial 84 points and is refined in three steps near high gradient error norm and is solved with 120, 181, and 284 nodes, respectively (see Figure 19). Figure 20(a) compares the displacement in -direction along the vertical edge for initial and refined configurations. In Table 6 the error norms of the node enrichment refinement procedure is used to plot the convergence curve and convergence rate in Figures 20(b) and 20(c).
(a) 84 initail nodes
(b) 84 refined nodes
(c) 138 initial nodes
(d) 138 refined nodes
(e) 218 initial nodes
(f) 218 refined nodes
(g) 299 initial nodes
(h) 299 refined nodes
(a)
(b)
(c)
(a) Initial nodes = 84
(b) First step enrichment nodes = 120
(c) Second step enrichment nodes = 181
(d) Third step enrichment nodes = 284
(a)
(b)
(c)
5. Conclusion
A mixed discrete least squares meshless method was extended for node moving and node enrichment adaptive refinements for efficient analysis of the planar elasticity problem. For the refinement procedures an error estimator based on least square residuals functional was formulated and used. Voronoi diagram was extended in the refinement procedures to find the neighbor nodes (node moving) and the position of the new nodes (node enrichment). For the moving node procedure, spring analogy was used to construct a system for computing the new place of each node after the refinement procedure. The efficiency and effectiveness of the proposed node moving and node enrichment adaptive refinement techniques in the MDLSM method by their application to the benchmark examples in the elasticity problems were verified. Results show that the proposed refinement methods are accurate and straightforward.
Conflict of Interests
The authors declare that there is no conflict of interest regarding the publication of this paper.