A geometrically exact formulation for thin multi-layered laminated composite plates: Theory and experiment

Composite Structures 93 (2011) 1649–1663

Contents lists available at ScienceDirect

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

A geometrically exact formulation for thin multi-layered laminated composite plates: Theory and experiment Walter Lacarbonara ⇑, Michele Pasquali Dipartimento di Ingegneria Strutturale e Geotecnica, Università di Roma La Sapienza, Via Eudossiana 18, Rome 00184, Italy

a r t i c l e

i n f o

Article history: Available online 12 January 2011 Keywords: Geometrically exact plate Nonlinear curvatures Nonlinear membrane stretching Multi-layered composite plates Experimental equilibrium paths

a b s t r a c t A geometrically exact approach is employed to formulate the equations of motion of thin multi-layered isotropic and laminated composite plates subject to excitations that cause large strains, displacements, and rotations. The linearization of the obtained semi-intrinsic theory leads to the Mindlin–Reissner theory while an ad hoc truncated kinematic approximation delivers, as a by-product, the Föppl–von Kármán theory of plates. An experimental validation is sought for fully clamped plates which are either of the isotropic single-layered type or of the multi-layered laminated composite type. To this end, nonlinear equilibrium paths are constructed both theoretically and experimentally when the plates are subject to a quasi-statically increasing central point load. The comparisons between the experimentally obtained results and those furnished by the geometrically exact theory as well as by the Föppl–von Kármán (FVK) theory show the high accuracy of the proposed nonlinear theory while the FVK theory becomes increasingly inaccurate at deflection amplitudes of the order of the plates thickness. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Composite laminated multi-layered plates with sufficiently small thickness-to-span ratios are often prone to deformations that involve large displacements, rotations, and, possibly, large strains. This geometric nonlinear phenomenology arises in various operating conditions, such as in the post-buckling or post-flutter range as well as when the plates are subject to time-varying excitations near resonances. The leading difference between single-layer plates and multilayer plates with thickness-to-span ratios approximately above 1/ 10 consists of the well-known zig-zag effect arising from the fact that different tangential elastic compliances of the layers cause the in-plane displacement components to exhibit a rapid change of their slopes in the thickness direction at the layer interfaces. In general, the in-plane stresses can be discontinuous at each layer interface while the transverse stresses are continuous. Multi-layered plates exhibit even more complex effects such as a high in-plane anisotropy and a high transverse flexibility when the laminates are made of anisotropic layers (in particular, orthotropic). Inspired by Koiter’s well-known remark [1,2], this instance has motivated the development of higher-order shear deformation models which describe, through refined kinematic representations, both thickness-wise and transverse shear deformations [3–5]. In ⇑ Corresponding author. Tel.: +39 64458 5293; fax: +39 64884 852. E-mail address: [email protected] (W. Lacarbonara). 0263-8223/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2010.12.005

addition, in the context of more accurate formulations, the continuity in the interlaminar transverse shear and normal stresses is also enforced [6,7]. Most of the higher-order shear deformation theories are linear both from a geometric and a constitutive point of view with the exception of the geometric stiffness terms associated with destabilizing forces. A large body of literature has addressed the task of constructing fully nonlinear two-dimensional reduced theories for isotropic plates while a less substantial literature has focused on multi-layered composite plates either in the context of intrinsic or semiintrinsic formulations (the Cosserat theory of plates, see [10–14]) or in the context of models which are fully deduced from the nonlinear elasticity theory. Podio-Guidugli [15] obtained the thin plate theory, based on the Kirchhoff–Love assumptions [16], from the three-dimensional linear theory for homogeneous, linearly elastic transversely isotropic, constrained materials. Similar attempts have been undertaken for the thin plate theory of Mindlin–Reissner (MR) [17] and its variants, also in the context of composite plates [18,19]. The Föppl–von Kármán (FVK) [20,21] plate equations are one of the popular nonlinear reduced descriptions of plates dynamics. The FVK model expresses the elastic energy of a deformed elastic plate as the sum of stretching (the stretches are expressed in terms of the first deflection gradient up to quadratic order) and bending energies (the curvatures are linear in the second-order gradient of the deflection). Other 2D reduced theories usually bear the same structure, i.e. their energy is given by the sum of

1650

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

a stretching term and a bending term. The FVK theory has been shown to be derivable from the three-dimensional elasticity theory by means of an ad hoc asymptotic expansion [22]. Due to its relative simplicity, the FVK theory of plates has been widely used in the study of static (post-buckling) and dynamic (post-flutter) problems. For example, the FVK thin plate theory was employed in [23] to predict the nonlinear behavior of plates subject to compressive in-plane loads via a semi-analytical approach. The pseudo-arclength path following allowed to trace the nonlinear equilibrium paths with limit points of various types when the considered rectangular plates were subjected to combined outof-plane and in-plane loads. The tremendous amount of literature based on the FVK plate equations is justified by the theoretical and practical need to compute the post-critical states/motions of plate structures not far from the static/dynamic bifurcation points where the kinematic and dynamic approximations inherent in the theory do not introduce in general appreciable errors. The physical soundness as well as the accuracy of the FVK plate equations have been often questioned. The theory is in fact based on: (i) an approximate geometry of deformation, (ii) assumptions about the way the stress varies across the thickness, (iii) linear constitutive relationships, (iv) neglect of some components of strain (normal and transverse shear strains), and (v) an apparent confusion of the referential and spatial descriptions. To overcome these limitations, Antman [24] proposed a nonlinearly elastic Cosserat (geometrically exact) plate theory to study the axisymmetric buckling of isotropic single-layered plates which can suffer thickness changes, as well as flexure, mid-plane extension, and shear. A variety of approximate nonlinear theories of plates have been proposed in the literature, often with assumptions of small strains but possibly large deflections and rotations [25,27,26]. In particular, in [26] both geometric and material nonlinearities exhibited by paper were accounted for in the attempt to describe its highly nonlinear responses. Only recently, a few works have proposed geometrically exact intrinsic theories for the dynamics of composite plates, seldom offering assessments of their accuracy through extensive numerical computations and experimental validations. One of the few exceptions is [28]. In particular, Yu et al. have proposed a variational-asymptotic approach [18] to decouple systematically the original three-dimensional problem into a linear onedimensional through-the-thickness analysis and a nonlinear twodimensional plate analysis, resulting into an equivalent singlelayer Reissner–Mindlin theory, obtaining an excellent accuracy (see [29]) which is comparable to that of higher-order, layer-wise theories. Extensive research efforts have been directed towards geometrically nonlinear finite element formulations of laminated composite plates and shells as [30–32], to cite but only a few. This work presents a geometrically exact theory of thin multilayered composite plates with general stacking sequences which accounts for mid-plane stretching, flexure, and transverse shear strains. The predictions of the theory, by considering fully clamped plates subject to a central point load, are compared with those of the FVK theory in the case of various multi-layered composite plates and isotropic single-layered plates. The theoretical predictions are contrasted with the experimental measurements in various points of the plates by means of a carefully designed experimental setup. An effective nonlinear geometrically exact model of composite plates is useful in various applications, including recent applications with carbon-nanotube-reinforced composite plates [33] and physics-based nonlinear system identification procedures of composite laminated plates [34]. As shown in [34] a fully nonlinear refined modeling approach delivers very accurate information through the dynamic response that can be unfolded and exploited for nonlinear identification purposes.

2. A geometrically exact formulation In the next sections we present the kinematic, dynamic, and constitutive formulation for thin multi-layered composite plates. To start with, the geometrically exact kinematics of the plate are presented. The kinematics are based on the fundamental assumption that transverse material fibers remain rigid (hence, they do not suffer stretching or bending) no matter what the loading conditions are. The assumed compatible nonlinear kinematic description clearly leads to a local violation of the interlaminar continuity of transverse stresses. However, in a large array of works, among which [8,9], it has been ascertained that this local violation does not introduce significant errors in the global response characteristics such as the elastic deflections, eigenfrequencies, and so on. 2.1. Notation We employ Gibbs notation for vectors and tensors. Vectors, which are elements of Euclidean 3-space E3 , and vector-valued functions are denoted by lower-case, italic, bold-face symbols. The dot product and cross product of (vectors) u and v are expressed as u  v and u  v, respectively. The value of tensor A at vector u is denoted by A  u. The transpose of tensor A is denoted by AT. A dyad formed by vectors u and v is expressed as uv instead of the more traditional notation u  v. The partial derivative of a function f with respect to the scalar argument s is denoted by either fs or @ sf. The time derivative of vector u is denoted by either _ Throughout the manuscript, (O, e1, e2, e3) indicates a fixed ut or u. (inertial) Cartesian frame with origin O. 2.2. Kinematics of thin plates The geometry of the undeformed (stress-free) configuration is described by the position vector of the material points of transvero o sal fibers collinear with unit vector b3 with b3  e3 . We let o r (xM) = x1e1 + x2e2 be the position vector of the material points of a base plane Bo of the reference configuration (see Fig. 1). The subscript M, here and henceforth, takes values 1 and 2, unless otherwise stated, and when it is repeated, the summation in its range must be considered. The position vector of the material points of the fiber through ro(xM) is then represented by o rðxM ; zÞ ¼ r o ðxM Þ þ zb3 where z  x3 denotes the position along the fiber and z 2 [z1, z2]. The superscript o characterizes all quantities referred to the base plane Bo . o o o The basis ðb1 ; b2 ; b3 Þ is a local basis employed to identify material fibers collinear with the coordinate lines (i.e., x1 = const, o x2 = const) in the base plane Bo (i.e., bM = eM) as well as those collinear with the thickness direction. The reference configuration of the plate is thus

 o B ¼ rðxM ; zÞ ¼ r o ðxM Þ þ zb3 ; ro ¼ x1 e1 þ x2 e2 ; x1 2 ½0; a; x2 2 ½0; b; z 2 ½z1 ; z2 g:

ð1Þ

The boundary of the reference configuration is denoted by @B and is assumed to be Lipschitz-continuous. It comprises the edges, expressed as @Bol ¼ @Bo  ½z1 ; z2  (@Bo is the boundary of the base

Fig. 1. Geometry of the stress-free plate (configuration Bo ) and deformed material lines through ro(x1, x2).

1651

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

plane), the lower and top boundary planes, denoted @Bo1 and @Bo2 , respectively. The internal kinematic constraint enforcing the rigidity of the transverse fibers implies that the plate nontrivial strains are five; for example, two strains denote the stretching of fibers collinear with the axes parallel to the base plane, one strain is the shear strain between in-plane fibers, while the other two strains are the shear strains between transverse fibers and in-plane fibers. We let po be the position vector in the actual configuration of the material point that occupies position ro in the reference configuration; we further denote by b3 the unit vector collinear at time t to the thickness-wise fiber, passing through ro in Bo . The current configuration of the plate is thus described by

 ¼ fpðxM ; z; tÞ ¼ po ðxM ; tÞ þ zb3 ðxM ; tÞ; B x1 2 ½0; a; x2 2 ½0; b; z 2 ½z1 ; z2 g:

ð2Þ

o

By letting u (xM, t) represent the displacement vector of the material point ro, we can express the actual position vector of the material points of the base plane as po(xM, t) := ro (xM) + uo(xM, t). It is convenient to introduce the unit vectors (b1, b2) which, together with b3 = b1  b2, make a suitable local basis for the current configuration. A key step of the kinematic formulation is the way the current orientation of the thickness-wise fibers is parametrized. The transverse fibers suffer finite rotations about a rotation axis orthogonal to the reference orientation of the fibers which implies that this o o axis is parallel to the base plane ðb1 ; b2 Þ. There are various possible parametrizations, for example, based on two Euler angles, Rodrigues parameters, or Quaternions. Since we are interested in parametrizing the full basis (b1, b2, b3) with the chosen angles, a straightforward approach is to conceive a sequence of two succeso o o sive rotations of the initial triad ðb1 ; b2 ; b3 Þ. The first rotation is, for o example, taken about b1  e1 of angle /2, expressed by the orthogð1Þ ð1Þ ð1Þ onal tensor R(1). The rotated basis ðb1 ; b2 ; b3 Þ is then subject to a ð1Þ rotation about b2 of angle /1, described by the orthogonal tensor R(2). The final unit vectors (b1, b2, b3) are given by o

bk ¼ R  bk ;

R ¼ Rð2Þ  Rð1Þ :

ð3Þ

The generalized strains can be calculated after the deformation gradient is determined according to T

o

T

F ¼ ð$pÞ ¼ ½ðe1 @ 1 þ e2 @ 2 þ e3 @ z Þðp þ zb3 Þ

¼ ð@ 1 po Þe1 þ ð@ 2 po Þe2 þ b3 e3 þ zð@ 1 b3 Þe1 þ zð@ 2 b3 Þe2 ;

(a)

where the notation @ M() := @()/@xM has been introduced for sake of notational simplicity. For admissible deformations we require that detF > 0. We let

moM :¼ @ M po

ð5Þ

denote the stretch vectors associated with the material fibers o through ro and collinear with bM (see Fig. 2). The stretch vectors o o o can be expressed as m M = @ Mp = @ Mro + @ Muo = bM + @ Muo. Moreover, it is straightforward to prove that

@ M b3 ¼ lM  b3 ;

ð6Þ

where (l1, l2) are the curvature vectors associated with the rotations of the transverse fibers. The curvature vector lM is the axial vector of the curvature tensor

CM ¼ ð@ M RÞ  RT :

ð7Þ

Consequently, by adopting the following componential form of the curvature vectors in the fixed basis (e1, e2, e3), lM = j1Me1 + j2Me2 + j3Me3, (7) yields j1M = CM(2,3), j2M = CM(1,3), j3M = CM(1, 2). These considerations allow to rewrite the deformation gradient in a more compact form as FðxM ; z; tÞ ¼ ðm o1 þ zl1  b3 Þ e1 þ ðm o2 þ zl2  b3 Þ e2 þ b3 e3 . To calculate the stretches of material fibers through r(xM, z) lying in a plane parallel to the base plane at a distance z from it, the above given deformation gradient F(xM, z, t) is o applied on the unit vector bM and obtain o

F  bM ¼ m oM þ zlM  b3 ¼: m M :

Hence, the vectors mM =: m + zlM  b3 represent the stretches of material fibers collinear with bMo lying at a distance z from the base plane Bo . These stretches thus vary linearly with z. As per kinematic assumption, transverse fibers do not suffer o stretching; in fact, m ¼ F  b3 ¼ b3 which implies jmj = jb3j = 1. The deformation gradient may thus be written in its most compact form as F = m1e1 + m2e2 + b3e3. By using this expression of the deformation gradient, we can calculate the Cauchy–Green strain tensor C – or the Green–Lagrange strain tensor E = (C  I)/2 – as

C ¼ FT  F ¼ m21 e1 e1 þ m22 e2 e2 þ ðm 1  m 2 Þðe1 e2 þ e2 e1 Þ þ ðm 1  b3 Þðe1 e3 þ e3 e1 Þ þ ðm 2  b3 Þðe2 e3 þ e3 e2 Þ þ e3 e3 :

ð4Þ

ð8Þ

o M

ð9Þ

Hence, the strain tensor C, by recalling m M ¼ m oM + z lM  b3, depends on four vectors, namely, ðm o1 ; m o2 ; l1 ; l2 Þ, which leads to 12

(b)

Fig. 2. (a) Undeformed and deformed material lines through ro(x1, x2) in the base plane Bo . The shaded areas denote infinitesimal material areas surrounded by the material lines through ro. (b) View of the stretch vectors (deformed material lines through ro(x1, x2) of the base plane) in the plane (b1, b2).

1652

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

generalized strain components. On the other hand, the independent kinematic parameters are five; for example, the position vector of the base plane po and the rotations (/1, /2) or the displacement vector uo and the rotations (/1, /2). The componential representation of m oM in the actual basis (b1, b2, b3) is given by

m oM ¼ mo1M b1 þ mo2M b2 þ mo3M b3 :

ð10Þ o

Therefore, the (actual) stretches of the fibers collinear with b1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi o and b2 are jm oM j ¼ m oM  m oM ; ðM ¼ 1; 2Þ while we let

mM :¼ m oM  bM

ð11Þ

be the (generalized) stretch of the fiber of the base plane collinear o with bM . According to (8), the stretches are linear functions of z through the multiplicative components of lM  b3. We thus conveniently let

l^ M :¼ lM  b3 ¼ l^ 1M b1 þ l^ 2M b2 ;

ð12Þ

^ 1M b1 þ l ^ 2M b2 Þ. Thus, the stretches of fibers iniso that m M ¼ m þ zðl tially collinear with e1 and e2 are given, respectively, by m1 ¼ m1 þ zl^ 11 ; m2 ¼ m2 þ zl^ 22 . If we let o M

lM ¼ l2M b1 þ l1M b2 þ l3M b3 ;

ð13Þ

^ 11 ¼ b2  l1 ; ^ M ¼ lM  b3 ¼ l1M b1  l2M b2 which entails l then l l^ 21 ¼ b1  l1 ; l^ 12 ¼ b2  l2 ; l^ 22 ¼ b1  l2 . For notational simplicity, we further let

l 1 :¼ l^ 11 ; l 2 :¼ l^ 22 ) l 1 ¼ b2  l1 ; l 2 ¼ b1  l2 ;

ð14Þ

moM  m 3 m oM  b3 m M  m 3 m M  b3 ¼ ; and sin c3M ¼ ¼ jm oM j jm oM j jm M j jm M j

ð15Þ

where, in force of the transverse inextensiblity, it is m3 = b3. The transverse shear strains turn out to be simply expressed as

sin c3M ¼

 12 is referred to as twist curvature. To make some later statewhere l ments more concise, henceforth, the algebraic vector u, collecting all the strain variables, will be used along with its rate u_ which lists the rates of the corresponding generalized strains. Besides the compatibility equations that relate the generalized strains, there is another set of compatibility equations, namely @ t(lM  bk) = @ M(x  bk) where the angular velocity vector x stems from @ tbk = x  bk. The kinematic boundary conditions can be prescribed by fixing the position and orientation of the fibers on the plate edges @Bo  ½z1 ; z2 . Let s be the arclength along the boundary of the base plane @Bo and let po(s, t) be the position vector of material points on @Bo in the reference configuration. The kinematic boundary conditions for shear deformable plates are expressed as

 3 ðs; tÞ; b3 ðs; tÞ ¼ b

s 2 @Bo ;

t 2 ½0; 1Þ;

ð21Þ

 3 ðs; tÞÞ represent the kinematic data. The initial conwhere ðp ðs; tÞ; b ditions are

 1; l  2 Þ are referred to as flexural curvatures. The scalar quantities ðl We next pause to discuss the shear strains of the plate. We express the transverse shear strains calculated at the base plane (z = 0) and those at a distance z from it as

g M ; jm oM j

ð20Þ

o

m1 ¼ m1 þ zl 1 ; m2 ¼ m2 þ zl 2 :

sin co3M ¼

l 12 :¼ l12  l21 ;

 o ðs; tÞ; po ðs; tÞ ¼ p

so that the stretches are accordingly expressed as

sin co3M ¼

On the other hand, the two angles that b3 makes with m o1 and m o2 (more precisely, the change in angles with respect to the initially orthogonal angles) measure the transverse shear strains. Since the curvature vectors l1 and l2 have altogether six components, it appears that the overall number of generalized plate strains is 11. However, it can be shown that, out of the six curvature components, they can be conveniently reduced to three independent curvatures, two flexural curvatures and one twist curvature (see Appendix B). Hence, the plate strain tensor can be effectively expressed in terms of eight generalized strains, out of which only five strains are truly independent. The generalized nonlinear strain measures coalesce into the well-known linear strains for shear deformable plates if we let

g M ; jm M j

ð16Þ

where

g M :¼ m oM  b3 ¼ @ M po  b3

ð17Þ

are referred to as the (generalized) transverse shear strains. We next calculate the shear strain between fibers of the base o o plane, collinear with ðb1 ; b2 Þ, as sin co12 ¼ ðm o1  m o2 Þ=ðjm o1 jjm o2 jÞ. If we let

g 12 :¼ mo1  m o2 ¼ @ 1 po  @ 2 po ;

ð18Þ

 12 =ðjm o1 jjm o2 jÞ. g  12 is referred to as the (generalized) then sin co12 ¼ g membrane shear strain. The above discussion shows that there are only five independent components of

m o1 ¼ m1 b1 þ m21 b2 þ g 1 b3 ; m o2 ¼ m12 b1 þ m2 b2 þ g 2 b3

ð19Þ o 12

12 ¼ m and where, for sake of notational consistency, we have put m m21 ¼ mo21 . This is in line with the notion that the magnitude of the vectors m o1 and m o2 represents the membrane stretching of material fibers in the plane and the change in the subtended angle (with respect to the original orthogonal configuration) measures the membrane shear strain.

 o ðxM Þ; po ðxM ; 0Þ ¼ p

 3;o ðxM Þ; b3 ðxM ; 0Þ ¼ b

ro ¼ xM eM 2 Bo : ð22Þ

The initial conditions must be compatible with the above kinematic data. 2.3. Equations of motion We proceed with the derivation of the equations of motion within the context of the geometrically exact formulation. The mechanical data of the problem are simply represented by f(xM, t) and c(xM, t) (the force and moment resultants per unit reference area of the base plane Bo ) and by fo(s, t) and co(s, t) (the force and moment resultants per unit reference length s of the edge) where s is the arclength along the curve represented by the boundary of the base plane @Bo . We consider a section plane, whose normal is the unit vector a, that cuts the plate in the reference configuration into two parts. The unit vector a is taken to be positive in the outward direction with respect to the plate part of interest. We consider its deformed image shown in Fig. 3. On the undeformed cross-sectional plane let us consider a material point r(xM, z) which occupies position p(r, t)  the at time t in the deformed configuration; further, let us denote a unit normal to the deformed region at p. Let t a ðp; tÞ be the Cauchy stress vector that quantifies the contact force at p and time t per unit actual area of the deformed region. In view of reducing the stress distribution along the thickness (by taking the pole of reduction on the base plane Bo ), we introduce generalized stress resultants by integrating the first Piola– Kirchhoff (or nominal) stress vector in the thickness domain [z1, z2] and computing its moment about the base plane Bo . The first Piola–Kirchhoff stress vector ta(r, t) quantifies the contact force at p(r, t) per unit area of the undeformed plane through r in the reference configuration having unit normal a.

1653

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

@ 1 m1 þ @ 2 m2 þ m o1  n1 þ m o2  n2 þ c ¼ qIo b3  uott þ qJ o xot ;

ð27Þ

where c denotes the external couples per unit of area and the moRz ment of inertia qJo is qJ o :¼ z12 qz2 dz. The mechanical boundary conditions on the edge boundary @Bo are na(s, t) = fo(s, t), ma(s, t) = co(s, t). 2.4. Componential form of the equations of motion

Fig. 3. Part of the undeformed plate obtained through a section plane of normal a and its deformed configuration with the Cauchy stress vector t a ðr; tÞ acting on the deformed plane at p(r, t).

In the adopted referential description, the stress resultants are

na ðro ; tÞ :¼

Z

z2

ma ðr o ; tÞ :¼ b3  ma ; ma :¼

t a dz;

z1

Z

z2

zt a dz:

z1

ð23Þ If T(r, t) denotes the first Piola–Kirchhoff stress tensor, the (nominal) stress vector is ta(r, t) = T(r, t)a(r, t). Consequently, Rz Rz na ðr o ; tÞ ¼ z12 ðT  aÞ dz; ma ðr o ; tÞ ¼ z12 zðT  aÞ dz. Vectors na and ma represent the contact force and the contact moment per unit reference length in the direction orthogonal to the thickness direction. These generalized stress resultants per unit reference length can be suitably expressed in the local basis (b1, b2, b3) with b1 col . The resulting componential form of na and ma is linear with a na = N1ab1 + N2ab2 + Qab3 and ma = M2ab1 + M1ab2 with

N1a :¼ M 2a

R z2

R z2

b1  t a dz; N2a :¼ z1 b2  t a dz; Q a :¼ R z2 Rz :¼  z1 zb2  t a dz; M 1a :¼ z12 zb1  t a dz: z1

R z2 z1

b3  t a dz; ð24Þ

(N11, N21, Q1) are the tension, the in-plane shear force, and the transverse shear force per unit reference length in the direction orthogonal to the thickness direction (i.e., by putting a = 1) as shown in Fig. 4. On the other hand, M21 and M11 are the twisting moment and the bending moment, respectively, per unit reference length. The statement of the local balance of linear momentum delivers the first equation of motion as o

@ 1 n1 þ @ 2 n2 þ f ¼ qh uott þ qIo b3;tt ;

ð25Þ

where o

qh :¼

Z

z2

z1

o

q dz; qI :¼

Z

It turns out that the most convenient componential form of the equations of motion (also for the successive numerical implementation of the model) is obtained by projecting the balance of linear momentum in the fixed basis (e1, e2, e3) while expressing the balance of angular momentum in the intrinsic basis (b1, b2, b3). In the fixed basis (e1, e2, e3), the generalized stresses are na = n1ae1 + n2ae2 + n3ae3 while the intrinsic representation of ma is ma ¼ b3  ma ¼ M 2a b1a þ M 1a b2a . If we choose the base plane to coincide with the plane of symmetry for the thickness-wise mass distribution, we have qIo = 0. The resulting componential form of the equations of motion reads: o

ð28Þ

o

ð29Þ

@ 1 n11 þ @ 2 n12 þ f1 ¼ qh uott  e1 ; @ 1 n21 þ @ 2 n22 þ f2 ¼ qh uott  e2 ; @ 1 q1 þ @ 2 q2 þ f3 ¼ qh

o

uott

zq dz:

ð26Þ

z1

On the other hand, enforcing that the resultant moments with respect to the origin O of the fixed frame (sum of the resultant moments of the applied forces with the moments of the contact forces and the contact moments) to be equal to the rate of change of angular momentum, by exploiting (25), yields the local balance of angular momentum as

ð30Þ

21 Q 1  g 2 Q 2  1 N21 þ m @ 1 M 21 þ @ 2 M22  l31 M 11  l32 M 12 þ m o o   g2 N22 þ c  b1 ¼ qJ xt  b1 ;

ð31Þ

1 Q 1 þ g  1 N11  m  2 N12 @ 1 M 11 þ @ 2 M12 þ l31 M 21 þ l32 M 22 þ g 12 Q 2 þ c  b2 ¼ qJo xot  b2 : m

ð32Þ

It can be shown (here it is omitted for sake of conciseness) that the linearization of equations of motion (28)–(32) leads to the equations of motion of the Mindlin–Reissner theory. At the same time, with ad hoc approximations for the gradients of the deflection and in-plane displacements, the geometrically exact theory can furnish the equations of motion of the FVK theory. 2.5. Constitutive equations for thermo-viscoelastic plates The constitutive equations for thermo-viscoelastic plates are deduced from the constitutive equations of the three-dimensional theory according to which

b ðF; F; _ D#; rÞ; Tðr; tÞ ¼ T

z2

 e3 ;

_ D#; rÞ; Tðr; tÞ ¼ F  b SðE; E;

ð33Þ

where D# := #(t)  #0 is the temperature difference (#(t) is the current temperature and #0 is the temperature of a reference state, for example, the nominally stress-free state). Thus, the constitutive equations for the stress resultants can be obtained as

R z2

_ D#; r o ; zÞ  a dz; b ðF; F; T Rz _ D#; r o ; zÞ  a dz; b ðF; F; ^ a ðu; u; _ r o Þ ¼ b3  z 2 z T ma ðro ; tÞ ¼ m 1

^ a ðu; u; _ ro Þ ¼ na ðro ; tÞ ¼ n

z1

ð34Þ

where the strain and strain rate vectors are

1 ; m 2 ; g  1; l  2; l  12 Þ;  12 ; g 1; g 2; l u ¼ ðm

u_

_ 1 ; m _ 2 ; g _ 1 ; l _ 2 ; l _ 12 Þ; _ 12 ; g _ 1 ; g _ 2 ; l ¼ ðm

Fig. 4. The generalized stress resultants: local components of the contact force n1 and contact moment m1 in the intrinsic basis (b1, b2, b3). (N11, N21, Q1) are the tension, in-plane shear and transverse shear forces while (M21, M11) are the twisting and bending moments, respectively.

and the tensor-valued functions F and F_ indicate the deformation gradient and its rate parametrized in terms of the strain variables _ u and strain rates u. The constitutive equations can be equivalently expressed, by virtue of the relationship between the first and second Piola–Kirchhoff stress tensors, T = F  S, as

1654

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

R z2

_ D#; ro ; zÞ  a dz; Fb SðE; E; R z _ D#; ro ; zÞ  a dz; ^ a ðu; u; _ ro Þ ¼ b3  z 2 zF  b m SðE; E; 1 ^ a ðu; u; _ ro Þ ¼ n

z1

ð35Þ

b S 11 ¼ 1Em2 ½11 þ m22   a 1E m D#;

where, in consonance with the above stated notation, the tensorvalued functions E and E_ indicate the Green–Lagrange deformation _ and its rate parametrized in ðu; uÞ. It is convenient to introduce a constitutive relationship for the first-order moment of the stress vector ma as

^ a ðu; u; _ ro Þ ¼ m

Z

z2

neering elastic constants (E, m, G) and by assuming S = SLMeLeM and E = IJeIeJ, are in the form

_ D#; ro ; zÞ  a dz b ðF; F; zT

b S 22 ¼ 1Em2 ½m11 þ 22   a 1E m D#; b S 13 ¼ 2G13 ;

Z

z2

2

_ D#; ro ; zÞ  a dz; zF  b SðE; E;

ð36Þ

z1

so that the constitutive equation for the stress couple per unit ref^ a ðu; u; ^ a ðu; u; _ r o Þ ¼ b3  m _ r o Þ. erence length becomes m The componential form of the constitutive equations is effectively obtained by letting the stress tensor be expressed in the b¼T b LM eL eM so that fixed basis T

^ 11 :¼ n

R z2

^ 11 :¼ m

R z2

z1

z1

Tb 11 dz; z Tb 11 dz;

^ 21 :¼ n

R z2 z1

^ 21 :¼ m

Tb 21 dz;

R z2 z1

^1 :¼ q

z Tb 21 dz;

^ 31 m

R z2

Tb 31 dz; Rz :¼ z12 z Tb 31 dz;

z1

ð37Þ ^ 12 :¼ n

R z2

^ 12 :¼ m

R z2

z1

z1

Tb 12 dz; z Tb 12 dz;

n22 :¼

R z2 z1

^ 22 :¼ m

Tb 22 dz;

R z2 z1

q2 :¼

z Tb 22 dz;

^ 32 m

16 E¼ 4 2

jm 1 j2  1

3

g12

g 1 7 2  2 5; jm j2  1 g

ð41Þ

0 ^ M þ z2 l ^ M , ^M  l where jm M j2 ¼ ½m oM  m oM þ 2zm oM  l Consequently,

h

for

i

11 ¼ 12 m21 þ g 21 þ 2zl 1 m1 þ z2 l 21 þ ðm21 þ zl^ 21 Þ2  1 h

M = 1,

Tb 32 dz; z1 Rz :¼ 2 z Tb 32 dz: z1

2.6. Constitutive equations for linearly elastic isotropic plates To gradually gain insight into the constitutive characterization of geometrically nonlinear composite plates, we first discuss materials whose constitutive behavior is linearly elastic and isotropic also at moderately large strains which can be described by the St. Venant–Kirchhoff constitutive equations as

ð39Þ

where (k, l) are Lamé’s constants. By taking into account also thermoelastic effects, (39) is modified into b S ¼ k trðE  ET ÞIþ T T 2lðE  E Þ, where E : = aD#I is the thermally induced strain tensor for isotropic materials, a is the thermal expansion coefficient and I is the identity tensor. The thermal expansion coefficient is assumed independent of the temperature and strain levels, an acceptable assumption provided that the temperature differences are not too severe and the strain levels are moderately large. For loading conditions that do not cause appreciable normal stresses (S33  0), constitutive equations (39), in terms of the engi-

2.

; i

22 ¼ 12 ðm2 þ zl 2 Þ2 þ g 22 þ m212 þ 2zl^ 12 m12 þ z2 l^ 212  1 ; 12 ¼ 12 g12 ¼ 12 m 1  m2 ¼ 12 ½g 12 þ zðm12 l 1 þ m2 l^ 21 þ m1 l^ 12 þ m21 l 2 Þ   1l ^ 21 l ^ 12 þ l  2 Þ ; M3 ¼ 12 g  M ¼ 12 m oM  b3 : þz2 ðl ð42Þ

R z2

ð38Þ

b S ¼ k trðEÞI þ 2lE;

ð40Þ

E where E is Young’s modulus, G ¼ 2ð1þ mÞ is the shear modulus, m is Poisson’s ratio, and the thermally induced stress reduces to b S T ¼ a 1E m D#I. The Green–Lagrange deformation tensor is

z1

¼

b S 12 ¼ 2G12 ;

b S 23 ¼ 2G23 ;

To obtain the constitutive equations, we express the first Piola– Kirchhoff stress tensor in terms of the components of the second b kM ¼ F kj b Piola–Kirchhoff stress tensor according to T S jM ; ðM ¼ 1; 2Þ. Consequently, the generalized stress and couple resultants are constitutively expressed by (37) and (38). Closed-form expressions for the constitutive equations can be obtained by accounting for the fact that F ¼ F ð0Þ þ zF ð1Þ ; b S¼ b S ð1Þ þ z2 b S T ðzÞ. S ð2Þ  b S ð0Þ þ z b 2.7. Constitutive equations for general multi-layered laminated plates We let the material axes of orthotropy of the kth layer (henceforth, for notational simplicity, we drop the superscript k) be collinear with the axes of a rotated basis denoted ð ek1 ;  ek2 ;  ek3 Þ. We further restrict our attention to the case of composite laminated plates whose mid-plane coincides with the (e1, e2)-plane and the ply material axes ð ek1 ;  ek2 Þ are rotated by an angle uk with respect to the fixed axis e1 (see Fig. 5). The resulting constitutive equations, in the fixed basis, are rearranged as follows:

3 2 k b L11 S 11 6b 7 6 4 S 22 5 ¼ 6 4 b S 12 2

Lk12 Lk22

32 3 3 2 bT bT S 1 þ S 2 cos 2uk 11 76 7 7 6 6 bT bT k7 Lk26 7 5 4 22 5  4 S 1  S 2 cos 2u 5 D#; g12 b Lk66 S T2 sin 2uk Lk16

ð43Þ

Fig. 5. Geometry of the multi-layered plate with the numbering of the layers, and the coordinates.

1655

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

"

b S 31 b S 32

#

" ¼

Lk55

Lk54 Lk44

#



g31 ; g32

ð44Þ

where

i 1h b S T1;2 ¼ Ek1 ak1 Ek2 ak2 þ Ek1 mk21 ðak2 ak1 Þ =ð1  mk12 mk21 Þ: 2

ð45Þ

The elastic coefficients, by exploiting some trigonometric identities, can be effectively parametrized in terms of the ply angle uk as n1

Lknn ðuk Þ ¼ C k1 þ ð1Þ

C k2 cos 2 k þ C k3 cos 4 k ; k k 1 k k Ln6 ð Þ ¼  2 C 2 sin 2 þ ð1Þn C k3 sin 4 k ; n ¼ 1; 2; Lk12 ð k Þ ¼ C k4  C k3 cos 4 k ; Lk66 ð k Þ ¼ C k5  C k3 cos 4 k ; m1 Lkmm ð k Þ ¼ C k6 þ ð1Þ C k7 cos 2 k ; m ¼ 4; 5; Lk45 ð k Þ

u u u

u

u

v

u

u

u

u u

u u ¼ 12 C k7 sin 2uk ; ð46Þ

where

   3Lk11 þ 2Lk12 þ 3Lk22 þ 4Lk66 ; C k2 :¼ 12 Lk11  Lk22 ;   C knþ2 :¼ 18 Lk11  2Lk12 þ Lk22  4ðd1n  d3n ÞLk66 ; n ¼ 1; 3;   C k4 :¼ 18 Lk11 þ 6Lk12 þ Lk22  4Lk66 ; C k1 :¼

C k6 :¼

1 8



1 ðLk55 2

þ Lk44 Þ;

C k7 :¼

1 ðLk55 2

ð47Þ

 Lk44 Þ:

The constants C kn , expressed in terms of the elastic coefficients Lkij in the material coordinates, can be directly given in terms of the five material constants of the kth (transversely isotropic) layer, lk ¼ ðEk1 ; Ek2 ; mk12 ; Gk12 ; Gk23 ; ak1 ; ak2 Þ. 2.8. Equations of motion for arbitrarily oriented multi-layered laminated composite plates In the preceding section, we have parametrized the constitutive equations for the kth layer through the angle uk according to (43)– (47) thus obtaining

b Sðu; r; lk ; uk Þ; Sk ¼ b

Fig. 6. The MTS 809 Axial-Torsional Test System and the plate specimen with its clamping apparatus.

1 ; m 2 ; g  1; l  2; l  12 Þ:  12 ; g 1; g 2; l u ¼ ðm

ð48Þ

To obtain the constitutive equations for the generalized stress and moment resultants, the first Piola–Kirchhoff stress tensor in bk ¼ F  b the kth layer is obtained according to T S k ; in particular, k b b k  e2 , that we are interested only in the stress vectors T  e1 and T is,

Fig. 7. Schematic of the experimental layout of the plate and its clamping apparatus. The geometric dimensions are: l1 = 0.0254 m, l2 = 0.304 m, l3 = 0.0699 m, l4 = 0.254 m, h1 = 0.0254 m, h2 = 0.0127 m, h3 = 0.0635 m, h4 = 0.0025 m, and h is the plate thickness.

1656

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Fig. 8. The aluminum plate response to an increasing central point load in A. The midpoint deflections predicted by the nonlinear theory (solid line) and by the FVK (dashed line) and MR (straight line) theories are shown in (a) and (b) while the equilibrium paths of a few control points are shown in (c) and (d). The deflected configurations are shown in (e) and (f) (increasing load) and in (g) and (h) (end load). The circles indicate the experimental measurements.

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

1657

Fig. 9. The steel plate response to an increasing central point load in A. The midpoint deflections predicted by the nonlinear theory (solid line) and by the FVK (dashed line) and MR (straight line) theories are shown in (a) and (b) while the equilibrium paths of a few control points are shown in (c) and (d). The deflected configurations are shown in (e) and (f) (increasing load) and in (g) and (h) (end load). The circles indicate the experimental measurements.

1658

S kj1 ; Tb kn1 ¼ F nj b

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

S kj2 ; Tb kn2 ¼ F nj b

ð49Þ

where F is again parametrized in the generalized strains. According to (37) and (38), we can calculate the stress and moment resultants as follows:

^ 11 ; m ^ 11 Þ :¼ ðn

NL Z X k¼1

^ 21 ; m ^ 21 Þ :¼ ðn

zk

zk

NL Z X k¼1

zkþ1 zk

ð1; zÞ Tb k21 dz;

ð1; zÞ Tb k12 dz; 4.1. Linearly isotropic metallic single-layered plates

zkþ1

zk

ð1; zÞ Tb k11 dz;

ð1; zÞ Tb k31 dz;

zkþ1

NL Z X k¼1

^ 32 Þ :¼ ðq2 ; m

zkþ1

NL Z X k¼1

^ 22 Þ :¼ ðn22 ; m

zkþ1

zk

NL Z X k¼1

^ 12 ; m ^ 12 Þ :¼ ðn

zk

NL Z X k¼1

^1 ; m ^ 31 Þ :¼ ðq

zkþ1

sitivity of the load cell is 0.1 N while the sensitivity of the displacement sensor is 105 m. A Honeywell VL7A-1000 LVDT sensor was used to record the plate deflections in selected locations when the load is applied. The details of the clamping fixture are shown in Fig. 7 which portrays the top view (left) and the two cross sections (right). The nonlinear equilibrium paths were recorded using four specimen: two isotropic plates (one made of aluminum and the other made of steel) and two multi-layered composite laminated plates (made of carbon fibers and epoxy resin with different stacking sequence). The motivation for exploring isotropic singlelayered metallic plates arises from the need to compare the relative ranges of nonlinear behavior in isotropic and composite plates given the generality of the developed geometrically exact plate model.

ð1; zÞ Tb k22 dz;

ð1; zÞ Tb k32 dz;

ð50Þ

where NL is the number of layers, zk and zk+1 are the coordinates of the lower and upper surfaces of the kth layer. The stress resultants (50) are then substituted into equations of motion (28)–(32). 3. Numerical implementation of the plate model The geometrically exact plate model was numerically implemented using the computational platform COMSOL Multiphysics. This FEM-based solver has proved to be particularly amenable to the objectives of the present work since it allowed to enter the strong form of the problem, namely, the equations of motion, while the software automatically defines the associated weak form and the FEM discretization procedure. In using the so-called PDE mode feature, the equations enforcing balance of linear and angular momentum in their componential form (28)–(32) were implemented. The nonlinear kinematic relationships and the constitutive equations for both the isotropic and the multi-layered plates were explicitly implemented after their pre-processing in MATHEMATICA. A general orientation of the fibers can be set for the model, thus allowing fast parametric investigations while the number of the layers is prescribed a priori. The boundary conditions as well as the geometric and constitutive parameters are further prescribed. Meshes with an increasing level of refinement have been adopted to have mesh-independent results. By following the same procedure, the Föppl–von Kármán and the Mindlin–Reissner plate models were also implemented in COMSOL Multiphysics for comparison purposes. While various kinds of static, modal, and transient response analysis were performed [34], here we report only the nonlinear static analysis results which are systematically contrasted with the experimental measurements for isotropic single-layered and multi-layered composite laminated plates. 4. Experimental validation The experimental campaign was conducted by Michele Pasquali during a mobility fellowship provided to visit and perform research at Clarkson University in collaboration with Mechanical and Aeronautical Engineering Professor Pier Marzocca. This section discusses the experimental test design, setup, procedures employed for data collection and analysis. The experimental rig consists of the plate test specimen with the supporting base fixed to an MTS 809 Axial-Torsional Test System with maximum static load of 250 kN (maximum dynamic load: 200 kN) and frequency range [0, 10] Hz (see Fig. 6). The sen-

The specimen consist of two homogeneous isotropic square metallic plates (provided by the manufacturer MACMASTER-CARR) with side length of 0.254 m (tolerance = 1.59  105 m) and thickness equal to 0.8128 mm (tolerance = 7.62  106 m). One of the plates is made of spring steel (tough ware resistant spring steel 1095) whose mechanical properties certified by the manufacturer are: E = 206.82 GPa, Poisson’s ratio m = 0.3, mass density q = 7833.41 kg/m3, thermal expansion coefficient a = 13.32  106C1, yield strength Sy = 551.58 MPa, and tensile strength equal to 979.06 MPa. The other plate is made of aluminum AA3003 with the same geometrical features. The mechanical properties of the aluminum plate are: E = [70, 80] GPa, m = 0.3, q = 2730 kg/m3, yield strength Sy = 125 MPa, and tensile strength equal to 130 MPa. The thickness-to-span ratio for both plates is 0.0032. The load test protocols feature a displacement control loading pattern with a monotonic loading ramp whose rate was kept at 0.1 mm/s. The sampling frequency was set to 100 Hz. The end displacement of the ramp was chosen to induce a maximum calculated stress equal to 90% of the yield stress. For the steel plate, the end displacement was set to 1.805 mm whose associated measured load was 132.93 N and the maximum estimated tensile stress was 496.42 MPa (=90%Sy). On the other hand, for the aluminum plate, the end displacement was fixed to 1.388 mm whose associated load was 27.94 N and the calculated tensile stress was 112.5 MPa (=90%Sy). The equilibrium load paths of the midpoint A are shown, for both the aluminum and the steel plate, in Figs. 8 and 9a, where the theoretical predictions of the nonlinear theory are compared with those of the FVK and MR theories. The filled circles denote the experimental results, where the radius of the circles is proportional to the sensitivity of the displacement sensor. A close agreement between the nonlinear theory and the experimental results can be observed. The FVK theory, as expected, predicts lower displacements due to the more constrained nature of the underlying approximate nonlinear model. In fact, the only nonlinearity included in the model is represented by the membrane stretching Table 1 Geometrical and mechanical properties of the Graphil 34-600/NTC301 and T600s/G91 plates. Composite plate

Graphil 34-600/NTC301

T600s/G91

Span Thickness Lamination sequence Young’s modulus E1 Young’s modulus E2 Poisson’s ratio m12 Poisson’s ratio m21 Density

0.254 m 2.54  103 m [0°/45°/90°/135°]s 137.137 GPa 9.308 GPa 0.304 0.017 1568 kg/m3

0.254 m 2.54  103 m [0°/15°/15°/0°]s 134.448 GPa 9.102 GPa 0.300 0.016 1501 kg/m3

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

1659

Fig. 10. The Graphil 34-600/NTC301 plate response to an increasing central point load in A. The midpoint deflections predicted by the nonlinear theory (solid line) and by the FVK (dashed line) and MR (straight line) theories are shown in (a) and (b) while the equilibrium paths of a few control points are shown in (c) and (d). The deflected configurations are shown in (e) and (f) (increasing load) and in (g) and (h) (end load). The circles indicate the experimental measurements.

1660

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Fig. 11. The T600s/G91 plate response to an increasing central point load in A. The midpoint deflections predicted by the nonlinear theory (solid line) and by the FVK (dashed line) and MR (straight line) theories are shown in (a) and (b) while the equilibrium paths of a few control points are shown in (c) and (d). The deflected configurations are shown in (e) and (f) (increasing load) and in (g) and (h) (end load). The circles indicate the experimental measurements.

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

effect while the flexural and twist curvatures are represented in their linearized form thus neglecting the nonlinear curvature contributions which exert a relaxing effect on the stiffness, and in general, on the motions of the plate. In Figs. 8b and 9b the nondimensional counterpart of the previous graphs for the aluminum and the steel plates are shown. In Figs. 8c, d and 9c, d, the equilibrium load paths of a few control points aligned along the mid-line and the diagonal line of the plates are depicted. As expected, the various points exhibit eminently increasing stiffnesses as they move closer to the clamping boundaries. In Figs. 8e, f and 9e, f the deflected configurations of the mid-line and diagonal line of the plates subject to a central point load of increasing magnitude are illustrated: for the aluminum plate the sequence for the increasing applied load is [5, 10, 15, 20, 25] N; on the other hand, for the steel plate, the sequence of the increasing applied load is [25, 50, 75, 100, 125] N. In Figs. 8g, h and 9g, h the deflected configurations of the mid-line and diagonal line of the plates subject to the maximum value of the load applied at the center of the plates are shown: for the aluminum plate this load is 27.94 N; for the steel plate the load is 132.93 N. In these plots the dashed lines indicate the predictions of the FVK theory while the solid lines denote the results of the geometrically exact theory. The filled circles denote the experimental results. Again, the agreement between the proposed nonlinear model and the experiment is remarkable. 4.2. Multi-layered laminated composite plates In this case the specimen are two multi-layered square plates (provided by Rock West Composites) with side length of 0.254 m (tolerance = 5.08  104 m) and thickness equal to 2.54 mm (tolerance = 2.54  104 m). One of the plates is made of Graphil 34-600 carbon fibers and Newport NTC301 epoxy matrix, symmetrically distributed in eight different layers: the stacking sequence of the plate is [0°/45°/90°/135°]s. The mechanical properties are: E1 = 137.137 GPa, E2 = 9.308, m12 = 0.304, m21 = 0.017, q = 1568 kg/ m3. The other plate is made of Torayca T600s carbon fibers and Estron G91 epoxy matrix, symmetrically distributed in eight layers: the stacking sequence of the plate is [0°/15°/15°/0°]s. The mechanical properties are: E1 = 134.448 GPa, E2 = 9.102, m12 = 0.3, m21 = 0.016, q = 1501 kg/m3. The thickness-to-span ratio for both composite plates is 1/100. A summary of the data of the plates is given in Table 1. The static tests have been run in displacement control with a monotonic loading ramp whose rate was 0.01 mm/s. The sampling frequency was set to 100 Hz. The end displacement of the ramp was chosen equal to 3 mm for both plates, so as to induce a clear geometrically nonlinear response while remaining within the linear elastic range. The end load for Graphil 34-600/NTC301 plate was 555.56 N while that for T600s/G91 was 488.06 N. The dimensional equilibrium load paths of the midpoint A are shown, for the Graphil 34-600/NTC301 and T600s/G91 plates, respectively, in Figs. 10a, b and 11a, b where the theoretical predictions of the nonlinear theory are compared with those of the FVK and MR theories. The filled circles denote the experimental results with the radius of the circles being proportional to the sensitivity of the displacement sensor. A close agreement between the proposed nonlinear theory and the experimental results can again be observed in the context of the multi-layered plates. The FVK theory, as expected, again predicts lower displacements due to the more constrained nature of the approximate nonlinear model. In Figs. 10c, d and 11c, d, the equilibrium load paths of a few control points aligned along the mid-line and the diagonal line of the plates are presented. Also in this case, the various points exhibit eminently increasing stiffnesses as we move closer to the clamping boundaries. In Figs. 10e, f and 11e, f the deflected configurations of the mid-line and diagonal line of the plates subject to a central

1661

point load of increasing magnitude are illustrated: for the Graphil 34-600/NTC301 plate the sequence of the increasing applied load is [100, 200, 300, 400, 500] N while, for the T600s/G91 plate, the load sequence is [90, 180, 270, 360, 450] N. In Figs. 10g, h and 11g, h the deflected configurations of the mid-line and diagonal line of the plates subject to the maximum central load are shown: for the Graphil 34-600/NTC301 plate the applied load is 555.56 N while for the T600s/G91 plate the load is 488.06 N. Also in these plots, the dashed lines indicate the predictions of the FVK theory, the solid lines denote the results of the geometrically exact theory, and the filled circles denote the experimental results. An interesting observation is that for the Graphil 34-600/NTC301 plate, the predictions of the FVK theory are closer to the geometrically exact theory in the whole range while this is no longer true for the T600s/ G91 plate. This is likely due to the higher stiffness induced by the stacking sequence in the Graphil 34-600/NTC301 plate which allows less nonlinear curvature effects to be exhibited in the deformation as it is clear by comparing Figs. 10g and 11g. 5. Comments The fidelity of two-dimensional theories of plates has been an outstanding question for over a century, first within the context of isotropic single-layered plates, more recently, also in relation to multi-layered composite laminated plates. While fully intrinsic theories à la Cosserat, which view the plate from the onset as a two-dimensional surface with additional extrinsic directors (in the sense of a Cosserat surface) and invoke concepts from differential geometry, have their own self-consistency, the difficulty arises in identifying the relationship between the constitutive laws of the thus reduced two-dimensional object and the constitutive properties of the three-dimensional physical counterpart. Semi-intrinsic approaches have been effectively employed by enforcing physically reasonable constitutive assumptions on the kinematics of the three-dimensional plate-like structures. In so doing, the constitutive relationships of the reduced two-dimensional theory can be obtained from the three-dimensional elasticity theory. This approach has been widely used in plate mechanics, especially in trying to establish nonlinear theories. Invariably, within the context of a yet rich literature on plate theories, there is a partial lack of efforts in the assessment of the reduced theories accuracy to be pursued by comparing their outcomes with experimental results, in particular, when the plates are forced in the nonlinear range. In this work, by resorting to a semi-intrinsic approach, we have devised a geometrically exact model of thin isotropic and multi-layered composite plates which relies fundamentally on one kinematic constraint, that of the rigidity of transverse material fibers. Other than this internal kinematic constraint and the assumption of linear elasticity in the constitutive behavior, no other restrictions have been enforced on the plate kinematics. The five reduced equations of motion have a form amenable to a numerical treatment within the available computational platforms which perform independently the discretization on the associated weak form. We have documented the nonlinear hardening-type response of fully clamped isotropic metallic and eightlayer composite laminated plates to a central point load of increasing magnitude. The selected plates are prominently thin; the isotropic plates exhibit a thickness-to-span ratio of the order of 3/ 1000 while the multi-layered plates are slightly thicker with a ratio equal to 1/100. A careful experimental campaign developed at Clarkson University, through measurements of the deflections of nine control points of the plate surface under various load levels, has shown the accuracy of the proposed theory in the full nonlinear range up to deflections of the order of one-hundredth of the plate span. We have observed clear deviations between the linear theory

1662

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663

(Mindlin–Reissner) and the nonlinear theories (both the present theory and the Föppl–von Kármán theory) at deflection amplitudes in the range 3  103–4  103 for isotropic metallic plates or 1.5  103–2  103 for the eight-layer composite plates. An important observation is that above the load levels where the response becomes clearly nonlinear, the Föppl–von Kármán theory predicts lower deflections than the present nonlinear theory or the experimental results. This is indeed expected since the FVK theory introduces the nonlinearity solely in the membrane stretching effects while it does not relax the deflections through the nonlinear curvatures and, less importantly, through the shear deformations. This is confirmed in the comparison between the two eight-layer laminated plates (Graphil 34-600/NTC301 with [0°/45°/90°/135°]s and T600s/G91 with [0°/15°/15°/0°]s). For the first plate the deviation between the FVK theory and the proposed geometrically exact theory is nearly negligible in the explored range likely due to the higher stiffness arising from the stacking sequence in the Graphil 34-600/NTC301 plate which allows less nonlinear curvature effects to be exhibited in the plate deformation. Of course, the nonlinear behaviors, here observed and documented in the fully clamped test specimen, depend significantly on the boundary conditions since different nonlinear load-carrying mechanisms may be activated in the presence of different constraining devices. Acknowledgments This work was partially supported by a FY 2009 AST Sapienza Grant and a 2009 Sapienza Fellowship for International Mobility granted to Michele Pasquali who visited Clarkson University (Potsdam, NY, USA) to conduct the experimental work here reported. Partial support provided to Prof. Marzocca by the ARMY Research Office, grant W911NF-05-1-0339 is also acknowledged. The authors would like to thank Prof. Pier Marzocca who helped with the experimental setup design and the protocols. The first author is deeply indebted to Prof. Stuart S. Antman for having shared with him some of the remarkable intricacies of the semi-intrinsic theories of nonlinear elasticity.

To obtain the orthogonal tensor R that transforms o o o ðb1 ; b2 ; b3 Þ ¼ ðe1 ; e2 ; e3 Þ into (b1, b2, b3), we consider the sequence of rotations /2 ? /1. We first rotate the basis about e1 by /2 and ð1Þ ð1Þ ð1Þ denote ðb1 ; b2 ; b3 Þ the obtained unit vectors. This choice is consistent with the traditional notation prevailed in the literature on ð1Þ plate theories. We then rotate this basis around b2 by the angle /1. The ensuing orthogonal matrix, R = R(2)  R(1), has the following entries:

R21 ¼ 0;

R12 ¼ sin /1 sin /2 ;

R22 ¼ cos /2 ;

R31 ¼  sin /1 ;

R23 ¼  sin /2 ;

ð51Þ

R33 ¼ cos /1 cos /2 :

The components of the unit vectors (b1, b2, b3) in the basis (e1, e2, e3) are

b1 ¼ cos /1 e1  sin /1 e3 ; b2 ¼ sin /1 sin /2 e1 þ cos /2 e2 þ cos /1 sin /2 e3 ;

ð52Þ

b3 ¼ cos /2 sin /1 e1  sin /2 e2 þ cos /1 cos /2 e3 : The actual position vectors of the material points at ro = x1e1 + x2e2 and at r = ro + ze3 are

po ¼ ðx1 þ uo1 Þe1 þ ðx2 þ uo2 Þe2 þ uo3 e3 ; p ¼ ðx1 þ uo1 þ z cos /2 sin /1 Þe1 þ ðx2 þ uo2  z sin /2 Þe2 þðuo3 þ z cos /1 cos /2 Þe3 :

ð54Þ ð55Þ

Appendix B. Compatibility and independent generalized strains for plates The four vector-valued fields m o1 ðxM ; tÞ; m o2 ðxM ; tÞ; l1 ðxM ; tÞ; l2 ðxM ; tÞ are not independent (in a differential sense) since they have to satisfy the compatibility conditions. To this end, recall the definitions: m oM = @ Mpo and @ Mbk = lM  bk. For continuity requirements of the vector-valued fields po(xM,t) and bk(xM, t) in Bo , according to the Schwarz Theorem, the following must hold: @ 2@ 1po = @ 1@ 2po and @ 2@ 1bk = @ 1@ 2bk. The above equalities imply that

@ 2 m o1 ¼ @ 1 m o2 ;

@ 2 ðl1  bk Þ ¼ @ 1 ðl2  bk Þ:

ð56Þ

Eqs. (56) yield six independent compatibility equations. The projection of (56)1 in the (b1, b2, b3)-basis yields

1  @ 1 m 12 þ l31 m 2  l32 m 21 þ l12 g  1  l11 g  2 ¼ 0; @2m       2 ¼ 0; @ 2 m21  @ 1 m2 þ l32 m1  l31 m12  l22 g1 þ l21 g 1  l21 m 2 þ l22 m 21 þ l11 m 12 ¼ 0: 1  @1g  2  l12 m @2g

ð57Þ ð58Þ ð59Þ

On the other hand, the compatibility equations for the curvatures are given by

@ 2 l11  @ 1 l12  l31 l22 þ l32 l21 ¼ 0; @ 1 l22  @ 2 l21 þ l11 l32  l12 l31 ¼ 0;

ð60Þ ð61Þ

@ 1 l32  @ 2 l31  l11 l22 þ l12 l21 ¼ 0:

ð62Þ

To make the generalized nonlinear strain measures coalesce into the well-known linear strains for shear deformable plates, we introduce the twist curvature as

l 12 :¼ l12  l21 :

ð63Þ

Eqs. (59) and (63) are used to obtain

ð64Þ ð65Þ

Eqs. (57) and (58) can be solved for (l31, l32). The outcomes, along with (64) and (65) are substituted into the curvature compatibility Eqs. (60)–(62) which turn out to be expressed in terms of the eight generalized strains

1 ; m 2 ; g  1; l  2; l  12 Þ:  12 ; g 1 ; g 2; l u ¼ ðm This proves that only five generalized strains are truly independent. References

R13 ¼ cos /2 sin /1 ;

R32 ¼ cos /1 sin /2 ;

to11 ¼ 1 þ @ 1 uo1 ; to21 ¼ @ 1 uo2 ; to31 ¼ @ 1 uo3 ; to12 ¼ @ 2 uo1 ; to22 ¼ 1 þ @ 2 uo2 ; to32 ¼ @ 2 uo3 :

l12 ¼ ð@ 2 g 1  @ 1 g 2 þ l11 m12 þ l22 m21 þ l 12 m2 Þ=ðm1 þ m2 Þ; l21 ¼ ð@ 2 g 1  @ 1 g 2 þ l11 m12 þ l22 m21  l 12 m1 Þ=ðm1 þ m2 Þ:

Appendix A. Special forms of the rotation tensor

R11 ¼ cos /1 ;

The stretch vectors can be calculated in a straightforward way as m oM = @ Mpo = tokM ek where

ð53Þ

[1] Koiter WT. A consistent first approximation in the general theory of thin elastic shells. In: Proceedings of the symposium on theory of thin elastic shells (IUTAM), Delft, August 24–28, North Holland, Delft, NL; 1959. p. 1233. [2] Koiter WT. Comment on: the linear and non-linear equilibrium equations for thin elastic shells according to the Kirchhoff–Love hypotheses. Int J Mech Sci 1970;12:663–4. [3] Reddy JN. Mechanics of laminated composite plates and shells. 2nd ed. Boca Raton, FL: CRC Press; 2004. [4] Batra RC. Higher order shear and normal deformable theory for functionally graded incompressible linear elastic plates. Thin Wall Struct 2007;45:974–82. [5] Fiedler L, Lacarbonara W, Vestroni F. A generalized higher-order theory for buckling of thick multi-layered composite plates with normal and transverse shear strains. Compos Struct 2010;92:3011–20. [6] Carrera E. Transverse normal stress effects in multilayered plates. J Appl Mech 1999;66:1004–12. [7] Carrera E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. Arch Comput Methods Eng 2002;9:87–140. [8] Lee J. Thermally induced buckling of laminated composites by a layer-wise theory. Comput Struct 1997;65:917–22.

W. Lacarbonara, M. Pasquali / Composite Structures 93 (2011) 1649–1663 [9] Fiedler L, Lacarbonara W, Vestroni F. A general higher-order theory for multilayered, shear-deformable composite plates. Acta Mech 2009;209:85–98. [10] Synge JL, Chien WZ. The intrinsic theory of elastic shells and plates. Theodore von Karman anniversary volume. California Institute of Technology; 1941. p. 103–20. [11] Naghdi PM. The theory of shells and plates. In: Truesdell S, Flügges S, editors. Encyclopedia of physics, vol. VI a/2. Springer-Verlag; 1972. p. 425–640. [12] Simo JC, Marsden JE, Krishnaprasad PS. The Hamiltonian structure of nonlinear elasticity: the material and convective representations of solids, rods and plates. Arch Ration Mech Anal 1988;104:125–83. [13] Antman SS. Problems of nonlinear elasticity. New York: Springer; 2005. [14] Nayfeh AH, Pai PF. Linear and nonlinear structural mechanics. New York: Wiley; 2004. [15] Podio-Guidugli P. An exact derivation of the thin plate equation. J Elasticity 1989;22:121–33. [16] Love AEH. The mathematical theory of elasticity. Cambridge: Cambridge University Press; 1906. [17] Mindlin RD. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. J Appl Mech 1951;38:31–8. [18] Yu W. Mathematical construction of a Reissner–Mindlin plate theory for composite laminates. Int J Solids Struct 2005;42:6680–99. [19] Wenbin Yu W, Kimb JS, Hodges DH, Chod M. A critical evaluation of two Reissner–Mindlin type models for composite laminated plates. Aerospace Sci Technol 2008;12(5):408–17. [20] Föppl A. Vorlesungen über technische Mechanik. Leipzig, Germany: B G Teubner, Bd. 5; 1907. [21] von Kármán T. Festigkeitsproblem im Maschinenbau. Encyk D Math Wiss 1910;IV:311–85. [22] Ciarlet PG, Destuynder P. A justification of a nonlinear model in plate theory. Comp Methods Appl Mech Eng 1979;17/18:227–58. [23] Shufrin I, Rabinovitch O, Eisenberger M. Elastic nonlinear stability analysis of thin rectangular plates through a semi-analytical approach. Int J Solids Struct 2009;46:2075–92.

1663

[24] Antman SS. Global properties of buckled states of plates that can suffer thickness changes. Arch Ration Mech Anal 1990;110:103–17. [25] Sugimoto N. Nonlinear theory for flexural motions of this elastic plate. J Appl Mech 1981;48:377–82. [26] Johnson Jr MW, Urbanik TJ. A nonlinear theory for elastic plates with application to characterizing paper properties. J Appl Mech 1984;51:146–52. [27] Hodges DH, Atilgan AR, Danielson DA. A geometrically nonlinear theory of elastic plates. J Appl Mech 1993;60:109–1126. [28] Hodges DH, Yu W, Mayuresh JP. Geometrically-exact, intrinsic theory for dynamics of moving composite plates. Int J Solids Struct 2009;46:2036–42. [29] Demasi L, Yu W. Assess the accuracy of the variational asymptotic plate and shell analysis (VAPAS) using the generalized unified formulation (GUF). In: Proceedings of the 50th structures, structural dynamics, and materials conference, Palm Springs, CA, 4–7 May; 2009. [30] Carrera E, Parisch H. An evaluation of geometrical nonlinear effects of thin and moderately thick multilayered composite shells. Compos Struct 1997;40(1):11–24. [31] Andrade LG, Awruch AM, Morsch IB. Geometrically nonlinear analysis of laminate composite plates and shells using the eight-node hexahedral element with one-point integration. Compos Struct 2007;79(4):571–80. [32] Han S-C, Tabiei A, Park W-T. Geometrically nonlinear analysis of laminated composite thin shells using a modified first-order shear deformable elementbased Lagrangian shell element. Compos Struct 2008;82(3):465–74. [33] Formica G, Lacarbonara W, Alessi R. Vibrations of carbon nanotube-reinforced composites. J Sound Vib 2010;329:1875–89. [34] Pasquali M, Lacarbonara W, Marzocca P. Nonlinear system identification of thin plates using higher-order spectra: numerical and experimental investigations. J Mech Syst Signal Pr, submitted for publication.