Strategy to Improve Edge-Based Smoothed Finite Element Solutions Using Enriched 2D Solid Finite Elements

Abstract

In this paper, we present an automatic procedure that enhances the solution accuracy of edge-based smoothed 2D solid finite elements (three-node triangular and four-node quadrilateral elements). To obtain an enhanced solution, an adaptive enrichment scheme that uses enriched 2D solid finite elements and can effectively improve solution accuracy by applying cover functions adaptively without mesh-refinement is adopted in this procedure. First, the error of the edge-based finite element solution is estimated using a devised error estimation method, and appropriate cover functions are assigned for each node. While the edge-based smoothed finite elements provide piecewise constant strain fields, the proposed enrichment scheme uses the enriched finite elements to obtain a higher order strain field within the finite elements. Through various numerical examples, we demonstrate the accuracy improvement and efficiency achieved.

1. Introduction

The finite element method (FEM) is one of the most successful numerical methods and has been widely used to solve mechanics problems in various engineering fields [1,2,3,4,5]. Despite its great success, the standard finite element method has some limitations: the analysis accuracy depends on the mesh used and the stress prediction accuracy is considerably poor when 3-node triangular (2D) and 4-node tetrahedral (3D) elements are used [5,6,7]. Many studies have been conducted to overcome these limitations of the standard FEM (ST-FEM), and several numerical methods, such as extended FEM (EX-FEM), enriched FEM (EN-FEM), and smoothed FEM (S-FEM), have been developed [7,8,9,10,11].

The S-FEM is a relatively recently proposed method among the several techniques mentioned above and has been successfully applied to various engineering problem [12,13,14,15,16,17,18,19,20,21,22,23,24]. In S-FEM, the strain smoothing technique originally developed for the Galerkin mesh-free method is applied to the FEM, and constant strain fields are constructed in newly defined smoothing domains [12,13]. Depending on how the smoothing domain is defined, node-based, cell-based, edge-based, and face-based smoothed FEM have been proposed. Each S-FEM has different characteristics and performance depending on the smoothing domain, and it is generally known that edge-based S-FEM (ES-FEM) is the most effective [19,20,21,22]. The solution accuracy of the ES-FEM is more accurate compared to the ST-FEM, especially when the same lower order element mesh is used, and unlike other methods such as EX-FEM and EN-FEM, no additional degrees of freedoms (DOFs) are required [6,8,11,25]. However, when high gradient occurs in the smoothing domain, the predictive ability of the ES-FEM with a lower order element mesh is not good enough to be used in engineering practice [25,26].

The enrichment scheme was developed as another method of improving finite element solution accuracy [5,9,11]. With the enrichment scheme, the solution space of the standard finite element can be expanded by introducing cover functions corresponding to enriched DOFs. The main advantage of the enriched FEM (EN-FEM) is that the solution improvement in the sense of the p-version of the ST-FEM is available without introducing additional nodes. In addition, adaptively applying the cover function to a local area where solution accuracy needs to be improved can enhance solution accuracy efficiently [5,11]. However, depending on the mesh used when applying the polynomial cover function, a linear dependence (LD) problem that leads to ill-conditioned global matrices can occur [27,28]. Many methods have been suggested to suppress the LD problem of enriched finite elements [3,29,30,31]. Among the suggested methods, the LD problem was resolved in a simple and effective way by adopting piecewise linear shape functions and suppressing enriched DOFs [32,33]. In addition, a procedure for automatically improving the analysis accuracy of the standard finite elements using enriched finite elements has been proposed [34,35]. In this procedure, the use of cover functions is determined automatically based on the analysis results of the standard finite elements. The mesh-refinement is not required, but the number of DOFs used increases.

In this paper, we propose an automatic procedure to effectively improve the solution accuracy of edge-based 2D solid finite elements (4-node quadrilateral and 3-node triangular) without mesh-refinement or introducing additional nodes. An error indicator that estimates the error of the finite element solution and cover function selection scheme based on the error indicator are adopted, and the enriched finite elements are employed to incorporate cover functions in the analysis results. The adaptive enrichment scheme can enhance the solution of edge-based finite elements and capture high gradients efficiently.

In the following sections, we firstly review the edge-based finite elements method (ES-FEM) and the enriched finite element method (EN-FEM). We then suggest an automatic procedure that provides solution improvement through adaptive enrichment including the error indicator and cover function selection scheme. Through several examples, the feasibility of improving the accuracy of edge-based 2D finite element solutions is demonstrated.

2. Enrichment Scheme for the Edge-Based Smoothed 2D Solid Finite Elements

In this section, we briefly review the edge-based smoothed three-node triangular and four-node quadrilateral elements (ES-QD4, ES-T3) and the enriched three-node triangular and four-node quadrilateral elements (EN-QD4, EN-T3), and suggest a strategy to obtain accurate solutions with an enrichment scheme.

It is not necessary to explicitly define the shape function to calculate the stiffness matrix of edge-based smoothed finite elements [6], but we present the formulation of edge-based smoothed finite elements (ES-FE) and enriched finite elements (EN-FE) with linear and piecewise linear shape functions [33]. These shape functions are adequate for ES-FE and EN-FE and also prevent the LD problem that can occur due to the enrichment scheme.

2.1. Edge-Based Smoothed 2D Finite Elements

The geometry of the ES-FE is given by

$$x={\displaystyle \sum _{i=1}^n{h}_i^n(r,s)x_i}\mathrm{with}{x}_i=[\begin{array}{cc}{x}_i& y_i{]}^T\end{array},$$

(1)

where $x_i$ is the position vector corresponding to node $i$ on the global Cartesian coordinate system $(x,y)$ shown in Figure 1a, $n$ is the number of nodes of element, and $h_i^n(r,s)$ are the 2D linear interpolation functions corresponding to node $i$ defined in the natural coordinate system $(r,s)$ in Figure 1b. The four-node quadrilateral and three-node triangular elements have center points ($r=s=0$ for the four-node element and $r=s=1/3$ for the three-node element) and are divided into four or three sub-triangles (cells), respectively, using its nodes and center points.

For the four-node quadrilateral element ($n=4$), the piecewise linear shape functions defined in each triangular cell shown in Figure 1b are [32]

$$\begin{gathered} h_1^4=(1+2r+s)/4,h_2^4=(1-2r+s)/4,h_3^4=(1-s)/4,h_4^4=(1-s)/4\mathrm{in}\mathrm{C}1, \\ {h}_1^4=(1+r)/4,h_2^4=(1-r+2s)/4,h_3^4=(1-r-2s)/4,h_4^4=(1+r)/4\mathrm{in}\mathrm{C}2, \\ {h}_1^4=(1+s)/4,h_2^4=(1+s)/4,h_3^4=(1-2r-s)/4,h_4^4=(1+2r-s)/4\mathrm{in}\mathrm{C}3, \\ {h}_1^4=(1+r+2s)/4,h_2^4=(1-r)/4,h_3^4=(1-r)/4,h_4^4=(1+r-2s)/4\mathrm{in}\mathrm{C}4. \end{gathered}$$

(2)

The linear shape functions for the three-node triangular element ($n=3$) are

$$h_1^3(r,s)=1+r+s,h_2^3(r,s)=r,h_3^3(r,s)=s,$$

(3)

where the linear shape functions are the same for all cells shown in Figure 1b.

The displacement of the ES-FE is interpolated by

$$\overline{u}={\displaystyle \sum _{i=1}^n{h}_i^n(r,s){\overline{u}}_i}\mathrm{with}{\overline{u}}_i=[\begin{array}{cc}{u}_i& v_i{]}^T\end{array},$$

(4)

in which ${\overline{u}}_i$ is the displacement vector of node $i$ in the global Cartesian coordinate system.

Employing the standard isoparametric finite element procedure with Equations (1) to (4), the strain field within a cell c of an element e is obtained using

$${\epsilon }^{(e,c)}=B^{(e,c)}{\overline{u}}^{(e)}\mathrm{with}{B}^{(e,c)}=[\begin{array}{cccc}{B}_1^c& B_2^c& \cdots & B_n^c]\end{array},{\overline{u}}^{(e)}=[{\begin{array}{cccc}{\overline{u}}_1& {\overline{u}}_2& \cdots & {\overline{u}}_n]\end{array}}^T,$$

(5)

where ${\epsilon }^{(e,c)}=\begin{array}{ccc}[{\epsilon }_{xx}& {\epsilon }_{yy}& 2{\epsilon }_{xy}{]}^T\end{array}$, $B^{(e,c)}$ is the displacement-strain relation matrix on the cell c in the element e, ${\overline{u}}^{(e)}$ is the nodal displacement vector of the element, and $B_i^c$ is the strain-displacement relation matrix corresponding to node i in cell c, as seen in Figure 2a. Since the geometry and displacement interpolations are based on the linear shape functions in Equations (2) and (3), the strain components in Equation (5) are constants in each cell.

For the ES-FEM, the smoothing domain is defined at each edge with its neighboring cells. Let us consider three-node triangular and four-node quadrilateral elements attached to the edge, colored in red (see, Figure 2). The edge-based smoothing domain is formed as the union of cells attached to the edge and belonging to different elements, and the smoothed strain corresponding to this smoothing domain is given by

$$\overline{\epsilon }=\frac{1}{A_{}^{(1,4)}+A_{}^{(2,4)}}(A_{}^{(1,4)}{\epsilon }^{(1,4)}+A_{}^{(2,3)}{\epsilon }^{(2,3)}),$$

(6)

in which $A_{}^{(1,4)}$ and $A_{}^{(2,3)}$ are the areas of cells attached to the target edge, and ${\epsilon }^{(1,4)}$ and ${\epsilon }^{(2,3)}$ are the strains of the neighboring cells. The strain field of the ES-FE is also constant in the smoothing domain, like the strain vector in Equation (5).

The edge-based smoothed strain within cell c of the target element e can be expressed in matrix and vector form as

$${\overline{\epsilon }}^{(e,c)}={\overline{B}}^{(e,c)}{\overline{u}}^{(e,c)}\mathrm{with}{\overline{B}}^{(e,c)}=[\begin{array}{cccc}{\overline{B}}_1& {\overline{B}}_2& \cdots & {\overline{B}}_{n_s}]\end{array},{\overline{u}}^{(e,c)}=[{\begin{array}{cccc}{\overline{u}}_1& {\overline{u}}_2& \cdots & {\overline{u}}_{n_s}]\end{array}}^T,$$

(7)

where $n_s$ is the number of nodes belonging to elements including cells constituting the smoothing domain (see Figure 2c). The edge-based smoothed elements also have a constant strain field and pass all basic tests (patch, isotropy, and zero energy mode tests).

2.2. Enriched 2D Finite Elements

Here we review the formulation of the EN-FEM. The enrichment scheme can expand the solution space of the FEM without mesh refinement or additional nodes, and the solution accuracy can be improved effectively by considering the adequate cover function for a particular problem.

The geometry interpolation of the EN-FE is the same as the interpolation of the standard and smoothed finite elements and is given in Equation (1). The interpolation of displacement is given by multiplying the shape functions with cover functions defined over cover area as follows:

$$\tilde{u}(r,s)={\displaystyle \sum _{i=1}^n{h}_i^n(r,s){\tilde{u}}_i}\mathrm{with}{\tilde{u}}_i={[\begin{array}{cc}{\tilde{u}}_i& {\tilde{v}}_i\end{array}]}^T,$$

(8)

where ${\tilde{u}}_i$ and ${\tilde{v}}_i$ are cover functions for the displacements u and v, respectively. The cover area consists of elements attached to node i as shown in Figure 3.

The polynomial cover functions are

$${\tilde{u}}_i=p_i(x)u_i^u,{\tilde{v}}_i=p_i(x)u_i^v\mathrm{in}{C}_i$$

(9)

with

$$p_i(x)=\left[\begin{array}{cccccc}1& {\xi }_i& {\eta }_i& {\xi }_i^2& \cdots & {\eta }_i^d\end{array}\right],{\xi }_i=\frac{(x-x_i)}{{\chi }_i},{\eta }_i=\frac{(y-y_i)}{{\chi }_i},$$

(10)

$$u_i^u={\left[\begin{array}{cccccc}{\overline{u}}_i& {\stackrel{⌢}{u}}_i^{\xi }& {\stackrel{⌢}{u}}_i^{\eta }& {\stackrel{⌢}{u}}_i^{{\xi }^2}& \cdots & {\stackrel{⌢}{u}}_i^{{\eta }^d}\end{array}\right]}^T,u_i^v={\left[\begin{array}{cccccc}{\overline{v}}_i& {\stackrel{⌢}{v}}_i^{\xi }& {\stackrel{⌢}{v}}_i^{\eta }& {\stackrel{⌢}{v}}_i^{{\xi }^2}& \cdots & {\stackrel{⌢}{v}}_i^{{\eta }^d}\end{array}\right]}^T,$$

(11)

where $p_i$ is the polynomial basis vector for node i, d is the degree of the polynomial bases, ${\chi }_i$ is the largest edge length of elements belonging to $C_i$, and $u_i^u$ and $u_i^v$ are the degrees of freedom vectors corresponding to the bases of the polynomial cover functions for the displacements in the x- and y- directions, respectively.

Substituting cover functions in Equation (9) into Equation (8), the displacement interpolation is represented by

$$\begin{array}{l}\tilde{u}=\overline{u}+\stackrel{⌢}{u}\\ ={\displaystyle \sum _{i=1}^n{h}_i^n(r,s){\overline{u}}_i}+{\displaystyle \sum _{i=1}^n{\stackrel{⌢}{H}}_i^n(r,s){\stackrel{⌢}{u}}_i}\end{array}$$

(12)

with

$${\stackrel{⌢}{u}}_i=\left[\begin{array}{c}{\stackrel{⌢}{u}}_i^u\\ {\stackrel{⌢}{u}}_i^v\end{array}\right],{\stackrel{⌢}{H}}_i^n(r,s)=\left[\begin{array}{cc}{\stackrel{⌢}{h}}_i^n(r,s)& 0\\ 0& {\stackrel{⌢}{h}}_i^n(r,s)\end{array}\right],$$

(13)

where ${\stackrel{⌢}{u}}_i$ and ${\stackrel{⌢}{H}}_i^n$ are the enriched DOFs vector and the corresponding interpolation matrix, respectively.

For the linear cover function used ($d=1$), the following components of the interpolation matrix and the enriched DOFs vector are employed

$${\stackrel{⌢}{h}}_i^n(r,s)=h_i^n(r,s)[\begin{array}{cc}{\xi }_i^{}& {\eta }_i]\end{array},{\stackrel{⌢}{u}}_i^u={\left[\begin{array}{cc}{\stackrel{⌢}{u}}_i^{\xi }& {\stackrel{⌢}{u}}_i^{\eta }\end{array}\right]}^T,{\stackrel{⌢}{u}}_i^v={\left[\begin{array}{cc}{\stackrel{⌢}{v}}_i^{\xi }& {\stackrel{⌢}{v}}_i^{\eta }\end{array}\right]}^T.$$

(14)

When the quadratic cover function is used ($d=2$), the components of the interpolation matrix and the enriched DOFs vector are

$$\begin{gathered} {\stackrel{⌢}{h}}_i^n(r,s)=h_i^n(r,s)\left[\begin{array}{ccccc}{\xi }_i& {\eta }_i& {\xi }_i^2& {\xi }_i{\eta }_i& {\eta }_i^2\end{array}\right], \\ {\stackrel{⌢}{u}}_i^u={\left[\begin{array}{ccccc}{\stackrel{⌢}{u}}_i^{\xi }& {\stackrel{⌢}{u}}_i^{\eta }& {\stackrel{⌢}{u}}_i^{{\xi }^2}& {\stackrel{⌢}{u}}_i^{\xi \eta }& {\stackrel{⌢}{u}}_i^{{\eta }^2}\end{array}\right]}^T,{\stackrel{⌢}{u}}_i^v={\left[\begin{array}{ccccc}{\stackrel{⌢}{v}}_i^{\xi }& {\stackrel{⌢}{v}}_i^{\eta }& {\stackrel{⌢}{v}}_i^{{\xi }^2}& {\stackrel{⌢}{v}}_i^{\xi \eta }& {\stackrel{⌢}{v}}_i^{{\eta }^2}\end{array}\right]}^T. \end{gathered}$$

(15)

For the cubic cover function used ($d=3$), the components of the interpolation matrix and the enriched DOFs vector become

$$\begin{gathered} {\stackrel{⌢}{h}}_i^n(r,s)=h_i^n(r,s)\left[\begin{array}{ccccc}{\xi }_i& {\eta }_i& {\xi }_i^2& \cdots & {\eta }_i^3\end{array}\right], \\ {\stackrel{⌢}{u}}_i^u={\left[\begin{array}{ccccc}{\stackrel{⌢}{u}}_i^{\xi }& {\stackrel{⌢}{u}}_i^{\eta }& {\stackrel{⌢}{u}}_i^{{\xi }^2}& \cdots & {\stackrel{⌢}{u}}_i^{{\eta }^3}\end{array}\right]}^T,{\stackrel{⌢}{u}}_i^v={\left[\begin{array}{ccccc}{\stackrel{⌢}{v}}_i^{\xi }& {\stackrel{⌢}{v}}_i^{\eta }& {\stackrel{⌢}{v}}_i^{{\xi }^2}& \cdots & {\stackrel{⌢}{v}}_i^{{\eta }^3}\end{array}\right]}^T. \end{gathered}$$

(16)

Using the displacement–strain relation, the strain vector within cell c of the target element e is given by

$$\begin{array}{cc}{\widehat{\epsilon }}^{(e,c)}& ={\widehat{B}}^{(e,c)}{\widehat{u}}^{(e)}\hfill \\ & =B^{(e,c)}{\overline{u}}^{(e)}+{\stackrel{⌢}{B}}^{(e,c)}{\stackrel{⌢}{u}}^{(e)}\hfill \\ & ={\epsilon }^{(e,c)}+{\stackrel{⌢}{\epsilon }}^{(e,c)}\hfill \end{array}$$

(17)

with

$${\stackrel{⌢}{B}}^{(e,c)}=[\begin{array}{cccc}{\stackrel{⌢}{B}}_1^c& {\stackrel{⌢}{B}}_2^c& \cdots & {\stackrel{⌢}{B}}_n^c]\end{array},{\stackrel{⌢}{u}}^{(e)}=[{\begin{array}{cccc}{\stackrel{⌢}{u}}_1& {\stackrel{⌢}{u}}_2& \cdots & {\stackrel{⌢}{u}}_n]\end{array}}^T,$$

(18)

in which ${\stackrel{⌢}{u}}^{(e)}$ is the enriched DOFs vector of the element e, and ${\stackrel{⌢}{B}}_i^c$ is the strain–displacement relation matrix corresponding to the enriched DOFs at node i in cell c of the target element e.

Note that the displacement of the ES-FE in Equation (12) is divided into two parts: Linear displacement $\overline{u}$ and higheorder displacement $\stackrel{⌢}{u}$. $\overline{u}$ is the standard finite element interpolation and $\stackrel{⌢}{u}$ is the additional enriched higher order displacement interpolation.

This polynomial enrichment scheme improves the solution accuracy in the sense of the p-version of the ST-FEM without adding extra nodes, and the cover functions with different degrees of polynomial can be used for each node. Thus, by applying cover functions to a local area where solutions need to be improved can increase the finite element solution accuracy effectively without mesh refinement [32,33].

2.3. The Procedure for Improving Edge-Based Finite Element Solutions

There have been studies to improve the accuracy of the standard finite element solutions using the enrichment scheme. Kim and Lee presented an error indicator for each node to estimate the standard finite element solution error and proposed a scheme to apply polynomial cover functions to the appropriate node using the calculated error indicator [34]. This procedure can automatically improve the stress and displacement predictions and is effective in terms of computational cost since no mesh-refinement or extra nodes are required.

Here, we present an automatic procedure that effectively improves the solution accuracy of edge-based smoothed 2D solid finite elements (ES-FE) with enriched 2D solid finite elements (EN-FE). Similar to the previous proposed procedures, the improvement of the solution accuracy is composed of two steps. First, the ES-FE solution is calculated, and then the error for each node is estimated through the error indicator. Based on this error indicator, appropriate cover functions are selected for each node, and re-analysis is performed to obtain an improved analysis result.

The error indicator that estimates errors at each node is defined by

$$M_i=\frac{1}{2}\left\{\frac{J_i}{\overline{J}}+\left(\frac{\overline{J}}{\overline{\tau }}\right)\frac{{\tau }_i}{\overline{\tau }}\right\}\left(\frac{{\chi }_i}{L_c}\right)$$

(19)

with

$$\overline{J}=\frac{1}{N}{\displaystyle \sum _{i=1}^N{J}_i},\overline{\tau }=\frac{1}{N}{\displaystyle \sum _{i=1}^N{\overline{\tau }}_i},J_i=\max ({\tau }_i)-\min ({\tau }_i),$$

(20)

where ${\tau }_i$ is the stress at node i, $J_i$ is the largest stress jump at node i, $N$ is the number of nodes, and ${\chi }_i$ is the largest edge length of elements attached to node i [35]. The error indicator in Equation (19) does not include any artificial constants and any stress quantity of interest can be used for this error indicator. In this paper, the von Mises stress quantity is applied to all examples.

Based on the error indicator in Equation (19), the cover selection scheme determines the nodes to which the cover function is applied and the degree of the cover function. This scheme should include minimum parameters and present the optimal nodes and degrees of the cover function to improve the solution accuracy.

In the cover selection scheme, index ($I_i$), which presents the arrangement of the nodal error indicator ($M_i$ with $i=1,2,\cdots ,N$) in ascending order, is defined. Then, the ascending index ($I_i$) and normalized index (${\widehat{I}}_i$) are given as follows:

$$I_i,I_j\in \{1,2,\cdots ,N\},M_{I_i}<M_{I_j}\mathrm{when}{I}_i<I_j$$

(21)

and

$${\widehat{I}}_i=I_i/N.$$

(22)

By using the normalized index, a weighted average of the error indicator is given by

$${\overline{M}}_w=\left({\displaystyle {\sum }_{i=1}^N{M}_i\cdot {\widehat{I}}_i}\right)/N.$$

(23)

With the weighted average of the error indicator in Equation (23) and artificial parameters (${\alpha }_i$), the cover selection is given as follows:

$$d(i)=\{\begin{array}{cc}0& if{\widehat{I}}_i<(1-{\alpha }_1{\overline{M}}_w)\hfill \\ 1& if(1-{\alpha }_1{\overline{M}}_w)\le {\widehat{I}}_i<(1-{\alpha }_2{\overline{M}}_w)\hfill \\ 2& if(1-{\alpha }_2{\overline{M}}_w)\le {\widehat{I}}_i<(1-{\alpha }_3{\overline{M}}_w)\hfill \\ 3& if(1-{\alpha }_3{\overline{M}}_w)\le {\widehat{I}}_i\hfill \end{array}.$$

(24)

For nodes with larger errors, a higher degree of cover function is applied, and artificial parameters (${\alpha }_i$) are the only values that should be determined. Kim and Lee suggested using $({\alpha }_1,{\alpha }_2,{\alpha }_3)=(21,14,7)$ as default values and demonstrated several examples that showed very effective improvement of the solution accuracy. These values for parameters (${\alpha }_i$) are adopted for all examples in the following section.

With the cover function selection scheme in Equation (24), the finite element model is updated to obtain a better solution and capture high gradients. The stiffness matrix of the updated FE model with smoothing and enrichment schemes can be calculated. As shown in Figure 4, smoothing and enrichment schemes are applied separately, so that they effectively provide low and high gradient solutions, respectively.

The presented procedure improves the solution accuracy of edge-based smoothed finite elements, which are known to have better performance than standard finite elements, while the automatic procedure proposed by Lee and Kim [35] improves the solution for standard finite elements. Iterative analysis is also possible by adopting a higher degree cover function or by applying a cover function to a node to which the cover function is not applied. Repeating this re-analysis requires parameter adjustment in Equation (24) at each iteration. The most accurate solution can be obtained when applying the highest order cover function to all nodes in a given finite element mesh. In this paper, we present a procedure that can effectively obtain an accurate solution through single re-analysis automatically without parameter adjustment.

We performed three basic tests (the isotropy, patch, and zero energy mode tests) with arbitrary cover function application to each case, as shown in Figure 4, and all cases passed these basic tests. With this cover selection scheme, the solution accuracy of the ES-FE can be improved effectively without mesh-refinement or introducing additional nodes and its performance is demonstrated through examples in the next section.

3. Numerical Examples

In this section, we demonstrate the feasibility of improving 2D ES-FE solutions with an adaptive enrichment scheme presented in the previous section. This procedure yields solutions with desired accuracy by applying appropriate cover functions for each node, and the higher-order displacement field can be obtained effectively. The Cook beam problem, Ad-hoc problem, Wheel problem, and L-shape structure problems are considered, and the finite element models are constructed using meshes of linear elements (four-node quadrilateral and three-node triangular elements) and quadratic elements (nine-node quadrilateral and six-node triangular elements).

Displacement, strain, von Mises stress, and energy norms were evaluated, and the relative errors of displacement ($v$) in the y-direction and von Mises stress ($\tau$) are given by

$$e_v=\frac{\left|v_h-v_{ref}\right|}{v_{ref}},$$

(25)

$$e_{\tau }=\frac{\left|{\tau }_h-{\tau }_{ref}\right|}{{\tau }_{ref}},$$

(26)

where the subscripts ‘h’ and ‘ref’ represent the finite element and the reference solutions, respectively. The energy norm E and its relative error are denoted by

$$e_E=\frac{\left|E_h-E_{ref}\right|}{E_{ref}}\mathrm{with}E={\displaystyle {\int }_{\Omega }{\epsilon }^T\sigma dV}.$$

(27)

The reference solutions are obtained using fine meshes of quadratic elements.

3.1. Cook’s Skew Beam Problem

The first numerical example is the Cook beam problem shown in Figure 5 [4]. A distributed force is applied on the right edge, and the left edge is clamped. The plane stress condition is employed with Young’s modulus $Y=3\times {10}^7{\mathrm{N}/\mathrm{m}}^2$, Poisson’s ratio $p=0.3$, and density $\rho =1{\mathrm{kg}/\mathrm{m}}^3$. The two quadrilateral meshes are shown in Figure 5, where edge-based smoothed four-node quadrilateral element (ES-Q4) are used, and the proposed scheme is applied for each mesh to obtain solutions with improved accuracy.

The relative errors of the vertical displacement in Equation (25) and the energy norm in Equation (27) are given in

Table 1. Relative errors (%) in the energy norm ($E$) and the vertical displacement ($v$) at the point A in the Cook’s skew beam problem.

Mesh		${\mathit{e}}_{\mathit{E}}[\%]$		${\mathit{e}}_{\mathit{v}}[\%]$		DOFs
		Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers
ES-Q4	I	12.277	0.244	13.348	0.889	40	360
	II	0.183	0.112	0.653	0.116	2112	3986
Reference solutions: $E_{ref}^{}=3.997\times {10}^{-7}\mathrm{N}\cdot \mathrm{m}$, $v_{ref}^{}=8.367\times {10}^{-7}\mathrm{m}$ (DOFs = 18,624)

. Figure 6 provides the distributions of the strain component ${\epsilon }_{xy}$ when the edge-based four-node quadrilateral elements (ES-Q4) and the proposed scheme are applied to the meshes shown in Figure 6. The distributions of the von Mises stress along line BC is given in Figure 7. The reference solution is calculated with a fine mesh of nine-node quadrilateral elements. The actual computational time for constructing the stiffness matrix and solving the linear system of equations are summarized in

Table 2. Solution times (in seconds) for constructing the stiffness matrix and solving linear equations for the Cook’s skew beam problem.

Mesh			DOFs	Solution Times [s]
				Stiffness Construction	Equation Solver	Total
ES-Q4	I	Linear elements	40	0.006	0.000	0.006
		Adaptive covers	360	0.372	0.001	0.373
	II	Linear elements	2112	1.102	0.001	1.103
		Adaptive covers	3986	10.202	0.004	10.206

.

It is shown that the calculated results of the strain energy norm, displacement, strain, and von Mises stress with the adaptive enrichment scheme are more accurate than those without cover functions. Furthermore, both the case of the adaptive cover in Mesh 1 and the case of using the linear element in Mesh 2 show sufficiently accurate results; however, the adaptive method consumes approximately three times less computational time.

Transient response at point A is also calculated, and Newmark time integration method with $\Delta t=0.001\mathrm{s}$ is used. The initial condition and harmonic excitation at the right edge are given by

$${}^{0}U{=}^0\dot{U}{=}^0\ddot{U}=0,$$

(28)

$$F(t)=\sin (4\pi t)\mathrm{N}.$$

(29)

ES-Q4 and the adaptive enrichment scheme with Mesh 1 are used for dynamic analysis. The adaptive use of cover functions is the same with static analysis shown in Figure 6, and Figure 8 shows the transient response at Point A. The solution of the proposed scheme almost converges to the reference solution, and it is confirmed that the suggested scheme can improve not only the static but also dynamic solution accuracy effectively.

3.2. Ad-Hoc Problem

Here, the ad-hoc problem shown in Figure 9 is considered. For the given in-plane displacements

$$u={(1-x^2)}^2{(1-y^2)}^2{e}^{my}\cos mx,$$

(30)

$$v={(1-x^2)}^2{(1-y^2)}^2{e}^{my}\sin mx\mathrm{with}m=5,$$

(31)

we can obtain the body forces satisfying the equilibrium equations

$$f_x^B=-\left(\frac{\partial {\tau }_{xx}}{\partial x}+\frac{\partial {\tau }_{xy}}{\partial y}\right),f_y^B=-\left(\frac{\partial {\tau }_{yy}}{\partial y}+\frac{\partial {\tau }_{xy}}{\partial x}\right).$$

(32)

The body forces are applied to the whole problem domain, and the fixed boundary condition is imposed on the bottom line ($y=-1$).

The plane stress condition is employed with $Y=2\times {10}^5{\mathrm{N}/\mathrm{m}}^2$ and $p=0.3$, and the solution is obtained using the three triangular meshes shown in Figure 9 with the edge-based three-node triangular elements (ES-T3) and the suggested adaptive enrichment scheme. The reference solutions are obtained using a fine mesh of the standard nine-node quadrilateral elements.

Figure 10 provides the von Mises stress distributions for the three triangular meshes shown in Figure 9 and the degrees of cover functions applied based on the adaptive enrichment scheme. The relative errors in the energy norm in Equation (27) and maximum von Mises stress that are averaged stress values at each node in Equation (26) are summarized in

Table 3. Relative errors (%) in the energy norm ($E$) and the maximum von Mises stress ($\tau$) in the ad hoc problem.

Mesh		${\mathit{e}}_{\mathit{E}}[\%]$		${\mathit{e}}_{\mathit{\tau }}[\%]$		DOFs
		Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers
ES-T3	I	2.851	0.036	23.549	8.169	144	1224
	II	7.145	0.002	10.275	0.081	544	2996
	III	2.174	0.045	4.524	0.157	2178	6140
Reference solutions: $E_{ref}^{}=7.513\times {10}^7\mathrm{N}\cdot \mathrm{m}$, ${\tau }_{ref}^{}=8.362\times {10}^6{\mathrm{N}/\mathrm{m}}^2$ (DOFs = 74,112)

. The von Mises stress fields obtained by applying the adaptive enrichment scheme almost converge to the reference solution for all the meshes (Mesh I, II, and III) considered.

3.3. Wheel Problem

We here consider the problem of an automotive wheel with a radius of 0.2 (see Figure 11). The lower part of the wheel is subjected to a constant pressure load and the fixed boundary condition is imposed on the inner circle. The plane stress condition with $Y=2\times {10}^5{\mathrm{N}/\mathrm{mm}}^2$ and $p=0.3$, and two different meshes (Mesh I and II) shown in Figure 11 are used.

Numerical results obtained employing the edge-based smoothed finite elements (ES-FE) are compared with the results of the adaptive enrichment scheme.

Table 4. Relative errors (%) in the energy norm ($E$) and the maximum von Mises stress ($\tau$) at point A in the wheel problem.

Mesh		${\mathit{e}}_{\mathit{E}}[\%]$		${\mathit{e}}_{\mathit{\tau }}[\%]$		DOFs
		Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers
ES-Q4&T3	I	6.522	2.066	36.424	1.5739	1532	4144
	II	1.690	1.289	13.230	0.429	5268	8614
Reference solutions: $E_{ref}^{}=3.997\times {10}^{-7}\mathrm{N}\cdot \mathrm{mm}$, ${\tau }_{ref}^{}=8.367\times {10}^{-7}{\mathrm{N}/\mathrm{mm}}^2$ (DOFs = 74,414)

gives the relative errors in the energy norm in Equation (27) and maximum von Mises stress that are averaged stress values at each node in Equation (26). The von Mises stress distributions are given in Figure 12 and the reference solutions are obtained employing a mesh of 9058 nine-node quadrilateral elements.

The proposed scheme applies the cover functions to the area where the solution accuracy needs to be improved and there is high stress gradient. Thus, the increase of the number of DOFs could be minimized, and the accuracy of the strain energy and the maximum von Mises stress prediction is improved.

3.4. L-Shape Structure Problem

The L-shape structure is subjected to pressure load varying along the edges, as shown in Figure 13 and the boundary conditions are imposed as: $u=0$ along AD and $v=0$ along CD. The plane stress condition is used with material properties $Y=200\times {10}^9{\mathrm{N}/\mathrm{m}}^2$ and $p=0.3$.

The finite element models are constructed using the standard four-node and three-node elements (ST-Q4, ST-T3) and the edge-based smoothed four-node and four-node elements (ES-Q4, ES-T3) for three different meshes (Mesh I, II, III), as shown in Figure 13. Triangular element meshes are generated from the quadrilateral mesh by dividing each quadrilateral element into two triangular elements. To obtain an improved solution, the enrichment schemes for the ST-FE and ES-FE are employed, and accuracy improvement and computational efficiency in terms of the number of DOFs used are compared. The reference solutions are calculated using a mesh of 3072 nine-node quadrilateral elements.

The von Mises stress fields and the strain component ${\epsilon }_{xx}$ obtained from the linear elements and adaptive enrichment scheme are presented in Figure 14, Figure 15 and Figure 16, and the relative errors of strain energy and maximum von Mises stress are summarized in

Table 5. Relative errors (%) in the energy norm ($E$) and the maximum von Mises stress ($\tau$) in the L-shape structure problem when the meshes of quadrilateral elements (ST-Q4, ES-Q4) are used.

Mesh		${\mathit{e}}_{\mathit{E}}[\%]$		${\mathit{e}}_{\mathit{\tau }}[\%]$		DOFs
		Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers
ST-Q4	I	6.285	4.954	34.647	5.711	120	704
	II	1.800	1.131	21.463	1.397	432	1442
	II	0.486	0.245	13.449	1.256	1632	3410
EN-Q4	I	0.787	1.403	39.683	4.337	120	660
	II	0.507	0.810	27.257	1.595	432	1324
	III	0.235	0.196	16.740	0.912	1632	3164
Reference solutions: $E_{ref}^{}=5.120\times {10}^{-9}\mathrm{N}\cdot \mathrm{m}$, ${\tau }_{ref}^{}=1.751\times {10}^0{\mathrm{N}/\mathrm{m}}^2$(DOFs = 6336)

and

Table 6. Relative errors (%) in the energy norm ($E$) and the maximum von Mises stress ($\tau$) in the L-shape structure problem when the meshes of triangular elements (ST-T3, ES-T3) are used.

Mesh		${\mathit{e}}_{\mathit{E}}[\%]$		${\mathit{e}}_{\mathit{\tau }}[\%]$		DOFs
		Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers	Linear Elements	Adaptive Covers
ST-T3	I	21.206	0.120	49.037	4.598	120	890
	II	6.980	1.483	33.617	1.075	432	1946
	II	1.977	0.884	18.077	0.174	1632	4050
EN-T3	I	3.846	0.123	38.008	4.541	120	836
	II	0.733	0.976	27.116	1.660	432	1716
	III	0.139	0.208	17.314	0.185	1632	3826
Reference solutions: $E_{ref}^{}=5.120\times {10}^{-9}\mathrm{N}\cdot \mathrm{m}$, ${\tau }_{ref}^{}=1.751\times {10}^0{\mathrm{N}/\mathrm{m}}^2$(DOFs = 6336)

.

As for the quadrilateral mesh, the ES-Q4 provides greatly improved strain energy prediction compared to the ST-Q4, whereas the ES-Q4 appears to be less accurate in predicting the maximum von Mises stress. By applying the suggested procedure, the solution of the ES-Q4 is significantly improved, especially in stress prediction, and the increase of DOFs used is less than the ST-Q4, as seen in Table 5. For the triangular mesh, the ES-T3 provides more accurate results than the ST-T3 for the same mesh. Therefore, when the same procedure is applied, the smaller number of DOFs is used for the ES-T3, but sufficiently accurate results are provided (see Table 6).

4. Conclusions

The objective of this paper was to propose an adaptive enrichment scheme for the edge-based smoothed 2D solid finite elements to effectively obtain a solution with desired accuracy. When the edge-based smoothed finite elements provide a solution that is not sufficiently converged, the suggested error indicator estimates the error for each node. Then, the cover function selection scheme assigns the polynomial cover functions of the appropriate degrees (up to a cubic degree) for each node. In order to apply the cover function, enriched finite elements free from the linear dependence problem were adopted, and the strain smoothing scheme was applied only to the area that is not subject to the applied cover function. Through this procedure, analysis accuracy can be automatically improved without mesh-refinement or introducing additional nodes, and a high-gradient solution can be effectively captured.

The performance of the proposed procedure is demonstrated through several numerical examples. In each example, the procedure selected the appropriate nodes located in an area where solution accuracy needs to be ameliorated, and the strain energy, stress, strain, and displacement predictions were enhanced effectively. In future work, it would be valuable to extend an adaptive enrichment scheme for nonlinear analysis to enhance the convergence performance in iterative analysis. In addition, the adaptive enrichment scheme can be easily applied to the 3D solid, shell, and beam finite elements [33,36,37,38,39,40].