# High-Precision Elastoplastic Four-Node Quadrilateral Shell Element
## 1 Introduction  
In engineering, various shell elements have been developed [1; 2; 3; 4; 5; 6], including the folded shell element, which is created by combining plane stress and plate bending elements, as well as the degenerate shell element, derived from the shell element hypothesis based on three-dimensional entity theory. The quadrilateral shell element, recognized for its simple structure and versatility, is widely utilized across different industries. However, the development of quadrilateral shell elements with enhanced accuracy presents a significant challenge, particularly in terms of addressing shear and membrane locking effects [7]. These issues can greatly affect the accuracy of solutions for plane stress and transverse shear problems, leading to substantial discrepancies between the calculated and actual results [8; 9]. The traditional quadrilateral shell element is limited by its inability to permit rotation around the normal vector of the neutral plane, which complicates the solution of plane stress problems and exacerbates the issue of film self-locking. The phenomenon of membrane self-locking has been extensively investigated in three-node curved beams and nine-node shell elements [10; 11; 12], leading to the proposal of various methods aimed at mitigating this issue. An improved four-node quadrilateral shell element [13; 14; 15] has been developed that incorporates a degree of freedom for rotation around the normal vector of the neutral plane. Although the Belytschko and Leviathan [16] units can function within a general six-DOF framework, the innovative projection scheme employed means that the sixth rotational degree of freedom does not represent a true degree of freedom. While certain specialized techniques can be used to address this limitation, the optimal solution remains a unit that allows for rotation around the normal. The methods developed by [PERSON] [17], [PERSON] [18], [PERSON] [19], and [PERSON] [20] have been employed to introduce quadratic functions for the interpolation of the element edges. This approach aims to eliminate the translational displacement of the edge midpoint and provide adjacent angular nodes with rotational degrees of freedom around the normal. Consequently, the four-node quadrilateral shell element mitigates the effects of membrane self-locking. However, the introduction of rotational degrees of freedom around the normal results in a redundant zero-energy mode, which does not possess physical significance. [PERSON] [21] addresses the zero-energy mode by incorporating additional stiffness through a penalty function, effectively simulating the strain induced by the rotational degrees of freedom around the normal. [PERSON] [22; 23] proposed four-node and eight-node shell elements (MITC4 and MITC8) utilizing the mixed tensor interpolation strain method [24; 25]. When calculating the transverse shear strain rate, the covariant component is expressed through tensor interpolation in the natural coordinate system and subsequently transformed to the element coordinate system, thereby eliminating shear self-locking. In the context of elastoplastic problems, the stress is updated according to the J2 flow theory [26; 27; 28; 29; 30]. The second invariant function of the skew stress tensor is employed to assess whether the yield condition has been met, thereby indicating the transition to the plastic stage. The cumulative equivalent plastic strain is derived by calculating the increment in the equivalent plastic strain using the radial return method, while simultaneously updating the back stress and yield stress. Additionally, the tentative deviatoric stress is proportionally reduced to return to the yield surface. It is important to note that the relationship between the yield stress and plasticity is expected to be linear (i.e., linear strengthening). By introducing a specific parameter, the unified representation of ideal elastic-plastic [31; 32], isotropichardening [33; 34], kinematic hardening [35; 36], and mixed hardening models [37] can be achieved.