A Simple Locking-Alleviated 3D 8-Node Mixed-Collocation C0 Finite Element with Over-Integration,for Functionally-Graded and Laminated Thick-Section Plates and Shells,with&without Z-Pins

Leiting DongAhmed S.El-GizawyKhalid A.Juhany and Satya N.Atluri

1 Introduction

It is known that primal finite elements of deformable solids,based on low-orderC0isoparametric displacement interpolations,suffer from shear locking for typical engineering structures with large length-to-thickness ratios,such as beams,plates,and shells.This is mainly because of the incompleteness of FEM displacement interpolations,as well as the incompleteness of the strains derived from the interpolated displacement fields,see[Atluri(2005)].Selective-reduced-order integrations[Hughes(1980)]are widely used in commercial software packages to improve the accuracy of FEM solutions by reducing the bending stiffness of elements.However,it relies on the decomposition of the element strain energy density into a dilatational part and a shear part,which is not straight-forward for anisotropic materials and nonhomogeneous composite materials.Moreover,for functionallygraded materials and laminated structures,over-integration is necessary in order to accurately evaluate the element stiffness matrix,see[Dong,El-Gizawy,Juhany and Atluri(2014)].

A more rational way of alleviating shear locking is to independently assume relatively complete stress/strain/displacement fields,and derive high-performance hybrid/mixed finite elements,see[Pian(1964);Atluri(1975);Atluri,Gallagher and Zienkiewicz(1983)].One of the most popular is the hybrid-stress type of element,see[Pian(1964);Pian and Chen(1983);Rubinstein,Punch and Atluri(1983);Punch and Atluri(1984);Pian and Sumihara(1984);Xue,Karlovitz and Atluri(1985);Pian and Wu(1988);Yuan,Huang and Pian(1993)].Various versions of assumed strains mixed elements were also developed in[Simo and Rifai(1990);Weissman and Taylor(1992)].These hybrid and mixed finite elements are all based on multi-field variational principles,using continuous Lagrange multiplier test functions to enforce the compatibility between the independently assumed stress/strain fields and those derived from mesh-based displacement interpolations.Such a stencilnotonly limitsthe optimization ofelementperformances,but also makes hybrid/mixed elements suffer from LBB conditions and saddle-pointproblem instabilities,see[Brezzi(1974),Rubinstein,Punch and Atluri(1983)].

In[Dong and Atluri(2011)],a new framework of developing hybrid/mixed elements was proposed,without using multi-field variational principles or continuous Lagrange multipliers.The essential idea was to enforce the compatibility between primitive and mixed variables by simple collocation at a set of pre-defined points within the element.This approach avoids LBB conditions,and provides great flexibility in selecting collocation points&methods to improve the accuracy and robustness of developed elements.The essential idea was successfully used to develop a series of computational grains,for direct numerical modeling of complex and random microstructures of heterogeneous materials,see[Dong and Atluri(2012a,b,c,2013);Bishay and Alturi(2012,2013)].In[Dong,El-Gizawy,Juhany and Atluri(2014)],a locking-alleviated,and almost-distortion-insensitive 4-node quadrilateralC0element(CEQ4)was also developed based on this approach,by defining a set of more rational constraints at each collocation point,to accurately model each deformation mode of tension,bending,and shear.By combining CEQ4 with over-integration along the thickness direction,it was shown that functionallygraded and laminated thick-section beamscan also be accurately modeled by CEQ4 in a very simple manner,without using higher-order theories[Lo,Christensen,and Wu(1977);Reddy and Robbins(1994)]or zig-zag displacement/stress assumptions[Carerra(2003)],and even without using theories of beams/plates/shells by Euler,Bernoulli,Timoshenko(1953),or of Reissner(1945)and Mindlin(1951).In this study,we extend the previous version of CEQ4 to a three-dimensional 8-node brickC0element,which we name as CEH8.The present element independently assumes an 18-parameter linearly-varying Cartesian strain field,which is then related to the Cartesian strains derived from mesh-based Cartesian displacement interpolations,by exactly enforcing 18 pre-defined constraints at 18 preselected collocation points.We then combine CEH8 with over-integration in the thickness direction,to model the deformation of functionally-graded or laminated composite plates and shells,with an arbitrary number of laminae.It is shown that,without using higher-order theories or zig-zag theories for laminated thick composite plates and shells,the present simple 3DC0element can reasonably capture the correct distributions as well as jumps of in-plane stresses,even if only one 3DC0element is used in the thickness direction.With a stress recovery approach using equilibrium equations of 3D elasticity,the transverse normal and shear stresses can also be computed easily,from the computed in-plane stresses and their variation in the thickness direction.Moreover,this simple eight-node 3DC0element is also used for the study of z-pinned laminated plates,by simply adding the stiffness of the z-pins to the stiffness matrix of each element.It was shown that by adding stiffer z-pins in the thickness direction,inter-laminar normal stresses and sheared stresses can be reduced,which is expected to alleviate the delamination of composite structures.

2 Detailed Formulation for the PresentLocking-Alleviated,Almost-Distortion-Insensitive,8-Node Three-Dimensional C0 CEH8 Element

Figure 1:The eight-node brick C0 element with non-dimensional coordinates.

2.1Independently assumed strain field

We firstly define a local Cartesian coordinate systemˉx1?ˉx2?ˉx3,with its origin located at the center of the element.The orthonormal base vectors of the local Cartesian coordinate system are defined as follows:

The procedure for determining such a strain-field assumption is described in detail in the first section of[Dong,El-Gizawy,Juhany and Atluri(2014)]:firstly derive strains from mesh-based displacement interpolations,then find out the locking-part between derived shear&normal strains,finally eliminate the locking-part of shear strains,leading to the locking-free strain-field assumption as given in Eq.(2).

We further rewrite Eq.(2)in a matrix-vector notation for convenience:

2.2Enforcing the compatibility between the independently assumed strain and displacement fields

Then in order to relate undetermined parametersγfor assumed strains to nodal displacements q,a set of 18 collocation points are adopted,as shown in Fig.3,which are the same as those used in[Bishay and Alturi(2012)].And the following rational collocation scheme is implemented following previous work of[Dong,El-Gizawy,Juhany and Atluri(2014)]:

Figure 2:(a)Stretch of an infinitesimal fiber(b)Change of the angle between two infinitesimal fibers.

where g1,g2,g3are covariant base vector evaluated at each collocation point.With these 18 equations,the 18 parameters ofγ1,···γ18are determined:

The strain fields are thus related to the nodal displacements q by:

The stiffness matrix is determined from the strain energy stored in element ?e:

We denote the presently developed 8-node brick element as CEH8.

2.3Some remarks on CEH8

Remark 1:Because of the assumption of linearly-varying strain fields,it is obvious that a 2×2×2 Gauss quadrature is necessary if each element is used to model a piece of a homogeneous material.However,if a non-homogeneous material within the element is considered,such as functionally-graded materials or thick-section laminated composites with an arbitrary number of laminae,we can use“overintegration”to accurately compute the stiffness matrix.For continuously graded materials,a 3×3×3 Gauss quadrature is good enough.However,for laminated plates and shells,it is more convenient to either use a layer-wise two-point Gauss quadrature in the thickness direction,or use a simple Trapezoidal rule in the thickness direction,with the number of sampling points depending on the number of plies in the thickness,to evaluate the stiffness matrix of the element.

Remark 2:If only one element of CEH8 is used in the thickness direction,the transverse normal and shear stresses directly computed by Eq.(2)may be inaccurate.In this study,we use a stress-recovery approach to compute the distribution of transverse stresses,by considering the equilibrium equations of 3D linear elasticity.For example,for FG and laminated plates,the distribution of transverse stresses can be obtained by numerically evaluating:

wherez=z0denotes the lower surface of plate.

For cylindrical shells,the distribution of transverse stresses can also be evaluated,by numerical solving the following 3 differential equations:

In Eq.(10),the left hand-side invovles stress components to be recovered,and the right-hand side are directly evaluated from the solutions of CEH8.Each equation is a first-order single-variable ODE,which can be solved with a variety of computational methods,see[Dong,Alotaibi,Mohiuddine,and Atluri(2014)].In this study,simple collocation of Eq.(10)is implemented at a variety of points along the thickness direction.Combined with the traction free condition at the inner surface of the cylindrical shell,stress componentsσrθ,σrz,σrrcan be efficiently recovered from the computed in-plane normal and shear stresses.

Remark 3:Higher-order theories[Reddy and Robbins(1994)]and zig-zag theories[Carerra(2003)]for beams,plates,and shells are popularized in the current literature for analyzing functionally-graded and laminated structures.For example,third-order theory of plates by[Reddy(1984)]adopts the following expansion of displacements in the thickness direction:

This not only complicates the problem by having 7 dependent variablesu0,v0,w0,?x,?y,?z,φzinstead of 3 variables ofu,v,w,but also requiresC1continuous trial functions forw0,?z,φz,which is extremely disadvantageous for the development of general-purpose finite elements of plates and shells.

Layer-wise theories express displacements in each layer of the laminated structure in terms of polynomial interpolations.For example,[Reddy(1987)]expresses displacements in thekthlayer of the laminate as:

where?ikare Lagrange interpolation functions(linear,quadratic,etc.)in the thickness direction.This further complicates the problem by having additional dependent variables in each layer,and may lead to ill-conditioned system of equations for laminated structures with a large number of plies.

In contrast to the above-mentioned higher-order or layer-wise theories for plates and shells,the currently-developed locking-alleviated 8-nodeC0brick element has the simplest topology with only 3 variables in each node.It automatically takes account of multi-layer effects(with an arbitrary number of laminae)by using overintegration along the thickness direction,without increasing the number of DOFs for each element.Moreover,the current framework of analyzing FG and laminated plates and shells with over-integration can also be combined with other 3DC0elements such as 20-node and 27-node bricks,which are already mature in most general-purpose FEM packages such as ANSYS and ABAQUS.Using lockingfree solid elements for direct numerical modeling of FG and laminated structures saves the trouble of developing specific theories for plates and shells,and thus provides an “one-size fit all”procedure for universal modeling of both bulk solids and engineering structures.

3 Numerical Examples

3.1Homogeneous beams

3.2Functionally-graded plates

Figure 4:A homogeneous cantilever beam(E=1,v=0)subjected to a bending load or a shear force at the free-end,modeled by 2 distorted elements.

Figure 5:Computed vertical displacement at point A of the homogenous material cantilever beam subjected to bending load at the free end.

Figure 6:Computed bending stress at point B of the homogenous material cantilever beam subjected to bending load at the free end.

Figure 7:Computed vertical displacement at point A of the homogenous material cantilever beam subjected to shear load at the free end.

Figure 8:Computed bending stress at point B of the homogenous material cantilever beam subjected to shear load at the free end.

Figure 9:A cantilever functionally graded plate(E=5z,v=0)subjected to tensile,bending,or shear load at the free end.

Figure 10:Computed tensile stress at x=5,y=0 for the square plate subjected to a uniformly distributed tensile load(N=1/length)at the free end.

Figure 11:Computed bending stress at x=5,y=0 for the square plate subjected to a uniformly distributed bending load(M=1/length)at the free end.

Figure 12:Computed bending stress at x=5,y=0.5 for the square plate subjected to a uniformly distributed shear load(P=1/length)at the free end.

Figure 13:Computed out-of-plane shear stress at x=5,y=1 for the square plate subjected to a uniformly distributed shear load(P=1/length)at the free end.

In this subsection,a cantilever unit-thickness 10×10 square plate is considered.Young’s modulus is exponentially varying in thezdirection,i.e.E=eβz,β=log5.Thus,we haveE=1 at the lower surface andE=5 at the upper surface.We also considerv=0 for illustration purposes.Three load cases are considered,where the free end of the plate is subjected to uniformly distributed tensile,bending,and shear loads respectively.Analytical solutions for this problem were given in[Kim and Paulino(2002);Zhong and Yu(2007)].For the case with tensile or bending load,the plate is modeled by only one element.And for cases with shear load,the plate is modeled by a 10×10 mesh withC0brick elements.Because of the exponentially varying material parameters,a 3×3×3 Gauss quadrature is used for evaluating the stiffness matrix.In Figs.10-13,computed in-plane normal stresses and out-of-plane shear stresses are compared to analytical solutions.Almost exact solution is obtained with the currently developed CEH8C0brick elements even though very coarse mesh is used.

3.3Laminated plates

Firstly,we consider a thick-section 50-ply([0?/90?]25)laminated 10 inches×10 inches square plate.Each layer of the laminate is composed of a Graphite/Epoxy composite,with the following material parameters:

whereLdenotes the fiber’s direction andTdenotes the transverse direction.

The thickness of the plate is 1 inch,so that each ply is 0.02 inch in thickness.The laminated plate is simply-supported at each edge.And it is subjected to a uniform lateral loadq=1 psi.

We solve this problem using a uniform 10×10 mesh with CEH8 elements,as well as using NASTRAN.The comparison between the meshes by CEH8 and by NASTRAN is given in Fig 14.Because of the large aspect ratio for each layer of the laminated plate,a very fine mesh is necessary for NASTRAN,withabout 1.5 million DOFs and about 2.5 hours of computational time on a regular PCwith i7 CPU.On the contrary,CEH8 only requires a very coarse mesh with726 DOFs and about 5 seconds of computational time.Computed in-plane and out-of-plane stresses by NASTRAN and CEH8 are also shown in Figs.15-16.It can be clearly seen that similar computational results are obtained even thoughCEH8 requires about 2000 times less computational time as compared to NASTRAN.

We also consider a different plate with a very-high aspect ratio.The same composite material,the same 50-ply laminate,the same thickness,and the same boundary conditions and loads are considered.The only difference for the current laminated plate is thata=b=1000 inches,so that it has an aspect ratio of 1000.We also solve this problem with 10×10 CEH8 elements,with computed stresses shown in Fig.17.This demonstrates that the current simple eight-nodeC0brick element can be used to tackle problems of both thick-section and thin-section plates,without having to resorting to theories of plates and shells.

Figure 14:Finite element model for the 50-ply laminated plate(a/h=10)by(a)NASTRAN and(b)currently-developed 3D C0 CEH8.

3.4Laminated shells

Figure 15:Computed σxx,σyy,σzz at x=y=4.5 inches,and computed σxz,σyz at x=y=1 inch,for the thick-section laminated plate(a/h=10),with CEH8.

Figure 16:x=y=1 inch,for the thick-section laminated plate(a/h=10),with NASTRAN.

Figure 17:x=y=100 inches,for the thin-section laminated plate(a/h=1000),with CEH8.

Figure 18:Finite element model for the 50-ply laminated shell by 400 CEH8 elements.

3.5Z-pin Reinforcements for Laminated Plates

In this section,we study the effect of z-pins on reducing inter-laminar stresses of laminated plates.The same simply-supported thick 50-ply laminated plate as studied in section 3.3 is considered here.A uniformly-distributed tensile load(p=1)is applied to the upper surface of the plate.100 uniformly distributed steel z-pins are used to reinforce the laminated plate in the thickness direction.The material parameters for steel z-pins areE=29×106psi,v=0.3.The diameter of each z-pin is 0.2523 inch,so that the volume fraction of z-pins is 5%.We solve this problem with 10×10 currently-developed 3DC0CEH8 elements,by simply adding the contributing stiffness of z-pins into the stiffness matrix of each CEH8 element:

Figure 20:Computed out-of-plane shear stress at tion in z direction constrained to simulate a plane strain case.

Figure 21:Computed bending stress at formation in z direction unconstrained to simulate a 3D case.

Figure 22:Computed out-of-plane shear stress at deformation in z direction unconstrained to simulate a 3D case.

Figure 23:A simply-supported 50-ply laminated plate(a/h=10)with 5%volume fraction of z-pins,modeled by 100 currently-developed 3D C0 CEH8

Figure 24:Computed out-of-plane shear stresses at x=y=1 inch,for the simplysupported thick-section laminated plate,with and without z-pin reinforcements

Figure 25:Computed out-of-plane normal stress at x=y=5 inches,for the simply supported thick-section laminated plate,with and without z-pin reinforcements

where Dzrepresents the elastic stiffness of the z-pin(?z)within each element.The computed out-of-plane normal and shear stresses are shown in Fig.24-25.It can be seen that z-pin reinforcements can help reducing inter-laminar stresses of laminated plates,thus reducing the possibility of structural failure caused by delamination.

4 Conclusion

A locking-alleviated 3D eight-nodeC0brick element is developed,following the previous work of[Dong,El-Gizawy,Juhany and Atluri(2014)].The present element independently assumes an 18-parameter linearly-varying Cartesian strain field.The 18 parameters for the assumed Cartesian strains are related to the Cartesian nodal displacements,by enforcing a set of prede fined constraints at 18 prede fined collocation points.The constraints are rationally de fined to capture the basic kinematics of the 3D 8-node element,and to accurately model each basic deformation mode of tension,bending,shear,and torsion.A scheme of over-integration is also used,for evaluating the stiffness matrices for functionally-graded materials or thick-section laminates with an arbitrary number of laminae.Through several numerical examples,it is clearly shown that,the current approach can obtain very accurate solutions for in-plane stresses of FG and laminated structures,even by using only one CEH8 3DC0element in the thickness direction.The out-of-plane normal and shear stress are also accurately recovered using equations of 3D elasticity.By adding the contributing stiffness of z-pins into the stiffness matrix of CEH8,it is also demonstrated that the presently developed method can be used to study the effect of using z-pin reinforcements to reduce the inter-laminar stresses of composite structures,in a very simple and computationally-efficient manner.In contrast to higher-order or layer-wise theories of plates and shells that are popularized in the current literature,the currently-developed locking-alleviated 8-nodeC0brick element saves the trouble of developing specific theories plates and shells,but simply uses the widely-available theories of elasticity for the modeling of FG and laminated structures.

Acknowledgement:This research is supported by the Mechanics Section,Vehicle Technology Division,of the US Army Research Labs,under a collaborative research agreement with UCI.The first author acknowledges the financial support of Natural Science Foundation Project of Jiangsu Province(Grant no.BK20140838).This work was funded by the Deanship of Scientific Research(DSR),King Abdulaziz University,under grant number(15-135-35-HiCi).The authors,therefore,acknowledge the technical and financial support of KAU.

Atluri,S.N.(1975):On hybrid finite element models in solid mechanics.Advances in Computer Methods for Partial Differential Equations,R.Vichnevetsky(eds.),AICA,pp.346-356.

Atluri,S.N.;Gallagher,R.H.;Zienkiewicz,O.C.(1983):Hybrid and mixedfinite element methods.John Wiley&Sons.

Atluri,S.N.(2005):Methods of Computer Modeling in Engineering and the Sciences,Tech Science Press.

Brezzi,F.(1974):On the existence,uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers.Revue Fran?aise D’automatique,Informatique,Recherche Opérationnelle,Analyse Numérique,vol.8,no.2,pp.129-151.

Bishay,P.L.;Atluri,S.N.(2012):High-Performance 3D Hybrid/Mixed,and Simple 3D Voronoi Cell Finite Elements,for Macro-&Micro-mechanical Model-ing of Solids,Without Using Multi-field Variational Principles.CMES:Computer Modeling in Engineering&Sciences,vol.84,issue 1,pp.41-97.

Bishay,P.L.;Atluri,S.N.(2013):2D and 3D Multiphysics Voronoi Cells,Based

on Radial Basis Functions,for Direct Mesoscale Numerical Simulation(DMNS)of the Switching Phenomena in Ferroelectric Polycrystalline Materials.CMC:Computers,Materials&Continua,vol.33,no.1,pp.19-62.

Carrera,E.(2003):Historical review of zig-zag theories for multilayered plates and shells.Applied Mechanics Reviews,vol.56,issue 3,pp.287-308.

Dong,L.;Atluri,S.N.(2011):A simple procedure to develop efficient&stable hybrid/mixed elements,and Voronoi cell finite elements for macro-µmechanics.CMC:Computers Materials&Continua,vol.24,no.1,pp.41-104.

Dong,L.;Atluri,S.N.(2012a):T-Trefftz Voronoi Cell Finite Elements with elastic/rigid inclusions or voids for micromechanical analysis of composite and porous materials.CMES:Computer Modeling in Engineering&Sciences,vol.83,no.2,pp.183-220.

Dong,L.;Atluri,S.N.(2012b):Development of 3D Trefftz Voronoi Cells with Ellipsoidal Voids&/or Elastic/Rigid Inclusions for Micromechanical Modeling of Heterogeneous Materials.CMC:Computers,Materials&Continua,vol.30,no.1,pp.39-82.

Dong,L.;Atluri,S.N.(2012c):Development of 3D Trefftz Voronoi Cells with Ellipsoidal Voids&/or Elastic/Rigid Inclusions for Micromechanical Modeling of Heterogeneous Materials.CMC:Computers,Materials&Continua,vol.30,no.1,pp.39-82.

Dong,L.;Gamal,S.H.;Atluri,S.N.(2013):Stochastic Macro Material Properties,Through Direct Stochastic Modeling of Heterogeneous Microstructures with Randomness of Constituent Properties and Topologies,by Using Trefftz Computational Grains(TCG).CMC:Computers,Materials&Continua,vol.37,no.1,pp.1-21.

Dong,L.;Alotaibi,A.;Mohiuddine,S.A.;Atluri,S.N.(2014):Computational methods in engineering:a variety of primal&mixed methods,with global&local interpolations,for well-posed or ill-Posed BCs.CMES:Computer Modeling in Engineering&Sciences,vol.99,no.1,pp.1-85.

Dong,L.;El-Gizawy,A.S.;Juhany,K.A.;Atluri,S.N.(2014):A Simple Locking-Alleviated 4-Node Mixed-Collocation Finite Elementwith Over-Integration,for Homogeneous or Functionally-Graded or Thick-Section Laminated Composite Beams.CMC:Computers,Materials&Continua,vol.40.no.1,pp.49-77.

Fung,Y.C.;Tong,P.(2001):Classicaland computationalsolid mechanics.World Scientific.

Hughes,T.J.(1980).Generalization ofselective integration procedures to anisotropic and nonlinear media.International Journal for Numerical Methods in Engineering,vol.15,no.9,pp.1413-1418.

Lo,K.H.;Christensen,R.M.;Wu,E.M.(1977):A high-order theory of plate deformation—part 2:laminated plates.Journal of Applied Mechanics,vol.44,issue 4,pp.669-676.

Mindlin,R.D.(1951):Influence of rotatory inertia and shear on flexural motions of isotropic,elastic plates.Journal of Applied Mechanics,vol.18,pp.31–38.

Pian,T.H.H.(1964):Derivation of element stiffness matrices by assumed stress distribution.A.I.A.A.Journal,vol.2,pp.1333-1336.

Pian,T.H.H.;Chen,D.(1983):On the suppression of zero energy deformation modes.International Journal for Numerical Methods in Engineering,vol.19,issue 12,pp.1741–1752.

Pian,T.H.H.;Sumihara,K.(1984):Rational approach for assumed stress finite elements.International Journal for Numerical Methods in Engineering,vol.20,issue 9,pp.1685–1695.

Pian,T.H.H.;Wu,C.C.(1988).A rational approach for choosing stress terms for hybrid finite element formulations.International Journal for Numerical Methods in Engineering,vol.26,issue 10,pp.2331-2343.

Punch,E.F.;Atluri,S.N.(1984):Development and testing of stable,invariant,isoparametric curvilinear 2-and 3-D hybrid-stress elements.Computer Methods in Applied Mechanics and Engineering,vol.47,issue 3,pp.331-356.

Reissner,E.(1945):The effect of transverse shear deformation on the bending of elastic plates.Journal of Applied Mechanics,vol.12,pp.69-77.

Reddy,J.N.(1984):A simple higher-order theory for laminated composite plates.Journal of Applied Mechanics,vol.51,issue 4,pp.745-752.

Reddy,J.N.(1987):A generalization of two-dimensional theories of laminated composite plates.Communications in Applied Numerical Methods,vol.3,issue 3,pp.173-180.

Reddy,J.N.;Robbins,D.H.(1994).Theories and computational models for composite laminates.Applied Mechanics Reviews,vol.47,issue 6,pp.147-169.

Ren,J.G.(1987):Exact solutions for laminated cylindrical shells in cylindrical bending.Composites Science and Technology,vol.29,issue 3,pp.169-187.

Rubinstein,R.;Punch,E.F.;Atluri,S.N.(1984):An analysis of,and remedies for,kinematic modes in hybrid-stress finite elements:selection of stable,invariant stress fields.Computer Methods in Applied Mechanics and Engineering,vol.38,issue 1,pp.63-92.

Simo,J.C.;Rifai,M.S.(1990):A class of mixed assumed strain methods and the method of incompatible modes.International Journal for Numerical Methods in Engineering,vol.29,issue 8,pp.1595-1638.

Timoshenko,S.(1953):History of strength of materials,McGraw-Hill New YorkTimoshenko,S.P.;Goodier,J.N.(1970):Theory of Elasticity,3rdedition,Mc-Graw Hill.

Weissman,S.L.;Taylor,R.L.(1992):A unified approach to mixed finite element methods:Application to in-plane problems.Computer Methods in Applied Mechanics and Engineering,vol.98,issue 1,pp.127-151.

Xue,W.;Karlovitz,L.A.;Atluri,S.N.(1985):On the existence and stability conditions formixed-hybrid finite elementsolutions based on Reissner’svariational principle.International Journal of Solids and Structures,vol.21,issue 1,1985,pp.97-116.

Yuan,K.Y.;Huang,Y.S.;Pian,T.H.(1993):New strategy for assumed stresses for 4-node hybrid stress membrane element.International Journal for Numerical Methods in Engineering,vol.36,issue 10,pp.1747-1763.

Zhong,Z.;Yu,T.(2007):Analytical solution of a cantilever functionally graded beam.Composites Science and Technology,vol.67,issue 3,pp.481-488.