*2.2. Adaptive Mesh Refinement*

To minimize expected errors after a finite element solution has been achieved, an adaptive mesh refinement technique is used. The method of adaptive mesh refinement measures the mesh's adequacy and refines the mesh wherever the estimated error is large. Until user-definable error tolerance is reached, the system iterates the mesh refinement and solution. Because the precision of the solution depends on these tolerance limits, it is important for the use of adaptive mesh generators to provide a good understanding of the FEM in an effective manner. The method is referred to as adaptive since at all times the process relies on previous results. The adaptive remeshing method was carried out on the basis of the posteriori stress error standard scheme to achieve the optimum mesh from [16]. The software adopted a frontal solver that is an effective direct solver used to solve a linear equation system. In h-type adaptive mesh refinement, the major point is to obtain the ratio of element normal stress error to the average normal stress error for the entire domain, which is also known as the relative stress norm error, and a new size can be predicted from this ratio for the refinement method. The mesh size is defined in the procedure of each element as:

$$h\_{\mathfrak{E}} = \sqrt{2A\_{\mathfrak{E}}} \tag{4}$$

where *Ae* is the area of the triangle element. The norm stress error for each element is defined by

$$\begin{aligned} \left\|\boldsymbol{\varepsilon}\right\|\_{\epsilon}^{2} &= \int\_{\Omega^{\boldsymbol{r}}} \left(\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\*}\right)^{\top} (\boldsymbol{\sigma} - \boldsymbol{\sigma}^{\*}) d\Omega \\ &= \int\_{\Omega^{\boldsymbol{r}}} \left(\begin{Bmatrix} \boldsymbol{\sigma}\_{\boldsymbol{x}} \\ \boldsymbol{\sigma}\_{\boldsymbol{y}} \\ \boldsymbol{\tau}\_{\boldsymbol{xy}} \\ \boldsymbol{\sigma}\_{\boldsymbol{z}} \end{Bmatrix} - \left\{\begin{array}{cc} \boldsymbol{\sigma}^{\*}\_{\boldsymbol{x}} \\ \boldsymbol{\sigma}^{\*}\_{\boldsymbol{y}} \\ \boldsymbol{\sigma}^{\*}\_{\boldsymbol{z}} \end{array}\right\} \right)^{\mathrm{T}} \left(\left\{\begin{array}{cc} \boldsymbol{\sigma}\_{\boldsymbol{x}} \\ \boldsymbol{\sigma}\_{\boldsymbol{y}} \\ \boldsymbol{\sigma}\_{\boldsymbol{z}} \\ \boldsymbol{\sigma}\_{\boldsymbol{z}} \end{array}\right\} - \left\{\begin{array}{cc} \boldsymbol{\sigma}^{\*}\_{\boldsymbol{x}} \\ \boldsymbol{\sigma}^{\*}\_{\boldsymbol{y}} \\ \boldsymbol{\sigma}^{\*}\_{\boldsymbol{x}} \\ \boldsymbol{\sigma}^{\*}\_{\boldsymbol{z}} \end{array}\right\} \right) d\Omega \end{aligned} \tag{5}$$

whereas the average norm stress error for the whole domain is

$$\begin{aligned} \left\|\boldsymbol{\ell}\right\|^{2} &= \frac{1}{m} \sum\_{\boldsymbol{\ell}^{x}=1}^{m} \int \boldsymbol{\sigma}^{\mathsf{T}} \boldsymbol{\sigma} d\boldsymbol{\ell} \Omega \\ &= \frac{1}{m} \sum\_{\boldsymbol{\ell}^{x}=1}^{m} \int \left\{ \begin{array}{c} \sigma\_{\boldsymbol{x}} \\ \sigma\_{\boldsymbol{y}} \\ \tau\_{\boldsymbol{xy}} \\ \sigma\_{z} \end{array} \right\}^{\mathsf{T}} \left\{ \begin{array}{c} \sigma\_{\boldsymbol{x}} \\ \sigma\_{\boldsymbol{y}} \\ \tau\_{\boldsymbol{xy}} \\ \sigma\_{z} \end{array} \right\} d\Omega \end{aligned} \tag{6}$$

where *m* indicates the total number of elements in the whole domain and *σ*<sup>∗</sup> is the smoothed stress vector. In the finite element treatment the integration with the isoparametric triangular element is converted by the summation of quadratures following the Radau rules [25] as follows:

e<sup>2</sup>*e* = 1 −1 1 −1 ⎛⎜⎜⎜⎝⎧⎪⎪⎪⎨⎪⎪⎪⎩ *<sup>σ</sup>*(*ξ*, *<sup>η</sup>*)*x <sup>σ</sup>*(*ξ*, *<sup>η</sup>*)*y <sup>τ</sup>*(*ξ*, *<sup>η</sup>*)*xy <sup>σ</sup>*(*ξ*, *<sup>η</sup>*)*z* ⎫⎪⎪⎪⎬⎪⎪⎪⎭ − ⎧⎪⎪⎪⎨⎪⎪⎪⎩ *<sup>σ</sup>*(*ξ*, *<sup>η</sup>*)<sup>∗</sup>*x <sup>σ</sup>*(*ξ*, *η*)∗*y <sup>τ</sup>*(*ξ*, *<sup>η</sup>*)<sup>∗</sup>*xy <sup>σ</sup>*(*ξ*, *<sup>η</sup>*)<sup>∗</sup>*z* ⎫⎪⎪⎪⎬⎪⎪⎪⎭⎞⎟⎟⎟⎠T⎛⎜⎜⎜⎝⎧⎪⎪⎪⎨⎪⎪⎪⎩ *<sup>σ</sup>*(*ξ*, *<sup>η</sup>*)*x <sup>σ</sup>*(*ξ*, *<sup>η</sup>*)*y <sup>τ</sup>*(*ξ*, *<sup>η</sup>*)*xy <sup>σ</sup>*(*ξ*, *<sup>η</sup>*)*z* ⎫⎪⎪⎪⎬⎪⎪⎪⎭ − ⎧⎪⎪⎪⎨⎪⎪⎪⎩ *<sup>σ</sup>*(*ξ*, *<sup>η</sup>*)<sup>∗</sup>*x <sup>σ</sup>*(*ξ*, *η*)∗*y <sup>τ</sup>*(*ξ*, *<sup>η</sup>*)<sup>∗</sup>*xy <sup>σ</sup>*(*ξ*, *<sup>η</sup>*)<sup>∗</sup>*z* ⎫⎪⎪⎪⎬⎪⎪⎪⎭⎞⎟⎟⎟⎠*t<sup>e</sup>*det*Jedξdη* = 3 ∑ *p*=1 ⎛⎜⎜⎜⎝⎧⎪⎪⎪⎨⎪⎪⎪⎩ *σξ p*, *<sup>η</sup>px σξ p*, *<sup>η</sup>py τξ p*, *<sup>η</sup>pxy σξ p*, *<sup>η</sup>pz* ⎫⎪⎪⎪⎬⎪⎪⎪⎭ − ⎧⎪⎪⎪⎨⎪⎪⎪⎩ *σξ p*, *<sup>η</sup>p*<sup>∗</sup>*x σξ p*, *<sup>η</sup>p*<sup>∗</sup>*y τξ p*, *<sup>η</sup>p*<sup>∗</sup>*xy σξ p*, *<sup>η</sup>p*<sup>∗</sup>*z* ⎫⎪⎪⎪⎬⎪⎪⎪⎭⎞⎟⎟⎟⎠T ⎛⎜⎜⎜⎝⎧⎪⎪⎪⎨⎪⎪⎪⎩ *σξ p*, *<sup>η</sup>px σξ p*, *<sup>η</sup>py τξ p*, *<sup>η</sup>pxy σξ p*, *<sup>η</sup>pz* ⎫⎪⎪⎪⎬⎪⎪⎪⎭ − ⎧⎪⎪⎪⎨⎪⎪⎪⎩ *σξ p*, *<sup>η</sup>p*<sup>∗</sup>*x σξ p*, *<sup>η</sup>p*<sup>∗</sup>*y τξ p*, *<sup>η</sup>p*<sup>∗</sup>*xy σξ p*, *<sup>η</sup>p*<sup>∗</sup>*z* ⎫⎪⎪⎪⎬⎪⎪⎪⎭⎞⎟⎟⎟⎠*t<sup>e</sup>*det*JeWp* (7)

where *WP* is a weighting factor and is *Je* is the Jacobian matrix. Similarly,

*e*<sup>ˆ</sup><sup>2</sup> = 1*m m*∑ *<sup>e</sup>*=1 3 ∑ *p*=1 ⎛⎜⎜⎜⎝⎧⎪⎪⎪⎨⎪⎪⎪⎩ *σξ p*, *<sup>η</sup>px σξ p*, *<sup>η</sup>py τξ p*, *<sup>η</sup>pxy σξ p*, *<sup>η</sup>pz* ⎫⎪⎪⎪⎬⎪⎪⎪⎭⎞⎟⎟⎟⎠T ⎛⎜⎜⎜⎝⎧⎪⎪⎪⎨⎪⎪⎪⎩ *σξ p*, *<sup>η</sup>px σξ p*, *<sup>η</sup>py τξ p*, *<sup>η</sup>pxy σξ p*, *<sup>η</sup>pz* ⎫⎪⎪⎪⎬⎪⎪⎪⎭⎞⎟⎟⎟⎠*t<sup>e</sup>*det*JeWp* (8)

where *te* is the element thickness for a plane stress condition and *te* = 1 for a plane strain condition. Therefore, the relative stress norm error *ζe* for each item is considerably less than some identified value [26]. Thus,

$$\mathcal{L}\_{\varepsilon} = \frac{||\varepsilon||\_{\varepsilon}}{||\mathbb{1}||} \le \zeta \tag{9}$$

And the relative stress error level of the new element is defined as permissible error as

$$\varepsilon\_{\mathcal{E}} = \frac{||\mathcal{E}||\_{\mathcal{E}}}{\mathcal{J}||\mathcal{E}||} \le 1 \tag{10}$$

This implies that any element with *εe* > 1 must be optimized and the new mesh size must be predicted. The asymptotic convergence rate criteria are used, which assumed

$$||e||\_{\mathcal{C}} \propto h\_{\mathcal{C}}^{p} \tag{11}$$

where *p* is the approximation of the polynomial order. In the analysis, *p* = 2 is used for the approximation of finite elements as a quadratic polynomial. The predicted sizes of the new element are stated as follows:

$$h\_N = \frac{1}{\sqrt{\varepsilon\_{t'}}} h\_{t'} \tag{12}$$

where *he* is the old element size and *p* is the order of the interpolation shape function.

Convergence of the mesh is dependent of the size of the new element, which defines how many elements in a model are required to ensure that the results of an analysis are not affected by changing the mesh size. System response (stress, deformation) converge with decreasing element size to a repeatable solution. Further refinement of the mesh does not affect results because the model and its results are now independent of the mesh.

The present mesh is known as the new background mesh and the advancing front method is replicated according to the amount of mesh refinements set by the user.

The mesh optimization is used in the final stage of the mesh generation in order to enhance the shape of the elements. The topological structure of the mesh is fixed in the process of mesh smoothing, i.e., the element's nodal connections are not changed, but the inner nodes are repositioned to create triangles with much improved shapes. The most effective computational smoothing algorithm is the well-known Laplacian smoothing [27], which repositions the inner node created by its neighboring nodes at the center of the polygon. The new position of an internal node *i* is computed as

$$(\mathbf{x}\_{i\prime}y\_i) = \frac{1}{N\_n} \sum\_{j=1}^{N} (\mathbf{x}\_{j\prime}y\_j) \tag{13}$$

where *Nn* is the number of nodes linked to node *i*. The mesh smoothing process consists of several iterations.
