Elastoplastic analysis of plane structures using improved membrane finite element with rotational DOFs

In this work, the small-strain elastoplastic behavior of structures is analyzed using an improved nonlinear finite element formulation. In this framework, an eight-node quadrilateral finite element denoted PFR8 (Plane Fiber Rotation) that belongs to the set of elements with rotational degrees of freedom is developed. Its formulation stems from the plane adaptation of the Space Fiber Rotation (SFR) concept that considers virtual rotations of nodal fiber within the element. This approach results in an enhancement of the displacement vector approximation. Von-Mises yield criteria have been applied for yielding of the materials along with the associated flow rule. Newton-Raphson method has been used to solve the nonlinear equations. To assess the performance of the proposed element, benchmark problems are addressed and the results are compared with some analytical and numerical solutions from the literature.


INTRODUCTION
onlinear problems arise in almost all disciplines of science and engineering practice. In structural mechanics, nonlinear analysis has emerged as a powerful platform for evaluation the performance of structural systems at the life safety and collapse prevention levels [1]. The literature on nonlinear structural behavior has highlighted several types of nonlinearities including geometric, material, kinematic and force nonlinearities [2]. In recent years, a large N and growing body of literature has investigated the nonlinear elastoplastic behavior of structures as it appears in a wide range of applications and industrial processes. For instance, crashworthy structures show a complex combination of structural nonlinearities in which plastic deformation is expected to contribute to the energy absorption [3]. Another example is drop test of some mobile devices such as mobile phones, personal digital assistants and laptops which end up generally with permanent deformation [4]. Due to complexity of geometry, boundary and loading conditions, there are limited practical problems in plasticity where analytical solution is available. Therefore, numerical analysis, and in particular, the finite element method has emerged as a viable alternative in addressing this kind of problems. Within this framework, accuracy and efficiency of finite elements have gained fresh prominence. To avoid the problem of singularity in the stiffness matrix and also to improve the accuracy of numerical results, there has been renewed interests in finite element with in-plane rotational degrees of freedom called also drilling rotations [5]. Since the first successful element with drilling variables proposed by Allman [6], a considerable amount of literature has been published on this subject [7][8][9][10][11][12][13]. Working along similar lines, many three-dimensional solid element formulations with rotational DOFs have been also proposed [14][15][16][17][18][19][20][21][22][23][24]. Ayad et al. [24], developed two eight-node hexahedral elements SFR8 and SFR8I based on the so-called Space Fiber Rotation concept (SFR). This approach, firstly introduced by Ayad [25], considers 3D virtual rotations of a nodal space fiber within the finite element that improves the displacement vector approximation. Similarly, many authors implemented this concept in different area of research [16,17,[23][24][25]27]. Despite the abundance of literature on the finite elements with rotational DOFs, only a few have been proposed to account for material and geometric nonlinear problems. These includes the works of Gruttmann et al. [28], Ibrahimbegovic et al. [29][30][31], Rebiai et al. [13] and Zouari et al. [32]. In the work of Zouari et al. [32], this approach has been successfully adopted to develop two four-node membrane elements, named PFR4 and PFR4I, to analyze linear and nonlinear geometric plane problems. This paper presents a new eight-node membrane finite element PFR8 in order to analyze elastoplastic nonlinear problems. The element's formulation is based on a planar adaptation of the SFR concept which results in a simpler and more concise formulation than the previous developed membrane elements with rotational DOFs. The remaining part of the paper proceeds as follows. In section 2, the formulation behind the proposed quadrilateral element is detailed. In section 3, the model is framed within the nonlinear analysis through the small strain elastoplastic constitutive model. Finally, and before highlighting the concluding remarks, numerical results relative to a series of benchmarks are presented.

BASIC CONSTITUTIVE EQUATIONS OF CONTINUUM ELASTOPLASTICITY
fter initial yielding the material behavior will be partly elastic and partly plastic. The phenomenon is characterized by an irreversible deformation where the total strain can be expressed: where   e ij  and   p ij  are the elastic and plastic components of the total strain respectively.
The elastoplastic behavior is governed by a scalar function called yield function of the form: where ij  denotes the indicial form of the stress tensor and  represents hardening parameters.
As shown in Fig. 1, the yield condition can be visualized as a surface in the stress space where the size of the surface depends on the value of the parameter  . In this paper, the elastoplastic model based on the Von-Mises associated yield criterion is adopted. Thus, the above yield condition can be expressed as: where e  is the Von-Mises effective stress and y  is the yield stress.
A Using the classical associative plastic flow rule, the plastic strain rate is assumed to be given by: Figure 1: Yield surface and normality rule in 2D stress space. [33][34] where   and P represent the plastic multiplier and the plastic flow direction, respectively.
The elastic part of strain field is linked to stress field through the relation: where e ijkl C is the tensor of elastic constants which for an isotropic material may be given as: in which G is the shear modulus;  is Poisson's ratio and ij  represents the Kronecker delta.
Using the consistency condition 0    , the plastic multiplier in Eq. (4) can be expressed as: Finally, by substituting the expression of the plastic multiplier into Eq. (5), the elastoplastic tangent modulus is derived as: where 1   for strict plastic loading and 0   for elastic loading/unloading. Therefore, the elastoplastic stress-strain relation can be expressed:

FINITE ELEMENT PROCEDURE
onsider a body of a domain  , in which the internal stresses    , boundary tractions T , and the distributed loads/unit volume v f form an equilibrium field; to undergo an arbitrary admissible displacement u  .
Therefore, the virtual work principle requires that: where    is the vector of associated virtual strains and S is the boundary on which tractions are applied.
The virtual linearized strain tensor  can be expressed in terms of the displacement vector as:

Kinematics of the PFR8 finite element
The Plane Fiber Rotation (PFR) formulation stems from the Space Fiber Rotational (SFR) concept. This concept, firstly introduced by Ayad [25], considers a 3D virtual rotations of nodal space fiber within the element that enhances the displacement vector approximation. Despite some differences (geometry, interpolation functions…), the PFR and SFR formulations proceed very much in the same way. In this section, the Plane Fiber Rotation approach is used to formulate an eight-node quadrilateral finite element called PFR8. As shown in Fig. 2, the PFR approach considers the rotation of the out of plane virtual nodal fiber iq within the element. This results in an enhancement in the displacement vector approximation. The PFR approximation of the displacement vector of a point q of the element takes the following form: is the additional displacement, i U is the vector of nodal displacements, and i N are the classical shape functions associated with the serendipity eight-node quadrilateral element. The additional displacements vector is given by: Therefore: where x and y are the Cartesian coordinates of q, i x and i y are the nodal DOFs (two displacements and one virtual rotation per node). In the matrix form, the approximation can be rewritten as: and is nodal DOFs vector of the PFR8 containing the mechanical displacements and rotations. The finite element discretization leads to write strains in the element level as: where   B is the discrete symmetric gradient which has the format: Since the shape function derivatives are expressed in terms of natural coordinates   ,   , a transformation to the physical domain   , x y should be performed. Accordingly, the shape function derivatives in terms of the physical domain coordinates   , x y are given by: where lk j are the terms of inverse Jacobian matrix. Based on Eqns. (9), (10), and (18), the virtual potential energy expression can be conveniently rearranged as: This equation results to the following expression: where e T K is the (24 × 24)-sized element stiffness matrix and e ext f is the external force vector. In this work, this system (Eq. 20) is solved by the classical frontal method [35,36].

NUMERICAL SOLUTION FOR ELASTOPLASTIC PROBLEMS
n order to solve elastoplastic problems, an algorithm based on the return mapping scheme [37,38] is adopted. This algorithm relies basically on a two-step procedure: Perform a predictor step in which we assume that the step   1 , n n t t  is elastic. Accordingly: The corresponding trial stress is given by: The flowchart in Fig. 3 presents an overall summary of the procedure. I Figure 3: Flowchart of the numerical procedure.

NUMERICAL EXAMPLES
n order to evaluate the accuracy, efficiency, and robustness of the proposed eight-node quadrilateral element PFR8, the element's formulation is implemented in the FORTRAN finite element code HYPLAS [39] to account for elastic and elastoplastic problems. In this section, representative benchmark problems are addressed. The obtained results are compared with the following reference elements:  Q8: the standard 8-node quadrilateral element with full numerical integration scheme.  HWQ8D: 8-node mixed quadrilateral element based on the Hu-Washizu functional [40].  CPS8: the plane stress eight-node quadrilateral Abaqus element with a full numerical integration scheme [41].  CPS8R: the plane stress eight-node quadrilateral Abaqus element with a reduced numerical integration scheme [41].  Q4: the standard 4-node quadrilateral element with full numerical integration scheme.  PFR4: 4-node quadrilateral element based on the concept of 'Plane Fiber Rotation' [32].  Q4PS: 4-node quadrilateral hybrid element [42].  Q4WT: 4-node quadrilateral non-conforming element [43,44].  Allman : 4-node quadrilateral element with drilling DOFs [45]. For each elastoplastic problem, the Von-Mises plasticity model is adopted where the Newton Raphson tangential stiffness method is used to solve the nonlinear resulting systems.

Cantilever beam under in-plane bending load
The first analyzed problem is a cantilever beam under pure plane bending. This example is a popular benchmark that has been studied by several authors [24,43,46] to evaluate the convergence rate and sensitivity to mesh distortion of the proposed elements. Geometry and material properties are shown in Fig. 4. As illustrated in Fig. 5, the beam is modeled with six different meshes: three structured and three distorted meshes. The maximum vertical displacement of point C, belonging to the right end side of the beam, can be obtained analytically by using the Timoshenko beam theory: 3 6 4.03 3 5 The normalized vertical displacement of point C is then determined and the obtained results are plotted in Fig. 6. For structured meshes, the PFR8 element has the best convergence rate when comparing to the other reference elements. On the other hand, the element exhibits very good performance for distorted meshes and it turns out to be less sensitive to mesh distortion.

Cook's membrane
This is a popular benchmark to evaluate the element's performance in bending dominated applications [7,47]. It consists of a trapezoidal cantilever plate completely constrained at the left-end and subjected to a uniform tangential load at the right free edge as shown in Fig. 7. The analysis is performed on the basis of six different structured meshes of 2 2 elements, with 1, … , 6. Material properties, boundary conditions and finite element mesh for 1 are illustrated in Fig. 7. The Reference solution of the vertical displacement of point A is reported by Bergan and Felippa in [7] and it is taken to be: The results of PFR8 element are compared with some reference elements in Tab. 1 where the convergence curves are depicted in Fig. 8. The findings demonstrate that the proposed PFR8 element for a given mesh density is better than the Q4 and Q8 elements and can be competitive when compared to the other reference elements. This example is also a well-known benchmark to evaluate the effectiveness of the finite element formulation in elastoplastic behavior [40,48,49]. Material properties adopted are: Poisson's ration 0.2

 
, Young modulus 70000 E  , Yield stress 0 243 y   and the isotropic hardening modulus 1000 H  . A total load 1600 F  is applied via an incremental scheme of 11 time steps where a mesh of (4 × 4) elements is adopted. Fig. 9 provides a comparison of the load-displacement curves of different elements. As can be seen, the PFR8 element results agree very well with the other reference solutions.

Cantilever beam
In this example, an elastoplastic cantilever beam with rectangular cross section is analyzed. As shown in Fig. 10, this beam is subjected to an in-plane load F at point C. Geometric as well as material properties are set out in Tab. 1. The material is assumed to be perfectly plastic where the Von-Mises model is implemented in the plane stress conditions. The cantilever is discretized into (2 × 50) regular quadrilateral finite elements. In this example, two cases of boundary conditions are studied as shown in Fig. 10. For case (a), the analytical limit load is given by Lubliner  . The evolution of the equivalent plastic strain for case (a) and (b) are depicted in Fig. 12 and Fig. 13. For case (a), it can be clearly seen that the plastic deformation starts at the node where the pinned boundary condition is applied and extends along the edge. On the other hand, the clamped edge doesn't show plastic deformation in case (b).

Thick cylinder
This example considers an infinitely long thick cylinder subjected to a gradually increasing internal pressure. This benchmark has been studied by several authors [39,51] and a reference solution was derived by Hill [52] where the limit pressure is expressed as:   The geometry of the problem as well as the adopted finite element mesh are shown in Fig. 13. The limit pressure for the present properties becomes: Owing to the symmetry of the problem, only a 30° segment of the cylinder is discretized assuming plane strain conditions along the axis of the cylinder. Material as well as geometric properties are summarized in Tab. 3.
The radial displacement at point B is plotted against the applied pressure in Fig. 14. The PFR8 results are compared to the Q8 standard element and the Hill's solution and a good correlation among the three set of results is obtained. In this example, yielding starts at the inner face of the cylinder and extends along the radius toward the outer face as depicted in Fig. 15.

CONCLUSION
n this paper, an eight-node quadrilateral element denoted PFR8 has been developed to account for elastic and elastoplastic problems. The element formulation relies on the plane adaptation of the Space Fiber Rotation (SFR) concept that considers an out of plane virtual rotation of a nodal fiber within the element which results in an enhancement in the displacement vector approximation. The present analysis focuses on small strain applications using the associative Von-Mises plasticity model. The proposed quadrilateral plane element PFR8 was implemented in HYPLAS, a Fortran finite element code developed by de Souza Neto et al. [39]. The element performance was evaluated through plane strain/stress benchmark examples. For all examples, the numerical results obtained with the PFR8 element showed good agreement with the reference solutions.