Strain smoothing for compressible and nearly-incompressible ﬁnite elasticity

.


Introduction
Low-order simplex (triangular or tetrahedral) finite element methods (FEM) are widely used because of computational efficiency, simplicity of implementation and the availability of largely automatic mesh generation for complex geometries.However, the accuracy of the low-order simplex FEM suffers in the incompressible limit, an issue commonly referred to as volumetric locking, and also when the mesh becomes highly distorted.
To deal with these difficulties various numerical techniques have been developed.A classical approach is to use hexahedral elements instead of tetrahedral elements due to their superior performance in plasticity, nearly-incompressible and bending problems, and additionally their reduced sensitivity to highly distorted meshes.However, automatically generating high-quality conforming hexahedral meshes of complex geometries is still not possible, and for this reason it is desirable to develop improved methods that can use simplex meshes.Significant progress has, however, been done in this direction [1].
Another option is to move to higher-order polynomial simplex elements.While they are significantly better than linear tetrahedral elements in terms of accuracy this is at the expense of increased implementational and computational complexity, and sensitivity to distortion.
Nodally averaged simplex elements [2,3] can effectively deal with nearly-incompressible materials, but they still suffer from an overly stiff behaviour in certain cases [4].
Meshfree (or meshless) methods [5][6][7] are another option because of their improved accuracy on highly-distorted nodal layouts, but the locking problem is still a challenging issue that needs careful consideration [8].To improve the non-mesh based methods, B-bar approach [9,10], which is appropriate not only to handle incompressible limits but also to model shear bands with cohesive surfaces, can be considered.Additionally, because they are substantially different to the FEM, they are not easily implemented in it existing software.
Isogeometric Analysis (IGA) is another high-order alternative and the interested reader is referred to [11,12].Moreover for the further studies for fractures undergoing large deformations, edge rotation algorithm can be an another option in large plastic strains [13,14].
Mixed and enhanced formulations are another popular remedy for volumetric locking [15,16], but they retain the sensitivity to mesh distortion of the standard simplex FEM [17].
Another approach, and the one that we employ in this paper, is the strain smoothing method developed by Liu et al. [18,19].The strain smoothing method has the advantage over the above methods that it improves both the behaviour of low-order simplex ele-ments with respect to both volumetric locking and highly distorted meshes, while being simple to implement within an existing finite element code.
The basic idea of strain smoothing is based on the stabilised conforming nodal integration (SCNI) proposed in the context of meshfree methods by Chen et al. [20,21].Later SCNI was extended to the natural element method (NEM) by Yoo et al. [22], and was shown to effectively handle nearly-incompressible problems.
In the smoothed finite element method (S-FEM), the domain is divided into smoothing domains where the strain is smoothed as shown in Fig. 1.Typically, the geometry of the smoothing domains is derived directly from the standard simplex mesh geometry.Then with the divergence theorem, numerical integration is transferred from the interior to the boundary of the smoothing domains [23,24].Critically, this procedure results in a discrete weak form without the Jacobian, the matrix used to map basis function derivatives from the reference element to the real element in the mesh.In the standard FEM the Jacobian is required to construct the derivatives of the basis functions.When distorted meshes are used in the standard FEM, the Jacobian becomes ill-conditioned, and this affects the accuracy of the method.Because the Jacobian is not required in S-FEM, the resulting method is significantly more robust than the standard FEM on highly distorted meshes.
It is also known that the S-FEM produces stiffness matrices that are less stiff than the standard FEM, and in certain cases this property can be used to overcome volumetric locking.Since S-FEM was introduced, its properties have been studied from a theoretical viewpoint [18,19,[25][26][27][28][29], extended to n-sided polygonal elements [30] and applied to many engineering problems such as plate and shell analysis [31][32][33][34].
Particularly, Bordas et al. [35] recalled the central theory and features of S-FEM and showed notable properties of S-FEM which depend on the number of smoothing domains in an element.Moreover, Bordas et al. [35] presented the coupling of strain smoothing and partition of unity enrichment, so called SmXFEM, with examples of cracks in linear elastic continua and arbitrary cracks in plates.
The contribution of this paper to the literature is to present a robust, efficient and stable form of the smoothed finite element methods to simulate both compressible and nearly-compressible hyperelastic bodies.We study two forms of smoothing (nodebased and edge-based) and compare their relative merits.A key ingredient of our method is to add cubic bubbles to each element to ensure inf-sup stability.Although bubbles have been suggested before in the context of linear elastic S-FEM by Nguyen-Xuan and Liu [36] here we make the non-trivial extension to deal with hyperelastic problems.Finally we present a rigorous testing procedure that demonstrates the superior performance of our approach over the standard FEM.
The outline of this paper is as follows; first, we briefly review the idea fundamentals of S-FEM.In Section 3 we formulate the non-linear S-FEM for hyperelastic neo-Hookean compressible materials.To demonstrate the accuracy and convergence properties of the proposed methods we present extensive benchmark tests in Section 4. Finally, conclusions and future work directions are summarised in Section 5.

Smoothed finite element method (S-FEM)
It was shown in numerous studies that S-FEM provides a higher efficiency, i.e. computational cost versus error than the conventional FEM for many mechanical problems.We list below some of the strengths and weaknesses of each variant: the cell-based smoothed FEM (CS-FEM), the edge-based smoothed FEM (ES-FEM), the node-based smoothed FEM (NS-FEM), and the face-based smoothed FEM (FS-FEM).
Volumetric locking.NS-FEM can handle effectively nearlyincompressible materials where Poisson's ratio v ?0.5 [37], while ES-FEM suffers from volumetric locking.Combining NSand ES-FEM gives the so-called the smoothing-domain-based selective ES/NS-FEM which also overcomes volumetric locking [38].In the case of CS-FEM, volumetric locking can be avoided by selective integration [39].Upper and lower bound properties.In typical engineering analysis with homogeneous Dirichlet boundary conditions the NS-FEM gives upper bound solution and FEM obtains lower bound solution in the energy norm.While, in the case of problem with no external force but non-homogeneous Dirichlet boundary conditions, NS-FEM and FEM provide lower and upper bounds in the energy norm, respectively [40,41].Static and dynamic analyses.ES-FEM gives accurate and stable results when solving either static or dynamic problems [42].I n contrast, although NS-FEM is spatially stable, it is temporally unstable.Therefore, to solve dynamic problems, NS-FEM requires stabilisation techniques [43,44].CS-FEM can also be extended to solve dynamic problems [45].Other features.In NS-FEM, the accuracy of the solution in the displacement norm is comparable to that of the standard FEM using the same mesh, whereas the accuracy of stress solutions  in the energy norm is superior to that of FEM [38].In terms of computational time, in general, ES-FEM is more expensive than conventional FEM on the same mesh [38].

Non-linear elasticity and S-FEM approximation
The principle of virtual work for finite elasticity can be written in the Galerkin weak form [46][47][48]: where the smoothed deformation gradient F ¼ I þ ru is written in terms of displacements u, v is the set of admissible test functions.The strain energy density function W for a compressible neo-Hookean material [49] is: 3 l, and the shear modulus l > 0 and the bulk modulus k > 0 are material parameters.
The smoothed deformation gradient F for the proposed technique is: where the deformation gradient F is given in Appendix A.
To find an approximate solution using Eq. ( 2) for the displacement field u, we employ the Newton-Raphson method.At iteration iter + 1, knowing the displacement u iter from iteration iter, find r iter that satisfies [46]: and i; j; k; l 2f1; 2g for two dimensional problems.The energy function Eq. ( 5) and its directional derivatives Eq. ( 6) become the following equivalent formulations, respectively: where i; j; k; l; p; s 2f1; 2g.The resulting algebraic system for the numerical approximation of Eq. ( 4) is assembled from the block systems: By taking v ¼ P I N I v I , we obtain the stiffness matrix Kiter with following components: and the components of the load vector are: The smoothed tangent stiffness Ktan ¼ Kmat þ Kgeo can be re-written using Eq. ( 10): where the smoothed strain-displacement matrices B0 and B can be expressed respectively as (also see in Fig. 2) and by Eq. ( 11) the load vector b is: where matrix S is: where the fourth-order elasticity tensor C is: Finally, the global system of equations Eq. ( 4) can be written as: and 3. Enriched strain smoothing method with bubble functions In the finite element method, to apply the Ritz-Galerkin method to a variational problem, a finite dimensional sub-space of space V is required.The space V defined on domain X is approximated by simple functions which are polynomials [50].
V ¼fu 2ðH1 ðXÞÞ where displacement u,b o u n d a r yC and a Hilbert space H 1 ðXÞ.I nt h i s space, we cannot avoid the locking phenomenon in the incompressible limit, and S-FEM may face this obstacle as well because in both FEM and S-FEM, the same low-order simplex elements are used.One popular technique to overcome the locking effects is employing bubble functions within mixed finite element approximation [51,52].Nguyen-Xuan and Liu [36] proposed a bubble enriched smoothed finite element method called the bES-FEM (see also [53]).In addition, further studies of bubble functions are used in mixed finite strain plasticity formulation with MINI element for quasi-incompressible plasticity fractures [54], and brittle and ductile models [14].
A bubble function supplements an additional displacement field at a node placed at centroid of triangle T. In contrast to the MINI element, ES-FEM constructs a displacement-based formulation.ES-FEM with a bubble function has only a linear displacement field as unknown which has value one at the centroid of triangle T and the pressure vanishes at the edges of triangle T. As shown in Fig. 3, and interior node is located at the geometric centre with an additional displacement field associated with the cube bubble.
The cubic bubble function introduced in [55] is used in this paper.Since the first three basis functions are not zero at the centroid (1/3, 1/3), a basis function Wðn; gÞ¼½1 À n À g; n; g; 27ngð1 À n À gÞ T is necessarily required transformation form gives as: Wðn; gÞ T ¼ Wðn; gÞ T B À1 S ¼½1 À n À g; n; g; 27ngð1 À ngÞ 100 0 010 0 001 0 and therefore the basis functions become as: Wðn; gÞ¼ The properties of renewed basis functions and cubic bubble function of a right 45°three-node triangular element are given as (also see in Fig. 4):

Simple shear deformation
For simple shear deformation, the deformation gradient takes the form: where k > 0. For this deformation, the strain invariants are: Thus the incompressibility condition is always satisfied regardless of the material characteristics (isochoric deformation).Substituting this in Eq. ( 2) gives the following strain energy function: Fig. 2. The integration is performed on Gauss points located at the mid-point of the boundaries C k of the smoothing domain X k .
The non-zero entries of the corresponding Cauchy stress tensor are [59,46]: Hence Eq. ( 27) can be written: For this section, the shear and bulk moduli used are l ¼ 0:6 and j ¼ 100, respectively.The higher value of j, the material is more incompressible.
Dirichlet boundary conditions.To obtain the simple shear of a square section as shown in Fig. 5, the following Dirichlet boundary conditions can be imposed: All edges: ðu 1 ; u 2 Þ¼ðkX 2 ; 0Þ.Fig. 6 illustrates the deformed shape of the standard FEM and the proposed technique for the simple shear deformation with Dirichlet boundary conditions when the deformation is k = 1 for both the FEM and the S-FEM.
The strain energies for the analytical, FEM and ES-FEM solutions are shown in Table 1.The analytical solution can be calculated by Eq. ( 26) and is such that W¼0:3.
Table 1 provides the values of the relative error in strain energy for FEM, ES-FEM and NS-FEM.The values of the proposed formulations are within machine precision for moderate and coarse meshes.

Pure shear deformation
In this section pure shear deformation is considered, the deformation of pure shear is given as [46,60]: and therefore the deformation gradient for pure shear F is: Therefore the left Cauchy-Green tensor B is: The Cauchy stress is: Mixed boundary conditions.The deformed shape of the approach for pure shear with the mixed Neumann and Dirichlet boundary conditions are shown in Fig. 7.

Uniform extension with lateral contraction
We deform a 3D sample of compressible material in Eq. ( 24) by the following triaxial stretch: where X =[X 1 , X 2 , X 3 ] T and x =[x 1 , x 2 , x 3 ] T denote the reference (Lagrangian) and current (Eulerian) coordinates, respectively, and k i > 0, i = 1, 2, 3, are positive constants.The corresponding deformation gradient is: and the left Cauchy-Green tensor is B = FF T .We can then calculate the strain invariants using the following formulae: For the triaxial deformation, the strain invariants are: In particular, if the deformation is isochoric (preserves volume), then I 3 =1.The biaxial deformation associated with a square section of the material is then obtained by setting k 3 ¼ 1.In this case, if the deformation is isochoric, then k 2 ¼ 1=k 1 , and the strain invariants are: Substituting these in Eq. ( 2) gives the following value for the strain energy function: By the Rivlin-Ericksen representation, the Cauchy stress takes the general form: where the elastic response coefficients are calculated as follows: In particular, for the biaxial deformation of the square material, the non-zero components of the Cauchy stress are: Hence, the non-zero components of the Cauchy stress tensor are: Dirichlet boundary conditions.To obtain the above biaxial stretch of a square section, assuming that the sides of the square are aligned with the directions X 1 and X 2 , and the bottom left-hand corner is at the origin O (0, 0), then the following Dirichlet boundary conditions can be imposed: Bottom edge: ðu 1 ; u 2 Þ¼ððk 1 À 1ÞX 1 ; 0Þ; Left-hand edge: ðu 1 ; u 2 Þ¼ð0; ð1=k 1 À 1ÞX 2 Þ; Top and right-hand edge: The deformed shapes for the uniform extension with lateral contraction with Dirichlet boundary conditions are illustrated in Fig. 8.The relative strain energy errors are shown in Table 2. 2  Mixed boundary conditions.Alternatively, Neumann boundary conditions can be imposed on some of the edges.Before we can do this, we need to recall the general formula for the first Piola-Kirchhoff stress tensor: Then, for the biaxial stretch with k 2 ¼ 1=k 1 and k 3 ¼ 1, we obtain the following non-zero components for this tensor: At the corners, if one of the adjacent edges is subject to Dirichlet conditions and the other to Neumann conditions, the Dirichlet conditions are essential and take priority over the Neumann conditions.If both edges are subject to Neumann conditions, these are to be imposed simultaneously at the corner.Fig. 9 represents the deformed shapes with mixed boundary conditions, and the relative errors for this problem are given in Table 3.Note that all methods provide, again, the exact results down to machine precision.

''Not-So-Simple" shear deformation
Consider now the non-homogeneous deformation: for which the deformation gradient is: where k >0.For clarity of presentation, denote K =2kX 2 .Then the strain invariants are:   Strain energy relative error is given by: WNumerical ÀWExact WExact Â 100%. 2 We observe that all methods provide the exact results at machine precision.
and substituting Eq. ( 49) in Eq. ( 2) gives the strain energy function: Note that this function is not constant.
Dirichlet boundary conditions.To obtain the simple shear of a square section, the following Dirichlet boundary conditions can be imposed (see Fig. 10): All edges: ðu 1 ; u 2 Þ¼ðkX 2 2 ; 0Þ The deformed shape of the ''Not-so-simple" shear deformation is illustrated in Fig. 11.The strain energy relative errors for FEM and S-FEM are given in Table 4.The results of the FEM are comparable to those of the S-FEM; however, errors for ES-FEM and NS-FEM are globally small, around À0.4% and À0.5% respectively.

Near-incompressibility
In this section, near-incompressibility tests are studied.For these examples, different bulk moduli are used, j =10 2 ,1 0 3 and 10 4 .With those bulk moduli, for which the Poisson's ratio is close to 0.5, the model becomes nearly-incompressible.The geometry of the structure is illustrated in Fig. 12.Because an analytical solution is not available for this problem we calculate a reference solution numerically using a mixed finite element method on a highly-refined mesh within the DOLFIN finite element software [61,62].As shown in Fig. 13, edge-and node-based S-FEM are proven to be accurate and reliable for both compressible and nearly-incompressible problems.The x-and y-directions represent logarithmic number of global degrees of freedom and logarithm of a fraction of numerical results and analytical solution, respectively.When the Poisson's ratio is close to 0.5, the convergence of the ES-FEM becomes slow.The NS-FEM provides here an upper bound solution.Tables 5-7 provide the strain energy relative errors for FEM, ES-FEM and NS-FEM.As shown in Table 7, S-FEM handles near-incompressibility excellently, with results provided by NS-FEM up to 140 times more accurate than the FEM.

Mesh distortion sensitivity
In this section, a mesh distortion sensitivity is considered.For this test, results of DOLFIN finite element software are compared with the gradient smoothing techniques.We use artificially distorted meshes which are given by [35]: where r c is a random number between À1.0 and 1.0, a is the magnitude of the distortion and Mx, My are initial regular element sizes in the x-and y-direction.The higher a the more distorted the mesh.
The geometry of the examples is given in Sections 2.2.6 and 5.2.4 of [58] (see also Fig. 14).Consider a rectangle in the reference Cartesian coordinates (X, Y) defined by:   Strain energy relative error is given by: WNumerical ÀWExact WExact Â 100%.
The deformation in cylindrical coordinates is: For implementation, the given cylindrical coordinates are rewritten in Cartesian form: Dirichlet boundary conditions.Dirichlet boundary conditions are imposed as following: Bottom edge (Y = ÀB): Left-hand edge (X = A 1 ): Parameters, a = 0.9, A 1 =2,A 2 = 3 and B = 2 for Dirichlet boundary conditions, the distortion factors a = 0.1, 0.2, 0.3, 0.4 and 0.45 for mesh distortion, and l = 0.6 and j = 1.95 (E % 1.6326, v % 0.3605) for neo-Hookean material, are used in this test.In addition, we can obtain an exact solution for this example [58].The deformation gradient F for this problem is:   Strain energy relative error is given by: WNumerical ÀWExact WExact Â 100%.where f, g, f 0 and g 0 are: The strain energy density can be rewritten as: where I 1 = f 02 +(fg 0 ) 2 + 1.Hence, Eq. ( 59) is: where a = 0.9 and then strain energy is WðXÞdYdX % 4:485618.Fig. 15 illustrates the deformed configurations of bending block with different distortion factors.When the distortion factor a is close to 0.5, the meshes become severely distorted.In this test, we only impose Dirichlet boundary conditions which means that applied external forces vanish and no body force acts on the domain.
Detailed values of strain energy relative error are given in Tables 8-11.The relative error of S-FEM is much less than that of the FEM: errors for ES-FEM are about À1.0% and À1.9%, those of NS-FEM are around À1.5% and À3.5% with finer meshes (2 Â 32 and 4 Â 32) and highly distorted meshes (a = 0.45) whilst errors for FEM are approximately À0.7% and 260%.Moreover, MINI element gives accurate results; however, when meshes are severely distorted, MINI element fails to converge.This indicates that the S-FEM can effectively alleviate the mesh distortion sensitivity.

Edge-based smoothing strain using bubble functions
Lastly, we provide the results of the enhanced strain smoothing method, implementing Cook's membrane with the larger bulk moduli j =10 5 ,10 6 and 10 7 .Parameters which are used in this section are exactly the same as in the previous section.Fig. 16 illustrates the convergence of the strain energy.DOLFIN finite element software based on mixed finite element formulation on highly refined meshes is used as a reference solution.
The strain energy convergence of given techniques are described in Fig. 16.As shown in Fig. 16, NS-FEM performs much better than ES-FEM and the classical FEM.However the bubble- enhanced ES-FEM produces more accurate results and higher convergence rates than NS-FEM.It is clearly shown that the bubble function within ES-FEM effectively improves the quality of T3 elements in the nearly-incompressible limit.
Relative errors in the strain energy for FEM, ES-FEM, NS-FEM and ES-FEM with the bubbles are given in .The relative errors of FEM and ES-FEM are around 50% for both methods with fine meshes, whereas NS-FEM and bES-FEM prevent volumet-Table 8 Strain energies relative error for the bending of a rectangle using the standard FEM with a = 0.0, 0.1, 0.2, 0.3, 0.4 and 0.45.The higher the value of a the more distorted the mesh is.Strain energy relative error is given by: WNumerical ÀWExact WExact Â 100%.

Table 9
Strain energies relative error of bending of a rectangle for the ES-FEM with a = 0.0, 0. Strain energy relative error is given by: WNumerical ÀWExact WExact Â 100%.

Table 10
Strain energies relative error of bending of a rectangle for the NS-FEM with a = 0.0, 0.1, 0.2, 0.3, 0.4 and 0.45.The higher the value of a the more distorted the mesh is.Strain energy relative error is given by: WNumerical ÀWExact WExact Â 100%.

Table 11
Strain energies relative error of bending of a rectangle for the MINI element with a = 0.0, 0.1, 0.2, 0.3, 0.4 and 0.45.The higher the value of a the more distorted the mesh is.
MINI ric locking in quasi-incompressible limit ðm !0:5Þ.Notable improvement of the bubble-enhanced ES-FEM is that its relative errors, À0.8% for the bulk moduli j =10 5 ,1 0 6 and 10 7 with 8 Â 8 elements, are smaller than those of NS-FEM, 0.9% for the bulk mod-uli j =10 5 ,1 0 6 and 10 7 with 16 Â 16 elements.In other words, bubble-enriched ES-FEM has more accurate results and faster convergence and overcomes the overestimation of the stiffness matrix and the locking problems.Strain energy relative error is given by: WNumerical ÀWReference WReference Â 100%.

Table 14
Strain energies relative error of the Cook's membrane for the standard NS-FEM with the higher bulk moduli j =10 5 ,1 0 6 and 10

Conclusions
In this work, we reviewed the basic theory of the smoothed finite element method in linear and finite elasticity.Through numerical examples, we showed the accuracy and convergence of the proposed method in hyperelasticity, and its ability to overcome locking and mesh distortion effects.
We also presented the analytical solutions for Simple Shear deformation with Dirichlet boundary conditions, Uniform Extension with lateral contraction with both Dirichlet and mixed boundary conditions, and ''Not-So-Simple" Shear deformation with Dirichlet boundary conditions.We analysed the accuracy of the proposed technique, compared to those analytical solutions and numerical results obtained with FEM.
To show the ability of the method to handle nearlyincompressible problems, bulk moduli j = 1.95, 10, 10 2 ,1 0 3 and 10 4 were used.For nearly-incompressible problems, FEM provides very slow convergence, whereas S-FEM is shown to be stable and accurate.When the bulk modulus is large, ES-FEM reveals relatively slower convergence than NS-FEM.Although NS-FEM itself is stable and reliable for near-incompressibility, enhanced ES-FEM, using the bubble functions, sufficiently improves the quality of lower-order simplex element and prevents the locking issue under large deformations.Lastly, to study mesh distortion sensitivity, artificially distorted meshes are constructed with various distortion factors.For heavily distorted meshes, FEM shows unreliable results, whilst S-FEM performs very well.
As shown in the numerical examples the S-FEM is able to alleviate the spurious effects of both shear locking and mesh distortion whilst requiring only simplex elements, meshes of which are easily generated.It is therefore apparent that these elements, which are easily implemented within existing FE codes offer an alternative to quadrilateral elements.We are currently extending this work to 3D hyperelastic problems and proceeding to GPU implementation for real-time applications [63]. ðA:1Þ where the undetermined coefficients a ij and b i , for i, j = 1, 2, are constant.
We here consider the smoothed deformation gradient F for ES-FEM.The deformation gradient on a triangle MABC for the standard FEM in Fig. A1 is: For the smoothed deformation gradient F in the smoothing domain 2 Þð A:2Þ Substituting Eq. (A.2) into Eq.(A.1), the displacement field on midpoint O 1 is given by: Similarly, the displacement fields on node B and C can be written as: Hence, the displacements on the mid-point O 1 are given by: 2 À a 22 hÞ ðA:5Þ From Eq. (A.5), the undetermined coefficient a ij are defined as follows: Similarly, the undetermined coefficient a ij for triangle MDCB in Fig.
A.1 are given by: The smoothed deformation gradient is given by Hu et al. [64]: where U is: 0 otherwise ðA:6Þ and then:

Fig. 4 .
Fig. 4. Renewed basis functions and the cubic bubble function associated the centroid of a right 45°three-node triangular (T3) element.

Fig. 16 .
Fig.16.Strain energy convergence of the Cook's membrane with the bulk moduli j =10 5 ,1 0 6 and 10 7 : DOLFIN finite element software is to be the reference solution.

Table 1
Strain energy relative error (Â10 À12 %) for the simple shear deformation with Dirichlet boundary conditions: FEM, edge-based smoothing and node-based-smoothing.

Table 2
Strain energy relative error (Â10 À12 %) for the uniform extension with lateral contraction with Dirichlet boundary conditions: FEM, edge-based smoothing and node-based smoothing.

Table 3
Strain energy relative error (Â10 À12 %) for the uniform extension with lateral contraction with mixed Dirichlet and Neuman boundary conditions: FEM, edgebased smoothing and node-based smoothing.

Table 4
Strain energy relative error (%) for the ''Not-so-simple" shear example: edge-based and node-based smoothing.

Table 5
Strain energies relative error of the Cook's membrane for the standard FEM, ES-FEM and NS-FEM with bulk modulus j = 100.

Table 6
Strain energies relative error of the Cook's membrane for the standard FEM, ES-FEM and NS-FEM with bulk modulus j = 1000.

Table 7
Strain energies relative error of the Cook's membrane for the standard FEM, ES-FEM and NS-FEM with bulk modulus j = 10,000.
Fig. 14.The geometry of bending of a rectangle.C.-K. Lee et al. / Computers and Structures 182 (2017) 540-555 1, 0.2, 0.3, 0.4 and 0.45.The higher the value of a the more distorted the mesh is.

Table 12
Strain energies relative error of the Cook's membrane for the standard FEM with bulk moduli j =10 5 ,1 0 6 and 10 7 .

Table 13
Strain energies relative error of the Cook's membrane for the standard ES-FEM with bulk moduli j =10 5 ,1 0 6 and 10 7 . 7.

Table 15
Strain energies relative error of the Cook's membrane for the standard ES-FEM with the bubbles with bulk moduli j =10 5 ,1 0 6 and 10 7 .