Thermo-mechanics of undamaged and damaged multilayered composite plates: assessment of the FEM sub-laminates approach

Composite Structures 81 (2007) 137–155 www.elsevier.com/locate/compstruct

Thermo-mechanics of undamaged and damaged multilayered composite plates: assessment of the FEM sub-laminates approach Marco Gherlone *, Marco Di Sciuva

1

Aerospace Engineering Department, Politecnico di Torino, Corso Duca degli Abruzzi, 24, 10129 Torino, Italy Available online 18 September 2006

Abstract Based on a previously developed higher-order plate theory and on the relative plate mixed finite element, a FEM approach based on the sub-laminates concept is proposed. A large number of numerical examples are presented and discussed in order to assess deeply the capabilities of the proposed approach to solve some non-classical problems. At first, it is shown the possibility to model plates subjected to tangential loads distributed on the top and bottom faces of the laminate. Moreover, the proposed sub-laminate strategy is particularly suitable for describing the effects of damaged interfaces, also in case of severe levels of delamination; particular attention is dedicated to explain the damage modelling strategy and its numerical implications. Finally, it is possible to study the effects of prescribed temperature distributions on the static response of multilayered composite plates.  2006 Elsevier Ltd. All rights reserved. Keywords: Multilayered plates; Tangential loads; Damaged plates; Thermal loads; Sub-laminates; Plate mixed finite element

1. Introduction The increasing use of composite and sandwich laminates as primary structural components stresses the need for models suitable to describe some typical behaviours of multilayered structures: transverse anisotropy and deformability, damage and environmental effects sensibility among the others [1–4]. The approaches usually adopted to face these modelling challenges are layerwise theories based on both displacement and stress unknowns and layerwise theories based on displacements only; the latter are a good compromise between accuracy and computational complexity when the number of unknowns does not depend on the number of layers as happens for zig-zag theories pioneered by Di Sciuva [5–8]. These models account for transverse shear deformability with continuity of transverse shear stresses (vanishing on top and bottom faces *

Corresponding author. Tel.: +39 0115646817; fax: +39 0115646899. E-mail addresses: [email protected] (M. Gherlone), marco. [email protected] (M. Di Sciuva). 1 Tel.: +39 0115646808; fax: +39 0115646899. 0263-8223/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2006.08.005

of the laminate) and are able to capture the so-called zig-zag behaviour, typical of thick and/or heterogeneous laminates. A recent development of the original zig-zag approach is the Hermitian Zig-Zag (HZZ) plate theory [9,10] that improves the modelling capabilities in particular in terms of transverse normal deformability. Moreover, the degrees of freedom of the model are the displacements and the transverse shear stresses of the top and bottom laminate faces; this allows coupling the HZZ plate theory with the sub-laminates approach. For a broad discussion on the HZZ plate theory and on the sub-laminates philosophy, see [10,11]; a summary of the main aspects of this modelling strategy is contained in Section 2 of this article. The main scope of this article is to assess the capabilities of the FEM methodology based on the Hermitian Zig-Zag plate finite element and on the sub-laminates approach also for non classical applications: • Tangential loads on the top and bottom faces: these loads can derive from friction due to contact or can simulate the effect of piezo-actuators attached to the laminate; in any case it is important that the transverse shear stresses

138

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

Nomenclature a, b, h in-plane dimensions and thickness of the plate L common value of the sides for a square plate N number of orthotropic layers n number of sub-laminates xi (i = 1, 2, 3) Cartesian orthogonal coordinates xa (a = 1, 2) in-plane orthogonal coordinates x3 thickness coordinate Vi displacement field eij (ca3) strain tensor (transverse shear components) rij (sa3) stress tensor (transverse shear components) DVa intensity of the in-plane displacements jump across an interface Tab sliding constants d sliding function giving the in-plane displacements jump as a function of the transverse shear stresses Ei (i = 1, 2, 3) Young moduli of the material Gij Shear moduli of the material mij Poisson coefficients of the material ai (i = 1, 2, 3) thermal expansion coefficients of the material

temperature variation with respect to the reference value u0i ; uhi ði ¼ 1; . . . ; 8Þ degrees of freedom of the bottom and top faces for in-plane displacement V1 v0i ; vhi ði ¼ 1; . . . ; 8Þ degrees of freedom of the bottom and top faces for in-plane displacement V2 w0i ; w0;xi ; w0;yi ; whi ; wh;xi ; wh;yi ði ¼ 1; . . . ; 4Þ degrees of freedom of the bottom and top faces for transverse displacement V3 and bending rotations s0xzi , shxzi ði ¼ 1; . . . ; 8Þ degrees of freedom of the bottom and top faces for transverse shear stresses s13 s0yzi , shyzi ði ¼ 1; . . . ; 8Þ degrees of freedom of the bottom and top faces for transverse shear stresses s23 [M], [K], {F} mass matrix, stiffness matrix and force vector of the plate finite element [V], {A} terms of the matrix equation governing the behaviour of the interfacial nodes with finite debonding {q} degrees of freedom of the plate finite element {k} Lagrange multipliers vector H

recover the value of the external loads on top and bottom faces. • Multiple delaminations: the original zig-zag models by Di Sciuva and the HZZ plate theory may be used to analyse plates with sliding type damages but it has been shown [11] that, for severe level of delamination, results may be inaccurate; the adoption of the sub-laminates approach allows modelling damages in a more correct way without any limitation on the damage pattern (see Section 3 for the whole discussion on damage modelling). • Thermal effects: the effect of temperature variations on a multilayered plate, in particular when exhibiting high transverse anisotropy, may lead to delaminations [4]; for this reason undamaged and damaged plates subjected to uniform or linear temperature field will be studied.

1. the in-plane displacements are through-the-thickness piecewise cubic (the use of shape functions that are a generalisation of the Hermite cubic polynomials is the origin of the ‘‘Hermitian’’ adjective) while the transverse displacement is linear; it is possible to model interlayer slip by means of jumps of axial displacements at the damaged interfaces; 2. transverse shear stresses are continuous along the thickness; 3. the model satisfies the traction equilibrium condition on both external surfaces, also in case of tangential loads; 4. the transverse normal deformability is taken into account adopting the complete constitutive equations and assuming the transverse normal stress to be constant along the laminate thickness; 5. the degrees of freedom of the model are the displacements and transverse shear stresses on the top and bottom faces.

Moreover, the applications will be performed for the case of plates with a clamped edge and the other edges free, being a more complicated case that the traditional simply supported panel. Due to the particular kind of applied mechanical and thermal loads, boundary conditions and damages, it has not been possible to find an analytical exact solution to be used as a reference result; for this purpose a 3D detailed FEM model has been used.

For the details about the formulation of the Hermitian Zig-Zag plate theory, see [9,10]. Based on the Hermitian Zig-Zag plate theory, a plate finite element has been formulated [10] with 8 nodes belonging to the top face and other 8 nodes belonging to the bottom one, according to the same scheme adopted by the traditional QUAD8 membrane element, but repeated for each face (Fig. 1). The degrees of freedom are:

2. Plate mixed finite element We briefly summarize the main properties of the Hermitian Zig-Zag plate theory on which the plate finite element formulation is based [9,10]:

• the in-plane displacements u and v and the transverse shear stresses sxz and syz; these degrees of freedom are defined in each of the 16 nodes of the element and the Serendipity parabolic shape functions are used for the interpolation;

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

139

3. Damage modelling One of the typical damages present in multilayered plates and due to manufacturing and/or operating conditions is delamination, i.e. a damage occurring at some interfaces between layers and characterised by a discontinuity of the in-plane displacements [3,12,13]; if one supposes that the delamination may not open, there is no discontinuity of the transverse displacement. Considering that delaminations are due to a weakening of the adhesive film between layers, some authors [14,15] adopt the socalled ‘‘compliant layer’’ method; it is supposed that a very thin layer is present at a damaged interface and very low mechanical properties are attributed to it. The Hermitian Zig-Zag plate theory allows modelling delaminations too [9–11] but without considering compliant layers. A linear slip law is adopted [9,10] that links the in-plane displacements jump (DVa) to the common value of transverse shear stresses at the damaged interface (sa3) DV a ðx1 ; x2 Þ ¼ T ab ðx1 ; x2 Þsb3 ðx1 ; x2 Þ

Fig. 1. Plate mixed finite element topology.

ð1Þ

It has been shown that this way of describing the effect of delaminations is correct for low levels of damage (moderate values of Tab) [11]; when we are studying severe damages the model exhibits a stiffening effect, called sliplocking phenomenon [11], leading to erroneous results.

• the transverse displacement w and its derivatives w,x and w,y; these degrees of freedom are defined only at the 8 corner nodes of the element and the Hermite cubic polynomials are adopted for their interpolation (these polynomials have no relations with the though-the-thickness shape functions of the Hermitian Zig-Zag plate theory). The mass matrix, the stiffness matrix and the mechanical force vector are calculated in a straightforward way following the classical FEM approach [10,11]. When thermal loads are prescribed, the equivalent thermal nodal forces are calculated, assuming that temperature, for each finite element, is through-the-thickness linear and distributed in x1 and x2 as well as the in-plane displacements [10]. The fact that nodes and degrees of freedom are located and defined at the laminate top and bottom faces, allows the introduction of the sub-laminates concept [10,11]. A sub-laminate is a group of adjacent layers or part of a single layer. In order to increase the accuracy of the FEM model of a laminate, it is possible to divide its thickness in some sub-laminates each of them corresponding to a finite element (the usual mesh in the (x1, x2) plane is obviously still present).

Fig. 2. Interfacial degrees of freedom.

Table 1 Material mechanical properties Material

E1 (MPa)

E2, E3 (MPa)

G12, G13 (MPa)

G23 (MPa)

m12, m13

m23

a1, a2 (C1)

a3 (C1)

a b

111,000 200,000

7900 8000

3300 5000

2000 2200

0.33 0.25

0.49 0.35

– 2e6

– 50e6

140

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

Table 2 Stacking sequences of the considered three-layered laminates (layers are considered along the positive x3-axis direction) Laminate

Thicknesses (mm)

Orientations ()

Materials

1 2

(0.25/0.25/0.25/0.25) (0.25/0.25/0.25)

(0/90/0/90) (0/90/0)

(a/a/a/a) (b/b/b)

The sub-laminates approach gives an alternative and effective way to model severe damage too; let us see how it is possible to model delaminations. Two elements – sub-laminates (a and b) share the 8 nodes of their common face (Fig. 2); we focus our attention on one of these nodes, for example node 2 of the lower face

Fig. 3. Mechanical loads considered for the numerical applications.

Fig. 4. Delamination geometries considered for the numerical applications.

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

of element b and node 2 of the upper face of element a (Fig. 1). There are three possible cases. • Undamaged interface: the in-plane displacements are continuous through the interface and the same is obviously valid for the transverse shear stresses: ðu02 Þb ¼ ðuh2 Þa

ð2Þ

ðv02 Þb ¼ ðvh2 Þa ðs0xz2 Þb

¼

ðs0yz2 Þb

¼ ðshyz2 Þa

141

• Damaged interface, partial debonding: transverse shear stresses are still continuous but not zero, thus giving a weak bonding between the two layers; in-plane displacements are obviously discontinuous: ( 0 )   ðu2 Þb ðuh2 Þa ð4Þ  ¼ dððs0xz2 Þb ; ðs0yz2 Þb Þ ðv02 Þb ðvh2 Þa ðs0xz2 Þb ¼ ðshxz2 Þa

ð5Þ

ðs0yz2 Þb ¼ ðshyz2 Þa

ðshxz2 Þa

ð3Þ

Conditions (2) and (3) may be imposed through the classical FEM assembling procedure.

While condition (5) is again satisfied by assembly, the relation between the in-plane displacements jump and the transverse shear stresses may be of several types

1

0.9

0.8

0.7

x3

0.6

0.5

0.4

0.3

0.2

HZZ4 NASTRAN

Fig. 5. Mesh of the MSC/NASTRAN models. 0.1

HZZ1 HZZ2

0 –1.22

–1.21

–1.2

–1.19

–1.18

–1.17

–1.16

–1.15

–1.14

Table 3 Number of elements and of degrees of freedom for the FEM models FEM model

Elements

Dofs

NASTRAN (HEXA8) HZZ1 HZZ2 HZZ3 HZZ4

19,200 (40 ·40 · 12) 64 (8 · 8 · 1) 128 (8 · 8 · 2) 192 (8 · 8 · 3) 256 (8 · 8 · 4)

65,559 2286 3429 4572 5715

–1.13 x 10–3

V3(L/2,L/2,x3)

Fig. 6. V3(L/2, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

1

0.9

0.8

0.7

0.6

x3

Table 4 List of numerical applications and of the relative figures Laminate (Table 2)

Mechanical load

Thermal load

Damages (Fig. 4)

Figures

1 2 2 1 1 1 1 1 2 2

TG – – TR TG TR TG TR – –

– UNI LIN – – – – – UNI UNI

– – – A A A+B A+B A C C+D

6–13 14–16 17, 18 19–25 26–33 34–39 40–47 48–50 51 52

0.5

0.4

0.3 HZZ4

0.2

HZZ2 0.1

HZZ1 NASTRAN

0 –9.1

–9

–8.9

–8.8

V3(L,0,x3)

–8.7

–8.6

–8.5 –3

x 10

Fig. 7. V3(L, 0, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

142

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

(depending on the nature of function d). Thus, this way of modelling damages is more general than the original one contained in the Hermitian Zig-Zag formulation and limited to linear slip law only [9]. However, and for sake of simplicity, let us consider again a linear law of delamination ( 0 )  )   ( 0 ðu2 Þb ðsxz2 Þb T 11 T 12 ðuh2 Þa  ð6Þ ¼ T 21 T 22 ðv02 Þb ðvh2 Þa ðs0yz2 Þb • Damaged interface, total debonding: in this case transverse shear stresses are zero at the interface and the inplane displacements exhibit a jump with no relation with the shear stresses

ðu02 Þb 6¼ ðuh2 Þa

ð7Þ

ðv02 Þb 6¼ ðvh2 Þa ðs0xz2 Þb ¼ 0 ðshxz2 Þa ¼ 0

ð8Þ

ðs0yz2 Þb ¼ 0 ðshyz2 Þa ¼ 0

Eq. (7) says that the in-plane degrees of freedom of the two elements at common node 2 are independent from one another and like that must be treated during the solution of the equations of motion; the vanishing of trans-

1

1 HZZ1 0.9

0.9

HZZ2 HZZ4

0.8

0.8

NASTRAN

0.7

0.7

0.6

x3

x3

0.6

0.5

0.5

0.4

0.4

0.3

0.3

HZZ4

0.2

0.2

NASTRAN HZZ1

0.1

0.1

HZZ2

0 1.5

0

1.6

1.7

1.8

1.9

2

2.1

2.2

V3(L,L,x3)

–2

–1

0

1

2

3

4

γ13(L/2,L/2,x3)

x 10–3

Fig. 8. V3(L, L, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

x3

0.6

0.5

0.5

0.4

0.4

0.3

0.3

HZZ4

0.2

HZZ1

0.2

NASTRAN

HZZ2

HZZ1

0.1

HZZ4

0.1

HZZ2

NASTRAN 0 –4

5 x 10–4

Fig. 10. c13(L/2, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

1

x3

3

–3

–2

–1

V1(L,L,x3)

0

1

0 –5.5

2 x 10–4

Fig. 9. V1(L, L, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

–5

–4.5

–4

–3.5

–3

–2.5

γ23(L/2,L/2,x3)

–2

–1.5

–1

–0.5 x 10–4

Fig. 11. c23(L/2, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

verse shear stresses may be considered as a normal boundary condition. When analysing a damaged plate, the equations that must be added to the equations of motion of the structure are like Eq. (6) and these appear only when the damage does not correspond to a total debonding. In all the other cases, the mathematical equations describing damage reduce to assembling operations (Eqs. (2), (3), (5) and (7)) or boundary conditions (Eq. (8)). Let us discuss how to manage the case of not totally debonded interfaces. Once we have written the whole set of equations like Eq. (6), valid for all the nodes where a damage is present, we can collect them in a compact matrix form

143

½V fqg ¼ fAg

ð9Þ

where {q} is the vector of degrees of freedom of the whole laminate, [V] and {A} are a constant matrix and a constant vector (in particular, in our case, {A} = {0}). In the case of a dynamic analysis, we should consider the system of Eq. (9) and the classical equation of motion ½Mf€qg þ ½Kfqg ¼ fF g

ð10Þ

where [M] is the mass matrix, [K] the stiffness matrix and {F} the force vector (sum of the mechanical and thermal contributions [10]). It is possible to solve the problem via the classical approach of the Lagrange multipliers, thus obtaining the modified equation of motion

1 0.7

0.9

0.8

0.6

0.7 0.5 HZZ3

x3

x3

0.6

0.5

0.4

0.4

HZZ1 NASTRAN

0.3

0.3 0.2 HZZ4

0.2

NASTRAN

0.1

HZZ1

0.1

HZZ2 0

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

0.8

1

0 –3.05

–3.045

–3.04

–3.035

τ13(L/2,L/2,x3)

–3.03

–3.025

–3.02

Fig. 12. s13(L/2, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

–3.015 x 10–4

V2(L,L,x3)

Fig. 14. V2(L, L, x3) through-the-thickness distribution (Laminate 2, UNI load, undamaged).

1 HZZ1 0.7

0.9

HZZ3 NASTRAN

0.8 0.6

0.7 0.5

0.5

x3

x3

0.6

0.4

0.4

0.3

0.3 0.2

HZZ3

0.2

NASTRAN 0.1

HZZ1

0.1

HZZ2 0

–1

–0.9

–0.8

– 0.7

– 0.6

–0.5

–0.4

–0.3

–0.2

–0.1

0

τ23(L/2,L/2,x3)

Fig. 13. s23(L/2, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, undamaged).

0 –47

–46

–45

–44

–43

σ33(0,L/2,x3)

–42

–41

–40

–39

Fig. 15. r33(0, L/2, x3) through-the-thickness distribution (Laminate 2, UNI load, undamaged).

144

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

HZZ3 0.7

0.7

HZZ1 NASTRAN

0.6

0.6

HZZ1 0.5

0.4

0.4

NASTRAN

x3

x3

HZZ3 0.5

0.3

0.3

0.2

0.2

0.1

0.1

0 –0.5

–0.4

–0.3

– 0.2

–0.1

0

0.1

0.2

0.3

0.4

0 0.9

0.5

1

1.1

T

½M ½V  ½V  ½0

#

f€ qg f€ kg



"

½K þ ½V 

T

½V  ½0

#

fqg fkg



 ¼

fF g fAg



where {k} is the vector of the Lagrange multipliers. This is the method used to analyse damaged plates in the present paper. We have limited our discussion to sliding type damages hence the numerical results presented for damaged plates concern this kind of damage. Nevertheless, the Hermitian Zig-Zag model, together with the sub-laminates approach,

4. Numerical results Some results of static analyses on undamaged and damaged plates subjected to mechanical and thermal loads and temperature variations are shown and discussed; the aim of these studies is to demonstrate the capabilities of the proposed model and the related plate finite element.

HZZ3 0.7

HZZ1 NASTRAN

0.6

x3

0.5

0.4

0.3

0.2

0.1

–20

–15

–10

–5

1.4 x 10–3

could be used to describe also opening mode damages (discontinuities of transverse displacement) or other more complicated fracture patterns.

ð11Þ

0

1.3

Fig. 17. V3(L, L/2, x3) through-the-thickness distribution (Laminate 2, LIN load, undamaged).

Fig. 16. s13(L/2, L/2, x3) through-the-thickness distribution (Laminate 2, UNI load, undamaged).

"

1.2

V3(L,L/2,x3)

τ13(L/2,L/2,x3)

0

5

10

15

20

σ33(0,L/2,x3)

Fig. 18. r33(0, L/2, x3) through-the-thickness distribution (Laminate 2, LIN load, undamaged).

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

4.1. Numerical data The mechanical properties of the adopted materials and the stacking sequence details for the considered laminates are listed in Tables 1 and 2, respectively. In all cases we have considered thick square plates, with a side length-tothickness ratio of 4. The plate is clamped along the edge x1 = 0 and free along the others.

145

Here below are listed the different mechanical load distributions used in the presented study, individually or combined: • transverse unit load uniformly distributed over the upper face and directed downwards (‘‘TR’’, Fig. 3); • tangential unit load uniformly distributed over the top and bottom faces, on the upper face directed according

Fig. 19. V1 overall distribution (Laminate 1, TR load, delamination A).

Fig. 20. V3 overall distribution (Laminate 1, TR load, delamination A).

146

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

to the positive x1-axis and on the lower face according to the positive x2-axis (‘‘TG’’, Fig. 3). For the thermal analyses, we have considered two through-the-thickness temperature distributions (no variations along x1 and x2): • uniform, H = + 100 C (‘‘UNI’’); • linear in x3, H = 50 C at the bottom face and H = +50 C at the top face (‘‘LIN’’). When analysing damaged plates, we will consider four possible delaminations, denoted by a capital letter (A, B,

C or D) in Fig. 4. If not specified, it is supposed that the delamination is characterised by a complete debonding at the damaged interface. Once again, we remind that only axial displacements exhibit a jump at the damaged interface, the transverse displacement being continuous. For a list of the numerical applications presented in this paper (stacking sequence, mechanical and/or thermal loads, damages) and of the relative figures, see Table 4. The results obtained with the present approach (Hermitian Zig-Zag model + sub-laminate approach) will be accompanied by the HZZn label, where n is the number of sub-laminates; if needed, a further description of how

1

1 HZZ2

0.9

HZZ2 0.9

HZZ4

HZZ4 NASTRAN

0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

NASTRAN

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 –1.5

–1

–0.5

0

0.5

1

1.5

2

2.5

V1(L,L,x3)

0 –15

3

Fig. 21. V1(L, L, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

–10

–5

0

5

10

15

σ11(3L/4,L/2,x3)

x 10–3

Fig. 23. r11(3L/4, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

1 1

HZZ2 0.9

HZZ4 0.9

NASTRAN 0.8

0.8

0.7

0.7

0.6

0.5

x3

x3

0.6

0.4

0.5

0.4

0.3

0.3

0.2

0.2

0.1

0.1

HZZ2 HZZ4 NASTRAN

0 –2

–1

0

1

2

3

V2(L,L,x3)

4

5

6

7

8 x 10–5

Fig. 22. V2(L, L, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

0 –0.0176 –0.0175 –0.0174 –0.0173 –0.0172 –0.0171

–0.017

–0.0169 –0.0168 –0.0167 –0.0166

V3(L,L,x3)

Fig. 24. V3(L, L, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

4.2. Undamaged plates

1 HZZ2 HZZ4

0.9

NASTRAN 0.8

0.7

x3

0.6

0.5

0.4

0.3

0.2

0.1

0

147

–2

–1.5

–1

–0.5

0

τ13(3L/4,L/2,x3)

Fig. 25. s13(3L/4, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

the laminate has been divided in sub-laminates will be added. The results will be shown both by means of through-the-thickness distributions of displacements, strains and stresses and by means of fringe representations of the displacement fields on the undeformed shape of the plate; in this second case, the representation is based on a sub-laminate per layer. Figures showing through-the-thickness distributions adopt the following units: (mm) for the thickness coordinate and for displacements, (MPa) for stresses (strains are non dimensional).

One of the capabilities of the Hermitian Zig-Zag plate model is that it allows to satisfy the traction equilibrium condition on the top and bottom laminate faces also when tangential loads are present. This is the case of the static analysis performed on undamaged laminate 1 (TG loads). Considering that there not exists an analytical solution for this problem, the reference results have been obtained by means of a detailed FEM analysis conducted with MSC/NASTRAN. In Fig. 5 the adopted mesh is represented: 40 · 40 · 12 HEXA8 elements have been used for the plates with Laminate 1 stacking sequence (three elements along the thickness of each of the four layers). As for the plates with Laminate 2 stacking sequence, the mesh is again 40 · 40 · 12 because four elements along the thickness of each of the three layers are used. As for the HZZ FEM solution, the mesh is 8 · 8 · n over the whole plate (see Table 3 for more details on the number of finite elements and of degrees of freedom). The plate exhibits a remarkable transverse displacement also if loaded in the x1- and x2-axis (Figs. 6–8). The HZZ approach with an adequate number of sub-laminates is able to predict the transverse deflection also when its through-the-thickness behaviour is far not only from a constant distribution but also from a linear one (Figs. 6–8). The in-plane displacement in the x1-direction is well captured also using one sub-laminate (Fig. 9). Figs. 10–13 show the capability of the approach to evaluate transverse shear strains and stresses when tangential loads are applied to the external top and bottom surfaces of the laminate; the values of s13 and s23 (Figs. 12 and 13) on the top and

Fig. 26. V1 overall distribution (Laminate 1, TG load, delamination A).

148

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

Fig. 27. V2 overall distribution (Laminate 1, TG load, delamination A).

1

1

HZZ2

HZZ2 0.9

HZZ4

HZZ4

0.9

NASTRAN

0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

NASTRAN

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 –8

–6

–4

–2

0

V1(L,L,x3)

2

4

6 x 10–4

0 1.5

2

2.5

3

3.5

V2(L,L,x3)

4

4.5

5

5.5 x 10–3

Fig. 28. V1(L, L, x3) through-the-thickness distribution (Laminate 1, TG load, delamination A).

Fig. 29. V2(L, L, x3) through-the-thickness distribution (Laminate 1, TG load, delamination A).

bottom faces correspond to the values of the external loads or are zero if no load is applied. As explained in Section 2, the formulated plate mixed finite element can also be used to analyse the static response of laminates under a known temperature profile. Some results concerning thermal loading of composite plates have already been presented [16]; herein we briefly summarize the main ones. At first we consider an undamaged square plate (Laminate 2) subjected to a uniform (UNI) and to a through-the-thickness linear (LIN) temperature variation. Also in this case the reference solution has

been obtained via a MSC/NASTRAN detailed 3D model (Table 3). Figs. 14–16 show some results in the UNI case. The inplane displacement V2(L, L, x3) (Fig. 14) is well captured also with one sub-laminate but only by using three sublaminates its through-the-thickness gradient is correct. As for transverse stresses, one sub-laminate gives the average value of the normal one (Fig. 15) while the shear one is incorrect (Fig. 16): only by using three sub-laminates the transverse shear stress is correctly evaluated. When considering the LIN temperature distribution (Figs. 17 and 18)

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155 1

1 HZZ2

HZZ2 HZZ4

0.9

HZZ4

0.9

NASTRAN

0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

NASTRAN

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 1.5

2

2.5

3

3.5

4

4.5

5

5.5

V3(L,L,x3)

0 –0.8

6

–0.6

–0.2

0

0.2

0.4

0.6

0.8

1

τ13(3L/4,L/2,x3)

Fig. 32. s13(3L/4, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delamination A).

1

1

HZZ2

HZZ2 0.9

HZZ4

0.9

HZZ4

NASTRAN

NASTRAN 0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

–0.4

x 10–4

Fig. 30. V3(L, L, x3) through-the-thickness distribution (Laminate 1, TG load, delamination A).

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

149

–10

–5

0

5

10

σ11(3L/4,L/2,x3)

0 1

–0.9

–0.8

–0.7

–0.6

–0.5

–0.4

0.3

0.2

0.1

0

τ23(3L/4,L/2,x3)

Fig. 31. r11(3L/4, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delamination A).

Fig. 33. s23(3L/4, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delamination A).

the transverse normal deformability is not negligible and the transverse deflection exhibits a relevant through-thethickness variation, well captured with a sub-laminate per layer (Fig. 17).

detailed description of how the HZZ approach models delaminations). The reference solution is obtained once again with MSC/NASTRAN. In this case the delamination is modelled meshing separately the two solids representing the adjacent layers divided by the damage, merging the coinciding nodes in the undamaged area of the interface and ‘‘connecting’’ the coinciding nodes of the damaged area with rigid elements enforcing the continuity of the transverse displacement only. Figs. 19–25 show some results for the TR load case and with the delamination A. In Fig. 19 the overall distribution of V1 displacement is shown while Fig. 20 gives the distribution of the transverse displacement. Note the presence of a discontinuity in the in-plane displacement where

4.3. Damaged plates Let us now consider the static response of damaged plates (Laminate 1, delaminations A and A + B, Fig. 4) subjected to mechanical loads (TR and TG). It is supposed that there is a complete debonding at the damaged interfaces; the damage is of the sliding type, so a jump in the axial displacements V1 and V2 is present while V3 is continuous through the debonded interface (see Section 3 for the

150

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

Fig. 34. V1 overall distribution (Laminate 1, TR load, delaminations A + B).

Fig. 35. V2 overall distribution (Laminate 1, TR load, delaminations A + B).

the delamination is located (Fig. 19). It can be seen that the use of two sub-laminates gives accurate results in terms of axial displacements and in-plane normal stress (Figs. 21–23) while four sub-laminates lead to a remarkable improvement for transverse quantities (Figs. 24 and 25). We highlight that the HZZ approach correctly predicts that at a damaged interface, where a total debonding has been reached, the transverse shear stress vanishes

(Fig. 25). Considering this modelling issue, some words must be spent on FEM commercial codes, MSC/NASTRAN in our case. At the damaged interfaces (Fig. 25) the model has two sets of coinciding nodes belonging to the two layers separated by the delamination. Thus, when plotting the transverse shear stresses, we obtain discontinuities that are a violation of interlaminar equilibrium. A common technique adopted to overcome this drawback,

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155 1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

151

HZZ3 HZZ4

0.5

x3

x3

NASTRAN

0.5

0.4

0.4

0.3

0.3

0.2

0.2 HZZ3

0.1

HZZ4

0.1

NASTRAN 0 1.5

1

0.5

0

0.5

1

1.5

2

2.5

3

V1(5L/8,L/2,x3)

0 –25

3.5 x 10–3

Fig. 36. V1(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delaminations A + B).

–20

–5

0

5

10

15

20

Fig. 38. r11(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delaminations A + B).

1 HZZ3

HZZ3 HZZ4

0.9

HZZ4

0.9

NASTRAN

NASTRAN

0.8

0.8

0.7

0.7

0.6

0.6

0.5

x3

x3

–10

σ11(5L/8,L/2,x3)

1

0.4

0.5

0.4

0.3

0.3

0.2

0.2

0.1

0 –0.0123

–15

0.1

–0.0122

–0.0121

–0.012

–0.0119

–0.0118

–0.0117

–0.0116

V3(5L/8,L/2,x3)

0

–3.5

–3

–2.5

–2

–1.5

τ13(5L/8,L/2,x3)

–1

–0.5

0

Fig. 37. V3(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delaminations A + B).

Fig. 39. s13(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delaminations A + B).

is to calculate an average value at the common location; looking at Fig. 25 one understands that also in this case the final result would not be accurate. Similar comments may be proposed for the case with delamination A but TG load (Figs. 26–33). In particular, we focus the attention on the transverse shear stresses (Figs. 32 and 33); their values on top and bottom laminate faces perfectly recover the value of the external tangential loads. Let us now consider what happens with TR load applied to a plate with two delaminations (A + B), Figs. 34–39. In this case, the solution obtained with three sub-laminates is based on the following thickness division: layer 1, layer 2, layers 3–4. The capability of the HZZ approach to describe the through-the-thickness distribution of in-plane

displacements is again confirmed (Fig. 36); at a (x1, x2) location where both delaminations are present, the two in-plane displacement jumps are well captured. As for stresses (Figs. 38 and 39) , results appear in good agreement with the reference ones when using four sub-laminates. When considering both delaminations (A + B) and tangential loads (TG), once again the use of more sublaminates (4 instead of 3) appears to be useful in particular for transverse quantities (Figs. 44, 46 and 47) while inplane displacements and stresses are well captured also with three thickness subdivisions (Figs. 42, 43 and 45). Figs. 40 and 41 show how the two delaminations affect the global response of the laminate in terms of in-plane displacements.

152

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

Fig. 40. V1 overall distribution (Laminate 1, TG load, delaminations A + B).

Fig. 41. V2 overall distribution (Laminate 1, TG load, delaminations A + B).

Having assessed that HZZ approach is able to describe both undamaged and damaged plates (where totally debonded interfaces are present), the same model can be used to study the effect of a more and more severe damage, i.e. a delamination in which a weakening adhesive effect is present (Figs. 48–50). Let us imagine that the usual square plate (Laminate 1, TR, HZZ4 scheme adopted) is damaged according to pattern A (Fig. 4), that Eq. (6) holds at the

damaged interface and that, in our case, Tab = Tdab. Figs. 48–50 show what happens for T ranging from 0 (undamaged interface) to 1 (corresponding to sa3 = 0, i.e. totally debonded interface). The in-plane displacement jump increases (Fig. 48), the transverse deflection becomes bigger due to the loss in overall plate stiffness (Fig. 49) and the transverse shear distribution becomes that of two separate two-layered plates (Fig. 50).

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155 1

153

1 HZZ3

HZZ3 0.9

HZZ4

0.9

HZZ4

NASTRAN

0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

NASTRAN

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 –4

–2

0

2

4

6

8

10

12

V1(5L/8,L/2,x3)

0 –2.8

14

–2.78

–2.76

–2.74

x 10–4

Fig. 42. V1(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delaminations A + B).

–2.72

–2.7

–2.68

–2.66

–2.6 x 10–3

1

HZZ3

HZZ3

HZZ4

0.9

HZZ4

0.9

NASTRAN

NASTRAN

0.8

0.8

0.7

0.7

0.6

0.6

0.5

x3

x3

–2.62

Fig. 44. V3(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delaminations A + B).

1

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

–2.64

V3(5L/8,L/2,x3)

1

1.5

2

2.5

3

V2(5L/8,L/2,x3)

3.5

4

4.5

5 x 10

–3

0 –15

–10

–5

0

5

σ11(5L/8,L/2,x3)

10

15

20

Fig. 43. V2(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delaminations A + B).

Fig. 45. r11(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delaminations A + B).

An interesting analysis is that concerning damaged plates subjected to a temperature variation ([16], Laminate 2, UNI, Figs. 51 and 52). In Fig. 51 it is shown what happens to the in-plane displacement in x1-direction when delamination C is present; the through-the-thickness distribution and the intensity of the jump due to debonding is correctly predicted by the HZZ3 model if compared with MSC/NASTRAN solution. In Fig. 52 it is shown the effect of both C and D delaminations on the V1(3L/8, L/2, x3) in-plane displacement.

laminates strategy is presented. The article aims to assess the modelling capabilities of the proposed approach also in presence of non-classical effects like: (i) tangential loads applied to the top and bottom faces, (ii) damaged interfaces and (iii) thermal effects (known temperature distributions). A lot of numerical examples are presented where damaged and undamaged thick laminates are subjected to mechanical or thermal loads. Different damage geometries and severities are tested, also considering multiple delaminations and complete debonding between layers. In each case the FEM approach has shown to be accurate and computationally effective with respect to FEM commercial codes. The main advantage of the proposed

5. Conclusions A FEM approach based on a higher order lamination theory (Hermitian Zig-Zag model) and on the sub-

154

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155 1

1 HZZ3

Undamaged

HZZ4

0.9

Damaged (T=1e–4)

0.9

Damaged (totally debonded)

0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

NASTRAN

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0 –0.8

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

0.8

0 –1.5

1

Damaged (T=5e–4)

–1

–0.5

τ13(5L/8,L/2,x3)

0

0.5

1

1.5

2

2.5

Fig. 46. s13(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delaminations A + B).

3 x 10–3

V1(L,L,x3)

Fig. 48. V1(L, L, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

1

1 HZZ3

0.9

HZZ4

0.9

0.8

0.8

0.7

0.7

0.6

0.6

x3

x3

NASTRAN

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

Undamaged Damaged (T=1e–4) Damaged (totally debonded)

0 –1

–0.9

–0.8

–0.7

–0.6

–0.5

–0.4

– 0.3

–0.2

–0.1

0

τ23(5L/8,L/2,x3)

0 –0.018

Damaged (T=5e–4)

–0.017

–0.016

–0.015

–0.014

–0.013

–0.012

–0.011

V3(L,L,x3)

Fig. 47. s23(5L/8, L/2, x3) through-the-thickness distribution (Laminate 1, TG load, delaminations A + B).

Fig. 49. V3(L, L, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

approach with respect to commercial codes is the reduced number of elements and of degrees of freedom needed to obtain a similar or even better description (vanishing of transverse shear stresses on external faces and on totally debonded interfaces, for example) of multilayered plates response; the reduced computational weight of the approach would allow to obtain also reduced CPU times when implemented in a more commercial programming environment. Some guidelines for modelling have been presented, being shown, as a general rule, that for in-plane displace-

ments and stresses a reduced number of sub-laminates is enough while for transverse quantities the use of more sub-laminates improves the model accuracy. The proposed approach is thus suggested for the analysis of very thick and/or highly anisotropic multilayered plates, especially when non-classical effects are present (tangential loads, damaged interfaces, thermal effects). The modelling freedom guaranteed by the sub-laminates philosophy gives some interesting perspectives for future further developments, among which: modelling of opening mode damages, formulation of a completely coupled

M. Gherlone, M. Di Sciuva / Composite Structures 81 (2007) 137–155

155

1 HZZ3

Undamaged 0.7

Damaged (T=1e–4)

0.9

NASTRAN

Damaged (totally debonded) 0.8

0.6

Damaged (T=5e–4)

0.7 0.5

x3

x3

0.6

0.5

0.4

0.3

0.4

0.3 0.2 0.2 0.1 0.1

0 –6

–5

–4

–3

–2

–1

0

0 –2.3

τ13(3L/4,L/2,x3)

–2.2

–2.1

–2

–1.9

V1(3L/8,L/2,x3)

Fig. 50. s13(3L/4, L/2, x3) through-the-thickness distribution (Laminate 1, TR load, delamination A).

–1.8

–1.7 x 10–4

Fig. 52. V1(3L/8, L/2, x3) through-the-thickness distribution (Laminate 2, UNI load, delaminations C + D).

References HZZ3 (damaged) 0.7

NASTRAN (damaged) HZZ3 (undamaged) NASTRAN (undamaged)

0.6

x3

0.5

0.4

0.3

0.2

0.1

0

–5.88

–5.86

–5.84

–5.82

V1(L,L/2,x3)

–5.8

–5.78

–5.76 x 10–4

Fig. 51. V1(L, L/2, x3) through-the-thickness distribution (Laminate 2, UNI load, delamination C).

thermo-mechanical problem, application to active control problems with adoption of piezo-layers and to the analysis of multilayered plates with functionally graded materials.

Acknowledgements This research is part of the Scientific Research Programs of National Interest, titled ‘‘Damage tolerance design of stiffened and sandwich composite panels for fuselages: modelling, experimental testing and monitoring’’ funded by the Italian Ministry of University and Research (MIUR), that the authors gratefully acknowledge.

[1] Librescu L. Elastostatics and kinetics of anisotropic and heterogeneous shell-type structures. Leyden (Netherlands): Noordhoff International Publications; 1975. [2] Di Sciuva M, Gherlone M, Librescu L. Implications of damaged interfaces and of other non-classical effects on the load carrying capacity of multilayered composite shallow shells. Int J Non-Linear Mech 2002;37:851–67. [3] Schmidt R, Librescu L. A general theory of geometrically imperfect laminated composite shells featuring damaged bonding interfaces. Quart J Mech Appl Math 1999;25(1):1–19. [4] Savoia M, Reddy JN. Three-dimensional thermal analysis of laminated composite plates. Int J Solids Struct 1995;32(5):593–608. [5] Di Sciuva M. Bending, vibration and buckling of simply supported thick multilayered orthotropic plates: an evaluation of a new displacement model. J Sound Vib 1986;105:425–42. [6] DiSciuvaM.Animprovedshear-deformationtheoryformoderatelythick multilayered anisotropic shells and plates. J Appl Mech 1987;54:589–96. [7] Di Sciuva M. Multilayered anisotropic plate models with continuous interlaminar stresses. Compos Struct 1992;22(3):149–67. [8] Di Sciuva M. Geometrically non-linear theory of multilayered plates with interlayer slips. AIAA J 1997;35(11):1753–9. [9] Di Sciuva M, Gherlone M. A global/local third-order Hermitian displacement field with damaged interfaces and transverse extensibility: analytical formulation. Compos Struct 2003;4:419–31. [10] Gherlone M, Di Sciuva M. Thermo-mechanics of undamaged and damaged multilayered composite plates: a sub-laminates finite element approach.ComposStruct,inpress,doi:10.1016/j.compstruct.2006.08.004. [11] Di Sciuva M, Gherlone M. A global/local third-order Hermitian displacement field with damaged interfaces and transverse extensibility: FEM formulation. Compos Struct 2003;4:433–44. [12] Cheng ZQ, Jemah AK, Williams FW. Theory for multilayered anisotropic plates with weakened interfaces. J Appl Mech 1996;63(4):1019–26. [13] Williams TO. Efficiency and accuracy considerations in a unified plate theory with delamination. Compos Struct 2001;52(1):27–40. [14] Lu X, Liu D. An interlaminar shear stress continuity theory for both thin and thick composite laminates. J Appl Mech 1992;59:502–9. [15] Averill RC. Static and dynamic response of moderately thick laminated beams with damage. Compos Eng 1994;4(4):381–95. [16] Gherlone M, Di Sciuva M. A higher-order FEM approach for the analysis of damaged composite multilayered plates subjected to thermal loads. In: Proceedings of the VI international congress on thermal stresses, Wien (Austria) 2005. p. 253–6.