Analysis of elasto-plastic thin-shell structures using layered plastic modeling and absolute nodal coordinate formulation

A new elasto-plastic thin-shell finite element of the absolute nodal coordinate formulation allowing for large deformation and finite rotation is proposed based on the Kirchhoff–Love theory and layered plastic model. The von Mises yield criterion of plane-stress with linear isotropic hardening is adopted in constitutive description of elasto-plastic material. Owing to the plane-stress constraint, special treatment should be given to the stress update algorithm for plasticity. To accommodate the plasticity formulation, the Gauss-point layered integration is inserted into the thickness of the element to produce the internal force. Then, the Jacobian of internal forces is deduced by deriving the consistent elasto-plastic tangent moduli. To accurately track the load–displacement equilibrium path in the buckling analysis of elasto-plastic thin shells, the arc-length method is used. The dynamics of the thin shells is also studied by using the generalized-alpha algorithm. Finally, several static and dynamic examples are presented to verify the accuracy of the proposed formulation.


Introduction
Thin-shell structures have found important applications in various fields, such as aerospace, mechanical, and civil engineering. When a thin-shell structure is subject to large deformations, any geometrically linear method does not work for the kinematics of the thinshell structure in general. As a consequence, it is necessary to develop geometrically nonlinear approaches to study the thin-shell structure subject to large deformations.
When the internal stress state of a thin-shell structure exceeds its yield stress, plastic strain begins to develop and gives rise to the problem of material nonlinearity. In this case, a nonlinear analysis becomes necessary to make the simulation of the structural behavior more realistic. Two different classical approaches have been available to deal with the material nonlinearity caused by plasticity in thinshell structures. One is the stress resultant approach with elasto-plastic constitutive models formulated directly in stress resultants [1][2][3][4][5][6], and the other is the layered approach, which integrates through the thickness of the shell to obtain the stress resultants [7][8][9][10][11]. The stress resultant approach has a lower computational cost since it does not deal with any throughthickness integration. However, it is very difficult to derive the plastic stress resultant constitutive models with this approach. Even the simplest yield function, such as the von Mises condition, has a complex functional form. Hence, it is difficult to implement such constitutive models in the finite element framework (see the work by Simo [3] and references therein). In addition, the information for the spread of plasticity through the shell thickness may be partially or completely lost in this approach. With the layered approach, the conceptual simplicity is inherited from the three-dimensional or plane-stress plasticity model, which enables the use of the standard return-mapping algorithm. Furthermore, it allows for more accurate prediction of the elastic or plastic state and more accurate calculation of the internal forces. However, this approach has a higher computational cost because of the integration through the thickness [9]. Despite these disadvantages, this study aims to develop the layered approach in order to obtain more accurate results and make it easier to shift from the problem of small strains to that of finite strains.
The absolute nodal coordinate formulation (ANCF), initially proposed for beam elements by Shabana [12], provides an accurate and non-incremental type of finite element to deal with the dynamics of flexible multibody systems subjected to both overall motions and large deformations. In the incremental procedure, such as a corotational procedure, the dynamic equations of an element are first defined in the frame of element coordinates and then transformed to get the dynamic equations in the global frame of reference. These equations are solved for the displacement increments, which are used afterward to update the displacement field of the element in the global frame [13]. To model the rotation and deformation field of a finite element, the gradients of global positions, which replace angles in classical finite element analysis, serve as the nodal coordinates in the ANCF. Because the nodal coordinates of a finite element of ANCF are defined in the global frame of coordinates, the dynamic equations of the system have a constant mass matrix and the absence of centrifugal and Coriolis forces. Therefore, the ANCF is one of the most important advances in the dynamics of flexible multibody systems [14][15][16]. Furthermore, the studies in [17][18][19] have shown that the B-spline geometry used in CAD systems is consistent with the geometric description of ANCF.
To model a thin-shell structure, Mikkola and Shabana [20] proposed a quadrilateral fully parameterized shell element of ANCF. However, Mikkola and Matikainen [21] showed that the shell element suffers from serious shear and thickness locking. To avoid shear and thickness locking, Dmitrochenko and Pogorelov [22] introduced a gradient-deficient thinshell element of ANCF based on the Kirchhoff-Love theory. This thin-shell element of ANCF drops off the gradients along the thickness direction and improves the computational efficiency. Dufva and Shabana [23] showed that this thin-shell element has good convergence. To prevent membrane locking in the gradientdeficient thin-shell elements of ANCF, Sanborn et al. [24] proposed a flat-mapped extension modeling method, which re-describes the displacement field of the shell elements. Liu et al. [25] proposed a thin cylindrical shell element of ANCF so as to model flexible curved parallelogram structures by introducing an angle between two base vectors of the local frame of coordinates on the shell element.
To study elasto-plasticity in flexible multibody systems, Ambrósio and Nikravesh [26] introduced the plastic constitutive law into the floating frame of reference formulation by using the updated Lagrangian formulation and rate-type constitutive equations. Sugiyama and Shabana [27] implemented the von Mises yield criterion with linear isotropic hardening by using the fully parameterized beam elements and shell elements of ANCF. These finite elements, however, surfer from both high computational costs and locking problems. Subsequently, Sugiyama and Shabana derived the consistent elasto-plastic tangent moduli for implicit integration methods [28]. Gerstmayr [29] employed the plasticity models of planestress in the planar four-node finite element of ANCF. Unfortunately, the increment of plastic strain computed was not accurate enough. Gerstmayr and Matikainen [30] used the beam elements of ANCF with 48 degrees of freedom to study the dynamics of an elasto-plastic pendulum so as to improve the accuracy of stress and strain quantities and prevent the locking behavior. Recently, Wang et al. [31] reported the contact dynamics of elasto-plastic spatial thin beams using the gradient-deficient beam elements of ANCF. As such, the studies on the dynamics of elastoplastic thin shells via the gradient-deficient shell elements of ANCF have yet to be performed.
The objective of this work is to extend the gradientdeficient shell element of ANCF to model the elastoplastic behavior of thin shells undergoing both overall motions and large deformations by using a layered plasticity approach. The work is based on the von Mises yield criterion of plane-stress with linear isotropic kinematic hardening based on the additive split of strains into elastic and plastic parts [32,33]. The backward Euler time integration scheme, required for the solution of the elasto-plastic constitutive equations, with the return-mapping algorithm is used. The plasticity incremental procedure is used to solve the equation of yield criterion, while the non-incremental procedure for describing the finite rotation is used to solve the equation of motion. To obtain the plastic multiplier, the yield equation need to be solved by using the local Newton iteration or the method proposed by Simo [34].
The remainder of this paper is organized as follows. In Sect. 2, the formulation of a gradient-deficient thinshell element of ANCF is outlined in view of the Kirchhoff-Love theory. In Sect. 3, the details of the layered plasticity approach and the constitutive models based on the von Mises yield condition are presented. They are then used as the bases to derive the vector of internal forces and the corresponding Jacobian of a shell element of ANCF. In Sect. 4, a computational procedure is presented for nonlinear static analysis, including buckling analysis. The arclength method is used to trace the nonlinear equilibrium path. A generalized-alpha algorithm for solving the dynamic equations is also applied. In Sect. 5, five numerical examples of both statics and dynamics of elasto-plastic thin-shells are presented, and simulations are conducted via the commercial software ABAQUS for comparisons. Finally, the major conclusions are drawn in Sect. 6.

Kinematics of the shell element
This subsection presents the gradient-deficient thinshell element of ANCF developed in previous studies [22,24,25] via the total Lagrangian formulation. Figure 1 shows a scaled view of an arbitrary layer p z in a thin-shell element with the nodes of the shell element defined on the mid-surface p. Based on the Kirchhoff-Love theory, the position vector of an arbitrary point P 0 ðn; g; zÞ on the surface p z reads r z 0 n; g; z ð Þ¼r 0 n; g ð Þþzn 0 n; g ð Þ where ''0'' indicates the initial state, 0 n 1 and 0 g 1 are the canonical parameters of the midsurface of the shell element, and z is the distance between the surface p z and the mid-surface p . In addition, r 0 is the location of the corresponding point Pðn; gÞ on the mid-surface p , n 0 is the unit normal vector of the mid-surface p at this point, and e 0 represents the vector of nodal coordinates of the shell element and yields [22] where ðr 0 Þ i is the position vector of the ith node, oðr 0 Þ i on and oðr 0 Þ i og are the slope vectors of the ith node, (i = 1, 2, 3, 4), respectively, for the shell element. The shape function matrix of the element yields where I 3 is a 3 Â 3 unit matrix. The shape functions S 1 ; :::; S 12 can be found in the research by Dmitrochenko and Pogorelov [22].

Deformation of the shell element
To accurately describe the deformation of the thinshell element, two local frames of coordinates are established at point P 0 ðn; g; zÞ, as shown in Fig. 1. The first frame consists of curved surface coordinates ðg z 0 Þ 1 À ðg z 0 Þ 2 À n 0 as follows: where the superscript z indicates that the curved surface coordinates are with respect to the surface p z . Furthermore, the frame of ðg z 0 Þ 1 À ðg z 0 Þ 2 À n 0 is not defined using unit vectors. Then, the second frame is the local frame of Cartesian coordinates ðe z 0 Þ 1 À ðe z 0 Þ 2 À ðe z 0 Þ 3 determined via the frame of curved surface coordinates as follows: The local frames of coordinate ðg z Þ 1 À ðg z Þ 2 À n and ðe z Þ 1 À ðe z Þ 2 À ðe z Þ 3 are also presented in the current configuration. Based on the Kirchhoff-Love theory, the relation between the infinitesimal arc segments dX z and dx z , defined in the initial and current configurations, respectively, can be expressed as where the matrix F is a 2 Â 2 matrix describing the deformation of the infinitesimal arc segment in the tangential plane of the surface p z . dX z and dY z are the components of dX z in the local frame of Cartesian coordinate ðe z 0 Þ 1 À ðe z 0 Þ 2 À ðe z 0 Þ 3 , and dx z and dy z are the components of dx z in the local frame of Cartesian coordinates ðe z Þ 1 À ðe z Þ 2 À ðe z Þ 3 . To obtain F, it is easy to decompose the infinitesimal arc segment dx z into the frame of coordinates ðg z 0 Þ 1 À ðg z 0 Þ 2 À n 0 and ðe z 0 Þ 1 À ðe z 0 Þ 2 À ðe z 0 Þ 3 , and get the two summations via two pairs of dummy indices where v 1 ¼ n; v 2 ¼ g; Furthermore, the covariant transformation matrix b j i enables one to derive the relation between the base vectors of the two local frame of coordinates ðg z 0 Þ 1 À ðg z 0 Þ 2 À n 0 and ðe z 0 Þ 1 À ðe z 0 Þ 2 À ðe z 0 Þ 3 as Because there is no distinction between the covariant and contravariant base vectors in the frame of Cartesian coordinates, the covariant transformation matrix b j i reads Therefore, the contravariant components of the infinitesimal arc segments dX z in the two frames of coordinates above satisfy where b i k is the contravariant transformation matrix and yields b i k b j i ¼ d j k . According to Eqs. (7), (8), and (10), one derives dn ¼ dn dg As the thickness in the Kirchhoff-Love theory is constant, F yields Taking the rigid rotation of the shell element into account, the matrix of displacement gradients can be obtained by is an orthogonal transformation matrix between the two local Cartesian coordinate frames ðe z 0 Þ 1 À ðe z 0 Þ 2 À ðe z 0 Þ 3 and ðe z Þ 1 À ðe z Þ 2 À ðe z Þ 3 . According to continuum mechanics [35] and the orthogonality relation between vectors e z 0 À Á 1 and e z 0 À Á 2 , one defines the Green-Lagrange strain as where g z To obtain the explicit expression of the Green-Lagrange strain, one can expand Eq. (4) as Finally, from Ref. [25], it is easy to derive the Green-Lagrange strain tensor as follows: where e mid and e j are, respectively, the mid-surface strain and bending strain of the shell element and yield Because of the assumption of plane-stress and Kirchhoff-Love, the 3-D plastic model generally used in flexible multibody systems cannot be used in the gradient-deficient thin-shell elements of ANCF. The detail for the plastic model used in this paper will be presented in Sect. 3.
3 Elasto-plasticity constitutive model for thin-shell elements

Three-dimensional elasto-plastic constitutive model
In this work, the description of the three-dimensional elasto-plastic constitutive model focuses on the case when the strains remain small. In this case, the total Green-Lagrange strain can be split into elastic and plastic parts under the assumption of additive decomposition [36], i.e., The elastic part described by the St. Venant-Kirchhoff constitutive model is as follows: where s is the second Piola-Kirchhoff stress tensor, and C denotes the elastic tangent moduli. If nonlinear elastic constitutive model is considered, the assumption of multiplicative decomposition is required to be used. The von Mises yield criterion is applied in this study, which can be represented as where dev s ð Þ :¼ s À 1 3 s : I ð ÞI is the deviatoric stress and dev s ð Þ k k:¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dev s ð Þ : dev s ð Þ p is its Euclidean norm. a is the effective plastic strain and RðaÞ is the radius of the yield surface.
In the case of full finite strain kinematics, the correct definition of the deviatoric part of the second Piola-Kirchhoff stress tensor is where C ¼ 2e þ I. Under the assumption of small strains, Eq. (25) can be recast as where o e k k ð Þ denotes the terms that tend to be zero. Therefore, when strains are restricted to be small, an approximate definition of the deviatoric stress can be employed. Accordingly, in situations where strains remain infinitesimal but have large deformations and finite rotations, the above von Mises model introduced by Eq. (24) is equivalent to the von Mises model of small strains. This makes all aspects of material identification in the von Mises model of small strains transferrable to the von Mises model of plane-stress.
The associated flow rule and hardening law for the von Mises yield criterion can be expressed as [32,33,37] and a ¼ c where c is the plastic multiplier that can be solved from the Kuhn-Tucker conditions as follows:

Plane-stress von Mises theory with linear isotropic hardening
The gradient-deficient thin-shell elements of ANCF are based on the Kirchhoff-Love theory and the assumption of plane-stress; thus, the 3-D plastic model generally used in flexible multibody systems should be replace by the model of plane-stress in this paper. This subsection presents the von Mises theory under the plane-stress condition in detail.
In accordance with the Voigt notation and planestress condition, the stress and strain tensors can be defined as The mapping, which connects the stress tensor s and its deviator dev(s), plays a crucial role in the von Mises theory of plane-stress as follows: where the matrix P can be expressed as In terms of the model described in Sect. 3.1, the von Mises model of plane-stress can be written as [38] s ¼ C e À e p ð Þ; where the radius of the yield surface RðaÞ is defined by the linear isotropic hardening, i.e., RðaÞ ¼ r Y þ Ka in this study. r Y and K denote the initial yield stress and plastic modulus, respectively. Matrix P has the form As stated by Kleiber [39], matrix P differs from the previous matrix P. Matrix C denotes the two-dimensional elastic tangent moduli of the St. Venant-Kirchhoff constitutive model in the form of Voigt notation as follows: The linear isotropic and kinematic hardening law can also be used. In this case, the center of the yield surface experiences a motion in the direction of the plastic flow. The yield criterion can be expressed as whereb is the back stress as follows: where H 0 a ð Þ is called kinematic hardening modulus. Nonlinear isotropic and kinematic hardening can also be accommodated, but the linearization of the yield function and the internal forces will become more complex. Furthermore, the implementation of the return-mapping algorithm is more difficult. In this paper, only the law of linear isotropic hardening is used.

Efficient formulations of internal forces and their Jacobians
Before calculating the vector of internal forces and the corresponding Jacobian to reduce the storage of the zeroes and thereby enhance computational efficiency, the position vector of an arbitrary point on the midsurface of the shell element can be recast as where e is a generalized coordinate matrix of size 3 Â 12 as follows: and where S 1 , …, S 12 are the nonzero entries of the shape functions. According to Eq. (19), the components of the mid-surface strain tensor e mid can be expressed as The subscripts in Eq. (41) are a = 1, 2; n = 1, 2; g = 1, 2;, r = 1, 2; l = 1, 2, 3; s = 1, 2, 3; t = 1, …, 12; and k = 1, …, 12. Similarly, according to Eqs. (20) and (21), the components of the bending strain tensor e j can be written as where rn H tgr are the third-order constant matrices. The subscripts in Eq. (42) are a = 1, 2; n = 1, 2; g = 1, 2; r = 1, 2; l = 1, 2, 3; s = 1, 2, 3; and t = 1, …, 12.
The vector of internal forces of the shell element can be expressed as [35] : where e denotes the vector of nodal coordinates of the shell element, and D¼ g 0 ð Þ 1 Â g 0 ð Þ 2 . The gradientdeficient thin-shell element of ANCF with a linear elastic constitutive law enables one to simplify the integration through the thickness by using analytical integration. Thus, the vector of the internal forces can be written as where S is the area of the regular master element. However, for the elasto-plastic constitutive law, owing to the different stress states in an arbitrary layer of the thin-shell element, such as the layer close to the midsurface, elasticity may be retained, but plastic deformation may appear on the layer far from the midsurface. This will require a numerical integration through the thickness to compute the internal forces. Therefore, Eq. (43) can be recast as According to Eq. (41), the term oe mid oe can be directly expressed in terms of its components as follows: where the subscript i ¼ 3 Â k À 1 ð Þþl: Similarly, from Eq. (21), after some complex mathematical manipulations, oj ab oe can be written as where S ;ab (a; b ¼ 1; 2) denotes the partial derivative of S with respect to n and g (e.g., S ;12 ¼ oS onog ), and t ¼ r n Â r g . The subscripts in Eq. (47) are f = 1, 2, 3; s = 1, …, 12; and a = 1, 2, 3. Then, the expression of oe j oe can be determined from Eqs. (42) and (47), respectively: Finally, the internal forces are obtained.
In terms of Eq. (45), the Jacobian of the internal forces can be written as where C ep denotes the consistent elasto-plastic tangent moduli, which will be introduced in Subsection 3.5. As can be seen from Eq. (48) where the subscript j ¼ 3 Â t À 1 ð Þþl: According to Eq. (47), o 2 j o 2 e can be expressed as in which and where 2 mnk is the component of the Eddington tensor in Cartesian coordinates, and L ¼ S ;n S ;g denotes the dyadic product of S ;n and S ;g . The subscripts in Eqs. (51), (52)  3.4 Return-mapping algorithm for plane-stress elasto-plasticity formulation From subsection 3.3, to obtain the internal forces of the thin-shell elements, the plastic strain and second Piola-Kirchhoff stress tensor must be updated at each time step. The standard return-mapping algorithm can be used to update the three-dimensional plastic models, but this violates the plane-stress condition for the gradient-deficient thin-shell element of ANCF. Hence, the algorithm is no longer applicable. For the plane-stress case, starting from Eq. (33), a backward-Euler difference scheme yields the following approximation of the plastic return-mapping: a nþ1 ¼ a n þ Dc where f nþ1 is defined as Then, the stress-strain relationship is expressed as where s trial nþ1 denotes the trial stress. In case of elastic state, s trial nþ1 yields From Eqs. (54) and (57), the relationship between s nþ1 and s trial nþ1 yields To facilitate the calculations, Eq. (59) can be recast as where NðDcÞ plays the role of a modified elastic tangent matrix defined as The stress updating in Eq. (60) is dependent on the plastic multiplier Dc, which is determined by solving the Kuhn-Tucker consistency condition at time t n?1 In the case of isotropic elasticity, Eq. (62) has a remarkably simple form to avoid complex linearization. The elastic tangent moduli C and the projection matrix P have the same characteristic subspaces. Hence, they can be easily diagonalized as follows: where the orthogonal matrix Q À1 ¼ Q T is expressed as The diagonal matrices K P and K C are given by the expressions From Eq. (66), the expression is obtained as where CðDcÞ is a diagonal matrix given by With this notation, Eq. (62) can now be expressed as With the simple expression for the Kuhn-Tucker consistency condition, complex linearization can be avoided. The differentiation of Eq. (73) yields Then, Newton iteration is used to update the consistency parameter Dc as follows: Remark In the study by Gerstmayr [29], an approximation of the flow rule was used to simplify the calculation. The updated plastic strain can be expressed as Then, Eq. (59) becomes It can be seen that the relationship between s nþ1 and s trial nþ1 is linear. Consequently, the Kuhn-Tucker consistency condition at time t n?1 can be written as a quadratic equation for Dc as follows: This method reduces the nonlinearity of the Kuhn-Tucker consistency equation to simplify the calculation, but the computational accuracy decreases. This will be verified later in the numerical examples.

Consistent elasto-plastic tangent moduli
Both the static and dynamic analyses of the thin-shell elements require linearization to complete the iterative procedure. To obtain the Jacobian of the vector of internal forces, the consistent elasto-plastic tangent moduli must first be solved. The differentiation of Eqs. (55) and (57) leads to where f nþ1 is defined by Eq. (56), and df nþ1 is computed by differentiating this expression as follows: Equation (81) can be recast from Eq. (61) as Differentiation of the consistency condition in Eq. (62) at t nþ1 using Eq. (82) yields the following expression Substituting Eq. (83) in Eq. (84), dDc can be solved as Finally, according to Eqs. (83), (85) and (86), the consistent elasto-plastic tangent matrix can be expressed as where x nþ1 is defined as 4 Computation strategy

Computation approach for nonlinear buckling analysis
The equilibrium equation that governs the nonlinear static or quasi-static problem can be expressed as where q is the vector of generalized coordinates of the system; F int and F ext are the vector of internal forces and vector of generalized external forces, respectively, k is the loading parameter, and f ext is the vector of external forces. The above equilibrium equation is usually solved by using the Newton-Raphson method. Therefore, the governing equilibrium equation can be further recast in an iterative form as where i denotes the ith step; and K t and G represent the tangent stiffness matrix and the vector of out-ofbalance forces, respectively. Normally, the loaddisplacement equilibrium path can be obtained by solving Eq. (89). However, the tangent stiffness matrix K t becomes singular and leads to the divergence of the Newton-Raphson solution when buckling phenomena occur. Analytical techniques that allow the limit points to be passed are required to trace the load-displacement equilibrium path. One of the most widely used approaches is the arc-length method. The main idea of the arc-length method is to increase the dimensions of the equilibrium equation by taking the loading parameter k an independent variable. The load increment is governed by an extra constraint equation with the general form [40] where Dq i and Dk i denote the vector of incremental displacements and an incremental load factor from the last converged point, respectively, A 0 is a scalar parameter that governs the relative contributions of displacement and load increments; and Dl represents the arc-length increment for the current load step. According to the widely used cylindrical arc-length method proposed by Crisfield [40], parameter A 0 is set as zero. The iterative process of the arc-length method is shown in Fig. 2.
There are difficulties when solving the quadratic constraint equation described by Eq. (91). The first is the choice of roots. Crisfield [40] proposed a method for evaluating the scalar product (Dq T iþ1 Dq i ) for each root and then choosing the root with the largest product, i.e., the smallest angle between Dq i and Dq iþ1 . Second, complex roots can sometimes occur. Lam and Morley [41] proposed a pseudo-line-search technique to solve this problem, which can be formulated without involving the concept of potential energy.
For all load steps, the initial estimate of the load increment is expressed as where d t denotes the tangential displacement, i.e., d t ¼ K t ð Þ À1 f ext . The success of the path-following technique depends on choosing the appropriate sign for the iterative load factor. In this study, the criterion proposed by Feng et al. [42] is used as follows: where Dq pr is the previously converged displacement increment. This criterion can accurately predict the path direction even in extreme situations, such as snap-backs, where some widely adopted methods may fail.

Computation procedure for dynamics
Based on the ANCF, the final dynamic equations for a system of thin shells can be expressed in a compact form as a set of differential algebraic equations (DAEs) with a constant mass matrix as follows [43]: where M is the constant mass matrix of the system; Fig. 2 Iterative process of the arc-length method F int q ð Þ is the vector of the internal forces; U q; t ð Þ denotes the vector of the constraint functions; U q is the derivative matrix of U with respect to the vector of generalized coordinates q; k represents the vector of Lagrange multipliers, and F ext q; _ q ð Þ is the vector of generalized external forces applied to the system, which can be obtained using the principle of virtual work. Several integration procedures [44][45][46] have been proposed to solve Eq. (94). In this study, the generalized-alpha algorithm is used to achieve the optimal combination of accuracy in the low-frequency range and numerical damping in the high-frequency range [47,48]. According to the generalized-alpha algorithm, the following set of linear algebraic equations need to be solved via Newton-Raphson iteration at each integration step: where Here, K int q ð Þ and K ext q; _ q ð Þ are the Jacobians of the vectors of internal forces and external force, respectively.ĉ q ð Þ andb q ð Þ are the algorithm parameters for the generalized-alpha algorithm, both of which are determined by the parameter q 2 0; 1 ½ of the algorithm.

Statics of a cantilever beam under tip loads
The example of concern is the linear isotropic elastoplastic analysis of a cantilever beam under two tip loads shown in Fig. 3, which was studied by Dvorkin et al. [49] and by Eberlein and Wriggers [50]. The major objective of this subsection is to make a comparison of the method proposed with previous studies. As in those studies, the length l, width w, and thickness h of the cantilever beam are 10 m, 1 m, and 0.1 m, respectively. To improve the convergence efficiency, both points A and B are subject to loading displacement u F , which increases from zero to the maximal value of 4 m. The elastic parameters are Young's modulus E = 1.2 9 10 7 N/m 2 and Poisson's ratio m ¼ 0:3, and the plastic parameters are the initial yield stress r Y = 2.4 9 10 4 N/m 2 and the plastic modulus of the von Mises criterion K = 1.2 9 10 4 N/ m 2 , respectively.
To check the deformation process in detail, this study uses 20 equidistant time increments and takes the error tolerance as 1 9 10 -6 . The study by Eberlein and Wriggers [50] showed that 20 finite elements led to a convergent solution. To achieve better results, 30 shell elements of ANCF are used in this example. Figure 4 presents the original and final deformed configurations of the cantilever beam, as well as the effective plastic strain. The figure indicates that the plastic deformation occurs mainly in the constraint region. Figure 5 shows the load-deflection curves  with three, five, seven, and nine Gauss points through the thickness and 3 9 3 Gauss points over the shell element. The curves look fairly invariant when the number of Gauss points increases from seven to nine. This implies that seven Gauss points through the thickness are sufficient to achieve a good approximation. Accordingly, a large number of Gauss points through the thickness will be necessary when the dominant deformation is bending. To validate the method proposed, Fig. 5 provides a comparison of load-deflection results between the method proposed and the study by Eberlein and Wriggers [50]. The comparison shows good agreement up to a displacement 2.5 m. For larger displacements, the increasing difference between the results of Eberlein and Wriggers [50] and the present study is due to the number of Gauss points through the thickness as Eberlein and Wriggers [50] used five Gauss points. To verify this observation, Fig. 5 also presents the numerical solution by using the S4R elements of ABAQUS with seven Gauss points through the thickness, exhibiting good agreement. Figure 6 shows a comparison between the linear hardening model and the perfect plasticity model for the same mesh and number of Gauss points. The figure indicates that the linear hardening model requires a larger force than the perfect plasticity model for the same deflection. The reason of the phenomenon is that the stress of the linear hardening model is larger than the perfect plasticity model for the same strain.

Statics of a pinched cylinder with end diaphragms
The example is a thin cylindrical shell, as shown in Fig. 7, bounded by two rigid diaphragms and subjected to two diametrically opposite loads P at the midpoint of the top and bottom surfaces. This is a classical and challenging benchmark to validate elasto-plastic shell formulations, as studied in previous works [3, 7-9, 51] using different approaches. As in previous studies, the radius, length, and thickness of the cylinder are R = 300 mm, L = 600 mm, and h = 3 mm, respectively. The elastic parameters are Young's modulus E = 3000 N/mm 2 and Poisson's ratio m ¼ 0:3, while the plastic parameters are the yield stress r Y = 24.3 N/mm 2 and the plastic modulus of the  The boundary condition of the rigid diaphragms, enforced at the two edges of the cylinder, is that the displacements of the edges in both X-direction and Zdirection are zeros. Because of the geometrical symmetry of the example, one only needs to study an octant of the shell as meshed in Fig. 7. The cylindrical arc-length method is used in this example, and the largest displacement is set to 299 mm. To display the deformation process in detail, the convergence tolerance is set as 1 9 10 -6 for the arc-length iteration, but reduced to 1 9 10 -8 in the iteration for incremental plastic multiplier.
In the numerical simulations, the cylinder is meshed via 16 9 16, 32 9 32 and 48 9 48 shell elements of ANCF, respectively, and then computed for five Gauss points through the thickness and 3 9 3 Gauss points over the shell element. Figure 8 shows the deformed configurations and effective plastic strain for different loading stages, where the structure  Figure 9 gives the load-deflection curve at point A of the shell. One can see the appearance of local buckles on the curve computed via a 16 9 16 mesh of shell elements. Increasing the number of shell elements leads to more accurate results. Accordingly, it is necessary to use a refined mesh of shell elements to accurately capture the wrinkles generated during the deformation process. As shown in Fig. 9, the mesh size of 32 9 32 is sufficient for convergent results. For comparison, Fig. 9 presents the result by Areias et al. [51], which the proposed formulation well agrees with. To demonstrate the advantage of the proposed method, Fig. 9 also gives the result by using the method of Gerstmayr [29], which looks almost identical with the results of this study when the deformation is small, but exhibits an obvious difference when the displacement under the point load reaches 200 mm. Therefore, the method proposed in this study works for large deformation problems.

Dynamics of a plane-stress cantilever beam with perfect plasticity
This subsection presents the dynamic analysis of a plane-stress cantilever beam under a concentrated load to validate the effectiveness of the proposed elastoplastic shell element of ANCF in a two-dimensional situation. As shown in Fig. 10, a cantilever beam of length L = 4 m, width w = 0.5 m, and thickness h = 0.01 m is initially at rest and then is suddenly subjected to a constant concentrated load P = 10 kN at the midpoint A on the right tip. The material density is 7800 kg/m 3 . The Young's modulus and Poisson's ratio are set to 2.1 9 10 11 N/m 2 and 0.3, respectively. The von Mises theory of perfect plasticity is used, and the initial yield stress is 1 9 10 8 N/m 2 . The generalized-alpha algorithm in Sect. 4 is used to solve the dynamic equations with the parameter q ¼ 0:8. The time step Dt of the integration is set as 1 9 10 -4 s. The convergence tolerance is 1 9 10 -6 for dynamic iteration, and 1 9 10 -8 when iterating for the incremental plastic multiplier. Because this is a twodimensional problem, zero Gauss points through the thickness and a 3 9 3 Gauss quadrature over the area of the element are chosen. For comparison, the same problem is solved using the commercial software ABAQUS. As shown in Fig. 11a, the convergent results with respect to the vertical displacements at  Fig. 11 Comparison of evolution curves of vertical displacement at point A with different element meshes and the same problem computed by ABAQUS: a elasto-plastic and b elastic point A by using perfect plasticity can be obtained using 24 9 12 mesh of shell elements. Good agreement is obtained between the results of the ANCF and ABAQUS. Figure 11b shows a comparison of the evolution curves of the vertical displacements at point A, which were computed using the elastic model with different element meshes. It can be seen that the elastic model required fewer elements to obtain a convergent result compared with the plastic model because different numbers of elements may cause different stress states at Gauss points. Accordingly, it is important to choose a suitable element size for the plastic model. Owing to the plastic dissipation of the cantilever beam, the displacement computed using the elasto-plastic shell elements of the ANCF is smaller compared with that using the elastic ones, as shown in Fig. 11a and b.

Dynamics of a free-falling plate-shaped pendulum
This subsection presents the dynamic simulation of a free-falling plate-shaped pendulum with elasto-plastic deformation. This problem was studied as a test example for nonlinear shell dynamics by several authors using various finite elements of the ANCF [20,22,23]. However, material nonlinearity, such as elasto-plasticity, has yet to be considered. As shown in Fig. 12, a heavy elasto-plastic plate with one fixed corner undergoes pendulum-like oscillations, which occur not only in large rectilinear displacements but also in large rotations under gravitational forces, and the plate is initially at rest. The length, width, and thickness of the plate are 0.3 m, 0.3 m, and 0.35 mm, respectively. The mass density is set to 1800 kg/m 3 . Because of the assumption of a small strain, a larger Young's modulus E = 2.4 9 10 9 N/m 2 compared with that in [20,22,23] and Poisson' s ratio of m = 0.3 are used. A linear isotropic hardening model with an initial yield stress of r Y = 6 9 10 6 N/m 2 and plastic modulus of K = 1.2 9 10 7 N/ m 2 is used to model the plasticity of the plate.
The gravitational acceleration is set as 9.81 m/s 2 . In the simulation, the parameter q of the generalizedalpha algorithm was chosen to be 0.8. The time step Dt of the integration was set as 1 9 10 -4 s. The convergence tolerance was 1 9 10 -6 for dynamic iteration, and 1 9 10 -8 when iterating for the incremental plastic multiplier. The plate-shaped pendulum is simulated for 0.5 s using a 10 9 10 mesh of shell elements and analyzed for five Gauss points through the thickness and 5 9 5 Gauss points over the area of the element. This is sufficient to obtain convergent results. Figure 13a and b shows the X-, Y-, and Z-displacement of the free tip of the pendulum for the elastic and elasto-plastic constitutive models. A comparison between the evolution curves of elasticity and elastoplasticity shows that the difference between the elastic and elasto-plastic models is not very obvious at approximately 0.13 s because the plastic strain and dissipation are small at the beginning of the simulation. Afterward, the plastic strain and dissipation increase over time in the dynamic simulation. The Z-displacement of the tip of the elasto-plastic material is smaller than that of elastic material. The maximum deflection of the elasto-plastic model is larger than that of the elastic model; therefore, the maximum strain is similar. To verify the validity of the proposed formulation, the same problem is solved using the S4R element of ABAQUS. As shown in Fig. 13a, b, good agreement is obtained between the results of the proposed formulation and ABAQUS.
The evolution curves of the tip Z-velocity are displayed in Fig. 14; these indicate that the elastic and elasto-plastic models are quite different. Owing to the velocity is the higher-order term of displacement, plastic strain has a stronger effect on velocity and the time of obvious difference between the elastic and elasto-plastic models decreased from 0.13 s to 0.07 s. Figure 15 shows the deformed shapes of the plateshaped pendulum as well as the effective plastic strain for different time steps. It can be seen that the plastic strain is large at the fixed end. The final case study focuses on the application of the shell elements of ANCF with a layered plastic model to a rigid-flexible multibody system undergoing large overall motions. As shown in Fig. 16, the system which consists of a rigid cube and a flexible panel arm is initially at rest. The initial configuration of the system is symmetrical with respect to the X-, Y-, and Z-axes. The side length and mass of the rigid cubes are 2 m and 800 kg, respectively. The two constant (a) (b) Fig. 13 Comparison of evolution curves of tip displacements of a plate-shaped pendulum for elasticity and elasto-plasticity: a Zdisplacement and b X-and Y-displacement The parameter q of the generalized-alpha algorithm and the convergence tolerance are set as 0.8 and 1 9 10 -6 , respectively. The time step Dt is 1 9 10 -4 s, and the total simulation time is 7 s. The rigid cube is modeled via the natural coordinate formulation [52]. A mesh of 20 9 4 shell elements of ANCF is sufficient to obtain a convergent result for the panel arms. It is analyzed for five Gauss points through the thickness and 5 9 5 Gauss points over the shell element. The global coordinates of the points constraining the rigid cube and panel arm can be obtained by using the local coordinates of the points in the rigid cube, and then determining the vector of the constraint function.
The evolution curves of the X-, Y-, and Z-displacement at midpoint A of the free edge are shown in Fig. 17 a, b, and Fig. 17 Comparison of evolution curves of displacements at point A with elasticity and elasto-plasticity: a X-displacement, b Ydisplacement, and c Z-displacement same problem is solved using elastic material. It is observed that the differences between the evolution curves of the elastic and elasto-plastic materials are not very obvious before approximately 3 s. Afterward, the two models exhibited different displacement values owing to the increase in the plastic strain and dissipation over time. Figure 18 shows the deformed configurations of the rigid-flexible multibody system as well as the effective plastic strain at some time steps during the simulation. In this case, bending and torsion occur in the panel arms, and the plastic strain first appears at the roots of the panel arm. With the propagation of the stress waves, the oscillation initially at the roots of the panel arm spread to all regions. Therefore, the plastic strain gradually moves toward the free edge of the panel arms. It is important to note that the plastic strain is very small in the early time step and is not shown in Fig. 18. In fact, the plastic strain is originally generated at 1.16 s. Fig. 18 Motion of the rigid-flexible multibody system at some time steps On the basis of the Kirchhoff-Love theory, an elastoplastic thin-shell element of ANCF is proposed to describe elasto-plastic thin-shell systems subject to the large deformations coupled with large overall motions. Both geometrical and material nonlinearities are taken into account in the thin-shell analysis. The layered approach and the von Mises plastic model of plane-stress with isotropic hardening are applied in the study. The analytical formulations of the vector of internal forces and the Jacobian of the elasto-plastic shell elements are derived from those of the elastic shell elements of ANCF. The static or quasi-static problems, including buckling with severe nonlinearities, are solved by using the cylindrical arc-length method and specific techniques. For the dynamic problems of flexible multibody systems, the dynamic equations in the form of DAEs are solved by using the generalized-alpha algorithm. The performance of the proposed elasto-plastic thin-shell elements of the ANCF is validated via five case studies. The first example shows that a large number of Gauss points through the thickness will be necessary when the dominant deformation is bending. From the pinched cylinder, the method using an approximation of the flow rule to update the plastic strain is found inaccurate for large deformations. The three dynamic examples indicate that the proposed elasto-plastic thin-shell elements of ANCF is suitable to model complex flexible multibody systems composed of elasto-plastic thin shells and reflect their dynamic behaviors appropriately. In future researches, hyperelastic-plastic constitutive models with multiplicative decomposition will be embedded into the proposed formulations to deal with large strain problems.