A robust triangular membrane element

To analyze the plane problem with irregular mesh and complicated geometry, it is helpful to utilize the triangular element. In this study, several optimization criteria will be elaborated. By utilizing these provisions and satisfying the equilibrium conditions, a novel triangular element, named SST, is developed. To demonstrate the high accuracy and efficiency of the new element, a variety of structures will be solved. The findings will prove that the presented element has a low sensitivity to the geometric distortion. Moreover, the parasitic shear error will not arise when this element is employed. In addition to these, the proposed element is rotational invariant. Comparison studies will reveal that the SST element is more robust than the other well-known triangular ones.

plates is quite complex, and requires innovations. The large number of the free parameters, symbolic processing and optimizing the entries of the matrix are among the difficulties, which are encountered in the mentioned process.
Two investigators, named Bergan and Nygard, developed the free formulation (Bergan and Nygard, 1984). In the related area, John Dow et al. proposed another scheme named strain gradient notation (Dow, 1999). Within this procedure, the displacement interpolation field is expressed in terms of strain states. In this way, the roots of many errors in the finite-element scheme can be explored and eliminated. It should be mentioned that strain gradient notation is an explicit and simple representation of the free formulation approach. It is possible to utilize this approach and write displacements and also strains in terms of strain states. Furthermore, some optimal criteria will generally be included to enhance the outcomes of the element. By using these optimality criteria the resultant element is rotational invariance and insensitive to geometry distortion. Moreover, the parasitic shear error is eliminated, and the equilibrium equations are satisfied in this formulation.
In this paper, strain states formulation is presented, and a new triangular element, named SST, is suggested. It is worth emphasizing that the recommended element is useful to analyze the plane problems with irregular mesh and complicated geometry. In fact, the triangular elements are able to model the geometry of the structure more properly. On the other hand, numerical experiences have shown that the low-order triangular elements are less accurate than the quadrilateral ones (Eom et al., 2009). To solve irregular meshes, various means of constructing triangular elements have been proposed by the researchers (Felippa, 2003a;Choo et al., 2006;Huang et al., 2010;Pan and Wheel, 2011;Caylak and Mahnken, 2011;). In the following sections; several numerical tests are performed to demonstrate the capabilities of the suggested element. The obtained results indicate good elemental performances, such as high accuracy, low sensitivity to geometry distortion, rotational invariance and rapid convergence.

STRAIN GRADIENT NOTATION FORMULATION
In strain gradient notation, Taylor's expansion of the strain field is utilized about the coordinates' origin (Dow, 1999). In the two-dimensional state, the strain field consists of three strain function as follows: (1) In these equations,  as a subscript denotes the strain gradient value at the origin. These values are called strain states. The magnitude of axial strain x e at the origin is denoted by ( ) (2) In the common finite element formulation, the displacements are interpolated using a polynomial field function similar to the below form: At the origin, strains and rotation will have the following values: The rigid body rotation is equal to ( ) r r  . The unknown parameters are written in terms of strains as below: The coefficients of the quadratic terms of the field interpolation function can be determined by a similar approach. For this purpose, the next derivatives of the strains are utilized: By replacing the coordinates of the origin and solving the resulting set of equations; the coefficients of the displacement interpolation functions are determined, as follows: The remaining coefficients of the displacement interpolation field are obtained in a similar manner. Let u  and v  be the rigid body translations along the x and y axes, based on this assumption; the element's field functions will have the below form:

OPTIMALITY CONDITIONS
In order to attain the efficient elements, the optimization condition can be inserted into the formulation. In this approach, the roots of many errors may be recognized. Several optimality constraints are included in the authors' formulation, which are discussed in the coming sections.

Satisfying the equilibrium equations
To solve the plane problems in the plane-stress or plane-strain state, the equilibrium relations have to be satisfied. For a homogeneous elastic continuum, these equations are written in the following shape: In these equalities, the functions ( , ) x F x y and ( , ) y F x y are the force field functions along the x and y directions in a rectangular coordinate system, respectively. It is obvious that in the planar problems, gradient of the force field, in the perpendicular direction (z) to plane of the element, is equal to zero. To include the stress fields, they are defined in the rectangular coordinate system. Utilizing the stress-strain relations, Eq. (10) can be rewritten in terms of the strain fields. Based on the Hooke's law for a homogeneous elastic condition, the following relations are valid for a plane problem: In these formulas, l for the plane-strain and plane-stress cases will be equal to respectively. The shear modulus, Poisson's ratio, and elastic modulus are denoted by G, n and E, correspondingly. The next equilibrium equations will be established by inserting Eq. (11) into Eq.

Planar pure bending test
This examination has been utilized for finding the optimal flexural template by Felippa (Felippa, 2003a(Felippa, , 2003b(Felippa, , 2006. In this work, Euler-Bernoulli beam has been applied to assess response of the template to the in-plane bending along the x and y axes. The energy ratios (r) are measured in this experiment. Two parameters, named x r and y r , denote the bending energy ratios of the beam's rectangular part in x and y directions, respectively. Fig. (1) shows a simple structure to be utilized in studying the bending behavior of a beam in the x direction. This structure has a cross section equal to b h , with a moment of x M at the free ends. Ignoring the complicated behavior of the two ends, the majority of the beam will undergo a constant bending moment of ( ) , which will produce a stress field equal to / It should be added that these are exact results obtained from the elastic theory of the beams.
Based on the aforementioned discussions, the internal elastic energy, which is stored in the a b section of this beam, has the following value: ture in the xy plane. If this section is composed of only one rectangular element, the potential energy can be written in the coming form: Figure 2: In-plane pure bending test along the y axes.
Where, bx D and by D are the nodal displacement vectors corresponding to relevant stress fields resulted from x M and y M moments, respectively. The stiffness matrix of the rectangular element is denoted by K . It should be mentioned that this test can also be used for triangular elements. In this case, the potential energy is equal to the sum of the potential energies stored in the two triangular elements, which is composed of section a b . The displacements bx D and by D are essential to evaluate this potential energy. For this purpose, stress fields are obtained for both cases of in-plane bending, which are shown in Figs. (1) and (2). These stress functions have the following relations: With stresses functions at hand, the strain field can be attained. Afterwards, the displacement fields are obtained by integration. The flexural energy ratios for the x and y directions are calculated in the coming shapes: Based on these parameters, an element will be capable of correctly represent the bending in any arbitrary directions, when 1 x r = and 1 y r = ( 1 r = ). It should be added that an element will be over stiff or over flexible while 1 r  or 1 r  . In the case of 1 r = , any aspect ratio of the element will have optimum behavior in bending. Furthermore, if (r) increases as the aspect ratio rises, the locking shear will appear in the element.
At this stage, this test will be investigated based on the strain states. In planar bending in x direction, the real stress field varies linearly along the y direction. According to the Hook's law, this strain state has the coming shape: In the last relationships, 1 b , 2 b , 3 b and 4 b are all constant coefficients. Equations (17) indicates that only ( ) are required for achieving the actual solution for an in-plane flexural problem. A necessary condition for convergence requires that the element assumed strain field should always include constant strain states and the rigid motions' ones. It is worth emphasizing that the planar bending test based on the strain states does not impose any limitations on the geometric shape of the element. In other words, not only it can be utilized for triangular and rectangular elements, but other elements as well.

Rotational invariance
In the suitable elements, properties of the element will not alter with rotation of the coordinate system. These types of elements fall into the category of rotational invariance ones. It is obvious that elements can be configured with different rotational orientation in the mesh of a structure. Therefore, the element should be rotational invariance. To have this important property, terms of the strain field with a complete order should be selected (Dow, 1999;Rezaiee-Pajand and Yaghoobi, 2012, 2013, 2014. For instance, the rotational mapping of a constant strain state has the next appearances: According to these relationships, only if the formulation takes into account all three cases of the constant strain states, then an element is capable of representing constant strains with respect to Latin American Journal of Solids and Structures 11 (2014) 2648-2671 any system of the coordinates. It should be stated that invariance constraint has been used by other researchers in the process of developing fine element templates (Felippa and Militello, 1999;Felippa, 2000).

Lack of the parasitic shear error
The existence of axial strain states in the shear strain interpolation function causes the parasitic shear effect. This bad property will lead to stiffen of the element (Dow, 1999). To grow the knowledge about parasitic shear error, formulation of the four-node rectangular element is examined in the coming lines. This element employs an incomplete polynomial, which contains neither 2 x nor 2 y terms. The strain gradient notation related to these terms can be represented in the following form: Based on the strain-deformation relations, which were given by Equations (2), the interpolation functions of the element can be written in the coming shapes: To find the advantages and deficiencies of this model, the Taylor's expansion of the polynomial should be considered. From the shear strain series, it is observed that the shear strains are independent of the axial strain. Two strain states, , improperly appear in the shear strain interpolation function. As a result, if the element undergoes flexural deformations, the strain states will be nonzero and will incorrectly represent a portion of the shear strain.
It is very important to know, if any element is formulated by using the strain gradient notation, the parasitic shear error can be easily eliminated by excluding the incorrect strain states from the shear strain polynomial. In fact, the parasitic shear error decreases as the finer meshes are utilized. As a result, even coarser meshes can produce good responses if elements are exempted from this error. Furthermore, utilizing the complete interpolation functions will prevent the appearance of this error. Including these important issues in the new formulation, the following strain states will be used for the model: The coming strain field is also utilized: It should be added that x-y coordinate is a generalized one, and can be set up in any optional places. Based on formula (23), Eq. (12) can be written as below: According to these relationships, the body forces are assumed to be zero. In the formulas (24) can be expressed as a function of other strain states. By using these equations, the number of the unknowns will be decreased to 10. Furthermore, the vector of the strain states has the next shape: Using this notion, the strain interpolation field can be represented in the below matrix form: In addition, the displacement interpolation field has the following matrix shape: Introducing nodal coordinates into relation (28) In the last relationship, D denotes the displacement vector. The q G matrix is obtained by inserting the nodal coordinates into q N . It is obvious that q G provides a relation between the strain state vector and the displacement one. The boundary conditions should be applied to the formulation. To have displacements and strains' interpolation fields in terms of nodal values, these functions can be expressed in the succeeding form: The potential energy function, which is in terms of the elastic matrix ( E ), should be minimized in the coming shape: The last formula leads to the structural stiffness matrix. This matrix, which is denoted by K , has the below form:

CREATION OF A SPECIAL ELEMENT
The effectiveness of the formulation is verified through the creation of a special element. For this purpose, the aforementioned model will be applied to a triangular element. As it was previously mentioned, ten nodal unknowns should be identified to formulate the element. The geometry of the proposed element is shown in Fig. (3). It is obvious that any adjacent two elements forms generalized four-sided elements. This is shown in Fig. (3-2). First, a mesh with quadrilateral elements is defined. By using diameter line, each four-sided element is divided into two triangular elements, as it is shown in Figure (3.2). All of these diameter lines have six degrees of freedom. In each element, the nodes located on the middle of diameter lines contain ninth and tenth degrees of freedom.
To calculate q G , the displacement in the direction perpendicular to the side of the element is expressed in terms of u and v. This process is carried out for the side of ij in Fig. (4).

Short cantilever beam under shear force
The first example is a short homogeneous cantilever beam, which is subjected to a parabolic shear force at its free end. In order to model the fixity at the other end of the beam, the nodal displacements on that side are set to zero. The exact deflection of the free end is equal to 0.35601 (Felippa, 2003a). As illustrated in Fig. (5), the beam has a length of 48, height of 12 and width of 1. The elastic modulus is equal to 30000, and the Poisson's ratio is 0.25. All the related values of this example are dimensionless. The total shear load acting on the beam is equal to 40. As it is shown in Fig. (5), the test was performed with regular meshes. These meshes are made of each rectangular unit divided into two-triangular elements. The x ý N N mesh indicates the number of subdivisions in the x-and y-directions. For example, Figs. (5-2), (5-3) and (5-4) show 2×2, 4×2 and 8×2 meshes, respectively. The outcomes of using the SST element and the results obtained by utilizing other triangular elements are given in Tables (1), (2) and (3). The effectiveness of the authors' formulation over the other elements is then examined using these tables. Deflection of the short cantilever beam under end shear for aspect ratio 4:1, 2:1 and 1:1 are given in Tables (1), (2) and (3) Table 3: Normalized deflection of the short cantilever beam for aspect ratio 1:1 shown in Fig. (5-4).
The SST element demonstrates insensitivity to the dimensions, and yields a uniform and rapid convergence towards the exact answer. Comparing the results obtained by using SST and the other elements in Tables (1), (2) and (3) indicates that the new element has its superiority for all cases of the meshing.

Slender cantilever beam under moment
The second benchmark problem is a homogeneous slender cantilever beam. This structure will be subjected to a bending moment at the free end. The nodes at the fixed end will have zero displacements. Based on the beam theory, the deflection of the free end of a beam with length L, moment of inertia z I and elastic modulus E , under a bending moment M will be equal to 2 / (2 ) z ML EI . In this equation, 3 / 12 z I hb = . The elastic modulus E and the Poisson's ratio ν are set to 768 and 0.25, respectively. As it is illustrated in Fig. (6), the length of the beam is 32, the height is 2, and the width is equal to 1 h = . The bending moment applied to this beam is equal to 100 M = . The exact deflection of the beam's free end is equal to 100 (Felippa, 2003a). Regular meshes ranging from 2×2 to 32×2 are used, each rectangle mesh unit consisting of four half-thickness overlaid triangles. The element aspect ratios vary from 1:1 to 16:1. The deflection of the homogeneous slender cantilever beam under a moment at its free end is given in Table (4). Table (4) indicates that SST and OPT elements are more efficient than the other ones. The error of SST is less than one percent for all cases of the meshing.

Thin cantilever beam
As it is shown in Fig. (7), the cantilever beam, which is under in-plane shear, has the length of 100, width of 1 and thickness of 1. The elasticity modulus and the Poisson ratio of this structure are equal to 1000000 and 0.3, respectively. It is important to note that this beam with a force of 1 at its free end is a severe test. This test was performed with regular meshes. These meshes for triangular elements are made of each rectangular unit divided into two-triangular elements. The y x N N mesh shows the number of parts in the y-and x-directions. The 1×100 mesh is utilized for analyzing the cantilever beam. This mesh demonstrates that the number of the subdivisions in the x and y directions are equal to 100 and 1, respectively. In another solution, the 2×100 mesh is used. Since the triangular elements have not been used by other researchers, for the sake of comparison, the responses of famous quadrilateral elements are taken to benefit. According to the reference, the exact free end displacements of the beam in the x-and ydirections are equal to 0.03 and 4, respectively (Wisniewski and Turska, 2009). To compare more accurately, the answers of other good elements are given in Table (5). Based on the results, the error of the answers belongs to the SST element is very low.

Elements
Mesh u x ×100 u y  (Paknahad et al., 2007) 3 4 Table 5: The displacement of the free end of the thin cantilever beam under in-plane shear.

MacNeal's beam
As it is shown in Fig. (8), the thin cantilever beam with three different meshes, including rectangular, parallelogram-shaped and trapezoidal meshes is considered. In this problem, it is intended to assess the loss of accuracy produced by utilizing elements shaped as parallelogram or trapezoid. In fact, this test was proposed by MacNeal as a benchmark to evaluate the sensitivity of quadrilateral elements to mesh distortion (MacNeal and Harder, 1985). The elasticity modulus, Poisson's ratio, and the beam's thickness are 10000000, 0.3 and 0.1, respectively. Two loading cases will be used, which are pure bending induced by unit bending moment at the free end and bending caused by unit shear force at the tip. Under these load cases, the exact displacements of the tip are respectively equal to 0.027 and 0.1081. The regular meshes will be first utilized to compare the accuracy of displacement for bending and shear load. This structure will be analyzed by using three rectangular meshes, 1 × 6, 2 × 12 and 4 × 24, with each rectangle mesh unit consisting of two triangles. The 1 × 6 mesh, shown in Fig. (8), is denoted by A. The related results for the proposed element and other researchers' wellknown elements are shown in Table (6). Based on the numerical outcomes, the high accuracies of the authors' element for all cases of the meshing are evidence.   Similarly, MacNeal's is solved by the SST triangular element for meshes A, B and C, which are shown in Fig. (8). It is important to note that a big aspect ratio and distortion in these meshes make the test difficult and reliable for evaluating the efficiency of the formulations. Table (7) demonstrates the results of meshes A, B and C for proposed element. Since in the cases of B and C meshing, the triangular elements have not been used by other researchers, for the sake of comparison, the responses of famous quadrilateral elements are utilized in Table (7). The outcomes for mesh A indicate that SST element performs well when the structure is subjected to the aforementioned load cases. For both meshes, namely B and C, SST has a low sensitivity to distortion. In contrast; other investigators' formulations are generally sensitive to the distortion of parallelogram-shaped and trapezoidal meshes. In fact, the error of response increases extensively under the aforementioned load cases for the trapezoidal mesh.

Cantilever shear wall
This structure is shown in Fig. (9-1). The elastic modulus, and the Poisson's ratio of the shear wall are assumed to be 20000000 kN/m2 and 0.2, correspondingly. The loading parameters q and p are equal to 500 kN and 100 kN, respectively. This shear wall will be analyzed by the new element using the different regular meshes of Fig. (9-2). Once again; these meshes are made of each rectangular unit divided into two-triangular elements. To provide a basis for comparison of the robustness and accuracy of the authors' element, an eight-node element is used to model the structure, as well.
The top-end lateral drift of this structure will be evaluated by taking advantage of the suggested and the eight-node isoparametric element with various meshing cases. Furthermore, the OPT* element will also be used for the comparison. The results of mentioned element for the different meshing cases are available (Paknahad et al., 2007). Fig. (10) demonstrates these outcomes and the ones obtained from the suggested element, and also the eight-node isoparametric element. It is important to note that the OPT* element is developed for analyzing shear walls. In other words, the OPT* element is a specialized and accurate element for shear wall analysis. The efficiency of the eightnode isoparametric element is another well-known fact. The advantage of the SST element over the OPT* element is clearly revealed in Fig. (10). The advantage of the SST element over the OPT* element is clearly revealed in Fig. (10).

Coupled shear wall
In this section, the new element will be used to analyze a coupled shear wall, which is illustrated in Fig. (11-1). The elastic modulus, and the Poisson's ratio are equal to 20000000 kN/m 2 and 0.25, respectively. The width of this shear wall is 0.4 m, and the magnitude of P is equal to 500 kN. Two different regular meshes, cases of a and b, will be utilized in this analysis. Both meshes are shown in Fig. (11-2). For the sake of comparison; all meshes are made of each rectangular unit divided into two-triangular elements. The lateral drifts of the 2nd, 4th, 6th and 8th stories for two meshing cases are calculated using the authors' element. In this study, the outcome of element OPT* is also utilized (Paknahad et al., 2007). Once again, the eight-node isoparametric element is also implemented. All of the results of suggested element, the OPT* element and the eight-node isoparametric one are given in Table (8). For better comparison, a very fine mesh of the eight-node isoparametric element is used to analyze the shear wall. In fact, the entire shear wall is modeled by using 1010 cm2 elements. This mesh will lead to 26880 strong eight-node isoparametric elements, which is denoted by c.

Thin curved beam
The geometry of the thin curved beam is shown in Fig. (12-1). One end of the structure is fixed, and the other one is under a shear force equivalent to 1. The elasticity modulus, the Poisson ratio and the beam thickness are equal to 10000000, 0.25 and 0.1, correspondingly. The near exact vertical displacement under the load is available and equals to 0.08734 . In this study, the beam will be analyzed by using three rectangular meshes, 6 × 1, 12 × 2 and 24 × 4, with each rectangle unit consisting of two triangles. The 6 × 1 mesh is exposed in Fig. (12-2). Table (

Thick curved beam
In this section, the structure of Fig. (13), will be solved. This thick carved beam is subjected to the shear force loading at the free end. The magnitude of the load is equal to 600. The elasticity modulus, Poisson's ratio and the structural thickness are 1000, 0 and 1, correspondingly. As before, all the unities are consistent. By employing the suggested triangular element, this thick carved beam is meshed according to the Fig. (13-2). It should be stated the other researchers have analyzed this structure by using quadrilateral elements. As it is shown in Fig. (13-1), they have utilized four elements. So far, the triangular elements have not been used by other investigators. Therefore, for the sake of comparison, the responses of famous quadrilateral elements are utilized here.  (Cen et al., 2007). According to the numerical results, the findings prove that the SST and Q6 elements perform much better than the well-known quadrilateral ones.

Cook's skew beam
Other investigators studied the performance of the quadrilateral elements by analyzing this structure (Cook et al., 1989). The geometry and loading of the beam is shown in Fig. (14-1). Shear displacements govern the structural behavior, and the distorted four-sided elements are utilized for meshing. One end of this beam is clamped, and the other end is subjected to a uniformly distributed shear load 1 P = . The elasticity modulus, Poisson's ratio and the structural thickness are 1, 1/3 and 1, respectively. In this study, two types of mesh will be used, as they are shown in Figs. (14-2) and (14-3). The quadrilaterals are generally divided into two triangles in the shortest-diagonal-cut of the Mesh A. To check the distortion insensitivity of the authors' element; the longest-diagonal-cut mesh is also used as Mesh B. Four types of mesh, including 2×2, 4×4, 8×8 and 16×16, are utilized for this analysis. Figs. (14-2) and (14-3) demonstrate the case in which 2×2 mesh, and the SST element are used for the solution of Cook's skew beam. The vertical displacement of point C is inserted in Table  ( 11). For the sake of comparison, the results of other researchers' well-known elements are considered too. It should be added that there is no known analytical solution for this problem, and the response of GT9M8 element, for a 64×64 mesh, is considered as the near exact answer (Long and Xu, 1994 (Long and Xu, 1994) 23.96

CONCLUSIONS
To take benefit from the optimality conditions along with the finite element scheme, free formulations as well as strain gradient notations were employed throughout this study. The number of nodal unknowns was decreased by satisfying the equilibrium equations. This also leads to the robustness of the new element. Since triangular elements are usually suitable for modeling the structures with complex shape, the authors' technique is quite useful. By utilizing a variety of numerical tests, the high accuracy, rapid convergence, low sensitivity to geometry distortion, rotational invariance and lack of parasitic shear error were specified in the SST element. Different benchmark problems were analyzed to confirm the efficiency and the accuracy of the suggested element. In this investigation, the new element and those from other researchers were compared. The findings of the paper confirmed that the authors' formulation was the most advantageous triangular element in solving the plane problems. In all the numerical examples, this element demonstrated satisfying features such as; small error in coarse meshes and rapid convergence to the exact solution.