Comparison of meshfree displacement and stress error recovery of finite element solutions using moving least squares interpolation

Abstract

The moving least square (MLS) interpolation based the recovery procedures have been successfully applied to recover the finite element solution errors in the analysis of elastic plates and pipes problems and can be advantageously applied for large deformation and fracture problems. The study presents the displacement and stress error recovery characteristics in the error estimation analysis employing Moving Least Squares interpolation approach. The study considers quartic spline, cubic spline, and exponential weight function with three different order of basis function in Moving Least Squares interpolation based error recovery analysis. The displacement/stress errors in finite element solution are quantified in energy norm. The cylinder and plate benchmark examples using triangular and quadrilateral elements are analyzed to compare the convergence, effectivity and adaptively improved meshes obtained using the various displacement/stress recovery procedures. The study shows that cubic spline weight function and quadratic basis function found to perform better in MLS based meshfree recovery technique for stress as well as displacement errors recovery of finite element solution. It is observed from the study that increasing the order of basis function will enhance the error estimation quality that is, rate of convergence become faster with improved effectivity of the results. The increase in convergence rate with the increase of the order of basis function is more in displacement recovery technique as compared to stress recovery technique. It is observed in the analysis of benchmark example with linear triangular meshing that the error reduction using meshfree MLS interpolation based displacement and stress recovery is about 10% and 150% respectively for the displacement and stress recovery over the mesh dependent least square based displacement and stress recovery. The study concludes that the effectiveness and efficiency of meshfree displacement/stress error recovery technique strongly depends on the weight and basis functions of MLS method to recover the errors.

Introduction

During recent years, the finite element technique has been used successfully in modeling complex static and dynamic problems in various areas of engineering and procedures are developed for error estimation in finite element solutions. Error estimation based adaptivity techniques for finite element method has significantly increased the reliability, safety, and cost effectiveness of finite element codes for solving engineering problems. A literature review on wide applicability of finite element methods is conducted by Cen et al. ¹ Chen et al. ² presents the review of last two decades literature in the mesh less method. The recovery type error estimation, relying upon recovery of higher order finite element solution using the mesh dependent node patches, is proposed by Zienkiewicz-Zhu. ³ In classical recovery methods such as Zienkiewicz-Zhu recovery methods, ⁴ interpolants of unknown parameters are related to the element’s geometry. However, in meshless recovery methods, such as meshfree Galerkin recovery technique based on Moving Least Squares (MLS) method, approximations of unknown parameters require distributed nodes in support domain without defining the element mesh and connectivity. In the MLS approximation, continuous field variables or their derivatives approximation are derived from a set of point values using weighted least squares approximation. As there is no dependency of error recovery on the domain mesh, the meshfree MLS recovery techniques is most suitable approach to situations in which the distortion of elements occur, such as for large domain changes and domain discontinuities such as for plasticity and fracture problems. The weight function in weighted least squares method for a specific given node is positive only on the patch of nodes close to a given node up to a determined distance, that is, its support comprises a patch of neighboring nodes around a vertex node. Oñate et al. ⁵ has discussed the completeness and continuity of MLS interpolation. They proved that MLS procedure interpolate random data with improved accuracy. The improvement in the moving least squares method is proposed by Most and Bucher ⁶ developing an interpolating non-singular weight function which lead to moving least squares shape function fulfilling the interpolation condition exactly. Ubertini ⁷ has proposed the recovery by patch equilibrium and recovery by patch compatibility to get improved stresses to the explicitly measured stresses. Ródenas et al. ⁸ have put forwarded the improvement of the super-convergent patch recovery (SPR) technique, called SPR-C technique (Constrained SPR) and found that the approach recovered the stress field more accurately and improve the elemental effectivity of the Zienkiewicz-Zhu error estimation.

The shape functions in meshfree Galerkin formulation is derived by Cao et al. ⁹ using MLS technique which avoids the moment matrix inversion. A hybrid finite element method (FEM) and the element-free Galerkin method (EFGM) based on local maximum entropy shape functions is proposed by Ullah et al. ¹⁰ for linear and nonlinear problems having material as well as geometrical nonlinearities. The Zienkiewicz and Zhu (ZZ) super-convergent patch recovery for strains and stresses is used in the FE region of the problem domain, while the Chung and Belytschko ¹¹ error estimates is used in the EFGM region. Mirzaei ¹² has conducted a comprehensive analysis of error estimation considering MLS approximation. The stability of the moving least squares (MLS) approximation is analyzed and discussed by Li and Li. ¹³ They found that the stability of the MLS approximation deteriorates severely with the decrease of nodal spacing. Parret-Fréaud et al. ¹⁴ presented an MLS recovery based procedure to obtain continuous stresses field in which the continuity of the smoothed field is provided by the shape functions of the underlying mesh. Liu et al. ¹⁵ developed a hybrid technique using wavelet-Galerkin method and finite element method and established its accuracy and stability for 2D to 3D elasticity problems. For 3D-finite elements analysis, more accurate post-processed stress field can be found by projecting the directly computed stress field onto the conveniently chosen space. ¹⁶ Ahmed et al. ¹⁷ have developed MLS based Meshfree recovery technique for elastic problems with super-convergent properties. The effect of recovery parameters such as support domain shape and size, polynomial expansion order in basis function, on solution error recovery using MLS based recovery approach in meshfree environment is investigated by Ahmed et al. ¹⁸ They have suggested optimal values of polynomial expansion order in basis function and dilation parameter for improving the performance of error estimation.

The development of meshfree and mesh dependent recovery methods for more accurate displacement and stress, motivate the author to compare the displacement/stress recovery characteristics obtained by different recovery techniques. The most of study in the mesh free recovery techniques is concerned with recovery of displacement errors of finite element solution. The research studied are also not devoted to compare the mesh free recovery techniques to popular and conventional mesh dependent recovery techniques. The present study aimed to use meshfree and mesh dependent recovery techniques in adaptive finite element analysis. The study investigates the impact of choice of element for discretization, weight functions, order of polynomial basis functions and dependency of node with mesh connectivity on finite element displacement and stress recovery. The techniques, for recovery of displacement and stress errors, employed in the study are recovery of displacement and stress field utilizing MLS procedure over meshfree node patch, recovery of displacement field employing least square approach over a mesh dependent element neighbor patches and recovery of stress field (ZZ recovery) employing least square approach over a mesh dependent node neighbor patch. The MLS interpolation procedure considers three weight functions, namely cubic spline, quartic spline, and exponential, and three set of polynomial order of basis functions. The errors in finite element solution are computed directly from the post-processed displacement/stress values. Adaptive finite element numerical experiment is carried out on cylinder and plate examples for error and convergence characteristics of the recovered solution at decreasing order of mesh density. The triangular/quadrilateral elements are used for domain discretization. The effectiveness of meshfree MLS interpolation recovery techniques is compared with effectiveness of mesh dependent displacement/stress recovery techniques employing different elements, weight function, and basis function. The efficiency of various recovery procedures is compared through error estimator guided adaptive improvement of the domain discretization.

Moving least squares interpolation method

The MLS method is presented by Lancaster and Salkauskas. ¹⁹ The MLS technique provide the field variables or their derivatives approximation at a point using the interpolation of the field variables or their derivatives at the nodes which lie in the support domain of the point in a weighted least square sense. The support domain or patch of arbitrarily distributed nodes is a zone up to a determined distance with influenced nodes (Figure 1). The procedure to get continuous approximation values at points may be used as a Meshfree error recovery that does not require element connectivity but only needs definition of nodes and description of boundaries. The MLS technique uses three components to express the function u(x) with the approximation u ^h (x) namely, a basis functions p(x), a set of coefficients a(x), and a weight function w(x) associated to each node. The many basis functions can be selected within a weighted least squares framework and the coefficients are functions of the coordinates of the nodes distributed within support domain. MLS approximation is desired should include enough supporting nodes in whose domain of influence the desired location resides. This number of nodes should be at least equal to the number of terms in the basis but should be significantly larger. It is also beneficial to construct the domain of influence to be smaller so that the local character of the approximation is maintained. Suppose the nodes in support domain are defined by x_i…x_n where x_i = (x_i, y_i) in two dimensions. In MLS, the approximated value of the field function or their derivatives (the contribution of point xI in the support domain of any point x to the field variable (or derivatives) at the point x) can be computed as follows.

$${\mathrm{u}}^h\left(x\right)=p^T\left(x\right)a\left(x\right)=\sum _{j=1}^m{p}_j\left(x\right)a_j\left(x\right)$$

(1)

Figure 1.

Meshfree node patches for MLS technique (displacement/stress recovery).

where m is the number of polynomial basis and a(x) is the vector of coefficient given by

$$a\left(x\right)={\{a_0\left(x\right)a_1\left(x\right)\dots \dots \dots a_{m-1}\left(x\right)a_m\left(x\right)\}}^T$$

(2)

where a₀(x) are the function of coordinates and p(x) is the vector of basis function that generally consist of monomial of the lowest order to ensure minimum completeness. Basis functions p(x) matrix and linear, quadratic, cubic orders of basis function are given below,

$$p=\left[\begin{array}{c}\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{1}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{1}\end{array}\begin{array}{c}\begin{array}{c}{x}_1\hfill \\ \multicolumn{1}{c}{x_2}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{x_n}\end{array}\begin{array}{c}\begin{array}{c}{y}_1\hfill \\ \multicolumn{1}{c}{y_2}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{y_n}\end{array}\begin{array}{c}\begin{array}{c}{x}_1{y}_1\hfill \\ \multicolumn{1}{c}{x_2{y}_2}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{x_n{y}_n}\end{array}\right]$$

(3)

$$p^T(x,y)=[1,x,y,]\mathrm{for}m=3$$

(4)

$$p^T(x,y)=[1,x,y,x^2,\mathrm{xy},y^2]\mathrm{for}m=6$$

(5)

$$p^T(x,y)=[1,x,y,x^2,x^3,x^{2}y,x{y}^2,y^2,y^3]\mathrm{for}m=9$$

(6)

The vector of coefficients a(x) can be obtained by minimizing a weighted residual as follows:

$$J=\sum _{j=1}^{n}w(x-x_I)[P^T\left(x_I\right)a\left(x\right)-u_I{]}^2$$

(7)

$$\frac{\partial J}{\partial a}=A\left(x\right)a\left(x\right)-B\left(x\right)u_s=0$$

(8)

Minimization of weighted residual leads to the following expression of the coefficient vector.

$$a\left(x\right)={\mathrm{A}}^{-1}\left(x\right)B(x)u_s$$

(9)

where A is called the MLS moment matrix given by

$$A\left(x\right)=\sum _{i=i}^n{w}_i(x-x_i)p^T\left(x_I\right)p\left(x_I\right)$$

(10)

where B(x) has the form

$$B\left(x\right)=[w_I(x-x_I)p\left(x_I\right),\dots \dots ,\dots .w_n(x-x_n)p\left(x_n\right)]$$

(11)

where u_s is the vector of nodal parameters of the field variables (or derivatives) for all nodes of the support domain and w(x) is the weight functions given as,

$$w(x)=\left[\begin{array}{c}\begin{array}{c}w(x-x_1)\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{w(x-x_2)}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{0}\end{array}\begin{array}{c}\begin{array}{c}\dots \hfill \\ \multicolumn{1}{c}{\dots }\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\dots }\end{array}\begin{array}{c}\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{w(x-x_n)}\end{array}\right]$$

(12)

Therefore, the approximated field variables can be computed as follows.

$$u^h\left(x\right)=\sum _{I=1}^n\sum _{j=1}^m{p}_j\left(\mathrm{x}\right){\mathrm{A}}^{-1}\left(x\right)B(x{)}_{\mathrm{jI}}{u}_s$$

(13)

MLS interpolation weight functions

The weight function plays an important role in performance of MLS interpolation method. The number of weight functions can be selected for MLS interpolation technique. However, the nodal weight functions should satisfy conditions that the weight functions should be positive only on the domain of influence that is, within the support domain. The weight functions should be equal to zero outside the support domain or in other words, the weight function vanishes when nodes do not lie in the support of point. The value of the weight functions increases as the distance between xI and x decreases, and vice versa. The domains of influence in the circular form that is, the circular support boundary is constructed using the radial distance $\overline{\mathit{d}}=\Vert \mathit{x}-{\mathit{x}}_{\mathit{i}}\Vert /{\mathit{d}}_{\mathit{m}}$ . The (x−x_i) is the distance from node x to point x_i and d_m is the size of influence domain of the point x_i, the support size of the Ith node, d_mI, is computed by d_mI = d_max c_I, in which d_max is a scaling factor called as dilation parameter, the distance c_I is determined by searching for sufficient neighbor nodes distances. The dilation parameter (d_max) is taken as 3.0 in the present study as an optimized value suggested in Ahmed et al. ¹⁸ The following weight functions are used in the present study.

Cubic Spline Weight Function:

$$w\left(\overline{d}\right)=\left\{\begin{array}{cc}\frac{2}{3}-4{\overline{d}}^2+4{\overline{d}}^3\hfill & \mathrm{for}\overline{d}\le \frac{1}{2}\hfill \\ \multicolumn{1}{c}{\frac{4}{3}-4\overline{d}+4{\overline{d}}^2-\frac{4}{3}{\overline{d}}^3}& \mathrm{for}1\le \overline{d}\le \frac{1}{2}\hfill \\ \multicolumn{1}{c}{0}& \mathrm{for}\overline{d}>1\hfill \end{array}\right\}$$

(14)

Quartic Spline Weight Function:

$$w\left(\overline{d}\right)=\left\{\begin{array}{cc}1-6{\overline{d}}^2+8{\overline{d}}^3-3{\overline{d}}^4\hfill & \mathrm{for}\overline{d}\le 1\hfill \\ \multicolumn{1}{c}{0}& \mathrm{for}\overline{d}>1\hfill \end{array}\right\}$$

(15)

where w( $\overline{d}$ ) is a weight function in 2-D associated to each node. The weights take a unit value in the vicinity of the point where the function and its derivatives are to be computed and vanishes outside a region Ωi covering the point x_i, and $\overline{\mathit{d}}=\Vert \mathit{x}-{\mathit{x}}_{\mathit{i}}\Vert /{\mathit{d}}_{\mathit{m}}$ .

Exponential Weight Function:

$$w\left(\overline{d}\right)=\left\{\begin{array}{cc}{e}^{-{(\frac{\overline{d}}{\alpha })}^2}\hfill & \mathrm{for}\overline{d}\le 1\hfill \\ \multicolumn{1}{c}{0}& \mathrm{for}\overline{d}>1\hfill \end{array}\right\}$$

(16)

where α is a shape parameter constant. The shape parameter constant is taken as 0.4 in the present study. The constant (α) controls the exponential weight function. The higher weights are assigned on points x_i close to x while lower weights on points far from x.

The smoothness level of cubic spline (W1), quartic spline (W2), and exponential (W3) weight functions are depicted in Figure 2.

Figure 2.

Smoothness level of cubic spline (W1), quartic spline (W2), and exponential (W3) weight functions. ²⁰

Stress error recovery technique (ZZ recovery)

The Zienkiewicz-Zhu super-convergent patch recovery technique ⁴ for finite element solution errors recovery assumes that the node values of field variable derivatives belong to a polynomial expansion of the same complete order as of the basis function and is valid over a patch of neighboring nodes of the vertex (Figure 3). Following polynomial expansion may be used for each component of stress.

$$\mathbf{\sigma }\mathbf{*}(\mathit{x})=\mathit{P}(\mathit{x}).\mathit{a}$$

(17)

Figure 3.

Mesh dependent node patches for least square (ZZ) technique (stress recovery): (a) triangular elements and (b) quadrilateral elements.

where P(x) is the basis function of the assumed polynomial, x = (x_i, y_i), are the coordinates of the sampling points and a is the unknown parameters vector. A least square fit of σ ^h values over the nodes patch, may be made by minimizing the following functional.

$${\pi }_f(a)=\frac{1}{2}\sum _{i=1}^{\mathrm{np}}[{\sigma }^h(x_i,y_i)-P(x_i,y_i).\mathbf{a}{]}^2$$

(18)

On simplification, it results into the following equation.

$$\mathit{A}a=\mathit{b}$$

(19)

The matrix A and b are given as.

$$\mathit{A}=\sum _{i=1}^{\mathrm{np}}[P^T(x_i,y_i)P(x_i,y_i)]$$

(20)

$$\mathit{b}=\sum _{i=1}^{\mathrm{np}}[P^T(x_i,y_i){\sigma }^h(x_i,y_i)]$$

(21)

where np is the number of nodes in a patch.

Displacement error recovery technique

The recovery of field variable (displacement) is obtained by least squares fit of the computed nodal field variable using a higher order polynomial over an element neighborhood patch consists of union of the elements surrounding an element (Figure 4). To perform least square fitting, the following functional is minimized.

$${\pi }_f(a)=\frac{1}{2}\sum _{i=1}^{\mathrm{np}}[u_i^h(x_i,y_i)-u(x_i,y_i){]}^2$$

(22)

Figure 4.

Mesh dependent node patches for least square technique (displacement recovery): (a) triangular elements and (b) quadrilateral elements.

where u_i (x_i, y_i) = P_i (x_i, y_i) ξ a

$${\mathit{u}}_{\mathit{i}}={\left[u_i{v}_i\right]}^{\mathrm{T}},\mathit{a}={\left[a_u{a}_v\right]}^{\mathrm{T}}$$

(23)

where u_i and v_i are the nodal parameters of field variables in x and y direction and a is the vector of unknown parameters a_u and a_v.

$${\mathit{P}}_i=\left[\begin{array}{cc}{\mathit{p}}_{\mathit{i}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\mathit{p}}_{\mathit{i}}\hfill \end{array}\right]$$

(24)

$${\mathit{P}}_{\mathit{i}}=[1,x_i,y_i,x_i^2,x_i{y}_i,y_i^2,.......]$$

(25)

where (x_i, y_i ) are the sampling points coordinates.

Minimization condition of π_f (a) implies that a satisfies the following relation.

$$\sum _{i=1}^{\mathrm{np}}{P}_i^T(x_i,y_i).P_i(x_i,y_i).a=\sum _{i=1}^{\mathrm{np}}{P}_i^T(x_i,y_i)\mathrm{x}{u}_i^h(x_i,y_i)$$

(26)

Solving for “a,” the following relation is obtained.

$$\mathit{a}={\mathit{A}}^{-1}\mathit{b}$$

(27)

Where

$$\mathit{A}=\sum _{i=1}^{\mathrm{np}}{P}_i^T(x_i,y_i)P_i(x_i,y_i)$$

(28)

$$\mathit{b}=\sum _{i=1}^{\mathrm{np}}{P}_i^T(x_i,y_i)u_i^h(x_i,y_i)$$

(29)

Quantification of error estimates

The error in computed state variable or state variable derivative that is, displacement (u) or stress (σ), e_u^* (or e_σ^*) is defined as the difference between the exact (or numerical) values of u (or σ) and respective computed values, u ^h (or σ^h) that is,

$${e_u}^{*}=u-u^h$$

(30)

$${e_{\sigma }}^{*}=\sigma -{\sigma }^h$$

(31)

The error is generally measured in “norm” representing the integral scalar quantities. For example, the energy norm is the strain energy contained in the difference between the discontinuous value σ^h and the recovered values σ*. The “energy norm” can be written as follows.

$$\Vert e{\Vert }_E={\left[\underset{\mathrm{\Omega }}{\int }{e}_{\sigma }^{*T}{D}^{-1}{e}_{\sigma }^{*}d\mathrm{\Omega }\right]}^{\frac{1}{2}}$$

(32)

where D is the elasticity matrix.

An estimator is asymptotically exact for a problem if the problem global (representing overall of all elements) and local (element) effectivity index (θ) that is, ratio of estimated error and actual error, converges to one when the mesh size approaches to zero.

$$\theta =\frac{\Vert e\Vert }{\Vert e_{\mathrm{ex}}\Vert }$$

(33)

where $\Vert e_{\mathrm{ex}}\Vert$ and $\Vert e\Vert$ denote the actual error and the estimated error estimate in energy norm.

The accuracy (η) of a finite element solution may be defined as follows.

$$\eta =\frac{\Vert e_{\sigma }\Vert }{\Vert {\sigma }^{*}\Vert }$$

(34)

$$\Vert {\sigma }^{*}{\Vert }^2=\Vert {\sigma }^h{\Vert }^2+{\Vert e{\Vert }_E}^2$$

(35)

$$\Vert {\sigma }^h{\Vert }^2=\sum _{i=1}^n\Vert {\sigma }^h{\Vert }_i^2$$

(36)

The solution is acceptable if η ≤ η _allow where η _allow is the allowable accuracy. If η > η _allow , modification of element size is needed.

Benchmark example 1: Cylinder subjected to pressure

The effectiveness and efficiency of the stress and displacement recovery in meshfree MLS and mesh dependent least square (LS) based error estimation is assessed through the finite element analysis of the thick cylinder under pressure example “P.” The analytical solution for the cylinder example is given by the following relations. ²¹

For a point (x, y), $\mathrm{c}=\mathrm{b}/\mathrm{a},\mathrm{r}=(\surd {\mathrm{x}}^2+{\mathrm{y}}^2)$ with “a” and “b” as the internal and external radius of cylinder, the displacement in radial direction is given by,

$${\mathit{u}}_{\mathit{r}}=\left[P(1+\nu )\right]/[E(c^2-1)][r(1-2\nu )+\left(\frac{b^2}{r}\right)]$$

(37)

The stresses in cylindrical coordinate are given by,

$${\mathbf{\sigma }}_{\mathit{r}}=[P/(c^2-1)][(1-\left(\frac{b^2}{r}\right)]$$

(38)

$${\mathbf{\sigma }}_{\mathit{\theta }}=[P/(c^2-1)][(1+(\frac{b^2}{r})]$$

(39)

where P, E, and ν are pressure, Modulus of elasticity, and Poisson’s ratio respectively with a value of 1.0, 1000, 0.3.

The linear triangular and quadrilateral elements have been used for the discretization of the cylinder domain. The cylinder domain, quarter part due to symmetry, with linear triangular and quadrilateral element having uniform and non-uniform discretization in distribution is shown in Figure 5. The test examples consider cubic spline, quartic spline, and exponential weight functions, linear, quadratic, cubic order of polynomial basis function (m = 3, 6, 9) and circular support domains for meshfree node distribution and, element neighborhood patch/node neighborhood patch for mesh dependent node distribution.

Figure 5.

Meshing schemes for cylinder problem domain: (a) triangular element meshes (uniform discretization), (b) triangular element meshes (non-uniform discretization), and (c) quadrilateral element meshes (uniform discretization).

Displacement and stress recovery with element type and patch configuration

The characteristics of the displacement and stress error estimation are obtained in the finite element analysis of benchmark thick cylinder problem with different elements discretization types that is, linear triangular, and linear quadrilateral element, and patch configurations that is, meshfree and mesh dependent uniform/non-uniform patch, in term of the error convergence and effectivity considering cubic spline weight function and quadratic order of basis function (m = 6). Tables 1 to 3 present error in energy norm and global effectivity for linear triangular and quadrilateral elements mesh, and meshfree and mesh dependent patches employed for displacement and stress finite element solution error recovery types. Table 1 also shows the percentage reduction in error in MLS based displacement and stress error recovery over the LS based displacement and stress error recovery for linear triangular meshing.

Table 1.

Error (energy norm) and global effectivity (θ) for cylinder problem using displacement/stress recovery and patch configuration, and percentage reduction of error using MLS over LS technique (linear triangular element, cubic spline, uniform mesh).

Mesh		FEM error (×10−3)	MLS (meshfree patch), m = 6				LS (mesh dependent patch)				%age reduction of error using MLS over LS
			Displacement		Stress		Displacement		Stress (ZZ)		Displacement	Stress
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
125	162	4.05	1.63	0.936	1.93	0.990	1.82	0.953	3.91	1.171	19	102
272	324	2.77	0.85	0.941	0.97	0.964	0.97	0.946	2.54	1.171	14	162
550	622	2.00	0.67	0.930	0.66	0.937	0.70	0.930	1.67	1.114	5	153

Table 2.

Error (energy norm) and global effectivity (θ) for cylinder problem using displacement/stress recovery and patch configuration (linear triangular element, cubic spline, non-uniform mesh).

Mesh		FEM error (×10−3)	MLS (mesh free patch), m = 6				LS (mesh dependent patch)
			Displacement		Stress		Displacement		Stress (ZZ)
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
136	162	3.85	1.48	0.969	1.86	1.009	1.70	0.964	3.59	1.194
291	344	2.73	0.92	0.958	1.04	0.977	1.00	0.942	2.46	1.173
587	758	2.02	0.72	0.934	0.72	0.937	0.74	0.923	1.75	1.128

Table 3.

Error (energy norm) and global effectivity (θ) for cylinder problem using displacement/stress recovery and patch configuration (linear quadrilateral element, cubic spline).

Mesh		FEM error (×10−3)	MLS (mesh free patch), m = 6				LS (mesh dependent patch)
			Displacement		Stress		Displacement		Stress (ZZ)
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
113	274	2.04	0.89	0.939	2.25	1.437	1.10	1.064	0.81	0.937
269	610	1.40	0.71	0.876	0.98	1.094	0.84	0.949	0.73	0.885
504	1102	1.16	0.73	0.791	0.65	0.885	0.80	0.844	0.76	0.803

Displacement and stress recovery with weight functions and order of basis functions

The sensitivity of finite element error estimation with weight functions and order of basic functions, in MLS technique is also studied. The rate of convergence and effectivity obtained using MLS based error recovery analysis for linear triangular elements are shown in Table 4 considering linear (m = 3) and cubic (m = 9) order of basis function. The error estimation results with quadratic order of basis function (m = 6) are given in Table 1. The rate of convergence and effectivity obtained using meshfree MLS technique for triangular elements considering quartic spline (QS) and exponential (E) weight function are shown in Table 5. The error estimation results obtained with cubic spline (CS) weight function for triangular elements are given in Table 1. Table 5 depicts the percentage reduction in error in meshfree MLS based displacement/stress recovery considering various order of basis functions and weight functions for linear triangular meshing.

Table 4.

Error in energy norm and global effectivity (θ) for cylinder problem using displacement/stress recovery and order of basis function (linear triangular element).

Mesh		FEM error (×10−3)	Linear basis (m = 3, cubic spline)				Cubic basis (m = 9, cubic spline)
			Displacement		Stress		Displacement		Stress
No. of elem.	DOF		Error (×10−3)	Θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
125	162	4.05	4.06	1.147	0.39	1.262	1.20	0.960	1.44	0.952
272	324	2.77	2.22	1.085	2.05	1.158	0.75	0.964	0.87	0.953
550	622	2.00	1.40	1.038	1.10	1.037	0.68	0.945	0.69	0.935

Table 5.

Error (energy norm) and global effectivity (θ) for cylinder problem using displacement/stress recovery and weight functions (linear triangular element).

Mesh		FEM error (×10−3)	Quartic spline (m = 6)				Exponential (m = 6)
			Displacement		Stress		Displacement		Stress
No. of elem.	DOF		Error (×10−3)	Θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
125	162	4.05	1.78	0.943	2.06	1.008	1.25	0.922	2.22	0.729
272	324	2.77	0.89	0.941	1.02	0.972	0.82	0.939	1.45	0.743
550	622	2.00	0.69	0.929	0.66	0.941	0.72	0.932	1.05	0.722

Table 6.

Percentage reduction of error using MLS technique for cylinder problem in displacement/stress recovery with various order of basis function and weight function (linear triangular element, uniform mesh).

Mesh		%age reduction of error using MLS with order of basis function and weight function
		m = 3, CS		m = 6, CS		m = 9, CS		m = 6, QS		m = 6, E
No. of elem.	DOF	Displacement	Stress	Displacement	Stress	Displacement	Stress	Displacement	Stress	Displacement	Stress
125	162	–	–	149	110	238	181	128	97	224	82
272	324	25	35	226	186	269	218	211	172	238	91
550	622	43	82	199	203	194	190	190	203	178	91

Adaptive refinement using displacement and stress recovery based error estimation

The efficiency and reliability of displacement and stress based recovery technique under meshfree, and mesh dependent environment is assessed through adaptive finite element analysis of thick sheet cylinder problem that is, refinement of finite element mesh under guidance of displacement and stress error estimators for satisfying the predefined error limit. The initial meshes are adaptively improved to bring the solution error below the predefined error limit. Table 7 shows the global error, that is, summing up of the error contribution of all elements, number of elements, and degree of freedom (DOF) in improved meshes at predefined error limit of 4% using the meshfree MLS based error recovery and mesh dependent least square error recovery employing linear triangular elements uniform and non-uniform mesh and linear quadrilateral elements uniform mesh. MLS based error recovery technique considers cubic spline (CS) weight function and quadratic order of basis function (m = 6) in displacement and stress solution error recovery. The adaptively improved meshes at 4% target error using meshfree MLS, and mesh dependent least square based displacement and stress error estimation with linear triangular and quadrilateral elements are portrayed in Figures 6 to 8.

Table 7.

Actual and computed global errors (energy norm, %) and displacement/stress error recovery based adaptive finite element analysis of cylinder problem using different element types for 4% target error (cubic spline, m = 6).

Recovery type	Linear triangle (initial uniform mesh: 125 elements & 162 DOF, initial non-uniform mesh: 136 elements & 162 DOF)				Linear quadrilateral (initial mesh: 113 elements & 274 DOF)
	FEM error	Error	Adaptive mesh		FEM error	Error	Adaptive mesh
			N	DOF			N	DOF
MLS (Meshfree-CS-Disp.)-uniform	14.7	14.1	593	666	7.4	7.0	525	1154
MLS (Meshfree-CS-Disp.)-non-uniform	14.0	13.8	733	818
MLS (Meshfree-CS-Stress)-uniform	14.7	15.1	607	686	7.4	11.0	1223	2618
MLS (Meshfree-CS-Stress)-non-uniform	14.0	14.7	863	958
LS (Mesh Dependent-Displacement)-uniform	14.7	14.4	639	716	7.4	7.9	653	1414
LS (Mesh Dependent-Displacement)-non-uniform	14.0	13.7	831	920
LS (Mesh Dependent-Stress)-uniform	14.7	17.3	694	780	7.4	7.0	503	1100
LS (Mesh Dependent-Stress)-non-uniform	14.0	16.7	925	1020

Figure 6.

Adaptively improved mesh in cylinder example using meshfree MLS and mesh dependent least square (LS) based displacement and stress recovery techniques (initial triangular mesh elements (uniform) = 125, 4% target error): (a) initial mesh, (b) MLS-displacement (CS, m = 6), (c) MLS-stress (CS, m = 6), (d) LS-displacement, and (e) LS-stress (ZZ).

Figure 7.

Adaptively improved mesh in cylinder example using meshfree MLS and mesh dependent least square (LS) based displacement and stress recovery techniques (initial triangular mesh elements (non-uniform) = 136, 4% target error): (a) initial mesh, (b) MLS-displacement (CS, m = 6), (c) MLS-stress (CS, m = 6), (d) LS-displacement, and (e) LS-stress (ZZ).

Figure 8.

Adaptively improved mesh in cylinder example using meshfree MLS and mesh dependent least square (LS) based displacement and stress recovery techniques (initial quadrilateral mesh elements = 113, 4% target error): (a) initial mesh, (b) MLS-displacement (CS, m = 6), (c) MLS-stress (CS, m = 6), (d) LS-displacement, and (e) LS-Stress (ZZ).

Table 8 shows the global error, number of elements, and degree of freedom DOF in improved meshes at 4% error limit using the MLS error recovery technique employing quartic spline and exponential weight functions, linear and cubic basis function. The adaptively improved meshes at 4% error limit using meshfree MLS based displacement and stress error estimation employing different weight function and order of basis function are displayed in Figures 9 and 10.

Table 8.

Global errors using various order basis function and weight function in recovery analysis and number of elements (N) with DOF in adaptively improved meshes for 4% target error (initial linear triangular mesh: 125 elements and 162 DOF).

Weight/basis function	FEM error	MLS-displacement			MLS-stress
		Error	Adaptive mesh		Error	Adaptive mesh
			N	DOF		N	DOF
Linear basis/cubic spline	14.7	17.7	700	786	20.4	821	910
Cubic basis/cubic spline	14.7	14.4	668	748	14.4	522	592
Quartic spline/quadratic basis	14.7	14.2	623	700	15.4	582	658
Exponential/quadratic basis	14.7	13.8	592	666	10.9	434	502

Figure 9.

Adaptively improved mesh in cylinder example using meshfree MLS based displacement and stress recovery technique using different order of basis function (initial triangular mesh elements = 125, 4% target error, *8% error): (a) MLS-displacement (CS, m = 3), (b) MLS-stress* (CS, m = 3), (c) MLS-displacement (CS, m = 9), and (d) MLS-stress (CS, m = 9).

Figure 10.

Adaptively improved mesh in cylinder example using meshfree MLS based displacement and stress recovery technique using different weight function (initial triangular mesh elements = 125, 4% target error): (a) MLS-displacement (quartic, m = 6), (b) MLS-stress (quartic, m = 6), (c) MLS-displace (exponential, m = 6), and (d) MLS-stress (exponential, m = 6).

Benchmark example 2: Infinite plate with hole subjected to traction

An example involving square portion from an infinite plate having a circular hole with a radius (a) of one unit is also analyzed to compare the performance of stress and displacement recovery technique in stress concentration condition. The plate domain with hole is discretized with linear triangular/quadrilateral elements as shown in Figure 11. The one quarter of the plate domain is taken for analysis because of symmetry of centrally holed plate. The normal displacement component and shear stress are zero along the symmetry line. A unit in-plane traction is applied in the x-direction. The analytical solution of plate example ⁴ is given in equations (40)–(42).

$${\sigma }_x={\sigma }_{\infty }[1-\frac{a^2}{r^2}(1.5\cos 2\theta -\cos 4\theta )-1.5\frac{a^4}{r^4}\cos 4\theta ]$$

(40)

$${\sigma }_y={\sigma }_{\infty }[0-\frac{a^2}{r^2}(0.5\cos 2\theta -\cos 4\theta )-1.5\frac{a^4}{r^4}\cos 4\theta ]$$

(41)

$${\sigma }_{\mathrm{xy}}={\sigma }_{\infty }[0-\frac{a^2}{r^2}(0.5\sin 2\theta -\sin 4\theta )-1.5\frac{a^4}{r^4}\sin 4\theta ]$$

(42)

Figure 11.

Triangular/quadrilateral elements mesh for centrally holed plate domain: (a) triangular elements and (b) quadrilateral elements.

where r² = y² + x² and σ_∞ is the uniaxial traction applied at infinity.

Displacement and stress recovery with element type and patch configuration

The displacement and stress recovery based error estimations are obtained in the finite element solution of holed plate example with different elements discretization types that is, linear triangular, and quadrilateral element, and patch configurations that is, meshfree and mesh dependent patch, considering cubic spline weight function and quadratic order of basis function (m = 6). The properties of error estimation that is, error in energy norm and global effectivity are tabulated in Tables 9 and 10 with increasing mesh density for linear triangular and quadrilateral elements.

Table 9.

Error (energy norm) and global effectivity (θ) for holed plate example using displacement/stress recovery and patch configuration (linear triangular element, cubic spline).

Mesh		FEM error (×10−3)	MLS (meshfree patch), m = 6				LS (mesh dependent patch)
			Displacement		Stress		Displacement		Stress (ZZ)
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
155	194	12.91	9.28	0.763	12.37	0.908	11.20	0.805	13.23	1.018
548	618	8.50	4.59	0.878	6.70	0.983	5.50	0.907	7.91	1.126
1585	1700	5.17	2.20	0.891	3.10	0.946	2.57	0.900	4.67	1.149

Table 10.

Error (energy norm) and global effectivity (θ) for holed plate example using displacement/stress recovery and patch configuration (linear quadrilateral element, cubic spline).

Mesh		FEM error (×10−3)	MLS (meshfree patch), m = 6				LS (mesh dependent patch)
				Displacement		Stress		Displacement		Stress (ZZ)
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
179	414	8.09	5.64	0.863	12.75	1.507	7.29	0.940	9.26	1.401
578	1254	4.29	2.49	0.902	5.36	1.356	2.95	0.999	2.11	0.899
1842	3854	2.88	1.45	0.846	2.15	0.991	1.70	0.946	1.65	0.912

Displacement and stress recovery with weight functions and order of basis functions

The matchability of weight functions and appropriate order of poly nominal basis function are investigated for effectiveness of meshfree MLS based error recovery technique. The rate of convergence with increasing mesh order and effectivity obtained using MLS based error recovery analysis considering linear (m = 3) and cubic (m = 9) order of basis function for linear triangular elements meshing are shown in Table 11. The convergence rate and effectivity with increasing number of elements obtained using meshfree MLS technique for linear triangular elements mesh employing quartic spline and exponential weight function are shown in Table 12. The error estimation results obtained with cubic spline weight function and quadratic order of basis function (m = 6) for triangular elements are given in Table 1.

Table 11.

Error in energy norm and global effectivity (θ) for holed plate example using displacement/stress recovery and order of basis function (linear triangular element).

Mesh		FEM error (×10−3)	Linear basis (m = 3, Cubic Spline)				Cubic basis (m = 9, Cubic Spline)
				Displacement		Stress		Displacement		Stress (ZZ)
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
155	194	12.91	12.77	0.911	16.46	1.128	8.20	0.776	10.56	0.824
548	618	8.50	7.05	0.975	9.72	1.216	4.07	0.886	5.28	0.908
1585	1700	5.17	3.96	0.998	5.26	1.182	1.84	0.91	2.45	0.907

Table 12.

Error (energy norm) and global effectivity (θ) for holed plate example using displacement/stress recovery and weight functions (linear triangular element).

Mesh		FEM error (×10−3)	Quartic Spline (m = 6)				Exponential (m = 6)
			Displacement		Stress		Displacement		Stress (ZZ)
No. of elem.	DOF		Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
155	194	12.91	9.71	0.772	12.85	0.935	9.33	0.831	9.18	0.580
548	618	8.50	4.81	0.883	7.05	1.008	4.05	0.848	5.15	0.647
1585	1700	5.17	2.31	0.893	3.32	0.963	2.24	0.874	2.89	0.742

Adaptive refinement using displacement and stress recovery based error estimation

The adaptive finite element analysis of holed plate example is carried to assess the effectiveness and efficiency of displacement/stress based error estimators under meshfree and mesh dependent environment. The initial meshes are adaptively improved to bring the solution error below the preassigned accuracy limit. The global error, number of elements, and degree of freedom (DOF) in improved meshes at preassigned accuracy limit of 1% using the meshfree MLS based error recovery and mesh dependent least square error recovery employing linear triangular elements mesh are shown in Table 13. MLS based error recovery technique considers cubic spline (CS) weight function and quadratic order of basis function (m = 6) in displacement and stress solution error recovery. The adaptively improved meshes at 1% target error using meshfree MLS, and mesh dependent least square based displacement and stress error estimation with triangular and quadrilateral elements are portrayed in Figure 12.

Table 13.

Actual and computed global errors (energy norm, %) and displacement/stress error recovery based adaptive finite element analysis of holed plate example for 1% target error (initial uniform triangular mesh: 155 elements and 194 DOF, cubic spline, m = 6).

Recovery type	FEM error	Displacement			Stress
		Error	Adaptive mesh		Error	Adaptive mesh
			N	DOF		N	DOF
MLS (meshfree)	10.38	7.99	560	628	9.47	521	588
LS (mesh dependent)	10.38	8.42	468	532	10.46	528	596

Figure 12.

Adaptively improved mesh in holed plate example using meshfree MLS and mesh dependent least square (LS) based displacement and stress recovery techniques (initial triangular mesh elements (uniform) = 155, 1% target error): (a) initial mesh, (b) MLS-displacement (CS, m = 6), (c) MLS-stress (CS, m = 6), (d) LS-displacement, and (e) LS-stress (ZZ).

Discussion

The effectiveness and efficiency of MLS method with different weight functions and polynomial order of basis functions to recover the displacement/stress error in meshfree environment is studied. The comparison of meshfree MLS based displacement/stress error recovery techniques with mesh dependent least square based displacement/stress recovery techniques is also carried out. The MLS technique uses the nodal values in the local support domain to interpolate the field variable or field variable derivatives. The cylinder under pressure and centrally holed plate under traction examples are analyzed using finite element method coupled with error recovery techniques considering meshfree and mesh dependent local domains. The MLS error recovery techniques employ three weight functions, namely cubic spline, quartic spline, and exponential, and three set of polynomial order of basis functions. The problem domain is discretized using triangular and quadrilateral shape elements. The error recovery analysis is further utilized for adaptive improvement of domain mesh. The quality of discretization error obtained using recovery techniques is compared in terms of error convergence characteristics, error effectivity and adaptively improved element meshes with preassigned accuracy limit. The computation results obtained with increasing fineness order, element types, patch configuration, weight functions, and order of basis functions, for cylinder problem under the pressure, are depicted in Tables 1 to 6. The results of adaptive refined meshes for target accuracy obtained for cylinder example is given in Tables 7 and 8 and Figures 6 to 10. The computation results of error convergence and effectivity index obtained with increasing fineness order, element types, weight functions and order of basis functions for holed plate example, are tabulated in Tables 9 to 12. Table 13 and Figure 12 present the details of errors and adaptive refined meshes for 1% target error limit for holed plate example.

The influence of choice of element for discretization, weight functions, order of polynomial basis functions, dependency of node with mesh connectivity on finite element displacement, and stress error recovery in meshfree and mesh based recovery technique are compared in this work. From the Tables 1, 2, 9, and 10, showing error convergence and effectivity with mesh size obtained with bench mark examples analyzed with different element for discretization and patch configuration, it is clear that the order of errors are lower and effectivity is better in displacement recovery based error estimators as compared to stress recovery based error estimators with triangular element discretization while with quadrilateral element discretization, the stress recovery based error estimators have better error convergence and effectivity for both meshfree and mesh dependent patch configuration. From Table 6, it is observed that that the error reduction using meshfree MLS interpolation based displacement and stress recovery is about 10% and 150% respectively for the displacement and stress recovery over the mesh dependent least square based displacement and stress recovery. The reduction in errors using MLS interpolation based displacement and stress recovery considering linear triangular meshing are about (35%, 70%), (200%, 175%), and (250%, 200%) respectively for linear, quadratic, cubic orders of basis function and for cubic spline, quartic spline, and exponential weight functions are (200%, 175%), (175%, 150%), and (200%, 90%) respectively. The meshfree MLS interpolation recovery method is found to perform better in comparison to mesh dependent least square recovery method with different domain discretization. The improved results of the MLS interpolation recovery method over mesh dependent recovery methods seem to be due to better error recovery on boundary nodes in meshfree MLS interpolation recovery method. From the Tables 3 to 6, 11, and 12, showing error convergence and effectivity with mesh size obtained on cylinder and plate examples analyzed with different order of basis function and weight function in meshfree MLS recovery method, it is evident that lower order of basis function will not be effective in improving the finite element error results. However, increasing the order of basis function will enhance the error estimation quality that is, rate of convergence become faster with improved effectivity of the results. The increase in convergence rate with the increase of the order of basis function is more in displacement recovery technique as compared to stress recovery technique. It is observed that employing cubic spline weight in MLS interpolation recovery method has better quality of error recovery, for displacement and stress error recovery, than those employing other weight functions. It can be concluded that the order of basis function and weight function employed in MLS error recovery method has pronounced influence on error recovery quality.

The efficiency of error estimation using meshfree MLS displacement/stress error recovery technique and mesh dependent least square recovery technique is demonstrated in terms of adaptively improved meshes with specified error limits. The adaptive finite element results for cylinder problem are presented in Tables 7 and 8 and Figures 5 to 8 at specified error limit of 4%. The adaptively improved meshes are the indicators of the error distribution pattern in the domain as the meshes become finer in areas of high errors to get a uniform accuracy throughout the domain. The details of adaptive improved meshes at specified error limit of 1% for hole plate example are presented in Table 13 and Figure 12, The change in the error limits will change the order of mesh fineness with similar elements density pattern as shown in adaptive meshes figures. ¹⁷ It is observed from the results that the number of elements required to achieve objective accuracy, in general, is smaller using meshfree MLS based error estimators as compared to using mesh dependent error estimators. The increasing the order of basis function in MLS bases error estimator, improve the effectiveness of adaptive finite element analysis. It concludes that MLS based error recovery technique is more effective than the least square based recovery technique. It is also concluded that in general, the displacement recovery based error estimators is more effective than the stress recovery based error estimators.

Conclusions

The present study evaluates the effectiveness and reliability of field variable and field variable derivatives recovery based meshfree Moving Least Squares (MLS) interpolation approaches and mesh dependent least square approaches in finite element analysis with different mesh type, weight, and basis function. The recovery of field variable and field variable derivatives through MLS methods is carried out over patch of Meshfree scattered nodes (the patch of influenced nodes or support close to a given node up to a determined distance). The mesh dependent element patch (a patch of elements surrounding the given element) is used to recover the field variable using least square methods. The mesh dependent nodes patch (a patch of elements surrounding the given node) is also used to recover the field variable derivatives using least square (ZZ) methods. Finite element numerical experiments on cylinder and plate examples are carried out employing MLS interpolation and least square recovery techniques using triangular and quadrilateral element to recover the displacement/stress error. The characteristics of error estimation with different recovery techniques is compared in terms of convergence characteristics, effectivity and adaptively refined meshes. It is observed in the analysis of benchmark example with linear triangular meshing that the error reduction using meshfree MLS interpolation based displacement and stress recovery is about 10% and 150% respectively for the displacement and stress recovery over the mesh dependent least square based displacement and stress recovery. It is also observed from the analysis results that employing cubic spline weight in meshfree MLS interpolation recovery method has better behavior, for displacement and stress error recovery, than other weight functions used in MLS interpolation recovery method. It is also observed that by increasing the order of basis function will enhance the error estimation quality that is, rate of convergence become faster with improved effectivity of the results. The increase in convergence rate with the increase of the order of basis function is more in displacement recovery technique as compared to stress recovery technique. It is concluded from the study that the effectiveness and efficiency of meshfree displacement/stress error recovery technique strongly depends on the weight and basis functions used in Moving Least Squares (MLS) interpolation method.