Through-the-thickness stress profiles in laminated composite and sandwich structure plates via unified formulation

Composites Part B 107 (2016) 29e42

Contents lists available at ScienceDirect

Composites Part B journal homepage: www.elsevier.com/locate/compositesb

Through-the-thickness stress profiles in laminated composite and sandwich structure plates via unified formulation Mauricio F. Caliri Jr. a, *, Antonio J.M. Ferreira b, Volnei Tita c a

Engineering Department, The University of Lavras, Mailbox 3037, Lavras, MG, Brazil Mechanical Engineering Department, Faculty of Engineering of Porto, University of Porto, Rua Dr Roberto Frias, Porto, Portugal c ~o Carlos School of Engineering, University of Sa ~o Paulo, Av. Joa ~o Dagnone, 1100, Sa ~o Carlos, SP, Brazil Aeronautical Engineering Department, Sa b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 31 August 2016 Received in revised form 13 September 2016 Accepted 17 September 2016 Available online 20 September 2016

This paper presents a discussion on the through-the-thickness stress profiles of thick laminated plates and sandwich structures calculated with a particular generalized unified formulation. It consists on a finite element Unified Formulation (UF) with C-1 continuity of the transverse displacement field implemented in a four node plate element (CGF- Caliri's Generalized Formulation). The current works on unified formulations attempt to simplify and settle the development, use and choice of plate and shell theories. However, even using unified approaches, no generalizations base on specific problems are possible. For instance, the present work demonstrates that a non-linear plate theory may yield worse results when compared to a linear solution. It was found that, for a selected thick sandwich structure, a cubic layerwise theory provided a set of results 27% less accurate than the linear one. Moreover, when the through-the-thickness stress profiles are closely evaluated, the accuracy of the stress profiles are subjected to the solution method and not only the theory. Particularly, the effects of poor mesh resolution and boundary conditions are demonstrated in this work. Nevertheless, CGF is able to predict pointwise accurate results (less than 1% of deviation) even for the for transverse stresses variables. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Unified formulations Sandwich structures Laminated composites Plates Finite element method CGF

1. Introduction Laminated composite and sandwich structural plates and shells comprise a class of structural components, which is applied on a daily basis in engineering. The structural simplification of 3D (three dimensional) structures into 2D (bi-dimensional) ones is very appealing, but not accomplished without overcoming a few hurdles. For example, wind turbine blade engineering is a top-notch application comprising such structures [1e3]. When designing these blades, the proper calculation of the internal stresses of these layered structures, especially those using thick laminates and sandwich structures [4] is crucial to the performance and life of the wind turbine. The simplification mentioned above is basically the removal of the thickness coordinate out of the problem, yielding the so-called plate/shell theories. Depending on how this simplification is performed, a different theory may emerge. There are, according to Carrera [5], three main approaches: continuum base, asymptotic or axiomatic theories. Each of these structural theories

* Corresponding author. E-mail address: [email protected]fla.br (M.F. Caliri). http://dx.doi.org/10.1016/j.compositesb.2016.09.055 1359-8368/© 2016 Elsevier Ltd. All rights reserved.

brings a different formulation to accurately estimate the displacement and/or stress fields. However, it is not a straightforward task to identify the origin of the lack of accuracy of a particular theory when compared to another one derived from a different approach. That is why the recent use Unified Formulations is very attractive [6]. By choosing one unified solution method, different plate/shell theories can be properly compared and tested, because the source of the variations in the results are controlled by the engineer [5e20]. The authors also recommend their lengthy review on plate and shell theories in Ref. [21]. It demonstrates, with over 100 papers, the actual infinite number of approaches to the plate/shell problem available in the literature. Thus, this paper presents the stress profiles obtained by the authors using a finite element based on unified formulation tagged as CGF (Caliri's Generalized Formulation) [22,23]. In fact, CGF is based on Carrera's Unified Formulation [5,6] and the Generalized Unified Formulation from Ref. [11]. The novelty of this solution method is the proposal of a unified finite element (4-node plate element) solution considering C-1 continuity of the transverse displacement field to mainly increase the accuracy of the transverse stresses without using mixed formulations. In the previous works

30

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

uk ðx; y; zÞ ¼ Ftk ðzÞukt ðx; yÞ t ¼ 0…Nt

Fig. 1. General form of acronyms used with unified formulations.

(2)

The Ft(z) variables are the transverse functions; the ut(x,y) are the in-plane functions. The subindex t is a recursive index for the current order of each component of the complete polynomial of order Nt assumed for the displacement variable u.

Table 1 Material properties of case study II. Ply

E1

E2

E3

G12

G13

G23

y12

y13

y23

Ply Core

25E2 0.04E2

1 GPa 0.04E2

1 GPa 0.5E2

0.5E2 0.016E2

0.5E2 0.06E2

0.2E2 0.06E2

0.25 0.25

0.25 0.1

0.25 0.1

of the authors [22,23] thin plates were investigated, as well as thick laminated and sandwich structure plates were studied after an improvement of the formulation. To appreciate the performance of CGF, regarding the stress profiles, several laminated composite plates are investigated and compared, considering the combination of different length-to-thickness ratios (a/h) and plate/shell theories. Additionally, two sandwich structures are investigated separately due to their higher material anisotropy and/or higher geometric aspect ratios.

uk ðx; y; zÞ ¼ Faku ðzÞukau ðx; yÞ; vk ðx; y; zÞ ¼ Fakv ðzÞvkav ðx; yÞ; wk ðx; y; zÞ ¼ Fakw ðzÞwkaw ðx; yÞ;

au ¼ 0…Nau av ¼ 0…Nav aw ¼ 0…Naw

(3)

In Ref. [11], Demasi generalized the unification process and decoupled the number of terms from each displacement field in equation (2) which produced equation (3). In Nau , for example, the subscript a has another subscript u to distinguish the order of the polynomials according to each displacement field. In the case of a laminated structure approached as one single monocoque ply (Equivalent Single Layer approach - ESL), the transversal functions F(z) can assume the form of the terms of a Taylor expansion, yielding:

2. The unified formulation The results obtained by CGF can be different from those obtained by Carrera and Demasi [5e8,10,11] (named CUF - Carrera's Unified Formulation and GUF - Generalized Unified Formulation) because the proposed finite element solution preserves the C-1 continuity of the transverse displacement field w(x,y,z) using Hermite polynomials. Unified formulations possess a kernel matrix from which the numerous theories stem. In CGF's kernel matrix, the finite element formulation yields in-plane rotations and a twist, which appear in the C-1 formulation process. The elementary equilibrium system for a steady-state case using CGF can be seen in equation (1).




P R

ðeÞ

 ¼

i

Kuu sym

Kuv Kvv

ðeÞ
ij

u w;

ðeÞ

 þ

j

Muu sym

0 Mvv

ðeÞ
ij

 u w€ ;

ðeÞ

FNk t ðzÞ ¼ FNt ðzÞ F0 ðzÞ ¼ 1 F1 ðzÞ ¼ z F2 ðzÞ ¼ z2 « FNt ðzÞ ¼ zNt

(4)

However, if a set of displacement functions is assumed to each ply, then the approach is termed Layerwise (LW) and the transverse functions must satisfy the compatibility requirements of the calculated variables at the interfaces. Carrera and Demasi [5e8,11] pointed to the use of Legendre polynomials to fulfill such task. In Ref. [23] this approach is also detailed. Thus, considering both ESL and LW cases, the finite element kernel of CGF, KCGF, is:

j

"

(1) Respectively, the upper case P, R and M variables stand for normal and moment loads and the mass. The variable K represents the matrix stiffness components. The sub-indices u and w represent the displacement fields and the curvatures, respectively, while the i and j are recursive indexes. The stiffness matrices in equation (1) are not in their unified form yet and neither of the unification or generalization parameters can be seen in this equation. To obtain an explicit unified kernel matrix, an axiomatic unification of the displacement fields needs to be equated into equation (1). Carrera's and Demasi's kernels stem from the respective unifications in equations (2) and (3). In these equations, the displacement fields are postulated. The superscript k refers to a particular ply of the layered structure in cases where the formulation is of the layerwise type [6,11]. Table 2 Material properties of case study III. Ply

E (MPa)

v

Top skin Core Bottom skin

8000 1000; 0.1 10000

0.34

KCGF ¼

~ krs K uu sym

~ krs K uw ~ krs K

# (5)

ww

where

0 ~ krs K uu

B ¼@

kru su Kuu

kru sv Kuv kru sv Kvv

0 sym kru sw K kru;xsw Kuw ;y B uw krs kru sw kru sw ~ B Kuw ¼ @ Kvw;x Kvw ;y kru sw kru sw Kww Kww ;x ;y 0 kru sw kru sw Kw K w w ;x w;y B ;x ;x kru sw ~ krs ¼ B K K ww w;y w;y @ sym

1 kru sw Kuw kru sw C A Kvw kru sw Kww 1 kru sw Kuw ;xy C kru sw C Kvw ;xy A kru sw Kww ;xy kru sw Kw ;x w;xy

(6) 1

C kru sw C Kw ;y w;xy A kru sw Kw ;xy w;xy

where the r,s,i and j letters are recursive indexes which come from the derivation of the finite element equilibrium equations based on the Principle of Virtual Displacements (PVD). krsij As an example, the Kuw component is shown in equation (7).

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42 Table 3 Coordinates of points used to evaluate the stress values.

x y z

Table 5 Comparison of shear stresses obtained with CGF and the analytical reference.

sxx

syy

szz

txz

tyz

txy

a/2 a/2 h/2

a/2 a/2 h/4

a/2 a/2 þh/2

0 a/2 0

a/2 0 0

a 0 h/2

kru sw;z

krsij Kuw ¼ Z13

Theory

ðzÞ9Ni;x ðx; yÞHj1 ðx; yÞ8þ

kru sw;z

ðzÞ9Ni;y ðx; yÞHj1 ðx; yÞ8þ kru;z sw 1 Z55 ðzÞ9Ni ðx; yÞHj;x ðx; yÞ8þ kru;z sw 1 Z54 ðzÞ9Ni ðx; yÞHj;y ðx; yÞ8

Z63

(7)

The triangles 9 8 indicate the in-plane integration, N are the linear Serendipity's interpolation polynomials whereas H are Herkru sw mite's. Also, Z13 ;z term in equation (8) refers to the thickness integration drill.

kru sw;z

Z13

k

k Zzt

¼ C 13

Fru ðzÞFsw;z ðzÞ dz

(8)

zkb k

The C 13 is the regular rotated orthotropic material constant. The full derivation, along with the 21 explicit components of this elementary symmetrical kernel, can be appreciated in Refs. [22,23] where it was introduced. Carrera and Demasi proposed two sets of acronyms which are summarized in Fig. 1. There are four blocks in this picture. The first one indicates whether the formulation is of the ESL or LW type. The second block tags whether the solution was derived via PVD or via a Mixed formulation. Third block refers to the interlaminar continuity constraint, C. Last block indicates, with the subscript Nui zi and the superscript Nsi zi , whether there are zig-zag refinements in a particular direction and what are the order of the polynomials used for each variable, respectively. However, in this paper, only the derivations via PVD are considered. Because of this, the tags are restricted to the displacement (D) tags only. Thus, an ED113 acronym means that the formulation is of the Equivalent Single Layer type, derived from the PVD method and does not account for interlaminar continuity nor zig-zag refinements. Lastly, the 113 number points to linear inplane displacement fields and a cubic one for the transverse displacement.

Table 4 Comparison of normal stresses obtained with CGF and the analytical reference. Theory

sxx

110 111 112 113 220 221 222 223 330 331 332 333 Abaqus 2D Abaqus 3D Analytical [24]

43.04 n/a 43.04 10.42 42.72 11.51 42.72 11.67 43.04 n/a 43.04 2.32 42.72 3.50 42.72 3.68 3.60 n/a 3.60 2.18 3.76 3.81 3.76 3.99 0.3705 0.7054 0.7200

ESL

syy LW

31

szz

Value

ESL

LW

ESL

LW

18.87 18.87 16.79 16.79 18.87 18.87 16.79 16.79 7.21 7.21 5.58 5.58 0.6593 0.6795 0.6630

n/a 10.11 8.91 8.57 n/a 2.58 1.27 1.07 n/a 1.83 0.30 0.11

73.10 73.10 125.48 125.48 73.10 73.10 125.48 125.48 58.41 58.41 102.70 102.70 n/a 0.9721 1.0000

n/a 0.59 14.88 16.98 n/a 3.61 11.28 17.03 n/a 3.84 10.26 16.00

%

Scaled

110 111 112 113 220 221 222 223 330 331 332 333 Abaqus 2D Abaqus 3D Analytical [24]

tyz

txz

txy

Value

ESL

LW

ESL

LW

ESL

LW

44.42 44.42 43.60 43.60 44.42 44.42 43.60 43.60 14.52 14.52 14.25 14.25 0 0.1134 0.292

n/a 35.34 31.23 30.41 n/a 3.53 0.96 1.85 n/a 11.95 9.18 8.22

45.11 45.11 42.37 42.37 45.11 45.11 42.37 42.37 2.74 2.74 0.59 0.59 0 0.0747 0.219

n/a 17.12 15.30 15.34 n/a 12.65 10.73 10.55 n/a 12.56 11.42 11.14

80.73 80.73 79.87 79.87 80.73 80.73 79.87 79.87 58.03 58.03 58.46 58.46 0.0083 0.0326 0.0467

n/a 34.90 33.62 32.12 n/a 4.93 10.28 12.42 n/a 12.42 20.99 23.34

%

Scaled

In Refs. [22,23], Caliri et al. focused on the introduction of the CGF solution method and showed the displacement and natural frequency results obtained with this formulation. The performance of the CGF solution method, regarding the magnitude of the stress in laminated structures and sandwich plates, is the subject of this paper. The same in-house FE code developed within MATLAB in Refs. [22,23] is used but a post-processing step of the displacement fields, to evaluate the stress profiles within the laminated structures, is investigated. To assess the accuracy of the proposed solution method, the commercial software Abaqus was employed to generate baseline results. This FE baseline comprises the results obtained from finite element models built with shell (S4) and brick elements (C3D8, C3D20). The validation herein proposed is divided into three case study sections. The first case analyzes symmetric laminated plates of different thicknesses with different plate/shell theories. The second case is a symmetric laminated sandwich structure. Lastly, the third case investigates a non-symmetric sandwich structure made of isotropic plies. 3. Results and discussions In this section, three case studies (I, II and III) are considered. All cases are submitted to the same simply-supported condition and the bi-sinusoidal pressure given by equation (9). The first case, case study I, is a one-meter four-ply square symmetric thick laminated plate with a stacking sequence of [0 /90 ]S. Its material properties can be found in the Ply row in Table 1.

pðx; yÞ ¼ p0 sinðpx=aÞsinðpx=bÞ

(9)

In equation (9), p0 is the magnitude of the pressure and the side lengths of the plate are a and b. For the second case, case study II, the sandwich structure is a one-meter wide and 0.1m thick square plate with stacking sequence of [0 /90 /core/90 /0 ]. The thickness of sandwich plate H is fractioned between skins and core. Since there are two plies in each skins, each skin is therefore, 0.05H thick and the core fills the remaining 0.8H of the thickness. The orthotropic material properties of the skins and the core are in Table 1. The last case study, case study III, comprises a rectangular plate (a/b¼900mm/300mm¼3) with three isotropic plies. For this sandwich structure, the thickness (H¼75mm) is fractioned so that the bottom, the core and the top ply are 0.1H, 0.7H and 0.2H thick, respectively. The material properties of this case can be seen in Table 2. In Table 2, the Young's modulus of the core may assume two

32

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

Fig. 2. Scaled stress field sx (a/h¼4).

values whose magnitudes are 104 times apart. They correspond to the two face-to-core stiffness ratios (FCSR¼10 and FCSR¼105), investigated via CGF. These ratios consider the moduli of the core and that of the bottom ply.

Fig. 3. Scaled stress field syy (a/h¼4).

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

Fig. 4. Scaled stress field szz (a/h¼4).

Before analyzing the stress results in next section, for CGF, the nodal stresses are calculated using the following steps: (i) The three displacement fields, the two rotations and the twist from CGF were built for each node.

33

Fig. 5. Scaled stress field tyz (a/h¼4).

(ii) These fields were evaluated throughout the thickness of the plate. (iii) The strains were then calculated at the centroid (x¼h¼0) of the elements using compatibility relations with the precalculated displacement fields in (ii)

34

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

Fig. 6. Scaled stress field txz (a/h¼4).

(iv) For each node, the nodal strain values were obtained from the average values from the surrounding elements. (v) Finally, the nodal stresses were evaluated directly from constitutive relations applied on the calculated strain fields in (iv).

Fig. 7. Scaled stress field txy (a/h¼4).

In addition, all finite element models were fully integrated and the 2D models (CGF and S4) have a 1010 mesh resolution whereas the 3D models (C3D8 and C3D20) were built with 10108 elements.

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

35

Table 6 Dimensionless normal stresses. a/h

4 10 20 100

sxx

syy

ESL

LW

Analytical [24]

C3D20

ESL

LW

Analytical [24]

C3D20

0.6929 0.5482 0.5284 0.5195

0.6926 0.5370 0.5272 0.5195

0.720 0.5590 0.5430 0.5390

0.7054 n/a n/a 0.5389

0.6260 0.3827 0.2964 0.2610

0.6610 0.4082 0.3021 0.2612

0.6630 0.4030 0.3090 0.2710

0.6795 n/a n/a 0.2680

Table 7 Dimensionless transverse normal stress and in-plane shear stress. a/h

4 10 20 100

szz

txy

ESL

LW

C3D20

ESL

LW

Analytical [24]

C3D20

0.0270 0.0124 0.0766 2.3560

0.8974 0.8504 0.7762 1.4949

0.9721 n/a n/a 1.572

0.0194 0.0163 0.0177 0.0203

0.0565 0.0172 0.0169 0.0202

0.0467 0.0275 0.0230 0.0214

0.0362 n/a n/a 0.0197

Table 8 Dimensionless transverse shear stresses. a/h

4 10 20 100

tyz

txz

ESL

LW

Analytical [24]

C3D20

ESL

LW

Analytical [24]

C3D20

0.2504 0.1667 0.1431 0.3934

0.2565 0.2235 0.1789 0.4175

0.2920 0.1960 0.1560 0.1390

0.1134 n/a n/a 0.4512

0.2177 0.2835 0.3057 0.4173

0.1946 0.3192 0.3400 0.4410

0.2190 0.3010 0.3280 0.3390

0.0747 n/a n/a 0.5785

3.1. Case study I Case study I provides two sets of results comparing different parameters. The first set investigates the stress results obtained for a thick laminate (a/h¼4) calculated by different plate theories. By choosing a thick laminate, it becomes easier to demonstrate the difference between the use of different plate theories. The second set of results compares the behavior of the stress profiles throughthe-thickness of the laminated structure for different length-tothickness ratios (a/h¼4,10,20,100). In both cases, the computational predictions are compared to the results obtained from the commercial software Abaqus. Two dimensional problems meshed with shell elements (S4) and three dimensional models meshed with brick elements (C3D8 and C3D20) as well as an analytical approach are used to evaluate the accuracy of CGF. The values of stresses at particular coordinate points are used for the comparative analyses. Equation (10) indicates, respectively, the multiplying scaling factors (sf) for the in-plane normal stresses (sxx,yy), transverse normal stress (szz) and the shear stresses (sxz,yz,xy).

sfxx;yy ¼

1 ðh=aÞ2 P0

sfzz ¼

1 P0

sfxz;yz;xy ¼

1 ðh=aÞ P0

(10)

The set of points assumes a coordinate system located at the mid-plane of the lower corner of the plate. Thus, the center of the plate is at x¼y¼a/2;z¼0. It is important to observe that the stress values at the top and bottom faces of the laminate have the same magnitude for equivalent single layer models, because the laminate is symmetric. In LW cases, even for symmetric laminates, this cannot be stated because there is one displacement field for each ply, unlike the ESL results, which come from one single displacement field shared among the plies. For brevity, only the bottom values are compared in Tables 4 and 5. However this is usually implicit in the references. In the original reference [24], one can check the values at the bottom and

Table 9 Stress comparison between CGF and the analytical reference for case study II. Theory ESL LW

111 333 111 333

Abaqus 2D Abaqus 3D Analytical [25] a

Comparison with Abaqus 3D.

sxx

syy

szza

txy

tyz

txz

Value

9.37 3.21 2.29 2.29 1.0799 1.1542 1.1417

5.22 3.32 4.75 5.70 0.0539 0.0615 0.0632

169.20 29.19 11.38 17.01 n/a 0.9403 n/a

55.84 49.89 51.59 50.53 0.0141 0.0045 0.0471

54.37 18.42 0.30 6.97 0.1659 0.1545 0.1694

52.42 16.48 0.22 6.85 0.1790 0.1649 0.1839

%

Scaled

36

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

at the top of the plate. For syy, the considered value is the one taken at the interface of the 1st to the 2nd ply and at the second ply. Based on Tables 4 and 5, it is possible to observe that the layerwise formulations are sensitive to every change in number of expansion terms. Nonetheless, it is worth mentioning that the transverse normal stress exhibits an opposite response for cubic polynomials. From LD22x to LD33x, the results worsen. This probably comes from the equilibrium of the transverse functions F(z), which are combinations of the Legendre polynomials in the case of layerwise formulation [23]. However, this should be interpreted as a mathematical result of the use of a unified formulation. If the polynomial order is increased (e.g. to 5th), the results are expected to improve, not to worsen. For the ESL solution, this is not verified. For the in-plane normal stresses, the results are improved in two cases. Firstly, when the formulation upgrades from EDxx1 to EDxx2, this yields a minor enhancement. Secondly, when the formulations changes from ED22x to ED33x. In this case, the change is substantially high for both sxx and syy. For the transverse normal stress szz, the same pattern is observed. The quadratic and cubic formulations seem to disturb the transverse normal stress and increase the deviations from the analytical value. Figs. 2e4 endorse these findings with detailed normal stress patterns through-the-thickness of the laminate. Figs. 5e7 give the same information for the shear components.

Fig. 9. Scaled stress field syy (a/h¼10).

Fig. 8. Scaled stress field sx (a/h¼10).

These graphics demonstrate how important it is for the designer to understand the plate theory used in these analyses. For this case study II, the geometry and the loading on this particular plate are symmetric. However, the material properties are not and this yields an orthotropic deformation in the global coordinate system xyz. Because of this, it is easy to verify from Figs. 2 and 3 that the use of a linear transverse function affects the sxx stress pattern more than the syy one, including the results obtained with Abaqus shown in Fig. 2c. Based on the transversal stresses in Fig. 5 through 7, it is shown that non-linear and layerwise theories for thick laminated plates are required. Moreover, when applying CGF, the discontinuities of the transverse stresses at the interfaces decreases and, for sufficiently non-linear theories, these discontinuities may actually vanish, numerically at least. In Fig. 5c, it is important to note the difference between the C3D8 and C3D20 results. Both models have 2 bricks in the thickness direction to represent each ply and even so, high discrepancies are observed. If a convergence test were to be performed, several additional elements would be required through-the-thickness of the plate. A proof of that is the non fulfillment of the top and bottom ply free transverse shear stress condition of the presented examples. Regarding the results calculated with the S4 elements, the nodal transverse shear stresses txz and tyz were obtained directly from Abaqus variables THSR13 and THSR23, respectively. These variables are available only for thick-shell elements and are calculated at different integration points through the section of the

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

37

Fig. 10. Scaled stress field szz (a/h¼10).

Fig. 11. Scaled stress field tyz (a/h¼10).

laminate. Simpson's integration rule was chosen because Gauss' in Abaqus does not output any values at the interfaces of the plies. Also, the respective curves in Figs. 5c and 6c, unlike the solid model responses, fulfill the free shear stress at the top and bottom plies faces. Moreover, the magnitude and shape of the curves are limited by the transverse shear strain energy and the ESL approach, respectively. In this work, CGF does not make use of shear correction factors to clearly show the result of using non-linear theories. However, for laminated shells in Abaqus, a shear correction factor for each ply is estimated. Based on a pre-calculated transverse shear stress profile for the shell, the respective strain energy is equated to that estimated from the forces and strains in each section. Concerning the szz stress, due to the directly calculation of stresses, qualitative deviations can be observed in the ESL results (Fig. 4a). It is once again evident that thick laminates cannot be handled with thin plate theories. Since CGF is a unified formulation, it englobes both thin and thick plate cases. In Table 6 through 8, the behavior of CGF for different plate thickness is investigated by using two plate theories (ED332 and LD332). From these results, it is easy to observe how the deviations between different theories diminishes as the plate gets thinner. Nonetheless, considerably high deviations remain in the transverse normal stress results szz and there are, in Table 7, two facts concerning this variable. First, there is no analytical reference for such stress, because it was disregarded

in reference [24]. This is usually the case for thin plate theories because the εzz strain is normally neglected. However, the respective stress may not be neglected. Hence, its calculation should not be performed directly. This can be related with the second fact, shown in the curves of Fig. 4a. The stress recovery process used with CGF, does not capture the unit value at the center of the top ply of the laminate. Actually, for the ESL case, the results deviate both qualitative and quantitatively. The LW case also struggles to recover this unit value. Abaqus' results can be explained due to poor mesh resolution, but CGF cannot. CGF computes the results directly and this result should be ignored because the transverse strain εzz is just too small to counterweight the other in-plane strains (equation (11)). The authors did keep these calculations for consistency of the results obtained via CGF.

k

k

k

skzz ¼ C 31 εxx þ C 32 εyy þ C 33 εzz

(11)

In Table 8, it shown that CGF's results agree well with both analytical and Abaqus ones. However, CGF agrees more with Abaqus, because it is a finite element solution method with C-1 continuity of the transverse displacement, which relaxes the equilibrium equations and increases the physical consistency of the plate problem.

38

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

Fig. 12. Scaled stress field txz (a/h¼10).

Fig. 13. Scaled stress field txy (a/h¼10).

3.2. Case study II For the coordinate points in Table 3, the stresses by CGF are compared to the analytical solution [25] and the results obtained by Abaqus. The 3D model built within Abaqus has five elements to represent the core and one element for each ply within the skins. For CGF, two theories were chosen. The linear one (LD/ED111), which is better represented by Mindlin's theory, despite the nonconstant transverse displacement field, and a cubic one (LD/ ED333), which is more physically consistent resembling Reddy's theory. The comparison of the magnitudes extracted from all methods are shown in Table 9. Inspecting Table 9, it is verified that the use of non-linear and/or LW theories increases the accuracy of the predictions for most of the stress values. Actually, for LW theories, the switch from linear to cubic has an overall negative impact on the results. The overall accuracy decreased 27% from LD111 to LD333 and the main contributions to such decrease come from the transverse stresses. It seems that for this particular case, the additional terms in the transverse functions could not predict the real physical deformation. Thus, the use of non-linear theories does not always represent a better option to a particular application. The particular high deviations (>50%) for the transverse stresses values should be interpreted with care. The reference solution is the analytical one. However, if Abaqus's 3D solution is assumed, this variation can be even higher. Particularly, the shear in-plane stress txy is the one

Fig. 14. Influence of mesh refinement on the txz stress profile of case study II.

with high deviations even for the non-linear layerwise approach in both quantitative and qualitative aspects (Fig. 13b). This again is explained by the boundary condition effects in the 3D models. In this same Fig. 13b, there are two results for the txy variable extracted from Abaqus's 3D model. The C3D20 through-the-

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

39

Table 10 Dimensionless transverse displacement and shear stress (zskin¼3h/10) of case study III.

txz

w FCSR Theory GUF [11] CGF Abaqus 2D Abaqus 3D Analytical [26]

10 ED225 2.0330 2.1154 3.1032 3.8885 3.0112

LD225 3.0098 3.8819

105 ED225 4.098E-4 2.664E-4 1.4457 0.1352 0.0132

Fig. 15. Scaled stress field sx (a/h¼12;b/h¼4).

thickness stress profile is the one according to Table 3. The other result, the C3D20diag, is the stress profile obtained at the opposite node of this element's diagonal (towards the center of the plate). Comparing these two curves, both the qualitative and quantitative results from CGF are once more endorsed by the 3D model because the CGF stress profiles are in between both results (C3D20 and C3D20diag). Despite these deviations, the overall quantitative results from Table 9 are good considering the sandwich structure. In addition, the stress profiles through-the-thickness of the structure agree considerably well with the results from Abaqus. The authors believe that higher order theories (e.g. LD444þ) should yield better results because these theories have more degrees of freedom to allow for a better fitting (or transition) of the

LD225 0.0132 0.1337

10 ED225 0.3250 0.641 0.3330 0.0055 0.3217

LD225 0.3214 0.256

105 ED225 0.3317 0.918 0.3203 0.0035 5.408

LD225 5.400E-4 0.056

Fig. 16. Scaled stress field syy (a/h¼12;b/h¼4).

boundary conditions effects into the bulk solution. It is important to highlight the fact that in Figs. 8e13, the curves were built from the nodal values of the C3D20 element, i.e. nine values were used. Because of this, these curves do not look smooth. On the other hand, the S4 results have five values per ply, including the interface heights, what smoothens the respective curves. In Fig. 10 through 13, neither the shear nor the transverse stresses agree very well due to different reasons. The transverse normal stress is influenced by the ESL assumptions and the direct post-processing. Nonetheless, the non-linear LW results are close to Abaqus'. It is important to note that the authors did not choose to build a refined 3D model to provide a more direct comparison of the FE approaches regarding the number of elements used. In addition, it highlights the need for a mesh convergence study for complex structures and

40

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

Fig. 17. Scaled stress field szz (a/h¼12;b/h¼4).

shows that CGF results can also be used for quick theory checks and assessments. To endorse this statement, Fig. 14 gives the results of txz with 9 and 18 elements in the thickness direction of the sandwich structure. This graphic shows that by doubling the number of elements, the free transverse shear stress condition at both bottom and top faces of the sandwich structure, is almost fulfilled. A further mesh refinement in either in-plane or thickness direction should correct this curve. 3.3. Case study III In this case, the analytical solution for a single point is selected to compare the transverse shear stress txz obtained via CGF, by Demasi [11] and via Abaqus. The selected point is not the one from Table 3 but rather: x ¼ 0; y ¼ b/2; z ¼ 3h/10. This point corresponds to the bottom face of the upper skin. The investigated sandwich structure is made of three different isotropic plies. The aspect ratio of b/a¼3 and the length-to-thickness ratio of b/h¼4 makes this structure an anisotropic plate, in geometric terms. Also, two values of the FCSR (Ebot_skin/Ecore) are investigated. This ratio is the quotient of the elastic modulus of the material of the bottom skin and the respective moduli of the core. The results are shown in Table 10. Since Demasi's solution [11] for this sandwich plate problem covers both ED225 and LD225 theories, such theories were also derived

Fig. 18. Scaled stress field tyz (a/h¼12;b/h¼4).

from CGF for direct comparison. All meshes have 10  10 elements in plane and the 3D model has 8 elements in the thickness direction (five in the core, one representing the lower skin and two for the upper skin). In Ref. [11], Demasi solved this sandwich problem via a Naviertype solution and reported the values of txz, shown in Table 10, calculated at the bottom face of the upper skin. From this table, there is no clear agreement between any pair of txz values calculated from different approaches but the pair GUF/[26]. Nevertheless, these two reference results can be interpreted as equivalent, because both methods (GUF and [26]) are based on double Fourier solution methods. On the other hand, taking the C3D20 result as the best reference, different conclusions can be drawn from Fig. 15 through 20. For the in-plane normal stresses, sxx and syy, the LW results agree quite well with the ones from Abaqus. Nonetheless, both the ESL results from CGF and S4 deviate for both FCSR cases and the FCSR¼100000 gives the worst set of results, because the peaks at the bottom and top skins interfaces with the core are not properly represented. Concerning the transverse normal stress calculated by the ESL methods, both qualitative and quantitatively results do not agree with Abaqus's 3D model. The LW theory provides the best alternative, even though these results also deviate from the expected null and unit values at the bottom and top heights of the sandwich plate, respectively. Regarding the transverse shear stresses, both Figs. 18 and 19 show that the LD225 theory best approaches the C3D20 results. The same applies to the

M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

Fig. 19. Scaled stress field txz (a/h¼12;b/h¼4).

sxy shear stress (20). One more time, the results from Abaqus demand a higher mesh resolution in all directions, especially the thickness one, to yield consistent through-the-thickness stress profiles of the transverse stresses. The present work does not seek to provide a mesh independent model of this test case via Abaqus, but rather to compare different FE procedures and solution methods for laminated plates and shells. Moreover, it does not seek a best theory derived from CGF for this sandwich plate problem. This can be accomplished with methods such as the Best Plate Diagram proposed in Ref. [27]. The choice of theories (ED225/LD225) had the purpose of offering a direct comparison with Demasi's works [26,11]. Finally, the influence of the additional C-1 degrees of freedom (DOFs) of CGF the transverse variable w (w,x,w,y and w,xy) were not fully explored in this paper. However, for the three structures herein studied, the simply supported condition assumed these DOFs as unknowns. In Table 10, the normalized transverse displacement variable w, (w ¼ wð100E2 h3 Þ=ðP0 a4 Þ), is given to show that, besides the stress values, the LD225 displacement results via CGF provides a better agreement with Abaqus result (99.8%) than the respective GUF result (79.8%). A detailed discussion on the displacement fields obtained through CGF the authors can be found in Refs. [22,23].

41

Fig. 20. Scaled stress field txy (a/h¼12;b/h¼4).

4. Conclusion This work showed, for three laminated plate structures, a complete set of the through-the-thickness stress profiles obtained using CGF. Based on these stress profiles, the following conclusions arise:  CGF and the current stress recovery method represent an accurate solution for sxx and syy for both thin and thick laminated and sandwich plates.  For thin to moderately thick plates, CGF can yield good transverse shear stress value for both ESL and LW methods.  The szz stress calculated by layerwise theories from CGF exhibited a good through-the-thickness profile, and for some cases, the null/unit scaled values at the bottom and top of the laminate were achieved.  Before switching to layerwise theories, the use of non-linear ESL theories is strongly advised.  The use of higher order theories, available in UFs, such as CGF, does not always give better results than linear or lower order theories.  Thick sandwich plates are complex structures which, need a 3D solution method for exact calculation of transverse stresses. However, LW plate/shell formulations, can be optimized, with a

42









M.F. Caliri Jr. et al. / Composites Part B 107 (2016) 29e42

UF tool such as CGF, to provide a faster alternative to full 3D models. The stress profiles shown in this work highlight the capability of the finite element method of including the boundary condition effects in the results and the need for mesh independence studies when designing complex structural regardless of the FE formulation. The stress profiles shown in this work were directly computed from constitutive equations. Further post-processing may improve these results. The use of a FE-based Unified Formulation is preferred over analytical solution methods because of the advantages that the Finite Element Method has for solving engineering problems. Further works on CGF include a refined post-processing of the stresses, development of a shell element, and inclusion of a damping kernel matrix in CGF.

[9]

[10] [11]

[12] [13] [14]

[15] [16] [17]

Acknowledgements The authors acknowledge the financial support of the National Council for Scientific and Technological Development (CNPq process number: 401170/2014-4 and 150492/2015-4). Volnei Tita acknowledges the financial support of Air Force Office of Scientific Research e AFOSR (Grant FA9550-16-1-0222).

[18]

[19]

[20]

References [21] [1] Hansen MOL, Sorensen JN, Voutsinas S, Sorensen N, Madsen HAa. State of the art in wind turbine aerodynamics and aeroelasticity. Prog Aerosp Sci 2006;42: 285e330. [2] Thomsen OT. Sandwich materials for wind turbine blades e present and future. J Sanwich Struct Mater 2009;11(1):7e26. [3] Ivanell SSA. Numerical computations of wind turbine wakes. Stokholm: Royal Institute of Technology; 2005. Ph.D. thesis. [4] Gibson LJ, Ashby M. Cellular solids: structures & properties. England: Pergamon Press-Headington Hill Hall; 1988. iSBN: 0-08-035910-8. [5] Carrera E. Theories and finite elements for multilayered anisotropic, composite plates and shells. Archives Comput Methods Eng 2002;9(2):87e140. [6] Carrera E. Theories and finite elements for multilayered anisotropic, composite plates and shells:a unified compact formulation with numerical assessment and benchmarking. Archives Comput Methods Eng 2003;10(3): 215e96. [7] Carrera E, Demasi L. Classical and advanced multilayered plate elements based upon pvd and rmvt. part1: derivation of finite element matrices. Int J Numer Methods Eng 2002;55:191e231. [8] Carrera E, Demasi L. Classical and advanced multilayered plate elements based

[22]

[23]

[24] [25]

[26] [27]

upon pvd and rmvt. part 2: numerical implementations. Int J Numer Methods Eng 2002;55:253e91. Demasi L. Treatment of stress variables in advanced multilayered plate elements based upon reissner's mixed variational theorem. Comput Struct 2006;84:1215e21. Demasi L. Hierarchy plate theories for thick and thin composite plates: the generalized unified formulation. Compos Struct 2008;84:256e70. Demasi L. An invariant model for any composite plate theory and fem applications: the generalized unified formulation. In: 50th AIAA/ASME/ASCE/ AHS/ASC structures, structural dynamics and materials conference; 2009. Demasi L. Mixed plate theories based on the generalized unified formulation. part i: governing equations. Compos Struct 2009;87:1e11. Demasi L. Mixed plate theories based on the generalized unified formulation. part ii: layerwise theories. Compos Struct 2009;87:12e22. Demasi L. Mixed plate theories based on the generalized unified formulation. part iii: advanced mixed high order shear deformation theories. Compos Struct 2009;87:183e94. Demasi L. Mixed plate theories based on the generalized unified formulation. part iv: zig-zag theories. Compos Struct 2009;87:195e205. Demasi L. Mixed plate theories based on the generalized unified formulation. part v: Results. Compos Struct 2009;88:1e16. Demasi L. Partially zig-zag advanced higher order shear deformation theories based on the generalized unified formulation. Compos Struct 2012;94: 363e75. Ferreira AJM, Carrera M, Cinefra M, Roque CMC, Polit O. Analysis of laminated shells by a sinusoidal shear deformation theory and radial basis functions collocation, accounting for through-the-thickness deformations. Compos Part B Eng 2011;42:1276e84. Ferreira AJM, Araujo AL, Neves AMA, Rodrigues JD, Carrera E, Cinefra M, et al. A finite element model using a unified formulation for the analysis of viscoelastic sandwich laminates. Compos Part B Eng 2012;45:1258e64. Mantari JL, Oktem AS, Soares CG. A new higher order shear deformation theory for sandwich and composite laminated plates. Compos Part B Eng 2012;43:1489e99. Caliri MFJ, Ferreira AJM, Tita V. A review on plate and shell theories for laminated and sandwich structures highlighting the finite element method. Compos Struct 2016 (in press). Caliri MFJ, Ferreira AJM, Tita V. New generalized unified solution method for thin laminated plates. Am Inst Aeronautics Astronautics J - AIAA J 2016;54(8): 2556e60. Caliri MFJ, Ferreira AJM, Tita V. A new finite element for thick laminates and sandwich structures using a generalized and unified plate theory. Int J Numer Method Eng 2016. http://dx.doi.org/10.1016/j.compstruct.2016.02.036 (in press). Pagano NJ, Hatfield SJ. Elastic behavior of multilayered bidirectional composites. Am Inst Aeronautics Astronautics J - AIAA J 1972;10(7):931e3. Pandit MK, Sheikh AH, Singh BN. Analysis of laminated sandwich plates based on an improved higher order zigzag theory. J Sandw Struct Mater 2010;12: 307e26. Demasi L. 2d, quasi 3d and 3d exact solutions for bending of thick and thin sandwich plates. J Sandw Struct Mater 2008;10(4):271e310. Carrera E, Miglioretti F, Petrolo M. Guidelines and recommendations on the use of higher order finite elements for bending analysis of plates. Int J Comput Methods Eng Sci Mech 2011;12:303e24.