A hierarchical multiple plate models theory for laminated composites including delamination and geometrical nonlinear effects

Composite Structures 93 (2011) 780–791

Contents lists available at ScienceDirect

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

A hierarchical multiple plate models theory for laminated composites including delamination and geometrical nonlinear effects S.L. Angioni, A. Visrolia, M. Meo ⇑ Department of Mechanical Engineering, University of Bath, BA2 7AY Bath, UK

a r t i c l e

i n f o

Article history: Available online 15 September 2010 Keywords: Laminated plates Unified plate theory Hierarchical modeling Mesh superposition Delaminations

a b s t r a c t Modeling and predicting the response and failure of laminated composite structures is a very challenging task. On one hand the designer is interested in the overall response of the structure, but on the other he is also interested in modeling phenomena such as damage which have a localized behavior. The ideal solution would be to carry out the analysis at a global level enhancing the model by adding localized information when needed. Plate models dimensionally reduce a finite element method (FEM) problem by imposing constraints on the through thickness variation of the displacement field. In this work, the displacement field is represented as the superposition of a number of plate fields, and a unifying strain field is derived. By appropriately defining boundaries to the ‘enhancing’ displacement fields, it is demonstrated that the superposition of displacement fields can be used to locally enrich the solution where more information is required. In this manner, an efficient global model can be used to determine gross displacements, and an enriching model can be used to determine stresses at lamina interfaces for the accurate prediction of localized phenomena. The proposed theoretical framework would allow an analyst to study both thin and thick walled laminated structures by introducing simple switch-over criteria between single and multiple layer theories. The model is implemented using an extended FEM (X-FEM) and superposition approach. Extra degrees of freedom are added to the model to represent the additional displacement fields, and the meshing process remains independent for each field. Both equivalent single layer theories and discrete layer theories are incorporated into the general unified framework. Geometric nonlinearities are also included in the form of the von Kármán equations. Delaminations are integrated by performing an extrinsic enrichment of the assumed displacement fields. A cohesive law between laminae was also incorporated to calculate the tensile and shear tractions. The present paper outlines the various general theories to describe in-plane and out-plane displacements fields, then a unified theoretical framework is formulated and finally the finite element formulation which is based on a multilevel mesh superposition approach is reported. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction 1.1. Motivation Key components of primary aircraft structures, such as wing skins, are currently being designed using laminated composites. Although the structural designer is interested in the global behavior of these structures, very often phenomena such as impact loading, varying geometrical features, such as tapering, etc. have localized structural effects. ⇑ Corresponding author. E-mail addresses: [email protected] (S.L. Angioni), [email protected] (A. Visrolia), [email protected] (M. Meo). 0263-8223/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2010.08.003

Plate models dimensionally reduce a FEM problem by imposing constraints on the through thickness variation of the displacement field. Unfortunately if a single plate model is used throughout the structure this will not be able to provide detailed information at all geometric scales. Furthermore, in large domain design problems using FEM increasing the number of nodes in the mesh in order to gain localized information is computationally inefficient. To accurately describe localized phenomena, different sub-regions of the structure can be analyzed using different plate models. Another methodology is to resort to hierarchical models that offer globallocal analysis capabilities. This is the approach adopted in this paper. In the following we will use the expressions plate model and assumed displacement field interchangeably.

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

1.2. Mathematical models for laminated composite plates A number of reviews [1–4] have been conducted into the various mathematical theories available for the analysis of composite laminate plates. In general, these models can be separated into two categories: (i) equivalent single layer (ESL) or smeared models and (ii) discrete layer (DL) or layerwise theory (LWT) models. ESL models assume a continuously differentiable (C1) displacement field in the through thickness direction (referred to in this paper as the z-direction). This plate-style analysis relies on homogenization (or smearing) of the laminate properties to determine overall structural response of the structure. Examples are classical lamination theory (CLT) and first-order shear deformation theory (FSDT). For thin walled structures, these models are sufficient for determining a global structural response, however, they are not capable of accurately resolving interlaminar stresses. Higher order or refined theories have been used in ESL models, such as the third order shear deformation theory (TSDT) [5], to improve stress resolution. However these models fail to address the cause of the models’ inability to accurately resolve interlaminar stresses. The alternative to ESL models are the LWT models. The layerwise model represents the displacement field as continuous but piecewise in the z-direction, i.e. C0. A number of authors present such a model, including Reddy [6]. Common approaches to discretizing the through thickness displacement field are to use piecewise Lagrangian interpolations to represent physical lamina [7–9], or to use a ‘zig-zag’ field [10]. The ESL models offer a computationally efficient method, whereby the number of degrees of freedom in the model is independent of the number of lamina, and, especially in cases of thin laminates, can offer accurate modeling of the global displacements. The layerwise models offer a solution for thick laminates, at the cost of increased complexity: the number of degrees of freedom in the model are proportional to the number of lamina in the material, often numbering in the hundreds. Carrera and Ciuffreda [11] conduct an extensive review of these models using various orders of polynomials, and assess their performance under a variety of transverse loading conditions. The review is extended in [12] to include analysis of functionally graded materials, using Carrera’s unified formulation. The ability to combine any arbitrary plate models (including multiple ESL and layerwise models) would enable modeling of complex geometries, such as homogenized periodic microstructure, and ‘smart’ components, such as piezo-electric inserts and shape-memory alloy wires. 1.3. Multiple plate models theories A number of authors have proposed multiple plate models theories for laminated composites by combining the various mathematical theories available to model laminates. Barbero and Reddy created the generalized plate theory [7,13] which combines an ESL theory with a layerwise expansion for the in-plane displacement components. This is defined as a partial layerwise theories by Reddy [5]. The theory can be augmented with a discontinuous field using a step function, as explained in [8], in order to model delamination. Robbins and Reddy [4,9] proposed a hierarchical model called variable kinematic model (VKM) that combines FSDT and layerwise displacement fields, as well as a discontinuous field for delamination. This is done by superposition of meshes in regions of the model where more accurate information is required. This VKM is expanded in [14] by making the model automatically adaptive. The use of Murakami’s zig-zag function [10] is always used as an enhancement to a simpler model, such as FSDT, and can itself be seen as an example of the unification of two different displacement fields.

781

Williams [15] proposes a model referred to as a two length scale or global-local approach, by introducing a coupling between the displacement fields used on the two different length scales. The through the thickness displacement function is of arbitrary selection. Delamination can be incorporated using an arbitrary choice of interfacial constitutive relationships. This global-local strategy is presented in its linear form in [16], and proves to be both accurate and efficient. It is worth noting that most authors (e.g. [17]) speak about multiscale modeling only when local information is present in the model along all the axis and not only the z-axis. 1.4. Hierarchical plate models The hierarchical modeling of composite structures is achieved through the use of sequence of plate models, each identified by specific parameters, of which the exact solutions constitute a converging sequence to the solution of the fully three-dimensional model. This modeling approach uses a predefined methodology for selecting the model with the appropriate level of sofistication from the sequence. Cho and Oden [18] applied the concept of hierarchical modeling to ESL and DL plates and shells identifying as relevant parameter the order of the displacement field in the through the thickness direction. However, when they applied this concept to laminated structures, the authors did not take into account the microstructure of the laminae. As pointed out by Actis et al. [19] the concept of hierarchical plate models must not be confused with that of hierarchical finite element spaces. Hierarchical finite element spaces produce a converging sequence of approximate solutions, the limit of which is the exact solution of the particular mathematical model of interest. Hierarchical plate models provide means for controlling modeling errors while hierarchical finite element spaces provide means for controlling discretization errors. Zienkiewicz and Taylor [20] present a review of hierarchical finite element spaces methods. 1.5. Mesh superposition techniques Of particular interest to this work is the concept of multilevel mesh superpositioning, which was first introduced by Fish [21] as the s-version or s-refinement of the finite element method, where the s stands for superpositioning. This technique was first developed for the two level case where a coarse mesh (global mesh) is superpositioned by an independent, refined mesh (overlay mesh or local mesh) in order to form a composite mesh. The total displacement field on the composite mesh is the sum of the displacement field interpolated on the global mesh (global displacement) and the displacement field interpolated on the local mesh (local displacement). The important aspects of this method is that the original mesh does not need to be changed, global and local meshes do not need compatible discretizations and the local mesh can enhance the global mesh regardless of the original global mesh topology. The mesh superpositioning is called structured when the superpositioning elements fit entirely within the domain of the superpositioned element. Otherwise the technique is called unstructured. The technique was extended to the multilevel case in [22–24] and to handle discontinuous fields in [25]. The method was applied to laminated composites in [26,27]. In these applications the displacement fields at the various levels are all based on the same assumed displacement field, hence the sub-regions differ only in the level of refinement of the interpolated solution. The approach used to tackle cracks in the s-version FEM is to mesh the superpositioning elements so that they follow the shape of the crack, and to double the number of nodes within the elements on the crack surface in order to represent the discontinuity.

782

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

1.6. Modeling of delamination using enrichment methods Another method for handling discontinuities and cracks is the extended finite element (X-FEM) or generalized finite element [28], which, using the partition of unity property of finite elements [29], allows local enrichment functions to be easily incorporated into a finite element approximation for modeling discontinuities without requiring remeshing. There are two ways of carrying out the enrichment [30]: enriching the shape functions (intrinsic enrichment) or adding to the problem additional unknowns or degrees of freedom associated to the enriched solution (extrinsic enrichment). As discussed by Belytschko et al. [28], the s-version FEM is closely related to the X-FEM. Lee et al. [31] combined X-FEM and s-version FEM for modeling cracks proposing a hybrid method called XS-FEM. Local enrichment can be used for modeling delamination in laminated composites, as delamination can be considered as a particular form of fracturing in which the planes of fracturing are known a priori [32].

to be h, as shown in Fig. 1. Moreover, the midplane is assumed to be an open and bounded region x 2 R2 with piecewise smooth boundary @ x, so that the total domain X = x  (h/2, h/2). The boundary oX then consists of the union of the plate’s top and bottom surfaces oXu = x  { h/2, h/2} and by the lateral surface oXt = @ x  (h/2, h/2). A typical point P in X is denoted by x ¼ ðx; x3 Þ ¼ ðx; zÞ, where x 2 x and h/2 6 z 6 h/2. 





2.2. Total displacement field The total displacement field u in the proposed hierarchical multiple plate models approach is given in (2.1) by the superposition [22–24] in X of all the individual plate models ua being considered

uðx; y; z; tÞ ¼

M X

ua ðx; y; z; tÞ

ð2:1Þ

a¼1

where M is the total number of plate models. A material point occupying the position (x, y, z) in the undeformed laminate moves to the position (x + u, y + v, z + w) in the deformed laminate, where u ¼ u^ex þ v ^ey þ w^ez .

1.7. Scope of this work

2.3. Dimensional reduction

While previous authors limit the analysis to only two levels of superposition [4,9] or only to specific plate models [26,27], this paper proposes a hierarchical multiple plate models approach for modeling structural response of thick laminated composite plates in which any number and type of plate models can be superimposed by using a multilevel mesh superposition approach. Both equivalent ESL and DL plate models are incorporated into the general framework. Geometric nonlinearities are also included in the form of von Kármán equations. Discontinuous through the thickness displacement fields, such as delaminations, are integrated by performing an extrinsic enrichment of the assumed displacement fields. A cohesive law between laminae is used where appropriate. In implementation, an adaptive method may be adopted to automatically enhance the model in the appropriate regions to the appropriate level of complexity. When a higher level of precision (accurate stress distribution) is required, then the approximation is enhanced in a hierarchical fashion by adding the extra degrees of freedom needed to represent the additional plate models in the regions of interest. The outline of the paper is the following. Section 2 outlines the general theory of the hierarchical multiple models theory, while Section 3 describes the finite element formulation of the theory.

As discussed in Section 1, the full, three-dimensional elasticity problem is computationally expensive to solve. For this reason the problem is dimensionally reduced by semi-discretization. Each plate model a is obtained by constraining the admissible solutions of the three-dimensional problem to have a specific dependence on the transverse variable z. In particular, the problem is discretized only in the transverse direction z [33], while the spatial x and y and time t co-ordinates are left unchanged. Based on the above assumptions, the components of assumed displacement field for any plate model a can be expressed as follows:

uai ðx; y; z; tÞ ¼ U ai ðx; tÞF ai ðzÞ 

ð2:2Þ

2. Theory 2.1. Notation Let us consider a open and bounded domain X 2 R3 with piecewise smooth boundary oX. The plate is the closure X of X. The global (laminate) orthonormal Cartesian co-ordinate system of the ^1 ; n ^2 ; plate is denoted as (x1,x2,x3) = (x,y,z) with unit vectors ðn ^ 3 Þ ¼ ðn ^x ; n ^y ; n ^ z Þ. We further assume the z-axis is pointing upward n and a right-handed (or positive) system. The boundary oX is decomposed as oX = oXu [ oXt, where oXu is the boundary where the prescribed displacements are imposed (Dirichlet boundary conditions) and oXt is the boundary where the prescribed tractions are imposed (Neumann boundary conditions). In general, oXt is a curved ^ n ¼ nx n ^ x þ ny n ^ y defined to be outward surface with a unit normal n with respect to X. We will assume the xy-plane (reference plane) of the problem to coincide with the midplane (z = 0) of X and the total thickness of X

Fig. 1. Geometry of a laminated plate with curved boundary, including co-ordinate systems and layer numbering.

783

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

where i = 1, 2, 3. U ai is function only of the variables associated with the midplane x and time t, while F ai is function only of the trans verse variable z. 2.4. Material layers The composite laminated plate is composed of Nm orthotropic material lamina, as shown in Fig. 1. Each lamina k is of uniform thickness h(k) = z(k+1)  z(k), where z(k) is the co-ordinate of the interface between the kth and the (k + 1)th material layers. The local (layer or material) Cartesian  orthonormal  co-ordinate system ðkÞ ðkÞ ðkÞ of the kth layer is denoted as x1 ; x2 ; x3 ¼ ðxðkÞ ; yðkÞ ; zðkÞ Þ, with x(k) oriented at an angle h(k) with respect to the global co-ordinate x. The angle h(k) is aligned with the direction of the fibers in the lamina k. The kth material layer is defined in z(k) 6 z 6 z(k+1). The global and local co-ordinate systems are collinear along the z directions. 2.5. Out-of-plane discretization In order to determine the through the thickness stresses, the microstructure of the laminate, e.g. fibers and matrix constituents and fibers/matrix interface of each lamina, has to be accounted for. We will take into account only of the homogenized mechanical properties of each individual lamina. For each model a the thickness h of the plate can be discretized along the z axis into a number of mathematical layers N aL that is independent from the number of material layers Nm. A mathematical layer may consist of several material lamina (sub-layers), of a single material lamina or a sub-region of a material lamina (sub-lamina). Under the above assumptions, each mathematical layer can be represented by the equivalent homogeneous mechanical properties of the laminae that compose it. The layers may be assumed to be perfectly bonded together or have displacement jumps across interfaces due to delaminations. Formally, the out-of-plane discretization can be viewed as a partition Paz of the z axis into N aL intervals through the thickness of the plate.

    Pz ¼ z 2 Z aI jZ aI ¼ zaI ; zaIþ1 ; Z aI 2 ðh=2; h=2Þ; I ¼ 1; s; NaL a

aT

a

a

a

a

ð2:3Þ

a

S

T

i ¼ x . where xi xj = ;, if xi 2 Px ; xj 2 Px and i – j, and x The number of elements of Mi is called the cardinality of Mi and is denoted jMi j. The greatest element Si of Mi is the element of Mi which is greater than or equal to any other element of Mi . We define the in-plane domain xa with boundary @ xa of a model a as a subdomain of x which is the union of all subdomains S xi where the model a is utilized. Formally we can write, xa = xi for all subdomains xi such that a 2 Mi . The elements of the in-plane domain are x 2 xa .  We define the out-of-plane domain of a model a as any subset n o of intervals Z a1 ; Z a2 ; . . . ; Z aNa within the partition Paz , with N aI 6 N aL , I

where the model a is utilized. Some of these intervals may be contiguous, but in general this is not a requirement for the intervals. The elements of the out-of-plane domain are za 2 n o Z a1 ; Z a2 ; . . . ; Z aNa and zaI 2 Z aI . I n o The domain of a model a then is xa  Z a1 ; Z a2 ; . . . ; Z aNa , which in I general is equal to one or more disjoint volumes. For example, in the case of the superposition of a FSDT with a linear LWT model displayed in Fig. 2, as the FSDT theory can be represented by the superposition of two assumed displacement fields {1,2} and we refer to the linear LWT as model 3, the set of available mathematical models in X is M ¼ f1; 2; 3g. In subdomain x1 only the FSDT model is available, so M1 ¼ f1; 2g. In subdomain x2 both FSDT and linear LWT models are available, so M2 ¼ f1; 2; 3g. The in-plane domain of the FSDT model coincides with x, while the in-plane domain of the linear LWT model equal to x2. 2.7. Out-of-plane approximation functions a

In each interval Z aI the displacement variables U I cana be i approximated using a predefined approximation function F I of a i degree L I . For simplicity, if we assume that the components of i the displacement are approximated using the same functions then we can drop the index i and write a

LI a  l   X F aI zaI ¼ F I zaI l¼0

l

for zaI 2 Z aI . Note that here l is the exponential operator.

where Z I Z J ¼ ;, if Z I 2 Pz and Z J 2 Pz and I – J, and S a Z I ¼ ½h=2; h=2.   We define the measure (length) of subdomains Z aI as k Z aI . The  a out-of-plane discretization step is then equal to k Z I . The exact solution of the semi-discretized problem converges asymptotically to the continuous one [34] as the discretization step along z tends to zero. 2.6. Domain of a model It’s useful at this stage to introduce the concept of domain of a model. This concept will be used later for example for imposing compatibility conditions on the assumed displacement fields, i.e. boundary conditions between different model domains. As we have seen in Eq. (2.2), the dimensionally reduced displacement field a is the product of an in-plane function and an out-of-plane function. The domain of a model a is then the product of an in-plane domain and an out-of-plane domain. If the set of available mathematical models in X is defined as M ¼ f1; . . . ; Mg, then the domain x can be partitioned into P subdomains xi with boundaries @ xi, where within each subdomain xi only a subset of models Mi # M are used. This partition is called Px , and can defined formally as follows:

Px ¼ fxi jxi # x; 1 6 i 6 Pg

ð2:4Þ

Fig. 2. Superposition of a FSDT model with a linear LWT model.

ð2:5Þ

784

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

2.8. Number of nodes through the thickness

eyy The number of nodes through the thickness of the plate in each point x for a model a is equal to 

a

Nap ¼

NI X 

LaI þ 1



ð2:6Þ

I¼1

if all N aI intervals are contiguous. Otherwise the expression in (2.6) is increased by one for every non contiguous intervals. For example, in the case of the superposition of a FSDT with a linear LWT model displayed in Fig. 2, there is one node through the thickness for the FSDT model in x, while there are Nm + 1 nodes through the thickness for the linear LWT model in x2, if a number of subdivisions through the thickness equivalent to the number of material layers Nm is used for the LWT. 2.9. Layer semi-discretized displacement fields Now, introducing the concept of mathematical layers into definition (2.2), we can write a

Np   X  a a uai x; z; t ¼ U I x; t F I ðzÞ 

i

I¼1



ð2:7Þ

i

n o where i = 1,. . .,3, x 2 xa and z 2 Z a1 ; . . . ; Z aNa . I  In Eq. (2.7) the variables U aI are called the generalized displacements or kinematic variables for the assumed displacement field a, and are the unknowns of the problem. In the global  co-ordinate system (x, y, z) the components of U aI are U aI ; V aI ; W aI . F aI are 1-D a-priori selected approximation functions that are continuous within one or more contiguous mathematical layers. Here, for simplicity, the   functions are assumed to have the same components F aI ; F aI ; F aI along the x-, y- and z-directions. Then the same number of points through the thickness N ap can be used for the components of displacement field a Eq. (2.7) can be written as a

a

u ¼

Np X

a

a a

UI F I

v

a

¼

I¼1

Np X

a

a a

V I FI

a

w ¼

I¼1

Np X

W aI F aI

ezz

1 N ap a X @V a I @ ¼ F A @y I a¼1 I¼1 20 a 10 b 13 Np Np M X @W bJ b 1 X @W aI a A@X 4 @ F F A5 þ 2 a;b¼1 @y I @y J I¼1 J¼1 0 a 1 Np a X dF a I @ A ¼ WI dz a¼1 I¼1 M X

0 a 1 Np Nap a X X @W aI a A a dF I @ cyz ¼ VI þ F dz @y I a¼1 I¼1 I¼1 0 a 1 Np N ap a a X X dF @W a a I I @ czx ¼ UI þ F A dz @x I a¼1 I¼1 I¼1 M X

2 a 3 Np  a a M X X @U @V I 4 cxy ¼ þ I F aI 5 @y @x a¼1 I¼1 20 a 10 b 13 Np Np a M X X X @W bJ b @W a I 4@ þ F A@ F A5 @x I @y J a;b¼1 I¼1 J¼1

0¼

Z

T

ðdU þ dV  dKÞdt

ð2:11Þ

where dU is the virtual strain energy, dV is the virtual work done by the applied forces, dK is the virtual kinetic energy and T is the time interval. The virtual strain energy can be expressed as

Z n

o

rxx dexx þ ryy deyy þ rzz dezz þ ryz dcyz þ rzx dczx þ rxy dcxy dX ð2:12Þ

! Np

ð2:9Þ

Z

qdwdX

ð2:13Þ

X

dK ¼

Z

2.11. Displacements and strains The nonlinear strains associated with the total displacement field in Eqs. (2.2) and (2.7) can be calculated using the von Kármán equations [35]. The nonlinear strains are given in Eqs. (2.10).

1

Np X @U aI a A 1 @ þ ¼ F 2 @x I a¼1 I¼1 20 a 10 b 13 Np Np b a M X X X @W @W J a b I 4@  F A@ F A5 @x I @x J a;b¼1 I¼1 J¼1

_ w _ ÞgdX _ u_ þ v_ dv_ þ wd fq0 ðud

ð2:14Þ

where q0 is the density of the plate material and a dot on the variable indicates time derivative. For laminated plate structures the following relation stands true



exx

dV ¼ 

X

and is in general a function of x.

a

The virtual work done by the applied forces, in the case of a distributed load, is

The virtual kinetic energy is

a

a¼1

M X

ð2:10fÞ

0

dU ¼

The functional degrees of freedom (FDOF) or nodal degrees of freedom of the multiple models theory can be very useful when comparing different combinations of models. This parameter has the following definition

0

ð2:10eÞ

The equations of motion can be derived using the dynamic version of the principle of virtual displacements

2.10. Functional degrees of freedom

FDOF ¼ 3

ð2:10dÞ

2.12. Equations of motion

X

M X

ð2:10bÞ

ð2:10cÞ

M X

ð2:8Þ

I¼1

0

M X

Z

fg dX ¼

Z (Z x

X

¼

h 2

) fg dz dxdy

h2

Z (X Nm Z x

k¼1

z

ðkÞ

zt ðkÞ b

) fgdz dxdy

ð2:15Þ ðkÞ

ð2:10aÞ

where Nm is the number of material layers in the laminate, zb ¼ zðkÞ ðkÞ and zt ¼ zðkþ1Þ are the co-ordinates of the bottom and top of the k-layer. Defining the force resultants, the moment resultants and the mass moments of inertia as

785

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791



a

a

NI

a

NI

1

2

¼

a



a

T ¼

5

a

T

a

b I N

ab

ab

ab

Q IJ

Q IJ

Q IJ

1

2

6

ab

I IJ ¼

b I N

ns

Nm Z X

¼

nz

¼

f rzz

ðkÞ

Nm X

Z

ðkÞ

zt

k¼1



ðkÞ z b

Nm Z X

U I ¼  nx U aI þ ny V aI

a

ryz

zb

k¼1

T

ðkÞ

zt

ð2:16Þ



r^ nn r^ ns r^ nz T F aI dz

ðkÞ

zt ðkÞ

Substituting the expressions of dU, dV and d K in Eqs. (2.17), (2.18) and (2.19) in Eq. (2.11), collecting the coefficients of each general  ized displacement dU aI ; dV aI ; dW aI and setting these coefficients to zero separately over x, we obtain the following Euler–Lagrange equations of the theory

ryy rxy gT F aI F bJ dz

f rxx

zb

a

dU aI :

ðkÞ a b 0 F I F J dz

q

ðkÞ

a

@N I

@N I

1

þ

@x a

8 a Z X Np M
a

a

a

a

a

a

@dU I @dV I @dW I @dW I dU ¼ NI þNI þNI þNI 5 1 @x 2 @y 4 @y @x x a¼1 : I¼1 ) 

a @dU a a a a @dV aI I e I dW a þ N e I dV a þ N e I dU a dxdy þN þ þN I I I I 5 6 3 4 @y @x 8 a b " Np Z X Np M




Z

Z

M X

2 4

X

@ Xt a¼1

dK ¼

Z





Z X M   qb dW a1 þ qt dW aNap dxdy x a¼1





^ ns dus þ r ^ nz dwÞdzds ^ nn dun þ r ðr

2h

@ Xt

dV ¼ 

h 2

Nap

3 a a a a a b I dU I þ N b I dU I þ N b I dW a 5ds N I

I¼1

nn

n

ns

s

nz

I¼1

ð2:18Þ

b¼1

þ

@N I

6

@x

a

eI¼ þN 4

0 @

b

Np X

ab

I IJ

@U bJ

J¼1

@t2

1 A

0 b 1 Np ab @V b X J @ A I IJ @t2 b¼1 J¼1

M X

a

@N I

a

5 e I q  q þ þN b t 3 @x 8 b " ! Np M
@y

þ

ab @W b ab @W b @ J J Q IJ þ Q IJ @y @y @x 2 6

!#)

b

¼

Np M X X b¼1 J¼1

ab

I IJ

@W bJ @t 2

a

a

nn

nn

ð2:23Þ

a

a

ns

ns

b I ¼0 N I N a

a

nz

nz

b I ¼0 N I N

ð2:24Þ

where

J¼1

ð2:19Þ where qb ðxÞ and dwðx; h=2Þ are respectively the distributed stress   and the component along z of the displacement on the bottom of   h the plate z ¼  2 ; qt ðxÞ and dwðx; h=2Þ are respectively the distrib  uted stress and the component along z of the displacement on the   h top of the plate z ¼ 2 and q = qb + qt. ds is the arc-length of an ^ nn ; r ^ ns ; r ^ nz Þ are infinitesimal line element along the boundary. ðr the specified stress components on the portion oXt of the boundary, where n is the normal direction and s is the tangential direction on oX, as shown in Fig. 1. Note that if the unit outward normal vector is oriented at an angle h from the x-axis with cosine directors nx = cosh and ny = sinh, then the transformation between the co-ordinate systems (x, y, z) and (n, s, r) is

^ n þ sin hn ^s ^ x ¼ cos hn n

a

a

a

nn

1

6

a

a

a

ns

6

2

N I ¼ N I nx þ N I ny N I ¼ N I nx þ N I ny 2 b ! Np M ab @W b ab @W b X X J J 4 N I ¼ N I nx þ N I ny þ Q IJ þ Q IJ nx nz 5 4 @x @y 1 6 b¼1 J¼1 ! # ab @W b ab @W b J J ny þ Q IJ þ Q IJ @y @x 2 6 a

a

a

The primary variables (displacements) and secondary variables (forces) for the hierarchical multiple plate models theory are a

a

Primary variables : U I ; U I ; W aI

^ n þ sin hn ^s ^ y ¼  cos hn n ^r ^z ¼ n n

@N I

a

dW aI :

5

M X

b I ¼0 N I N

2 a b 3 Np Np M  ab  X X X b b b a a a _ dW _ 5dxdy 4 I IJ U_ I dU_ J þ V_ I dV_ J þ W I J

x a;b¼1

@y

eI¼ þN

6

@y a

@N I

2

a

Note that in the equations above integration-by-parts (spatial and time) was used to relieve the generalized displacements of any differentiation of the variational operator d, in order to use the fundamental lemma of calculus. Also the terms in x evaluated at t = 0 were set to 0 because the generalized displacements are zero there. The essential (geometric) boundary conditions for the theory are imposed on U aI ; V aI and W aI . The natural (force) boundary conditions for the theory are

Z h     i    dV ¼  qb x dw x; h=2 þ qt x dw x; h=2 dxdy x

a

dV I : a

a

ð2:22Þ

s

dF rzx g I dz dz T

and integrating through the thickness of the laminate, we can write

Z

ð2:21Þ

n

ðkÞ

zt

zb

k¼1

Nm Z X k¼1

b I N



U I ¼nx U aI þ ny V aI

ryy ryz rzx rxy gT F aI dz

a

eI N

4

nn

6

a

a

eI N

3

5

f rxx

ðkÞ b

a

eI N

T

NI

ðkÞ

zt

z

k¼1

a

NI

4

Z

Nm X

a

NI

 a a In virtue of Eq. (2.20) the displacements U I ; U I are related to  a a n s U I ; V I by

n

ð2:20Þ

ð2:25Þ

s

a

a

a

nn

ns

nz

Secondary variables : N I ; N I ; N I

ð2:26Þ

786

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

8 b " !#9 Np b M
2.13. Laminate constitutive equations

a

For the kth (orthotropic) mathematical layer we have the following 3-D stress–strain relationships

8 ðkÞ C 11 > > > > ðkÞ > > C 21 > > > > < C ðkÞ

8 ðkÞ rxx 9 > > > > > > > > > ryy > > > > > > > < =

rzz 31 ¼ > > r yz > > > > 0 > > > > > > > > > > > rzx > > > > 0 > > > : > ; > : ðkÞ rxy

ðkÞ

C 13

ðkÞ

0

0

C 22

ðkÞ

C 23

ðkÞ

0

0

ðkÞ C 32

C 33

ðkÞ

0

0

0

C 45 0

0

0

ðkÞ C 44 ðkÞ C 54

ðkÞ C 62

ðkÞ C 63

0

0

C 61

ðkÞ 9ðkÞ 8 9ðkÞ C 16 > > exx > > > > > > > > ðkÞ > > > C 26 > > eyy > > > > > > > > > > > ðkÞ > < = = ezz > C

C 12

36

> > cyz > > > 0 > > > > > > >c > > > > > zx > > > > 0 > > > > > :c > ; > ; ðkÞ xy C 66

ðkÞ

ðkÞ C 55

2 b !3 Np M ba ba @W b ba ba @W b X X J J b b e I ¼ 4 5 N A JI V J þ A JI þ A JI U J þ A JI i¼4;5 i4 i4 @y i5 i5 @x b¼1 J¼1 a

ð2:27Þ ðkÞ

where C ij are the transformed elastic coefficients in the global (laminate) co-ordinates, which are related to the elastic coefficients ðkÞ in the local (layer) co-ordinates C ij by the following relationship

C ¼ TCT T

IJ

¼

i¼1;2;6

ð2:28Þ

where T is the transformation tensor with components

2 6 6 6 6 6 T ij ¼ 6 6 6 6 4

2

cos2 hðkÞ 2 sin hðkÞ

sin hðkÞ cos2 hðkÞ

0 0

0

0

1

0

0

0

0 ðkÞ

sinh cosh

sin2hðkÞ sin2hðkÞ

0

0

0

coshðkÞ

sinhðkÞ

0

0 sinhðkÞ

0

ðkÞ

0 0

0 0

ðkÞ

sinh cosh

ðkÞ

0

coshðkÞ

0

0

cos2 hðkÞ  sin hðkÞ

0

2

ab

A IJ ¼

Z

k¼1

ab

Nm Z X

ij

ðkÞ

ðkÞ

zt

ij

k¼1

ðkÞ

Z

ðkÞ

zt ðkÞ

zb

b

ðkÞ

C ij

zb

k¼1

ab

Nm X

C ij F aI F bJ dz A IJ ¼

ðkÞ zb

ij

A IJ ¼

ðkÞ

zt

b dF J ðkÞ C ij F aI

a Nm abc X dF I dF J dz B IJK ¼ dz dz ij k¼1

dz

Z

c

B IJK

3 7 7 7 7 7 7 7 7 7 5

ð2:29Þ

Sometimes it is useful to express the equations of motion (2.23) in   terms of generalized displacements U aI ; V aI ; W aI . This can be achieved by substituting the expressions for the force and moment resultants from Eqs. (2.32) in to Eq. (2.23).

The laminate stiffnesses are defined as follows: Nm X

abc

c

abc @V abc @U K c K þ B IJK þ B IJK W K : @x @y i1 i2 i3 c¼1 K¼1 2 c c  Np Np  c  c c M abc @U c X X X @V K 1 abcd @W K @W L K 4 þ B IJK þ þ D IJKL @y @x @x @x 2 i6 i1 c;d¼1 K¼1 L¼1 c c c c a b c d a b c d 1 @W K @W L @W K @W L þ D IJKL ð2:31eÞ þ D IJKL @y @y @x @y 2 i2 i6

ab

Q

8 c Np M
ð2:31dÞ

dz

ðkÞ

zt ðkÞ

zb

C ij F aI F bJ F cK dz ðkÞ

3. Finite element formulation 3.1. In-plane discretization

abc

B IJK ¼ ij

Nm Z X k¼1

¼

ðkÞ

zb

Nm Z X k¼1

ðkÞ

zt

c

C ij F aI F bJ ðkÞ

abcd dF K dz D IJKL dz ij

ðkÞ

zt ðkÞ

zb

C ij F aI F bJ F cK F dL dz ðkÞ

ð2:30Þ

The laminate constitutive equations can be written as

N

a

I

i¼1;2;6

8 b " Np M
ab @V b ab @U bJ J þ A IJ þ A IJ W bJ : J¼1 i1 @x i2 @y i3 b¼1 !#) b b ab @U @V J J þA IJ þ @y @x i6 2 b c N N b b p p X M X X 1 abc @W J @W cK 1 abc @W J @W cK 4 B IJK þ þ B IJK 2 i1 @x @x 2 i2 @y @y b;c¼1 J¼1 K¼1 !# b c abc @W @W J K þB IJK ð2:31aÞ @x @y i6

¼

ab

A IJ

2 b !3 Np M ab ab @W b ab ab @W b X X J J b b 4 5 A IJ V J þ A IJ þ A IJ U J þ A IJ N I ¼ i¼4;5 i4 i4 @y i5 i5 @x b¼1 J¼1 a

ð2:31bÞ

The in-plane discretization is based on the multilevel mesh superposition approach [21], which is shown in Fig. 3. In the multilevel mesh superposition finite element analysis, the domain  is discretized into a collection of elements such that each inx plane model domain xa is discretized independently with its own mesh a

a  x

Ne [ e¼1

 ae x [

 a ¼ xa @ xa x [  ae ¼ xae x @ xae

ð3:1Þ ð3:2Þ ð3:3Þ

  We define the measure (area) of subdomain xae as k xae . In each subdomain xi there will be jMi j non overlapping meshes, where each level of mesh than the previous ones in the  is more refined  hierarchy, i. e. k xeaþ1 6 k xae . In general the boundaries of the models in Mi do not coincide. The highest level of mesh in xi is the one for model Si. This mesh will be referred to [36] as the top mesh, while the meshes of the remaining plate models within Mi are referred to as nontop meshes. The representative element domain xe of this mesh will be used to derive the weak forms of the equations for the theory, under the

787

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791 a

assumptions that each element has a known shape and associated interpolation functions.

model a and interface I and w e ðxÞ are two-dimensional interpolaim  tion functions in xae .

3.2. Weak forms

3.4. Order of a model

Multiplying the three equations in (2.23) by dU aI ; dV aI and dW aI , respectively, then integrating over the representative element domain xe and finally integrating per parts using the component forms of the Green–Gauss theorem (gradient or divergence theorem) to weaken the differentiability of U aI ; V aI and W aI results in the following equations

We will assume the models within M to be defined in a hierarchical fashion. In particular, we will assume that the properties of the models determining their hierarchal nature are the following:

3 Nbp M X ab @U b X @dU aI a @dU aI a e a a J5 a 4 dxdy 0¼ I IJ NIþ N I þ N I dU I  dU I 5 1 6 @x @y @t 2 xe b¼1 J¼1 I  a a  N I nx þ N I ny ds Z

2

@ xe

1

6

ð3:4aÞ 3 N bp a a a a M X ab @V b X a @dV @dV J a a I I e I dV  dV 4 5dxdy 0¼ I IJ NIþ N I þN I I 2 6 4 @y @x @t 2 xe b¼1 J¼1 I  a a N I nx þ N I ny ds  Z

2

@ xe

6

2

ð3:4bÞ Z



a

@dW I a @dW I a e a 0¼ NIþ N I þ N I dW aI  qb dW aI  qt dW aI dxdy 5 4 3 @y @x xe 8 2 b Np Z


m¼0

im

im



L1I 6    6 LaI 6    6 LM I       k x1e P    P k xae P    P k xM e L1e 6    6 Lae 6    6 LM e

ð3:6Þ

The order of a model a within the hierarchy M is identified by       parameters k Z aI ; LaI ; k xae ; Lae . A higher order model can be added to a sequence of models by creating a new model where any one of the above parameters is increased with respect to greatest element in the set M. 3.5. Finite element model Substituting Eq. (3.5) in to the weak forms in Eqs. (3.4) we obtain the semidiscrete finite element model of the hierarchical multiple plate models theory. e

e

e

e

e

IJ

I

IJ

I

I

€ a¼Fa K ab U a þM ab U

a

      k Z 1I P    P k Z aI P    P k Z M I

ð3:5Þ

where i = 1, . . . , 3 and Lae þ 1 is the number of nodes per 2-D element domain xae used to approximate the generalized displacements for

ð3:7Þ

e

where the matrix K ab is the direct element stiffness matrix, the IJ

e

e

matrix M ab is the element mass matrix, the column vector F a is I e

IJ

the element force vector and the column vector U a ¼ I

e  e e U aI; U aI; . . . ; U a I is the vector of the nodal degrees of freedom 1 2

e Lae þ1 e e e per element xe, where U aI ¼ U aI; U aI; U aIg ¼ j ij ij ij

e  P  e e a U aI; V aI; W aI ; I ¼ 1; . . . ; a2Mi N p ; i ¼ 1; . . . ; 3 and j ¼ 1; . . . ; j

j

j

Lae þ 1. Eq. (3.7) has 2ne unknowns, where

ne ¼ 3

X

!

N ap

 a  Le þ 1

ð3:8Þ

a2Mi

Because of the adopted modeling approach, both discretization and modeling errors can be reduced within a certain region of xi by selecting, in a hierarchical fashion, the effective number of models from Mi that are used within that region. Every time a new model a is superimposed on the previous models in certain region, the nodal degrees of freedom of that region are enriched by the degrees of freedom of new model. As consequence this will increase the total number of nodal degrees of freedom ne per element domain xe in that region, so other elements where less precision is required should not be enriched is order to keep the computational cost to the minimum required for the analysis. For example, in the case of the superposition of a FSDT with a linear LWT model displayed in Fig. 2, the set of available mathematical models in X is M ¼ f1; 2; 3g. In subdomain x2 the subset of available models is M2 ¼ f1; 2; 3g. During normal operation only models 1 and 2 are used, but when discretization and modeling errors need to be reduced then model 3 can be superimposed in x2 on models 1 and 2, at the expense of increasing the number of nodal degrees of freedom and consequently the computational efficiency. The direct element stiffness and element mass matrices can be expressed as follows:

788

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

2

ab

ab

ab

3

6 K IJ K IJ K IJ 7 6 11 12 13 7 7 6 ab ab ab 7 e 6 7 K ab ¼ 6 K IJ K IJ K IJ 7 6 IJ 22 23 7 6 21 7 6 ab a b a b 5 4 K IJ K IJ K IJ 31 32 33 3 2

ð3:9Þ

ab

6 M IJ 6 11 6 e 6 M ab ¼ 6 6 0 IJ 6 6 4 0

0 ab

IJ 22

0

07 7 7 7 07 7 7 7 ab 5 IJ 33

where a; b 2 Mi . The direct stiffness matrix is symmetric only in the linear case, while the mass matrix is always symmetric. Substituting the laminate constitutive Eqs. (2.32) in to the weak forms (3.4) we can express the element stiffness matrix and element mass matrix in terms of the laminate stiffnesses (2.30) as follows: ab

ab

ab

ab

abc

abc

abcd

!

K IJ ¼ f A IJ ; A IJ ; A IJ ; B IJK ; B IJK ; D IJKL ij

ij

 ab ab M IJ ¼ f I IJ

ij

ij

ij

ij

Fig. 3. Multilevel mesh superposition finite element model.

ð3:10Þ

ij

ð3:11Þ

ii

where i,j = 1, 2, 3. e It’s important to note that the diagonal terms in K ab present IJ linear and nonlinear self coupling contributions plus nonlinear coupling contributions from all the other models. The non diagonal terms instead present linear and nonlinear coupling contributions between each two models. This is summarized as follows (dropping indices IJ): e

e

L

NL

e

e

LC

NLC

K eaa ¼ K aa þK aa þ

M X

b¼1

e

K ab NLC

b–a

ð3:12Þ

K eab ¼ K ab þK ab Moreover, the general form of the above contributions is independent of any two models being considered. After assembly, the semi-discretized finite element equation can be given in the following form

€ a ¼ Fa K IJab U aI þ M IJab U I I

ð3:13Þ

where a; b 2 Mi . The fully discretized finite element model is obtained using a time approximation, e.g. the Newmark scheme. 3.6. Multiple models compatibility and uniqueness conditions Certain degrees of freedom (DOF) should be suppressed to ensure compatibility of the total displacement field and uniqueness of the solution of the multiple assumed displacement fields. For this reason additional boundary conditions need to be imposed on top of the physical boundary conditions before the finite element Eq. (3.13) can be solved. The additional boundary conditions are generally referred to [24,36] as compatibility and uniqueness conditions, and are shown in Fig. 4. 3.6.1. Compatibility conditions The total displacement field should be continuous throughout the composite mesh, which is the requirement of compatibility or C0 continuity. For this reason, homogeneous essential boundary

Fig. 4. Compatibility and uniqueness conditions for the case of structured mesh superposition.

conditions must be enforced on incompatible regions, setting to zero the DOF of the plate models in those regions. Note that the compatibility conditions are not imposed on boundaries that coincide with physical boundaries, as displayed in Fig. 4. Here DOF need only to be suppressed in accordance with the physical boundary conditions. Given two neighboring subdomains xi and xj, if an assumed T displacement field a is in Mi , but not in Mj , then the region @ xi @ xj is defined to be an incompatible region. Homogeneous essential boundary conditions must be enforced on incompatible regions, setting to zero the DOF of model a in those regions. An example is displayed in Fig. 4 for the case of structured mesh superposition. 3.6.2. Uniqueness conditions The other issue is the singularity of the assembled stiffness matrix in Eq. (3.13), which is caused by redundant rigid body modes. In fact there may be more than one set of assumed displacement fields that can be summed to yield the total displacement field,

789

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

as shown in Fig. 4 for the case of structured mesh superposition. To guarantee uniqueness of the solution, homogeneous essential boundary conditions must be enforced in coinciding nodes of different mesh levels by setting to zero the redundant DOF. These should be selected between the DOF of the highest level mesh, as shown in Fig. 4. For the case of the superposition of a FSDT with a linear LWT model there are five redundant variables that need to be set to zero to permit a unique solution for the remaining variables. Fig. 5 shows the in-plane deformation u(x, y) = u1(x, y) along a transverse material line A  B for a 4 layer plate. Two variables need to be set to zero. In this case the particular pair of variables is arbitrary [5], and although the same deformation is achieved in Fig. 5a and b, the numerical values of the remaining five nonzero values are different. 3.7. Calculation of the coupling terms The coupling terms in Eq. (3.12) are the result of an integral whose integrands contain products of in-plane interpolation functions defined for different mesh levels. A mapping between the different co-ordinate systems of the overlaid meshes is necessary to perform the numerical integration. This can be quite complex when unstructured mesh superposition is used. In the case of structured mesh superposition displayed in Fig. 4, a top/non-top transformation matrix [36] can be calculated that allows expressing the interpolation functions of any mesh level a 2 Mi to those of mesh level Si. A specific form of this matrix exists for any pair of mesh levels considered. 3.8. Inclusion of delaminations The occurrence of delaminations is common in composite laminates [37,38], so it is important to incorporate the kinematics of single and multiple delaminations in to the general theory presented here. For the purpose of the following discussion it is assumed that the domain X is crossed by D discontinuities @ XdI , where I = 1, . . . , D, between the lamina, i. e. one or more delaminations. As shown in Fig. 6, the two resulting parts of the domain are called ^ dI is the unit vector along the normal to @ XdI pointing XIþ and XI . n into XIþ . The following relation holds for all discontinuities: S XI XIþ ¼ X; 8I ¼ 1; . . . ; D. Furthermore, it is assumed that the delaminated interface @ XdI has co-ordinate zdI . Recalling that x ¼ ðx; zÞ, then x 2 XIþ , if x 2 x   and z P zdI . Otherwise, if x 2 XI , then x 2 x and z < zdI .

Fig. 6. The domain X is crossed by a delamination @ XdI at co-ordinate zdI .

Discontinuous through the thickness displacement fields, such as delaminations, are incorporated in the theoretical framework by performing an extrinsic enrichment of one or more of the assumed displacement fields. This can be achieved when, for the set of nodes that belong to the domain of the discontinuity, the functions F aI are multiplied by discontinuous enrichment functions, such as the Heaviside unit step functions [30] defined as follows:



HI z 

zdI



( ¼

1 if z P zdI ; 0

if z < zdI :

Another possible choice is the Heaviside sign-functions [30]. The additional generalized displacements udI are physically meaningful and they represent jump discontinuities in the displacement components of the total displacement field across the delamination. A cohesive law between laminae is incorporated where appropriate. This approach is outlined in Section 3.9. For example, Fig. 7 shows the in-plane deformation u(x, y) = u1(x, y) along a transverse material line A  B for the case of a 4



(a) u31 = u35 = 0

ð3:14Þ

(b) u32 = u33 = 0

Fig. 5. Suppression of redundant in-plane deformations for the case of superposition of FSDT, u1 and u2, and linear LWT models, u3.

790

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791

T d ¼ LT T d L

ð3:18Þ

where the orthogonal transformation matrix L performs the transformation from the co-ordinate system aligned to the discontinuity to the global co-ordinate system. Therefore there is a new contribution dVczm to the virtual work energy given in Eq. (2.11) performed by cohesive forces amongst the layers, which has the following form

dV czm ¼ ¼

1

2

Fig. 7. In-plane deformation for the case of superposition of FSDT, u and u , linear LWT, u3, and delamination models, u4.

We assume that the tractions between layers can be modeled using a cohesive zone model (CZM) [39–43]. Note that only the presence, and not the particular form, of the CZM is assumed. Considering the discontinuities @ XdI as an internal boundary condition, the following relationships can be added to the boundary conditions of the problem

x 2 @ XdI

ð3:15Þ

where t dI are the tractions at the internal boundaries @ XdI . In the co-ordinate system aligned with the orientation of the discontinuity, which is displayed in Fig. 8, the tractions t dI are related to the displacement jumps v dI via tangent stiffness matrix of the traction-separation law [41]

t dI ¼ T d v dI

x 2 @ XdI

ð3:16Þ

Note that normally in the co-ordinate system aligned with the orientation of the discontinuity the z-direction is chosen to coincide with the direction normal to the discontinuity. In general, the displacement jumps v dI are related to the additional generalized displacements udI via the magnitude of the jump in the step function, which in the present case is equal to one.

v dI ¼ udI

x 2 @ XdI

ð3:17Þ

The tangent stiffness matrix in the global co-ordinate system is equal to

Fig. 8. Orientation configuration.

of

the

I¼1

@ XdI

I¼1

XdI

"Z D X



dudI

dxdy

#   T dxx U dI dU dI þ T dyy V dI dV dI þ T dzz W dI dW dI dxdy

d

d

NI

NI

1

T ¼

2

Nm Z X k¼1

z

ðkÞ

zt ðkÞ b

f rxx

n

ryy gT HI dz þ T dxx T dyy

oT ð3:19Þ

d

e I ¼ Td N zz 3

3.9. Cohesive zone model

^ dI r ¼ t dI n

# t dI

If tractions between layers are modeled using a CZM, then the force resultants in Eqs. (2.16) can be modified as follows for the delamination field

layer plate and the superposition of FSDT, linear LWT and delamination models. A single delamination (field u4) occurs at z = 0.

"Z D X

discontinuity

within

the

laminate

deformed

The terms introduced by the cohesive forces in Eq. (3.19) will have to be taken into account on top of the contributions given in (3.11) when computing the element stiffness and mass matrices for the delamination field. 4. Conclusions A hierarchical multiple plate models approach has been proposed by which any number and type of plate models may be superimposed in order to perform a global-local analysis of laminated composite structures. This is particularly useful for thick laminates when localized phenomena such as impact loading or damage may arise. First the general theory and then the finite element formulation based on a multilevel mesh superposition method was derived. The capability of the theory to model delaminations was also shown. A cohesive zone model was included in the development. The theory presented in this paper can be extended to incorporate thermal effects, shell structures, mixed variational formulations and multifield problems. The theoretical framework will be validated against numerical results in a follow up paper by the authors. References [1] Rohwer K, Friedrichs S, Wehmeyer C. Analyzing laminated structures from fibre-reinforced composite material—an assessment. Tech Mech 2005;25(1):59–79. [2] Mohite PM, Upadhyay CS. Region-by-region modeling of laminated composite plates. Comput Struct 2007;85:1808–27. doi:10.1016/j.compstruc.2007.04.005. [3] Carrera E. Evaluation of layerwise mixed theories for laminated plates analysis. AIAA J 1998;36(5):830–9. [4] Reddy JN, Robbins Jr DH. Theories and computational models for composite laminates. Appl Mech Rev 1994;47(6):147–69. doi:10.1115/1.3111076. [5] Reddy J. Mechanics of laminated composite plates and shells: theory and analysis. CRC; 2004. [6] Reddy JN. A generalization of two-dimensional theories of laminated composite plates. Commun Appl Numer Methods 1987;3:173–80. [7] Barbero EJ, Reddy JN. Nonlinear analysis of composite laminates using a generalized laminated plate theory. AIAA J 1990;28:1987–94. doi:10.2514/ 3.10509. [8] Barbero EJ, Reddy JN. Modeling of delamination in composite laminates using a layer-wise plate theory. Int J Solids Struct 1991;28(3):373–88. doi:10.1016/ 0020-7683(91)90200-Y. [9] Robbins Jr DH, Reddy JN. Variable kinematic modelling of laminated composite plates. Int J Numer Methods Eng 1996;39:2283–317. [10] Demasi L. Refined multilayered plate elements based on Murakami zig-zag functions. Compos Struct 2005;70(3):308–16. doi:10.1016/ j.compstruct.2004.08.036.

S.L. Angioni et al. / Composite Structures 93 (2011) 780–791 [11] Carrera E, Ciuffreda A. A unified formulation to assess theories of multilayered plates for various bending problems. Compos Struct 2005;69:271–93. doi:10.1016/j.compstruct.2004.07.003. [12] Carrera E, Brischetto S, Robaldo A. Variable kinematic model for the analysis of functionally graded material plates. AIAA J 2008;46(1):194–203. doi:10.2514/ 1.32490. [13] Barbero EJ. On a generalized laminated plate theory with application to bending vibration and delamination buckling, Ph.D. thesis, Virginia Polytechnic Institute and State University, Blacksburg, VA; October 1989. [14] Robbins Jr DH, Reddy JN. Adaptive hierarchical kinematics in modeling progressive damage and global failure in fiber-reinforces composite laminates. J Compos Mater 2008;42(2):143–72. doi:10.1177/ 0021998307086210. [15] Williams TO. A generalized multilength scale nonlinear composite plate theory with delamination. Int J Solids Struct 1999;36:3015–50. doi:10.1016/S00207683(98)00138-3. [16] Mourad HM, Williams TO, Addessio FL. Finite element analysis of inelastic laminated plates using a global-local formulation with delamination. Comput Methods Appl Mech Eng 2008;198:542–54. doi:10.1016/j.cma.2008.09.006. [17] Kanouté P, Boso D, Chaboche J, Schrefler B. Multiscale methods for composites: a review. Arch Comput Methods Eng 2009;16(1):31–75. [18] Cho J, Oden J. Priori error estimations of hp-finite element approximations for hierarchical models of plate- and shell-like structures. Comput Methods Appl Mech Eng 1996;132(1):135–77. [19] Actis RL, Szabo BA, Schwab C. Hierarchic models for laminated plates and shells. Comput Methods Appl Mech Eng 1999;172(1–4):79–107. doi:10.1016/ S0045-7825(98)00226-6. [20] Zienkiewicz O, Taylor R. The finite element method for solid and structural mechanics. Butterworth-Heinemann; 2005. [21] Fish J. The s-version of the finite element method. Comput Struct 1992;43(3):539–47. [22] Fish J, Markolefas S. Adaptive s-method in linear elastostatics, in: 33rd AIAA/ ASME/ASCE/AHS/ASC structures, structural dynamics and materials conference, Dallas, TX; 1992. p. 313–6. [23] Fish J, Markolefas S, Guttal R, Nayak P. On adaptive multilevel superposition of finite element meshes for linear elastostatics. Appl Numer Math 1994;14(1– 3):135–64. [24] Park J, Hwang J, Kim Y. Efficient finite element analysis using mesh superposition technique. Finite Elem Anal Des 2003;39(7):619–38. [25] Fish J. Hierarchical modelling of discontinuous fields. Commun Appl Numer Methods 1992;8(7):443–53. [26] Fish J, Markolefas S. The s-version of the finite element method for multilayer laminates. Int J Numer Methods Eng 1992;33(5):1081–105.

791

[27] Fish J, Guttal R. The s-version of finite element method for laminated composites. Int J Numer Methods Eng 1996;39(21):3641–62. [28] Belytschko T, Gracie R, Ventura G, Ogata S, Umeno Y, Kohyama M, et al. A review of extended/generalized finite element methods for material modeling. In: Proceedings of the 7th World congress of computational mechanics on computational bridging of length scales in complex materials, vol. 20; 2006. p. 21. [29] Melenk J, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996;139(1–4):289–314. [30] Mohammadi S. Extended finite element method. Oxford: Blackwell Publishing; 2008. [31] Lee S, Song J, Yoon Y, Zi G, Belytschko T. Combined extended and superimposed finite element method for cracks. Int J Numer Methods Eng 2004;59(8):1119–36. [32] Crisfield M, Alfano G. Adaptive hierarchical enrichment for delamination fracture using a decohesive zone model. Int J Numer Methods Eng 2002;54(9):1369–90. [33] Schwab C. Hierarchic modelling in mechanics in wavelets. Multilevel methods and elliptic PDEs. Oxford University Press; 1997. [34] Dauge M, Faou E, Yosibash Z. Plates and shells: asymptotic expansions and hierarchical models. Encyclopedia of computational mechanics, vol. 1; 2004. p. 199–236. [35] Ciarlet P. A justification of the von Karman equations. Arch Ration Mech Anal 1980;73(4):349–89. [36] Yue Z, Robbins D. Adaptive superposition of finite element meshes in nonlinear transient solid mechanics problems. Int J Numer Methods Eng 2007;72(9):1063–94. [37] Meo M, Thieulot E. Delamination modelling in a double cantilever beam. Compos Struct 2005;71(3–4):429–34. [38] Barbieri E, Meo M. A meshfree penalty-based approach to delamination in composites. Compos Sci Technol 2009;69(13):2169–77. [39] Needleman A. A continuum model for void nucleation by inclusion debonding. J Appl Mech 1987;54:525. [40] Tvergaard V. Effect of fibre debonding in a whisker-reinforced metal. Mater Sci Eng: A 1990;125(2):203–13. doi:10.1016/0921-5093(90)90170-8. [41] Camanho P, Davila C, De Moura M. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003;37(16):1415. [42] de Borst R, Gutiérrez M, Wells G, Remmers J, Askes H. Cohesive-zone models, higher-order continuum theories and reliability methods for computational failure analysis. Int J Numer Methods Eng 2004;60(1):289–315. [43] Xie D, Waas A. Discrete cohesive zone model for mixed-mode fracture using finite element analysis. Eng Fract Mech 2006;73(13):1783–96.