Formulation of enriched macro elements using trigonometric shear deformation theory for free vibration analysis of symmetric laminated composite plate assemblies

Composite Structures 119 (2015) 38–49

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Formulation of enriched macro elements using trigonometric shear deformation theory for free vibration analysis of symmetric laminated composite plate assemblies Rita F. Rango a,⇑, Liz G. Nallim a, Sergio Oller b,c a b c

Facultad de Ingeniería, INIQUI (CONICET), Universidad Nacional de Salta, Av. Bolivia 5150, 4400 Salta, Argentina CIMNE International Center for Numerical Method Engineering, Spain UPC, Technical University of Catalonia (Barcelona Tech), Edif. C1, Campus Nord, Jordi Girona 1–3, 08034 Barcelona, Spain

a r t i c l e

i n f o

Article history: Available online 23 August 2014 Keywords: Composite plates Macro element Free vibration Trigonometric Shear Deformation Theory

a b s t r a c t The formulation of an enriched macro element suitable to analyze the free vibration response of composite plate assemblies is presented in this article. Based on the Trigonometric Shear Deformation Theory (TSDT) and the use of Gram–Schmidt orthogonal polynomials as enrichment functions a finite macro element is developed. In the TSDT framework, shear stresses are vanished at the top and bottom surfaces of the plates and shear correction factors are no longer required. The Principle of Virtual Work is applied to derive the governing equations of motion. A special connectivity matrix is obtained; so that hierarchically enriched global stiffness matrix and mass matrix of general laminated plate structures are derived, allowing to study generally coplanar plate assemblies by combining two or more macro elements. This procedure gives a matrix-eigenvalue problem that can be solved with optimum efficiency. Results of free vibration analysis for symmetric laminated plates of different thickness ratios, geometrical planform shapes and boundary conditions are presented. The accuracy of the formulation is ensured by comparing some numerical examples with those available in the literature. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Currently, composite materials have become indispensable in several applications, such as high-performance structures in many fields of civil, marine and aerospace engineering, among others. While composite materials have many advantages over conventional materials, they also present complex and challenging problems for structural analysis and design engineers and researchers. Many structural elements, such as cylinders, beams, plates and shells, can potentially be used for the analysis of composite laminates. Particularly, laminated plates of fiber reinforced composite materials, with different shapes, have great advantages and are widely used in high-performance structural components. During the last few decades the use of composite plates has increased in various engineering applications. The laminated plates are attractive structural components in many industries, because of the high stiffness-weight ratio along with the possibility to tailor

⇑ Corresponding author. E-mail address: [email protected] (R.F. Rango). URLs: http://www.unsa.edu.ar, http://www.conicet.gob.ar (R.F. Rango). http://dx.doi.org/10.1016/j.compstruct.2014.08.012 0263-8223/Ó 2014 Elsevier Ltd. All rights reserved.

the lamination scheme, which could be adapted to the requirements of design. The global deformation of laminated composite plates is, in the general, characterized by complex couplings between extension, bending, torsion and shear modes. Furthermore, due to their low transverse shear stiffness, laminated composite plates exhibit a much more significant transverse shear deformation than homogeneous isotropic plates with the same geometric dimensions, even for low thickness-to-length ratios. In order to consider these aspects into the analysis and design of laminated plates, and to exploit the potential advantages of these materials, it is necessary to develop methodologies that include these effects. Furthermore, it is necessary to have accurate analysis tools that allow arrive to appropriate and versatile designs, according to increasingly stringent requirements. For the study of plates with any thickness-to-length ratio, the formulations based on Equivalent Single Layer (ESL) theories must include the higher-order effects. This is accomplished by applying theories which take into account the shear deformation in the cinematic expressions, with the advantage of not requiring the use of shear correction factors and reproducing more accurately the distribution of interlaminar stresses in thick plates. Interesting Higher

39

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

Order Shear Deformation Theories (HOSDT) were proposed by authors like Reddy and Liu [1], Touratier [2], Soldatos [3], Kant and Swaminathan [4], Karama et al. [5], among others, and more recently, by Mantari et al. [6], Mantari and Guedes Soares [7] and Grover et al. [8]. These theories satisfy the boundary conditions corresponding to the free surfaces of the plates and represent approximately the parabolic distribution of shear stress in the thickness. Most of the HOSD existing theories have polynomial expressions for the shear deformation. For example, in the general formulations presented by Carrera [9], Carrera et al. [10] and Demasi [11–15] the unknown variables are represented throughout the thickness by polynomial functions. However, in relation to the ESL theories, in accordance with the reviews made by the authors and also with Mantari et al. [16], it is quite important to explore the behavior of other functions in the implementation of new shear deformation theories. It can be said that there is evidence of the demand generated by the higher order trigonometric theories [17], fundamentally due to the fact that they are much richer than polynomial functions, which are at the same time, simpler and more precise, and the boundaries conditions in the free surfaces of the plate are guaranteed at priori. The above mentioned trigonometric theories were applied by Shimpi and Ainapure [18] and later by Arya et al. [19] for studying laminated beams. Subsequently, Ferreira et al. [20] used for the first time a trigonometric shear deformation theory for modeling symmetric laminated square plates by a meshless method, based on global multiquadric radial basis functions, obtaining very good results. Then, Roque et al. [21] used this trigonometric theory for laminated plates but incorporating the concept of multilayer laminates, obtaining very good results for the static analysis of symmetric laminated square plates. Starting from these studies, Xiang and Wang [22] presented the analysis of free vibrations of square laminated plates, using the trigonometric shear deformation theory and inverse multiquadric radial basis functions, arriving to very good results of natural vibration frequencies for different material and geometric parameters. Recently, Mantari et al. [17,23] presented a trigonometric shear deformation theory to model laminated composite and sandwich plates, with square or rectangular planform, by formulating a discrete finite element. Then, Mantari and Guedes Soares [7,24,25] completed these studies with the analysis of graded plates and advanced composite plates, respectively. Regarding the analysis of thick plates of general geometries, Ramesh et al. [26] presented a higher-order triangular plate element based on the Third-order Shear Deformation Theory and a layer-wise plate theory of Reddy for the bending analysis of laminated composite plates. Zamani et al. [27] presented a transformation of coordinates combined with the differential equations from First Order Shear Deformation Theory, to model the problem of free vibration of laminated plates with trapezoidal and skew planform with different geometrical parameters, various aspect ratios and boundary conditions. Using a Higher Order Shear Deformation Theory, Fazzolari et al. [28] present an exact dynamic stiffness method for free vibration analysis of composite plate assemblies. Previously, Houmat and Rashid [29] presented a method for coupling isoparametric cubic quadrilateral h-elements and straight sided serendipity quadrilateral p-elements for the free vibration analysis of plates with curvilinear planforms. In previous papers, the authors have presented the formulation of hierarchical finite macro elements (h-p version of FEM), enriched with Gram–Schmidt orthogonal polynomials, using the Classical Laminated Plates Theory (CLPT) [30] and the First Order Shear Deformation Theory (FSDT) [31–33]. In this paper, the concept of macro element formulated by the authors is extended, so as to incorporate the kinematic corresponding to the Trigonometric Shear

Deformation Theory (TSDT), which allows the study of thick laminated plates and plate assemblies, due the incorporation of mapping and assembly techniques. 2. Displacement field A general quadrilateral thick laminated plate element, as shown in Fig. 1, is considered. A laminate of uninform thickness h with Nl layers is adopted for the analysis. Each layer consists of unidirectional fiber reinforced composite material. The fiber angle of kth layer counted from the surface z = h/2 is b and it is measured from the x axis to the fiber direction. Symmetric lamination of plies is considered in this work (see Fig. 1A and C). Based on the Trigonometric Shear Deformation Theory (TSDT) and taking into account the corresponding hypothesis [20], the displacement field can be described as:

@w0 ðx; y; tÞ pz þ sin /x ðx; y; tÞ @x h @w0 ðx; y; tÞ p z v ðx; y; z; tÞ ¼ z þ sin /y ðx; y; tÞ @y h  y; z; tÞ ¼ w0 ðx; y; tÞ wðx;

ðx; y; z; tÞ ¼ z u

ð1Þ

; v ; w  ¼ w0 are the displacements of a generic point on the where u mid plane (z = 0) along (x, y, z) and (/x , /y ) are the rotations of the transverse normal about y and x axis respectively. During free vibration, the displacements are assumed split in the spatial and temporal parts, being the last one periodic in time: w0(x, y, t) = w (x, y) sin(xt) where x is the radian natural frequency. The linear strains associated with the displacement fields (Eq. (1)) are given by:

8 @/ x > > @x pz < @/

8 2 9 > @ w > > < @x2 2 > = y eyy ¼ sin þ z  @@yw2 @y : ; > > > h > > > > : @/x þ @/y > : ; cxy @2 w ; 2 @x@y @y @x     cyz p pz /y ¼ cos cxz h /x h 8 9 < exx =

9 > > =

ð2Þ

3. Equations for dynamic analysis The governing equations of the problem come from the dynamic version of the Principle of Virtual Work.

Z

t2

ðdLÞdt ¼ 0

ð3Þ

t1

where L is the lagrangian and is defined as L = T  (U + V) where U, V, T are the strain energy, work done by applied forces and kinetic energy, respectively. In this article the term V is omitted as the analysis is limited to free vibration response. The virtual strain energy dU is given by:

dU ¼

Z (Z R

h=2 h

h=2

i

)

rxx dexx þ ryy deyy þ sxy dcxy þ sxz dcxz þ syz dcyz dz dxdy ð4Þ

where R is the mid-surface area of the plate (see Fig. 1). Replacing Eq. (2) into Eq. (4), the following expression is obtained:

! ! @ 2 dw0 @d/x pz @ 2 dw0 r sin þ r z  þ xx yy h @x2 @x @y2 h=2 R !   @d/y pz @ 2 dw0 pz @d/x @d/y þ sxy sin sin þ sxy z 2 þ þryy h h @y @x@y @y @x p pz p pzi o þsxz d/x cos þ syz d/y cos dz dxdy ð5Þ h h h h

dU ¼

Z (Z

h=2

"

rxx z 

40

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

Fig. 1. (A) General planform plate macro element (B) Mapped macro element (C) N-layered symmetric laminate.

The virtual kinetic energy dT is given by:

dT ¼

Z (Z R

)

h=2

_ u_ þ v_ dv_ þ wd _ w _ dz dxdy q½ud

ð6Þ

h=2

where q is the material density, which is considered uniform through the volume of the laminate. Replacing Eq. (1) into Eq. (6) the following expression is obtained: Z (Z

  _0 _0 @w pz @dw pz þ sin /_ x þ sin d/_ x z h h @x @x h=2 R      _0 _0 @w pz _ @dw pz _ _ 0 dw _ 0 dz dxdy ð7Þ þ z þ sin /y þ sin d/y þ w z h h @y @y

dT ¼

h=2

q



z

Replacing Eqs. (5) and (7) in the expression (3) of the Principle of Virtual Works the following expression can be obtained: ! ! Z t2 (Z Z h=2 " @ 2 dw @d/x pz @ 2 dw 0¼ rxx z  2 þ rxx sin þ ryy z  @x h @y2 @x h=2 t1 R !   @d/y pz @ 2 dw pz @d/x @d/y þryy þ sxy sin sin þ sxy z 2 þ h @x@y h @y @y @x i p pz p pz dzdxdy þ syz d/y cos þsxz d/x cos h h h h   Z Z h=2  _ _ @w pz @dw pz  q z þ sin /_ x z þ sin d/_ x h h @x @x h=2 R      _ _ @w pz @dw pz _ w _ dzdxdy dt þ z þ sin /_ y þ sin d/_ y þ wd z h h @y @y

ð8Þ Considering the constitutive relationship [34] and the maximum kinetic energy corresponding to a complete cycle of vibration, Eq. (8) can be written as:

Z (Z

"

2  11 @/x sin pz  Q  11 z @ w þ Q  12 @/y sin pz Q h @x2 h @x @y h=2 R !   2 2  12 z @ w þ Q  16 @/x þ @/y sin pz  Q  16 2z @ w Q @y2 h @x@y @y @x ! 2 2 @d/x pz @ dw  12 @/x sin pz  Q  12 z @ w þ Q sin  z  2 h @x h @x2 @x @x   2  22 @/y sin pz  Q  22 z @ w þ Q  26 @/x þ @/y sin pz þQ h @y2 h @y @y @x ! !  2 2 @d/ @ w p z @ dw @/ y  26 2z  16 x sin pz þ Q sin  z Q @x@y h @y2 h @y @x   2 2  16 z @ w þ Q  26 @/y sin pz  Q  26 z @ w þ Q  66 @/x þ @/y sin pz Q 2 2 @x h @y h @y @y @x !  !  2 2 @d/x @d/y pz @ dw  66 2z @ w sin  2z þ Q @x@y h @x@y @y @x  p p z p p z p pz  45 / cos þ Q  55 / cos þ Q d/x cos y x h h h h h h  i o  44 / p cos pz þ Q  45 / p cos pz d/ p cos pz dz dxdy þ Q y x y h h h h h h    Z Z h=2 @w pz @dw pz  x2 q z þ sin /x z þ sin d/x @x h @x h R h=2     @w pz @dw pz z þ sin /y þ sin d/y þ wdw dzdxdy þ z @y h @y h

0¼

h=2

ð9Þ  ij are the mechanical reduced rigidities referred to x,y axes where Q [34]. The integration of Eq. (9) throughout the thickness (z coordinate) is developed in Appendix A. After integration, Eq. (9) can be written in a matrix form as:

41

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

Fig. 2. Assembly procedure (A) general planform geometry (B) adopted mesh (C) mapped macro element (D) degrees of freedom and generalized displacement.

0¼

Z h R



h

@d/y @y

@d/x @x

@/y @y

@/x @x

2

2 S A1 i6 11 @d/y 6 @d/x S þ 4 A112 @y @x @/x @y

þ

@/y @x

3

D11 6  4 D12

D12

D16

D22

7 D26 5

D16

D26

D66



Z h R



h

@/x @x

2



@ 2 dw @x2

@/y @y

@/x @y

iT

3 þ

7 A1S26 7 5

A1S16 A1S26 A1S66 Z h 2  2 i @ 2 dw @ dw dw 2 @@x@y dxdy þ @x2 @y2

h



@2 w @x2

@2 w @y2

2

@/y @x

A1S22

A1S16



 2 iT @ w 2 @x@y dxdy

BS12

BS16

ð10Þ S

S

S

where A1 , A2 , B , D are given in Appendix A.

BS11

 2 i6 dw 6 BS 2 @@x@y 4 12 iT

dxdy 

BS16 Z h R

@d/x @x

BS12 BS22 BS26 @d/y @y

BS16

4. Approximating functions

3

7 BS26 7 5 BS66

@d/x @y

þ

@d/y @x

i

3  2 iT 6 S 7h @ 2 w @2 w S S 7 @ w 6 2 @x@y dxdy 4 B12 B22 B26 5 @x2 @y2 S S S B16 B26 B66 " # Z

 A2S44 A2S45

T d/y d/x /y /x dxdy þ S S R A245 A255 " 3 Z 2 2 h @w @dw 2h @w 2h @dw h  x2 q  2 d/x  2 /x þ /x d/x @x 2 12 @x @x p @x p R BS11

# 3 2 2 h @w @dw 2h @w 2h @dw h  2 d/y  2 /y þ /y d/y þ hwdw dxdy @y 2 12 @y @y p @y p

R

@ 2 dw @y2

þ

A1S12

Applying the mapping technique, an arbitrarily shaped quadrilateral plate in Cartesian coordinates (Fig. 1A) may be expressed by a square plate defined in the natural coordinates (Fig. 1B) by the simple boundary equations n = ± 1 and g = ± 1. The mapping of the Cartesian coordinate system is given by [35]:

x¼

4 X M i ðn; gÞxi i¼1

4 X y¼ M i ðn; gÞyi

ð11Þ

i¼1

where (xi, yi), i = 1, . . ., 4 are the coordinates of the four corners of the quadrilateral region R and Mi(n,g) are the interpolation functions of the serendipity family given by:

42

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

M i ðn; gÞ ¼

1 ð1 þ gi gÞð1 þ ni nÞ 4

ð12Þ

The transformation Eq. (11) maps a point (n, g) in the master ~ onto a point (x, y) in the real plate domain R and vice versa plate R if the Jacobian determinant of the transformation given by:

jJj ¼

@x @y @x @y  @n @ g @ g @n

ð13Þ

is positive. Applying the chain rule of differentiation it can be shown that the first derivatives of a function in both spaces are related by:

"

@ @x @ @y

#

" 1

¼J

@ @n @ @g

2

#

¼4

J 22 jJj

 Jj12Jj

 Jj21Jj

J 11 jJj

3" 5

@ @n @ @g

# ð14Þ

where J is the Jacobian given by:

J¼



J 11

J 12

J 21

J 22



P xi Mi;n ¼ P xi M i;g

P P

yi Mi;n

 ð15Þ

yi M i;g

and by Bardell et al. [38,39]. Briefly, the convergence in th h-p version of FEM is sought by simultaneously refining the mesh and increasing the degree of the elements. As has been demonstrated in a previous work [33] a very good convergence results can be obtained increasing the amount of hierarchy Gram–Schmidt orthogonal polynomials and using a single quadrilateral element. The unknown functions w, /x , /y in Eq. (10) are approximated by the product of the shape functions in natural coordinates n, g, by the respective generalized displacements: T

wðn; gÞ ¼ fNw ðn; gÞg fcw g T

/x ðn; gÞ ¼ fN/ ðn; gÞg fc/x g

@2 w @x2

3

2

@2 w @n2

3

7 h 6 2 7 h i6 i 6@ w7 ð1Þ 6 2 7 ð2Þ 6 @y2 7 ¼ Op 6 @@ gw2 7 þ Op 5 5 4 4

where: {cw}, fc/x g, fc/y g are the vectors of unknown generalized displacement and h ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ fNðÞ ðn; gÞg ¼ p1 q1 ; p1 q2 ; p1 q3 ; p1 q4 ; p1 q5 ; . . . ;

@2 w @x@y

@2 w @n@ g

ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ ðÞ

ðÞ ðÞ

 p1 qðÞ n ; p2 q1 ; p2 q2 ; p2 q3 ; p2 q4 ; p2 q5 ; . . . ;  p2 qðÞ n ; p3 q1 ; p3 q2 ; p3 q3 ; p3 q4 ; p3 q5 ; . . . ; ðÞ

 p3 qðÞ n ;...;

" @w # @n @w @g

ð17Þ

/y

/y ðn; gÞ ¼ fN ðn; gÞg fc g

The elemental area dxdy in the Cartesian domain R is transformed into jJjdndg. Applying again the chain rule of differentiation [37] in Eq. (14), results:

2

T

/

ðÞ

ðÞ

ðÞ

ðÞ

ðÞ

ðÞ

ðÞ ðÞ ðÞ ðÞ ðÞ ðÞ ¼ w; /x ; /y pðÞ n q1 ; pn q2 ; pn q3 ; pn q4 ; pn q5 ; . . . ; pn qn

ð16Þ

¼ NðÞ

where the elements of the matrices [Op(1)] and [Op(2)] are given in Appendix B. The h-p formulation developed in this article follows the general guidelines to those presented elsewhere by the authors [33]

are the vectors that include the shape functions. The first polynomials in N(), () = w, /x , /y are the Hermite w cubic polynomials pw i ðnÞ; qj ðgÞði; j ¼ 1 . . . 4Þ and Hermite linear / / polynomials p/i x ðnÞ; q/j x ðgÞ; pi y ðnÞ; qj y ðgÞði; j ¼ 1; 2Þ.

Table 1 qffiffiffiffi  ¼ x ah2 Eq of square plate (0/90/90/0). Non-dimensional first frequency x 2 a/h

E1/E2

5

Present macro element

5

10

20

30

40

6.482 6.397 6.328 6.306

8.490 8.467 8.330 8.319 8.421 8.270

10.035 10.022 9.834 9.828 9.671 9.528

11.007 10.983 10.773 10.769 10.416 10.279

11.718 11.675 11.458 11.454 10.938 10.773 10.807 10.729

7.240

9.915 9.876 9.806 9.777 9.912 9.847

12.458 12.438 12.308 12.292 12.316 12.225

14.261 14.250 14.064 14.054 13.943 13.987

15.653 15.644 15.412 15.405 15.213 15.112 15.101 15.166

GS = 2 GS = 3 GS = 4 GS = 5

7.522 7.515 7.512 7.503

10.619 10.612 10.608 10.600

13.879 13.873 13.864 13.859

16.487 16.480 16.465 16.460

18.716 18.710 18.688 18.684 18.659 18.902

GS = 2 GS = 3 GS = 4 GS = 5

7.533 7.530 7.530 7.527

10.644 10.640 10.639 10.637

13.934 13.929 13.927 13.926

16.578 16.573 16.569 16.568

18.851 18.845 18.840 18.839 18.822 19.157

GS = 2 GS = 3 GS = 4 GS = 5

7.536 7.536 7.536 7.536

10.651 10.650 10.650 10.650

13.950 13.948 13.948 13.948

16.609 16.604 16.604 16.604

18.896 18.891 18.891 18.891

GS = 2 GS = 3 GS = 4 GS = 5

Xiang and Wang [22] Liu et al. [41] Ferreira and Fasshauer [42] Zhen and Wanji [43] 10

Present macro element

6.557

GS = 2 GS = 3 GS = 4 GS = 5

Xiang and Wang [22] Liu et al. [41] Ferreira and Fasshauer [42] Zhen and Wanji [43] 50

Present macro element

7.197 7.115 7.083 7.029

Ferreira and Fasshauer [42] Zhen and Wanji [43] 100

Present macro element

Ferreira and Fasshauer [42] Zhen and Wanji [43] 1000

Present macro element

i

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

43

Fig. 3. Geometry of laminated plate (A) and mesh of macro elements (B).

Then, an adequate number of Gram–Schmidt polynomials are added to formulate a polynomially-enriched plate macro element: / / w pw ði; j ¼ 5 . . . nÞ and p/i x ðnÞ; q/j x ðgÞ; pi y ðnÞ; qj y ðgÞ i ðnÞ; qj ðgÞ ði; j ¼ 3 . . . mÞ. The degree of the first Gram–Schmidt polynomial in both natural coordinates is 4 for the transversal displacement (w) and 2 for the rotations (/x , /y ). The members of the set of characteristic orthogonal polynomials are obtained following the procedure described in references [33,36,37]. This way, the

hierarchical modes contribute only to the internal displacement of the element, and do not therefore affect the displacement along the element edges or at the element nodes. Nevertheless, products obtained between any of the Gram–Schmidt (GS) characteristics orthogonal polynomials and the Hermite polynomials will constitute what amounts to edge freedoms along the element boundaries. Replacing the expressions of Eq. (17) into Eq. (10) and applying the mapping technique (Eqs. (13)–(16)) results:

 ¼ x ah2 Fig. 4. Non-dimensional frequencies x

qffiffiffiffi q

E2

case 1.

44

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

 ¼ x ah2 Fig. 5. Non-dimensional frequencies x

2 3 2 3 0 0 0 h 0 A1 0 i  6 7 6 7 S 0¼ dc w dc /x dc /y 4 A1 0 A2 5 A1ij 4 0 0 A2 5 33 1 1 0 A2 A1 0 A2 A1 2 w3 2 3 c A3 A4 2A5 Z 1 Z 1

 6 7 6 7   4 c/x 5jJj dn dg þ 0 0 5 Dij 33 dc w dc /x dc /y 4 0 1 1 c/y 0 0 0 2 32 w 3 A3 0 0 c Z 1 Z 1

w  6 76 7  4 A4 0 0 54 c/x 5jJj dn dg þ dc dc /x dc /y 1 1 c/y 2A5 0 0 2 3 2 32 w 3 A3 A4 2A5 h 0 A1 0 c i 6 7 6 76 7 S 0 0 5 Bij  4 0 0 A2 54 c/x 5jJj dn dg 4 0 3 3 /y c 0 0 0 0 A2 A1 2 3 0 0 0 h Z 1 Z 1 i

w 6 7 S þ dc dc /x dc /y 4 A1 0 A2 5 Bij  3 3 1 1 0 A2 A1 2 32 w 3 A3 0 0 c Z 1 Z 1

w  6 76 7  4 A4 0 0 54 c/x 5jJj dn dg þ dc dc /x dc /y 1 1 c/y 2A5 0 0 2 3 2 3 " # cw 0 0 h i 0 0 N/ 6 /x 7 6 S /7  4 0 N 5 A2ij  4 c 5jJj dn dg  x2 qhfdcg 2 2 0 N/ 0 / N 0 c/y Z

1

Z

1

qffiffiffiffi q

E2

case 2.

h i 2 2 2 h2  p2h2 N/ A10 ðNw Þ þ 12 ðA10 Þ þ ðA20 Þ 1 6 6 2 1 6  p2h2 N/ A10 ðN/ Þ 2 1 4 2





Z

1

1

Z

 p2h2 N/ A20

0

 p2h2 N/ A20 0 2 1 ðN/ Þ 2

3 7 7 7 5

 jJj dn dgfdcg

where:

J 22 @N/ J 12 @N/ J @N/ J 11 @N/  ; A2 ¼  21 þ jJ j @n jJ j @ g jJ j @n jJ j @ g  w w 3 2 w 2 w 2 w X @ N @ N @ N @N 0 @N 0 0 0 A3 ¼ a01 þ a02  a þ a a þ b 3 i i i 2 @ g2 @n@ g i¼1 @n @g @n  w w 3 2 w 2 w 2 w X 0 @ N 0 @ N 0 @ N 0 0 @N 0 @N þ b  b þ b a þ b A4 ¼ b1 i 2 3 @ g2 @n@ g i¼1 i i @n @g @n2  w w 3 2 w 2 w 2 w X @ N 0 @N 0 @ N 0 @ N 0 0 @N A5 ¼ c01  c þ c  c a þ b i 2 3 @ g2 @n@ g i¼1 i i @n @g @n2 w w J @N J @N  12 A10 ¼ 22 jJ j @n jJ j @ g J 21 @Nw J 11 @Nw 0 A2 ¼  þ jJ j @n jJ j @ g A1 ¼

ð19Þ

45

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

 ¼ x ah2 Fig. 6. Non-dimensional frequencies x

Finally, cancelling the virtual nodal displacement, Eq. (13) results: E

E

2

E

f½K   x ½M gfc g ¼ 0

ð20Þ

E

qffiffiffiffi q

E2

case 3.

h i 3 2 2 2 h2 ðNw Þ þ 12 ðA10 Þ þ ðA20 Þ  p2h2 N/ A10  p2h2 N/ A20 7 16 6 7 2 1 6 7jJjdndg  p2h2 N/ A10 ðN/ Þ 0 2 5 1 4 2

h

i M E ¼ qh

Z

1 Z

1

 p2h2 N/ A20

where [K ] is the stiffness matrix of mapped macro element according to TSDT and it is given by:

Z

1

Z

1



½B1½A1S ½B1T þ ½B3½D½B3T þ ½B3½BS ½B1T þ½B1½BS ½B3T þ ½B2½A2S ½B2T jJjdndg

½K E  ¼

1

0 6 ½B1 ¼ 4 A1 2

0

0 0

3 0 7 A2 5;

A2 A1

A3 A4 2A5

6 ½B3 ¼ 4 0

0

0

0

E

2 1 ðN/ Þ 2

ð23Þ

5. Assembly and global equations

1

ð21Þ

with:

2

0

3

2

0 6 ½B2 ¼ 4 0 N/

3 0 7 N/ 5; 0

ð22Þ

7 0 5 0

and [M ] is the mass matrix of macro element in natural coordinates given by:

For the purposes of modeling plates with more general planform, it is here necessary to address how to combine the individual macro finite elements obtained in the previous section. The approach adopted here follows the general outlines of previous works [33,40]. In order to obtain the global matrices, the additional complexity generated by the use of hierarchical modes has to be taken into account. So, it is necessary to identify in the element stiffness matrix and in the element mass matrix (Eqs. (21) and (23)) the row/column corresponding to particular degrees of freedom. These are nodal (N), edge (E) or purely internal (I) degrees of freedom. After the identification is carried out, they are separated and then rearranged according to the following equation [40]:

46

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

Fig. 7. Geometry of triangular plate with an internal hole and adopted mesh.

2

½NN ½NE ½NI

6 ½K  or ½M  ¼ 4 ½EN ½IN E

E

½EE ½IE

3

7 ½EI 5 ½II

ð24Þ

Each {cE} is formed by sub vectors {cw}, {c/x}, {c/y} that allow obtaining the displacement functions of the whole macro element. 6. Numerical results

In Eq. (24) N is associated to Hermite polynomials, E is associated to hierarchical GS and Hermite polynomials, finally I is associated to hierarchical GS polynomials. The assembly is carried out by equating generalized displacements of two adjacent macro elements (MEi, with i: number of elements). The assembly procedure is schematically shown in Fig. 2, where the mesh is obtained using the lower number of quadrilateral macro elements. Inter-element compatibility is achieved simply by matching the appropriate generalized coordinates at common element nodes and along common edges (Fig. 2). This methodology ensures the elements are fully conforming and leads to the global stiffness and mass matrices, [KG] and [MG]. Different boundary conditions may be applied to the laminated plate assembly, removing from the global [KG] and [MG] matrices, the rows and columns which correspond to the degrees of freedom associated with the corresponding support condition. Assuming simple harmonic motion, the governing equation of motion for free vibration can be obtained as:

f½K G   x2 ½M G gfcG g ¼ f0g G

ð25Þ

where {c } is the global vector of generalized nodal displacement, that allows obtaining the vector {cE} (Eq. (20)) of each macro element of the adopted mesh, from correct vectors of indexation.

6.1. Convergence study The formulation developed in this work has been implemented in a computer program that allows working with laminated anisotropic plates (built by a single layer or symmetric laminates), with different combinations of boundary conditions, material properties, thickness ratios, and also to assembly two or more macro elements in order to include the analysis of plates with complex geometry. The convergence study and verification of results are presented for the fundamental frequency corresponding to a square fourlayer laminated plate with stacking sequence (0/90/90/0) and with the four edges simply-supported. The material properties of each lamina are G12 = G13 = 0.6E2, G23 = 0.5E2, m12 = 0.25. The non-dimensional fundamental frequencies for five thickness ratios a/h and for different orthotropic relations E1/E2 are depicted in Table 1. The values obtained are summarized for the first frequency, using from two (GS = 2) to five (GS = 5) polynomials of Gram–Schmidt in a single macro element, and it is shown the references values of Liu et al. [41] and Xiang and Wang [22], who employed a higher order theory and multiquadratic radial basis function, Ferreira and Fasshauer [42] who presented results of free vibration frequencies using an RBF-pseudospectral method together with FSDT and Zhen and Wanji [43] who used global–local higher order theory.

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

 ¼ x ah2 Fig. 8. Non-dimensional frequencies x

It can be seen that from the use GS = 3 Gram–Schmidt polynomials, the corresponding frequencies have a tendency to stabilize, showing also a stable convergence without oscillations, and with a great concordance with the results of others authors. Furthermore at a fixed length-to-thickness ratio, the dimensionless fundamental frequency increases when increasing the orthotropic ratio for these boundary conditions. A similar behavior can be observed when increasing the length-to-thickness ratio but by fixing the orthotropic ratio. Clearly, Table 1 demonstrates that the computer results with the proposed finite macro element are in any case well in agreement with the solutions reported by others researchers, whether thin or thick laminates are concerned. In fact results so far suggested that this formulation is insensitive to shear locking. 6.2. Results for thick laminated plate assemblies To show the applicability of the proposed formulation for studying thick laminated plates with arbitrary geometrical

qffiffiffiffi q

E2

47

of triangular plate.

planforms and general boundaries conditions, examples of assembly of macro elements are presented. The results shown here and in section 6.3 for thick plate structures with irregular planform can serve as a benchmark comparison for further investigations. In addition they provide clear evidence of the versatility and capability of the current methodology. Vibration analyses are presented here for two laminated plate structures with mixed boundary conditions; trapezoidal plates with different internal continuity conditions and a triangular plate with an internal hole. For each plate, the first six natural frequencies and modes shapes are presented. The analyzed plates are symmetric three-layers laminates with the following material properties of each lamina: E1 = 60.7 GPa, E2 = 24.8 GPa, G12 = 12 GPa, G13 = 0.5E2, G23 = 0.2E2, m12 = 0.23. The first example concerns the trapezoidal plate shown in Fig. 3A, with stacking sequence (60/30/60). The plate is clamped along the x axis, while the other edges are free, and the thickness ratio is a/h = 10. This laminate can be studied using just one finite

48

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

macro element. However, in order to showing the potentiality of the present formulation, the assembly of the two macro elements (Fig. 3B) is considered and used to study three different cases: Case 1: plate with an interface within elements that require continuity of displacement and rotations. Case 2: plate with internal discontinuity (along the shared side of the macro elements of the mesh) and joint point in x = 0.5 a, y = 0.6 a. Case 3: plate with a simply support line coinciding with the common side of the elements of the mesh. The first six non-dimensional frequencies of free vibration of these three cases are presented in Figs. 4–6, respectively, where the associated mode shapes have been also included. The second example shows the application of a simple mesh, with a minimum number of macro elements, for the dynamic study of a triangular plate with an internal triangular hole. The plate and the adopted mesh are presented in Fig. 7. The stacking sequence is (45/-45/45) and thickness ratio is a/h = 10. The plate is simply supported along one side, free along another one and clamped along the last one, and the hole-sides are all free (Fig. 7). The non-dimensional frequencies and the corresponding mode shapes are presented in Fig. 8. 7. Conclusions An enriched plate macro element, based on Trigonometric Shear Deformation Theory is developed in this work. The formulation is implemented in a computer program to carry out free vibration analysis of composite structures modeled as plate assemblies. The proposed macro element is a significant refinement over previously developed macro elements using classical and first order shear deformation theories. The present formulation is particularly useful when analyzing thick composite plates and is free from shear locking phenomenon, because of the use of higher order macro elements in conjunction with higher order shear deformation theory. This represents an important difference with the employment of the classical finite element method. The macro element has been tested through numerical examples, producing very good results. It is possible to achieve very good accuracy in the results using a low number of polynomials and a mesh with a minimum number of macro elements in plates with general planform. The assembly examples provide evidence of the versatility and capability of the current method. Acknowledgments

      @/y @/x @d/x @d/y @d/x @d/y þ A1S26 þ þ @x @y @x @y @y @x    @/x @/y @d/x @d/y S dxdy þ þ þA166 @y @x @y @x " ! ! ! ! Z @2w @ 2 dw @2w @ 2 dw þ D11 þ D12 @x2 @x2 @y2 @x2 R ! ! ! ! @2w @ 2 dw @2w @ 2 dw þ2D16 þ D12 @x@y @x2 @x2 @y2 ! ! ! ! 2 2 2 @ w @ dw @ w @ 2 dw þ D22 þ 2D26 @y2 @y2 @x@y @y2 ! ! ! ! @2w @ 2 dw @2w @ 2 dw þ 2D26 þ2D16 @x2 @x@y @y2 @x@y ! !# @2w @ 2 dw dxdy þ4D66 @x@y @x@y ! !   Z " @ 2 w @d/x @ 2 w @d/x S  þ B BS11 12 @x2 @y2 @x @x R ! !   2 2 @ w @d/x @ w @d/y þ BS12 þ2BS16 @x@y @x2 @x @y ! !   2 2 @d/y @ w @d/y @ w þ 2BS26 þ BS22 2 @y @x@y @y @y ! !   2 2 @ w @d/x @d/y @ w @d/x @d/y S S þ B26 þ þ þ B16 2 2 @x @y @y @x @y @x ! #  2 ! Z "  2 @ w @d/x @d/y @/x @ dw þ BS11 þ2BS66 dxdy  @x@y @x2 @y @x @x R   2 !   2 ! @/y @ dw @/x @/y @ dw þ BS16 þ þBS12 @x2 @x2 @y @y @x   2 !   2 ! @/y @/x @ dw @ dw þ BS22 þ BS12 @y2 @y2 @x @y   2 !   2 ! @/x @/y @ dw @/x @ dw S þ 2B þ þ BS26 16 @y2 @x@y @y @x @x   2 !   2 !# @/y @ dw @/x @/y @ dw þ2BS66 dxdy þ þ 2BS26 @x@y @x@y @y @y @x Z h þ A2S44 ð/y Þðd/y Þ þ A2S45 ð/y Þðd/x Þ þ A2S45 ð/x Þðd/y Þ R " 3 Z 2 i h @w @dw 2h @w S 2  2 d/ þA255 ð/x Þðd/x Þ dxdy  x q 12 @x @x p @x x R

þ A1S16

The present investigation has been sponsored by the CONICET Project PIP 0105/2010, CIUNSA Project 1903 and REDES V Project, Supported by SPU. Appendix A The integration of Eq. (11) throughout the thickness (z coordinate) is:

      Z  @/y @d/x @/x @d/x þ A1S12 A1S11 @x @x @y @x R       @d/y @/x @/y @d/x @/x þA1S16 þ A1S12 þ @y @x @x @x @y       @/ @d/ @/ @d/ @/ y y y y x þ A1S22 þ A1S26 þ @y @y @y @x @y

0¼

2h

2

3

2

@dw h h @w @dw 2h @w þ /x d/x þ  2 d/ @x 2 p2 12 @y @y p @y y # 2 2h @dw h  2 /y þ /y d/y þ hwdw dxdy @y 2 p 

/x

where now the stiffness are given by:

pz p 2 Z h=2 pz 2 Q ij sin dz Q ij cos2 ; A2Sij ¼ dz h h h h=2 h=2 Z h=2 Z h=2 pz Q ij z sin Q ij z2 dz dz; Dij ¼ BSij ¼ h h=2 h=2

A1Sij ¼

Z

h=2

Appendix B The matrices [Op(1)] and [Op(2)] that appear in Eq. (18) are as follow:

R.F. Rango et al. / Composite Structures 119 (2015) 38–49

2

2 0 a1 h i 6 0 Opð1Þ ¼ 4 b1 c01

3 X a0 a0 6 6 i¼1 i i 3 6 a03 3 h i 6 6 X 0 7 0 ð2Þ b3 5 Op ¼6 bi a0i 6 6 i¼1 0 c3 6 6 X 3 4  c0i a0i

a02 0

b2 c02

i¼1

3

3 X a0i b0i 7 7 i¼1 7 3 X 0 0 7 7 bi bi 7 7 7 i¼1 7 3 X 07 5 0  ci bi i¼1

where:

a01 ¼

0

b1 ¼

c01 ¼

J 222 jJj

2

J 221 jJj

2

;

a02 ¼

;

b2 ¼

J 21 J 22 jJj2

0

;

J 212 jJj2 J 211 jJj2

c02 ¼

;

a03 ¼ 2

;

b3 ¼ 2

J 11 J 12 jJj2

J 11;n J 22 þ J 12;n J 21 jJj J J  J J a03 ¼ 11;g 22 22;n 21 jJj

a01 ¼

0

;

J 12 J 22 jJj2 J 11 J 21 jJj2

c03 ¼

a02 ¼

J 11 J 22 þ J 12 J 21 jJj2

J 21;g J 22 þ J 22;g J 21 ; jJj

J 21;g J 12  J22;g J 11 J 11;n J 12  J 12;n J 11 ; b02 ¼ ; jJj jJj J 11;g J 12 þ J 22;n J 11 b03 ¼ jJj

b01 ¼

References [1] Reddy JN, Liu CF. A higher-order shear deformation theory of laminated elastic shells. Int J Eng Sci 1985;23:319–30. [2] Touratier M. An efficient standard plate theory. Int J Eng Sci 1991;29(8):901–16. [3] Soldatos KP. A transverse shear deformation theory for homogeneous monoclinic plates. Acta Mech 1992;94:195–220. [4] Kant T, Swaminathan K. Analytical solutions for the static analysis of laminated composite and sandwich plates based on a higher order refined theory. Compos Struct 2002;56:329–44. [5] Karama, M, Afaq, KS, Mistou, S. A new theory for laminated composite plates. In: Proc IMechE, vol.223. Part L: J. Mater: Des Appl., 2009. [6] Mantari JL, Oktem AS, Guedes Soares C. Static and dynamic analysis of laminated composite and sandwich plates and shells by using a new higherorder shear deformation theory. Compos Struct 2011;94:37–49. [7] Mantari JL, Guedes Soares C. Bending analysis of thick exponentially graded plates using a new trigonometric higher order shear deformation theory. Compos Struct 2012;94:1991–2000. [8] Grover N, Singh BN, Maiti DK. Analytical and finite element modeling of laminated composite and sandwich plates: an assessment of a new shear deformation theory for free vibration response. Int J Mech Sci 2013;67:89–99. [9] Carrera E. Theories and finite elements for multilayered plates and shells: a unified compact formulation with numerical assessment and benchmarks. Arch Comput Meth Eng 2003;10:215–96. [10] Carrera E, Brischetto S, Cinefra M, Soave M. Effects of thickness stretching in functionally graded plates and shells. Compos: Part B 2011;42:123–33. [11] Demasi L. Mixed plate theories based on the generalized unified formulation. Part I: governing equations. Compos Struct 2009;87:1–11. [12] Demasi L. Mixed plate theories based on the generalized unified formulation. Part II: layerwise theories. Compos Struct 2009;87:12–22. [13] Demasi L. Mixed plate theories based on the generalized unified formulation. Part III: advanced mixed high order shear deformation theories. Compos Struct 2009;87:183–94. [14] Demasi L. Mixed plate theories based on the generalized unified formulation. Part IV: zig – zag theories. Compos Struct 2009;87:195–205.

49

[15] Demasi L. Mixed plate theories based on the generalized unified formulation. Part V: results. Compos Struct 2009;88:1–16. [16] Mantari JL, Oktem AS, Guedes Soares C. Bending and free vibration analysis of isotropic and multilayered plates and shells by using a new accurate higherorder shear deformation theory. Compos: Part B 2012;43(8):3348–60. [17] Mantari JL, Oktem AS, Guedes Soares C. A new trigonometric shear deformation theory for isotropic, laminated composite and sandwich plates. Int J Solids Struct 2012;49:43–53. [18] Shimpi RP, Ainapure AV. A beam finite element based on layerwise trigonometric shear deformation theory. Compos Struct 2001;53:153–62. [19] Arya H, Shimpi RP, Naik NK. A zigzag model for laminated composite beams. Compos Struct 2002;56:21–4. [20] Ferreira AJM, Roque CMC, Jorge RMN. Analysis of composite plates by trigonometric shear deformation theory and multiquadrics. Comput Struct 2005;83:2225–37. [21] Roque CMC, Ferreira AJM, Jorge RMN. Modelling of composite and sandwich plates by a trigonometric layerwise deformation theory and radial basis functions. Compos: Part B 2005;36:559–72. [22] Xiang S, Wang KM. Free vibration analysis of symmetric laminated composite plates by trigonometric shear deformation theory and inverse multiquadric RBF. Thin-Walled Struct 2009;47:304–10. [23] Mantari JL, Oktem AS, Guedes Soares C. A new trigonometric layerwise shear deformation theory for the finite element analysis of laminated composite and sandwich plates. Comput Struct 2012:45–53. [24] Mantari JL, Guedes Soares C. A trigonometric plate theory with 5-unknowns and stretching effect for advanced composite plates. Compos Struct 2014;107:396–405. [25] Mantari JL, Guedes Soares C. Four-unknown quasi-3D shear deformation theory for advanced composite plates. Compos Struct 2014;109:231–9. [26] Ramesh SS, Wang CM, Reddy JN, Ang KK. A higher-order plate element for accurate prediction of interlaminar stresses in laminated composite plates. Compos Struct 2009;91:337–57. [27] Zamani M, Fallah A, Aghdam MM. Free vibration analysis of moderately thick trapezoidal symmetrically laminated plates with various combinations of boundary conditions. European J Mech A/Solids 2012;36:204–12. [28] Fazzolari FA, Boscolo M, Banerjee JR. An exact dynamic stiffness element using a higher order shear deformation theory for free vibration analysis of composite plate assemblies. Compos Struct 2013;96:262–78. [29] Houmat A, Rashid MM. Coupling of h and o finite elements: application to free vibration analysis of plates with curvilinear plan-forms. Appl Math Model 2012;36:505–20. [30] Rango RF, Nallim LG, Oller S. Formulation and assembly of hierarchical finite elements to the static and dynamic analysis of laminated quadrilateral plates. Revista Sul-Americana de Engenharia Estrutural 2012;9:4–21. [31] Rango RF, Nallim LG, Oller S. Static analysis of thick laminated plates using enriched macro elements. In: Proceedings of the 16th International Conference on Composite Structures – ICCS. Portugal, 2011. [32] Rango RF, Bellomo FJ, Nallim LG. A general Ritz algorithm for the static analysis of arbitrarily laminated composite plates using first order shear deformation theory. J Eng Res (TJER) 2013;10(2):1–12. [33] Rango RF, Nallim LG, Oller S. Static and dynamic analysis of thick laminated plates using enriched macroelements. Compos Struct 2013;101:94–103. [34] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. 2nd ed. Boca Raton, USA: CRC Press; 2003. [35] Zienkiewicz OC, Taylor RL. The finite element method for solid and structural mechanics, 6th ed., vol. 2. Great Britain: Elsevier; 2005. [36] Nallim LG, Oller S, Grossi RO. Statical and dynamical behaviour of thin fibre reinforced composite laminates with different shapes. Comput Meth Appl Mech Eng 2005;194:1797–822. [37] Nallim LG, Oller S. An analytical–numerical approach to simulate the dynamic behaviour of arbitrarily laminated composite plate. Compos Struct J 2008;85:311–25. [38] Bardell NS, Dunsdon JM, Langley RS. Free vibration analysis of thin rectangular laminated plate assemblies using the h-p version of the finite element method. Compos Struct J 1995;32:237–46. [39] Bardell NS, Dunsdon JM, Langley RS. Free vibration analysis of coplanar sandwich panels. Compos Struct J 1997;38(1-4):463–75. [40] Bardell NS, Dunsdon JM, Langley RS. Free vibration analysis of thin coplanar rectangular plate assemblies – Part I: theory, and initial results for specially orthotropic plates. Compos Struct J 1996;34:129–43. [41] Liu L, Chua LP, Ghista DN. Mesh-free radial basis function method for static, free vibration and buckling analysis of shear deformable composite laminates. Compos Struct 2007;78:58–69. [42] Ferreira AJM, Fasshauer GE. Analysis of natural frequencies of composite plates by an RBF-pseudospectral method. Compos Struct 2007;79(2):202–10. [43] Zhen W, Wanji C. Free vibration of laminated composite and sandwich plates using global-local higher-order theory. J Sound Vib 2006;298:333–49.