A Quasi-Conforming Point Interpolation Method (QC-PIM) for Elasticity Problems

It is well known that a high-order point interpolation method (PIM) based on the stan- dard Galerkin formations is not conforming, and thus the solution may not always be convergent. This paper proposes a new interesting technique called quasi-conforming point interpolation method (QC-PIM) for solving elasticity problems, by devising a novel scheme that smears the discontinuity. In the QC-PIM, the problem domain is ﬁrst discretized by a set of background cells (typically triangles that can be automatically generated), and the average displacements on the interfaces of the two neighboring cells are assumed to be equal. We prove that when the size of background cells approaches to zero, all the additional potential energy coming from the discontinuous displacement ﬁeld becomes zero, which ensures the pass of the standard patch test and hence the convergence. Numerical experiments verify that QC-PIM can produce the convergent solutions with higher accuracy and convergent rate that is in between fully conforming linear and quadratic models.


Introduction
The development of meshfree methods has achieved remarkable progress in recent years. Methods and techniques developed so far include the smooth particle hydrodynamic (SPH) method, the element-free Galerkin (EFG) method, the meshless On the other hand, Xu and Liu et al. proposed a novel conforming point interpolation method (CPIM) [Xu et al. (2010)] for generating higher order PIM shape functions without using the strain smoothing techniques. A solution with very high accuracy and convergent rate has been obtained using the CPIM. Recently, Xu proposed the ES-CPIM [Xu et al. (2012)] and NS-CPIM [Xu et al. (2011)] using the techniques combined by CPIM and smoothing technique. Using these schemes, we got the numerical solutions of some real problems with high accuracy and convergent rate.
In this paper, we extend the ideas of CPIM, and propose a novel quasiconforming point interpolation method (QC-PIM) without using the generalized strain smoothing technique. We develop a technique for shape function construction and guarantee the conformability when the size of back ground approaches to zero.
The paper is outlined as follows. Section 2 briefs the linear elasticity, and Sec. 3 gives a briefing on PIM and T-Scheme. The ideas of the QC-PIM are presented in Sec. 4. In Sec. 5, numerical examples and further investigations on QC-PIM are presented. Conclusions are drawn in Sec. 6.

Basic Equations for Planar Elastic Problems [Zienkiewicz and Taylor (2000)]
Let Ω be a planar region occupied by the elastic body, and Γ = ∂Ω is its boundary. There are three groups of state variables: the stress σ = (σ 11 , σ 22 , σ 12 ) T , the strain ε = (ε 11 , ε 22 , γ 12 ) T and the displacement u = (u 1 , u 2 ) T . Then σ, ε and u satisfy the following equations: where b = (b 1 , b 2 ) T is the vector of body forces, L d is a matrix of differential operators defined as For the homogeneous and isotropic materials, the elasticity constant matrix D is: where E is Young's modulus and v is Possion's ratio. Suppose boundary Γ = Γ u + Γ t and Γ u ∩ Γ t = Φ. A displacement boundary condition u = u 0 is given on Γ u , and a natural boundary conditions are L T n σ = t on Γ t , where L n is the matrix of unit outward normal which can be expressed as Eliminating σ and ε from (1) to (3) yields a system equation of displacement u: where operator A = L T d (DL d ).

Briefing on the PIM and T-scheme
In the PIMs, the problem domain is first discretized by a set of background cells as shown in Fig. 1. Since triangular background cells can always be generated efficiently and automatically, PIMs usually use triangular background cells for node selection and numerical integration. This is called T-Schemes [Liu and Zhang (2008a)]. In the T-schemes, a home cell refers to the cell which holds the point of interest.
x Interested point Support points for T3-scheme Support points for T6-scheme  Neighboring cells refer to the cells which share one edge with a cell. T3-scheme is used to produce linear PIM shape function by using the linear basis function, and T6-scheme is used to produce quadratic PIM shape function. In the T3-scheme (see Fig. 1), the point of interest point (x) is located in a cell; only the three vertexes of the home cell are used for shape function construction. In the T6scheme, the six nodes are selected to interpolate an interested point located in a cell: three vertexes of the home cells and the other three vertices of the three neighboring cells.
The details for T6-schemes are as follows [Liu (2002)]. Consider a displacement component u(x) for the solid mechanics problems. It can be approximated as the following formulation where p i (x) is the monomial basis function of x = (x, y) T , a i is the corresponding coefficient, and n is the number of nodes in the local support domain. The complete polynomial basis of order 2 can be written as In the original quadric PIM, the coefficients in Eq. (9) can be determined by enforcing u(x) to displacements at three vertexes and the other three vertices of the three neighboring cells: for i = 1, 2, . . . , 6. Following the standard FEM process, we can obtain solution for a as: where and Note that the linear independence of three vertexes and the other three vertices of the three neighboring cells. Therefore, matrix P n must be nonsingular and Eq. (12) has a unique solution. It is clear from the above discussion that when the high-order interpolation is used for shape functions reconstruction, the displacement along interface edge is not continuous in between two adjacent cells. Hence, the original high-order PIM is not conforming. If only the linear basis function p T (x) = (1, x, y) and T3-scheme are used for shape function construction, then it is called linear PIM. The linear PIM is conforming.

Quasi-conforming PIM
In the QC-PIM, according to the Hellinger-Reissner's two-field variational principle [Liu and Zhang (2008b)], we construct the following energy functional for the discretized system equations: where ε = L d U is the assumed strain vector, stresses σ are dependent on the strain ε through the stress-strain relationship σ = Dε, and t is the traction on the boundary that is dependent on the stress σ in the form of L T n σ = t. Using Green's divergence theorem, the second term in the right-hand of Eq. (15) can be further expressed as Substituting Eq. (16) into Eq. (15) gives After the problem domain is discretized by a set of background triangular cells Ω i (i = 1, 2, . . . , n), Eq. (17) can be written as the discretized form: where t i is the tractions on the cells boundary

between a cell and his neighbor cells. Denotes
as the sum of additional potential energy of all the cells boundary, which is produced by the discontinuous of displacement field. For a conforming methods, displacement field satisfy u i =ȗ i on all the cells boundary ∂Ω i . Then for any traction t i we have H P = 0. Note that displacement field in the higher order PIMs is not continuous. However, the original PIM ignores this additional potential energy H P for the sake of simplicity, which may result in incompatibility.
In order to solve the incompatibility, this paper proposes a quasi-conforming PIM (QC-PIM). QC-PIM relaxes the conditions of displacement continuous between the adjacent background cells, and just satisfies some generalized conforming conditions H P = 0 when the size of cells approaches to zero. Here, for the quadric PIM, we give a simple scheme to satisfy the above generalized conforming conditions.
In a similar way to Sec. 2, assume that a component of displacement field u(x) can be approximated as Eq. (9), and the complete quadric polynomial basis in Eq. (10). There are six coefficients in Eq. (9) can be determined. We first enforce u(x) to displacements at three vertexes in the home cell of x: We now introduce the other three conditions for determining coefficients a whereȗ is the same component of the displacement obtained for a neighboring cell, b 1 b 2 b 3 are illustrated in Fig. 1. The displacement field constructed using (20) and (21) is not continuous. However, when the size of cells approached to zero, the traction t i on the boundary of a cell trends to a constant, and then Eq. (19) can be written as In deriving Eq. (22), we have used the conditions (21). It is clear from [Zienkiewicz and Taylor (2000)] the system can pass the standard patch test when the addition potential energy H P = 0. However, the quartic displacements in (20) and (21) satisfy generalized conforming conditions when the size of cells approached to zero. This is the difference between the proposed method and fully conforming methods, and it is the reason that the proposed method is called the QC-PIM.

Numerical Examples
In this section, a number of numerical examples will be examined using the QC-PIM with piecewise quadratic displacement field. To investigate quantitatively the numerical results, the displacement error and energy norms are defined as follows, where the superscript "ref " denotes the reference or analytical solution, "num" denotes a numerical solution obtained using a numerical method.

Standard path test
This section investigates if the QC-PIM can pass the standard path test. A rectangular patch of 10 × 50 is considered, and the displacements are prescribed on all outside boundaries by the following linear function The patch is represented using nodes shown in Fig. 2 and displacement errors obtained using QC-PIM are illustrated in Table 1. It is clear that when the size S of the background cells approaches to zero, the displacement error is also approach to zero. This example verifies numerically that the QC-PIM can pass the standard path test when the size of background cells approach to zero.

Cantilever 2D solid
A 2D cantilever solid with length L = 50 m and height D = 10 m is now studied. The solid is subjected to a parabolic traction at the right end as shown in Fig. 3. The cantilever beam is studied as a plane stress problem with E = 3.0 × 10 7 Pa, P = −1000 N and v = 0.3. Analytical solutions of this problem can be found in [Timoshenko and Goodier (1970)].   Fig. 3. A 2D cantilever solid subjected to a parabolic traction on the right edge.
Using Eqs. (23) and (24), errors in displacement and energy norms are computed and plotted in Figs. 4 and 5. When QC-PIM is used, the convergence rate of 2.86 in displacement norm is found, this is much higher than the original quadratic PIM. The convergence rate in energy norm is found as high as 1.52 for QC-PIM as shown in Fig. 5, which is much higher than that of original quadratic PIM. These numerical results indicate that QC-PIM is of convergence with high accuracy and high convergence rates.

Infinite 2D solid with a circular hole
An infinite 2D solid with a central circular hole (a = 1 m) and subjected to a unidirectional tensile (T x = 10 N/m) is studied. Owing to its two-fold symmetry, one quarter is modeled with b = 5 m (as shown in Fig. 6). Symmetry conditions are imposed on the left and bottom edges of the quarter model, and the inner boundary of the hole is traction free. For this benchmark problem, the analytical solution can be found in [Timoshenko and Goodier (1970)]. Using Eqs. (23) and (24), errors in displacement and energy norms are computed and plotted against the average nodal spacing (h) as shown in Figs. 7 and 8. When QC-PIM is used, the convergence rate of 2.75 in displacement norm is found, this is much higher than the original quadratic PIM. The convergence rate in energy norm is found as high as 1.41 for QC-PIM as shown in Fig. 8, which is much higher than that of original quadratic PIM. These numerical results indicate that QC-PIM is of convergence with high accuracy and high convergence rates.

Convergence rates
By the theory of standard FEM, the convergence order of a computational scheme is determined by the order of trial functions [Zienkiewicz and Taylor (2000)].
Quadratic QC-PIM satisfies generalized conforming conditions and corresponding interpolation functions are the piecewise quadratic polynomials which indicates that the convergence rate is higher than that in linear FEM. It should be noted that QC-PIM has not introduced any additional nodes compared to the linear FEM; and the number of nodes in quadratic QC-PIM is obviously less than those in quadratic FEM when the same background cells are used. Considering the irregular nodal distribution in real computation, the rough estimation of the convergence rate in quadratic QC-PIM is in between linear FEM and quadratic FEM: namely, convergence rate in displacement norm is 2 < r d < 4 and 1 < r e < 2 in energy norm. The convergence of CPIM has been verified by the numerical examples in this paper.

Band width of stiffness matrix
The band width of stiffness matrix is very important for the computational efficiency [Xu et al. (2010)]. Here, we give an estimation of band widths of stiffness matrix for different methods. Consider a simple problem domain and corresponding triangular background cells as shown in Fig. 9. The influence domain of a node includes usually some closely related nodes around this interested node. These nodes in influence domain are of very importance to band width of stiffness matrix. In Fig. 9, • denotes an interested node i; • and is the related nodes of node i for linear FEM; •, and are related nodes for quadric QC-PIM, ES-PIM (T3) and the original quadric PIM; •, , and are related nodes for quadratic CPIM. There are more nodes for the other high-order S-PIM including quadric NS-PIM, ES-PIM, ES-CPIM, NS-CPIM. Figure 9 illustrates the influence domains of a node for the different schemes. It is important that the band width of stiffness matrix in quadratic QC-PIM is exactly the same with that in original quadric PIM, and far less that in CPIM. In addition, considering the extra computing cost in constructing shape functions, QC-PIM will take a little more computation time compared to original quadric PIM, and is very similar to that in ES-PIM (T3) with the same nodes.
Note that quadratic QC-PIM assumes a piecewise quadratic polynomials displacement field. So, under the same computational accuracy, computational cost of QC-PIM is far less than that in CPIM; under the same computational cost, computational accuracy is better than that in ES-PIM (T3). Therefore, the performance of QC-PIM is better than that of ES-PIM (T3) and quadric CPIM.

Conclusion
This paper proposes a new technique to produce a quasi-conforming point interpolation method (QC-PIM) for solving elasticity problems. In the QC-PIM, the problem domain is first discretized by a set of background cells, and the average displacements on the boundary of a cell are assumed to be equal to that on its neighboring cells. We prove theoretically that when the size of background cells approach to zero, the additional potential energy coming from the boundary force of all the background cells are zero which ensures the pass of patch test, and thus the convergence of QC-PIM. Numerical experiments verify the QC-PIM can pass the standard patch test, and produce the solutions with higher accuracy and high convergent rate.