# Study on the Vector Quadrilateral Isoparametric Shell Element with Shear Deformation

Authors: Zhen Wang; Yang Zhao; Zhi Ding; Huiping Tang
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Based on the linear superposition of the quadrilateral isoparametric membrane element and quadrilateral Mindlin plate element, basic vector formulas for the quadrilateral isoparametric shell element with the shear deformation is derived. Two different coordinate modes for updating the particle displacement and the element node internal force are put forward. The uniform numerical solution for an irregular quadrilateral element is realized by the isoparametric integral scheme along the element plane. Then, the bilinear elastoplastic constitutive of material is introduced into the vector finite shell element, and the nonlinear effect of the material is realized based on the integral method along the thickness. On this basis, the nonlinear calculation and analysis program of the quadrilateral isoparametric shell element is developed. Numerical example results show that the static and dynamic analysis, large deformation and large rotation analysis, and buckling analysis can all be well performed for the shell structures by the developed program, which verifies the validity of theoretical derivation and computer program.


---

## Body

## 1. Introduction

The shell structure is a two-dimensional element structure, which is widely used in engineering structures, including hyperbolic flat shells, spherical shells, cylindrical shells, and surface structures with complex surface structures and boundaries [1–5]. Under the action of external loading, the shell structure has both in-plane tensile deformation and out-of-plane bending deformation, which can be regarded as the linear superposition of membrane elements and plate elements.The shell elements include planar shell elements and surface shell elements, and the former is relatively simple and has good precision. Green et al. [6] first adopted a triangular flat shell element and applied it to a shell structure of arbitrary shape. Zienkiewicz [7] proposed the analysis theory that combined membrane element and plate element to simulate planar shell element. Bathe and Ho [8] developed the basic theory for isoparametric thin shell element obtained from the surface of higher-order shape function, which was analyzed and compared with the calculation accuracy of thin shell element composed by linear superposition film element and the plate element. It was shown that the latter method had good accuracy and its theoretical formula was simple and easy to apply. Yang et al. [9] pointed out that by assuming the value of normal rotation stiffness to avoid stiffness matrix singularity, the obtained convergence result might be poor. It is suggested that the rotation degree of freedom should be considered separately as one of the degrees of freedom for the shell element to obtain more effective and accurate calculation results. In terms of the shear deformation of shell elements, there are some relevant literature studies. Wang et al. [10] proposed a simple first-order shear deformation shell theory (S-FSDST) for free and transient vibration analysis of composite laminated open cylindrical shells with general boundary conditions. Lezgy–Nazargah and Salahshuran [11] presented a novel mixed-field theory with a relatively low number of unknown variables for bending and vibration analysis of multilayered composite plates. Lezgy–Nazargah [12] developed a high-performance finite element model for bending and vibration analysis of thick plates based on the parametrized mixed variational principle. Lezgy–Nazargah and Meshkani [13] proposed a four-node quadrilateral partial mixed plate element with low degrees of freedom (dofs) for static and free vibration analysis of functionally graded material (FGM) plates rested on Winkler–Pasternak elastic foundations.Vector Form Intrinsic Finite Element (VFIFE) is a new dynamic calculation analysis method based on particle mechanics and vector operation mode [14–16]. This method is based on the second law of motion of a particle and belongs to a class of discrete element method. It does not involve mathematical continuous function and has a simple theory and no iteration. By means of vector operation mode, the relations among variables such as force (moment), displacement, and acceleration are expressed by particles, while the numerical calculation results of multi-particle discretization approach the approximate value of the real unique solution, and each particle is calculated separately and cyclically step by step. Therefore, the VFIFE method is more convenient for analyzing the dynamic behavior of complex structures such as large deformation, large rotation, stability, collision, and fracture. At present, some research achievements have been made on two-dimensional plate and shell element of vector finite element. Wu et al. [17] studied the vector formula theory of triangular membrane element and applied it to large deformation and large rotation analysis of membrane structure. Wang et al. [18] studied the vector formula theory of triangular thin shell element and the treatment and application of elastic-plastic materials constitutive model. However, the existing literature of VFIFE mainly focuses on the research of triangular element, and there is no relevant study on the theory and application of quadrilateral isoparametric shell element, and the influence of shear deformation on thick shell is not involved.In this paper, by linear superposition of quadrilateral isoparametric membrane element and quadrilateral Mindlin plate element, the basic vector formula of quadrilateral isoparametric shell element with shear deformation is derived firstly. Two different coordinate modes for solving particle displacement and internal force of element node and the integration of internal force for element node along element thickness and plane are then proposed. The bilinear elastoplastic material constitutive model is used to deal with the nonlinear effect of structure in vector calculation. On this basis, the nonlinear calculation and analysis program of quadrilateral isoparametric shell element is developed, and the validity of derived formulas and program codes is verified by analysis of typical examples.


## 2. Theory of Vector Quadrilateral Isoparametric Shell

The VFIFE method discretizes structure into two parts including mass particles and massless elements. And, the solving process includes three steps: the total displacement of particles, the pure deformation displacement of nodes, and the internal force of nodes. The shell element is superimposed of the quadrilateral membrane element and Mindlin plate element with shear deformation in this paper.In this section, equation (1) is obtained by using the relevant formula in literature [19]. Equations (2)–(9) are the formulas of quadrilateral isoparametric shell element under the VFIFE theoretical framework derived in this paper, without reference to other literatures. Equations (10)–(32) are the relevant formulas of quadrilateral isoparametric shell element under the VFIFE theoretical framework derived with reference to the literature [19, 20], which belong to the original content, and have been cited in Sections 2.2 and 2.3, respectively.


### 2.1. Total Displacement of Particles

The particle motion equation is obtained through Newton’s second law, which is given by the following equation:(1)Mx¨+αMx˙=F,Iθ¨+αIθ˙=Fθ,where M and I are the mass matrix and moment of inertia matrix, respectively, x is the linear coordinate vector, θ is the angular coordinate vector, F=fext+fint is the resultant force vector, Fθ=fθext+fθint is the resultant moment vector, and the superscript “ext” and “int”, respectively, represent the external force and internal force, α is the damping parameters.During each update substep, the central difference formula is used for numerical calculation, and the total particle displacement is gradually obtained through a loop substep update. For the quadrilateral isoparametric shell element in this paper,fint and fθint in equation (1) are obtained from coordinate definition mode I, x and θ in equation (1) are obtained from coordinate definition mode II, which are described in Section 3.1.


### 2.2. Pure Deformation Displacement of Element Nodes

By introducing the inverse motion concept, the rigid displacement is deducted from the total displacement of a particle to get pure deformation displacement for element nodes. The actual spatial quadrilateral of the element is transformed into a planar quadrilateral by the projection method, which is that of coordinate definition modeI. The actual coordinate vectors of the four element nodes i (i = a, b, c, d) are denoted in this section as xi, and the unit edge vector of the spatial quadrilateral ekl is given by the following equation:(2)ekl=xl−xkxl−xkkl=ab,bc,cd,da,where k and l are two adjacent points of the spatial quadrilateral.The unit normal vector of the plane corresponding to two adjacent edges of the spatial quadrilateralnj0 is given by the following equation:(3)nj0=ekl×emnekl×emnkl,mn=ab,bc,cd,da;j=1,2,3,4,where kl and mn are two adjacent edges of the spatial quadrilateral.The unit mean normal vector of the projection planen0 is given by the following equation:(4)n0=nn=14∑j=14nj0/14∑j=14nj0j=1,2,3,4,where n=1/4∑j=14nj0 is the mean normal vector of the projection plane.The projection coordinate vectors on the projection plane of element nodei (i = a, b, c, d) are denoted as x¯i, which is given by the following equation:(5)xi¯=xc+dxi−dxi⋅n0n0,where xC is the coordinate vector of the centroid point C, dxi=xi−xC is the distance vector for element nodes i relative to centroid point C.The projection distance vectors of projection points of element nodesi relative to centroid point C are denoted as dx¯i, which is given by the following equation:(6)dx¯i=x¯i−xc,where x¯i is the coordinate vector of the projection point of the element node i.The actual angular coordinate vectors of the four element nodesi (i = a, b, c, d) are denoted as θi, then the rotational angle of the element node caused by projection dθi (that is, the rotation angle from dxi to dx¯i with unit vector npi0 as the rotation axis) is given by the following equation:(7)dθi=cos−1dxi⋅dx¯idxidx¯i,dθi∈0,π,where npi0=dxi×dx¯i/dxi×dx¯i.The rotational angular displacement vector of element nodedθi is given by the following equation:(8)dθi=dθinpi0.The projection angular coordinate vector of element nodeθ¯i is given by the following equation:(9)θ¯i=θi+dθi.After projection method processing, the spatial quadrilateral element is transformed into the planar quadrilateral element, corresponding to the coordinate definition modeI (x¯i,x¯i′,θ¯i,θ¯i′,n,n′), where x¯i and x¯i′ are the coordinate vectors for projection points of element nodes i corresponding to the beginning and end of particle motion time substep, θ¯i and θ¯i′ are the angular coordinate vectors for projection points of element nodes i corresponding to the beginning and end of particle motion time substep, n and n′ are the mean normal vectors of the projection planes corresponding to the beginning and end of particle motion time substep. The total linear displacement vectors and total angular displacement vectors of element nodes i are, respectively, ui=x¯i′−x¯i and uθi=θ¯i′−θ¯i.The spatial motion for the projection plane (i.e., planar quadrilateral) includes rigid body translation, out-of-plane rigid body rotation, in-plane rigid body rotation, and pure deformation displacement. The specific derivation is similar to the triangle plate element theory in reference [19], and only the main derivation formulas are given in this paper.Rigid body translation is taken as the total linear displacements vectorua for reference node a, then the relative linear displacements Δηi and relative angular displacements Δηθi for element node i after deducting the rigid body translation are given by the following equation:(10)Δηa=0,Δηi=ui−ua,i=b,c,d,Δηθi=uθi,i=a,b,c,d,Rigid body rotation is estimated based on the projection plane, including out-of-plane rigid body rotation and in-plane rigid body rotation. The opposite direction of rotation is inverse rotation and the linear displacementΔηi−opr and angular displacement Δηθi−opr of the out-of-plane reverse rigid body rotation are given by the following equation:(11)θ1Δηa−opr=0,Δηi−opr=RopT−θ1−Irai″=Rop∗−θ1⋅rai″,i=b,c,d,Δηθi−opr=−θ1=−θ1nop0,i=a,b,c,d,where θ1 is the out-of-plane rotational angle of projection plane (that is, the rotational angle from n to n″ with unit vector nop0=n×n″/n×n″=lopmopnopT as the rotational axis) which is given by equation (12), −θ1 is the out-of-plane inverse rotational angle of the projection plane, θ1 is the out-of-plane rotational angle vector of the projection plane, Rop∗−θ1 is the out-of-plane inverse rotation matrix which is given by equation (13), rai″=x¯i″−x¯a″=x¯i′−x¯a′ is the edge vector of the projection planar element (i.e., planar element consisting of projection points of element nodes) after deducting the rigid body translation, xi″ is the coordinate vector of the projection point of element node i after deducting the rigid body translation.(12)n″θ1=cos−1n⋅n″nn″,θ1∈0,π,where n″ are the mean normal vectors of the projection planes after deducting the rigid body translation which can be obtained by substituting x¯i″ for xi in equation (2)∼equation (4).(13)Rop∗−θ1=1−cos−θ1Aop2+sin−θ1Aop,where Aop=0−nopmopnop0−lop−moplop0.The linear displacementsΔηi−ipr and angular displacements Δηθi−ipr for the in-plane reverse rigid body rotation are given by the following equation:(14)Δηa−ipr=0,Δηi−ipr=RipT−θ2−Irai″=Rip∗−θ2⋅rai″,i=b,c,d,Δηθi−ipr=−θ2=−θ2nip0,i=a,b,c,d,where θ2 is the in-plane rotational angle of projection plane (that is, the average in-plane rotational angle of the projection plane with the unit vector nip0=n0=n/n=lipmipnipT as the rotation axis) which is given by equation (15), −θ2 is the in-plane inverse rotational angle of the projection plane, θ2 is the in-plane rotational angle vector of the projection plane, Rip∗−θ2 is the in-plane inverse rotation matrix which is given by equation (18), rai‴=x¯i‴−x¯a‴=x¯i′−x¯a′+Δηi−opr is the edge vector of the projection planar element (i.e., planar element consisting of projection points of element nodes) after deducting the rigid body translation and the out-of-plane rigid body rotation, xi‴ is the coordinate vector of the projection point of element node i after deducting the rigid body translation and the out-of-plane rigid body rotation.(15)θ2=∆φa+∆φb+∆φc+∆φd4,θ2∈−π,π,where Δφi (i = a, b, c, d) is the in-plane rotational angle of element i which is given by the following equations:(16)Δφi=cos−1ei⋅ei‴,i=a,b,c,d,where ei=x¯i−x¯C/x¯i−x¯C, ei‴=x¯i‴−x¯C/x¯i‴−x¯C, x¯C is the coordinate vector of the projection point of centroid point C.(17)∆φi=Δφi,whenei×ei‴andn0areinsamedirection,∆φi=Δφi,whenei×ei‴andn0areinoppositedirection,(18)Rip∗−θ2=1−cos−θ2Aip2+sin−θ2Aip,where Aip=0−nipmipnip0−lip−miplip0.After deducting the rigid body translation, the out-of-plane rigid body rotation, and the in-plane rigid body rotation of the projection plane, the pure deformation linear displacementΔηid and pure deformation angular displacement Δηθid of element node i are given by the following equation:(19)Δηad=0,Δηid=ui−ua+Δηi−opr+Δηi−ipr,i=b,c,d,Δηθid=uθi+Δηθi−opr+Δηθi−ipr,i=a,b,c,d,


### 2.3. Internal Force of Element Nodes

When solving the internal forces of element nodes, the VFIFE method transforms the spatial element problem into a planar element problem by defining the deformation coordinate system. The reference element node a is defined as the origin point in this section for the deformation coordinate system, and thex^ axis is defined as the direction along the edge ab of the element, which results in u^a=v^a=0 and w^i=0 (i = a, b, c, d).The transformational relation between the global coordinate system and the deformation coordinate system is given by the following equation:(20)x^=Q^x−xa,where x=xyz is the linear coordinate vector in the deformation coordinate system, Q^ is the transformation matrix from the global coordinate system to the deformation coordinate system which is given by equation (21). Since z^=0, x^ could be simplified to x^=x^y^T.(21)Q^=e^1e^2e^3T,where e^1=x¯b−x¯a/x¯b−x¯a,e^3=n0,e^2=e^3×e^1.The pure deformation linear displacements vectoru^i and pure deformation angular displacements vector u^θi in the deformation coordinate system for the element node i are given by the following equation:(22)u^i=Q^∆ηid,i=a,b,c,d,u^θi=Q^∆ηθid,i=a,b,c,d,where u^i=u^iv^iw^iT, u^θi=u^θiv^θiw^θi. Since w^i=0 and w^θi can be not considered in the solving process of the internal forces of element nodes, u^i and u^θi can be simplified as u^i=u^iv^iT and u^θi=u^θiv^θiT, respectively.The pure deformation linear displacementsu^i,v^i of element nodes are applied for calculating the node internal forces for membrane element part, the pure deformation linear displacements w^i and pure deformation angular displacements u^θi,v^θi of element nodes are applied for calculating the node internal forces and node internal moments for plate element part. Then, the deformation virtual work equation of the quadrilateral shell element is given by the following equation:(23)∑iδu^m_iTf^m_i+δu^θb_i′Tf^θb_i′=∫VδΔε^mTσ^mdV+δΔε^bTσ^bdV,where u^m_i=u^iv^iT and u^θb_i′=w^iu^θiv^θiT are the pure deformation linear displacement vector of node i for membrane element part and the pure deformation linear (angular) displacement vector of node i for plate element part in the deformation coordinate system, respectively.By introducing the generalized form function of the quadrilateral isoparametric element in traditional finite element, the pure deformation linear displacements vectoru^ and pure deformation angular displacements vector u^θ are shown as equation (24) at any position of elements.(24)u^=∑i=a,b,c,dΝiu^i,u^θ=∑i=a,b,c,dΝiu^θi,where Nix^,y^,z^ is the element form function, whose deformation form is the same as the traditional finite element.For the quadrilateral isoparametric element with four nodes, the square element coordinatesξ,η in the element coordinate system are introduced to solve the relevant variable values in the deformation system, and the expression of the element shape function Nix^,y^,z^ is shown as the following equation:(25)Nix^,y^,z^=Niξ,η=141+ξiξ1+ηiη,where ξi and ηi are the square element coordinates of the four nodes in the element coordinate system, and the value is ±1.0. The transformation of the element coordinate system and deformation coordinate system is completed in this paper by Jacobian matrix J, whose expression is given in reference [20].The node pure deformation linear displacement combines the vector of the membrane element partu^m∗ and the node pure deformation linear (angular) displacement combines the vector of the plate element part u^′b∗ in the deformation coordinate system are given by the following equation:(26)u^m∗=u^bTu^cTu^dTT,u^′b∗=u^′θaTu^′θbTu^′θcTu^′θdTT,where u^i=u^iv^iT, u^θi′=w^iu^θiv^θiT=0u^θiv^θiT.The tensile strain vector of the membrane element partΔε^m and bending (including shear) strain vector of the plate element part Δε^b are given by the following equation:(27)Δε^m=Bm∗u^m∗,Δε^b=B′b∗u^′b∗=z^Bb∗Bs∗u^′b∗,where Δε^m=Δε^mxΔε^myΔγ^mxyT, Δε^b=Δε^bxΔε^byΔγ^bxyΔγ^byzΔγ^bzxT, Bm∗ and B′b∗ are the tensile strain-displacement relationship matrix of the membrane element part and bending (including shear) strain-displacement relationship matrix of the plate element part, respectively, Bb∗ and Bs∗ are bending strain-displacement relationship matrix and shear strain-displacement relationship matrix of plate element part. The detailed expressions of Bm∗, B′b∗, Bb∗, and Bs∗ are given in reference [20].The material stress-strain relationship matrices (i.e., constitutive matrices) of the membrane element part and plate element part are denoted asD and D′, respectively, which could be linear elastic or elastoplastic. Then, stress vectors for the membrane element part and plate element part are given by the following equation:(28)Δσ^m=DΔε^m=DbBm∗u^m∗,Δσ^b=D′Δε^b=DbDsz^Bb∗Bs∗u^′b∗,where Δσ^m=Δσ^mxΔσ^myΔτ^mxyT, Δσ^b=Δσ^bxΔσ^byΔτ^bxyΔτ^byzΔτ^bzxT, Db and Ds are the constitutive matrices of plane strain (including tensile and bending) and shear strain, respectively.For the case of elastic materials,Db and Ds are given by the following equation:(29)Db=E1−υ21υ0υ10001−υ2Ds=kG1001,,where E is Young modulus, G=E/21+υ is shear modulus, υ is Poisson's ratio, k is the shear correction coefficient, and k=5/6 is generally taken to consider the shear non-uniform change in the direction of shell thickness.Considering the influence of transverse shear stressτ^byz and τ^bzx in quadrilateral isoparametric shell element, the right-hand side formula for deformation virtual work equation (23) can be expressed in a more detailed form, which is given by the following equation:(30)δU=∫VδΔε^mTσ^m0+Δσ^mdV+∫VδΔε^bTσ^b0+Δσ^bdV=δu^m∗T∫VBm∗Tσ^m0dV+∫VBm∗TDbBm∗dVu^m∗+δu^′b∗T∫Vz^Bb∗Bs∗Tσ^b0dV+∫Vz^Bb∗Bs∗TDbDsz^Bb∗Bs∗dVu^′b∗,where σ^m0 and σ^b0 are the initial stresses of the membrane element part and plate element part.According to the left-hand side formula of the deformation virtual work equation (23), the internal force vector of the membrane element part and plate element part f^m∗ and internal force (moment) vector of plate element part f^′θb∗ at the end of particle motion time substep can be obtained by the following equation:(31)f^m∗=f^m_bTf^m_cTf^m_dTT=f^m0∗+Δf^m∗=∫VBm∗Tσ^m0dV+∫A^a∫−ta/2ta/2Bm∗TDbBm∗dz^dA^au^m∗,f^′θb∗=f^′θb_bTf^′θb_cTf^′θb_dTT=f^′θb0∗+Δf^′θb∗=∫Vz^Bb∗Bs∗Tσ^b0dV+∫A^a∫−ta/2ta/2z^2Bb∗TDbBb∗dz^dA^a+∫A^a∫−ta/2ta/2Bs∗TDsBs∗dz^dA^au^′b∗,where f^m_i=f^ixf^iyT is the internal force of node i for membrane element part and plate element part, f^θb_i′=f^izf^θixf^θiyT is the internal force (moment) of node i for plate element part, ta is the thickness of ehe shell element, f^m0∗ is the internal force vector of the membrane element part and plate element part at the beginning of the particle motion time substep, f^′θb0∗ is the internal force (moment) vector of plate element part at the beginning of particle motion time substep, Δf^m∗ is the internal force increment vector of the membrane element part and plate element part during the particle motion time substep, Δf^′θb∗ is the internal force (moment) increment vector of the plate element part during the particle motion time substep. In addition, the internal force f^m_a=f^axf^ayT at the reference node a could be obtained from the static equilibrium condition under planar conditions.The internal forces (moments) for the nodei (i = a, b, c, d) in the above-given formula are superimposed and extended to obtain the internal force f^i=f^ixf^iyf^izT and internal moment f^θi=f^θixf^θiyf^θizT of 3D shell element, where f^θiz=0. Then, the internal forces (moments) for the node i in the global coordinate system is obtained system through coordinate system transformation from the deformation coordinate. And, the internal force vector fi and internal moment fθi for the element node i at the time t are obtained through forward motion transformation, which is given by the following equation:(32)fi=Ripθ2Ropθ2Q^Tf^i,fθi=Ripθ2Ropθ2Q^Tf^θi,where Ropθ2=Rop∗θ2+IT, Ripθ2=Rip∗θ2+IT, Rop∗θ2 and Rip∗θ2 are the out-of-plane and in-plane forward rotational matrices, which are obtained according to equation (13) and equation (18). The internal force fi and internal force (moment) fθi of element node i are applied to particle i in reverse order to obtain the internal force fiint and internal force (moment) fθiint transferred from shell element to particle.It is worth noting that the research topic of this paper is to develop a quadrilateral shell element considering shear deformation under the theoretical framework of the new analysis method VFIFE, rather than to develop a quadrilateral shell element considering shear deformation with better performance based on the traditional finite element method. VFIFE is a new analysis method based on particle motion, which has great advantages in structure discontinuity behavior, dynamic behavior, and other aspects (see Section1 for details). However, there is no relevant literature on the quadrilateral shell element theory under the VFIFE method. In this paper, considering from the relatively simple quadrilateral shell element (that is, the linear superposition of the quadrilateral isoparametric membrane element and quadrilateral Mindlin plate element), so as to promote and improve the theoretical development and application of VFIFE, and give full play to the advantages of VFIFE method compared with the traditional FE method. Other refined and high-order shell element models considering shear deformation will be further studied in the future.


## 2.1. Total Displacement of Particles

The particle motion equation is obtained through Newton’s second law, which is given by the following equation:(1)Mx¨+αMx˙=F,Iθ¨+αIθ˙=Fθ,where M and I are the mass matrix and moment of inertia matrix, respectively, x is the linear coordinate vector, θ is the angular coordinate vector, F=fext+fint is the resultant force vector, Fθ=fθext+fθint is the resultant moment vector, and the superscript “ext” and “int”, respectively, represent the external force and internal force, α is the damping parameters.During each update substep, the central difference formula is used for numerical calculation, and the total particle displacement is gradually obtained through a loop substep update. For the quadrilateral isoparametric shell element in this paper,fint and fθint in equation (1) are obtained from coordinate definition mode I, x and θ in equation (1) are obtained from coordinate definition mode II, which are described in Section 3.1.


## 2.2. Pure Deformation Displacement of Element Nodes

By introducing the inverse motion concept, the rigid displacement is deducted from the total displacement of a particle to get pure deformation displacement for element nodes. The actual spatial quadrilateral of the element is transformed into a planar quadrilateral by the projection method, which is that of coordinate definition modeI. The actual coordinate vectors of the four element nodes i (i = a, b, c, d) are denoted in this section as xi, and the unit edge vector of the spatial quadrilateral ekl is given by the following equation:(2)ekl=xl−xkxl−xkkl=ab,bc,cd,da,where k and l are two adjacent points of the spatial quadrilateral.The unit normal vector of the plane corresponding to two adjacent edges of the spatial quadrilateralnj0 is given by the following equation:(3)nj0=ekl×emnekl×emnkl,mn=ab,bc,cd,da;j=1,2,3,4,where kl and mn are two adjacent edges of the spatial quadrilateral.The unit mean normal vector of the projection planen0 is given by the following equation:(4)n0=nn=14∑j=14nj0/14∑j=14nj0j=1,2,3,4,where n=1/4∑j=14nj0 is the mean normal vector of the projection plane.The projection coordinate vectors on the projection plane of element nodei (i = a, b, c, d) are denoted as x¯i, which is given by the following equation:(5)xi¯=xc+dxi−dxi⋅n0n0,where xC is the coordinate vector of the centroid point C, dxi=xi−xC is the distance vector for element nodes i relative to centroid point C.The projection distance vectors of projection points of element nodesi relative to centroid point C are denoted as dx¯i, which is given by the following equation:(6)dx¯i=x¯i−xc,where x¯i is the coordinate vector of the projection point of the element node i.The actual angular coordinate vectors of the four element nodesi (i = a, b, c, d) are denoted as θi, then the rotational angle of the element node caused by projection dθi (that is, the rotation angle from dxi to dx¯i with unit vector npi0 as the rotation axis) is given by the following equation:(7)dθi=cos−1dxi⋅dx¯idxidx¯i,dθi∈0,π,where npi0=dxi×dx¯i/dxi×dx¯i.The rotational angular displacement vector of element nodedθi is given by the following equation:(8)dθi=dθinpi0.The projection angular coordinate vector of element nodeθ¯i is given by the following equation:(9)θ¯i=θi+dθi.After projection method processing, the spatial quadrilateral element is transformed into the planar quadrilateral element, corresponding to the coordinate definition modeI (x¯i,x¯i′,θ¯i,θ¯i′,n,n′), where x¯i and x¯i′ are the coordinate vectors for projection points of element nodes i corresponding to the beginning and end of particle motion time substep, θ¯i and θ¯i′ are the angular coordinate vectors for projection points of element nodes i corresponding to the beginning and end of particle motion time substep, n and n′ are the mean normal vectors of the projection planes corresponding to the beginning and end of particle motion time substep. The total linear displacement vectors and total angular displacement vectors of element nodes i are, respectively, ui=x¯i′−x¯i and uθi=θ¯i′−θ¯i.The spatial motion for the projection plane (i.e., planar quadrilateral) includes rigid body translation, out-of-plane rigid body rotation, in-plane rigid body rotation, and pure deformation displacement. The specific derivation is similar to the triangle plate element theory in reference [19], and only the main derivation formulas are given in this paper.Rigid body translation is taken as the total linear displacements vectorua for reference node a, then the relative linear displacements Δηi and relative angular displacements Δηθi for element node i after deducting the rigid body translation are given by the following equation:(10)Δηa=0,Δηi=ui−ua,i=b,c,d,Δηθi=uθi,i=a,b,c,d,Rigid body rotation is estimated based on the projection plane, including out-of-plane rigid body rotation and in-plane rigid body rotation. The opposite direction of rotation is inverse rotation and the linear displacementΔηi−opr and angular displacement Δηθi−opr of the out-of-plane reverse rigid body rotation are given by the following equation:(11)θ1Δηa−opr=0,Δηi−opr=RopT−θ1−Irai″=Rop∗−θ1⋅rai″,i=b,c,d,Δηθi−opr=−θ1=−θ1nop0,i=a,b,c,d,where θ1 is the out-of-plane rotational angle of projection plane (that is, the rotational angle from n to n″ with unit vector nop0=n×n″/n×n″=lopmopnopT as the rotational axis) which is given by equation (12), −θ1 is the out-of-plane inverse rotational angle of the projection plane, θ1 is the out-of-plane rotational angle vector of the projection plane, Rop∗−θ1 is the out-of-plane inverse rotation matrix which is given by equation (13), rai″=x¯i″−x¯a″=x¯i′−x¯a′ is the edge vector of the projection planar element (i.e., planar element consisting of projection points of element nodes) after deducting the rigid body translation, xi″ is the coordinate vector of the projection point of element node i after deducting the rigid body translation.(12)n″θ1=cos−1n⋅n″nn″,θ1∈0,π,where n″ are the mean normal vectors of the projection planes after deducting the rigid body translation which can be obtained by substituting x¯i″ for xi in equation (2)∼equation (4).(13)Rop∗−θ1=1−cos−θ1Aop2+sin−θ1Aop,where Aop=0−nopmopnop0−lop−moplop0.The linear displacementsΔηi−ipr and angular displacements Δηθi−ipr for the in-plane reverse rigid body rotation are given by the following equation:(14)Δηa−ipr=0,Δηi−ipr=RipT−θ2−Irai″=Rip∗−θ2⋅rai″,i=b,c,d,Δηθi−ipr=−θ2=−θ2nip0,i=a,b,c,d,where θ2 is the in-plane rotational angle of projection plane (that is, the average in-plane rotational angle of the projection plane with the unit vector nip0=n0=n/n=lipmipnipT as the rotation axis) which is given by equation (15), −θ2 is the in-plane inverse rotational angle of the projection plane, θ2 is the in-plane rotational angle vector of the projection plane, Rip∗−θ2 is the in-plane inverse rotation matrix which is given by equation (18), rai‴=x¯i‴−x¯a‴=x¯i′−x¯a′+Δηi−opr is the edge vector of the projection planar element (i.e., planar element consisting of projection points of element nodes) after deducting the rigid body translation and the out-of-plane rigid body rotation, xi‴ is the coordinate vector of the projection point of element node i after deducting the rigid body translation and the out-of-plane rigid body rotation.(15)θ2=∆φa+∆φb+∆φc+∆φd4,θ2∈−π,π,where Δφi (i = a, b, c, d) is the in-plane rotational angle of element i which is given by the following equations:(16)Δφi=cos−1ei⋅ei‴,i=a,b,c,d,where ei=x¯i−x¯C/x¯i−x¯C, ei‴=x¯i‴−x¯C/x¯i‴−x¯C, x¯C is the coordinate vector of the projection point of centroid point C.(17)∆φi=Δφi,whenei×ei‴andn0areinsamedirection,∆φi=Δφi,whenei×ei‴andn0areinoppositedirection,(18)Rip∗−θ2=1−cos−θ2Aip2+sin−θ2Aip,where Aip=0−nipmipnip0−lip−miplip0.After deducting the rigid body translation, the out-of-plane rigid body rotation, and the in-plane rigid body rotation of the projection plane, the pure deformation linear displacementΔηid and pure deformation angular displacement Δηθid of element node i are given by the following equation:(19)Δηad=0,Δηid=ui−ua+Δηi−opr+Δηi−ipr,i=b,c,d,Δηθid=uθi+Δηθi−opr+Δηθi−ipr,i=a,b,c,d,


## 2.3. Internal Force of Element Nodes

When solving the internal forces of element nodes, the VFIFE method transforms the spatial element problem into a planar element problem by defining the deformation coordinate system. The reference element node a is defined as the origin point in this section for the deformation coordinate system, and thex^ axis is defined as the direction along the edge ab of the element, which results in u^a=v^a=0 and w^i=0 (i = a, b, c, d).The transformational relation between the global coordinate system and the deformation coordinate system is given by the following equation:(20)x^=Q^x−xa,where x=xyz is the linear coordinate vector in the deformation coordinate system, Q^ is the transformation matrix from the global coordinate system to the deformation coordinate system which is given by equation (21). Since z^=0, x^ could be simplified to x^=x^y^T.(21)Q^=e^1e^2e^3T,where e^1=x¯b−x¯a/x¯b−x¯a,e^3=n0,e^2=e^3×e^1.The pure deformation linear displacements vectoru^i and pure deformation angular displacements vector u^θi in the deformation coordinate system for the element node i are given by the following equation:(22)u^i=Q^∆ηid,i=a,b,c,d,u^θi=Q^∆ηθid,i=a,b,c,d,where u^i=u^iv^iw^iT, u^θi=u^θiv^θiw^θi. Since w^i=0 and w^θi can be not considered in the solving process of the internal forces of element nodes, u^i and u^θi can be simplified as u^i=u^iv^iT and u^θi=u^θiv^θiT, respectively.The pure deformation linear displacementsu^i,v^i of element nodes are applied for calculating the node internal forces for membrane element part, the pure deformation linear displacements w^i and pure deformation angular displacements u^θi,v^θi of element nodes are applied for calculating the node internal forces and node internal moments for plate element part. Then, the deformation virtual work equation of the quadrilateral shell element is given by the following equation:(23)∑iδu^m_iTf^m_i+δu^θb_i′Tf^θb_i′=∫VδΔε^mTσ^mdV+δΔε^bTσ^bdV,where u^m_i=u^iv^iT and u^θb_i′=w^iu^θiv^θiT are the pure deformation linear displacement vector of node i for membrane element part and the pure deformation linear (angular) displacement vector of node i for plate element part in the deformation coordinate system, respectively.By introducing the generalized form function of the quadrilateral isoparametric element in traditional finite element, the pure deformation linear displacements vectoru^ and pure deformation angular displacements vector u^θ are shown as equation (24) at any position of elements.(24)u^=∑i=a,b,c,dΝiu^i,u^θ=∑i=a,b,c,dΝiu^θi,where Nix^,y^,z^ is the element form function, whose deformation form is the same as the traditional finite element.For the quadrilateral isoparametric element with four nodes, the square element coordinatesξ,η in the element coordinate system are introduced to solve the relevant variable values in the deformation system, and the expression of the element shape function Nix^,y^,z^ is shown as the following equation:(25)Nix^,y^,z^=Niξ,η=141+ξiξ1+ηiη,where ξi and ηi are the square element coordinates of the four nodes in the element coordinate system, and the value is ±1.0. The transformation of the element coordinate system and deformation coordinate system is completed in this paper by Jacobian matrix J, whose expression is given in reference [20].The node pure deformation linear displacement combines the vector of the membrane element partu^m∗ and the node pure deformation linear (angular) displacement combines the vector of the plate element part u^′b∗ in the deformation coordinate system are given by the following equation:(26)u^m∗=u^bTu^cTu^dTT,u^′b∗=u^′θaTu^′θbTu^′θcTu^′θdTT,where u^i=u^iv^iT, u^θi′=w^iu^θiv^θiT=0u^θiv^θiT.The tensile strain vector of the membrane element partΔε^m and bending (including shear) strain vector of the plate element part Δε^b are given by the following equation:(27)Δε^m=Bm∗u^m∗,Δε^b=B′b∗u^′b∗=z^Bb∗Bs∗u^′b∗,where Δε^m=Δε^mxΔε^myΔγ^mxyT, Δε^b=Δε^bxΔε^byΔγ^bxyΔγ^byzΔγ^bzxT, Bm∗ and B′b∗ are the tensile strain-displacement relationship matrix of the membrane element part and bending (including shear) strain-displacement relationship matrix of the plate element part, respectively, Bb∗ and Bs∗ are bending strain-displacement relationship matrix and shear strain-displacement relationship matrix of plate element part. The detailed expressions of Bm∗, B′b∗, Bb∗, and Bs∗ are given in reference [20].The material stress-strain relationship matrices (i.e., constitutive matrices) of the membrane element part and plate element part are denoted asD and D′, respectively, which could be linear elastic or elastoplastic. Then, stress vectors for the membrane element part and plate element part are given by the following equation:(28)Δσ^m=DΔε^m=DbBm∗u^m∗,Δσ^b=D′Δε^b=DbDsz^Bb∗Bs∗u^′b∗,where Δσ^m=Δσ^mxΔσ^myΔτ^mxyT, Δσ^b=Δσ^bxΔσ^byΔτ^bxyΔτ^byzΔτ^bzxT, Db and Ds are the constitutive matrices of plane strain (including tensile and bending) and shear strain, respectively.For the case of elastic materials,Db and Ds are given by the following equation:(29)Db=E1−υ21υ0υ10001−υ2Ds=kG1001,,where E is Young modulus, G=E/21+υ is shear modulus, υ is Poisson's ratio, k is the shear correction coefficient, and k=5/6 is generally taken to consider the shear non-uniform change in the direction of shell thickness.Considering the influence of transverse shear stressτ^byz and τ^bzx in quadrilateral isoparametric shell element, the right-hand side formula for deformation virtual work equation (23) can be expressed in a more detailed form, which is given by the following equation:(30)δU=∫VδΔε^mTσ^m0+Δσ^mdV+∫VδΔε^bTσ^b0+Δσ^bdV=δu^m∗T∫VBm∗Tσ^m0dV+∫VBm∗TDbBm∗dVu^m∗+δu^′b∗T∫Vz^Bb∗Bs∗Tσ^b0dV+∫Vz^Bb∗Bs∗TDbDsz^Bb∗Bs∗dVu^′b∗,where σ^m0 and σ^b0 are the initial stresses of the membrane element part and plate element part.According to the left-hand side formula of the deformation virtual work equation (23), the internal force vector of the membrane element part and plate element part f^m∗ and internal force (moment) vector of plate element part f^′θb∗ at the end of particle motion time substep can be obtained by the following equation:(31)f^m∗=f^m_bTf^m_cTf^m_dTT=f^m0∗+Δf^m∗=∫VBm∗Tσ^m0dV+∫A^a∫−ta/2ta/2Bm∗TDbBm∗dz^dA^au^m∗,f^′θb∗=f^′θb_bTf^′θb_cTf^′θb_dTT=f^′θb0∗+Δf^′θb∗=∫Vz^Bb∗Bs∗Tσ^b0dV+∫A^a∫−ta/2ta/2z^2Bb∗TDbBb∗dz^dA^a+∫A^a∫−ta/2ta/2Bs∗TDsBs∗dz^dA^au^′b∗,where f^m_i=f^ixf^iyT is the internal force of node i for membrane element part and plate element part, f^θb_i′=f^izf^θixf^θiyT is the internal force (moment) of node i for plate element part, ta is the thickness of ehe shell element, f^m0∗ is the internal force vector of the membrane element part and plate element part at the beginning of the particle motion time substep, f^′θb0∗ is the internal force (moment) vector of plate element part at the beginning of particle motion time substep, Δf^m∗ is the internal force increment vector of the membrane element part and plate element part during the particle motion time substep, Δf^′θb∗ is the internal force (moment) increment vector of the plate element part during the particle motion time substep. In addition, the internal force f^m_a=f^axf^ayT at the reference node a could be obtained from the static equilibrium condition under planar conditions.The internal forces (moments) for the nodei (i = a, b, c, d) in the above-given formula are superimposed and extended to obtain the internal force f^i=f^ixf^iyf^izT and internal moment f^θi=f^θixf^θiyf^θizT of 3D shell element, where f^θiz=0. Then, the internal forces (moments) for the node i in the global coordinate system is obtained system through coordinate system transformation from the deformation coordinate. And, the internal force vector fi and internal moment fθi for the element node i at the time t are obtained through forward motion transformation, which is given by the following equation:(32)fi=Ripθ2Ropθ2Q^Tf^i,fθi=Ripθ2Ropθ2Q^Tf^θi,where Ropθ2=Rop∗θ2+IT, Ripθ2=Rip∗θ2+IT, Rop∗θ2 and Rip∗θ2 are the out-of-plane and in-plane forward rotational matrices, which are obtained according to equation (13) and equation (18). The internal force fi and internal force (moment) fθi of element node i are applied to particle i in reverse order to obtain the internal force fiint and internal force (moment) fθiint transferred from shell element to particle.It is worth noting that the research topic of this paper is to develop a quadrilateral shell element considering shear deformation under the theoretical framework of the new analysis method VFIFE, rather than to develop a quadrilateral shell element considering shear deformation with better performance based on the traditional finite element method. VFIFE is a new analysis method based on particle motion, which has great advantages in structure discontinuity behavior, dynamic behavior, and other aspects (see Section1 for details). However, there is no relevant literature on the quadrilateral shell element theory under the VFIFE method. In this paper, considering from the relatively simple quadrilateral shell element (that is, the linear superposition of the quadrilateral isoparametric membrane element and quadrilateral Mindlin plate element), so as to promote and improve the theoretical development and application of VFIFE, and give full play to the advantages of VFIFE method compared with the traditional FE method. Other refined and high-order shell element models considering shear deformation will be further studied in the future.


## 3. Key Problems and Corresponding Solutions in Numerical Calculation

### 3.1. Coordinate Definition Modes

The coordinate definition mode I refers to the coordinates of the element nodes when the pure deformation linear (angular) displacement and internal force (moment) of shell element nodes are solved from timet0 to t, so they need to be on the same reference plane. The four-node shell element is formed into a spatial quadrilateral after undergoing motion and deformation. The projection coordinates (x¯i,x¯i′ and θ¯i,θ¯i′) for the element nodes on the reference planes which are perpendicular to the mean normal vector (n,n′) of the element and passing through the centroid points (C,C′) of the original element can be used as the coordinate definition mode I to realize the transformation of the spatial quadrilateral to the planar quadrilateral.Coordinate definition mode II refers to the coordinates of each particle when the central difference formulas are applied for solving the particle motion equation (1), so it needs to be a unique value. The following two methods can be adopted: (1) the actual coordinates of the particle corresponding to the spatial quadrilateral after the element movement and deformation; (2) I element nodes corresponding to the same single particle in coordinate definition mode I obtained by plane projection do not coincide with each other, so the average coordinate value (x=1/k∑i=1kx¯i,x′=1/k∑i=1kx¯i′ and θ=1/k∑i=1kθ¯i,θ′=1/k∑i=1kθ¯i′) of each element nodes corresponding to the same single particle in coordinate definition mode I can be taken as coordinate definition mode II, where k is the total number for that corresponding element nodes connected to each element by the same single particle.


### 3.2. Integral of Internal Forces for Element Nodes

When solving the node internal forces and node internal moments for the membrane element part and plate element part of the quadrilateral isoparametric shell element according to equation (31), the numerical integration operation will be involved since the internal stresses and strains of the elements change with their positions. Considering the case of general elastic-plastic material, the integrations along the thickness and element plane are solved by Newton–Cotes integral scheme and 2D Gaussian integral scheme in the traditional finite element method, respectively, [21]. The numerical integration formulas of Δf^m∗ and Δf^′θb∗ are given by the following equation:(33)Δf^m∗=∫A^a∫−ta/2ta/2Bm∗TDbBm∗dz^dA^au^m∗=∫−11∫−11∫−ta/2ta/2Bm∗ξ,ηTDbξ,η,z^Bm∗ξ,ηdz^Jξ,ηdξdηu^m∗=∫−11∫−11∑k=1n1HkFm1ξ,η,z^kJξ,ηdξdηu^m∗≅∑i,j=1n2HijFm2ξi,ηju^m∗,Δf^′θb∗=∫A^a∫−ta/2ta/2z^2Bb∗TDbBb∗dz^dA^a+∫A^a∫−ta/2ta/2Bs∗TDsBs∗dz^dA^au^′b∗=∫−11∫−11∫−ta/2ta/2z^2Bb∗ξ,ηTDbξ,η,z^Bb∗ξ,ηdz^Jξ,ηdξdη+∫−11∫−11∫−ta/2ta/2Bs∗ξ,ηTDsξ,η,z^Bs∗ξ,ηdz^Jξ,ηdξdηu^′b∗=∫−11∫−11∑kn1Hkz^k2Fb1ξ,η,z^kJξ,ηdξdη+∫−11∫−11∑kn1HkFs1ξ,η,z^kJξ,ηdξdηu^′b∗≅∑i,j=1n2HijFb2ξi,ηj+∑i,j=1n2HijFs2ξi,ηju^′b∗,where Fm1ξ,η,z^, Fm2ξ,η, Fb1ξ,η,z^, Fb2ξ,η, Fs1ξ,η,z^, and Fs2ξ,η are given by equation (34), Jξ,η is the Jacobian determinant, n1 is the total integral points along the direction of element thickness, and n2 is the total integral points along the direction of each natural coordinate axis (i.e., square element coordinate axis) in the element reference plane.(34)Fm1ξ,η,z^=Bm∗ξ,ηTDbξ,η,z^Bm∗ξ,η,Fm2ξ,η=∑k=1n1HkFm1ξ,η,z^kJξ,η,Fb1ξ,η,z^=Bb∗ξ,ηTDbξ,η,z^Bb∗ξ,η,Fb2ξ,η=∑k=1n1Hkz^k2Fb1ξ,η,z^kJξ,η,Fs1ξ,η,z^=Bs∗ξ,ηTDsξ,η,z^Bs∗ξ,η,Fs2ξ,η=∑k=1n1HkFs1ξ,η,z^kJξ,η.


## 3.1. Coordinate Definition Modes

The coordinate definition mode I refers to the coordinates of the element nodes when the pure deformation linear (angular) displacement and internal force (moment) of shell element nodes are solved from timet0 to t, so they need to be on the same reference plane. The four-node shell element is formed into a spatial quadrilateral after undergoing motion and deformation. The projection coordinates (x¯i,x¯i′ and θ¯i,θ¯i′) for the element nodes on the reference planes which are perpendicular to the mean normal vector (n,n′) of the element and passing through the centroid points (C,C′) of the original element can be used as the coordinate definition mode I to realize the transformation of the spatial quadrilateral to the planar quadrilateral.Coordinate definition mode II refers to the coordinates of each particle when the central difference formulas are applied for solving the particle motion equation (1), so it needs to be a unique value. The following two methods can be adopted: (1) the actual coordinates of the particle corresponding to the spatial quadrilateral after the element movement and deformation; (2) I element nodes corresponding to the same single particle in coordinate definition mode I obtained by plane projection do not coincide with each other, so the average coordinate value (x=1/k∑i=1kx¯i,x′=1/k∑i=1kx¯i′ and θ=1/k∑i=1kθ¯i,θ′=1/k∑i=1kθ¯i′) of each element nodes corresponding to the same single particle in coordinate definition mode I can be taken as coordinate definition mode II, where k is the total number for that corresponding element nodes connected to each element by the same single particle.


## 3.2. Integral of Internal Forces for Element Nodes

When solving the node internal forces and node internal moments for the membrane element part and plate element part of the quadrilateral isoparametric shell element according to equation (31), the numerical integration operation will be involved since the internal stresses and strains of the elements change with their positions. Considering the case of general elastic-plastic material, the integrations along the thickness and element plane are solved by Newton–Cotes integral scheme and 2D Gaussian integral scheme in the traditional finite element method, respectively, [21]. The numerical integration formulas of Δf^m∗ and Δf^′θb∗ are given by the following equation:(33)Δf^m∗=∫A^a∫−ta/2ta/2Bm∗TDbBm∗dz^dA^au^m∗=∫−11∫−11∫−ta/2ta/2Bm∗ξ,ηTDbξ,η,z^Bm∗ξ,ηdz^Jξ,ηdξdηu^m∗=∫−11∫−11∑k=1n1HkFm1ξ,η,z^kJξ,ηdξdηu^m∗≅∑i,j=1n2HijFm2ξi,ηju^m∗,Δf^′θb∗=∫A^a∫−ta/2ta/2z^2Bb∗TDbBb∗dz^dA^a+∫A^a∫−ta/2ta/2Bs∗TDsBs∗dz^dA^au^′b∗=∫−11∫−11∫−ta/2ta/2z^2Bb∗ξ,ηTDbξ,η,z^Bb∗ξ,ηdz^Jξ,ηdξdη+∫−11∫−11∫−ta/2ta/2Bs∗ξ,ηTDsξ,η,z^Bs∗ξ,ηdz^Jξ,ηdξdηu^′b∗=∫−11∫−11∑kn1Hkz^k2Fb1ξ,η,z^kJξ,ηdξdη+∫−11∫−11∑kn1HkFs1ξ,η,z^kJξ,ηdξdηu^′b∗≅∑i,j=1n2HijFb2ξi,ηj+∑i,j=1n2HijFs2ξi,ηju^′b∗,where Fm1ξ,η,z^, Fm2ξ,η, Fb1ξ,η,z^, Fb2ξ,η, Fs1ξ,η,z^, and Fs2ξ,η are given by equation (34), Jξ,η is the Jacobian determinant, n1 is the total integral points along the direction of element thickness, and n2 is the total integral points along the direction of each natural coordinate axis (i.e., square element coordinate axis) in the element reference plane.(34)Fm1ξ,η,z^=Bm∗ξ,ηTDbξ,η,z^Bm∗ξ,η,Fm2ξ,η=∑k=1n1HkFm1ξ,η,z^kJξ,η,Fb1ξ,η,z^=Bb∗ξ,ηTDbξ,η,z^Bb∗ξ,η,Fb2ξ,η=∑k=1n1Hkz^k2Fb1ξ,η,z^kJξ,η,Fs1ξ,η,z^=Bs∗ξ,ηTDsξ,η,z^Bs∗ξ,η,Fs2ξ,η=∑k=1n1HkFs1ξ,η,z^kJξ,η.


## 4. Realization of Nonlinear Material Constitutive

The bilinear elastoplastic constitutive model is widely used in nonlinear analysis of structures considering the linear plastic hardening effect, which is used in this paper for theoretical derivation and numerical analysis.


### 4.1. Bilinear Elastoplastic Constitutive Model

Considering the case of isotropic material constitution, the dynamic yield stress is obtained by scaling static yield stress through plastic factor, and its constitutive equation is given by the following equation:(35)σd=σ0+Epε¯p,where Ep is the plastic hardening modulus, Ep=EtE/E−Et, E is the Young modulus, Et is the tangent modulus, σ0 is the static yield stress, ε¯p=2εij,pεij,p/3 is the equivalent plastic strain, εij,p is the plastic strain tensor.This constitutive model is suitable for bilinear elastic-plastic materials. The case ofEt=0 corresponds to the ideal elastic-plastic model, and the case of σ0≫σ corresponds to the linear elastic model.The von Mises yield criterion is adopted in this paper as a dynamic yield criterion, that is, the plastic stress state (i.e., yield condition equation) is given by the following equation:(36)ϕ=σ¯−σd=0,where σ¯ is the equivalent von Mises stress, which is given by equation (37), σd is the dynamic yield stress.(37)σ¯=32sijsij=32σij−δijσmσij−δijσm,where σm=σ11+σ22+σ33/3 is the mean principal stress, σij is the stress component, sij is the deviatoric stress component, δij=1i≠j, δij=0i=j.For shell element, the principal stress along the thickness direction of elementσ33 is always zero, that is σ33=0. Accordingly, the transverse principal strain along the thickness direction of element ε33 is an dependent variable, which is given by the following equation:(38)ε33=−υ1−υε11+ε22.In process of solving the elastic-plastic incremental matrix constitutiveDep, the 6 × 6 elastic constitutive matrix De containing all stress-strain components is used for theoretical derivation. The elastic stress-strain relationship used for trial calculation is given by the following equation:(39)Δσ^=DeΔε^,where De is given by equation (40), Δσ^=Δσ^11Δσ^22Δσ^33Δτ^12Δτ^23Δτ^31T is the stress increment vector, Δε^=Δε^11Δε^22Δε^33Δγ^12Δγ^23Δγ^31T is the strain increment vector.(40)D=E1−υ21υ0000υ100001−υυ1−2υ1−υυ1−2υ1−υ21−2υ0000001−υ2000000k1−υ2000000k1−υ2,where E is Young modulus, k is the shear correction coefficient, and υ is Poisson's ratio, which are given in equation (29).


### 4.2. Plastic Incremental Analysis Steps

Given the initial stress and strain at the beginning of the time substep and the incremental strain during the time substep, the elastic-plastic incremental analysis is used to obtain the updated stresses and strains under the condition for satisfying the material constitutive and yield criteria at the end of time substeps. Detailed steps are as follows:Step 1: elastic prediction. The stress componentσij,t and the equivalent plastic strain ε¯p,t at the time t (i.e., at the beginning of time substep) and the incremental strain component Δε^ij from the time t to the time t+Δt (i.e., during the time substep) are given. Assuming that the time substep is completely elastic, the predicted elastic stress matrix σtr,t+Δt and the predicted plastic hardening parameter matrix κtr,t+Δt=Ep,t+Δtε¯p at the time t+Δt are obtained.Step 2: yield judgment. Substituteσtr,t+Δt and κtr,t+Δt at the time t+Δt into the yield condition equation (36), and judge whether it is an elastic or plastic incremental time substep according to the positive or negative of the predicted yield function ϕtr,t+Δt. If ϕtr,t+Δt>0, it is a plastic time substep; Otherwise, it is an elastic time substep.Step 3: plastic correction. If it is an elastic incremental time substep, then skip this step. If it is a plastic incremental time substep, the plastic flow factor incrementΔλ is calculated on the basis of normal flow rule and dynamic yield criterion. The case of Δλ>0 corresponds to the plastic incremental time substep, and the case of Δλ=0 corresponds to the elastic incremental time substep.Step 4: stress and strain updating. If it is an elastic incremental time substep, thenσt+Δt=σtr,t+Δt. If it is a plastic incremental time substep, the stress and strain updating is carried out by Δλ. The updated variables include the equivalent plastic strain ε¯p,t+Δt, the dynamic yield stress σd,t+Δt, the deviatoric stress component sij,t+Δt, the von Mises stress σ¯t+Δt and the actual stress σt+Δt at the time t+Δt.Step 5: solution of the elastic-plastic matrix. If it is an elastic time substep, there is no plastic matrixDp. If it is a plastic time substep, the elastic-plastic matrix Dep needs to be solved. Through dynamic yield criterion and normal flow rule, the plastic flow factor increment dλ and plastic matrix Dp are obtained, and then Dep=De−Dp. Then, the internal forces (moments) of element nodes are solved by integral formulas.


## 4.1. Bilinear Elastoplastic Constitutive Model

Considering the case of isotropic material constitution, the dynamic yield stress is obtained by scaling static yield stress through plastic factor, and its constitutive equation is given by the following equation:(35)σd=σ0+Epε¯p,where Ep is the plastic hardening modulus, Ep=EtE/E−Et, E is the Young modulus, Et is the tangent modulus, σ0 is the static yield stress, ε¯p=2εij,pεij,p/3 is the equivalent plastic strain, εij,p is the plastic strain tensor.This constitutive model is suitable for bilinear elastic-plastic materials. The case ofEt=0 corresponds to the ideal elastic-plastic model, and the case of σ0≫σ corresponds to the linear elastic model.The von Mises yield criterion is adopted in this paper as a dynamic yield criterion, that is, the plastic stress state (i.e., yield condition equation) is given by the following equation:(36)ϕ=σ¯−σd=0,where σ¯ is the equivalent von Mises stress, which is given by equation (37), σd is the dynamic yield stress.(37)σ¯=32sijsij=32σij−δijσmσij−δijσm,where σm=σ11+σ22+σ33/3 is the mean principal stress, σij is the stress component, sij is the deviatoric stress component, δij=1i≠j, δij=0i=j.For shell element, the principal stress along the thickness direction of elementσ33 is always zero, that is σ33=0. Accordingly, the transverse principal strain along the thickness direction of element ε33 is an dependent variable, which is given by the following equation:(38)ε33=−υ1−υε11+ε22.In process of solving the elastic-plastic incremental matrix constitutiveDep, the 6 × 6 elastic constitutive matrix De containing all stress-strain components is used for theoretical derivation. The elastic stress-strain relationship used for trial calculation is given by the following equation:(39)Δσ^=DeΔε^,where De is given by equation (40), Δσ^=Δσ^11Δσ^22Δσ^33Δτ^12Δτ^23Δτ^31T is the stress increment vector, Δε^=Δε^11Δε^22Δε^33Δγ^12Δγ^23Δγ^31T is the strain increment vector.(40)D=E1−υ21υ0000υ100001−υυ1−2υ1−υυ1−2υ1−υ21−2υ0000001−υ2000000k1−υ2000000k1−υ2,where E is Young modulus, k is the shear correction coefficient, and υ is Poisson's ratio, which are given in equation (29).


## 4.2. Plastic Incremental Analysis Steps

Given the initial stress and strain at the beginning of the time substep and the incremental strain during the time substep, the elastic-plastic incremental analysis is used to obtain the updated stresses and strains under the condition for satisfying the material constitutive and yield criteria at the end of time substeps. Detailed steps are as follows:Step 1: elastic prediction. The stress componentσij,t and the equivalent plastic strain ε¯p,t at the time t (i.e., at the beginning of time substep) and the incremental strain component Δε^ij from the time t to the time t+Δt (i.e., during the time substep) are given. Assuming that the time substep is completely elastic, the predicted elastic stress matrix σtr,t+Δt and the predicted plastic hardening parameter matrix κtr,t+Δt=Ep,t+Δtε¯p at the time t+Δt are obtained.Step 2: yield judgment. Substituteσtr,t+Δt and κtr,t+Δt at the time t+Δt into the yield condition equation (36), and judge whether it is an elastic or plastic incremental time substep according to the positive or negative of the predicted yield function ϕtr,t+Δt. If ϕtr,t+Δt>0, it is a plastic time substep; Otherwise, it is an elastic time substep.Step 3: plastic correction. If it is an elastic incremental time substep, then skip this step. If it is a plastic incremental time substep, the plastic flow factor incrementΔλ is calculated on the basis of normal flow rule and dynamic yield criterion. The case of Δλ>0 corresponds to the plastic incremental time substep, and the case of Δλ=0 corresponds to the elastic incremental time substep.Step 4: stress and strain updating. If it is an elastic incremental time substep, thenσt+Δt=σtr,t+Δt. If it is a plastic incremental time substep, the stress and strain updating is carried out by Δλ. The updated variables include the equivalent plastic strain ε¯p,t+Δt, the dynamic yield stress σd,t+Δt, the deviatoric stress component sij,t+Δt, the von Mises stress σ¯t+Δt and the actual stress σt+Δt at the time t+Δt.Step 5: solution of the elastic-plastic matrix. If it is an elastic time substep, there is no plastic matrixDp. If it is a plastic time substep, the elastic-plastic matrix Dep needs to be solved. Through dynamic yield criterion and normal flow rule, the plastic flow factor increment dλ and plastic matrix Dp are obtained, and then Dep=De−Dp. Then, the internal forces (moments) of element nodes are solved by integral formulas.


## 5. Programming and Example Verification

Based on the above-given theoretical formulas of vector quadrilateral isoparametric shell element with shear deformation, a computational analysis program is developed by Matlab software to realize the calculation and analysis of shell structure with shear deformation. The program implementation process is shown in Figure1. Then, the rationality and practicability of derived formulas and program codes are verified by the following examples.Figure 1
Analysis flow chart for shell structures.


### 5.1. Static Bending of Annular Thick-Plate (Shear Effect along Thickness)

Taking an annular thick-plate with circular line load as an example [22], the influence of the shear effect along thickness for shell element on small bending deformation of thick plate structure is studied. Figure 2 shows the typical annular thick-plate structure, including the geometric model and mesh model. The inner radius of plate is R1 = 1.4 in, the outer radius of the plate is R2 = 2.0 in and the thickness of plate is ta = 0.5 in (corresponding width-thickness ratio is (R2 — R1)/ta = 1.2 < 5, that is, thick plate). The plate is placed horizontally and initially at rest. The boundary conditions for the inner circular edge and outer circular edge are respectively simply supported and free. The line loading q = 800 lb/in is loaded on the circular line at the radius R = 1.8 in of the annular thick-plate, and the static results are obtained by applying damping to eliminate the dynamic oscillation effect. The mass density is ρ = 0.1 lbs/in3, the Young modulus is E = 1.8 × 107 lb/in2 and the ’Poisson’s ratio is υ = 0.3. The quadrilateral isoparametric shell element is used for numerical simulation. The time substep is h = 6.0 × 10−6 s and the damping parameter α = 2 × 104.Figure 2
Annular thick-plate model. (a) Geometric model. (b) Mesh model.
(a)(b)In order to investigate the influence of different mesh sizes on the VFIFE method in this paper, the mesh sizesl are set as 0.05 in, 0.10 in, and 0.20 in, respectively, for solving calculation. Figure 3 shows the comparison of vertical displacement convergence values of the outer-edge nodes with different mesh sizes. It can be seen that with the decrease of element mesh size, the vertical displacement decreases and the decreasing speed becomes slow. Therefore, the element mesh size l of 0.10 in is selected for in-depth analysis in the following, which has a relatively small calculation scale and high calculation accuracy.Figure 3
Comparison of vertical displacements of outer-edge nodes with different mesh sizes (ta = 0.5 in).In order to investigate the influence of different thicknesses on the VFIFE method in this paper, the shell thicknessesta are set as 0.1∼1.5 in (the corresponding length-to-thickness L/ta is 6∼0.4), respectively, for solving calculation. Figure 4 shows the comparison of vertical displacement convergence values of the outer-edge nodes with different shell thicknesses. It can be seen that when ta ≥ 0.5 in (that is ta/l ≥ 5, which belongs to medium-thick shell, corresponding to L/ta = 1.2), the error between the VFIFE method and ANSYS calculation results is relatively small, and the calculation accuracy of the VFIFE method in this paper is improved with the increase of shell thickness. When ta < 0.5 (that is, belonging to thin-shell), the error between the VFIFE method and ANSYS calculation results is relatively large. When ta ≤ 0.2, the thick-shell element superimposed by the VFIFE method is not applicable, so there is no corresponding calculated value. Therefore, for the medium-thick shell element with ta/l ≥ 5, the VFIFE method has a high accuracy.Figure 4
Comparison of vertical displacements of outer-edge nodes with different thicknesses (mesh size is 0.10 in).Figure5 shows the computational convergence process for the vertical displacements of outer-edge nodes of the annular thick-plate with the mesh size l of 0.10 in. The result shows that it has good convergence and the vertical displacement convergence result is about 0.00533 in. Figure 6 shows the calculated results of the vertical displacements and von Mises stresses for radial nodes along the radial direction (corresponding to radial distance) of the annular thick-plate, respectively, and the comparison with ANSYS calculated results including that with shear effect (i.e., ANSYS-shear) and without shear effect (i.e., ANSYS-no-shear). Among them, I results of the cases for ANSYS-shear and ANSYS-non-shear are obtained by applying the shell43 element (considering Mindlin–Ressner first-order shear effect) and shell63 element (considering no shear effect), respectively.Figure 5
Convergence process of vertical displacement for outer-edge node.Figure 6
Results of the radial nodes along the radial direction. (a) Vertical displacement. (b) Von mises stress (SNEG surface).
(a)(b)The result from Figure6 shows that the vertical displacements of the outer circular edge node of the annular thick-plate calculated by the method of this paper (i.e., VFIFE), ANSYS-shear, and ANSYS-no-shear are 0.00533 in, 0.00534 in and 0.00521 in respectively. The errors of the displacement results of VFIFE are about 0.18% compared with the results of ANSYS-shear and about 2.31% compared with the results of ANSYS-no-shear considering the influence of shear deformation. The calculated result is basically consistent with the result in literature [22]. The case of results for von Mises stress is similar, that is, the maximum errors of the von Mises stress results of VFIFE are about 0.37% compared with the results of ANSYS-shear and about 2.58% compared with the results of ANSYS-no-shear. The calculated results of the proposed method in this paper (i.e., VFIFE method) are close to those of ANSYS-shear, which verifies the correctness and effectiveness of the VFIFE method. For the thick-plate structure, the shear effect has a certain influence, and the VFIFE method presented here has good accuracy for analyzing the thick-plate structure with shear deformation. With the increase of plate thickness, the shear effect cannot be ignored gradually.The quadrilateral shell element of the VFIFE method in this paper is based on the linear superposition of the quadrilateral isoparametric membrane element and quadrilateral Mindlin isoparametric plate element, which is mainly applied for medium-thick shell structures. For thin shell structures, in addition to the low calculation accuracy, there is also a shear self-locking phenomenon. Due to the complexity of this problem, this paper did not make a detailed analysis and considered further research.


### 5.2. Dynamic Bending of Long Cantilever Plate (Large Displacement and Rotation)

Taking a long cantilevered plate with a concentrated load at its cantilevered end as an example [23]. Figure 7 shows the typical cantilevered plate structure, including the geometric model and mesh model. The width of the plate is B = 0.2 m, the length of plate is L = 1.0 m and the thickness of plate is ta = 0.015 m. The plate is placed horizontally and initial at rest. The boundary conditions for one end and other end are, respectively, fixed and free cantilever. The concentrated load P = 100 kN is loaded on the midpoint of the cantilevered end for the plate. The ramp-platform dynamic loading mode is adopted, and the convergence value is achieved through damping. The mass density is ρ = 7850 kg/m3, the Young modulus is E = 201 GPa and the ’Poisson’s ratio is υ = 0.3. The quadrilateral isoparametric shell element is used for numerical simulation. The time substep is h = 2.0 × 10−6 s and the damping parameter α = 200.Figure 7
Long cantilevered plate model. (a) Geometric model. (b) Mesh model.
(a)(b)The four cases of quadrilateral shell element (i.e., the superposition of quadrilateral membrane elements and Mindlin plate elements) in this paper with shear correction coefficientsk = 5/6, 2/3, 1/3, 0 and the case of triangular shell element in literature [23] (i.e., the superposition of CST membrane element and DKT plate element) without shear deformation are considered for analysis and comparison. Figure 8 shows the vertical displacement-time curves for a long cantilevered plate at the central node of the free end, which shows good convergence under quasi-static calculation.Figure 8
Vertical displacements-time curves at the central node of the free end.It can be observed from Figure8 that the vertical displacement convergence values of the free end central node for the long cantilevered plate calculated by the VFIFE method of this paper with shear correction coefficients k = 5/6, 2/3, 1/3, and 0 are about 0.186 m, 0.227 m, 0.369 m, and 0.983 m, respectively, and that of literature [23] without shear deformation is 0.800 m. The vertical displacement at the central node for free end with k = 5/6 (i.e., general shear) is only about 19% and 23% of that with k = 0 (i.e., no shear) and that of literature [23], respectively. For long cantilever shell structures, the shear effect of shell element is extremely significant and cannot be ignored. In this case, the triangular shell element without shear deformation given in literature [23] is no longer applicable, and the quadrilateral isoparametric shell element with shear deformation given in this paper should be adopted.In addition, the vertical displacement-time variation curves of the quadrilateral shell element withk = 0 (i.e., no shear) and the nonshear triangle shell element given in reference [23] are basically consistent before the time of 0.02 s, which also verifies the practicability of the proposed method (VFIFE) in a certain extent.Figure9 shows the displacement deformation diagrams of the long cantilever plate with shear correction coefficients k = 5/6, 2/3, 1/3, and 0 at the time of 0.04 s. It can be observed from Figure 9 that the long cantilever plate with a shear deformation effect presents a bending-shear bending form, which is a process of large displacement and large rotation. Taking the case of k = 5/6 as an example, Figure 10 shows the von Mises stress contours (MATLAB self-programming is used to draw) of the long cantilever plate at the time of 0.04 s. It turned out from Figure 10 that maximum stress appears in the root region of the long cantilever plate.Figure 9
Displacement deformation diagrams. (a)k = 5/6 (b) k = 2/3 (c) k = 1/3 (d) k = 0.
(a)(b)(c)(d)Figure 10
Von mises stress contours of the long cantilevered plate at the time of 0.04 s.In the above, Figure9 is not a curve, but a side view of the long cantilever plate with large displacement and rotation, which reflects the applicability of the VFIFE method in this paper under the condition of large deformation and rotation. As for the quantitative analysis of deformation under different shear coefficients, it has been described in Figure 8 and above.


### 5.3. Nonlinear Stability Analysis of Shallow Spherical Shell (Elastic-Plastic Material Constitutive)

Taking a jump instability of the shallow spherical shell with the central concentrated load as an example [24], the influence of elastic-plastic material constitutive for shell element on nonlinear stability is studied. Figure 11 shows the typical shallow spherical shell structure, including the geometric model and mesh model. The spherical radius of shell is R = 2.54 m, the projected side length of shell is L = 1.5698 m and the thickness of shell is ta = 99.45 m. The shell is placed horizontally and initial at rest. The boundary conditions are that four sides are simply supported. The concentrated load P is loaded on the central node of shell. The mass density is ρ = 2500 kg/m3, the Young modulus is E = 68.95 MPa and the Poisson’s ratio is υ = 0.3. The quadrilateral isoparametric shell element is used for numerical simulation. The time substep is h = 2.0 × 10−4 s and the damping parameter α = 20. The whole process of elastic-plastic instability of the shallow spherical shell is tracked by the displacement control method.Figure 11
Shallow spherical shell model. (a) Geometric model. (b) Mesh model.
(a)(b)Figure12 shows the concentrated load-vertical displacement curves for the central node of the shallow spherical shell. It can be observed that with the increases of the vertical displacement for the central nodes, there is an obvious load extreme point of the shallow spherical shell. After crossing the extreme point, the vertical displacement deformation turns over greatly and the load drops sharply. At this time, the structure cannot be used normally due to the excessive deformation, that is, the jump buckling occurs. After complete overturning, the load continues to increase with the increase of structural vertical displacement.Figure 12
Load-displacement curves for center nodes of the shallow spherical shells.There are some differences in the load-displacement process between the VFIFE method of this paper and the ABAQUS method, but the overall changing trend is consistent, and the load extreme points are about the same. The maximum extreme point loading from the VFIFE method is about 16.45 kN when the vertical displacement is about 210 mm, and that from the ABAQUS method is about 18.76 kN, with an error of about 12.3%, which is relatively small. The calculation errors mainly come from the error of solving the jump instability problem and the number of shell elements. The error can be further reduced with the increase of the number of shell elements. As for the large difference between VFIFE and ABAQUS results in Figure12, it may be related to the loading path and displacement accumulation effect of the buckling problem involving large deformation and abrupt deformation. However, the buckling extreme points are basically the same, which verifies the effectiveness of this method to some extent.Figure13 shows the von Mises stress contours (MATLAB self-programming is used to draw) of the shallow spherical shell when the maximum extreme point load (corresponding to the vertical displacement 210 mm) is reached by the VFIFE method. It can be observed from Figure 13 that the region near the central node is in a state of plastic hardening due to the concentrated load.Figure 13
Von mises stress contours of the shallow spherical shell.


## 5.1. Static Bending of Annular Thick-Plate (Shear Effect along Thickness)

Taking an annular thick-plate with circular line load as an example [22], the influence of the shear effect along thickness for shell element on small bending deformation of thick plate structure is studied. Figure 2 shows the typical annular thick-plate structure, including the geometric model and mesh model. The inner radius of plate is R1 = 1.4 in, the outer radius of the plate is R2 = 2.0 in and the thickness of plate is ta = 0.5 in (corresponding width-thickness ratio is (R2 — R1)/ta = 1.2 < 5, that is, thick plate). The plate is placed horizontally and initially at rest. The boundary conditions for the inner circular edge and outer circular edge are respectively simply supported and free. The line loading q = 800 lb/in is loaded on the circular line at the radius R = 1.8 in of the annular thick-plate, and the static results are obtained by applying damping to eliminate the dynamic oscillation effect. The mass density is ρ = 0.1 lbs/in3, the Young modulus is E = 1.8 × 107 lb/in2 and the ’Poisson’s ratio is υ = 0.3. The quadrilateral isoparametric shell element is used for numerical simulation. The time substep is h = 6.0 × 10−6 s and the damping parameter α = 2 × 104.Figure 2
Annular thick-plate model. (a) Geometric model. (b) Mesh model.
(a)(b)In order to investigate the influence of different mesh sizes on the VFIFE method in this paper, the mesh sizesl are set as 0.05 in, 0.10 in, and 0.20 in, respectively, for solving calculation. Figure 3 shows the comparison of vertical displacement convergence values of the outer-edge nodes with different mesh sizes. It can be seen that with the decrease of element mesh size, the vertical displacement decreases and the decreasing speed becomes slow. Therefore, the element mesh size l of 0.10 in is selected for in-depth analysis in the following, which has a relatively small calculation scale and high calculation accuracy.Figure 3
Comparison of vertical displacements of outer-edge nodes with different mesh sizes (ta = 0.5 in).In order to investigate the influence of different thicknesses on the VFIFE method in this paper, the shell thicknessesta are set as 0.1∼1.5 in (the corresponding length-to-thickness L/ta is 6∼0.4), respectively, for solving calculation. Figure 4 shows the comparison of vertical displacement convergence values of the outer-edge nodes with different shell thicknesses. It can be seen that when ta ≥ 0.5 in (that is ta/l ≥ 5, which belongs to medium-thick shell, corresponding to L/ta = 1.2), the error between the VFIFE method and ANSYS calculation results is relatively small, and the calculation accuracy of the VFIFE method in this paper is improved with the increase of shell thickness. When ta < 0.5 (that is, belonging to thin-shell), the error between the VFIFE method and ANSYS calculation results is relatively large. When ta ≤ 0.2, the thick-shell element superimposed by the VFIFE method is not applicable, so there is no corresponding calculated value. Therefore, for the medium-thick shell element with ta/l ≥ 5, the VFIFE method has a high accuracy.Figure 4
Comparison of vertical displacements of outer-edge nodes with different thicknesses (mesh size is 0.10 in).Figure5 shows the computational convergence process for the vertical displacements of outer-edge nodes of the annular thick-plate with the mesh size l of 0.10 in. The result shows that it has good convergence and the vertical displacement convergence result is about 0.00533 in. Figure 6 shows the calculated results of the vertical displacements and von Mises stresses for radial nodes along the radial direction (corresponding to radial distance) of the annular thick-plate, respectively, and the comparison with ANSYS calculated results including that with shear effect (i.e., ANSYS-shear) and without shear effect (i.e., ANSYS-no-shear). Among them, I results of the cases for ANSYS-shear and ANSYS-non-shear are obtained by applying the shell43 element (considering Mindlin–Ressner first-order shear effect) and shell63 element (considering no shear effect), respectively.Figure 5
Convergence process of vertical displacement for outer-edge node.Figure 6
Results of the radial nodes along the radial direction. (a) Vertical displacement. (b) Von mises stress (SNEG surface).
(a)(b)The result from Figure6 shows that the vertical displacements of the outer circular edge node of the annular thick-plate calculated by the method of this paper (i.e., VFIFE), ANSYS-shear, and ANSYS-no-shear are 0.00533 in, 0.00534 in and 0.00521 in respectively. The errors of the displacement results of VFIFE are about 0.18% compared with the results of ANSYS-shear and about 2.31% compared with the results of ANSYS-no-shear considering the influence of shear deformation. The calculated result is basically consistent with the result in literature [22]. The case of results for von Mises stress is similar, that is, the maximum errors of the von Mises stress results of VFIFE are about 0.37% compared with the results of ANSYS-shear and about 2.58% compared with the results of ANSYS-no-shear. The calculated results of the proposed method in this paper (i.e., VFIFE method) are close to those of ANSYS-shear, which verifies the correctness and effectiveness of the VFIFE method. For the thick-plate structure, the shear effect has a certain influence, and the VFIFE method presented here has good accuracy for analyzing the thick-plate structure with shear deformation. With the increase of plate thickness, the shear effect cannot be ignored gradually.The quadrilateral shell element of the VFIFE method in this paper is based on the linear superposition of the quadrilateral isoparametric membrane element and quadrilateral Mindlin isoparametric plate element, which is mainly applied for medium-thick shell structures. For thin shell structures, in addition to the low calculation accuracy, there is also a shear self-locking phenomenon. Due to the complexity of this problem, this paper did not make a detailed analysis and considered further research.


## 5.2. Dynamic Bending of Long Cantilever Plate (Large Displacement and Rotation)

Taking a long cantilevered plate with a concentrated load at its cantilevered end as an example [23]. Figure 7 shows the typical cantilevered plate structure, including the geometric model and mesh model. The width of the plate is B = 0.2 m, the length of plate is L = 1.0 m and the thickness of plate is ta = 0.015 m. The plate is placed horizontally and initial at rest. The boundary conditions for one end and other end are, respectively, fixed and free cantilever. The concentrated load P = 100 kN is loaded on the midpoint of the cantilevered end for the plate. The ramp-platform dynamic loading mode is adopted, and the convergence value is achieved through damping. The mass density is ρ = 7850 kg/m3, the Young modulus is E = 201 GPa and the ’Poisson’s ratio is υ = 0.3. The quadrilateral isoparametric shell element is used for numerical simulation. The time substep is h = 2.0 × 10−6 s and the damping parameter α = 200.Figure 7
Long cantilevered plate model. (a) Geometric model. (b) Mesh model.
(a)(b)The four cases of quadrilateral shell element (i.e., the superposition of quadrilateral membrane elements and Mindlin plate elements) in this paper with shear correction coefficientsk = 5/6, 2/3, 1/3, 0 and the case of triangular shell element in literature [23] (i.e., the superposition of CST membrane element and DKT plate element) without shear deformation are considered for analysis and comparison. Figure 8 shows the vertical displacement-time curves for a long cantilevered plate at the central node of the free end, which shows good convergence under quasi-static calculation.Figure 8
Vertical displacements-time curves at the central node of the free end.It can be observed from Figure8 that the vertical displacement convergence values of the free end central node for the long cantilevered plate calculated by the VFIFE method of this paper with shear correction coefficients k = 5/6, 2/3, 1/3, and 0 are about 0.186 m, 0.227 m, 0.369 m, and 0.983 m, respectively, and that of literature [23] without shear deformation is 0.800 m. The vertical displacement at the central node for free end with k = 5/6 (i.e., general shear) is only about 19% and 23% of that with k = 0 (i.e., no shear) and that of literature [23], respectively. For long cantilever shell structures, the shear effect of shell element is extremely significant and cannot be ignored. In this case, the triangular shell element without shear deformation given in literature [23] is no longer applicable, and the quadrilateral isoparametric shell element with shear deformation given in this paper should be adopted.In addition, the vertical displacement-time variation curves of the quadrilateral shell element withk = 0 (i.e., no shear) and the nonshear triangle shell element given in reference [23] are basically consistent before the time of 0.02 s, which also verifies the practicability of the proposed method (VFIFE) in a certain extent.Figure9 shows the displacement deformation diagrams of the long cantilever plate with shear correction coefficients k = 5/6, 2/3, 1/3, and 0 at the time of 0.04 s. It can be observed from Figure 9 that the long cantilever plate with a shear deformation effect presents a bending-shear bending form, which is a process of large displacement and large rotation. Taking the case of k = 5/6 as an example, Figure 10 shows the von Mises stress contours (MATLAB self-programming is used to draw) of the long cantilever plate at the time of 0.04 s. It turned out from Figure 10 that maximum stress appears in the root region of the long cantilever plate.Figure 9
Displacement deformation diagrams. (a)k = 5/6 (b) k = 2/3 (c) k = 1/3 (d) k = 0.
(a)(b)(c)(d)Figure 10
Von mises stress contours of the long cantilevered plate at the time of 0.04 s.In the above, Figure9 is not a curve, but a side view of the long cantilever plate with large displacement and rotation, which reflects the applicability of the VFIFE method in this paper under the condition of large deformation and rotation. As for the quantitative analysis of deformation under different shear coefficients, it has been described in Figure 8 and above.


## 5.3. Nonlinear Stability Analysis of Shallow Spherical Shell (Elastic-Plastic Material Constitutive)

Taking a jump instability of the shallow spherical shell with the central concentrated load as an example [24], the influence of elastic-plastic material constitutive for shell element on nonlinear stability is studied. Figure 11 shows the typical shallow spherical shell structure, including the geometric model and mesh model. The spherical radius of shell is R = 2.54 m, the projected side length of shell is L = 1.5698 m and the thickness of shell is ta = 99.45 m. The shell is placed horizontally and initial at rest. The boundary conditions are that four sides are simply supported. The concentrated load P is loaded on the central node of shell. The mass density is ρ = 2500 kg/m3, the Young modulus is E = 68.95 MPa and the Poisson’s ratio is υ = 0.3. The quadrilateral isoparametric shell element is used for numerical simulation. The time substep is h = 2.0 × 10−4 s and the damping parameter α = 20. The whole process of elastic-plastic instability of the shallow spherical shell is tracked by the displacement control method.Figure 11
Shallow spherical shell model. (a) Geometric model. (b) Mesh model.
(a)(b)Figure12 shows the concentrated load-vertical displacement curves for the central node of the shallow spherical shell. It can be observed that with the increases of the vertical displacement for the central nodes, there is an obvious load extreme point of the shallow spherical shell. After crossing the extreme point, the vertical displacement deformation turns over greatly and the load drops sharply. At this time, the structure cannot be used normally due to the excessive deformation, that is, the jump buckling occurs. After complete overturning, the load continues to increase with the increase of structural vertical displacement.Figure 12
Load-displacement curves for center nodes of the shallow spherical shells.There are some differences in the load-displacement process between the VFIFE method of this paper and the ABAQUS method, but the overall changing trend is consistent, and the load extreme points are about the same. The maximum extreme point loading from the VFIFE method is about 16.45 kN when the vertical displacement is about 210 mm, and that from the ABAQUS method is about 18.76 kN, with an error of about 12.3%, which is relatively small. The calculation errors mainly come from the error of solving the jump instability problem and the number of shell elements. The error can be further reduced with the increase of the number of shell elements. As for the large difference between VFIFE and ABAQUS results in Figure12, it may be related to the loading path and displacement accumulation effect of the buckling problem involving large deformation and abrupt deformation. However, the buckling extreme points are basically the same, which verifies the effectiveness of this method to some extent.Figure13 shows the von Mises stress contours (MATLAB self-programming is used to draw) of the shallow spherical shell when the maximum extreme point load (corresponding to the vertical displacement 210 mm) is reached by the VFIFE method. It can be observed from Figure 13 that the region near the central node is in a state of plastic hardening due to the concentrated load.Figure 13
Von mises stress contours of the shallow spherical shell.


## 6. Conclusions

In this paper, the theory and application for the quadrilateral isoparametric shell elements based on the VFIFE method are researched, and the following conclusions are obtained:(1)
Based on the linear superposition of the quadrilateral isoparametric membrane elements and the quadrilateral Mindlin plate elements, the basic vector formula of the quadrilateral isoparametric shell element with shear deformation is derived. Two different coordinate modes for updating the particle displacement and the element node internal force are proposed. The uniform solution of the irregular quadrilateral element is realized by the isoparametric integral scheme along the element plane.(2)
The bilinear elastoplastic constitutive model for the material is introduced into vector finite shell elements, and the nonlinear effect of the material is realized based on the integral method along the thickness of the element.(3)
Based on the established basic theory of quadrilateral isoparametric shell elements with vector finite element, the corresponding calculation, and analysis program is developed. Example analysis shows that the vector finite element program can obtain good computational accuracy in the static and dynamic analysis, large deformation, and large rotation analysis for shell structure, which verifies the rationality and practicability for derived formulas and program codes.(4)
VFIFE method is applied to the stability analysis for shell structures with double nonlinearity including geometric nonlinearity and material nonlinearity. By means of displacement control, the jump instability process of plate and shell structure can be obtained without special treatment, and the critical loading can be obtained, which reflects the advantages of this method in analyzing the stability problems of material with large deformation and rotation and elastic-plastic.



---

*Source: 4412788-2022-11-16.xml*