A Simple Locking-Alleviated 4-Node Mixed-Collocation Finite Element with Over-Integration,for Homogeneous or Functionally-Graded or Thick-Section Laminated Composite Beams

Leiting Dong,Ahmed S.El-Gizawy,Khalid A.Juhany and Satya N.Atluri

1 Introduction

1.1　Locking of low-order isoparametric primal elements

It is known that primal finite elements,based on low-order isoparametric displacement interpolations,suffer from shear locking for beam-shaped structures.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.Consider a four-node quadrilateral element as an example(see Fig.1),wherein the physical Cartesian coordinates as well as the Cartesian displacement components are interpolated using the same Lagrange shape functions:

where the subscript i(=1,2)denotes the index of the Cartesian coordinates,and the superscripts 1 to 4 denote the nodes of the 4-node isoparametric element.

From such a mesh-based interpolation,we have:

The Cartesian strain components are derived by using the following chain rule:

where repeated indices indicate summation.

1.2　Selective reduced-order integrations

Equations for primal FEMs are derived from the symmetric Galerkin weak-form or the equivalent principle of minimum potential energy,see[Atluri(2005);Dong,Alotaibi,Mohiuddine and Atluri(2014)]:

where ke,fe,qeare element stiffness matrix,load vector,and displacement vector respectively.N is the shape function for displacement interpolations,B is the derived matrix for strain interpolations,D is the matrix of elastic stiffness.Andˉt,ˉb represent applied tractions and body forces respectively.

For a homogeneous-material 4-noded quadrilateral element,it is clear that at least a 2 by 2 Gauss quadrature is necessary to exactly-evaluate ke.However,because fully-integrated 4-node elements are too stiff forbending problems,selectivereduced-order integrations are typically used in commercial FEM codes to improve such a decomposition is simple,by expressing the strain energy density function using Lamé constants.For anisotropic materials,the decomposition of strain energy density into dilatational and shear parts is not straight-forward,see[Hughes(1980)].Moreover,if each element represents a piece of an inhomogeneous material,such as the functionally-graded or a thick-section laminated composite beam considered in this study,over-integration is necessary instead of under-integration.In this sense,rigorously formulated hybrid and mixed finite elements are more favorable as compared to the numerical tricks of reduced-order integrations,to alleviate shear-locking in 4-noded quadrilaterals.

1.3　Hybrid&mixed finite elements

In order to overcome the shear locking,hybrid and mixed elements independently assume a stress/strain/displacement field,and derive FEM stiffness matrices using multi-field variational principles,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);Pian and Sumihara(1984);Punch and Atluri(1984);Xue,Karlovitz and Atluri(1985);Pian and Wu(1988)].Among the many variants of four-node hybridstress elements developed by Pian,the Pian-Sumihara(1984)element is currently considered to be among the best-performing four-noded elements.Through a geometrical perturbation analysis,[Pian and Sumihara(1984)]concluded that it is rational to assume the following distribution of Cartesian stress components within the element:

However,this is exactly equivalent to assuming the following distribution of the contra-variant components of the stress tensor,as independently given in[Xue,Karlovitz and Atluri(1985)]:

Rewriting the stress-field assumption asσ=Pγ,the stiffness matrix of the hybridstress element can be obtained by using Reissner’s variational principle[Reissner(1950)]:

Numerical examples have shown that the Pian-Sumihara element performs excellently for homogeneous isotropic beams if undistorted rectangular elements are used.But its accuracy is reduced when the mesh is significantly distorted.Many later works have tried to improve the Pian-Sumihara element by using different independently assumed fields,and different variational set-ups,see[Simo and Rifai(1990);Weissman and Taylor(1992);Yuan,Huang and Pian(1993)].However,very similar results were obtained by all these researchers for a four-node quadrilateral element.

Moreover,if we consider varying material properties within each element,such as functionally graded materials or thick-section laminated composite materials,the assumption of linear stress distribution is generally invalid if a coarse mesh is used.It is shown in the numerical examples of this study,that the Pian-Sumihara element fails for functionally graded materials and laminated structures,even in the most simple problem of pure tension.With this being understood,it is clear that assumed linearly-varying strains are more favorable than assumed stresses.

Another disadvantage of the current variational frameworks for hybrid/mixed elements is their questionable stability,because continuous Lagrange multiplier test functions are used to enforce the compatibility between the independently assumed stress/strain fields and those derived from mesh-based displacementinter-polations.Brezzi(1974)analyzed the existence,uniqueness,stability and convergence of saddle point problems and established the so-called LBB conditions.Inability to satisfy LBB conditions a-priori,in general would plague the solvability and stability of hybrid/mixed finite element equations.[Rubinstein,Punch and Atluri(1983);Punch and Atluri(1984);Xue,Karlovitz and Atluri(1985)]used sophisticated group theories to develop guidelines for selecting least-order stress interpolations which satisfy the LBB conditions for an undistorted element.For distorted elements,there is generally no rational approach to satisfy the LBB condition a-priori.In this study,we develop a new type of four-node quadrilateral element,which we denote as“CEQ4”.We demonstrate that,without using any multi-field variational principle or selective reduced-order integration,the CEQ4 gives much more accurate locking-alleviated and distortion-insensitive solutions than the Pian-Sumihara element for the modeling of homogeneous beams.We then combine CEQ4 with over-integration in the thicknessdirection,to modelthe deformation offunctionallygraded or laminated thick composite beams.It is shown that,without using higherorder theories[Lo,Christensen,and Wu(1977);Reddy and Robbins(1994)]or zig-zag displacement/stress assumptions[Carerra(2003)],the present CEQ4 can reasonably capture the correct distributions as well as jumps of in-plane stresses in the thickness direction,for functionally-graded or thick-section laminated beams,even if only a few elements are used.By using the equilibrium equations of elasticity,the transverse normal and shear stresses for thick-section laminates can be computed easily,from the computed in-plane stresses and their variation in the thickness direction.Detailed formulations and numerical examples are presented in the next 2 sections.

2 Detailed Formulation for the PresentLocking-Alleviated,Almost-Distortion-Insensitive,4-Node Planar CEQ4 Element:

2.1　Independently assumed strain field

It is postulated in the present paper that,instead of using assumed stresses,assumed strains are of fundamental interest for the development of FEMs.Thus,the local Cartesian components of the strain tensor are independently assumed as:

In Eq.(9),ˉx1-ˉx2is a local Cartesian coordinate system with its origin located at the center of the element.Thus,the local direct strains are assumed to be varying linearly with respect to Cartesian coordinates,in order to capture the basic bending deformation modes of the element,and the local shear strain is assumed to be a constant.

We further rewrite(9)in a matrix-vector notation for convenience:

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

It is also understood that if a two-field variational principle is used to enforce the methods were used in the context of Mesheless-Local Petrov Galerkin approaches in[Atluri,Han and Rajendran(2004);Avila,Han and Atluri(2011)].

In the work of[Dong and Atluri(2011)],the following 5 collocation points are preselected to relateγ1,...,γ5to nodal displacements(see Fig.1):

Figure 1:CEQ4:enforces 5 pre-defined constraints at 5 preselected collocation points.

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

The fundamental idea is to capture the basic kinematics of the 4-node element,to accurately model each deformation mode of tension,bending,and shear.In order to do this,we first study theinfinitesimal deformation of an infinitesimalfiberABin Fig.2(a).As illustrated in many textbooks of solid mechanics,such as[Fung and Tong(2001)],the ratio of stretch in the direction of the fiber’s axis can be calculated as:

Thus in order to model the basic kinematics modes of the four-node element,the following scheme of collocation is used in this study:

At pointA:

With these 5 equations,the five parameters ofγ1,···γ5are determined as:

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

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

The presently developed four-node quadrilateral element is denoted as CEQ4.

2.3　Some remarks on CEQ4

Remark 1:For the assumed linearly-varying strain filed,one can find an equivalent displacement field:

Without much derivation,one can verify that this method is entirely equivalent to the assumed strain formulation with the collocation scheme presented in the last subsection.

Remark2:Because ofthe assumption oflinearly-varying strain fields,itis obvious that a 2 by 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,we can use “over-integration”to accurately compute the stiffness matrix.For continuously graded materials,3 by 3 Gauss quadrature is good enough.However,for very-thick laminates,it is more convenient to either use a layer-wise 2 by 2 Gauss quadrature,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 3:If only a few elements of CEQ4 are used to model thick-section beams,the transverse normal and shear stresses directly computed by Eq.(9)may be inaccurate.In this study,following the work of[Timoshenko and Goodier(1970)],we use a stress-recovery approach to compute the distribution oftransverse stresses,by considering the equilibrium equations of linear elasticity.With bending stressesσxxcomputed by Eq.(9),the distribution of transverse stresses can be obtained by numerically evaluating:

wherey=y0denotes the lower edge of the beam.By several numerical examples in section 3,it is demonstrated that this approach gives excellent solutions for transverse stresses,even for very-thick-section laminated beams.

3 Numerical Examples

3.1　Homogeneous cantilever beams

Figure 3:A homogeneous cantilever beam subjected to a bending load or a shear force at the free-end,modeled by 2 distorted elements.

Figure 4:Computed vertical displacement at point A of the homogenous cantilever beam subjected to bending load.

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

Figure 6:Computed vertical displacement at point A of the homogenous cantilever beam subjected to shear load.

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

Figure 8:Computed transverse shear stress of the homogenous material cantilever beam subjected to shear load.

3.2　Functionally-graded materials

In this subsection,we study structures composed of functionally-graded materials.The first example is a functionally-graded square plate subjected to a tensile load along the upper side.The analytical solution for this problem is given in[Kim and Paulino(2002)].As shown in Fig.9,the plate has a unit height,width,and thickness.A plane stress condition is considered.The Young’s modulus is exponentially varying in thexdirection,i.e.E=eβx,β=log 5.Thus,we haveE=1 at the left side,andE=5 at the right side.We also considerv=0 for illustration purposes.This problem is solved by using the primal quadrilateral element,the Pian-Sumihara element,and the mixed-collocation element(CEQ4)presented in this study.The plate is modeled by only one element.Because of the exponentially varying material parameters,3 by 3 Gauss quadrature is used for evaluating the stiffness matrix.The computed vertical displacement along the upper side,and computed tensile stress along the lower side,are given in Figs.10-11.It is shown that,even for this most simple problem,the stress distribution is not linearly varying.Thus the Pian-Sumihara element yields large computational errors for the stress distribution.On the other hand,exact solution is obtained by the assumed strain mixed collocation element presented in this study.

We also consider a functionally graded cantilever beam subjected to a unit bending load or a unit shear force at the free end.As shown in Fig.12,the length and thickness of the beam are 5 and 1 respectively.Young’s modulus is exponentially varying in theydirection,i.e.E=eβy,β=log 5.Thus we haveE=1 at the lower side,andE=5 at the upper side.We also consider thatv=0 for this problem.

Figure 9:A functionally-graded square plate subjected to a uniform tensile load,modeled by 1 element.

Figure 10:Computed vertical displacement along the upper edge of the square plate.

Figure 11:Computed tensile stress along the lower edge of the square plate.

Figure 12:A functionally-graded cantilever beam subjected to a bending load(by 1 element)or a shear force(by 5 elements).

Figure 13:Computed bending stress along the left side of the functionally-graded cantilever beam(pure bending case).

Figure 15:Computed shear stress of the functionally-graded cantilever beam(endshear case).

For the pure bending case,we use merely one element to solve the problem,with 3 by 3 Gauss quadrature to evaluate the stiffness matrix.The computed bending stress along the left side is given in Fig.13,as compared to the analytical solution of[Kim and Paulino(2002)].Because of shear-locking,primal FEM gives meaningless solutions.Moreover,the Pian-Sumihara element cannot accurately capture the stress distribution because of its linear stress assumption.On the other hand,an accurate solution is obtained by using the mixed-collocation FEM(CEQ4)presented in this study.

For the shear-load case,we use 5 even-sized elements along the length direction to solve the problem.The computed bending stress along the linex=0.5(midspan of the first element)is compared to the analytical solution of[Zhong and Yu(2007)],see Fig.14.And computed transverse shear stress,which is invariant with respect tox,is plotted againstyin Fig.15.Because of the inaccuracy of primal FEM and the Pian-Sumihara element,only the results for CEQ4 are presented in Fig.14-15.Excellent agreement is found between the solution by CEQ4,and the exact solution given in[Zhong and Yu(2007)].

3.3　Thick-section laminated composite beams

In this subsection,we consider examples of thick-section laminated beams.In all of these examples,the Pian-Sumihara element and the primal 4-node element give very poor solutions unless a very fine mesh is used,thus only results by the presently developed mixed collocation FEM are presented.

Figure 16:A simply-supported thick laminated beam subjected to sinusoidal lateral load.

The first example is the classic problem of a 2-ply(0°/90°)laminated beam.Each layer of the laminate is composed of a Graphite/epoxy composite,with the following material parameters:

whereLdenotes the fiber direction andTdenotes the transverse direction.

The length and thickness of the beam are 12 inches and 1 inch respectively,so in the length direction,and only 1 element in the thickness direction,using the mixed-collocation FEM presented in this study.In Fig.17-19,we compare the computed normal stresses the shear stress with the analytical solution of Pagano(1969).Excellent agreements are obtained between the numerical and analytical solutions.

Figure 17:Computed in-plane normal stress at the mid-span of the 2-ply laminated beam.

Figure 18:Computed transverse normal stress at x=3.5(mid-span of the third element),of the 2-ply laminated beam.

Figure 19:Computed shear stress at x=1.25(right-side of the first element),of the 2-ply laminated thick beam.

Figure 20:Computed bending stress at the mid-span of the 50-ply laminated beam(l/h=6)with 8 elements by the mixed-collocation FEM(CEQ4)presented in this study.

Figure 21:Computed transverse shear stress at x=0.75(right-side of the first element),of the 50-ply laminated beam(l/h=6)with 8 elements by the mixedcollocation FEM(CEQ4)presented in this study.

Figure 22:Computed bending stress at the mid-span of the 50-ply laminated beam(l/h=12)with 8 elements by the mixed-collocation FEM(CEQ4)presented in this study.

Figure 23:Computed transverse shear stress at x=1.5(right-side of the first element)of the 50-ply laminated beam(l/h=12)with 8 elements by the mixedcollocation FEM(CEQ4)presented in this study.

Figure 24:Computed bending stress at the mid-span of the 50-ply laminated beam(l/h=6)with 375,000 Q4 elements by NASTRAN.

In the last example,we study a complex 50-ply([0°/90°]25)simply-supported laminated beam.We considerh=1 for this laminated beam,thus each ply is 0.02 inch in thickness.And we consider two different aspect ratios:(a)l/h=6 as a very-thick-section laminated beam;(b)l/h=12 as a slightly-thick-section laminate beam.The same graphite/epoxy material as in the last example is considered for each lamina,and the same sinusoidal load is applied on the upper edge of the laminated beam.We solve both cases using only 8 elements presented in this study.In Figs.20-23,the computed bending stresses and transverse shear stresses are compared to the analytical solution of Pagano(1969).It is shown that,reasonably-accurate bending and transverse shear stresses are obtained by the mixed-collocation FEM presented in this study.The only slight discrepancy is that,the present method gives a linear distribution of stresses in each of the 0°and 90°lamina,while analytical solution of Pagano(1969)demonstrates a slightlynonlinear trend of the bending stress,for the very-thick-section laminate beam(l/h=6).Extension of the present 4-node element to 8-node or 9-node planar elements should be able to give a better approximation of stress distributions,see[Bishay and Atluri(2012)].

Moreover,for the very-thick-section laminate beam(l/h=6),the computed bending and transverse shearstressesobtained by NASTRANwith 375,000 Q4 elements are also plotted in Fig.24-25.The results obtained by using NASTRAN compare well with those obtained by using the CEQ4 presented in this study,although the number of elements used in obtaining the results from NASTRAN is about five orders of magnitude more than the number of CEQ4 elements.

Figure 25:Computed transverse shear stress at x=0.75 of the 50-ply laminated beam(l/h=6)with 375,000 Q4 elements by NASTRAN.

4 Conclusion

A four-node mixed finite element is developed,following the initial work by[Dong and Atluri(2011)].The present element independently assumes a 5-parameter linearly-varying Cartesian strain field.The 5 parameters for the assumed Cartesian strains are related to the Cartesian nodal displacements,by enforcing a set of predefined constraints at 5 pre-defined collocation points.A scheme of overintegration is also proposed,for evaluating the stiffness matrices if functionallygraded or thick-section laminated composite materials are considered.Through several numerical examples,it is clearly shown that,the CEQ4 is much more accurate than the Pian-Sumihara element for the modeling of homogeneous beams.For functionally-graded materials,the presently-developed element can accurately capture not only the in-plane stress distribution and its variation in the thickness direction,but also the transverse shear stresses,even when very few elements are used.However,the Pian-Sumihara element fails because the assumption of linearly-varying stress-field is generally invalid for such materials,if a very coarse mesh is used.For thick laminated composite beams,reasonably accurate solutions are obtained even by using only 1 element of CEQ4 in the thickness direction.The methodology presented in this paper thus provides a very simple alternative to the complex higher-order theories and zig-zag assumptions of displacements/stresses in the thickness direction,for thick laminates.The approach presented in this paper can be very easily extended to study the effects of Z-Pins and stitching on interlaminar stresses in laminated composites,for structural applications.These,and the extension of CEQ4 to C0elements of higher-order,for plates and shells as well as to multi-physics will be pursued in future studies.

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.;Han,Z.D.;Rajendran,A.M.(2004):A new implementation of the meshless finite volume method,through the MLPG “mixed”approach.CMES:Computer Modeling in Engineering&Sciences,vol.6,no.6,pp.491—514.

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

Avila,R;Han,Z.D.;Atluri,S.N.(2011):A novel MLPG-finite-volume mixed method for analyzing Stokesian flows&study of a new vortex mixing flow.CMES:Computer Modeling in Engineering&Sciences,vol.71,no.4,pp.363-396.

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-µ-mechanical modeling of solids,without using multi-field variational principles.CMES:Computer Modeling in Engineering&Sciences,vol.84,no.1,pp.41-97.

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.;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.

Fung,Y.C.;Tong,P.(2001):Classical and Computational Solid Mechanics.World Scientific.

Hughes,T.J.(1980).Generalization ofselective integration proceduresto 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.

Pagano,N.J.(1969):Exact solutions for composite laminates in cylindrical bending.Journal of composite materials,vol.3,issue 3,pp.398-411.

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.(1950):On a variational theorem in elasticity.Journal of Mathematical Physics,vol.29,pp.90-95.

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

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.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 conditionsformixed-hybrid finite elementsolutions based on Reissner’s variational 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.