Mesh-Free MLS-Based Error-Recovery Technique for Finite Element Incompressible Elastic Computations

Abstract

The finite element error and adaptive analysis are implemented in finite element procedures to increase the reliability of numerical analyses. In this paper, the mesh-free error-recovery technique based on moving least squares (MLS) interpolation is applied to recover the errors in the stresses and displacements of incompressible elastic finite element solutions and errors are estimated in energy norms. The effects of element types (triangular and quadrilateral elements) and the formation of patches (mesh-free patch, mesh-dependent element-based patch, and mesh-dependent node-based patch) for error recovery in MLS and conventional least-square interpolation-error quantification are also assessed in this study. Numerical examples of incompressible elasticity, including a problem with singularity, are studied to display the effectiveness and applicability of the mesh-free MLS interpolation-error recovery technique. The mixed formulation (displacement and pressure) is adopted for a finite element analysis of the incompressible elastic problem. The rate of convergence, the effectivity of the error estimation, and modified meshes for desired accuracy are used to assess the effectiveness of the error estimators. The error-convergence rates are computed in the original FEM solution, in the post-processed solution using mesh-free MLS-based displacement, stress recovery, mesh-dependent patch-based least-square-based displacement, and stress recovery (ZZ) as (0.9777, 2.2501, 2.0012, 1.6710 and 1.5436), and (0.9736, 2.0869, 1.6931, 1.8806 and 1.4973), respectively, for four-node quadrilateral, and six-node triangular meshes. It is concluded that displacement-based recovery was more effective in the finite element incompressible elastic analysis than stress-based recovery using mesh-free and mesh-dependent patches.

1. Introduction

The discretization of problem domains generates errors in the numerical technique. The finite element error and mesh adaptive analysis are implemented in finite element procedures to increase the reliability of finite element analyses. The recent research directions in finite element methods include the advancement of the finite element technique by overcoming its drawbacks and the enhancement of the method’s reliability and efficiency. Cen et al. [1] provide a survey of the finite element methods’ vast range of applications. Several techniques have been proposed to recover displacements or their gradients and to improve the finite element solution’s accuracy. A critical review of various error-estimation techniques to obtain the practical finite element results of linear and non-linear problems is summarized by Gratsch and Bathe [2]. New finite element techniques are developed to decrease the reliance on mesh for the analysis of problems. Chen et al. [3] summarized the surveys of mesh-free method developed to address the weaknesses of the finite element approach. Some the recent developments in alternative finite element methods (FEM) are the generalized finite element method [4], the extended finite element technique [5], the smoothed finite element method [6], the smoothed-particle-hydrodynamics method [7], the reproducing kernel particle method [8], and the moving-particle semi-implicit method [9]. A number of materials used in industry have very high bulk moduli in comparison with the shear modulus, i.e., nearly incompressible materials, and materials with large deformations may also be treated as incompressible. The displacement-based finite element approach has difficulties in the analysis of nearly incompressible materials, such as stiffness-matrix ill-conditioning, the locking phenomenon and spurious stresses. A summary of various methods used in the finite element analysis of incompressible and nearly incompressible elasticity is presented by Brink and Stein [10]. The developments over the past ten years in the application of various approaches in finite element methods to deal with volumetric locking are reviewed by Zeng and Liu [11]. Mixed Kirchhof stress–displacement–pressure formulation [12], a displacement-pressure formulation, and a simple pressure-projection-stabilized method [13] can be used to deal with incompressibility constraints. Doll et al. [14] investigated the effectiveness of selective reduced integration in overcoming the volumetric locking phenomenon in 2D and 3D solid elements. Numerous solutions to the problem of incompressible locking in displacement-based finite element approaches were been proposed by Boffi and Stenberg [15]. Nemer et al. [16] applied a stabilized finite element procedure using a mixed formulation in which the momentum equation of the continuum is extended by a pressure equation that treats the incompressibility constraints, in order to solve transient linear and nonlinear solid dynamics in compressible and incompressible materials. The finite-strain-incompressible-elasticity problem was analyzed by Gültekin et al. [17] using a variational technique based on a finite element method with a solution for volumetric constraints. The mixed meshless local Petrov–Galerkin approach was applied by Jarac et al. [18] to gradient elasticity, and they found high accuracy with a low order in the meshless approximation functions. The neural networks (ANN) technique was applied by Saikia et al. [19] to recuperate the stress errors in finite element solutions.

The numerical errors due to discretization in the finite element method can be minimized using the error estimation and adaptive-mesh-optimization techniques. In incompressible and large deformations such as those for rubber materials and metal-forming processes, the accuracy is also lost due to volumetric locking, in addition to errors due to discretization. Various recovery- and residual-based error estimators have been developed for utilization in adaptive-mesh-improvement strategies and for enhancing the accuracy and reliability of finite element method results [20]. Cai and Cai [21] proposed explicit residuals and an improved Zienkiewicz–Zhu (ZZ) error-estimator-based hybrid a posteriori error estimator, applicable in tandem with the finite element technique. Gabriel et al. [22] presented reliable and efficient residual-based a priori and a posteriori error estimators for mixed finite element methods, and the good performance of the proposed method was confirmed by solving the linear elasticity problem. Li et al. [23] proposed a mesh-free MLS approach with greater computational efficiency for the static and quasi-static analysis of thin-shell fractures. They utilized error estimations based on gradients to add linear reproducing points in locations needing refinement. A hp-adaptive computational analysis applicable to a variety of solid mechanics problems was carried out by Bird et al. [24]. Lancaster and Salkauskas [25] presented a MLS interpolation approach with high computational stability and accuracy, in which shape functions at nodes perform the interpolation activity. The effectiveness of the discrete stationary moving least-squares fitting approach, the continuous stationary moving least-squares approach, and the moving least-squares approach, suggested for the element-free Galerkin method, were evaluated by Lee and Zhou [26]. They found that that the performance of the discrete stationary moving least-squares approach was better than the other recovery approaches. The mixed-formulations-based finite element approach in conjunction with error estimators is utilized for the adaptive finite element simulation of linear incompressible elasticity problems [27,28]. Some a posteriori error estimators for nearly incompressible linear elasticity were developed by Kumor and Rademacher [29] using goal-oriented estimations based on the dual weighted residual approach. The application of the MLS-based mesh-less recovery technique for gradient-error recovery in elastic analysis was proposed by Ahmed et al. [30]. Karvonen et al. [31] performed a probabilistic error analysis on a Gaussian-process-based statistical finite element method (statFEM). The adaptive simulation of magnetized plasma transport in fusion reactors was carried out by Vogl et al. [32]. They employed Zienkiewicz–Zhu error estimation and quadrilateral elements in their study. The accuracy and efficiency investigations of finite element formulation and the non-polynomial trigonometric higher-order shear deformation theory (HSDT) used in the vibration analysis of bi-directional functionally graded plates with geometric imperfection and porosity was carried out by Katiyar et al. [33]. The accuracy-and-robustness study of finite element formulation and the improved first-order shear-deformation theory used in the bending/buckling analysis of bi-directional functionally graded plates with porosity was performed by Vinh et al. [34]. Kahla et al. [35] proposed a mesh-free radial point interpolation technique to recuperate the field-variable error in the finite element analysis of incompressible elastic problems and showed that error estimators can be successfully implemented for reliable adaptive discretization.

According to the review of the relevant literature, it is clearly evident that mesh-free recovery-based error estimation is a relatively new interest in finite element analysis, and few studies have investigated the different mesh-free procedures in mesh-free recovery-based error estimation specifically for incompressible elasticity and large domain changes. Therefore, there is a need to develop mesh-free recovery techniques using different mesh-free procedures, especially for problems of incompressible elasticity and large domain changes. The efficient and reliable error estimator based on the recently presented mesh-free radial point interpolation technique for the adaptive discretization of incompressible elastic problems motivated us to apply other mesh-free techniques for error-recover-based solutions and to present a comparison of mesh-free error-recovery results with mesh-dependent error-recovery results. In this study, a mesh-free interpolation error recovery technique is explored for incompressible elastic finite element analysis. Moving least squares (MLS) interpolation considering radial weights over circular support domains was implemented for the recovery of solution errors in incompressible elastic finite element analysis and the errors were quantified in an energy norm. The MLS interpolation technique utilizes the finite element method’s solutions in weighted least squares to create a continuous displacement/stress approximation. The displacement/pressure-based mixed formulation was used in the finite element analysis. The most common element types, six-node triangular and four-node quadrilateral elements, were employed for the mesh generation. Analyses of numerical examples, including problems with singularity, were performed to illustrate the effectiveness and efficiency of the meshless MLS interpolation approach for displacement/stress-solution-error recovery in incompressibility conditions. The precision of the error estimator was measured by its convergence rate and effectiveness in error estimation, and the meshes were updated to ensure the desired accuracy. The quality of the mesh-less MLS interpolation approach for displacement/stress error recovery was compared with the least-squares interpolation approach for displacement/stress-solution-error recovery in incompressibility conditions. The node-based and element-based node patches were employed for the displacement and stress, respectively, in the least-squares-interpolation-recovery approach. The mesh-dependent patches of elements for displacement recovery consider all the surrounding elements of the specified element and mesh-dependent patches of node for stress recovery consisted of the union of the nodes surrounding the specified node [20].

2. Finite-Element Incompressible Elastic Formulation

The mixed formulation is used for incompressible elastic problem analyses [35]. The basic governing equations for incompressible elastic formulation can be written as:

$$\begin{array}{l}\mathrm{Equilibrium}\mathrm{equation}\nabla \sigma +f=0in\mathsf{\Omega }\end{array}$$

(1)

$$\mathrm{Boundary}\mathrm{conditions}\sigma n=\stackrel{-}{t}on{\mathsf{\Gamma }}_t,\mathrm{and}u=\stackrel{-}{u}on{\mathsf{\Gamma }}_u,$$

(2)

where Ω is a problem domain, and $\overline{t}$ and $\overline{u}$ are given natural and essential boundary conditions on Γ_t, Γ_u boundaries respectively.

The stress (σ)- and strain (ε)-related equation for mixed formulations can be written as:

$$\mathrm{Constitutive}\mathrm{equation}\sigma \left(u\right)=2\mu \epsilon \left(u\right)+\lambda tr\left[\epsilon \right(u\left)\right]\mathit{I},\mathrm{or}$$

(3)

$$\sigma (u,p)=2\mu \epsilon +p\mathit{I},\mathrm{or}$$

(4)

$$\frac{p}{\lambda }=div\left(u\right),$$

(5)

where u = displacement, p = pressure.

Lame parameters:

$$\mathsf{\mu }=E/2(1+\nu )],\mathsf{\lambda }=E.\nu /\left[\right(1-2\nu \left)\right(1+\nu \left)\right],$$

(6)

The strain-displacement equation and the approximation of displacements (u) and pressure (p) in terms of nodal values ($\stackrel{-}{d},\stackrel{-}{p}$) of element using interpolation function (N) can be written as [36]:

$$\epsilon =Lu,u={\mathrm{N}}_u\stackrel{-}{d},p={\mathrm{N}}_p\stackrel{-}{p},$$

(7)

The system of equations obtained by applying the Galerkin method are given as:

$$\left[\begin{array}{cc}A& B\\ B^T& 0\end{array}\right]\left\{\begin{array}{c}\stackrel{-}{d}\\ \stackrel{-}{p}\end{array}\right\}=\left[\begin{array}{c}{f}_1\\ 0\end{array}\right]$$

(8)

where $\mathit{B}={\int }_{\mathit{\Omega }}^{}{\mathit{G}}^{\mathit{T}}\mathit{I}{\mathit{N}}_{\mathit{P}}\mathit{d}\mathit{\Omega },\mathit{A}={\int }_{\mathit{\Omega }}^{}{\mathit{G}}^{\mathit{T}}2\mathsf{\mu }\mathit{G}\mathit{d}\mathsf{\Omega },$ ${\mathbf{f}}_1={\int }_{{\Gamma }_{\mathbf{i}}}^{}{\mathbf{N}}_{\mathbf{u}}^{\mathbf{T}}{\stackrel{-}{\mathbf{t}}}_{\mathbf{i}}\mathbf{d}\Gamma +{\int }_{\mathsf{\Omega }}^{}{\mathbf{N}}_{\mathbf{u}}^{\mathbf{T}}\mathbf{f}\mathbf{d}\mathsf{\Omega }$, $\mathbf{G}=\mathbf{L}{\mathbf{N}}_{\mathbf{u}},$

3. Finite-Element-Solution Effectivity and Accuracy

The finite-element-analysis solution errors are the variation in the field variables or gradients in the finite-element=analysis results from those in the field variables or gradients’ post-processed results (or exact results). It is standard practice to measure the solution errors in terms of a norm such as energy or L₂ norms, giving the required information in scalar quantity. The error-estimation reliability is found by estimating problem effectivity. The ratio between the projected error and the exact error can be used to define the effectivity (θ). The effectivity for an asymptotically exact error estimator converges to one as the mesh size approaches zero [20].

$$\mathsf{\theta }=\frac{‖\mathrm{e}‖}{‖{\mathrm{e}}_{\mathrm{e}\mathrm{x}}‖},$$

(9)

where $‖e‖$ ($={\left[{{\int }_{\Omega }^{}{e}_{\sigma }^{*T}{D}^{-1}{e}_{\sigma }^{*}d\Omega }_{}^{}\right]}^{\frac{1}{2}}$) represents the evaluated error, and $‖e_{es}‖$ is the exact error (in energy norms).

The FEM solution’s accuracy (η) is as follows.

$$\mathsf{\eta }=\frac{‖{\mathit{e}}^{\mathit{*}}‖}{‖{\mathit{\sigma }}^{\mathit{*}}‖},$$

(10)

$$\mathrm{where}{‖{\sigma }^{*}‖}^2={‖{\sigma }^h‖}^2+{{‖e‖}_E}^2$$

(11)

4. MLS-Procedure-Based Mesh-Free Error Recovery

Due to its completeness and continuity [37], the moving least squares (MLS) approach can successfully interpolate data with reasonable accuracy. The MLS approach provides the displacement or stress estimate at a node via the interpolation of the displacement or stress in a weighted least-squares sense at the mesh-free local domains. The nodes can be expressed in 2-D as x1–xn, where x1 = (x1,y1). The displacement (or stress) estimate u^h(x) can be expressed as the product of the polynomial basis, P(x), and a coefficients vector, a(x).

$$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)$$

(12)

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

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

(13)

In this study, for linear element, m is taken to be six, and for quadratic elements, m is taken to be nine in the basis function P(x).

It is possible to obtain the vector of coefficients a(x) by minimizing a weighted residual in the following manner.

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

(14)

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

(15)

When the weighted residual is minimized, the relation of the coefficient vector is developed.

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

(16)

The MLS moment matrix, A(x), and B(x) are given as

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

(17)

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

(18)

The displacement (or stress) estimate may be expressed as:

$$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{\left(x\right)}_{jI}{u}_I$$

(19)

where u_s = nodal parameter vector of the displacement or stress, A(x) = MLS moment matrix, n = number of nodes i, and w(x − x_i) = weighting function in 2D associated with each node (weight takes a value of one in the vicinity of the node where the function and its derivatives are to be computed and becomes zero outside a region Ωi surrounding the node x_i)

The weights are created using the distance $\stackrel{-}{d}$ =‖x − x_i‖/d_m. The domain of influence for the radial weight is in the circular form, i.e., the support boundary is circular (Figure 1). The ‖x − x_i‖ is the node x’s distance from node x_i and d_m is the influence domain size of the node xi. The support size of the Ith node, d_mI, is computed by d_mI = d_max c_I. The value of c_I is calculated from distances of the neighbour nodes. The dilation parameter (d_max) is used as 3.0 in the analysis. The cubic spline weight function for MLS interpolation given below is considered in the present study.

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

(20)

The mesh-free MLS-based recovery estimation and mesh-improvement strategies were employed in two-dimensional FEM-based computer program. The mesh-dependent patch-based conventional error-recovery estimations and incompressible elastic formulation were also implemented. The program was run on i7-configuration computer with processor of 2.6 GHz and RAM of 16 GB to obtain the computational results. Flow chart for MLS-based recovery technique coupled with adaptive analysis of incompressible elastic problems is shown in Figure 2.

5. Illustrative Applications

5.1. Incompressible Infinite Square Elastic Plate

The effectiveness of the proposed mesh-free displacement or stress-recovery scheme is examined through the error-convergence rate of the error recovery of the displacement or stress in case of an incompressible elastic body under self-weights. The problem is used by Zienkiewicz et al. [36]. The example does not present singularity, so the theoretical convergence rate of the recovery scheme can be compared with other recovery schemes. The exact solutions for the example are included in Equations (21)–(25).

Domain and boundary conditions: Ω [0 × 0] × [1 × 1], u = v = 0 on Γ

Exact solutions:

$$u=2x^{2}y{\left(1-x\right)}^2(1-y)(1-2y)$$

(21)

$$v=-2{xy}^2(1-x)(1-2x){\left(1-y\right)}^2,$$

(22)

$$p=x^2-y^2,$$

(23)

Body forces:

$$b_x=4y\left(1-6x+6x^2\right)\left(1-3y+2y^2\right)+12x^2\left(1-2x+x^2\right)\left(-1+2y\right)-2x$$

(24)

$$b_y=-4x\left(1-6y+6y^2\right)\left(1-3x+2x^2\right)+12y^2\left(1-2y+y^2\right)\left(-1+2x\right)+2y,$$

(25)

The meshes generated using triangular and quadrilateral elements are shown in Figure 3. The plate problem was analyzed considering one-point reduced integration with quadrilateral elements (four-node) and two-point reduced integration with triangular elements (six-node) for volumetric strain term. The desired error limit for four-node quadrilateral and six-node triangular elements was kept as 2%. The finite-element-analysis results for convergence rate and effectivity of solution with different error-recovery approaches in energy norms considering meshfree patches and mesh-dependent patches are tabulated in

Table 1. Displacement/stress errors’ convergence and effectivity (θ) obtained with displacement/stress-recovery scheme considering mesh-free and mesh-dependent patches (four-node quadrilateral regular mesh).

Mesh Size (1/h)	Mesh-Free Patch					Mesh-Dependent Patch
	FEM Error (×10−3)	Displacement Recovery		Stress Recovery		Displacement Recovery		Stress Recovery (ZZ)
		Error (×10−3)	θ	Error (×10−3)	Θ	Error (×10−3)	θ	Error (×10−3)	Θ
1/4	29.44	15.87	0.8174	22.77	0.8890	18.23	0.8788	22.74	0.9573
1/16	7.69	0.66	0.9838	2.35	2.1200	2.05	0.9992	3.36	1.0264
1/32	3.85	0.15	0.9956	0.36	1.1855	0.57	1.0004	0.92	1.0080
Convergence Rate	0.9777	2.2501		2.0012		1.6710		1.5436

,

Table 2. Displacement/stress errors’ convergence and effectivity (θ) obtained with displacement/stress recovery scheme considering mesh-free and mesh-dependent patches (six-node triangular regular mesh).

Mesh Size (1/h)	Mesh-Free Patch					Mesh-Dependent Patch
	FEM Error (×10−3)	Displacement Recovery		Stress Recovery		Displacement Recovery		Stress Recovery (ZZ)
		Error (×10−3)	θ	Error (×10−3)	Θ	Error (×10−3)	θ	Error (×10−3)	Θ
1/4	40.20	23.89	1.0135	20. 52	1.0022	18.88	0.9304	21.58	0.8843
1/12	14.11	2.27	0.9879	2.98	0.9815	2.19	0.9655	4.22	0.9712
1/24	7.02	0.57	0.9929	0.99	0.9904	0.65	0.9715	1.48	0.9870
Convergence Rate	0.9736	2.0869		1.6931		1.8806		1.4973

,

Table 3. Displacement/stress errors’ convergence and effectivity (θ) obtained with displacement/stress recovery scheme considering mesh-free and mesh-dependent patches (four-node quadrilateral irregular mesh).

Elements	DOF	Mesh-Free Patch					Mesh-Dependent Patch
		FEM Error (×10−3)	Displacement Recovery		Stress Recovery		Displacement Recovery		Stress Recovery (ZZ)
			Error (×10−3)	θ	Error (×10−3)	Θ	Error (×10−3)	θ	Error (×10−3)	Θ
64	162	18.00	8.97	0.8602	8.60	0.8841	11.04	0.9555	14.32	1.0206
702	1514	5.02	0.96	0.9763	0.66	0.7432	1.25	1.0002	1.94	1.0004
1212	2566	3.67	0.47	0.9882	0.48	0.7812	0.69	1.0012	1.37	1.0134

and

Table 4. Displacement/stress errors’ convergence and effectivity (θ) obtained with displacement/stress recovery scheme considering mesh-free and mesh-dependent patches (six-node triangular irregular mesh).

Elements	DOF	Mesh-Free Patch					Mesh-Dependent Patch
		FEM Error (×10−3)	Displacement Recovery		Stress Recovery		Displacement Recovery		Stress Recovery (ZZ)
			Error (×10−3)	θ	Error (×10−3)	Θ	Error (×10−3)	θ	Error (×10−3)	Θ
28	146	35.23	23.54	0.5717	29.64	0.7046	22.62	0.5817	21.08	0.4905
340	1470	6.31	1.97	0.9077	2.47	0.9008	1.57	0.8710	2.04	0.8347
580	2462	2.20	1.13	0.8750	1.16	0.8889	1.21	0.8118	1.31	0.7814

. The updated mesh plots for desired accuracy of 2% with various recovery procedures and meshing schemes is given in Figure 4 and Figure 5. The element numbers (N) and degrees of freedom (DOF) after mesh update for desired accuracy of 2% in adaptive analysis with displacement/stress-error-recovery techniques for quadrilateral/triangular discretization are presented in

Table 5. Updated mesh DOF and element numbers (N) in incompressible elastic plate adaptive analysis with displacement/stress-recovery techniques considering quadrilateral/triangular elements and mesh-free/mesh-dependent patches (2% desirable error rate).

Meshing Scheme	Initial Mesh		Mesh-Free Patch				Mesh-Dependent Patch
			Displac. Recovery		Stress Recovery		Displac. Recovery		Stress Recovery
	N	DOF	N	DOF	N	DOF	N	DOF	N	DOF
Quadrilateral	64	162	9209	19,766	2222	4584	11,283	22,926	14,658	29,730
Triangular	28	146	1439	5936	791	3280	2577	10,540	1767	7264

.

Effect of Patch Formation on Error Recovery

In order to study the effect of the patch formation on the MLS-based recovery, the mesh-dependent-element patch for the MLS-based displacement-recovery scheme was also considered with four-node and six-node elements in the incompressible-elastic-finite-element analysis. The analysis results for the MLS interpolation using mesh-dependent patches are given in

Table 6. Displacement errors; convergence and effectivity (θ) in incompressible plate with displacement recovery considering MLS-interpolation recovery with mesh-dependent patches (four-node quadrilateral and six-node triangular mesh).

Mesh-Dependent Patch-MLS (4-Node)				Mesh-Dependent Patch-MLS (6-Node)
Mesh Size (1/h)	FEM Error (×10−3)	Displacement Recovery		Mesh Size (1/h)	FEM Error (×10−3)	Displacement Recovery
		Error (×10−3)	Θ			Error (×10−3)	Θ
1/4	29.44	16.54	0.7234	¼	40.20	18.71	0.9409
1/16	7.69	0.99	0.9677	1/12	14.11	2.94	0.9811
1/32	3.85	0.24	0.9909	1/24	7.02	0.89	0.9859
Convergence Rate	0.9777	2.0253			0.9736	1.7001

.

5.2. Infinite Incompressible Elastic Plate with Rigid Inclusion

The characteristics of the mesh-free displacement/stress-recovery schemes were also investigated with the help of another example of an infinite incompressible elastic plate with rigid inclusion. An adaptive analysis and automatic mesh-update strategies were used in the example. The exact solution to the incompressible elasticity example is known [38]. The solution derivatives are continuous on the boundary and in the inside domain but with a single point or points lying outside of the domain. The exact displacements and stresses are given by the following equations.

$$u_r=\left(\frac{T_x}{8Gr}\right)\left\{\right(k-1\left)r^2\right]+2\gamma R^2+[\beta \left(k-1\right)R^2+2r^2-2\delta \left(\frac{R^4}{r^2}\right)]\}cos2\theta ,$$

(26)

$$u_{\theta }=-\left(\frac{T_x}{8Gr}\right)[\beta \left(k-1\right)R^2+2r^2-2\delta \left(\frac{R^4}{r^2}\right)]sin2\theta ,$$

(27)

$${\sigma }_{rr}=(T_x/2)[1-\gamma (R^2/r^2\left)\right]+(T_x/2)[1-2\beta (R^2/r^2)+3\delta (R^4/r^4\left)\right]cos2\theta ,$$

(28)

$${\sigma }_{\mathrm{\theta \theta }}=\left(\frac{T_x}{2}\right)\left[1+\gamma \left(\frac{R^2}{r^2}\right)\right]-(T_x/2)[1-3\delta (R^4/r^4\left)\right]cos2\theta ,$$

(29)

$${\tau }_{r\theta }=-(T_x/2)[1+\beta \left(\frac{R^2}{r^2}\right)+3\delta (R^4/r^4\left)\right]sin2\theta ,$$

(30)

where r = y² + x², T_x is the uniaxial traction applied at infinity, and, for the plain strain case, the constants k, β, γ, and δ can be written in terms of the Poisson ratio υ as k = 3 − 4υ, β =−2/(3 − 4υ), γ =−(2 − 4υ)/2, δ = 1/(3 − 4υ).

The problem domain is defined by radius (R) = 1 unit, sided (w,b) = 4 unit, as shown in Figure 6. The rigid inclusion center is a singular point. The arc AE boundary has both x- and y-displacement as zero. The lines of the CB and DE boundary have zero shear stress and normal displacement components. The tractions computed from Equations (28)–(30) are imposed along the boundaries of AB and CD. The corresponding load vectors were computed by numerical quadrature using twelve Gauss points per element side. The plate problem was analyzed using quadrilateral elements with one-point reduced integration and six-node triangular elements with two-point reduced integration for the volumetric strain term (Figure 7). The desired accuracy in the energy norms for the four-node and six-node elements was kept as 1%. The computational results for the convergence rates and the effectivity for the different error-recovery techniques for the mesh-free patch and standard patch are presented in

Table 7. Displacement/stress errors’ convergence and effectivity (θ) in incompressible plate with rigid inclusion in displacement/stress-recovery scheme considering mesh-free and mesh-dependent patches (four-node quadrilateral mesh).

Elements	DOF	Mesh-Free Patch					Mesh-Dependent Patch
		FEM Error (×10−3)	Displacement Recovery		Stress Recovery		Displacement Recovery		Stress Recovery
			Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	θ
266	594	157.20	90.10	0.8989	123.04	1.7373	123.63	1.1371	148.44	1.1276
793	1698	91.17	31.40	0.9212	43.18	1.7426	51.56	1.0264	65.40	1.0214
1842	3854	67.81	22.71	0.8972	19.45	0.9901	26.84	0.9587	37.62	0.9270

and

Table 8. Displacement/stress errors’ convergence and effectivity (θ) in incompressible plate with rigid inclusion in displacement/stress-recovery scheme considering mesh-free and mesh-dependent patches (six-node triangular irregular mesh).

Elements	DOF	Mesh-Free Patch					Mesh-Dependent Patch
		FEM Error (×10−3)	Displacement Recovery		Stress Recovery		Displacement Recovery		Stress Recovery
			Error (×10−3)	θ	Error (×10−3)	θ	Error (×10−3)	Θ	Error (×10−3)	θ
124	558	150.98	104.88	0.7295	93.95	0.6296	133.92	0.9637	137.51	0.8937
373	1604	73.53	40.66	0.7791	40.66	0.6626	65.97	1.0426	67.65	0.9273
874	3666	40.52	36.41	1.0168	25.23	0.6304	21.14	0.8994	23.00	0.8181

. The updated mesh plots for the desired accuracy of 1% with various recovery procedures and meshing schemes are shown in Figure 8 and Figure 9. The degrees of freedom and element numbers (N) after the mesh improvement to ensure the desired accuracy of 1% using displacement/stress-recovery techniques and quadrilateral/triangular mesh in the adaptive analysis of the incompressible plate with rigid inclusion are presented in Table 8.

6. Discussion

The use of the moving least squares (MLS) interpolation-based error-recovery approach to increase the reliability of incompressible elastic analysis is introduced in this study. The computational results were obtained for error quality, i.e., convergence rates, effectivity, and updated meshing for the intended accuracy by analyzing the incompressible plate problems. Moving least squares (MLS) interpolation considers radial weights over circular support domains and errors are estimated in energy norms. The incompressible constraints are implemented through a displacement/pressure-based mixed approach. The plate domains are discretized using triangular and quadrilateral elements. The incompressible elastic benchmark examples, including problems with singularity, were studied to demonstrate the effectiveness and efficiency of mesh-free MLS interpolation for displacement/stress-error recovery, and of mesh-dependent MLS interpolation for displacement error recovery. The MLS interpolation-based error-recovery results were also compared with the least-squares interpolation-based error-recovery results. The mesh-dependent patch formation for displacement recovery considers all the elements surrounding the specified element, and mesh-dependent patch formation for stress recovery considers the union of the nodes surrounding the specified node. The target error limit in the energy norm is kept as 2% in the benchmark elastic plate and 1% for the benchmark plate with rigid circular inclusion. The incompressible elastic analysis results for error convergence and effectivity for various recovery procedures considering mesh-free and mesh-dependent patches are tabulated in Table 1, Table 2, Table 3, Table 4 and Table 6.

For both the four-node and the six-node elements, the convergence obtained with the help of the mesh-free MLS-based error recovery was found to be better than that obtained with the mesh-dependent error recovery techniques. This indicates the higher efficiency of the mesh-free error-estimation scheme. In the mesh-dependent stress-recovery technique, there is a chance of a reduction in the precision of the recovery for the boundary nodes, since fewer points are available. The mesh-free recovery technique eliminates such difficulties, since more points are available in the support domain (mesh-free patch). It is evident from the tables that the error-recovery performance of the four-node-element discretization of the incompressible elastic problem was better than that of the performance of the six-node-element discretization. The tables also show that the performance of the displacement-based recovery was better than that of the stress-based recovery for both mesh-free and mesh-dependent patch formation. The error-convergence rates were computed in the original FEM solution, the post-processed solution using the mesh-free MLS-based displacement/stress-recovery scheme, and the mesh-dependent patch LS-based displacement/stress-recovery scheme, yielding values of (0.9777, 2.2501, 2.0012, 1.6710, and 1.5436) and (0.9736, 2.0869, 1.6931, 1.8806, and 1.4973) for the four-node quadrilateral and six-node triangular meshes, respectively.

The mesh-dependent patch (mesh-dependent-element patch) for the MLS-based displacement-recovery scheme was also considered with four-node and six-node elements in order the study the effect of patch configuration on MLS-based recovery. The analysis results for the standard patch for the MLS interpolation are given in Table 6. It is evident from the tables that the recovery-based error estimations were affected by the patch formation and that the mesh-dependent MLS interpolation had a lower convergence rate than and was effective compared to the mesh-free MLS interpolation for solution-error recovery. The convergence rates in the mesh-less and mesh-dependent MLS interpolation for field-variable-error recovery were 2.2501 and 2.0253 and 2.0869 and 1.7001, respectively, for the four-node quadrilateral and six-node triangular meshes However, the MLS interpolation-based error recovery using mesh-dependent patches provided better error-recovery quality than the least-squares interpolation-based error recovery using mesh-dependent patches. Table 7 and Table 8 show the error convergence and effectiveness of the displacement/stress-recovery procedures considering mesh-free and mesh-dependent patches in the incompressible elastic plate problem with rigid circular inclusion. The computational results showed similar trends for the incompressible elastic plate problem, i.e., the performance of the displacement-based recovery was better than that of the stress-based recovery for both mesh-free and mesh-dependent patch formation. The error-recovery performance was better for the quadrilateral elements than for the triangular elements selected for the meshing of the problems.

In all the adaptive finite element analyses, irregular meshes were considered. The robustness of the mesh-free error-estimation technique was also tested using irregular meshes, so initial irregular meshes were considered in the adaptive finite element analysis. The updated mesh plots for the desired accuracy with various recovery procedures and meshing schemes are shown in Figure 4, Figure 5, Figure 8 and Figure 9. The DOF and element numbers after the mesh improvement for the desired accuracy with the displacement and stress-recovery techniques and quadrilateral/triangular meshes are presented in Table 5 and

Table 9. Updated mesh-element numbers (N) and degrees of freedom (DOF) in incompressible plate with rigid circular inclusion adaptive analysis with displacement/stress-recovery techniques considering quadrilateral/triangular elements and mesh-free/mesh-dependent patches (1% desirable error).

Meshing Scheme	Initial Mesh		Mesh-Free Patch				Mesh-Dependent Patch
			Displac. Recovery		Stress Recovery		Displac. Recovery		Stress Recovery
	N	DOF	N	DOF	N	DOF	N	DOF	N	DOF
Quadrilateral	266	594	694	1490	419	926	760	1622	609	1314
Triangular	124	558	151	676	121	548	218	962	188	842

. It was observed from the mesh plots that the mesh-free error-estimation scheme was more efficient than the mesh-dependent patch-based error estimation. It was inferred that in general, the number of elements (N) necessary to attain the desired accuracy with the displacement-error-recovery scheme was lower than with the stress-error-recovery scheme with both mesh-free and standard patches and with all the types of element mesh. However, the lowest number of elements was required for the desired accuracy in the mesh-free stress-recovery technique, which suggests that this technique is highly effective in the recovery of stress. It can be concluded that the application of the mesh-modification procedure under the guidance of an effective error estimator can predict the distribution of errors in the domains and zones of high gradients.

7. Conclusions

The present study explored the effectiveness and applicability of the mesh-free MLS interpolation-approach-based error recovery of displacement/stress in incompressible elastic finite element analysis. The moving least-squares (MLS) interpolation technique considering radial weights over circular support domains was employed for the recovery of solution errors and the errors were quantified in energy norms. The incompressible elastic plate problems, with six-node triangular and four-node quadrilateral discretization schemes, were solved to investigate the order of the errors, the rate of the error convergence, the effectiveness (θ), and the updated meshes for the desired error limit. It was found from the computational results that the recovered field-variable type, the types of elements, and the patch formation had strong effects on the quality of the solution-error recovery. The error convergence obtained with the help of the mesh-free recovery was found to be better than that obtained with the mesh-dependent error-recovery technique. It was also found that the error-recovery performance was better for the quadrilateral elements than for the triangular elements selected for the meshing of the problems. The error-convergence rates were computed in the original FEM solution, in the post-processed solution using mesh-free MLS based displacement, and in the mesh-dependent patch-based least-squares-based displacement, and the stress-recovery values were (0.9777, 2.2501, 2.0012, 1.6710, and 1.5436) and (0.9736, 2.0869, 1.6931, 1.8806, and 1.4973) for the four-node quadrilateral and six-node triangular meshes, respectively. These illustrative examples demonstrate the mesh-free MLS based error recovery technique’s ability to extract errors effectively and efficiently in the finite element analysis of incompressible elastic problems.