Characterization of anisotropic, non-homogeneous plates with piezoelectric inclusions

Computers and Structures 83 (2005) 1171–1190 www.elsevier.com/locate/compstruc

Characterization of anisotropic, non-homogeneous plates with piezoelectric inclusions Pierangelo Masarati *, Gian Luca Ghiringhelli Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, via La Masa 34, 20156 Milan, Italy Accepted 18 October 2004 Available online 2 March 2005

Abstract A formulation for the determination of the stiffness, the dielectric and the piezoelectric properties of arbitrarily anisotropic and non-homogeneous laminated composite plates is presented. The linear elastic problem of a sample fibre, an ideal line through the thickness of the plate, is solved in terms of influence coefficients for both the generalized deformations (i.e. linear and angular strains of the plate) and the warping unknowns. The generalized deformations and the warping unknowns, as well as their plane-wise spatial derivatives, are preserved in analytical form. The in- and outof-plane warping displacements are described by means of arbitrary one-dimensional finite elements, accounting for the full three-dimensional material constitutive law; as a consequence, the desired level of accuracy can be achieved in the determination of interlaminar and intralaminar stresses and strains. The primary result is represented by the stiffness parameters of a general plate, which implicitly account for plane stress effects and for both interlaminar and intralaminar equilibrium and compatibility. Secondary results are represented by the full three-dimensional strain and stress fields due to the local internal force and moment fluxes. The extension to the inclusion of laminæ of piezoelectric materials in the lamination sequence of the plate is addressed. Numerical examples related to the mechanical and piezoelectric characterization of laminated plates are discussed.
2005 Elsevier Ltd. All rights reserved. Keywords: Laminated plates; Structural analysis; Plate characterization; Active materials; Piezoelectric structures

1. Introduction The plate model represents an important structural analysis tool. It can be efficiently used to model structures with a significant two-dimensional characterization, saving computational cost while the implied accuracy loss is acceptable. The elastic properties of a *

Corresponding author. Tel.: +39 2 2399 8365; fax: +39 2 2399 8334. E-mail addresses: [email protected] (P. Masarati), [email protected] (G.L. Ghiringhelli).

plate made of isotropic material can be easily calculated by hand, but when the properties of the material and the layout of the lamination sequence are a bit odder, adequate analysis tools are required to obtain reliable parameters for general purpose finite element (FE) plate or shell elements. This is particularly true when highly heterogeneous materials concur in the lamination sequence, and significantly when smart materials, like piezoelectric layers, are embedded in passive host structures. In these cases, the elastic interaction between different materials, including the bonding layers represented by conducting laminæ and glue, highly influences

0045-7949/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2004.10.017

1172

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

Nomenclature Symbols linear elastic constitutive law matrix C(E) D dielectric displacement E electric field F distributed forces per unit surface F plate fibre compliance matrix K plate fibre stiffness matrix dLe, dLi virtual external and internal work Mij, Ci, Lij, E, Ri, Aij, i,j=x,y internal work matrices N(z), B(z) generic shape function and shape function gradient Pi, i=x,y generalized warping forces Q generalized free charge matrix Ti, i=x,y arm matrices V electric potential X; Y; Z collocation matrices Z rigid translation and rotation matrix e linear piezoelectric constitutive law matrix ei, i=x, y, z coordinate vectors f distributed forces per unit volume h plate thickness pi, i=x,y stress projections in the coordinate directions qk, k=1,n electric conductor free charge r(x, y) linearized rigid displacement and rotation vector s displacement of an arbitrary point u(x, y) warping discrete unknowns v rigid part of the displacement of an arbitrary point w warping part of the displacement of an arbitrary point x, y, z coordinates of an arbitrary point H0i, i=x,y reference internal force and moment fluxes P position of an arbitrary point

the overall elastic and piezoelectric properties of the assembly, and an integrated and accurate analysis approach is mandatory. The conventional laminate theory (CLT) represents such an analysis tool, since it allows to determine reasonably accurate values for the membrane and the bending properties of a plate made of anisotropic materials. It consists in calculating the in-plane, bending and coupling properties as zero, first and second order moments of the thickness of each lamina, weighed by the elastic modulus of the materials, sometimes under the assumption that the bending stiffness of the single lamina is negligible. But the CLT fails when calculating the transverse shear coefficients and in general it is unable to determine the three-dimensional stress and strain state through the thickness; moreover, it fails to comply with the equilibrium and the compatibility at

generic rigid displacement and rotation multipliers
(e) linear dielectric constitutive law matrix ij, i=x, y, z; j=x,y pseudo-strain components eij, i=x, y, z; j=x,y strain components ei, cij, i=x,y; j=y,z engineering strain components #i, i=x,y internal force and moment fluxes ji, jij, i,j=x,y curvature components li, lij, i=x,y; j=y moment flux components acting on face i (about direction j) rij, i=x, y, z; j=x,y stress components /i, i=x, y, z components of a rigid rotation ui, uij, i=x,y; j=y,z force flux components acting on face i (in direction j) vi, i=x, y, z components of a rigid displacement wi, i=x,y linear and angular strains of the plate a

Operators partial derivative with respect to i (Æ)/i (Æ)· vector product o(Æ) compatibility operator that yields the strains from the displacements $(Æ) gradient operator (for electric potential) (Æ) matrix pseudo-inverse Acronyms AFC active fibre composite ATR active twist rotor CLT conventional laminate theory DAE differential algebraic equations FE, FEA, FEM finite elements, FE analysis, FE method PDE partial differential equations RHS right-hand side (of an equation)

the boundaries between the laminæ. The transverse shear properties can be quite important for different reasons. In fact, if the transverse shear is neglected, the displacements of a relatively thick plate can be underestimated, as well as its natural frequencies, buckling loads, and so on. Moreover, if damage criteria for composite materials must be used, the accurate determination of the three-dimensional strains and stresses is of paramount importance. Different approaches to the problem have been proposed, generally based on higher-order, through-the-thickness and plane-wise series expansion of the displacement unknowns. The literature on the subject is so broad that even an essential review would be out of the scope of this paper. Among the others, the works of Krishna and Vellaichamy [1], Kant and Pandya [2], Barboni et al. [3],

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

Di Sciuva [4], and Carrera [5] are mentioned as wellknown modeling techniques; a good reference may be found in [6]. In fact the domain of a composite plate sample fibre is quite simple, typically consisting in a straight line segment through the thickness, divided into slices where the material and the lamination properties change. If a linear displacement function is assumed through the whole thickness, as in the CLT, usually a saving in computational effort is achieved due to the limited number of degrees of freedom, at the cost of a rough modeling of strains and stresses, since, in general, both equilibrium and compatibility at the interfaces between the laminæ are violated. A further enhancement consists in using independent functions to model the displacements through each lamina, with compatibility constraints at the interfaces (maybe implicitly accounted for by the shape functions). In this case, if a C0 formulation is considered, both equilibrium and compatibility can be satisfied if appropriate functions are chosen. The formulation presented in this paper lies in the mainstream of this approach, but it operates an important distinction: the analytical expressions of the global displacements and deformations of the plate are considered separately from the warping, which is discretized following the finite element approach. On one hand, the global deformations are expressed in terms of the rigid displacements and rotations of the reference surface, so they can be directly used as global unknowns in view of a finite element discretization of the plate. On the other hand, the warping allows the displacement field in the plate to satisfy the natural boundary conditions both at the top and bottom surfaces of the plate and at the interfaces between the laminæ. The solution of the equilibrium problem of a sample fibre subject to unity internal loads that implicitly satisfy the five indefinite equilibrium equations (assuming that the drill equilibrium equation is degenerate, as discussed later, and thus implicitly satisfied) leads to the direct determination of the compliance properties of the plate, which are the inverse of the stiffness properties. At the same time, the local solution in terms of strains and stresses related to both the global deformation of the plate and to the warping effects is achieved. This formulation is an extension of the semi-analytical characterization of a composite beam section presented by Giavotto et al. in 1983 [7] and subsequently reviewed in a more general form. The resulting stiffness properties can be used in a conventional finite element formulation to determine the three-dimensional solution in terms of localized internal force fluxes, which can be subsequently substituted into the fibre solution to determine the strain and stress fields. Throughout this paper, the strains are assumed to be infinitesimal and only linear constitutive laws are considered, since the elastic characterization of the plate is addressed. The formulation can be easily generalized to a non-linear, pre-stressed and pre-strained configura-

1173

tion; in this sense the present work discusses a local linearization about an arbitrary configuration.

2. Problem description The notation used throughout the paper for most of the entities has been inherited from the semi-analytical formulation for the analysis of beam sections revisited by Ghiringhelli and Mantegazza [8], that is here applied to the characterization of plates. In fact, while the beam section characterization essentially was essentially based on the separation of the three-dimensional beam into a two-dimensional beam section manifold lying on a one-dimensional structured support, in this work the three-dimensional plate formulation is separated into a one-dimensional plate fibre manifold lying on a twodimensional structured support. The dependence on the in-plane position is thus treated analytically, while the dependence on the transverse position is discretized in a finite element way. A ‘‘sample fibre’’ of a plate is represented by a segment on a line, oriented as the axis, z, that is locally transverse to the plate, bounded by two points: the intersections of the transverse axis with the lower and upper faces of the plate. Two directions, x and y, mutually orthogonal and orthogonal to the fibre, identify the reference plane of the plate. The thickness of the plate is assumed to be constant with respect to the plane-wise directions in the reference configuration; in case of a laminated plate, the thickness of each lamina is assumed to be constant. The fluxes of internal forces and moments, collected in the 6 · 1 vectors #x, #y, respectively defined as 8 8 9 9 uxy > uxx > > > > > > > > > > > > > > > > > uyy > uyx > > > > > > > > > > > > = zy zx #x ¼ ; #y ¼ ; > > lxy > lxx > > > > > > > > > > > > > > > > > > > lyy > l > > > > > > > > : yx > : ; ; l lzx zy act on the sides of the fibre normal to the x and y axes. They must satisfy the local indefinite equilibrium equation, #x=x þ #y=y  T Tx #x  T Ty #y ¼ 0;

ð1Þ

where matrices Ti, i = x,y, defined as   0 ex  0 ey  ; Ty ¼ ; Tx ¼ 0 0 0 0 result from the differentiation of the moment transport operator  I 0 ; Uðx; yÞ ¼ xex  þyey  I

1174

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

namely T i ¼ U T=i ; ei are the unit vectors in the coordinate directions, while the operator (Æ)· performs the vector product in matrix notation: being a and b two vectors, a· is the matrix that yields the vector product a · b when it pre-multiplies b. Matrices Ti represent the gradient of the arm by which the shear fluxes on the side faces of the fibre participate in the moments equilibrium. A finite contribution of the shear flux to the equilibrium of the in-plane bending moments exists, since the thickness of the fibre is finite and thus it is able to counteract a bending moment; but it has infinitesimal plane-wise dimension, so it cannot react a concentrated moment in the transverse (drill) direction. As a consequence, the reaction moment about the z-axis, ez, is null, as well as its derivatives, and the moment equation in the z direction simply states the symmetry of the in-plane shear fluxes. Eventually, the contribution of distributed loads to the equilibrium equation will be taken into account, and the equilibrium equation itself will be derived by means of the Virtual Work Principle. At present, without any loss in generality, the fibre is assumed to bear internal loads only. In a classical displacement approach to the solution of the elastic problem of a plate, the internal force and moment fluxes are expressed in terms of the linear and angular strains by means of an elastic constitutive law. These, in turn, can be written in terms of displacements and displacement rates by means of a compatibility relationship. The linearization of an arbitrary displacement of a reference point along the fibre is represented by the vector r(x, y) = [vx, vy, vz, /x, /y, /z]T. Its gradient represents a measure of the deformation of the reference plane. The generalized deformations wx, wy of the plate, i.e. the linear and angular strains, are defined as: 8 8 9 9 xx > xy > > > > > > > > > > > > > yx > yy > > > > > > > > > > > > > > > > > < < zx = zy = wx ¼ ; wy ¼ ; > > / > / > > > > x=x > > x=y > > > > > > > > > > > > > /y=x > /y=y > > > > > > > > > : : ; ; /z=x /z=y where
results from the perturbation of the gradient of the displacement (i.e. the so-called deformation gradient) with respect to a reference configuration; as a consequence, in general it is not symmetric; in fact ii corresponds to the axial strain eii, while ij + ji = 2eij. The above presented generalized deformations depend on the displacement of the reference point and on its gradient wx ¼ r=x þ T x r;

wy ¼ r=y þ T y r:

ð2Þ

They must be null in case of rigid motion of the reference point; by considering the rigid linearized translation and rotation r = (I  xTx  yTy)a, the wi are

identically null regardless of the multiplier a if the property TiTj = 0 is accounted for. This implies that the generalized deformations defined by Eq. (2) are a meaningful measure of the deformation of the reference plane. A compatibility problem arises from the consideration that the 12 generalized deformations descend from six global displacements and rotations. This implies some degree of uncertainty in the determination of the deformations defined by the wi, which must be removed by introducing a compatibility condition, based on the consideration that the strain field is non-rotational, so the curl of the deformations must satisfy the following equation: wy=x  wx=y ¼ T x wy þ T y wx ;

ð3Þ

where the equivalences Txr/y = Txwy and Tyr/x = Tywx have been exploited, and the vanishing of the curl of a gradient has been considered. When the wi are used as the generalized deformations of the plate, the constitutive characterization consists in determining the constitutive law that links the 12 internal force fluxes #i, respectful of the six equilibrium equations of Eq. (1), to the 12 deformations wi, under the six constraints stated by Eq. (3). It is worth stressing that all 12 fluxes and deformations are considered mainly for ease of notation, although it must be clearly understood that those related to the drill degree of freedom are totally unnecessary.

3. Kinematic characterization The displacement vector s of an arbitrary point P, whose reference position lies on the fibre in position P = {0, 0, z}T, can be decomposed in a ‘‘rigid’’ part v, due to the rigid displacement and rotation r of the sample fibre, and a ‘‘warping’’ part w, that accounts for the rest of the total displacement, namely: s ¼ v þ w: In the present work, the warping w is arbitrary, and will be later described by a general finite element discretization. The rigid part of the displacement is given by v ¼ Zr; where the matrix Z, describing a linearized translation and rotation, is Z ¼ ½I

ðP  OÞT ;

O = {0, 0, 0}T is the reference point of the fibre, namely the intersection of the fibre with the reference plane of the plate. The warping displacement w in general depends on a finite set of unknowns u, e.g. nodal displacements, if usual finite element shape functions are considered, through arbitrary trial functions N that are

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

defined along the thickness of the fibre, the ez-axis, oriented from the lower to the upper surface: wðx; y; zÞ ¼ N ðzÞuðx; yÞ: The dependence of the warping on the transverse position along the fibre is accounted for by the shape functions N(z), while the unknowns u(x, y) depend on the in-plane position in the reference surface. If the shape functions are able to describe a rigid displacement, the problem will be under-determined because the rigid displacement is already accounted for by the unknowns that describe the reference position and rotation of the fibre. These rigid modes related to the warping unknowns will require adequate constraints, and the true warping displacements will be determined by decoupling the resulting warpings from the rigid modes. The strains e = [exx, eyy, ezz, cyz, czx, cxy]T at an arbitrary point are related to the displacements s by means of the differential matrix operator o(Æ): 2 3 0 0 ð Þ=x 6 7 0 7 ð Þ=y 6 0 6 7 6 0 0 ð Þ=z 7 6 7 ›ð Þ ¼ 6 7: ð Þ=z ð Þ=y 7 6 0 6 7 6 ð Þ 0 ð Þ=x 7 4 =z 5 ð Þ=y

ð Þ=x

0

The operator o(Æ) can be divided in three parts, each accounting for the differentiation in one of the coordinate directions, ›ð Þ ¼ Xð Þ=x þ Yð Þ=y þ Zð Þ=z . They are defined as: 2 3 2 3 2 3 1 0 0 0 0 0 0 0 0 60 0 07 60 1 07 60 0 07 6 7 6 7 6 7 6 7 6 7 6 7 60 0 07 60 0 07 60 0 17 7; Y ¼ 6 7; Z ¼ 6 7 X¼6 60 0 07 60 0 17 6 0 1 0 7: 6 7 6 7 6 7 6 7 6 7 6 7 40 0 15 40 0 05 41 0 05 0 1 0

1 0 0

0 0 0

The strains due to the reference motion of the fibre are 8 9 vx=x þ z/y=x > > > > > > > > vy=y  z/x=y > > > > > > > > < = 0 : ð4Þ er ¼ > > vz=y  /x > > > > > > > > > > vz=x þ /y > > > > : ; vy=x  z/x=x þ vx=y þ z/y=y These terms are usually accounted for in the theory of thick, or Mindlin, plates (C0). If the transverse shear strain is neglected, i.e. /x = vz/y, /y = vz/x, it Kirchhoffs plate theory is obtained (C1). The transverse axial strain due to the reference motion is null: this means that the plate is assumed to be in a plane strain state; in other

1175

words, the plate is infinitely stiff thickness-wise. This assumption is unrealistic; it should be relaxed at least to the plane stress state. In a CLT-like formulation, the plane stress condition can be obtained by using plane stress two-dimensional properties for the materials of the laminæ. This relaxation is implicit in the proposed formulation thanks to the transverse warping terms. It is apparent from Eq. (4) that the rigid rotation about the transverse axis, z, does not participate in the straining of the plate. This introduces a drilling singularity in the stiffness matrix of the plate when all the six rigid degrees of freedom of the fibre reference point are taken into account. It is worth noting, from the last row of Eq. (4), that the mixed curvature terms and the in-plane shear derivatives are always paired in the strains, i.e. the constant and the linear terms of the in-plane shear strain, respectively vy/x + vx/y and /y/y  /x/x; This implies that their constituents are not independent, suggesting that extra singularities are likely to appear, as will be discussed later. The strains due to the rigid displacement of the fibre can be expressed as: er ¼ XZwx þ YZwy : The warping strains are eu ¼ XN u=x þ YN u=y þ ZN =z u: Matrix B ¼ ZN =z is made of the derivatives of the shape functions with respect to z. The complete strains are e ¼ XZwx þ YZwy þ XN u=x þ YN u=y þ Bu:

ð5Þ

The stresses at an arbitrary point are obtained from the strains by means of a constitutive law r ¼ rðe; zÞ

ð6Þ

that accounts for the anisotropy and non-homogeneity of the material. Eq. (6) allows to compute accurate stress predictions as soon as the strains are reconstructed from the solution in terms of linear and angular strains and warping unknowns.

4. External work The thickness-wise integral of the divergence of the work done by the stresses acting on the side surface of the fibre describes the transmission of the internal forces through the fibre itself. The external work results in: Z n Z   T    o o2 dLe ¼ ds px =x þ dsT py =y dz þ dsT f dz þ dsT F oh ; oxoy h h where px ¼ XT r; py ¼ YT r are the stresses that act on the fibre on the sides normal to the x and y axes, respectively. The work done by the distributed forces is temporarily neglected; it will be discussed in the section

1176

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

devoted to the solution of local problems. When the dependence of the virtual displacements ds on the discrete unknowns of the problem is exploited, the work becomes: Z   o2 dLe ¼ drT Z T px=x þ py=y dz oxoy hZ   þ duT N T px=x þ py=y dz Zh Z T T T þ dr=x Z px dz þ du=x N T px dz Zh Zh T T T Z py dz þ du=y N T py dz; þ dr=y h

h

or, by rearranging the terms, and recalling that the linear rigid displacement matrix Z and the shape functions N do not depend on the in-plane coordinates: Z  Z  ! o2 dLe T T T ¼ dr Z px dz þ Z py dz oxoy h h =x =y Z  Z  ! þ duT

þ

drT=x

N T px dz

h

Z

T

Z px dz þ h

Z

N T py dz

þ

=x

drT=y

Z

Z

h

=y

T

duT=x

Z py dz þ

h

N T px dz þ duT=y N T py dz: h R R After defining #i ¼ h Z T pi dz; Pi ¼ h N T pi dz, it is worth noticing that #i represents the forces and moments resulting from the stress distribution on the side of the fibre in the ith direction, namely the flux of internal forces and moments in that direction, while Pi is the work done by the same stresses for the warping displacement du. If the warping unknowns u are nodal displacements, the Pi are the corresponding generalized nodal forces. By recalling the definition of the generalized deformations of the plate in Eq. (2), the gradient of the reference displacements becomes: 

h

r=x ¼ wx  T x r;

r=y ¼ wy  T y r:

The final form of the external work is   o2 dLe ¼ drT #x=x þ #y=y  T Tx #x  T Ty #y þ dwTx #x oxoy   þ dwTy #y þ duT Px=x þ Py=y þ duT=x Px þ duT=y Py :

deT r dz: h

The local stresses depend on the local strains by means of a three-dimensional linear constitutive law C(E) = C(E)(z), and pre-stresses r0 and imposed strains e0 can be accounted for as follows: r ¼ r0 þ C ðEÞ ðe þ e0 Þ: The effect of a pre-stress r0 is important in the evaluation of the bending properties of a plate subjected to membrane loads, while the imposed strains e0 can be useful in evaluating the behavior of thermally loaded plates. In the following sections, these terms are not considered, but they can be included in a straightforward manner. In fact, they contribute to the right-hand side of the problem in a way that resembles that of the external loads, so they can have the same transverse and inplane distribution properties that will be addressed in a subsequent section dedicated to local problems. When both stresses and strains are expressed in terms of the linear and angular strains and of the nodal warping displacements, the internal work is 8 9T 2 u=x > Mxx Mxy Cx Lxx > > > > > > > 6 > u=y > Myy Cy Lyx > > 6 > > < = 6 2 o dLi 6 E Rx ¼d u 6 > > 6 oxoy > > > > 6 > > w Axx > x> > 4 > > : > ; wy sym:

h

Cx ¼

Z

N T XT C ðEÞ B dz;

Cy ¼

h

Lxy ¼

Z Z

N T XT C ðEÞ XZ dz; h

N T XT C ðEÞ YZ dz;

Z

N T YT C ðEÞ XZ dz;

h

Lyy ¼

Z

N T YT C ðEÞ YZ dz;

h

The strain work inside the fibre is represented by the thickness-wise integral of the work done by the stresses for a virtual variation of the strains:

E¼

Z

B T C ðEÞ B dz; h

Z h

h

Lyx ¼ 5. Internal work

38 9 Lxy > u=x > > > 7> > > > Lyy 7> > u=y > > < > = 7> 7 u Ry 7 ; 7> > > > > > > > w Axy 7 5> x > > > > : > ; wy Ayy

where the submatrices are Z N T XT C ðEÞ XN dz; Mxx ¼ h Z Mxy ¼ N T XT C ðEÞ YN dz; h Z Myy ¼ N T YT C ðEÞ YN dz;

Lxx ¼ It is worth noticing that the expression that multiplies the virtual variation of the reference displacements r is the differential equilibrium equation, Eq. (1), that was initially written based on equilibrium considerations.

Z

o2 dLi ¼ oxoy

N T YT C ðEÞ B dz;

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

Rx ¼

Z

BT C ðEÞ XZ dz;

Ry ¼

h

Axx ¼ Axy ¼ Ayy ¼

Z

Z

2

BT C ðEÞ YZ dz;

6 6 4

h

Cy  CTy

Z T XT C ðEÞ XZ dz;

Zh Zh

skw: 2

Z T XT C ðEÞ YZ dz; T

T

Z Y C

ðEÞ



E

6 þ4

Rx Axx

sym:

YZ dz:

1177

38 9 Lyy > > u > > 7< = 7 wx 0 0 5> > > :w > ; y 0 =y 8 9 8 9 3 u Ry > 0 > > < > = > < > = 7 w Axy 5 ¼ #x : x > > > > > : > ; A :w ; # Lyx

yy

y

ð8Þ

y

h

In case of typical material properties, i.e. thickness-wise constant, the integration of matrices Aij can be performed analytically, leading to a set of coefficients that resemble those of the CLT; actually they would be the CLT coefficients if plane stress material coefficients were used. The other matrices will generally require numerical integration, although a closed form solution may be found in the simplest cases.

6. Plate fibre problem By equating the internal and the external work, the problem becomes: 9 2 38 9 8 u=x > Mxx Mxy Cx Lxx Lxy > Px > > > > > > > > > > > > > > 6 7 > > > > u Myy Cy Lyx Lyy 7> Py > > 6 =y > < < = = 6 7 6 7 u E R R P þ P ¼ : x y x=x y=y 6 7> > > > > 6 7> > wx > > > > > sym: Axx Axy 5> # > > > 4 x > > > > > > : > : ; > ; wy Ayy #y ð7Þ The indefinite equilibrium equation should be added for completeness, but it will be implicitly satisfied during the solution phase by an appropriate choice of the internal forces #i. The divergence of the first and the second block rows of the system allows to eliminate the unknown vectors Pi, whose divergence appears at the right-hand side of the third block equation. The system becomes partial differential, with respect to x and y, first order in the generalized deformations wi, and second order in the discrete warping unknowns u: 2 38 9 2 38 9 Mxx 0 0 > Myy 0 0 > < u > < u > = = 6 7 6 7 4 4 0 0 5 wx 0 0 5 wx > > : > : > ; ; wy =xx wy =yy sym: 0 sym: 0  38 9 2 u > 0 0 > Mxy þ MTxy 7< = 6 7 6 4 0 0 5> w x > : ; wy =xy sym: 0 8 9  2 3 Lxx Lxy > Cx  CTx < u > = 6 7 4 0 0 5 wx > : > ; wy =x skw: 0

In terms of differential/algebraic equations (DAE), it is an index 1 problem, that can be cast in partial differential form (PDE) by computing wx and wy from the last two block equations, and by substituting them and their partial derivatives in the first block equation. This operation is not recommended, though, because it destroys the sparsity of the problem. The right-hand side contains the internal force fluxes in the in-plane directions. The eigensolution of Eq. (8) is self-balanced by definition. It leads to the so-called ‘‘extremity problem’’, which allows one to solve a plane-wise boundary problem in terms of stresses and warpings. The extremity problem will be addressed in further research; in this paper only the solution related to non-null internal force fluxes, which has been called ‘‘central solution’’, will be addressed. The nomenclature for the solutions first introduced by Giavotto et al. in 1983 in the work on the characterization of composite beam sections [7] is here preserved.

7. Right-hand term The central solution requires the determination of indicial solutions for the internal fluxes. They result from impulsive loads and satisfy the indefinite equilibrium of the plate, but they are not unique, so they are determined by imposing compatibility conditions to the corresponding strain and displacement fields. The fluxes are assumed to take the following form, depending on their unknown value H0x, H0y at the reference point (i.e. x = y = 0):     #x ¼ I þ xT Tx H0x ; #y ¼ I þ yT Ty H0y : ð9Þ These are beam-like internal force fluxes that implicitly satisfy the indefinite equilibrium of Eq. (1); Fig. 1 illustrates the transverse shear case. There are other combinations of linear internal force and moment fluxes that satisfy Eq. (1) and at the same time are self-balanced, i.e. their value is null in (0, 0); they can be combined to those in Eq. (9) to enforce compatibility. The compatibility problem basically arises from the transverse shear solution, since it involves a linear internal moment flux to balance the constant transverse shear flux. Two self-balanced in-plane internal bending moments and two self-balanced in-plane internal forces are added.

1178

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

z

z

Fzx Fxy/y

Fyx/y y x

M yx/x

y x Fxx/x

Fig. 1. Main transverse shear on the face normal to the x-axis of the fibre.

Fig. 2. Self-balanced in-plane forces in the x direction.

Formally, three moment and three force fluxes are added, but those related to the drill degree of freedom will be eliminated later. They are " #sx ¼

ez 

0

0

ez 

" þ " #sy ¼

# x 0

0

I þ ex  ex  0

0

I þ ey  ey  # !

"

# !

I þ ey  ey 

I þ ex  ex 

þ

z

ez 

0

0

ez 

y

y H0s ;

# x

y H0s ;

x

M yy/y M xx/y

M yx/x

Fig. 3. Self-balanced in-plane bending moments in the x direction.

corresponding to: 2

0 x

6 6x 6 6 60 #sx ¼ 6 60 6 6 60 4 0 2 x 6 6 y 6 6 6 0 #sy ¼ 6 6 0 6 6 6 0 4 0

0 0

0

y

0 0

0

0

0 0

0

0

0

y

x

0

0

x

0

0

0 0

0

0

0

0

0 0

0

0

0 0

0

0

0 0

0

y

y

0 0 y

x

0 0

0

0

0

3

7 07 7 7 07 7H0s ; 07 7 7 07 5 0 3 0 7 07 7 7 07 7H0s : 07 7 7 07 5 0

These internal forces implicitly satisfy Eq. (1) and, as they have no constant part, their value is zero in the reference point (0, 0). The self-balanced force terms are considered first. In the second column, for instance, the normal force flux on the face normal to the x-axis,

in direction x and varying linearly with x, Fxx/x, is balanced by a shear force flux on the face normal to the y-axis, in direction x and varying linearly with y, Fxy/y. A corresponding shear force flux on the face normal to the x-axis, in direction y and varying linearly with y, Fyx/y, is required to satisfy the symmetry of the in-plane shear stress as shown in Fig. 2. The same considerations hold for the moment terms. The two self-balanced bending moment fluxes are also null at (0, 0). In the case of column 4, the bending moment flux on the face normal to the x-axis of the fibre, in direction y and varying linearly with x, Myx/x, is reacted by a moment flux on the side normal to the y-axis, in direction y and varying linearly with y, Myy/y. A corresponding moment flux on the side normal to the x-axis, in direction x and varying linearly with y, Mxx/y, is required to satisfy the symmetry of the shear stress, as shown in Fig. 3.

8. Solution The complete right-hand side of the characterization problem,

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190 2 3 3 0 0 I 0 6 7 6 7 6 0 6 xex  I 7 07 7H0y 7H0x þ 6 ¼6 6 I 6 0 07 07 4 5 4 5 yey  I 0 0   2 3 0 xez  þy I þ ey  ey  6 7 6 0 xez  þy ðI þ ex  ex Þ 7 7H0s þ6 6 xðI þ e  e Þ  ye  7 0 4 5 x x z   0 x I þ ey  ey   yez  2

(

#x #y

)

suggests that a linear solution in x and y of Eq. (8) must be sought. By rearranging the constant and the linear terms of the internal force fluxes, the following expression is obtained:   #x ¼ #0x þ x#x1x þ y#y1x H 0 ;   #y ¼ #0y þ x#x1y þ y#y1y H 0 ; where H 0 ¼ ½H0x ; H0y ; H0s T . In the same fashion, the displacement and strain unknowns are written in linear form, depending on the influence coefficients, as follows:   wx ¼ w0x þ xwx1x þ ywy1x H 0 ; ð10Þ   wy ¼ w0y þ xwx1y þ ywy1y H 0 ;

ð11Þ

  u ¼ u0 þ xux1 þ yuy1 H 0 :

ð12Þ

The problem parts, which problem in x 2 E Rx 6 A 4 xx sym:

is decomposed in its linear and constant are solved consecutively. The first step yields: 32 x 3 2 3 u1 0 Ry 76 x 7 6 x 7 Axy 54 w1x 5 ¼ 4 #1x 5; Ayy

wx1y

#x1y

where the right-hand side is composed by the x-dependent part of the assumed internal forces. It should be noted that matrices Tx and Ty, which represent the linear parts of the right-hand term, have only two non-null coefficients, and thus only two solutions at each first step, plus the four non-null terms out of the six self-balanced load terms added for compatibility enforcement, are actually required. The problem is written in matrix notation for the sake of notation clarity. Moreover, one of the non-null coefficients of each of the matrices Tx, Ty represents the contribution of the in-plane shear to the equilibrium of the drilling moment, which can be discarded since, as shown later, that degree of freedom will need to be constrained. Similarly, the first step problem in y yields: 2 32 y 3 2 3 0 u1 E Rx Ry 6 76 y 7 6 y 7 Axx Axy 54 w1x 5 ¼ 4 #1x 5: 4 #y1y wy1y sym: Ayy

1179

The second step problem requires the two previous solutions to build the right-hand side: 2 32 u 3 2 3 E Rx Ry 0 0 6 7 6 7 w 7 6 7 Axx Axy 56 4 4 0x 5 ¼ 4 #0x 5 w0y sym: Ayy #0y 3 3 2 2  ux1 Cx  CTx Lxx Lxy 76 x 7 6 7 6 þ6 0 0 7 54 w1x 5 4 x w1y skw: 0  2 2 3 3 uy1 Cy  CTy Lyx Lyy 6 76 y 7 7 76 þ6 0 0 54 w1x 5: 4 skw:

0

wy1y

The strains and the stresses in the general case can be obtained in a straightforward manner from these solutions as functions of the imposed fluxes of internal forces and moments; the former can be obtained from their definition, Eq. (5), and the latter from the former, by means of the constitutive law. It is worth remarking that the same matrix must be solved for all three problems. It is intrinsically singular, due to nine different displacement mechanisms. The first five singularities are related to the warping shape functions, that allow the three rigid translations and two rotations, respectively about the axes orthogonal to the fibre, x and y. Since the rigid displacements of the fibre are already accounted for by ad hoc rigid unknowns, the redundant rigid displacements introduced by the warping must be eliminated. This result can be achieved with little effort by noticing that a statically determined constraint is required. Because the domain of the fibre is a straight line, it is sufficient to constrain one node in the z direction, and two distinct nodes in both the x and y directions. The procedure can be automated, e.g. by constraining the top and bottom nodes in the in-plane directions, and an arbitrarily chosen node in the transverse direction. There are four other singularities related to the rigid degrees of freedom. The first two are represented by the two components of the gradient of the rotation about the z-axis, /z/x and /z/y. Since the drill singularity is intrinsic in the one-dimensional nature of the fibre model, it can be fixed by eliminating that degree of freedom, or, to preserve the structure of the generalized deformations with six degrees of freedom, the singularity can be avoided by means of a penalty function, namely an arbitrary coefficient on the diagonal of matrix A corresponding to the two above mentioned degrees of freedom. There are two more singularities related to the definition of the in-plane shear strain in the generalized deformation vectors, as previously anticipated. Since the symmetry of the strains has not been imposed, the kinematic model still contains two rigid rotation terms about the z-axis one is constant,

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

given by the two contributions of the reference displacement to the in-plane strain, assuming the value vx/y = vy/x, which corresponds to a rigid rotation about the z-axis; the other is thickness-wise linear, given by the twist rotations about the in-plane axes, when they assume the value /x/x = /y/y. A synthetic way to eliminate the latter two singularities, which preserves the structure of the generalized deformation unknowns, is presented later.

9. Compatibility enforcement The compatibility of the solution obtained so far must be enforced. The generalized deformations due to constant internal forces intrinsically satisfy the compatibility equation, Eq. (3). Only the deformations due to the transverse shear loads still need some care. These deformations violate the compatibility only in the inplane directions for both the forces and the moments, due to the in-plane linearly varying moments that are required to satisfy the constant transverse shear equilibrium equations. The generalized deformations due to the (nominally) 12 beam-like plus the (nominally) 6 self-balanced internal force fluxes that result from the solution of the two-step problem will be condensed to 12 solutions related to (nominally) 12 independent and compatible unit internal force and moment fluxes. The generalized deformations in Eqs. (10) and (11) are partitioned as follows: ( ) " # wxe0 þ xwxxe1 þ ywyxe1 wxs0 þ xwxxs1 þ ywyxs1 wx ¼ wye0 þ xwxye1 þ ywyye1 wys0 þ xwxys1 þ ywyys1 wy ! " H0 ; ð13Þ  H0s where the subscripts e and s respectively stand for ‘‘equilibrium’’ and ‘‘self-balanced’’, while the H0 are the 12 indicial internal force and moment fluxes. After applying Eq. (3) to the resulting generalized deformations of Eq. (13), a linear compatibility problem in x and y is obtained, resulting in three separate sub-problems:   wxy1e  wyx1e þ T x wy0e  T y wx0e H0   þ wxy1s  wyx1s þ T x wy0s  T y wx0s H0s ¼ 0; ð14Þ 

   T x wxy1e  T y wxx1e H0 þ T x wxy1s  T y wxx1s H0s ¼ 0;

while the first, Eq. (14), can be solved to determine the six H0s as functions of the twelve imposed unit internal forces H0, yielding:  y H0s ¼  wxy1s  wyx1s þ T x wy0s  T y wx0s    wxy1e  wyx1e þ T x wy0e  T y wx0e H0 : ð17Þ The matrix in Eq. (17) cannot be inverted because two of the self-balanced solutions refer to the drill degree of freedom, so the pseudo-inversion, denoted by the (Æ) operator, is required. Fig. 4 shows how the stresses of the non-compatible transverse shear solution combine with those of the non-null self-balanced solutions to yield the corresponding compatible transverse shear solution in a 4 ply (45/90/0/45) laminate. Notice that the plots related to the non-compatible and to the compatible solution enclose the same area, while those related to the self-balanced solutions enclose a null area. The solutions related to the warping and the generalized deformation must be corrected to account for the compatibility solution, to yield the (nominally) twelve compatible solutions that describe the behavior of the sample fibre for unit internal forces. This is accomplished by substituting the multipliers H 0 with H 0 "

¼ 

wxy1s  wyx1s þ T x wy0s  T y wx0s

I y 

wxy1e  wyx1e þ T x wy0e  T y wx0e

#  H; 0

where Eq. (17) has been used. The reference displacement can be obtained from the compatible solution as follows: a solution that is quadratic in the in-plane coordinates is expected, namely 1 1 2 y r ¼ r0 þ xrx1 þ yry1 þ x2 rx2 þ xyrxy 2 þ y r2 : 2 2

stress yz h/2

non-compatible transverse shear self-bal. x moment self-bal. y moment self-bal. x force self-bal. y force

thickness

1180

compatible transverse shear

0

ð15Þ 

   T x wyy1e  T y wyx1e H0 þ T x wyy1s  T y wyx1s H0s ¼ 0:

-h/2

ð16Þ The last two problems, Eqs. (15) and (16), vanish due to the intrinsic structure of the solution, i.e. Tiwj1 = 0,

-1

0

1

2

3

4

5

6

stress Fig. 4. Transverse shear stresses in a 0/90/45/45 laminate: five non-compatible solutions combine to yield the compatible one.

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

By recalling the definitions of the generalized strains, one obtains:   y x wx ¼ rx1 þ xrx2 þ yrxy 2 þ T x r0 þ xr1 þ yr1 þ   y x wy ¼ ry1 þ yry2 þ xrxy 2 þ T y r0 þ xr1 þ yr1 þ ; where the second order terms have been neglected, resulting in the constraints T i rx2 ¼ 0;

T i rxy 2 ¼ 0;

T i ry2 ¼ 0

for i = x, y. By collecting the constant and the linear terms, and by considering that the constant displacement r0 can be arbitrarily set equal to zero without any loss in generality, the following relations can be obtained: rx1 ¼ wx0 ; ry1 ¼ wy0 ; rx2 ¼ wxx1  T x wx0 ; ry2 ¼ wyy1  T y wy0 ; y rxy 2 ¼ wx1  T x wy0 ;

¼ wxy1  T y wx0 : Notice that the last expression results in rxy 2 being overdetermined by the two definitions of wx and wy; this should not be surprising since the compatibility equation, Eq. (3), directly states that the two expressions found for rxy 2 must be equal for the solution to be compatible.

10. Characterization of the sample fibre The elastic properties of the fibre can be obtained indifferently by the definition of the internal or of the external work, by substituting all the unknowns with the previously determined solutions that depend on the imposed fluxes of internal forces H0, where the face indices x, y have been dropped for simplicity, namely: dHT0 Fðx; yÞH0 ¼

o2 dL : oxoy

Matrix F is the compliance matrix of the fibre. Because the plate is assumed to have constant properties regardless of the position in the reference plane, the compliance matrix must be independent of x, y, as can be easily proven. The inverse constitutive relation of the fibre is w ¼ F#:

ð18Þ

The direct constitutive relation is obtained by inverting the compliance matrix: Kw ¼ #; which yields the stiffness matrix K. It is worth remembering that if all the 12 generalized strain degrees of freedom have been considered, the 12 · 12 compliance

1181

matrix F is four times singular due to the missing imposition of the symmetry of the strains tensor and to the drill degrees of freedom, and thus cannot be inverted. A pseudo-inversion can be used, but here a degree of freedom condensation that leads to a more intuitive definition of the generalized deformations and internal force fluxes is favored. Let the condensed generalized ~ ¼ fex ; ey ; c ; c ; deformations vector be defined by w xy xz cyz ; jx ; jy ; jxy g that includes the in-plane axial and the transverse shear strains unchanged from the previous definition, the in-plane shear strain being defined as the sum of the two parts that were used in the wi vectors. The bending curvatures are maintained, together with a twist cross-coupling term made of the difference between the derivatives of the in-plane twist rotations, as shown in the transformation matrix S: 2 3 0 0 0 0 0 0 0 0 0 1 0 0 60 0 0 0 0 0 0 1 0 0 0 07 6 7 6 7 60 1 0 0 0 0 1 0 0 0 0 07 6 7 6 7 60 0 1 0 0 0 0 0 0 0 0 07 6 7: S¼6 0 0 0 0 0 1 0 0 07 60 0 0 7 6 7 60 0 0 7 1 0 0 0 0 0 0 0 0 6 7 6 7 40 0 0 0 1 0 0 0 0 0 0 05 0 0 0 1 0 0 0 0 0 0 1 0 The pseudo-inverse transformation is 2 0 0 0 0 0 0 0 0 6 0 0 1=2 0 0 0 0 0 6 6 60 0 0 1 0 0 0 0 6 6 6 0 0 0 0 0 0 0 1=2 6 6 0 0 0 0 0 0 1 0 6 6 6 0 0 0 0 0 0 0 0 6 Sy ¼ 6 6 0 0 1=2 0 0 0 0 0 6 6 0 0 0 0 0 60 1 0 6 60 0 0 0 1 0 0 0 6 6 6 0 0 0 0 0 1 0 0 6 6 4 0 0 0 0 0 0 0 1=2 0 0 0 0 0 0 0 0

3 7 7 7 7 7 7 7 7 7 7 7 7 7 7; 7 7 7 7 7 7 7 7 7 7 7 5

which can be interpreted as the inverse of S in an energetic sense, that satisfies the virtual external work bal~ ¼ dwT #, with w ¼ fwT ; wT gT and # ¼ f#T ; ~ T# ance dw x y x T T ~ ¼ f/x ; /y ; /xy ; /xz ; /yz ; lx ; ly ; lxy g. The #y g , where w condensation is based on the consideration that the inplane strains at an arbitrary point due to the generalized deformations have the expression cxy ¼ ðvy=x þ vx=y Þ þ zð/y=y  /x=x Þ;

ð19Þ

as results from the last row of Eq. (4). The two strain terms in Eq. (19) represent the symmetric part of the

1182

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

gradient of the reference position; this is consistent with the consideration that the anti-symmetric part of the gradient describes a rigid rotation and thus does not produce strain. As a consequence, these terms represent a measure of the in-plane shear strain and of the in-plane twisting. The conjugated condensed force and moment fluxes ~ ¼ frx ; ry ; sxy ; szx ; szy ; lx ; ly ; lxy g. The following are # expressions apply: ~ ¼ Sw; w ¼ S y w; ~ w ~ ¼ S yT #; #

~ # ¼ S T #:

By means of the above reported transformations, matrices Mij, Ci, Lij, Ri and Aij can be easily reduced to the new set of generalized deformations, to intrinsically eliminate the singularities related to the rigid rotations. This reduction can be performed even before the system is solved, to obtain an invertible matrix, provided the singularities related to the warping are also eliminated. By substituting the reduced deformations and internal fluxes in Eq. (18), ~ ¼ SFS T #; ~ w ~ ¼ SFS T is obtained. the reduced compliance matrix F ~ ¼F ~ 1 , and the comThe reduced stiffness matrix is K yT ~ y plete stiffness matrix is K ¼ S KS . It is worth recalling that the stiffness matrix could not have been directly computed by inverting the compliance matrix, because they are structurally singular. Notice that, if matrices S, S are defined in such a manner that the transverse shear strains cxz, cyz, i.e. rows 4 and 5 of matrix S and columns 4 and 5 of matrix S, are eliminated, the Kirchhoff plate stiffness properties can be directly obtained.

11. Warping decoupling As previously noticed, the warping can contain some rigid displacement modes: the three rigid displacements and the two rigid rotations about the in-plane axes. In fact, the statically determined constraints of the warping unknowns that are required to eliminate the singularities of the matrix, and the arbitrary choice of the reference point along the fibre, usually imply that the resulting displacements contain some rigid motion. The warping w can be written as the sum of a rigid term wr and of some true, unknown warping wd, namely w = wr + wd, provided a definition of the decoupling criterion is given. A rigid displacement v is given by v = Zr as a function of the six rigid displacements and rotations of the reference point. The nodal displacements that correspond to r can be written as follows: ur ¼ Z node r ¼ Hr; where matrix H is made of a stack of matrices Z evaluated at the nodes. Then the warping at an arbitrary

point, due to a rigid motion, can be written as wr = N(z)Hr. In fact, the shape functions assume null value at every node, except for the one they refer to, where they take unit value. Provided the shape functions are able to describe a constant and a linear z-wise displacement, the rigid rotation-translation matrix can be written as Z(z) = N(z)H. The warping displacements can be uncoupled from the rigid ones in a least squares sense as follows: H T ðu  HrÞ ¼ 0: The problem can be separated in two linear parts, respectively in x and y, and in a constant part, namely:  1  1 rx1 ¼ H T H H T ux1 ; ry1 ¼ H T H H T uy1 ;  1 r0 ¼ H T H H T u0 : A more sound decoupling can be obtained in an energetic sense, following the idea that was introduced by Borri and Merlini for the warping of a beam section [9], which consists in imposing that the decoupled warping be orthogonal to the generalized nodal forces: Z   wTd px þ py dz ¼ 0; h

which, after substituting wd = w  Zr, and after discretizing, results in:     uT Px þ Py ¼ rT #x þ #y ; ð20Þ where the property #i = HTPi has been exploited. The resulting quantity is quadratic in x and y, since all the quantities in Eq. (20) are linear, so only the linear part will be decoupled. The constant and the x, y linear problems are  T  T r0 ¼ #x0 þ #y0 Px0 þ Py0 u0 ;   T  T rx1 ¼ #x0 þ #y0  #xx1 þ #xy1 r0   T   x x x þ Px1 þ Py1 u0 þ Px0 þ Py0 u1 ;   T  T  #yx1 þ #yy1 r0 ry1 ¼ #x0 þ #y0   T   þ Pyx1 þ Pyy1 u0 þ Px0 þ Py0 uy1 : In both cases, the same matrix needs to be inverted for the three problems. Notice that, in the second case, the problem is over-determined, because 6 (actually 5) rigid displacements are required to decouple the warpings due to 12 (actually 8) load conditions. In this case the warping and the nodal forces are not exactly orthogonal, but their product is constant in a least squares sense. The decoupled warpings are obtained by subtracting the rigid decoupling part from the coupled solution: uxd1 ¼ ux1  Hrx1 ;

uyd1 ¼ uy1  Hry1 ;

ud0 ¼ u0  Hr0 :

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

Attention must be paid to the generalized deformations, since they depend on the rigid displacements and their derivatives; as a consequence, the derivatives of the rigid part of the warping that has been subtracted from the solution must contribute to the generalized deformations. Recalling their definition, the uncoupled generalized deformations become: wxd1 wd0



 Tx x Tx y r1 ; wyd1 ¼ wy1 þ r; ¼ wx1 þ Ty Ty 1  x  Tx r r0 : ¼ w0 þ 1y þ Ty r1

It is worth noting that the characterization of the sample fibre, i.e. the determination of the flexibility and the stiffness matrices, does not depend on the mentioned decoupling, which needs to be performed only to correctly visualize the deformed shapes. Even the recovery of strains and stresses does not depend on the decoupling. Also the decoupling problem should be condensed by using matrices S, S, because it is singular, as described for the main problem.

12. Distributed loads The external work done by distributed loads can depend on the in-plane and transverse (thickness-wise) position in an arbitrary fashion. Such work, done by volume and surface loads, is given by: o2 dLd ¼ oxoy

Z h

 dsT f dz þ dsT F oh ;

which, after discretization, becomes: o2 dLd ¼ drT oxoy

Z h

T

þ du

  Z T f dz þ Z T F oh

Z h



   N f dz þ N F oh : T



T

1183

8 > > > > > > <

9 Px > > > > > Py > = RHS (7) ¼ Px=x þ Py=y þ F u > > > > > > #x > > > > > > : ; #y and that of Eq. (8) becomes: 8 9 > < Fu > = RHS (8) ¼ #x : > : > ; #y From now on, the distributed loads f, F are assumed to have the form F(x, y, z) = V(z)U(x, y), or, in other words, the thickness-wise variable z is separable from the planewise ones, x and y; if this is not the case, no particular integral can be computed. The indefinite equilibrium equation, Eq. (1), becomes: #x=x  T Tx #x þ #y=y  T Ty #y þ F r ¼ 0: The load vector can be written as: F r ¼ F r Uðx; yÞ; where U(x, y) is a scalar function depending only on the in-plane position, while F r is a constant vector (a ‘‘load mode’’), i.e. Fr when U = 1. Among the cases of distributed loading of interest, consider the effects of distributed inertia forces due to rigid body motion acceleration field, or the effects of constant pressure on either or both sides of the plate. Specialized solutions for these cases can be easily found. A general solution for distributed loads can be found only at the cost of considering a combination of recursive solutions, resulting in a rather complicated formulation that will be investigated in future works. The case of constant pressure is addressed here; given a load function U(x, y) = p0, with thickness-wise load vector F r ¼ feTz ; 0T gT , the following function is a solution in terms of internal force and moment fluxes:    p 1 2 #i ¼  0 x þ y 2 T Ti þ iI F r : 2 2

The dependence of matrix Z and of N on the transverse position is considered in the integration, as well as that of the loads. The boundary of the domain h, indicated as oh, is represented by the intersection of the lower and upper surfaces of the plate with the axis of the fibre, z. In such points, a node should be defined, so the surface forces F (pressures) act as external nodal loads on the boundary nodes. The work can be written as:

This is not the only solution for U(x, y) = p0; it is the one that propagates isotropically from the local origin, i.e. the point that is being characterized. This solution does not account for any specific boundary conditions.

o2 dLd ¼ drT F r þ duT F u ; oxoy R R where F r ¼ h Z T f dz þ Z T Fjoh and F u ¼ h N T f dz þ T N T Fjoh with Fr = H Fu, H being the rigid nodal displacement matrix. The right-hand side of Eq. (7) becomes:

The extension to the analysis of plates incorporating layers of piezoelectric materials is relatively straightforward. In practical cases, the topology of the piezoelectric inclusions is very simple, consisting in one or more layers of piezoelectric material intermixed with conventional layers. The electric potential in the piezoelectric

13. Piezoelectric plate analysis

1184

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

material is usually applied by means of conductive laminæ, annealed on the faces of the active layers, and the active and the host layers are separated by thin films of epoxy glue. As a consequence, the most significant component of the electric potential is thickness-wise, so there is little interest in highlighting the plane-wise component. In analogy with the warping part of the displacement, the electric potential V(x, y, z) is discretized by means of the shape functions introduced in previous sections, yielding V ¼ NðzÞdðx; yÞ; and the electric field, derived from the potential by means of the gradient operator $(Æ) is E ¼ $V ¼ ex Nd =x  ey Nd =y  ez N =z d: After discretization, the nodal unknowns related to the electric potential are formally treated in the same manner of those related to the warping displacement; as a consequence, the solution procedure is not altered when considering the piezoelectric problem. However, the electric potential unknowns related to the conductive laminæ are collected in the set of generalized degrees of freedom of the plate. The linear piezoelectric constitutive law is !

r D

"

" ¼

C ðEÞ e

eT
ðeÞ

#!

e E

" :

The internal work becomes Z  T  o dLi de r  dE T D dz; ¼ oxoy h 2

where the minus before the dielectric contribution is used to preserve the overall symmetry of the problem; this operation is allowed by the arbitrariness of the virtual variations d(Æ). The external work must also be expanded to account for the dielectric contributions; formally: Z !    T    o2 dLe ¼ ds px =x þ dsT py =y  dV eTx D =x oxoy h Z   "  dz þ dsT f dz þ dsT F   dV eT D y



X

=y

h

oh

dV i qi

8 > > > > > > <

9 Px > > > > > Py > = RHS (7) ¼ Px=x þ Py=y þ F u þ Q > > > > > > #x > > > > > > : ; #y and that of Eq. (8) becomes: 9 8 > = < Fu þ Q > #x RHS (8) ¼ > > ; : #y Note that, since the electric charges qi are applied on conductors, their electric potential must be constant, so in general Q does not depend on the plane-wise position. One consequence of this, which results in a modification of the first step of the solution procedure for the plate fibre characterization, is that the plane-wise derivatives of the electric potential on the conductors must be grounded in a sort of non-holonomic constraint, to make the electric potential on the conductors constant; this has been clearly discussed in [10]. Unless the dielectric properties of the overall laminate, and possibly those of the surrounding medium, are considered, each piezoelectric lamina represents an isolated dielectric domain. In any case, exactly one electric potential degree of freedom needs to be grounded for each dielectric domain; this physically expresses the fact that the difference between the electric potential at the two conductors of each piezoelectric patch is retained as generalized electric unknown. The generality of the approach allows to consider special cases, in which electrostatic circuitry components are used to connect the conductors of the piezoelectric devices; e.g. shunts, capacitors and more.

14. Numerical results The case of an isotropic, homogeneous thick plate is considered first. The analytical characterization of such a structural component is quite easy, so the problem is used to draw some initial considerations on the performances of the different discretization models. When isoparametric shape functions are used, in two, three and four-node elements, the following considerations are valid:

i

The last term on the right-hand side is the electric work made to charge the conductors that bound the piezoelectric inclusions. Notice that formally it resembles the work done by the external pressure on the faces of the laminate. As a consequence, after discretization, an additional dielectric nodal force term appears; the right-hand side of Eq. (7) becomes:

• The two-node element is able to describe a constant thickness-wise strain and stress condition. This means that, as the number of elements increases, the two node element can asymptotically reach the correct solution, but its behavior is quite poor. Nonetheless, a single element is able to determine the membrane properties of a plate.

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190 0.5 0.4 0.3

Thickness position (%thk)

• The three-node element is able to model a constant and a linear thickness-wise strain and stress condition. It is thus able to describe both the membrane and the bending behavior of a plate, and the transverse contraction due to the bending load, yielding the correct stiffness properties for the bending case. • The four-node element is able to model up to a parabolic thickness-wise strain and stress condition. As a consequence, it describes the shear strain and stress condition related to a transverse shear load, as given by the approximate shear theory. A single four-node element is thus able to determine all the stiffness characteristics of an isotropic plate.

1185

Present FEM Layered Elm. FEM 2D FEM 3D

0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.06

-0.04

-0.02

0

0.02

0.04

0.06

Sy

Fig. 6. 0/90 lateral (y) stress compared to FEM solutions.

0.5 Present

0.4

Thickness position [% thk]

The propagation of the effects of local boundary conditions can usually be described by means of transcendent functions [11]. These effects can be captured only asymptotically by the piece-wise polynomial interpolation performed by the finite elements. The use of thickness-wise finite elements instead of regular shape functions spanning the entire thickness, as accomplished by other methods, allows to comply with interlaminar continuity requirements. In fact, compatibility and equilibrium issues respectively require the in-plane strains and the transverse stresses to be continuous, while the transverse strains and the in-plane stresses can be discontinuous. Provided inter-element boundaries fall at each interlaminar boundary, the above described discontinuities can occur freely, as illustrated by Figs. 5–12, discussed later.

0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5

14.1. Isotropic homogeneous plate

0

0.5

1

1.5

2

2.5

Txz

Fig. 7. 0/90 transverse (xz) shear stress.

The stiffness properties of a plate made of isotropic homogeneous material have been calculated, both analytically and numerically, by means of different 1

Thickness position [% thk]

0.75 0.5

Present FEM Layered Elm. FEM 2D FEM 3D

0.4

Thickness position (%thk)

0.3 0.2 0.1 0

0.25

Sx Sy Txy

0 -0.25 -0.5 -0.75

-0.1 -0.2

-1 -0.8

-0.3

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

Stress

-0.4 -0.5 -3

0.5

Fig. 8. Five-layer laminate unit axial (x) force flux. -2

-1

0

1

2

3

4

5

6

Sx

Fig. 5. 0/90 longitudinal (x) stress compared to FEM solutions.

order elements in different mesh patterns. The plate is made of Aluminum Alloy, with elastic modulus

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190 1

1

0.75

0.75

Thickness position [% thk]

Thickness position [% thk]

1186

0.5 0.25

Sx Sy Txy

0 -0.25 -0.5

-1 -0.5 -0.4 -0.3 -0.2 -0.1

0

0.1

0.2

0.3

0.4

0.5

0 -0.25 -0.5

1 Sx Sy Txy

0.75 0.5 0.25 0 -0.25 -0.5 -0.75 0

0.1

0.2

0.3

0.4

0.5

-1 -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

Fig. 12. Five-layer laminate unit transverse (yz) force flux.

Fig. 9. Five-layer laminate unit in-plane (xy) force flux.

Thickness position [% thk]

0.25

Stress

Stress

0.6

0.7

0.8

0.9

1

Stress

Fig. 10. Five-layer laminate unit axial (y) force flux.

1 Txz Tyz

0.75

Thickness position [% thk]

0.5

-0.75

-0.75

-1 -0.1

Txz Tyz

0.5 0.25

while in the analytical computations the plane stress material properties have been used, in the numerical analysis the three-dimensional properties of the material have been considered. In fact, the proposed formulation implicitly accounts for the transverse deformation of the plate due to the imposition of plane stress boundary conditions. The comparisons are reported in Tables 1– 3. The single two-node element model gives the same shear coefficient in both the in-plane and transverse directions. In fact, this element is able to model only constant strain, so the parabolic shear strain due to transverse shear is interpolated by means of a constant strain evaluated at midpoint, giving the same result as the in-plane strain, which is actually constant and thus exact. The single three-node element model gives the exact values for the bending coefficients, while the single four-node element model gives the exact value for all the coefficients, as expected. Good convergence can be appreciated for the eight two-node element model in modeling the bending coefficients, as well as for the eight three-node element model in approximating the shear coefficients. 14.2. Anisotropic laminated plate

0 -0.25 -0.5 -0.75 -1 -0.2

0

0.2

0.4

0.6

0.8

1

1.2

Stress

Fig. 11. Five-layer laminate unit transverse (xz) force flux.

E = 74.5 GPa and Poisson modulus m = 0.30, with a h = 2.0 · 103 m thickness. It is important to note that,

The properties of a plate composed of a stack of laminæ of anisotropic material with different orientations have been analysed. The material presented by Pagano has been considered, with properties E1/E2 = 25, G12/ E2 = 0.5, G23/E2 = 0.25, m12 = 0.25 [11]. Comparison is made with the CLT, with constant stress in each lamina. The constant strain effect is achieved by means of twonode elements for each lamina. Accurate results, obtained by means of four-node elements, are presented in order to underline the improvement both in the determination of the stiffness properties and in the recovery of stresses and strains. Figs. 5 and 6 show the longitudinal and the lateral stresses due to unit axial (x) load and

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

1187

Table 1 Membrane stiffness Analytical (N/m)

2 n, 1 el. Ratio

3 n, 1 el. Ratio

4 n, 1 el. Ratio

8

1.63736 · 10

1.000

1.000

1.000

mE K 1;2 ¼ 1m 2 h

4.91209 · 107

1.000

1.000

1.000

E h K 3;3 ¼ 2ð1þmÞ

5.73077 · 107

1.000

1.000

1.000

K 1;1 ¼ K 2;2 ¼

E 1m2

h

Table 2 Bending stiffness Analytical (N/m) K 6;6 ¼ K 7;7 ¼

E h3 1m2 12

3

mE h K 6;7 ¼ 1m 2 12 3

E h K 8;8 ¼ 2ð1þmÞ 12

2 n, 1 el. Ratio

2 n, 8 el. Ratio

3 n, 1 el. Ratio

4 n, 1 el. Ratio

5.45788 · 10

1.225

1.004

1.000

1.000

1.63736 · 101

1.750

1.012

1.000

1.000

1.91026 · 101

1.000

1.000

1.000

1.000

1

Table 3 Shear stiffness Analytical (N/m) K 4;4 ¼ K 5;5 ¼

E 5 2ð1þmÞ 6 h

7

4.77564 · 10

2 n, 1 el. Ratio

3 n, 1 el. Ratio

3 n, 8 el. Ratio

4 n, 1 el. Ratio

1.200

1.197

1.000

1.000

unit bending (y) moment. The FEM layered element solution has been obtained with a model composed of ABAQUS 3D layered elements, while the 2D solution has been obtained by using both NASTRAN and ABAQUS 2D layered elements (the results match exactly), and the 3D solution has been obtained by stacking ABAQUS 3D elements in a 3D model. Finally, Fig. 7 shows the transverse shear stress due to unit transverse force flux. Four third-order elements per lamina have been used to improve the quality of the plot, but the same results can be obtained with one single element. Consider now a more sophisticated laminate made of five layers of the previously considered unidirectional material, with the stacking orientation of +45/45/90/ 45/+45, and with thicknesses 0.25/0.25/1.0/0.25/0.25, resulting in a laminate that is representative of a composite rotor blade structure. The stresses through the

thickness for unit axial and transverse force fluxes are presented in Figs. 8–12. Notice how the in-plane axial stress is carried by the mid-part of the laminate only in the (y) case (Fig. 10), in which the mid-part is oriented as the load; in the other cases (Figs. 8 and 9), there is a significant participation of the shearing of the ±45 laminæ. The transverse stresses are dominated by the cross-coupling induced by the ±45 laminæ, while the compatibility requirement causes the unidirectional mid-part to give a very limited contribution. 14.3. Comparison with NASTRAN PCOMP analysis This section presents a comparison between the equivalent material properties computed for NASTRAN finite element analysis code as resulting from its internal lamination procedure and from the proposed

Fig. 13. NASTRAN PCOMP card for 0/90 laminate referred to bottom surface, with corresponding MAT2 material properties compared with present ones (ANPA).

1188

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

Fig. 14. NASTRAN PCOMP card for ±45/0/90 laminate referred to mid surface, with corresponding MAT2 material properties compared with present ones (ANPA).

Table 4 Pressure transducer: charge per unit surface

Charge

[mm] 2

[C/mm ]

0.25

1.0

2.0

5.0

10.0

0.574e3

0.422e3

0.388e3

0.365e3

0.356e3

formulation. Three different case are considered: 0/90, ±45 and ±45/0/90 laminates, referred both to the bottom and to the mid surface, to assess the invariance with respect to the reference surface (Figs. 13 and 14). As expected, the equivalent properties match quite well in all cases except for those involving the transverse shear. The current analysis results are labeled as ANPA (Anisotropic Non-homogeneous Plate Analysis). 14.4. Piezoelectric pressure transducer A pressure transducer has been analyzed by modeling a piezoelectric lamina, made of PZT5H, bonded on the top surface of an aluminum alloy host structure, with elastic modulus E = 70e3 N/mm2 and m = 0.33. The transducer is loaded, on the side of the piezoelectric, with uniform pressure p0. The resulting polarization charge per unit surface can be computed by multiplying the constitutive matrix of the smart plate by the generalized deformations and the electric potential on the conductors. The charge resulting from different host/piezo thickness ratios is reported in Table 4; the thickness of the piezo is 0.25 mm, and a 0.01 mm thick glue layer (E = 3e3 N/mm2, m = 0.3) has been considered as well. Fig. 15 shows the normal in-plane and thickness-wise stresses of the model with 1 mm thick host structure; the electric voltage is also shown, which varies almost linearly through the thickness of the piezoelectric layer of the plate. This example shows the versatility of the proposed formulation in characterizing as well as in analysing local problems for arbitrarily loaded, anisotropic composite plates incorporating piezoelectric layers.

0.6

0.4

Thickness, mm

Thickness

0.2

0

-0.2

-0.4

-0.6

Sx (0., 0.) Sz (0., 0.) Txz (1., 0.) V (0., 0.) -1

-0.5

0

0.5

Stress, N/mm^2; Electric Potential x 10^5, V

Fig. 15. Pressure transducer stresses and electric potential due to unit pressure. Table 5 ATR material properties Property

Glass–epoxy

Graphite–epoxy AFC

c11 c12 c22 c66 q

14.8e3 1.3e3 13.7e3 1.9e3 1.8e3

181.8e3 2.9e3 10.3e3 7.2e3 1.6e3

83.5e3 14.8e3 42.4e3 11.4e3 1.8e3

[MPa] [MPa] [MPa] [MPa] [kg/m3]

Thickness h

0.130

0.125

0.140

[mm]

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

1189

Table 6 ATR blade active skin stiffness properties (units are N, C and mm)

ux uy uxy uxz uyz lx ly lxy Qp

ex

ey

cxy

cxz

cyz

jx

jy

jxy

2.57e5 2.29e4 0.00e0 0.00e0 0.00e0 9.34e3 1.36e5 1.08e3

2.29e4 6.33e4 0.00e0 0.00e0 0.00e0 1.93e4 9.34e3 1.08e3

0.00e0 0.00e0 2.57e4 0.00e0 0.00e0 1.08e3 1.08e3 9.35e3

0.00e0 0.00e0 0.00e0 7.62e3 8.65e0 0.00e0 0.00e0 0.00e0

0.00e0 0.00e0 0.00e0 8.65e0 9.20e3 0.00e0 0.00e0 0.00e0

9.34e3 1.93e4 1.08e3 0.00e0 0.00e0 4.38e4 1.48e4 1.36e3

1.36e5 9.34e3 1.08e3 0.00e0 0.00e0 1.48e4 1.92e5 1.36e3

1.08e3 1.08e3 9.35e3 0.00e0 0.00e0 1.36e3 1.36e3 1.66e4

0.00e0

0.00e0

6.50e1

0.00e0

0.00e0

7.36e0

7.36e0

4.08e1

14.5. Active fibre composite rotor blade Active fibre composite (AFC) materials represent a very interesting and innovative material for adaptive structures [12]. Their application in distributed strain induced actuation of helicopter rotor blades has been investigated at model scale since the mid 90s at MIT and at the NASA Langley Research Center in the Active Twist Rotor (ATR) concept [13]. An analytical model has been investigated in [14] by means of the piezoelectric beam characterization formulation presented in [10]. The same model is here characterized as a plate fibre. The laminate is made of a substrate of unidirectional carbon fibre on the inner side, with 4 layers of ±45 fibreglass/+45 AFC/±45 fibreglass/45 AFC, all wrapped by bottom and top ±45 fibreglass plies. The AFC is laid on the front part of a ‘‘D’’-shaped

web made of T-300 unidirectional graphite-epoxy. Data is taken from [15]; it refers to a generic medium weight helicopter rotor that is representative of a broad class of helicopters. The properties of the materials are reported in Table 5; the AFC properties result from the finite element analysis of a specimen [14]. The AFC uses interdigitated electrodes to exploit the larger direct e33 piezoelectric coupling coefficients, instead of the transverse e31 coefficients. Since the proposed formulation only considers actuation by means of thickness-wise electric field, an equivalent material has been considered such that a transverse electric field generates the corresponding amount of direct and transverse strain in the active fibres. The fibre stiffness matrix resulting from the proposed analysis method is reported in Table 6. Notice the high piezoelectric coupling of the in-plane shear and twist. Fig. 16 shows the internal stresses in the active laminate related to unit pressure applied on the top face.

1.5

15. Concluding remarks 1

Thickness, mm

0.5

0

-0.5

-1

-1.5 -5

Sx (0., 0.) Sy (0., 0.) Sz (0., 0.) Txz (10., 0.) -4

-3

-2

-1

0

1

2

3

4

5

Stress, N/mm^2

Fig. 16. Active fibre composite: stresses due to unit pressure.

A formulation for the semi-analytical characterization of generally anisotropic, non-homogeneous plates incorporating piezoelectric inclusions has been presented. A deep insight into the kinematics and the equilibrium of a plate fibre has been given, to justify the assumptions made in developing the model and to illustrate its capability of capturing the inmost aspects of inter-laminar equilibrium and compatibility. A thorough characterization of the overall properties of the plate fibre has been presented, and the detailed stress recovery has been addressed. The formulation is best suited for analyzing the elastic properties of sophisticated laminated plates, which can be subsequently modeled in a finite element manner in conventional analysis codes. The results of such analyses, in terms of internal force and moment fluxes, can be used to recover the laminar stresses and the strains with a high level of detail. The analogies with widely used finite element codes laminate analysis capabilities, and the possible improvements that

1190

P. Masarati, G.L. Ghiringhelli / Computers and Structures 83 (2005) 1171–1190

the formulation can give, have been discussed. The extension to the modeling of piezoelectric laminæ included in the lamination sequence has been discussed, and results concerning advanced applications have been presented.

References [1] Krishna Murty AV, Vellaichamy S. On higher order shear deformation theory of laminated composite panels. Compos Struct 1987;8:247–70. [2] Kant T, Pandya BN. A simple finite element formulation of a higher-order theory for unsymmetrically laminated composite plates. Compos Struct 1988;9:215–46. [3] Barboni R, Gaudenzi P. A class of C0 finite elements for the static and dynamic analysis of laminated plates. Comput Struct 1992;44(6):1169–78. [4] Di Sciuva M. Multilayered anisotropic plate models with continuous interlaminar stresses. Compos Struct 1992;22(3):149–67. [5] Carrera E. C0 Reissner–Mindlin multilayered plate elements including zig-zag and interlaminar stress continuity. Int J Numer Meth Eng 1996;39(11):1797–820. [6] Narasimha Reddy J. Theory and analysis of elastic plates. Philadelphia: Taylor & Francis; 1999. [7] Giavotto V, Borri M, Mantegazza P, Ghiringhelli GL, Caramaschi V, Maffioli GC et al. Anisotropic beam

[8]

[9] [10]

[11]

[12]

[13]

[14]

[15]

theory and applications. Comput Struct 1983;16(1–4): 403–13. Ghiringhelli GL, Mantegazza P. Linear straight and untwisted anisotropic beam section properties from solid finite elements. Compos Eng 1994;4(12):1225–39. Borri M, Merlini T. A large displacement formulation for anisotropic beam analysis. Meccanica 1986;21:30–7. Ghiringhelli GL, Masarati P, Mantegazza P. Characterisation of anisotropic, non-homogeneous beam sections with embedded piezo-electric materials. J Intell Mater Syst Struct 1997;8(10):842–58. Pagano NJ. Exact solutions for rectangular bidirectional composites and sandwich plates. J Compos Mater 1970;4:20–34. Bent AA, Hagood NW, Rodgers JP. Anisotropic actuation with piezoelectric fiber composites. J Intell Mater Syst Struct 1995;6:338–49. Wilkie WK, Wilbur ML, Mirick PH, Cesnik CES, Shin SJ. Aeroelastic analysis of the NASA/Army/MIT active twist rotor. In: American Helicopter Society 55th Forum, Montreal, Canada, May 25–27 1999. Ghiringhelli GL, Masarati P, Mantegazza P. Analysis of an actively twisted rotor by multibody global modelling. Compos Struct 2001;52/1:113–22. Wilkie WK, Park KC, Belvin WK. Helicopter dynamic stall suppression using piezoelectric active fiber composite rotor blades. In: AIAA/ASME/AHS structures, structural dynamics and materials conference, Long Beach, CA, April 20–23, 1998. AIAA-98-2002.