1 Introduction

The finite element method (FEM) has been rapidly developed since the 1950s when it was first introduced and successfully commercialized. Compared with other numerical analysis methods, the finite element method is more widely used in actual engineering fields. A key factor in its success is the use of element meshes, which allow the ease of handling complex geometries and boundary conditions. In addition, as the mesh becomes finer, the finite element solutions become more accurate; hence, this convergence characteristic is a crucial feature for ensuring the reliability of solutions.

Ideally, the convergence of finite element (FE) solutions should only depend on the mesh density. However, actual FE solutions rely on domain geometry and material properties, and under certain conditions of the foregoing, the solution convergence can deteriorate; this is called locking. There are various types such as shear locking, membrane locking, volumetric locking and so on. Over the past 50 years, efforts have been devoted to develop various methods, such as reduced integration [1,2,3,4], incompatible mode [5, 6], assumed strain, and enhanced strain [7,8,9] techniques, for alleviating locking. However, the effectiveness of these methods rapidly diminishes as the element meshes become more distorted. The influence of this distortion becomes even more severe on low-order elements, which are generally preferred in engineering practice.

This study proposes a completely different approach for improving FE solutions. We devise a new type of finite element called “self-updated finite element” (SUFE) and a method that considerably increases the accuracy of FE solutions through iterative analysis without mesh refinement. In conventional FE analysis, the stiffness matrix of an element is determined by element geometry and material properties. In the proposed technique, the key idea is to consider element deformation when calculating the stiffness matrix. First, the deformation is estimated through an initial analysis; then, the improved stiffness matrix reflecting this deformation is calculated. Thereafter, through a second analysis, more accurate deformation is derived. By repeating the foregoing process, the FE solution accuracy improves.

The proposed concept is implemented for a two-dimensional four-node plane-stress finite element. The objective is to improve the solutions by minimizing the shear locking effect on the four-node element through iterative analysis. Accordingly, a mode-based FE formulation is adopted. The four-node element has three rigid-body modes and five deformation modes. Among the deformation modes, two different bending modes can be chosen according to the bending direction. For the bending modes, analytical bending strains are enforced using the assumed strain method. For a given element geometry, material properties, and element deformation, the optimal bending directions are determined such that the element strain energy is minimized. To resolve the time-consuming problem of searching the optimal bending directions, deep learning is employed. The four-node self-updated finite element passes the patch and zero-energy mode tests. The main advantage of this approach is that highly accurate FE solutions can be derived even when the meshes are coarse and severely distorted.

In Sect. 2, the basic formulation of the standard four-node finite element is briefly reviewed. The mode-based finite element formulation used for the SUFE is introduced in Sect. 3. The iterative solution procedure for the SUFE is explained in Sect. 4. In Sect. 5, the use of deep neural network to improve the numerical efficiency of the iterative procedure is described. The results of various numerical tests are elaborated in Sect. 6. Finally, the conclusions are presented in Sect. 7.

2 Isoparametric finite element formulation

The basic formulation of the standard four-node finite is reviewed in this section [10].

Consider a mesh of four-node finite elements in the global Cartesian coordinate system defined by \({\mathbf{i}}_{x}\) and \({\mathbf{i}}_{y}\), as shown in Fig. 1a. The geometry of a four-node finite element \(m\) is interpolated by

$$ {\mathbf{x}} = \sum\limits_{i = 1}^{4} {h_{i} (r,s){\mathbf{x}}_{i} } \,{\text{with}}\;{\mathbf{x}} = \left[ {\begin{array}{*{20}c} x & y \\ \end{array} } \right]^{{\text{T}}} \,{\text{and}}\,{\mathbf{x}}_{i} = \left[ {\begin{array}{*{20}c} {x_{i} } & {y_{i} } \\ \end{array} } \right]^{{\text{T}}}, $$
(1)

where \(r\) and \(s\) are the natural coordinates in Fig. 1b, \(h_{i} (r,\,s)\) is the two-dimensional interpolation function of the standard isoparametric procedure corresponding to node i, and \({\mathbf{x}}_{i}\) is the position vector of node i in the global Cartesian coordinate system.

Fig. 1
figure 1

Four-node finite element in a global Cartesian coordinate system and b natural coordinate system

Full size image

The corresponding displacement interpolation of the element is

$$ \left[ {\begin{array}{*{20}c} u \\ v \\ \end{array} } \right] = \sum\limits_{i = 1}^{4} {h_{i} (r,s)\left[ {\begin{array}{*{20}c} {u_{i} } \\ {v_{i} } \\ \end{array} } \right]} , $$
(2)

where \(u\) and \(v\) are the displacements in the directions of \({\mathbf{i}}_{x}\) and \({\mathbf{i}}_{y}\), respectively, and \(u_{i}\) and \(v_{i}\) are the nodal displacements at node i, respectively.

The derivatives of displacements with respect to the global coordinates are calculated using the Jacobian matrix, \({\mathbf{J}}\):

$$ \left[ {\begin{array}{*{20}c} {\frac{\partial u}{{\partial x}}} \\ {\frac{\partial u}{{\partial y}}} \\ \end{array} } \right] = {\mathbf{J}}^{ - 1} \left[ {\begin{array}{*{20}c} {\frac{\partial u}{{\partial r}}} \\ {\frac{\partial u}{{\partial s}}} \\ \end{array} } \right]\,{\text{and}}\,\left[ {\begin{array}{*{20}c} {\frac{\partial v}{{\partial x}}} \\ {\frac{\partial v}{{\partial y}}} \\ \end{array} } \right] = {\mathbf{J}}^{ - 1} \left[ {\begin{array}{*{20}c} {\frac{\partial v}{{\partial r}}} \\ {\frac{\partial v}{{\partial s}}} \\ \end{array} } \right]\,{\text{with}}\, \, {\mathbf{J}} = \left[ {\begin{array}{*{20}c} {\frac{\partial x}{{\partial r}}} & {\frac{\partial y}{{\partial r}}} \\ {\frac{\partial x}{{\partial s}}} & {\frac{\partial y}{{\partial s}}} \\ \end{array} } \right]. $$
(3)

The derivatives in Eq. (3) are utilized to obtain the strain-displacement matrix, \({\overline{\mathbf{B}}}_{m}\):

$$ {{\varvec{\upvarepsilon}}} = {\overline{\mathbf{B}}}_{m} (r,s){\mathbf{u}}, $$
(4)
$$ {{\varvec{\upvarepsilon}}} = \left[ {\begin{array}{*{20}c} {\varepsilon_{xx} } & {\varepsilon_{yy} } & {\gamma_{xy} } \\ \end{array} } \right]^{{\text{T}}} \,{\text{with}}\,\gamma_{xy} = 2\varepsilon_{xy} , $$
(5)
$$ {\varvec{\mathbf{u}}} = \left[ {\begin{array}{*{20}c} {u_{1} } & {u_{2} } & {u_{3} } & {u_{4} } & {v_{1} } & {v_{2} } & {v_{3} } & {v_{4} } \\ \end{array} } \right]^{\text{T}} , $$
(6)

where \({{\varvec{\upvarepsilon}}}\) is the strain vector, and \({\mathbf{u}}\) is the nodal displacement vector of element \(m\).

The stiffness matrix of the standard four-node finite element \(m\) is then obtained as

$$ \overline{{\mathbf{K}}}_{m} = \int_{V} {\overline{{\mathbf{B}}}_{m}^{\text{T}} {\mathbf{C}}\overline{{\mathbf{B}}}_{m} dV} , $$
(7)

where \(V\) is the element volume, and \({\mathbf{C}}\) is the material law matrix.

The equilibrium equation of an FE model is given by

$$ {\mathbf{\overline{K}U = R}}\, \;{\text{with}}\,\,\,{\overline{\mathbf{K}}} = {\text{A}}_{m = 1}^{N} {\overline{\mathbf{K}}}_{m} , $$
(8)

where \({\overline{\mathbf{K}}}\), \({\mathbf{U}}\), and \({\mathbf{R}}\) denote the stiffness matrix, displacement vector, and load vector for the entire FE model, respectively; \({\text{A}}_{m = 1}^{N}\) is the FE assembly operator; and \(N\) is the number of finite elements in the FE model.

3 Mode-based finite element formulation

The four-node finite element has eight nodal degrees of freedom (DOFs). The number of kinematic modes that can be represented by the finite element is equal to the number of DOFs. The kinematic modes shown in Fig. 2 can be categorized into three groups, as follows:

  • Three rigid body (called zero strain) modes: two translational (\({\mathbf{\varphi }}_{1}\) and \({\mathbf{\varphi }}_{2}\)) and one rotational (\({\mathbf{\varphi }}_{3}\)) modes

  • Three constant strain modes: two stretching modes (\({\mathbf{\varphi }}_{4}\) and \({\mathbf{\varphi }}_{5}\)) and one shearing (\({\mathbf{\varphi }}_{6}\)) mode

  • Two linear strain modes: two bending modes (\({\mathbf{\varphi }}_{7}\) and \({\mathbf{\varphi }}_{8}\))

Fig. 2
figure 2

Kinematic modes of the four-node finite element in Fig. 1: a three rigid body modes (\({\mathbf{\varphi }}_{1}\), \({\mathbf{\varphi }}_{2}\), and \({\mathbf{\varphi }}_{3}\)), b three constant strain modes (\({\mathbf{\varphi }}_{4}\), \({\mathbf{\varphi }}_{5}\), and \({\mathbf{\varphi }}_{6}\)), and c two linear strain modes (\({\mathbf{\varphi }}_{7}\) and \({\mathbf{\varphi }}_{8}\)); solid blue lines indicate kinematic mode shapes

Full size image

A kinematic mode-based formulation for the four-node finite element incorporated with the assumed strain method is presented in this section.

3.1 Modal representation of strain–displacement matrix

First, two translational and one rotational rigid body modes can be computed through the following form:

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{\varphi }}_{1} } & {{\mathbf{\varphi }}_{2} } & {{\mathbf{\varphi }}_{3} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 \\ { - y_{1} } & { - y_{2} } & { - y_{3} } & { - y_{4} } & {x_{1} } & {x_{2} } & {x_{3} } & {x_{4} } \\ \end{array} } \right]^{{\text{T}}} , $$
(9)

where \(x_{i}\) and \(y_{i}\) are the coordinates of node i in the global Cartesian coordinate system. Note that the nodal components of the translational mode \({\mathbf{\varphi }}_{1}\) are obtained when the displacement fields are u = 1 and v = 0. Similarly, the translational mode \({\mathbf{\varphi }}_{2}\) is obtained when u = 0 and v = 1, and the rotational mode \({\mathbf{\varphi }}_{3}\) is obtained when u =  − y and v = x.

Three constant strain modes (two stretching modes and one shearing mode) can be similarly obtained as

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{\varphi }}_{4} } & {{\mathbf{\varphi }}_{5} } & {{\mathbf{\varphi }}_{6} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {x_{1} } & {x_{2} } & {x_{3} } & {x_{4} } & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & {y_{1} } & {y_{2} } & {y_{3} } & {y_{4} } \\ {y_{1} } & {y_{2} } & {y_{3} } & {y_{4} } & 0 & 0 & 0 & 0 \\ \end{array} } \right]^{{\text{T}}} , $$
(10)

where the components of the three constant strain modes are obtained when the displacement field is given by \(u = x\) and \(v = 0\) for the stretching mode \({\mathbf{\varphi }}_{4}\); \(u = 0\) and \(v = y\) for the stretching mode \({\mathbf{\varphi }}_{5}\); \(u = y\) and \(v = 0\) for the shearing mode \({\mathbf{\varphi }}_{6}\).

Two linear strain modes (two bending modes) can be calculated by

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{\varphi }}_{7} } & {{\mathbf{\varphi }}_{8} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {x_{1} y_{1} } & {...} & {x_{4} y_{4} } & { - \frac{{x_{1}^{2} + \upsilon y_{1}^{2} }}{2}} & {...} & { - \frac{{x_{4}^{2} + \upsilon y_{4}^{2} }}{2}} \\ { - \frac{{\upsilon x_{1}^{2} + y_{1}^{2} }}{2}} & {...} & { - \frac{{\upsilon x_{4}^{2} + y_{4}^{2} }}{2}} & {x_{1} y_{1} } & {...} & {x_{4} y_{4} } \\ \end{array} } \right]^{{\text{T}}} , $$
(11)

where \(\upsilon\) is Poisson’s ratio. The bending mode components are obtained when the displacement field is given by

$$ u = xy\,{\text{and}}\,v = - \frac{1}{2}\left( {x^{2} + \upsilon y^{2} } \right)\,{\text{for the bending mode}}\,{\mathbf{\varphi }}_{7}; $$
(12)
$$ u = - \frac{1}{2}\left( {\upsilon x^{2} + y^{2} } \right)\,{\text{and}}\,v = xy\,{\text{for the bending mode}}\,{\mathbf{\varphi }}_{8}. $$
(13)

The nodal displacement vector of the four-node finite element can be expressed using the eight kinematic modes in Eqs. (9)–(11)

$$ {\mathbf{u}} = {{\varvec{\Phi}}}\,{\mathbf{q}} $$
(14)
$$ {\text{with}}\,\,{{\varvec{\Phi}}}\, = \left[ {\begin{array}{*{20}c} {{\mathbf{\varphi }}_{1} } & {{\mathbf{\varphi }}_{2} } & {...} & {{\mathbf{\varphi }}_{8} } \\ \end{array} } \right]\,{\text{and}}\,{\mathbf{q}} = \left[ {\begin{array}{*{20}c} {q_{1} } & {q_{2} } & {...} & {q_{8} } \\ \end{array} } \right]^{\text{T}} , $$
(15)

where \({{\varvec{\Phi}}}\) is the kinematic mode matrix, and \({\mathbf{q}}\) is the generalized coordinate vector. By substituting Eqs. (14) into (4), the strain vector can be expressed as \({\mathbf{q}}\):

$$ {{\varvec{\upvarepsilon}}} = {\overline{\mathbf{B}}}_{m} {\mathbf{u}} = {\overline{\mathbf{B}}}_{m} {\mathbf{\Phi q}} = {\hat{\mathbf{B}}}_{m} {\mathbf{q}}, $$
(16)

in which \({\hat{\mathbf{B}}}\) is the matrix for relating the strains with the kinematic modes, and

$$ {\hat{\mathbf{B}}}_{m} = \left[ {\begin{array}{*{20}c} {{\mathbf{e}}_{1} } & {{\mathbf{e}}_{2} } & {...} & {{\mathbf{e}}_{7} } & {{\mathbf{e}}_{8} } \\ \end{array} } \right]\,{\text{with}}\,{\mathbf{e}}_{i} = {\overline{\mathbf{B}}}_{m} {\mathbf{\varphi }}_{i} \,{\text{for}}\,i = 1,2, \ldots ,8. $$
(17)

The strain–displacement relationship is given by

$$ {{\varvec{\upvarepsilon}}} = {\hat{\mathbf{B}}}_{m} \,{\mathbf{q}} = {\hat{\mathbf{B}}}_{m} \,{{\varvec{\Phi}}}^{ - 1} {\mathbf{u}} = {\overline{\mathbf{B}}}_{m} {\mathbf{u}}. $$
(18)

3.2 Assumed modal strain

It has been well known that shear locking occurs in bending. To improve the solutions by minimizing shear locking, the strain vectors, \({\mathbf{e}}^{7}\) and \({\mathbf{e}}^{8}\) in \({\hat{\mathbf{B}}}_{m}\), corresponding to the two bending modes (\({\mathbf{\varphi }}_{7}\) and \({\mathbf{\varphi }}_{8}\)) in the four-node finite element \(m\) have to be assumed.

In Eq. (11), two bending modes are calculated in the global Cartesian coordinate system. Noting that the bending modes can be chosen differently according to the bending directions is important.

Consider a Cartesian coordinate system defined by two orthonormal base vectors, \({\mathbf{i^{\prime}}}_{x}\) and \({\mathbf{i^{\prime}}}_{y}\) (called bending directions), whose origin is at the same location as that of the global Cartesian coordinate system, as shown in Fig. 3a. The basis vector, \({\mathbf{i^{\prime}}}_{x}\), has an angle, \(\alpha_{m}\), from \({\mathbf{i}}_{x}\). The corresponding bending modes are obtained as follows:

$$ \left[ {\begin{array}{*{20}c} {{\mathbf{\varphi }}_{7}^{\prime } } & {{\mathbf{\varphi }}_{8}^{\prime } } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {x_{1}^{\prime } y_{1}^{\prime } } & {...} & {x_{4}^{\prime } y_{4}^{\prime } } & { - \frac{{x_{1}^{\prime 2} + \nu y_{1}^{\prime 2} }}{2}} & {...} & { - \frac{{x_{4}^{\prime 2} + \nu y_{4}^{\prime 2} }}{2}} \\ { - \frac{{\nu x_{1}^{\prime 2} + y_{1}^{\prime 2} }}{2}} & {...} & { - \frac{{\nu x_{4}^{\prime 2} + y_{4}^{\prime 2} }}{2}} & {x^{\prime}_{1} y^{\prime}_{1} } & {...} & {x_{4}^{\prime } y_{4}^{\prime } } \\ \end{array} } \right]^{\text{T}} , $$
(19)

where \(x^{\prime}\) and \(y^{\prime}\) are the coordinates in the new Cartesian coordinate system defined by \({\mathbf{i^{\prime}}}_{x}\) and \({\mathbf{i^{\prime}}}_{y}\), respectively.

Fig. 3
figure 3

Bending directions and assumed deformed shapes: a bending direction angle, \(\alpha_{m}\), and b assumed deformed shapes for bending modes \({\mathbf{\varphi^{\prime}}}_{7}\) and \({\mathbf{\varphi^{\prime}}}_{8}\)

Full size image

The nodal components of the bending modes in Eq. (19) are obtained as follows:

$$ u^{\prime} = x^{\prime}y^{\prime}\,{\text{and}}\,v^{\prime} = - \frac{1}{2}\left( {x^{{\prime}{2}} + \upsilon y^{{\prime}{2}} } \right)\,{\text{for bending mode}}\,{\mathbf{\varphi^{\prime}}}_{7}; $$
(20)
$$ u^{\prime} = - \frac{1}{2}\left( {\upsilon x^{{\prime}{2}} + y^{{\prime}{2}} } \right)\,{\text{and}}\,v^{\prime} = x^{\prime}y^{\prime}\,{\text{for bending mode}}\,{\mathbf{\varphi^{\prime}}}_{8}, $$
(21)

where \(u^{\prime}\) and \(v^{\prime}\) are the displacements in the \(x^{\prime}\) and \(y^{\prime}\) directions, respectively. The assumed deformed shape of the element domain obtained using Eqs. (20) and (21) is shown in Fig. 3b. Note that the deformed shape cannot be represented by the displacement interpolation of the standard four-node finite element in Eq. (2).

The corresponding analytical bending strains are employed as the assumed strains (\({\mathbf{e^{\prime}}}_{7}\) and \({\mathbf{e^{\prime}}}_{8}\)) for the bending modes (\({\mathbf{\varphi^{\prime}}}_{7}\) and \({\mathbf{\varphi^{\prime}}}_{8}\)):

$$ {\mathbf{e}}_{7}^{\prime } = {\mathbf{T}}\left[ {\begin{array}{*{20}c} {y^{\prime } } & { - \upsilon y^{\prime } } & 0 \\ \end{array} } \right]^{{\text{T}}} \,{\text{and}}\,{\mathbf{e}}_{8}^{\prime } = {\mathbf{T}}\left[ {\begin{array}{*{20}c} { - \upsilon x^{\prime } } & {x^{\prime } } & 0 \\ \end{array} } \right]^{{\text{T}}} $$
(22)
$$ {\text{with}}\,{\mathbf{T}} = \left[ {\begin{array}{*{20}c} {\cos^{2} \alpha_{m} } & {\sin^{2} \alpha_{m} } & { - \sin \alpha_{m} \cos \alpha_{m} } \\ {\sin^{2} \alpha_{m} } & {\cos^{2} \alpha_{m} } & {\sin \alpha_{m} \cos \alpha_{m} } \\ {2\sin \alpha_{m} \cos \alpha_{m} } & { - 2\sin \alpha_{m} \cos \alpha_{m} } & {\cos^{2} \alpha_{m} - \sin^{2} \alpha_{m} } \\ \end{array} } \right], $$
(23)

where \({\mathbf{T}}\) is the strain transformation matrix [11].

Finally, the strain-displacement matrix with the assumed modal strains in Eq. (22) is obtained as follows:

$$ {\mathbf{B}}_{m} {\mathbf{ = \hat{B}^{\prime}}}_{m} \,{\mathbf{\Phi^{\prime}}}^{ - 1} \, $$
(24)
$$ {\text{with}}\,{\mathbf{\hat{B}^{\prime}}}_{m} = \left[ {\begin{array}{*{20}c} {{\mathbf{e}}_{1} } & {...} & {{\mathbf{e}}_{6} } & {{\mathbf{e^{\prime}}}_{7} } & {{\mathbf{e^{\prime}}}_{8} } \\ \end{array} } \right]\,\,{\text{and}}\,{\mathbf{\Phi^{\prime}}}\, = \left[ {\begin{array}{*{20}c} {{\mathbf{\varphi }}_{1} } & {...} & {{\mathbf{\varphi }}_{6} } & {{\mathbf{\varphi^{\prime}}}_{7} } & {{\mathbf{\varphi^{\prime}}}_{8} } \\ \end{array} } \right]. $$
(25)

The stiffness matrix of the four-node finite element adopting the assumed modal strains is obtained as follows:

$$ {\mathbf{K}}_{m} = \int\limits_{V} {{\mathbf{B}}_{m}^{\text{T}} {\mathbf{CB}}_{m} dV}. $$
(26)

The element stiffness matrix in Eq. (26) depends on the bending directions determined by angle \(\alpha_{m}\).

The stiffness matrix of the FE model is assembled using the direct stiffness method:

$$ {\mathbf{K}} = {\text{A}}_{m = 1}^{N} {\mathbf{K}}_{m} . $$
(27)

4 Iterative solution procedure

The solution accuracy of the standard four-node finite element deteriorates because of shear locking when the displacements contain bending modes. Furthermore, the locking effect becomes more severe when the element mesh is distorted [12,13,14,15,16,17]. In this section, an iterative procedure to alleviate shear locking is introduced.

The entire solution procedure that adopts the SUFE can be divided into two stages. The preliminary stage determines whether or not the iterative solution procedure proceeds. Once the FE solution does not contain bending modes, the iterative procedure is unnecessary. The determination process is implemented as follows.

The equilibrium equation in Eq. (8) is solved using the standard four-node finite elements to obtain the displacement solution, \({\overline{\mathbf{U}}}\), and the corresponding strain energy,

$$ \overline{E} = 0.5 \cdot {\overline{\mathbf{R}}}^{{\text{T}}} {\overline{\mathbf{U}}}. $$
(28)

Using the stiffness matrix in Eq. (27) and the solution (\({\overline{\mathbf{U}}}\)), the strain energy is calculated as

$$ E = 0.5 \cdot {\overline{\mathbf{U}}}^{{\text{T}}} {\mathbf{K\overline{U}}}, $$
(29)

where \({\mathbf{K}}\) is calculated using \(\alpha_{m} = 0\) for all finite elements (\(m = 1\) to \(N\)).

When \(\overline{E} = E\), the solution displacement, \({\overline{\mathbf{U}}}\), only contains six kinematic modes, \({\mathbf{\varphi }}_{1}\)–\({\mathbf{\varphi }}_{6}\), except for the bending modes. Therefore, any treatment for shear locking is unnecessary, and the iterative solution procedure is not required.

Otherwise, in the iterative stage, for \(i = 1,\,\,2,\,\,3 \ldots\), the following equilibrium equation is repeatedly solved:

$$ {\mathbf{K}}^{(i)} {\mathbf{U}}^{(i)} = {\mathbf{R}}\,{\text{until}}\;\frac{{E^{(i)} - E^{(i - 1)} }}{{E^{(i - 1)} }} < \delta , $$
(30)

where \({\mathbf{K}}^{(0)} = {\mathbf{K}}\), \(E^{(i)} = 0.5 \cdot {\mathbf{R}}^{{\text{T}}} {\mathbf{U}}^{(i)}\), and \(\delta\) is the threshold value.

In iterating Eq. (30), the values of angle \(\alpha {}_{m}\) for all finite elements in the FE model are determined by solving the following optimization problem:

$$ \alpha_{m} = \mathop {\min }\limits_{{\alpha_{m} }} \,{\mathbf{u}}^{{(i - 1)}{\text{T}}} {\mathbf{K}}_{m}^{(i)} {\mathbf{u}}^{(i - 1)} \,{\text{for}}\,m = 1\,{\text{to}}\,N $$
(31)
$$ {\text{with}}\,{\mathbf{K}}_{m}^{(i)} = \int_{V} {{\mathbf{B}}_{m}^{{(i){\text{T}}}} {\mathbf{CB}}_{m}^{(i)} dV} , $$
(32)

where angle \(\alpha_{m}\) (\(0^\circ \le \alpha_{m} < 90^\circ\)) determines the optimal bending directions of element m, and \({\mathbf{u}}^{(i - 1)}\) is obtained from the solution displacement, \({\mathbf{U}}^{(i - 1)}\). This means that the optimal bending directions minimize the strain energy stored in the element under bending; consequently, shear locking is minimized.

The strain-displacement matrix in Eq. (32) is calculated using Eq. (24). Note that when the kinematic mode matrix, \({\mathbf{\Phi^{\prime}}}\), in Eq. (24) is singular, the strain–displacement matrix, \({\mathbf{B}}_{m}\), cannot be assembled. This can easily be resolved by adding a small value (e.g., 0.001) to \(\alpha_{m}\).

The entire solution procedure using the self-updated four-node finite element is shown in Fig. 4.

Fig. 4
figure 4

Flowchart of the iterative solution procedure for SUFE

Full size image

5 Deep learning for estimating the optimal bending directions

Solving the optimization problem described in Sect. 4 is the most time-consuming step in the iterative solution procedure for SUFE. To considerably reduce the computational time required to solve the optimization problem, a neural network can be used to approximate the result [18]. Here, the neural network is constructed to determine a suitable bending direction angle, \(\alpha {}_{m}\). Recently, applications of neural networks to the finite element method have been studied [19,20,21].

This section introduces the methodology for constructing the neural network model in detail including data generation, network configuration, and training. The optimal bending direction angle is approximated through the trained neural network, and the angle is employed to update the stiffness matrix of SUFE.

5.1 Generation of training datasets

Numerous training datasets, including random geometries, displacements, and material properties, are required to train the neural network to accurately approximate the optimal bending direction angle (\(\alpha {}_{m}\)) of arbitrary four-node finite elements. Because the processing of all data spaces is extremely difficult, a previously developed methodology is applied [21] to reduce the amount of data required to achieve efficient network training. Note that deep learning is used to generate stiffness matrices of finite elements in Ref. [21].

The geometry of a four-node element is determined by its four nodal positions. A normalized geometry is adopted to represent all similar shapes [21]. Here, the normalized geometry refers to a quadrilateral where the two nodes (node 1 and 2) at either end of the longest side are fixed at (0, 0) and (1, 0) in the two-dimensional (2D) Cartesian coordinate system. To generate the nth geometry, two other nodes (nodes 3 and 4) are placed randomly, as shown in Fig. 5. In the geometry generation, extremely distorted geometries are excluded. Excluded geometries are quadrilaterals whose interior angle is less than 1° or greater than 179° and those whose ratios between the maximum and minimum side lengths exceed 100. Then, the normalized geometry is generated by four values (i.e., the coordinates of nodes 3 and 4): \(x_{3}\), \(y_{3}\), \(x_{4}\), and \(y_{4}\) [21].

Fig. 5
figure 5

Random generation of the nth normalized element geometry

Full size image

The optimal bending directions depend on the nodal displacements representing the deformed shapes of the finite element. Normalized nodal displacements are considered to represent all similar deformed shapes. Consider an element subjected to the minimum number of constraints to prevent rigid body motions, as shown in Fig. 6a. Then only five nodal displacements (i.e., \(u_{2}\), \(u_{3}\), \(v_{3}\), \(u_{4}\), and \(v_{4}\)) can represent the entire deformed shape. The nth normalized nodal displacement is generated by

$$ u_{2} = \frac{{u^{\prime}_{2} }}{{\delta^{\prime}_{\max } }},\,u_{3} = \frac{{u^{\prime}_{3} }}{{\delta^{\prime}_{\max } }},\,v_{3} = \frac{{v^{\prime}_{3} }}{{\delta^{\prime}_{\max } }},\,u_{4} = \frac{{u^{\prime}_{4} }}{{\delta^{\prime}_{\max } }},{\text{and}}\,v_{4} = \frac{{v^{\prime}_{4} }}{{\delta^{\prime}_{\max } }}, $$
(33)

where \(u^{\prime}_{2}\), \(u^{\prime}_{3}\), \(v^{\prime}_{3}\), \(u^{\prime}_{4}\), and \(v^{\prime}_{4}\) are uniformly generated random numbers ranging from −1.0 to 1.0, and \(\delta^{\prime}_{\max }\) is the largest among the absolute values of the five generated random numbers.

Fig. 6
figure 6

Random generation of normalized nodal displacements for the nth element geometry: a randomly generated displacements for an element with minimum number of constraints to remove rigid body motions; b normalization of generated displacements

Full size image

Poisson’s ratio, \(\upsilon\), is also randomly generated with a uniform distribution in the range 0–0.499999999, and Young’s modulus is arbitrarily chosen because it does not affect the determination of the optimal bending directions.

The optimal bending direction angle, \(\alpha^{(n)}\), is required as labeled data for the generated nth geometry, displacements, and material properties. The angle is calculated by solving the optimization problem in Eq. (31). The well-known Nelder–Mead method [22] is employed, and 30 initial seeds evenly distributed in the range 0°–90° of angle \(\alpha\) are considered to find the global optimum.

Finally, the nth training dataset, \({\mathbf{D}}^{(n)}\), is configured as

$$ {\mathbf{D}}^{(n)} = \left[ {\begin{array}{*{20}c} {{\mathbf{G}}^{(n)} } & {\alpha^{(n)} } \\ \end{array} } \right] $$
(34)

with

$$ {\mathbf{G}}^{(n)} = \left[ {\begin{array}{*{20}c} {x_{3} } & {y_{3} } & {x_{4} } & {y_{4} } & {u_{2} } & {u_{3} } & {v_{3} } & {u_{4} } & {v_{4} } & \upsilon \\ \end{array} } \right]^{(n)}. $$
(35)

5.2 Network configuration and training

The neural network model architecture with 10 inputs (four nodal coordinates, five nodal displacements, and one Poisson’s ratio), one output (\(\alpha ({{\varvec{\uptheta}}})\)), and ten hidden layers is shown in Fig. 7. The angle \(\alpha^{(n)}\) of the training dataset is used for the cost function \(C({{\varvec{\uptheta}}})\):

$$ C({{\varvec{\uptheta}}}) = \frac{1}{M}\sum\limits_{n = 1}^{M} {\left| {\alpha ({{\varvec{\uptheta}}}) - \alpha^{(n)} } \right|^{\frac{1}{2}} } , $$
(36)

where \({{\varvec{\uptheta}}}\) denotes the network weights, and \(M\) is the total number of training datasets.

Fig. 7
figure 7

Network configuration used for determining the bending direction angle (FC: fully connected layer; BN: batch normalization; ELU: exponential linear unit)

Full size image

The deep neural network shown in Fig. 7 is modeled using a fully connected neural network with 10 layers; The neural network depth is 320, and batch normalization is applied to each layer. The effect of the network structure and data size on the cost function values of the training and test datasets is shown in Table 1. As the depth and width of the network and training data size increase, the cost function values decrease. In this study, the network structure shown in Fig. 7 was adopted in consideration of efficiency.

Table 1 Cost function values of the training and test datasets according to the depth and width of the neural network and training data size; The batch size is used as half the training data size
Full size table

To account for the nonlinear behavior of the training datasets, an exponential linear unit used as the activation function is added to each neuron. The effect of the activation function type on the cost function values of the training and test datasets is shown in Table 2. The results show that using the activation function results in smaller cost function values, but there is no significant difference in the results for the activation function type.

Table 2 Cost function values of the training and test datasets according to the activation function type using the network in Fig. 7; The network is trained with a total of 100,000 datasets (M = 100,000), and the batch size is used as half of the training data size
Full size table

To correctly infer the value of \(\alpha_{m}\), 3,000,000 training datasets and 50,000 test datasets are utilized. Histograms of \(\alpha^{(n)}\) in the training and test datasets are shown in Fig. 8. The distribution of the two datasets is similar.

Fig. 8
figure 8

Histograms of \(\alpha\) in the a training and b test datasets

Full size image

TensorFlow [23] is used as the application programming interface to model and train the neural network. Adam optimization [24] is adopted as the optimizer, and Xavier initializer [25] is applied to randomize the initial weights of the neural network. Training is executed for 30,000 epochs, and a batch size of 50,000 is used. The learning rate linearly converged from 0.01 to 0 as the epoch progressed to approach better results. As a result of training, the average errors of the training and test datasets are 0.274° and 0.438°, respectively. All training and evaluation calculations are performed using NVIDIA Tesla V100 (5120 CUDA cores and 32 GB memory).

5.3 Estimation of the bending direction angle from the trained network

Normalized geometries, normalized displacements, and Poisson’s ratios were employed for efficient network training. To use the trained neural network for the iterative solution procedure proposed in Sect. 4, preprocessing the geometry and displacements of the four-node finite element \(m\) is required to obtain network inputs, as shown in Fig. 9.

Fig. 9
figure 9

Preprocessing to obtain network inputs: a original and b normalized element geometries

Full size image

The geometry normalization shown in the figure is performed for the four-node element used in the FE model. The FE connectivity is assigned such that the side length between nodes 1 and 2 is the largest. Then, the normalized nodal coordinates \({\overline{\mathbf{x}}}_{i}\) (\(\overline{x}_{i}\) and \(\overline{y}_{i}\)) as inputs of the neural network are obtained by translating, rotating, and resizing the nodal coordinates, \({\mathbf{x}}_{i}\), of the finite element \(m\) [21]:

$$ {\overline{\mathbf{x}}}_{i} = \frac{{{\mathbf{R}}^{{\text{T}}} \left( {{\mathbf{x}}_{i} - {\mathbf{x}}_{1} } \right)}}{{l_{\max } }}\,{\text{for}}\,i = 1,2,3,4 $$
(37)

with \({\mathbf{R}} = \left[ {\begin{array}{*{20}c} {\cos \beta } & { - \sin \beta } \\ {\sin \beta } & {\cos \beta } \\ \end{array} } \right],\) where \({\mathbf{R}}\) is the rotation matrix, \(\beta\) is the angle between the longest side and x-axis, and \(l_{\max }\) is the largest side length.

Consider the eight nodal displacements (\(u_{i}\) and \(v_{i}\)) of the four-node finite element \(m\). The nodal displacements are rotated through angle \(\beta\) to match the displacements with the normalized geometry, as shown in Fig. 10b. Translational and rotational rigid-body motions are eliminated because the optimal bending directions are not affected by rigid body motions, as shown in Fig. 10c and d. Then the resulting nodal displacements are normalized, as shown in Fig. 10e.

Fig. 10
figure 10

Preprocessing of nodal displacements for network inputs: a initial arbitrary nodal displacements, b rotation of element geometry and displacements, c elimination of translational rigid body modes, d elimination of rotational rigid body mode, and e nomalization of nodal displacements

Full size image

According to the explained procedure, the normalized nodal displacements, \(\overline{u}_{i}\) and \(\overline{v}_{i}\), can be calculated by the following equations:

$$ \overline{u}_{i} = \frac{{\hat{u}_{i} }}{{\delta_{\max } }}\,{\text{and}}\,\overline{v}_{i} = \frac{{\hat{v}_{i} }}{{\delta_{\max } }}\,{\text{for}}\,i = 1,2,3,4 $$
(38)

with

$$ \hat{u}_{i} = \tilde{u}_{i} - \tilde{u}_{1} - \overline{y}_{i} \left( {\tilde{v}_{1} - \tilde{v}_{2} } \right), $$
(39)
$$ \hat{v}_{i} = \tilde{v}_{i} - \tilde{v}_{1} + \overline{x}_{i} \left( {\tilde{v}_{1} - \tilde{v}_{2} } \right), $$
(40)
$$ \left[ {\begin{array}{*{20}c} {\tilde{u}_{i} } \\ {\tilde{v}_{i} } \\ \end{array} } \right] = {\mathbf{R}}^{{\text{T}}} \left[ {\begin{array}{*{20}c} {u_{i} } \\ {v_{i} } \\ \end{array} } \right], $$
(41)

where \(x_{i}\) and \(y_{i}\) are the nodal coordinates of element \(m\), and \(\delta_{\max }\) is the largest among the calculated absolute values of \(\hat{u}_{i}\) and \(\hat{v}_{i}\) for all nodes. As a result, \(\overline{u}_{1}\), \(\overline{v}_{1}\), and \(\overline{v}_{2}\) become zero, and the remaining non-zero displacements, \(\overline{u}_{2}\), \(\overline{u}_{3}\), \(\overline{v}_{3}\), \(\overline{u}_{4}\), and \(\overline{v}_{4}\), are used as network inputs.

Poisson’s ratio (\(\upsilon\)) is used as an input without any modification into the neural network. By using the network inputs given by

$$ {\mathbf{G}}_{input} = \left[ {\begin{array}{*{20}c} {\overline{x}_{3} } & {\overline{y}_{3} } & {\overline{x}_{4} } & {\overline{y}_{4} } & {\overline{u}_{2} } & {\overline{u}_{3} } & {\overline{v}_{3} } & {\overline{u}_{4} } & {\overline{v}_{4} } & \upsilon \\ \end{array} } \right], $$
(42)

the output, \(\alpha_{out}\), of the trained network is derived. Finally, the bending direction angle, \(\alpha_{m}\), for element \(m\) is calculated by

$$ \alpha_{m} = \beta + \alpha_{{{\text{out}}}} . $$
(43)

6 Numerical results

In this section, three basic tests are performed for the SUFE to verify its performance through several numerical problems: Cook’s skew beam, MacNeal’s cantilever beam, cantilever beam modeled using two elements with distortion parameter, cantilever beams modeled using distorted elements, thick and thin curved beams, and strap plate. To thoroughly investigate the performance of SUFE, the results are compared with those of several existing elements listed in Table 3. All compared elements pass the patch test. In this section, the self-updated four-node element is called “SU4.” Note that the Q8 element in Table 3 is an eight-node element unlike the others.

Table 3 List of finite elements used for comparison. The Q8 element is an eight-node element unlike the others
Full size table

The smaller \(\delta\), the higher the accuracy, but the more iterations, see Table 5. In the following numerical examples, the iterations are performed using a threshold \(\delta\) of 0.001 by default.

The computations are performed using Python scripts on an AMD Ryzen 3900X @ 3.80 GHz, 64-GB RAM, 64-bit Windows 10 operating system. The linear equations are solved using a sparse solver in the Python package, NumPy [26].

6.1 Basic numerical tests

In this section, we perform three basic tests: patch test, zero-energy mode test, and rotation dependency test. The patch test is passed if constant stress fields are represented exactly considering shearing and stretching in the finite element model in Fig. 11a [10]. The self-updated four-node finite element (SU4) passes the standard patch tests.

Fig. 11
figure 11

Geometries used for a patch test (\(q = 1.0\), \(E = 3.0 \times 10^{7}\), \(\upsilon = 0.3\), and thickness = 1) and b zero-energy mode test (E = 1.5 \(\times\) 103, \(\upsilon\) = 0.3, and thickness = 1)

Full size image

The zero-energy mode test has been widely used to verify that a finite element can correctly represent rigid body motions; the four element geometries shown in Fig. 11b are considered. In calculating the stiffness matrix of an SUFE, the element deformation is used as the input. Accordingly, the zero-energy mode test is performed considering various element deformations. The SU4 element correctly produces only three rigid body modes (two translations and one rotation mode) without any spurious energy mode in all geometry cases regardless of element deformation. That is, the SU4 element passes the zero-energy mode test.

The rotation test is intended to investigate the angular dependency of the finite element solution. The cantilever beam shown in Fig. 12 is considered [27]. It is discretized into two distorted elements with lengths L1 = 8 and L2 = 2; the width is H = 5 and thickness = 1. The plane stress condition is used, and the material properties are Young’s modulus (E) = 100 and Poisson’s ratio (\(\upsilon\)) = 0.3. A moment with P = 0.2 is applied at the free end. The beam geometry is rotated counterclockwise by up to 90° at 10° increments with respect to the x-axis. Table 4 lists the tip displacements at point A as the rotation angle varies. The SUFE provides almost the same results regardless of the rotation angle.

Fig. 12
figure 12

Cantilever beam for rotational frame invariance test: a problem description (L1 = 8, L2 = 2, H = 5, P = 0.2, \(E = 100\), \(\upsilon = 0.3\), and thickness = 1) and beam geometries rotated by b 30° and c 60° relative to x-axis

Full size image
Table 4 Tip displacements at point A of the cantilever beam for rotation dependency test shown in Fig. 11; values in parentheses are final numbers of iterations
Full size table

6.2 Cook’s skew beam

To further examine the SUFE, widely known benchmark problems are solved. The first benchmark test is Cook’s skew beam problem [28]. The geometric dimensions and boundary conditions used are shown in Fig. 13a. The plane stress condition is assumed, and the material properties are Young’s modulus (E) = 1.0 and Poisson’s ratio (\(\upsilon\)) = 1/3. A shear force, P = 1 (1/16 force per unit length), is applied at the right end. The solutions of the skew beam are obtained using N × N element meshes (N = 2, 4, 8, and 16). The reference solution of the vertical displacement at point A is 23.965, as given in Ref. [28].

Fig. 13
figure 13

Cook’s skew beam problem: a problem description (P = 1, E = 1, \(\upsilon = 1/3\), and thickness \(= 1\)) and b FE meshes used

Full size image

To investigate the predictive capability of SUFE in detail, convergence studies are performed. The solution accuracy is measured using the relative strain energy error, which is given as

$$ E_{e} = \left| {\frac{{E_{{{\text{ref}}}} - E_{h} }}{{E_{{{\text{ref}}}} }}} \right|, $$
(44)

where Eref and \(E_{h}\) denote the strain energies stored in the entire structure obtained from the reference and FE solutions, respectively. The optimal convergence of the relative strain energy error is estimated as \(E_{e} \cong ch^{2}\), where c is a constant, and h is the element size. In the convergence curves, h = 1/N is used. The reference solution, \(E_{{{\text{ref}}}}\), is obtained using the 100 × 100 element mesh of standard nine-node quadrilateral elements.

Table 5 lists the deflections at point A (\(v_{A}\)) and their normalized values. The proposed element yields accurate solutions even when the meshes are considerably coarse (N = 1 and 2); Fig. 14a shows the normalized value of \(v_{A}\) with respect to the number of iterations. Within a few iterations, the solution rapidly converges to the reference. The normalized value of \(v_{A}\) according to the number of layers (N) are shown Fig. 14b. The SU4 element has an excellent convergence behavior.

Table 5 Displacements at point A (\(v_{A}\)) for Cook’s skew beam shown in Fig. 13a; values in parentheses are final numbers of iterations for the SU4 element
Full size table
Fig. 14
figure 14

Normalized vertical displacement in Cook’s skew beam with respect to a number of iterations and b number of element layers, \(N\)

Full size image

The convergence curves are shown in Fig. 15a; even when only one iteration is executed, the SU4 element outperforms the other elements. The convergence behavior of strain energy according to the number of iterations is illustrated in Fig. 15b. At the beginning of the iteration, the strain energy rapidly converges to a certain level. As the mesh becomes finer, the convergence with respect to the iterations becomes faster.

Fig. 15
figure 15

Convergence curves for Cook’s skew beam: a relative strain energy error with respect to element size (bold line represents the optimal convergence rate); b relative strain energy error with respect to number of iterations

Full size image

The shear stress contour plots obtained using SU4 elements with N = 8 layers are shown in Fig. 16; the iterations are 0 and 1. The reference stress distribution shown in Fig. 16c is calculated using a 64 × 64 mesh of standard nine-node quadrilateral elements [10]. In the first iteration, the shear stress distribution in the SU4 element rapidly improves.

Fig. 16
figure 16

Shear stress distributions calculated in Cook’s skew beam using the SU4 element with N = 8: a iteration = 0 and b iteration = 1; c reference solution obtained using standard nine-node quadrilateral elements with N = 64

Full size image

6.3 MacNeal’s cantilever beam

Consider MacNeal’s cantilever beam problem [29], which has been frequently employed to test the sensitivity of quadrilateral elements to mesh distortion. The cantilever beam shown in Fig. 17 has a length (L) = 6, width (H) = 0.2, and thickness = 0.1. The plane stress condition is assumed, and the material properties are Young’s modulus (E) = 1.0 × 107 and Poisson’s ratio (\(\upsilon\)) = 0.3. The beam is clamped at the left end, and two loading cases are considered: (I) shear force P = 1 and (II) moment M = 0.2 at the free end. Three different mesh patterns, including regular (Mesh I), parallelogram (Mesh II), and trapezoidal (Mesh III) meshes, are considered.

Fig. 17
figure 17

MacNeal’s beam subjected to two load cases: tip shear force P and tip bending moment M (L = 6, H = 0.2, P = 1, M = 0.2, \(E = 1 \times 10^{7}\), \(\upsilon = 0.3\), and thickness \(= 0.1\)); three different mesh patterns: a regular (Mesh I), b parallelogram (Mesh II), and c trapezoidal (Mesh III) meshes

Full size image

Table 6 summarizes the normalized vertical displacements at point A. The reference solutions are −0.1081 for the shear load case and −0.0054 for the moment load case. Even with only one iteration, the SU4 element exhibits highly accurate solutions in all three mesh patterns regardless of mesh distortion. The SU4 element produces considerably more accurate solutions than the other elements with parallelogram (Mesh II) and trapezoidal (Mesh III) meshes.

Table 6 Normalized vertical displacements at point A for MacNeal’s cantilever beam subjected to either tip shear force or tip moment; values in parentheses are final numbers of iterations for SU4 element
Full size table

6.4 Cantilever beam modeled using two elements with distortion parameter

To further test the mesh distortion sensitivity, modeling the cantilever beam using two elements is considered as shown in Fig. 18. The geometry of the two elements varies with respect to the distortion parameter, e, which reflects the degree of element distortion. The cantilever beam has a length (L) = 10, width (H) = 2, and thickness = 1. The plane stress condition is considered, and the material properties are Young’s modulus (E) = 1500 and Poisson’s ratio (\(\upsilon\)) = 0.25. The cantilever is subjected to a bending moment M = 2000 at the right end. The exact solution of the vertical displacement at point A is 100, as given in Ref. [30].

Fig. 18
figure 18

Cantilever beam modeled using two elements with distortion parameter, e (L = 10, H = 2, M = 2000, \(E = 1500\), \(\upsilon = 0.25\), and thickness = 1)

Full size image

The cantilever is modeled using two four-node elements. As parameter e varies from 0 to 4.9, the two finite elements become more distorted. The normalized displacements at point A are listed in Table 7 and Fig. 19. The SU4 element yields the exact solution regardless of the magnitude of e even with only one iteration, whereas the other elements suffer the effects of mesh distortion. The stress distributions that are calculated using elements Q4, QM6, and SU4 when e = 4.9 are shown in Fig. 20.

Table 7 Normalized vertical displacements at point A in the cantilever modeled using two elements with distortion parameter, e; values in parentheses are final numbers of iterations for the SU4 element
Full size table
Fig. 19
figure 19

Normalized vertical displacements at point A according to distortion parameter, e, in the cantilever beam modeled using two elements; reference solution is \(v_{A}\) = 100 [30]

Full size image
Fig. 20
figure 20

Stress distributions calculated using elements Q4, QM6, and SU4 for the cantilever beam shown in Fig. 19 with distortion parameter (e) = 4.9; reference solutions obtained using standard nine-node quadrilateral elements with 4 × 20 regular mesh

Full size image

6.5 Cantilever beams modeled using distorted elements

As shown in Fig. 21, a cantilever beam is modeled with five distorted quadrilateral elements. Two loading cases are considered: (I) moment M = 2000 and (II) shear force P = 150 at the free tip. The plane stress condition is considered, and the material properties are Young’s modulus (E) = 1500 and Poisson’s ratio (\(\upsilon\)) = 0.25. Table 8 lists the vertical displacements at point A and stresses at point B. The SU4 element provides the exact solutions for load case (I) and highly accurate results for load case (II).

Fig. 21
figure 21

Cantilever beam modeled using five distorted elements (P = 150, M = 2000, \(E\) = 1500, \(\upsilon = 0.25\), and thickness \(= 1\))

Full size image
Table 8 Normalized vertical deflection at point A (\(v_{A}\)) and normalized stress at point B (\(\sigma_{B}\)) for the cantilever beam modeled using five distorted elements; values in parentheses are final numbers of iterations for the SU4 element
Full size table

6.6 Thick and thin curved beams

Thick and thin curved beams fully clamped at the bottom are shown in Fig. 22. For the thick beam case, the geometry is given by R = 10 and H = 5, and the material properties are Young’s modulus (E) = 1000 and Poisson’s ratio (\(\upsilon\)) = 0. A shear force, P = 600, is applied to the free tip.

Fig. 22
figure 22

Thick and thin curved beams: a thick beam (R = 10, H = 5, P = 600, E = 1000, \(\upsilon\) = 0, and thickness \(= 1\)) and b thin beams (P = 1, \(E\) = 365,010 (for case H/R = 0.03) and E = 44,027,109 (for case H/R = 0.006), \(\upsilon\) = 0, and thickness \(= 1\))

Full size image

For the thin beam case, two ratios of thickness to radius (H/R = 0.03 and 0.006) are considered. Poisson’s ratio (\(\upsilon\)) = 0; when H/R = 0.03, Young’s modulus (E) = 365,010, and when H/R = 0.006, E = 44,027,109. A vertical load, P = 1.0, is applied at the free tip. The curved beam is modeled using only four or five elements, as shown in Fig. 22. For both thick and thin beam problems, the plane stress condition is assumed.

Tables 9 and 10 summarize the normalized vertical displacements at point A for the thick and thin beams, respectively. Regardless of the thickness-to-radius ratio, the SU4 element exhibits excellent solutions. However, the solution accuracies of the other elements except for the Q8 element are extremely inadequate, especially those for the thin beam case.

Table 9 Normalized vertical deflection at point A in the thick curved beam shown in Fig. 22a; values in parentheses are final numbers of iterations for the SU4 element
Full size table
Table 10 Vertical displacements at point A in the thin curved beam shown in Fig. 22b; values in parentheses are final numbers of iterations for the SU4 element
Full size table

The bending directions calculated using the SU4 elements for the thin beam case (H/R = 0.03) are also investigated. The bending directions and effective stress (von Mises stress) distributions are shown in Fig. 23. Ideally, the bending direction should be toward the center of curvature of the thin beam. The differences between the calculated and ideal bending directions are less than 0.4°. In addition, the SU4 elements provide a highly accurate stress distribution.

Fig. 23
figure 23

Bending directions calculated by the SU4 elements and effective stress distributions for the thin beam case (H/R = 0.03); green arrows show bending directions, and values right next to arrows indicate differences between the calculated and ideal bending directions; reference stress distribution obtained using standard nine-node quadrilateral elements

Full size image

6.7 Strap plate problem

The strap plate problem is shown in Fig. 24. The strap plate has three holes, two of which are fixed. A uniformly distributed load, q = 10 (force per unit length), in the negative y-direction is applied to the right hole of the plate. The plane stress condition is considered with Poisson’s ratio (\(\upsilon\)) = 0.3 and Young’s modulus (E) = 2 \(\times\) 1011. The solutions are obtained using coarse, medium, and fine meshes. The reference solutions are calculated using the reference mesh of standard nine-node quadrilateral elements with the number of element layers \(N = 50\) (109,582 DOFs), as shown in Fig. 24b.

Fig. 24
figure 24

Strap plate problem: a problem description (q = 10 force per unit length, \(E = 2 \times 10^{11}\), \(\upsilon\) = 0.3, and thickness = 0.01), b reference mesh and c element meshes used

Full size image

The computational efficiency of the SU4 element is investigated. Figure 25a shows the relation between the computation time versus the relative strain energy error in Eq. (44) for the strap plate problem. The meshes in Fig. 24c are used. The SU4 element shows better computational efficiency among the tested elements. In other words, the SU4 element requires less computation time compared to other elements when obtaining similar solution accuracy.

Fig. 25
figure 25

Performance of the SUFE: a computational efficiency and b convergence curves for the strap plate problem; Computation time is measured in seconds, and the number of element layers (\(N\)) is displayed near each mark

Full size image

The convergence curve is shown in Fig. 25b. The proposed element provides accurate solutions with few DOFs. Figure 26 displays the effective stress (von Mises stress) distributions obtained using the proposed SU4 element.

Fig. 26
figure 26

Effective stress distributions obtained using the SU4 element for the strap plate problem

Full size image

7 Conclusions

In this paper, an SUFE and an iterative solution procedure were proposed. This new concept affords the capability to improve the accuracy of FE solutions without mesh refinement. The assumed modal strain was employed based on the mode-based FE formulation, leading to the SUFE. To calculate the element stiffness matrix for a certain element deformation, the search for the optimal bending directions is required. To reduce the search time, deep learning is employed. The performance of SUFE was investigated through various numerical examples. The accuracy and computational efficiency of the proposed SU4 element were compared with those of other existing four-node elements. The SUFE produced highly accurate solutions with few iterations, especially when the meshes are coarse and distorted.

This study only examined the four-node finite element. However, the proposed concept can be extended to various types of finite elements, including solid, plate, and shell elements [31,32,33,34,35,36]. The conduct of research on such elements in the future will certainly be advantageous.