JIA Hong－xue(賈紅學(xué))，LIU Xi－la(劉西拉)

School of Naval Architecture，Ocean and Civil Engineering，Shanghai Jiao Tong University，Shanghai 200240，China

Force－Based Quadrilateral Plate Bending Element for Plate Using Large Increment Method

JIA Hong－xue(賈紅學(xué))*，LIU Xi－la(劉西拉)

School of Naval Architecture，Ocean and Civil Engineering，Shanghai Jiao Tong University，Shanghai 200240，China

A force-based quadrilateral plate element(4NQP13)for the analysis of the plate bending problems using large increment method(LIM)was proposed．The LIM，a force-based finite element method(FEM)，has been successfully developed for the analysis of truss，beam，frame，and 2D continua problems．In these analyses，LIM can provide more precise stress results and less computational time consumption compared with displacement-based FEM．The plate element was based on the Mindlin-Reissner plate theory which took into accountthe transverse shear effects．Numerical exampleswerepresented to study itsperformance including accuracy and convergence behavior，and the results were compared with the results have been obtained from the displacementbased quadrilateral plate elements and the analytical solutions．The 4NQP13 element can analyze the moderately thick plates and the thin plates using LIM and is free from spurious zero energy modes and free from shear locking for thin plate analysis．

large increment method(LIM);displacement-based finite element method(FEM);Mindlin-Reissner plate theory;spurious zero energy modes;shear locking

Introduction

Plate elements are widely used in structural analysis in engineering．A number of research works have been proposed over several decades about the development of simple and efficient plate finite elements．The thin plate element based on the Kirchhoff theory is required that the deflections and their derivatives(C1continuity) should be continuous on the interfaces between elements．Because the C1continuity is rather difficult to guarantee for the displacement－based thin plate element，the focus has switched to the Mindlin－Reissner plate theory．The plate elements based on Mindlin－Reissner plate theory become the dominant model to calculate the bending of thin and moderately thick plate because only C0continuity is required．All existing Mindlin－Reissner plate elements can be generally divided into two groups，such as the displacementbased finite element method(FEM)［1－3］and the mixed/hybrid element method［4－7］．But the displacement－based Mindlin－Reissner plate elements are commonly numerical over－stiffness and suffer from shear locking when applied in thin plate problems，thereby giving inaccuracy for the analysis of thin plate problems．Many relatively successful methods have been proposed to circumvent shear locking，such as the reduced integration［8］， selective integration［9］， penalty－parameter modification［10］，assumed shear strain method［11］，and element boundary interpolation using Timoshenko beam function［12］．The selective and reduced integration can avoid the shear locking effectively in displacement－based model plate elements，but may result in singularity of the system stiffness matrix and such an element often exhibits spurious zero energy modes．As an alternative to the FEM，the large increment method(LIM) is extended to analyze the thin and moderately thick plate problems．

Zhang and Liu［13］first presented the LIM for material nonlinear problem，based on the force method and the generalized inverse matrix theory．In LIM，the governing equations are setin three parts: equilibrium equations，compatibilityequations，and the local physical equations．Unlike the traditional force method，the independent elemental (generalized)forces variables are adopted as system unknowns，and the basic determinate structure is no longer needed，so the LIM becomes a systematic method and it is easy to implement in the program and computer．The advantage of LIM is that it separates the global equilibrium and compatibility equations from the local constitutive equations(the system matrix does not change even for nonlinear material problems)，and the parallel computing in the time and spatial domain．For the traditional force method，the main unknowns are the redundant forces．However，there is no unique way to choose the redundant force in the statically indeterminate and a bad selection may cause computational instability，so it is difficult to implement in the program and computer．

In recent years，LIM has been extended for solving elastic perfectly plasticanalysisofplane framestructuresunder monotonic and cyclic loading［14－15］，2D continuum elastoplastic problems［16］，and the elastic plate problems［17］．In these papers，LIM has showed considerable computational savings and accuracy．In this paper，a new force－based plate element is proposed for the analysis of the thin plate and moderately thick plate problems．The new plate element(4NQP13)shows the advantages of avoiding the shear locking for thin plate problems and the spurious zeroenergymodes．The efficiencyand accuracy of the LIM are illustrated with a few benchmark problems and the results are compared with the analytical solution and the existing displacement－based quadrilateral plate elements．Only linear elastic material and small deformation are considered．

1 Governing Equations of Plate Element

In LIM，on a solid body Ω，its boundary is represented as S．The boundary is divided into two parts:the surface force intensity boundary Sσand the displacement boundary Su，where Sσ∩Su=0．The stress－resultant M and displacement u at each point in any element can be given by the element generalized force variables and element nodal displacement variables．The stress－resultant fields of the plate element are written as

where M is the generalized force vector in element，Z is the shape function of generalized force parameter vector in element，and Fedenotes the elemental generalized inner force vector．

The displacement fields in element are presented as

and the strain vector in element is given by

where u is the displacement vector in element，N is the shape function of node displacement vector in element，dedenotes the nodal displacement vector of the element，κ is the strain vector of the element，L is the differential operator，and B is the strain－displacement matrix in element．

If the displacement set in the element is compatible，the weak form of the equilibrium equations of the element can be obtained by using the equation of virtual work

Substituting Eq．(1)and Eq．(3)into Eq．(4)，Eq．(4) can be rewritten as

and the element equilibrium equations can be expressed as

where

In Eq．(7)，Ceis the element equilibrium matrix，Peis the element equivalent nodal force vector．

Thereby，theequilibrium equations ofthebody are obtained by assembling the element equilibrium equations and are presented as

where C is an m×n equilibrium matrix of the structure，and it is assembled from the element equilibrium matrices．For the body，a statically indeterminate structure，the system coefficient of equilibrium C is an m×n non－square matrix with m＜n．F is an n×1 generalized inner force vector of the body，and P is an m×1 nodal load vector．

Similarly，if the force set in the body is in equilibrium，using Eq．(1)，Eq．(3)，and the principle of complementary virtual work，it can be obtained that

Thus，the assembled system compatibility equations of the body can be expressed as

where D is the nodal displacement vector of the body and δ is the generalized deformation vector of the body．

The elemental constitutive equations can be written as

where

In Eq．(12)，δedenotes the element generalized deformation vector，Φeis the element flexibility matrix，and φ represents the relationship between stress and strain of the material．

Here the generalized inverse theory will be utilized to solve the linear equilibrium equation denoted by Eq．(8)．The system equilibrium matrix C is an m×n non－square matrix with m＜n for the statically indeterminate structure and the C－1does not exist．The generalized inverse of matrix C is defined as

then

where Im×mis an identify matrix．

The general solution of Eq．(8)can be written as

and the orthogonal projection operator β can be presented as

where the orthogonal projection operator α=

By eliminating the nodaldisplacementD using the generalized inverse theory，Eq．(10)can be rewritten as

then，

and Eqs．(17)and(18)become the generalized compatibility equations．

Based on the discussion above，the stress－resultant fields are used to form the equilibrium and constitutive equations and play an important role in obtaining the precise stress results，and the displacement fields are used to form the compatibility equations in the LIM．The stress－resultant and displacement fieldswillbe presented in subsections 1．1 and 1．2，respectively．

1．1 Assumed stress-resultant fields

In Eq．(7)and Eq．(12)，matrix Z appears in both the element equilibrium matrix Ceand the element flexibility matrix Φe．Therefore，it is important to properly assume stressresultant interpolation polynomials for obtaining the accurate results．The 4－node quadrilateral element has 12 degrees of freedom(DOF)，three of which are related by three equations of equilibrium．Therefore，a stress field with a minimum of 9 independent forces is needed to describe the stress field．The stress field should be selected in a rational methodology for avoiding the spurious zero energy mode．The number of independent generalized forces m should satisfy

where n is the total number of nodal displacement in an element and r is the number of kinematic DOF．Equation(19)is the only necessary condition but not sufficient for the absence of the zero energy modes(ZEM)．The generalized inverse matrix theory is applied to detecting the zero energy modes and the necessary and sufficient condition are given as follows

where*Ceis the element equilibrium coefficient matrix of after the element kinematic DOF is eliminated from Ce，the details will be presented in Section 3．

The assumed stress－resultant fields should be satisfied with the equilibrium equations

Based on the discussion above，the reasonable stressresultant fields are given in Table 1．The 4－node plate element is named 4NQP13．Here4N means4 nodes，QP means quadrilateral plate element，and 13 means the number of the unknown forces．

Table 1 Assumed stress－resultant fields for the plate element

1．2 Assumed displacement fields

The main assumption of the plate bending element based on Mindlin－Reissner plate theory states that a straight line originally normal to the mid－surface of the plate remains straight but not necessarily normal to the deformed mid－surface．The node numbering for the element in isoparametric coordinates is shown in Fig．1．Therefore，the displacement fields for the transverse displacement and two rotations are interpolated in the interior of the element in terms of nodal values and assumed as

where wi，θxi(θyi)，and Niare the corresponding displacements，rotation values，and the shape functions of node i，respectively．The shape functions are given below

Fig．1 Node numbering of the plate element

2 Basic Theory of LIM

In LIM，the governing equations of the body are global equilibrium equations， compatibility equations， and local constitutive equations．In the case of small deformation，the global equilibrium equations and the compatibility equations are constantlinear algebraic equations． However， the local constitutive equations are nonlinear if the constitutive relations are material nonlinear．Unlike the FEM，LIM separates the linear global equilibrium and compatibility equations from the local constitutive equations．For the indeterminate structures，its equilibrium matrixis non－square．The generalized inverse matrix method is employed to obtain the set of solutions of the non－square matrix equations directly．First the initial forces are obtained and the constitutive equations are then utilized to obtain the deformations．After that the compatibility equations are used to check whether the deformations satisfy the deformation compatible condition．Now the main goal of LIM is to find an optimized solution F that makes‖βδ‖ to approach zero．For solving the optimization problem，an iterative procedure is employed with the conjugate gradient method to find the correct solution．The governing equations of the plate element are given in Section 1．For the linear elastic material plate bending problem，the procedure of LIM is as follows．

Step 1 Using the theory of generalized inverse of matrix to solve the body equilibrium Eq．(8)and find out a special solution F0as the initial value of the iteration

Step 2 Iterate from Fnto Fn+1

1)Calculate the deformation vector δncorresponding to Fn:

2)Check whether the deformations are sufficiently compatible or not:

If ε(δn) ≤ e0(e0is the given error tolerance)，then the deformation are sufficiently compatible，no more iterations are required，and go to Step 3．Otherwise，the deformation vector δnis not compatible and Fnshould be modified to improve the compatibility of the deformation vector．

3)Calculate the search direction Sn:

where K(δn)is the current stiffness matrix，which can be obtained by calculating the inverse of the flexibility matrix Φ(Fn)．

4)Calculate the search step length:

5)Fn+1can be calculated as:

Step 3 Obtain the final results

3 The Spurious Zero Energy Modes Test

The stress－resultant and displacement fields within an element are assumed，respectively．The stress－resultant fields cannot be chosen arbitrary since the spurious zero energy modes (ZEM)should be avoided in the plate element．In order to check the ZEM，the generalized inverse matrix theory is applied to detecting the ZEM．

Using Eq．(1)，the complementary virtual work equation can be written as

when

from Eqs．(12)and(31)，the left of Eq．(31)can be written as

where Ceis the element equilibrium matrix，δFeis the virtual generalized element inner force，deis the nodal displacement vector of an element．By eliminating the kinematic DOF from Ceand de，Eq．(33)can be represented as

When

it can be found that de*has non－trivial solution by using the theory of generalized inverse of matrix and it can be presented as

where X is arbitrary vector and has the same dimension as the dimension of Fe．

It can be shown that the existence of the non－trivial solution can lead to rank deficiency of the element equilibrium matrix*Ce．The rank deficiency of the element equilibrium matrix is a sign of the appearance of spurious zero energy modes in the plate element．

Equation(35)is represented as

where r is the number of rigid body modes in an element．It is recognized that the number of independent generalized inner force parameters in the assumed stress－resultant fields has to satisfy

and the necessary and sufficient conditions are given as

where*Ceis the element equilibrium coefficient matrix after the element kinematic DOF is eliminated from Ce．

The rank of an elemental flexibility matrix Φehas to satisfy

The correct rank of the element equilibrium coefficient matrix guarantees that the plate element can avoid from spurious zero energy modes．

4 Numerical Results and Discussion

Several standard numerical examples are presented in this section．The square plate is used for assessing the convergence and the accuracy．The shear locking of the proposed plate elements is tested．The circular plate is used to investigate the sensitivity to distorted mesh．The results of the proposed elements are compared with existing 4－node quadrilateral displacement－based plate elements，namely the Q4 element［8］，the Q4－R element［18］，the RDKQM element［19］，the MQP4 element［20］， the MITC4 element［21］， and the analytic solutions［22－24］．

4．1 Convergence tests

A square moderately thick plate with simply supported/ clamped boundary conditions subjected to a uniform load．The parameters are the plate side L=0．4 m，the thickness t=0．08 m，the elasticity modulus E=2．1×1011Pa，Poisson's ratio ν= 0．3，and the uniform loading q=1010Pa．Due to symmetries both in geometry and the loading，a quarter of the plate is analyzed．The 4×4 regular meshes are shown in Fig．2．The corresponding convergence trends are shown in Figs．3－6．

Fig．2 4×4 mesh for quarter of the square plate

Fig．3 Central deflection for a clamped square plate with uniform load

Fig．4 Central deflection for a simply supported square plate with uniform load

Fig．5 Central moment coefficient for a clamped square plate with uniform load

Fig．6 Central moment coefficient for a simply supported square plate with uniform load

In Figs．3－6，it is clear that the mesh convergence of the proposed quadrilateral force－based 4NQP13 element is as good as the famous displacement－based Q4－R and RDKQM element，and is better than the displacement－based Q4 element and the force－based MQP4 element in estimating the central deflection and central moment coefficient．

4．2 Shear locking tests

A simply supported square plate with various thickness/span ratios t/L subjected to uniform load is used to test shear locking phenomenon of the proposed plate element using LIM，and the mesh size is 6×6 in a quarter of plate．The analytical solutions for the centraldeflection coefficientand centralmoment coefficient are obtained by Timoshenko and Krieger［22］and Liu et al．［23］for thin plate and moderately thick plate problems．

The results of central deflection coefficient and central moment coefficient presented in Tables 2 and 3 show that the new force－based plate element can avoid the shear locking phenomenon．However，the Q4 element is not free from shear locking and cannot get the correct solution for the thin plate．

Table 2 Central deflection coefficient for a simply supported plate with uniform load

Table 3 Central moment coefficient for a simply supported plate with uniform load

4．3 Mesh distortion test

A circular plate with simply supported boundary conditions is subjected to a uniform loading．The dimensionless quantities parameters are:the radius R=5，the thickness t=1 and 2．5 for the moderately thick plate，t=0．1 for the thin plate，the elasticity modulus E=10．92，the Poisson's ratio ν=0．3，the uniform loading q=1．Due to symmetry of the plate，a quarter of the plate analyzed with 12 and 48 elements is shown in Fig．7．The analytical solutions for the central displacement and central moment are obtained by Ayad and Rigolot［24］: simply supported plate:

where

Fig．7 Two kinds of mesh for quarter of a circular plate

In Tables 4 and 5，it should be noted that the 4NQP13 element can lead to accurate results for the moderately thick plate and thin plate，and it is insensitive to the mesh distortion．The 4NQP13 element has good monotonic convergence towards the exact solutions and the accuracy of the 4NQP13 element for the central deflection is quite well like the Q4－R element．It also shows that the 4NQP13 element gives more accurate stressresultant than the Q4 element．

Table 4 Central deflection for a simply supported circular plate with uniform load

Table 5 Central moment for a simply supported circular plate with uniform load

5 Conclusions

The LIM has been extended for the analysis of thin and moderately thick plates．A 4－node quadrilateral plate bending element(4NQP13) based on Mindlin－Reissnertheory is proposed in this paper． Numerical comparisons with displacement－based plateelements show thatthe 4NQP13 element for the analyses of the thin and moderately thick plates using LIM provides excellent convergence and good prediction for both deflection and stress resultants，and it is insensitive to the mesh distortion．The main advantages are that the 4NQP13 element is free from spurious zero energy modes and do not exhibit shear locking in the thin plate problems．A rational selection of the stress－resultant fields which can avoid shear locking and spurious zero energy modes is very important．

The stress distributions within the plate are particularly important in nonlinear analysis．As a force－based FEM，LIM can provide more precise stress results with high efficiency and low cost．Considering both the spatial and time domain parallel computation，the LIM has great potential for solving the material nonlinearity problems with significant computational savings．As a result，LIM will be a powerful alternative method to solve the elastoplastic problems of plate problems．

［1］Brezzi F，Evans J A，Hughes T J R，et al．New Rectangular Plate Elements Based on Twist－Kirchhoff Theory［J］．Computer Methods in Applied Mechanics and Engineering，2011，200(33): 2547－2561．

［2］Castellazzi G，Krysl P．A Nine－Node Displacement－Based Finite Element for Reissner－Mindlin Plates Based on an Improved Formulation of the NIPE Approach［J］．Finite Elements in Analysis and Design，2012，58(1):31－43．

［3］Serpik I N．Development of a New Finite Element for Plate and Shell Analysis by Application of Generalized Approach to Patch Test［J］．Finite Elements in Analysis and Design，2010，46 (11):1017－1289．

［4］Hu B，Wang Z，Xu Y C．Combined Hybrid Method Applied in the Reissner－Mindlin Plate Model［J］．Finite Elements in Analysis and Design，2010，46(5):428－437．

［5］Choo Y S，Choi N，Lee B C．A New Hybrid－Trefftz Triangular and Quadrilateral Plate Elements［J］．Applied Mathematical Modelling，2010，34(1):14－23．

［6］Darilmaz K，Kumbasar N．An 8－Node Assumed Stress Hybrid Element for Analysis of Shells［J］．Computers＆Structures，2006，84(29/30):1990－2000．

［7］Carstensen C，Xie X P，Yu G Z，et al．A Priori and a Posteriori Analysis for a Locking－Free Low Order Quadrilateral Hybrid Finite Element for Reissner－Mindlin Plates［J］．Computer Methods in Applied Mechanics and Engineering，2011，200(9/10/11/12): 1161－1175．

［8］Zienkiwicz O C，Taylor R C，Too J M．Reduced Integration Technique in General Analysis ofPlates and Shells［J］．International Journal Numerical Methods in Engineering，1971，3 (2):275－290．

［9］Hughes T J R，Cohen M，Haron M．Reduced and Selective Integration Techniques in the Finite Element Analysis of Plates［J］．Nuclear Engineering and Design，1978，46(1):203－222．

［10］Tessler A．A Priori Identification of Shear Locking and Stiffening in Triangular Mindlin Elements［J］．Computer Methods in Applied Mechanics，1985，53(2):183－200．

［11］Hinton E，Huang H C．A Family of Quadrilateral Mindlin Plate Elements with Substitute Shear Strain Fields［J］．Computers＆Structures，1986，23(3):409－431．

［12］Zhang H X，Kuang J S．Eight－Node Reissner－Mindlin Plate Element Based on Boundary Interpolation Using Timoshenko Beam Function［J］．International Journal Numerical Methods in Engineering，2007，69(7):1345－1373．

［13］Zhang C J，Liu X L．A Large Increment Method for Material Nonlinearity Problems［J］．Advances in Structural Engineering，1997，1(2):99－109．

［14］Barham W，Aref A，Dargush G．Development of the Large Increment Method for Elastic Perfectly Plastic Analysis of Plane Frame Structures Under Monotonic Loading［J］．International Journal of Solids and Structures，2005，42(26):6586－6609．

［15］Barham W，Aref A，Dargush G．On the Elastoplastic Cyclic Analysis of Plane Beam Structures Using a Flexibility－Based Finite Element Approach［J］．International Journal of Solids and Structures，2008，45(22/23):5688－5704．

［16］Long D B，Liu X L．Development of 2D Hybrid Equilibrium Elements in Large Increment Method［J］．Journal of Shanghai Jiaotong University，2013，18(2):205－215．

［17］Jia H X，Long D B，Liu X L．The Force－Based Quadrilateral Plate Elements for Plate Analysis Using Large Increment Method［C］．Proceeding the 2nd International Conference on Civil Engineering and Building Materials，Hong Kong，China，2012: 221－225．

［18］Malkus D S，Hughes T J R．Mixed Finite Element Methods－Reduced and Selective Integration Techniques:a Unification of Concepts［J］．Computer Methods in Applied Mechanics and Engineering，1978，15(1):63－81．

［19］Chen W J，Cheung Y K．Refined Quadrilateral Element Based on Mindlin/Reissner Plate Theory［J］．International Journal for Numerical Methods in Engineering，2000，47(1/2/3):605－627．

［20］Dhananjaya H R，Nagabhushanam J，Pandey P C．Bilinear Plate Bending Element for Thin Moderately Thick Plate Using Integrated Force Method［J］．Structural Engineering and Mechanics，2007，26(1):43－68．

［21］Bathe K J，Dvorkin E H．A Four－Node Plate Bending Element Based on Mindlin－Reissner Plate Theory and Mixed Interpolation［J］． International Journal for Numerical Methods in Engineering，1985，21(2):367－383．

［22］Timoshenko S P，Krieger S W．Theory of Plates and Shells［M］．New York:McGraw－Hill，1959:181－226．

［23］Liu J，Riggs H R，Tessler A．A Four Node Shear－Deformable Shell Element Developed via Explicit Kirchhoff Constraints［J］．International Journal Numerical Methods in Engineering，2000，49(8):1065－1086．

［24］Ayad R，Rigolot A．An Improved Four－Node Hybrid－Mixed Element Based upon Mindlin's Plate Theory［J］．International Journal Numerical Methods in Engineering，2000，55(6):705－731．

O34;TU311．4

A

1672－5220(2015)03－0345－06

date:2014－01－08

National Natural Science Foundation of China(No．10872128)

* Correspondence should to be addressed to JIA Hong－xue，E－mail:jiahongxue211@126．com