Multi-model finite element scheme for static and free vibration analyses of composite laminated beams 1

A transition element is developed for the local global analysis of laminated composite beams. It bridges one part of the domain modelled with a higher order theory and other with a 2D mixed layerwise theory (LWT) used at critical zone of the domain. The use of developed transition element makes the analysis for interlaminar stresses possible with significant accuracy. The mixed 2D model incorporates the transverse normal and shear stresses as nodal degrees of freedom (DOF) which inherently ensures continuity of these stresses. Non critical zones are modelled with higher order equivalent single layer (ESL) theory leading to the global mesh with multiple models applied simultaneously. Use of higher order ESL in non critical zones reduces the total number of elements required to map the domain. A substantial reduction in DOF as compared to a complete 2D mixed model is obvious. This computationally economical multiple modelling scheme using the transition element is applied to static and free vibration analyses of laminated composite beams. Results obtained are in good agreement with benchmarks available in literature.


INTRODUCTION
Laminated composites are finding varied engineering applications as these materials possess a good environmental resistance, strength and are light in weight.Failure modes of laminated composites include the delamination failure leading to separation of the layers and loss of integrity of the structure.Due to the heterogeneous material properties of the different layers of the composite laminate high interlaminar stresses develop at the interfaces.For the kinetic equilibrium in the thickness direction, these interlaminar stresses need to be continuous between the layers.Pipes and Latin American Journal of Solids and Structures 12 (2015) 2061-2077 Pagano (1970); Rybicki (1971) brought out that delamination failure can be attributed to the interlaminar transverse shear and normal stresses.For a sound design, an accurate evaluation of the interlaminar stresses becomes inevitable.Different analytical and finite element equivalent single layer (ESL) and layerwise (LW) models have been proposed for analysis of composites.The ESL models in the literature predict the global parameters with a reasonable accuracy but fail to predict the continuity of interlaminar stresses.The displacement based LW models also fail to predict continuity of transverse stresses.These are more accurate relative to ESL models for global parameters but need higher computational effort.For estimation of the transverse stresses, a separate stress function or integration of stress equilibrium equations is essential.LW mixed model with transverse stresses as nodal DOF are accurate for global parameters and also for local interlaminar stresses.
Accurate elastic analysis of interlaminar stresses in composite laminates was presented by Srinivas and Rao (1970); Pagano (1969Pagano ( , 1970)); Pagano and Hatfield (1972).As an offshoot of the plate formulation, many researchers have presented theories for analysis of composite laminated beams.Kant (1982) presented a comprehensive higher order theory for analysis of thick plates.Kant et al. (1997) applied the fully cubic theory for the dynamic analysis of beams.Various forms of higher order sub-theories were evolved by eliminating some terms from the comprehensive cubic theory (Marur and Kant, 1996;Manjunatha and Kant, 1993a;1993b).Reddy (1987) presented higher order ESL as well as displacement based LW formulations.Contributions were also made by Lo et al. (1978); Spilker (1982), and others in the form of higher order theories.The formulation of Spilker had a stress based hybrid element having separate through thickness distributions for stress and displacements.To overcome the limitation of displacement based ESL, Rao et al. (2001) proposed an analytical solution based on mixed formulation in which the transverse stresses are invoked as DOF ensuring their continuity.This model has been employed for static and dynamic analyses of laminated plates and beams.Shimpi and Ghugal (1999) presented a trigonometric shear deformable LW theory capable of maintaining the continuity of transverse stresses.Desai and Ramtekkar (2002) presented a mixed finite element LW model with transverse stresses also included in the set of nodal DOF.The issue of continuity is inherently resolved and the results are seen to be in good agreement with the exact solutions.The formulation was applied to static analysis of composite beams.Ramtekkar et al. (2002) extended the mixed formulation for free vibration analysis of beams.Adyodgu (2006) presented a pair of shear deformable beam theories by employing parabolic and exponential shape functions respectively in conjunction with the classical beam theory for free vibration analysis of angle ply beams.Marur and Kant (2007) presented a higher order theory for free vibration analysis of angle ply beams.
Literature indicates that LW approach is the most suitable finite element formulation for accurate estimation of the interlaminar stresses.Within the LW models reported in literature, it is observed that the mixed formulation (Desai and Ramtekkar, 2002) yields very accurate solution and also satisfies the continuity requirements.However, such formulation is computationally expensive.Thus, layerwise formulation may not be suitable at all locations in the domain of a laminate.Amongst the ESL formulations the higher order theory (Kant et al., 1997) with 8 DOF(HOSTB8) has been reported to be comprehensive and computationally economical for the estimation of global parameters.Structures 12 (2015) 2061-2077 In this work, the positive traits of the higher order ESL theory and the mixed LWT are used together to obtain global as well as the local transverse stress parameters with a good accuracy.Laminate is modelled using a stack of 2D elements having mixed DOF in the zone where the transverse stresses are critical.Remaining zone is modelled using higher order ESL (HOSTB8).A unique transition element is developed to bridge the interface of ESL and mixed LWT based elements.

THEORETICAL FORMULATION FOR BEAM ANALYSIS
Three models have been formulated for analysis of transversely loaded laminated composite beams consisting of several orthotropic layers each having different properties; (c) Model 3: This model is based on a local global finite element procedure to take advantage of computational efficiency of higher order ESL theory and accuracy of the LW mixed model.

Model 1: Development of ESL based theory (HOSTB8)
Displacements in two principal directions of the laminated beam as a fully cubic function of the thickness co-ordinate are The above displacement field eliminates any requirement of shear correction factor and chances of shear locking.Here ( 0 u and 0 w ) are deformations in X , Z (laminate co-ordinate) directions at the mid-plane.( x θ and z θ ) are rotations at mid-plane about principal directions of laminated beam and ( * 0 u , * x θ , * 0 w and * z θ ) are higher order terms from Taylor's series.By using material property, strain displacement relationship and the principle of minimum potential energy, the stiffness matrix for the laminate is developed.By using shape functions similar to those used for stiffness evaluation, the mass matrix is also developed.Detailed formulation can be seen in Kant et al. (1997).A fournode isoparametric line element of Lagrange family is used to discretize a laminated beam.
Integration is performed by employing 5x5 Gauss quadrature rule for the extension, bending, mass component and 3x3 Gauss rule for the shear part.Structures 12 (2015) 2061-2077 2.2 Model 2: Development of mixed LW model A 6-node two-dimensional plane stress element based on mixed formulation is used by considering displacement fields ( ) , u x z , and ( ) , w x z having quadratic variation along the length of beam and cubic variation in the transverse direction.The cubic variation has been adopted to invoke the transverse stresses as the nodal parameters in addition to the nodal deformations.The displacement field is expressed as

Latin American Journal of Solids and
where Further, mik a ( m = 0, 1, 2, 3; i = 1, 2, 3) are the generalized coordinates.Variation of displacement fields has been assumed to be cubic through the thickness of element, although there are only two nodes along ' z ' axis of an element.Derivative of displacement with respect to the thickness coordinate has also been included in the displacement field.Such a inclusion is required for invoking transverse stress components z σ , and xz τ as nodal DOF in the formulation.Further, it also ensures parabolic variation of the transverse stresses through the thickness of an element.
By making use of the elasticity relationship, displacement field ( ) (2) can be shown to be Here, i = 1, 2, 3 for the nodes with 1 ξ = − , 0 ξ = and 1 ξ = , respectively; q = 1, 2 and p = 3, 4 for the nodes with 1 η = − and 1 η = , respectively for node numbers 1 to 6 and; f and 4 f correspond to derivative of displacements with respect to thickness co-ordinate, whereas, Application of principle of minimum potential energy is used to develop element property matrix.Detailed formulation can be seen in Desai and Ramtekkar (2002).Numerical integration of system matrices (Stiffness property matrix, Mass property matrix and load vector) has been performed by using a 3x3 integration scheme in the length and 5x5 integration scheme in the thickness direction.Solids and Structures 12 (2015) 2061-2077

Model 3: Development of transition between ESL (HOSTB8) and 2D mixed LW model
Compatibility between two differently modelled sub-domains (by using Model 1 and Model 2) is enforced by degenerating a 2D mixed element through kinematic constraints compatible with deformations predicted by ESL HOSTB8 element.Cook et al. (2003) has presented a methodology to connect dissimilar elements.
A 2D-to-1D transition element has one or two edges of 2D element that are kinematically restrained to enforce compatibility with adjacent ESL elements.Such a edge is denoted as a transition edge in the sequel.The 2D element on the transition edge is conditioned for compatibility with the DOF of the ESL (HOSTB8) element to ensure continuity of the combined model.Such an element acts as a transition element to connect two independently modelled sub-domains.Transition is achieved by placing a stack of such transition elements used in different layers of a laminate at the transition edge.
A pair of incompatible mesh formulations is shown in Fig. 1 wherein a four-node ESL element with eight DOF per node (node numbers denoted with a prime) is connected to a stack of 2D mixed elements with four DOF per node (two translations and two transverse stresses). where Similar transformation can be performed for mass matrix of an element.The transformation in Eq. ( 9) degenerates the transition edge of the 2D element which becomes follower to the corresponding HOSTB8 leader node.All elements in the interior of the local transition edge are 6node elements with all the nodes modelled with the mixed formulation.The stress DOF at the 2D nodes on the transition edge are condensed prior to imposition of the restraint for complete compatibility.Considering the stiffness matrices of the ESL elements, transition elements and the interior LW mixed elements the global equation can be obtained in the following form after assembly.
Latin American Journal of Solids and Structures 12 (2015) 2061-2077 Here are the global stiffness property and mass matrices, respectively; M are the element property and mass property matrices of th i element, respectively, formed by using mixed LWT; The displacement vector of the transition element ˆtr q is composed of the DOF of the ESL node on the transition edge, and the DOF of the 2D nodes on the interior of element.The transition element that is developed by the application of the restraints consists of 5 nodes and 24DOF in case of an one edge transition.

Analysis of plate under static loading and vibration analysis
Standard finite element procedures are used to assemble and formulate global system matrices.For the static analysis the load deformation relationship is Gauss elimination method is used for determination of DOF in Eq. (9).By using Hamilton's variational principle as brought out in Bathe (1997) and solution to the equation of motion of laminate, the characteristic equation of the eigenvalue problem is formed  σ is estimated through Mesh 2 by using a stack of 2D elements with mixed DOF at and in the vicinity of ( )

ILLUSTRATIVE EXAMPLES AND DISCUSSIONS
and the remaining part of the laminate meshed with ESL elements.
For the static and free vibration analyses, symmetrical loading on the beam and simply supported end conditions are assumed.An example with different boundary conditions (Example 6) for an four layered angle ply beam is also considered for free vibration analysis.However, the same meshing patterns are used for discretization.Implementation of this novel hybrid finite element mesh is done on full length of beam.Free vibration analysis is also performed for two meshing arrangements to further validate the soundness and efficacy of the combined model.Stress estimation in static analysis for transverse stresses has been done by using selective meshing pattern according to suitability.Results obtained for static and free vibration analyses are compared with those from higher order ESL theories and LW mixed and displacement based finite element formulations available in literature.
A combined mesh size of (10)/(2x30) implies that full length of beam is discretized with 10 elements and out of which a portion of 2 elements is modelled with 30 2D elements with mixed DOF in the thickness direction.The remaining domain is modelled with ESL (HOSTB8) elements.Thus, there will be 68 elements (8#HOSTB8+60#2D) in the entire domain.If the same laminate is modelled using only 2D elements, the total number of elements for a mesh (10x30) will be 300.Proposed combined mesh reduces total number of DOF due to reduction in number of elements required to model the domain and also ensures continuity of transverse stresses.
Beams with different layups, span to thickness ratios (S), and material properties with various end conditions have been considered for static and free vibration analyses.Boundary conditions are listed in Table 1.Material properties and normalisation factors used for static and free vibration analyses are mentioned with the results.

Example 1
A simply supported cross ply (0 °/90 °) beam with layers of equal thickness and subjected to unidirectional sinusoidal loading is considered.The beam is analysed using both combined mesh patterns for stresses ( xz τ and z σ ), longitudinal and transverse displacements ( u and w ).The transverse shear and transverse normal stresses has been selectively estimated by using Mesh 1 and Mesh 2 respectively.The results for S = 4 and 10 are compared in Table 2. Exact elasticity solution by Pagano (1969), and FE results obtained by Desai and Ramtekkar (2002), Liou and Sun (1987), Lu and Liu (1992), Manjunatha and Kant (1993b), and Shimpi and Ghugal (1999) have been presented for proper comparison.Results from the present analysis have been found to be in agreement with the established results.The variation of normalized longitudinal normal stress ( x σ ), transverse normal and shear stresses ( z σ and xz τ ) through the thickness of the beam with S = 4, have been shown in Fig. 3. Proximity of present results with the benchmarks validate the suitability of the combined model.Essential requirement of continuity of transverse stresses is also fulfilled with a substantial reduction in the computational effort.
Simply supported beams (SS) (with ESL HOSTB8 Elements at ends in Mesh 2)

Example 2
Static analysis of a three-layered cross-ply (0°/90°/0°) simply supported beam with equal thickness of each layer, has been considered in this example.The beam is subjected to uni-directional sinusoidal loading applied at top edge.The transverse shear and transverse normal stresses has been selectively estimated by using Mesh 1 and Mesh 2 respectively.Normalized maximum transverse displacement ( w ), in-plane normal stress ( x σ ) and transverse shear stress ( xz τ ) for the simply supported laminated beam (with S = 4 and 10) have been presented in Table 3. Results obtained through present analysis are compared with elasticity solution by Pagano (1969) and various analytical/ FE solutions given by Desai and Ramtekkar (2002); Lo et al. (1978); Spilker (1982); Engblom and Ochoa (1985); Toledano and Murakami (1987); Manjunatha and Kant (1993b).Present results shows good agreement with the elasticity solutions, which underlines the usefulness of the present model.Variation of normalized longitudinal normal stress ( x σ ) and transverse stresses ( z σ and xz τ ) through thickness of beam with S = 4 is presented in Fig. 3.It is observed that the results from the present multi-model scheme are in very good agreement with the elasticity solution (Pagano, 1969).It is also noted that the displacement based models (Lo et al., 1978;Engblom and Ochoa, 1985), have not been able to predict the in-plane and the transverse stresses.The continuity of the transverse stresses is ensured and the integration of equilibrium equations is advantageously avoided.the present model are compared in Table 4 with FOBT, HOBT by Marur and Kant (1996), mixed theory by Rao et al. (2001) and FEM solution by Ramtekkar et al. (2002).The through thickness variation of non dimensional inplane displacements, transverse normal and shear stresses for first three modes is shown Fig. 5.It can be observed that the results from the present studies are in good agreement with HOBT and the mixed theory with a substantial reduction in the total DOF.

Example 5
A symmetric cross-ply (0°/90°/90°/0°) thin beam under simple support condition is considered for free vibration analysis.The non-dimensional natural frequencies ( ω ) obtained from the present investigation are compared in Table 6, with the first order beam theory FOBT (4a) by Marur and Kant (1996), HOBT by Kant et al. (1998) and the mixed theory by Rao et al. (2001) and FEM solution by Ramtekkar et al. (2002).Results have been observed to be in very good agreement with the HOBT and the mixed theory.

Example 6
A four layered angle ply beam with layup ( θ °/θ °/θ °/ θ °) is considered in this example.To study the influence of ply angle, fundamental natural frequencies of laminated beam are determined  for θ = 0°, 15°, 30°, 45°, 60°, 75°, 90°.Simply supported at both ends (SS), clamped at both ends (CC) and clamped free (CF) boundary conditions are considered.CF condition with Mesh 1 meshing pattern has a stack of 2D mixed elements at the clamped end and remaining domain is modelled with HOSTB8 elements.An end condition, free clamped (FC) is also considered with Mesh 1 meshing pattern having the stack of 2D mixed elements at the free end.The non-dimensional fundamental frequencies are compared with those presented by Adyodgu (2006) in Table 7. Variation of the fundamental frequency with the ply angle for boundary condition (SS) is shown in Fig 6 .Except for θ = 15° and 30°, natural frequencies for all other orientation of plies are in close proximity of those presented by Adyodgu (2006).The consistency between results of the two meshing patterns under all boundary conditions strengthens the efficacy of the combined mesh model.The results of above examples illustrate that the present multi-modelling scheme is apt for complete analysis of laminated beams.The accuracy requirements are fulfilled without sacrificing the computational efficiency.Continuity of transverse stresses is also ensured.It is observed that Ramtekkar et al. (2002) presented few frequencies which were not reported earlier in the literature.Present combined mesh model misses out on few of those additional frequencies reported in Ramtekkar et al. (2002) but on the other hand some new frequencies are encountered.The genuineness of these new frequencies is under investigation by the authors.

CONCLUSION
A finite element formulation is presented which combines the higher order ESL sub-domain to 2D mixed LW sub-domain using a transition element.It is observed that a reasonably good accuracy is achieved in the estimation of the deformations, stresses.andnatural frequencies.to reduction in number of elements rquired for the entire domain.This makes the procedure computationally economical.Increase in the size of the local LW mixed sub-domain would make the mesh to tend towards full LW modeling and the methodology may loose its advantage.An accuracy above 90% in comparison with elasticity solution is seen to be achieved by using LW mixed model on 25% to 40% of the entire domain.The present combined formulation overcome the shortcoming of ESL in prediction of continuity of the transverse stresses.The method suggested is confirming to the elasticity principles as both the models used on the sub-domains are derived from the elasticity equations and two dimensional state of stress and strain.The results given by the developed model show reasonably good accuracy as compared with the benchmarks from the literature.

(a) Model 1 :
This model adopts a cubic displacement field in the thickness direction for displacements (U , W ) and has 8DOF per node.The theory has been identified as HOSTB8.The model is based on the plane stress state of stresses and strains.(b) Model 2: In this model, mixed finite element LWT which has two displacements (U , W ) and the transverse stresses ( , xz zτ σ ) as the nodal DOF is used.The theory is based on elasticity relationships.Therefore, requirement of any additional parameters/stress variation functions are advantageously avoided.
the nodal transverse stress variables.

Figure 1 :Figure 2 :
Figure 1: The configuration of connection between 2D mixed and HOSTB8 elements.
By developing the restraint sub-matrices kj R     for all pairs of ESL and 2D nodes, the transformation matrix R     for the entire element can be formulated by appropriately populating sub-matrices kj R     corresponding to every pair.Finite element stiffness property, mass property matrices and internal force vector for the transition element are obtained by matrix transformations using the constructed stiffness property matrix, force vector ,mass matrix of 2D element and associated transformation matrix as follows, property and mass matrices of th j transition element, respectively; stiffness and mass matrices of th l ESL element, respectively; m , n and k in Eq. (10) represent number of LWT, transition and ESL elements.

D
are global nodal load and DOF vectors respectively.{ } D is the modal vector and ω is natural frequency of laminate.Solution of Eq. (10) after imposition of boundary conditions yield natural frequency ω and corresponding modal eigenvector { } D .Sub-space iteration technique is used in the current work for solution of eigenvalue problem.

Following
two mesh patterns are possible for a symmetrically loaded simply supported laminated beams; (a) Peak xz τ is estimated through Mesh 1 by using a stack of 2D elements with mixed DOF at and in the vicinity of ( ) 0 x = and the remaining part of the laminate meshed with ESL elements.Latin American Journal of Solids and Structures 12 (2015) 2061-2077 (b) Peak z −' indicates no boundary condition imposed on that degree-of-freedom at that location '*' 0 0 q = for free vibration analysis

Figure 5 :
Figure 5: Through thickness variation of (a) Normal displacement, (b) transverse normal stress, (c) transverse shear stress, for first three modes of vibration.
Parallely a substantial reduction in the global DOF can be observed as compared to full 2D LWT model due Latin American Journal of Solids and Structures 12 (2015) 2061-2077

Table 1 :
Boundary conditions for composite beams.

Table 2 :
Comparison of the maximum transverse displacement, the in-plane normal and the transverse shear stresses for simply supported laminated beam under sinusoial loading (lamination scheme: 0 °/90 °)

Table 4 :
Comparison of non-dimensional natural frequencies

Table 5 :
has been done.The results show good agreement with the earlier reported results in the literature.Comparison of non-dimensional natural frequencies

Table 6 :
Comparison of non-dimensional natural frequencies

Table 7 :
Variation of non-dimensional fundamental frequency