- Email: [email protected]
On the linear three-dimensional behaviour of composite beams
Composites Part B 2,8B (1997) 613-626
N;" I
PII: S1359-8368(96)00075-3
ELSEVIER
© 1997 Elsevier Science Limited Printed in Great Britain. All rights reserved 1359-8368/97/$17.00
On the linear three-dimensional behaviour of composite beams
G.L. Ghiringhelli Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, via C. Golgi, 40, 20133 Milano, Italy (Received 11 September 1995; accepted 20 November 1996) This paper describes a method of evaluating the linear stiffness matrix of a three-dimensional beam which accounts for all the possible couplings involved in the use of composite laminations appropriate to satisfy design requirements. The starting point of the procedure is represented by a complete evaluation of the beam crosssection properties by use of a two-dimensional finite element analysis. A mass matrix of the beam is also evaluated as well as the external forces, to account for overall thermal deflections. Some examples are presented and the results are compared either with their theoretical counterparts or with numerical results obtained by means of a full three-dimensional finite element analysis. © 1997 Elsevier Science Limited. (Keywords: B. mechanical properties; B. numerical modelling; composite beams)
INTRODUCTION Notwithstanding the availability of tools of ever-increasing power and reliability, the beam continues to be a model used widely in many branches of structural analysis, such as in the design of helicopter or aeolian generator blades, or the study of aircraft dynamics. Indeed, the advantages of an approach based on this model are well known. The development of tools able to predict the complex behaviour of composite beams has been driven by the increasing usage of these materials; there are several reasons for this, such as the significant improvement in fatigue life and damage tolerance, and the emphasis set on the possibilities offered from modem techniques of design, e.g. aeroelastic tailoring. In this context particular relevance is given to correct prediction of the elastic couplings between different load conditions. They can be introduced intentionally but can also represent an unwanted consequence of the design of composite laminations. Yet another problem is associated with both thermal load conditions and the anisotropic thermal characteristics of the materials, which can be responsible for global beam deflections. For example, in extreme environmental situations, owing to non-balanced thermal deformation, an aeroelastic tailored helicopter blade could show overall deflections that interfere with the designed structural couplings; the effort of the designer can then be nullified, possibly leading to a serious deterioration in performance. Since closed-form solutions are not general enough, the structural designer/analyst has to resort to numerical approximations which, when carried out on three-dimensional
finite element models, can be too costly. This is evident during the initial design phase, when the structural lay-out is not yet finalised and many analyses have to be carried out. The importance of the subject is confirmed by the available literature. In 1990 Hodges ~ carried out an exhaustive survey of the state-of-the-art on composite rotor blade modelling, citing almost 30 works. He pointed out that two approaches are mainly adopted, the first one being an analytical development of the problem while the second is based on some kind of finite element formulation. With regard to modelling generic, non-homogeneous and/or anisotropic beams showing any kind of thermoelastic coupling, the first category does not seem to be general enough in face of the design needs for composite beams and in particular for rotor blades. A further classification of finite-element-based methods can be made in terms of the in-plane and out-of-plane deformation of the beam cross-section. The availability of both these strain modes is strictly connected to a correct analysis of a solid beam. Several works allow out-of-plane cross-section warpings only facing the problem of thin-walled anisotropic beams; e.g. Giavotto et al. 2, Banchau 3 and Bauchau and Hong4 all adopted a two-dimensional (2D) approach to characterise the cross-section, while Stample and Lee 5 preferred a threedimensional (3D) formulation with warping nodes placed over the cross-section to describe bending, torsional and warping deformation. Similarly, Cofie and Bank 6 developed a 3D two-node beam element capable of accounting for warping and anisotropic coupling in thin-walled composite beams. These are suitable for a direct analysis of the stiffness matrix. Also papers by Smith and Chopra 7, Barbero et al. 8
613
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli and Bicos and Springer 18 suffer from some kind of limitation in blade topology. Other works finally remove this limitation: e.g. Giavotto et al. 9 and Kosmatka 1° both employ a similar finite element approach, and Hodges et al. 11 who instead exploit a variational-asymptotical beam sectional analysis. The work of Kosmatka is devoted to curved blades while the method of Hodges et al. is suitable for the analysis of straight, slender beams. Giavotto et al. consider straight beams that not are necessarily slender, taking into account a full cross-section characterisation. In fact, this work presents a completely general formulation for the evaluation of the cross-section elastic and mass properties and the stress distribution in a non-homogeneous and/or anisotropic beam. Irrespective of the geometrical shape of the crosssection (i.e. solid or hollow), thick or thin walls, single or multiple cells and type of materials adopted, the de Saint Venant approach is extended to straight, untapered and untwisted composite beams. The problem is solved by using a 2D finite element model of the cross-section without any undue simplification for geometry and material properties; the elastic characterisation of the cross-section is complete and the shear stiffnesses and related couplings are also furnished by the method. In this way great computational efficiency is gained during all analysis phases. Following the same basic idea of a semi-discretisation, the formulation was rewritten to extend the approach to curved and twisted beams by Borri et al. 12. A final extension was presented by the author t3, the basic formulation being completed to take into account strains due to the thermal anisotropy of the material and/or a complex temperature distribution on the beam cross-section. The present paper is a natural completion of these activities: it resumes previous works and describes a procedure to evaluate the stiffness and mass matrices of a two-node, three-dimensional beaml4; this element, which is suitable for use within the framework of a linear analysis, shows all the possible couplings introduced by the use of appropriate composite laminations. The starting point of the procedure is the characterisation of the beam cross-section performed with the cited 2D finite element analysis 9"15, the main points of which are briefly recalled. A theoretical section is then devoted to evaluation of the thermal equivalent loads that can be applied to a beam element to obtain a given thermal strain. The results obtained with the described procedures are compared either with their theoretical counterparts or with numerical results obtained by a full 3D finite element analysis. A final example compares the experimental frequencies of a helicopter rotor blade with results from an eigenvalue analysis on a beam model.
Tx2~_~~ M2
7j"
X3
M 3~
T1
XI
Figure 1 Beam reference frame and section load convention
to Giavotto et al. 9 and Ghiringhelli and Mantegazza 15 for an exhaustive discussion of the formulation. The classical theory of elasticity defines a beam as a prismatic, slender and straight body having both geometric and elastic properties that are constant along each generatrix. A section is obtained by cutting the beam along a plane normal to its axis in the undeformed configuration. According to the de Saint Venant hypotheses, we assume that loads applied along the axis produce negligible effects with respect to those due to extremity loads. Furthermore, the effects of their distribution at the beam extremities vanish far from the end at a distance that is small with respect to the beam length. To define the deformation of each cross-section we assume small displacements and strains. A linear behaviour of material mechanical properties is also assumed. Arbitrarily distributed anisotropy or non-homogeneity can characterise the cross-section but, according to the previous hypothesis, they must be constant along the beam axis. If these hypotheses are matched the elastic problem for a beam is greatly simplified, being reduced to a monodimensional problem. So the stress distribution on each cross-section can be analysed independently from each other, once the distribution of internal forces is known. It must be noted that linearity assumptions, in both displacements-strain and stress-strain laws, do not prevent a nonlinear analysis of a body exhibiting large displacements but small strains; in fact, the linearity assumptions refer only to the local behaviour of a section of the beam t6. The beam geometry, depicted in Figure 1, is defined in an orthogonal Cartesian reference frame, {x} [XlX2X3] T, with the axis x3 parallel to the beam generatrix. Let us define the strain and stress columns in the following form: =
{/3} ~- [/311 /322 2/321 2/331 2/332 /333] T
{ o / = I°,1 °22 o21 °3, ANISOTROPIC BEAM THEORY In this section some basic elements of the formulation of the elastic problem for a non-homogeneous/anisotropic beam are presented. They are fundamental for the arguments that will be presented in the following text. The reader is referred
614
so that the stress components acting on a section, i.e. tY3i, i ---1,3, can be easily isolated by formally extracting the second half of the stress vector with an appropriate matrix operator: {p} = [s]{~}
where [S] T = [0 I3], [I3] being a third-order identity matrix.
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
x2l TI,2
M T
Ml:~
~ X3
a
1,2
~
x,
Tl,i
T2.1
Figure 2 Three-dimensional elemental frame
In a beam model the displacement of a generic point is described by means of rigid displacements and rotations of the origin of the reference frame, while the stress state is reduced to resultants of internal forces and moments. If {T} and {M} are the arrays of resultant forces and moments, respectively, the conventions adopted being depicted in Figure 2, we can define the cross-section load vector as {0} T = [{T}T,{M}T]. It can be evaluated from: {0} =
~A[Z]T{p} dA
(1)
i.e. by integrating the elementary forces and moments due to stresses acting on the cross-section. The matrix [Z] is defined in such a way that the cross-product matrix related to the point with coordinate {x} is included:
[z] = [I (XA)T] The displacements, {X}, and rotations, {~o}, of a reference point belonging to the corresponding cross-section can be organised into a single array, {r}: {r} = [{x}TI{so}T] T When significant warping is present, i.e. in composite beams, the displacement of a generic point of the crosssection is profitably described as the sum of a term depending on the rigid movement of the cross-section and a warping term straining it, {g}: {s} = {g} + [Z]{r}
(2)
The sixth-order matrix [O] has only two non-null terms according to eqn (4):
[O]=
,
0
0
0
{g(xl,x2,x3) } = [Ni(xl,x2)] {u(xil,xi2,x3) }
(4)
(5)
eqn (5) states that nodal warpings are discrete parameters on the cross-section, still depending on the position along the beam axis; i.e. a semi-discretisation of warping is used. The following corresponding approximations for {s } and {e} are obtained: {s} = [Nl{u} + [Z]{r} {el = [B]{u} +
[S][NI{uI/3[S][ZI{d/}
(6) (7)
where [B] = [©][N] and
(3)
It takes into account the strain due to differential rigid displacement of adjacent cross-sections. The subscript /3 indicates differentiation with respect to the x3 direction; this notation will be used in the following developments.
1
The first two elements of {~b}, i.e. ~b, = Xl/3 - ~o2 and ~2 = X2/3 + ~°1, represent the transverse differential displacement, i.e. shear strain of the cross-section; the third term, i.e. ~b3 = X3z3, is the axis elongation; while the three last elements are associated with bending, ffl = ¢t/3 and ~b2 = ~o2/3, and torsion, ~b3 = ~o3/3, of the reference line, respectively. After discretising the beam section with appropriate twodimensional finite elements (Figure 2), we approximate the warping function {g} by interpolating the nodal values {u} with usual shape functions [N] with:
By writing strains as differentiation of displacements allows us to introduce the cross-section strain array, {~b}: {f} = {r}/3 q- [O]{r}
[t]=
0
"/1
0
0-
0
/2
0
12
/1
0
0
0
/1
0
0
/2
0
0
0
[D] =
615
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli The virtual work principle per unit length of the beam is written as
fA(~{s}T{p})/3dA=
IA~{8}T{0r} dA
(8)
By substituting eqns (6) and (7) into eqn (8) and carrying out the explicit differentiation of the term appearing in the lefthand side, we can write, recalling eqn (1):
{u/3/TI il{u/3/ LTCcT u { U/3}T/i 3 + {rrlT({o}13 -- [oIT{0})
(9)
where the following definitions are introduced: [M] =
JA[N]T[s]T[D][S][N]dA
{0} = ([I] +X3[O]T}{00}
(12)
where {00} is the resultant array at the origin. We then try to find a particular solution for {u}, {P} and {~b} having an equivalent linear form, i.e.
~A[B]T[D][B]dA
[El =
[C] =
(10)
depends on the way the beam ends are constrained and is generally negligible after a length of the order of the largest linear dimension of the section itself. Beyond this length at both the beam extremities, i.e. in its central portion, only the particular integral, which will be seen to depend only upon {0}, is significant. This is in accordance with Saint Venant's principle and consequently the general and particular integral were called by Giavotto et al.9 the extremity and central solution, respectively. The central solution is thus related to the determination of the usual global elastic properties and stresses of a beam section. The reader can refer to 9,a4for a comprehensive explanation of the complete solution of equations (11). The central solutions must satisfy eqn (1 ld) for generic extremity loads, so we directly look for linear internal resultants {0} of the type:
{u}= ([.o] +x3[ulJ){0o}
~A[B]T[D][S][N]dA
(13)
{P}= ([P0] +x3[P1]){0o} [R] =
~A[B]T[D][S][Z]dA
[L] =
~A[N]T[s]T[D][S][Z]dA
[A] =
fA[Z]T[s]T[D][S][Z] dA
{P} =
{P}/3 =
=
where [Uo], [ud, [Po], [P1], [~bo] and [~bl] are solution matrices; each of their six columns represents the solution for the corresponding unitary term of the related constant and linear part of {0 } respectively. After substituting the above equations into equations (11) and equating the terms having the same power of x3, we end up with the following system of linear equations with unknowns [u0], [ui], [Po], [P,], [¢'o] and [~bl]:
fA[N]T{p} dA fA[N]T{p}/3 dA
These matrices are evaluated successfully by using isoparametric elements and then assembled by the usual procedure in a finite element code. For arbitrarily varying {u}, {U/3}, {~} and {r}, eqn (9) leads to the following linear system of algebraic equations in {~b} and ordinary, first-order differential equations in {u }, {P} and {0}:
(14)
[pl] = tc]T[u,] + u0
-- C
(15) (16)
{P} = [M]{u}/3 + [c]T{u} -{- [L]{~}
(1 la)
[Po] = [M] [ul] + [c]T[uo] + [L] [~bo]
[C]{u}/3 -t- [E]{u} + [R] {~b}
(1 lb)
[L]T{u}/3-q-[R]T{u} -4- [A]{~b}
(1 lc)
The solution of the above sets of equations is naturally carried out by solving eqns (14), (15), (16) and (17). The right-hand side of eqn (14) reduces to only two columns because of the structure of [®] [see eqn (4)]; it must be noted that eqns (14) and (16) have the same coefficient matrix, thus a single factorisation is required for their solution. This is carded out after removing, in any convenient way, the six singularities due to the unknown warping {u } which, because of the properties of the finite element shape functions, already contains the rigid motions explicitly accounted for by {r}. Because the warping {g} contains
{P}/3 = {0} =
{0}/3 = [®]T{0}
(lld)
It can be seen that the equilibrium equations of a slice of infinitesimal length, i.e. Eqn (1 ld), are obtained. The solution of the above system is composed of a general and a particular integral. For a sufficiently long and slender beam, the general integral is related to the diffusion of stresses from the terminal sections of the beam itself. It
616
+x3
(17)
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli arbitrary rigid contributions owing to the applied constraints, a procedure to remove them is needed, as explained previously 14,16.
u/3/aNlS/j
additional term to eqn (9), the expression of the virtual work principle:
/
CROSS-SECTION STIFFNESS/COMPLIANCE MATRIX To recover the elastic properties of the beam we note that the central solution can be written by using eqn (13), so that the external virtual work of the beam slice, i.e. the first term on the right-hand side of eqn (9), can be written in a very concise form:
~{00}T([~0] T q-[u0]T [PI] "1=[t/lIT [P0]) {00} =6{00}T[F]{00}
(18)
The analysis of this equation allows us to define the compliance matrix of the section in the following form:
[F]=[J/o]TW[uo]T[P,]+[ul]T[po]
(19)
This can be shown to be symmetric by taking into account the first of eqn (16), eqn (17) and the expression for the identity matrix [/] arising from the second of eqn (16). Furthermore, it is independent of any rigid motion of the section, as demonstrated by Ghiringhelli and Mantegazza 15. In particular, it does not depend on the way in which the partition of eqn (2) is performed; in fact, constraints applied to nodal warpings to remove singularities from eqns (14) and (16) affect only displacements while stresses and the compliance matrix are not influenced. This agrees with the definition of that matrix as an intrinsic characterisation of the overall elastic properties of the beam. Conventional section parameters, used in standard engineering beam theory, can be obtained from the compliance matrix of the section. If the [F] matrix is diagonal or if couplings are restricted to axial/bending and/ or shear/torsion load conditions, engineering parameters, such as stiffness, shear centre and normal stress centre locations, can be evaluated so that the stiffness matrix of a beam is built with known techniques. But when more complex couplings are present, e.g. in the case of composite beams, the matrix could be fully populated. Some terms can be physically interpreted: by means of coefficients (4,6) and (5,6) the orientation of the torque moment that does not couple with bending moments can be defined as:
=
~
fA[ZIT[s]T[D]{¢T}dA
The temperature behaviour along the beam axis can be modelled with a polynomial expansion across the analysed cross-section, but higher order contributions are not of interest because, according to the Saint Venant principle, local effects due to distributed loads should be negligible. So a linear development is assumed, i.e. using nodal temperatures {t} and nodal temperature gradient {t'} vectors, and the thermal strain is expressed by: {eT} = {~}[N]T({t} +x3{t' })
(21)
where {o~} is the vector form of the thermal expansion tensor, i.e.
331T
=
In eqn (21) the distribution on the cross-section of the temperature and its axial gradient are obtained by the solution of the thermal problem 13. By using this equation the following matrices of nodal thermal influence coefficients can be defined:
IHT] fa [N]T [S]T[D]{~ }[N]T dA =
~A[B]T[DI{ot}[N]TdA
IVy] =
=
(22)
fA[Z]T[s]T[D]
{or}[U] T dA
and eqn (20) becomes:
~£Therm=
aU ~¢
[gl({t} q-x3{/'}) [W] ({1} +x3{t'})
(23)
Coherently with eqn (21), we still look for a linear {u}, {~b} and {P}:
° t x l = t a n - l ( -F46~F44j' ° t x 2 = t a n - l ( -F56~F55j A general procedure is needed to obtain a model showing the same behaviour of the original beam; this is the goal of a subsequent section.
THERMAL STRAIN The presence of thermal strain can determine the global deflection of a beam and this problem has great importance because when some constraints are applied, a severe tension state can be introduced. Thermal strain results in an
{u} = {u 0} + x 3 {u, }
(24a)
{~b} = {~b0} + x3{ ~bl }
(24b)
{P} = {Po} +x3{P, }
(24c)
In the absence of external loads, applying the same procedure as used previously leading to eqns (14)-(17), we obtain the following form of these equations:
Ul
617
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
{P1 } = [c]W{ul } + [L] {¢~ } -- [H] {t'}
U0
-- C
(25b)
the forces acting at the two ends of the beam, i.e. {0~} and {02}, are related by the following equation:
u1
{0, } = [H]{02 }
(28)
+ [:]{t}
(25C)
In addition, remembering that {0} is relative to cross-sections having the outgoing normal in the direction of x3, the nodal forces of the three-dimensional element can be written as:
{Po}=[M]{u,} +[c]T{u0} +[L]{¢o}-[H]{t}
(25d)
{S}~'{ -{01} }{02}
(29)
{s} = [L] {o2}
(30)
where the unknowns {u0}, {ul}, {P0}, {P1}, {¢0} and {¢1} are related to a specific temperature distribution given by {t} and {t' }. While eqns (14)-(17) are solved only once for unit external loads, in combining these solutions to obtain actual loads, the analysis of several thermal load conditions needs forwards/backwards substitution processes to be repeated for each set of temperature distributions. A restart procedure can be used successfully in this case. It can be appreciated that the assumption of linear behaviour along the axis enables us to keep the way of solving eqns (14)-(17) with the only modification being evaluation of the right-hand term. Nonetheless, higher terms could be accounted in a very simple way leading to further computational steps of the same kind as in eqns (25c) and (25d).
THE EQUIVALENT THREE-DIMENSIONAL BEAM
or, by eqn (28):
wherein the following matrix has been introduced: [LI=
From previous considerations, the internal forces in a genetic point, i.e. ~', can be interpolated from the end values, and from eqn (29) we obtain: {0} = [M]{S}
(26)
where load and displacement arrays, {S} and {W} respectively, are defined with the usual assumptions These terms are defined in the elemental frame so that the equilibrium equations expressed by the stiffness matrix must be rewritten in the basic frame by means of an appropriate transformation. To develop the formulation let us consider a generic beam, of length l, whose cross-section is characterised by the compliance matrix [F]. A natural local frame is conveniently assumed, the origin of which is placed at the middle of the beam. A normalised abscissa ~"is used for convenience, it ranges from - 1 to 1 when spanning from the first to the second node. According to beam theory, the equilibrium equations of a slice of infinitesimal length are
(Figure2).
{0}/3=[o]T{0} {0}/33=0
[M]=I~(~-l)[16] 1(~+1)[16]J Remembering equations (25), the virtual variation of the internal complementary work of the beam can be written as:
Thus the internal forces vary linearly along the beam and their end values are correlated; in fact, by integrating eqn (27) and defining a matrix [H] in the following form: [H]=
618
l[t]
[I1
6Li= ~lfl _ 1{60}y{F]{O}d~
(32)
while the term related to external forces is given by: 6Le = {rs}T{ W}
(33)
By introducing eqn (31) into eqn (32) and integrating we obtain the following form of the internal complementary work: ~Li = {c5S}T[p] {S }
(34)
where: I I2[F] [F]=6[-[F]
-[F]] 2[F] J
The complementary work equation is written by using eqn (30) to obtain an equation showing independent virtual variations only, i.e. Because of the arbitrariness of these virtual variations we obtain:
{602}.
{02} = ([L]T[pl[L]) -I[L]T{w} (27)
(31)
with
The aim of this section is to formalise a procedure to evaluate the 12th-order stiffness matrix of a simple, linear, 3D beam element with two nodes, in the form: {S} = [K]{ W}
[I] J
(35)
Finally, premultiplying by [L] and remembering eqn (30), we obtain eqn (26) where the stiffness matrix of the beam is defined as: [K] = [L] ([L]r [p] [L]) - 1[L]T
(36)
It must be stressed that this matrix accounts for all the global effects of the cross-section strain; so there is no need to
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli append other generalised degrees of freedom to nodal displacements accounting for cross-section warpings. During the evaluation of elastic matrices, i.e. Eqn (10), the mass matrix of the cross-section [m] can also be evaluated by a simple summation of inertia products, until the second order, integrated over each element. As far as a beam slice is concerned, once the terms due to the second time derivatives of warpings are neglected, the inertial forces per unit length related to the motion of the cross-section can be approximated by the inertial forces associated with rigid motion of the cross-section only: {FI} = [m]{}:}
with prescribed boundary conditions on the lateral surface. The effects due to the actual boundary conditions on the extremity cross-sections are neglected. This approach furnishes a temperature distribution throughout the beam volume: the behaviour is linear along the beam generatrix while the distribution on the cross-section depends on the thermal properties. This temperature distribution defines the initial thermal strain which leads to the thermal load defined by eqn (23) and the elastic problem is solved by equations (25). The displacement solution is still obtained in linear terms of strain cross-section parameters and nodal warpings. In particular, in order to describe the global deflection of the beam, we are interested in the first element expressed in the form of eqn (24). Once the section strains given by the solution of equations (25) are known, eqn (3) can be easily integrated to obtain the displacement components of the reference line at any point; i.e. in case of constant cross-section strain parameters it can be verified that:
(37)
In order to allow dynamic analyses these distributed loads must be concentrated at the beam nodes: a consistent inertial load array can be evaluated by introducing the interpolation of the rigid displacement of the reference frame at any position along the beam axis. Appropriate shape functions [N] can be applied to nodal displacements and rotations:
{ r} = [U]{ W}
(38) {r}=
According to a possible organisation of nodal unknowns and forces arrays, the [IV] operator can assume the form:
['o
{rl} + x 3 I
['o
{¢1}
(41)
where the subscript 1 indicates the boundary conditions at
"N 1
0
0
0
N2
0
N3
0
N1
0
N2
0
0
0
0
0
N5
0
0
0
0
0
Nl/x3
0
N2/x3
0
0
NI/x3
0
0
0
0
0
0
0
0
0
0
N4
0
0
N4
0
0
0
N6
0
0
0
0
N3/x3
0
N4/x3
0
0
N3/x3
0
0
0
N4/x3
0
0
0
0
0
0
N6
~
[N] =
N2/x3 0 0
A well-known linear polynomial can be applied to axial and torsional displacement, i.e. N5 and N6, while Hermitian cubic functions and their derivatives are needed for remaining terms, i.e. N1 to N 4. Then, the equivalence of the virtual work exploited by the inertial forces can be forced by: 6{W}'r[91//]{ff "} = IL(6{r}T[m]{~'}) dL
(39)
where the 12 × 12 mass matrix of the two node beam is introduced: f [Yff] = JL[N]T[m][N] dL
(40)
In eqns (26) and (37) the nodal displacement and force arrays are defined in the local reference frame of the beam, and generally a transformation should be used to write equilibrium equations in the basic reference frame.
Thermal equivalent loads A finite element approach to the thermal problem for a beam element 13 leads to the particular solution associated
N5
the first beam end. The application of eqn (41) is trivial, but a more general and practical procedure can be useful in a three-dimensional frame analysis. If the elastic and thermal problems are uncoupled, the presence of thermal strain simply leads to an additional term depending on cross-section strain parameters. The following thermal contribution to the complementary virtual work is identified:
l I' _ i {~0}T{~T} d~
6Li, Therm = ~
(42)
The stress distribution resulting from thermal anisotropy or non-homogeneity is self balanced. So, the term of the internal complementary work related to the cross-section force resultants is not affected by their presence and the contribution of eqn (42) is simply added to the complementary work equation. Remembering eqn (31) we obtain:
6Ci,zhorm={rS} T
l[Mlr ¢0}+5(~+1){¢, } d (43)
619
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli x A
EXAMPLES
V 10t
l
S.C.
Table 1
~" xI
N.C.
In this final section some examples are presented covering a large range of applications and proving the capability and efficacy of the two-dimensional approach to the analysis of composite beams.
b
A=200mm 2 t=-lmm E=70 GPa
I;:>
<:~-q>
<3
Figure 3
V t
A
Hollow isotropic two-cells beam
G=2300 GPa
The first example refers to an hollow, isotropic, two-cells beam. Theoretical solutions under the semi monocoque hypothesis are available. The geometry is presented in Figure 3: it can be considered unreasonable but it outlines
2b Hollow isotropic test
Stiffness properties of the cross-section ANBA
Axial stiffness, EA (N) Bending stiffness, E1 (N nlrfl 2) Torsional stiffness, G J 0 (N m m 2) Shear stress coordinates (mm) Normal stress coordinates ([ram)
Table 2
y
x
y
0.84 × 107 0.56 X 10 u 0.1877 X 10 II - 71.43 0
0.21 x l 0 II -0 0
0.84 × 107 0.56 × l 0 II 0.1877 × 1011 - 71.43 0
0.21 × l 0 II -0 0
Comparison of end displacements
Hollow beam
Shear stiffness
Theory
Included No Included No Included
3D FEM Present
End displacement ([mm) 1.669 1.614 1.669 1.614 1.669
The structure of this equation allows the identification of a displacement array which, when integrated, assumes the form:
{Or}=[Ho]{~bo}+[Ht]{~bl}
(44)
where:
[Ho]=21I- [I6]
-[•6]
{Or} is a nodal displacement array corresponding to thermal cross-section strain parameters. Thus the loads able to determine them are simply evaluated by premultiplying with the stiffness matrix of the element:
{Sr}=[KI([Ho]{~o}+[H1]{~bl})
(45)
However, this load array can be derived in a straightforward fashion by writing the complementary virtual work equation, which leads to
{S}=[K]({W}-{OT}) and then to eqn (45). These self-balanced nodal loads are 'thermal equivalent' in the sense that they determine a beam deflection characterised by {~bT} cross-section strain parameters and can be merged with the external forces in static problems.
620
Theory
x
× x x X x
10 -2 10 -2 10 -2 10 2 10 -2
End slope (rad)
End twist (rad)
3.805 X 10 -6 -3.805 x 10 -6 __ 3.805 x 10 -6
2.381 x 10 -5 2.381 x 10 -5 2.381 x 10 -s
the effects of the non-alignment of shear and normal stress centres. Two thin elements, working in plane stress state, were employed to model each panel, while stiffeners were modelled as concentrated axial elements. In Table 1 the section properties obtained with the computer code ANBA, developed by Giavotto et al.9, are compared with theoretical counterparts in terms of stiffness (axial, bending and torsional) and coordinates of the normal and shear stress centres. A clamped beam having this cross-section is loaded with a unit transverse force at its tip. In Table 2 the displacements and rotations at the free end are reported. The results obtained with the present formulation assume that only out-of-plane cross-section warpings are present. The values from the theoretical approach are equal to those obtained by using a BEAM element in the MSC/NASTRAN code once the gap between neutral and elastic axis is furnished. It can be seen that the present results automatically take into account all the elastic properties and also the shear stiffness contribution.
Thin-walled composite beam This example matched the experiment of Smith and Chopra7: it refers to thin-walled composite beams with rectangular cross-section, the geometric dimensions of
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli Table 3
Physical geometry of the box-beam
Dimensions (in) L d c Ply thickness Wall thickness (six plies)
Table 4
Lid = 56
L/d = 29
30 0.537 0.953 0.005 0.030
30 1.025 2.060 0.005 0.030
l/
3D FEM Model L/d=29 Specimen
/"
/
Properties of ply material (unidirectional graphite/epoxy)
Ell (msi) E22 (msi)
20.59 1.42 0.87 0.42
G 12.23,12 ( m s i )
u12
Y Z ,[.... X Table 5
Summary of reference cases
Specimen #
Lay-up
L/d
Top and bottom
Sides
Antisymmetric Antisymmetric Antisymmetric Symmetric Symmetric Symmetric Symmetric Symmetric Symmetric
56 56 56 29 29 29 56 56 56
(15)0 (0/30)3 (0/45)3 (15)6 (30)6 (45)6 (15)6 (30)6 (45)6
(15)6 (0/30)3 (0/45)3 (15/ -- 15)3 (30/-- 30)3 (45/--45)3 ( 1 5 / - 15)3 (30/-30)3 (45/-45)3
Symmetric layup
Figure 5
Anti-Symmetric layup
J d
Elastic Couplings: Bending/Torsion Extensio~'Shear
Elastic Couplings: Extension/Torsion Bending/Shear
Figure 4 Composite lay-ups and elastic couplings
which are summarised in Table 3. The structure is clamped at one end and loaded with transverse force, torque and axial load at the free end. Material properties are shown in Table 4; the units adopted in the reference paper are maintained in this comparison. The test allowed Smith and Chopra to investigate various methods to describe experimentally verified elastic couplings. In fact, several stacking sequences, summarised in Table 5, were analysed showing different coupling behaviour according to the scheme depicted in Figure 4. For the sake of brevity only the results relating to the coupled cases are compared. Specimens are classified into symmetric and antisymmetric configurations (specimens # 1 - 3 and #4-9, respectively). Two aspect ratios were analysed in the symmetric case, L/d = 56 (#4-6) and L/d = 29 (#7-9), while only L/d = 56 was analysed for antisymmetric specimens (#1-3). The present results are obtained by using a 2D mesh
3D FEM model for L/d = 29 specimen
of the cross-section, realised with thin laminate panels, to evaluate the stiffness matrix: isoparametric, three-noded elements were used, upper and lower walls were modelled with four elements while lateral walls were modelled with two elements. Then, a frame analysis was performed by means of a single beam element built using eqn (36). Because of some discrepancies found during the comparisons of these results, 3D meshes were prepared to confirm the results of the present approach. The commercial code I-DEAS from SDRC was used. In Figure 5 the mesh related to the L/d = 29 specimen is shown: 600 isoparametric, eight-noded laminate elements were adopted. The model relating to the slender case (L/d = 56) has the same topology and is obtained by simply scaling the node coordinates. The boundary conditions at the restrained end of the beam consist of perfect clamps on all the nodes. In this way, according to the Saint Venant principle, a stress state arises owing to the inconsistency between central and extremity solutions but, as the beam is slender enough, these effects are negligible. The results obtained by the author are referred to as 'Present' and '3D FEM' in figure legends; they are compared with data gathered from Smith and Chopra 7 and related to tests (Experimental), analytical tools (Chopra) and a finite element approach (Beam FEM) using the formulation presented by Stample and Lee 5. According to the reference paper, results are compared in terms of rotations because of their sensitivity to the precision of the method. Figures 6-8 show some comparisons in terms of the twist and bending slope at the middle of the beam due to tip torque, and twist due to axial load. The results obtained with the present formulation and the 3D FEM analysis overlap and agree very well with the cited Beam FEM. The comparison is also generally good with respect to experimental and analytical counterparts, but exceptions were found in some cases. Significant
621
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli 2.0
(a)
[] Present • 3D FEM
1.8 -- ' 7 " ~
• Experiment [] Chopra • Beam FEM
1.6 "
1.4
O
[] [] • [] •
1.4 - (b) 1.2
g
-[
1.0
Present 3D FEM Experiment Chopra Beam FEM
1.2
2 08
1.0
o
~
0.8 .~ 0.6 [0.4
0.6
"~ 0.4 0.2
0.2 Spec #1
Spec #2
0 ....d.__
Spec #3
Spec #2
Spec #1
Spec #3
Twist for specimen with L/d = 56 (asymmetrical lay-up): (a) unit tip torque and (b) unit axial load
Figure 6
1.6
[] • • • •
1.4 1.2 O
1.0
0.8
Present 3D FEM Experiment Chopra Beam FEM
g 0.6
0.8
o
[] Present • Chopra • 3D FEM • Beam FEM [] Experiment
(b)
eL O
=~ 0.6
0.4
~0
.=. 0.2
",~ 0.4 0.2 0
Spec #7
Spec #8
Spec #9
'
0
Spec #7
Spec #8
Spec #9
Twist (a) and bending slope (b) for unit tip torque in specimen with L/d = 56 (symmetrical lay-up)
Figure 7
1.0-
2.0(
[] [] • • •
(a)
1.75
~1.50 o
1.25
0)
Present 3D FEM Experiment Chopra Beam FEM
0 . 9 - (b) --, 0 . 8 -
[] Present [] Chopra [ ] 3D FEM • Beam FEM • Experiment
0.7 0.6 0.5-
1.00
O
0.4
0.75
08
0.3 "~ 0.50
¢~ 0.2
0.25
0.1 Spec #4
Figure 8
Spec #5
Spec #6
Spec #4
Spec #5
Spec #6
Twist (a) and bending slope (b) for unit tip torque in specimen with L/d = 29 (symmetrical lay-up)
discrepancies are present with respect to Experimental and Smith & Chopra's results for specimen #7 and #6, Figure 6(a) and Figure 7(b) respectively. So, for symmetric cases, the coupling parameters were evaluated by stepping the stacking sequence every 5 °. These results are depicted in Figures 9 and 10 from which it can be seen that the mismatched points
622
0
are quite far from a regular curve describing the phenomenon. It must be noted that the shape of these curves is characteristic of the change in properties due to a progressive orientation of fibre. Usually, a significant variation first occurs with lamination angle between 0 and 2 0 - 3 0 °, while little changes arise from 30 to 45 ° .
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
i (a) i 0 0 0 0
I
....
'
~x 0,6
&
I
~
;
0,6
!
underio.ueJ
........................................................................................
I!Symm. Config. I/d=56
0 0
4
:
,
Q-
~
..................... -----U-,I-+7-~7
(a)
I
I
0,8
0,8
+ :i: ~ ~ : ; i _ Symm. Conflg. L/d=29 under torque I
0,4 +
- !
!
•
,
m ~ 0,2
0,4
03
Present
¢-
• - Present
Exper
~5 c 0,2
•
Exper
Chopra
El
•
Chopra
o
Beam F E M
Beam FEM 0
/
..............................
0
5
10
15
20
25
0
-
.
30
35
40
0
45
5
10
i
2,5 (b)
I i
_
. . . . . . Symm. Conflg. I/d=29 under torque I i
i
2
20
25
30
35
40
!
i
__
~___
45
]
Symm. Config. I/d=56 under torque J :
I
~
15
Fibre orientation (deg)
Fibre orientation (deg)
I
1
-
g O O
o
o°° i
p. . . . .~. . . . ._. . . . .~. . . ~
~.
+l
'- ~
.
•
°
-~|
~- 0,5
r~ Chopra
I
+
-
~
-
" BeamF EMI 0
5
10
0
15
20 25 30 F i b r e odentation (deg)
35
40
-
, • tJ
•
I
0
;. :
• Exper 0,5
1
I 5
45
Present ~ l Exper
1 --
Chopra
/
'
o
+
+
•
i
Beam FEM I 10
15
20
25
30
35
40
45
Fibre o r i e n t a t i o n (deg)
Figure 10 Bending slope (a) and twist (b) for unit tip torque versus fibre orientation for symmetric lay-up (x/L = 0.5, L/d = 56)
Figure 9 Bending slope (a) and twist (b) for unit tip torque versus fibre orientation for symmetric lay-up (x/L = 0.5, L/d = 29)
Finally, Figure l l shows a comparison of the axial behaviour of bending slope in the case of symmetric lay-up for specimens #1 and #2. These results demonstrate that the present beam model automatically includes the shear flexibility properties of the beam and that no specialised
Table 6 Mechanical properties of isotropic materials E (kg mm 2) v (mm mm ~ °C ~)
Steel
Aluminium
21 000 0.31 12 × 10 -6
7000 0.3 24 × 10-6
Table 7 Mechanical properties of graphite fibre E u (kg mm -2) 14074
E22,33
(kg mm 2)
v12,31
923
1478
0.209
0.33
formulations are needed to account for this aspect of the problem. Thermal deflection o f a solid beam
This example refers to a beam 18 mm long, with solid rectangular cross-section, made by two layers with different mechanical properties (Figure 12) subject to uniform heating. Two situations were examined: a bimetallic realisation, made of steel and aluminium alloy, and a composite set-up using unidirectional graphite laminae in a cross-ply configuration, i.e. (0°/90°). The mechanical properties assumed for the simulation are
GI2,23,12 (kg mm 2) c~z (ram mm -I °C-l) 598
0.0
Ol2, 3
(mm mm I oC ]) 20 × 10 -6
summarised in Table 6 and Table 7. In Table 7 a thermal expansion coefficient typical for an epoxy resin is assumed along the transverse direction. No corrections were applied to take into account the actual graphite expansion coefficient to outline the effect of this kind of anisotropy. The specimen section was 6 mm width and 1 mm thick and was subjected to a +20°C uniform heating. The dimensions employed are the minimum required to avoid any interference from free-edge effects at the middle of the specimen. The results presented are related to a meshes having 12 two-dimensional, isoparametric, parabolic elements along the specimen width and eight elements along the thickness
623
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
~
(a) 16 /' [] - .....
~, 12 --
Present Smith & Chopra experiment Smith & Chopra analitycal without shear effect Smith & Chopra with shear effect []
X2
[]
zX []
A X3
Figure 13
The structural analysis with thermal strain gave the following results in terms of cross-section strain parameters for the bimetallic and composite specimen:
SS ~ s J ssSS SSSSSS
4
Bimetallic Composite
s SSS
~
(b) 1 2 -
I I 0.4 0.6 Normalised abscissa
0.2
t,
[] - 10.....
I 0.8
I 1.0
Present Smith & Chopra experiment Smith & Chopra analitycal without shear effect Smith & Chopra with shear effect [] []
/x
~ 6 ,~
[]
s SS
ssS
4
ssss
I
0.2
Solid specimen: two-dimensional mesh
I
I
0.4 0.6 Normalised abscissa
I
0.8
I
1.0
Figure 11 Effect of the stiffness for L/d = 56: (a) specimen #1; (b) specimen #2
{~br} = [000.3413 × 10-30.3316 X 10-300] T {~br} = [00 0.1230 × 10 -3 - 0 . 4 0 9 6 x 10-30 0] T
Corresponding to this thermal strain are a uniform axial force (bimetallic, 25.1876 kg; composite, 9.6201 kg) and a bending moment (bimetallic, - 1.2624 kg m; composite, -2.7578 kg m), these values being determined from eqn (42). As far as displacements of the free end are concerned, the results shown in Table 8 were obtained by a static analysis performed using the stiffness matrix of eqn (36) and the external loads given previously. Because of the symmetry of the problem, the overall deflection of the equivalent beam was evaluated by solving a frame model with a single element, 9 mm long, clamped at one end. These results are compared with the outputs of a threedimensional reference analysis in Table 8. It has been carried out by using a mesh composed of 216 isoparametric, parabolic, brick elements leading to a total of 3675 degrees of freedom. It represents a quarter of the complete beam (Figure 12), appropriate boundary conditions being applied. The 3D displacements were evaluated by computing the least-squares fit of axial and transverse displacements at the middle of the specimen end. The results well agree in the bimetallic case. Some discrepancies occur in the composite test but they are small enough, a few per cent, and are probably a consequence of specimen dimensions. In fact, a test with a shorter beam length, i.e. 3 mm only, gave worse results.
Modal frequencies of a helicopter rotor blade
Figure 12
The solid specimen and three-dimensional mesh
(Figure 13). Mesh statistics refers to 96 elements and 993 degrees of freedom. Because of the very small mesh dimensions, the whole specimen was modelled without paying any attention to structural or load symmetry.
624
The last example refers to evaluation of the frequencies and modes of a helicopter rotor blade, data for which are available by courtesy of Agusta. The blade, depicted in Figure 14, is constructed from an aluminium alloy spar wrapped in an extemal sheet also made of an aluminium alloy. A rear spar closes the cross-section, while the back of the profile is honeycombed. In this figure solid lines indicate the 14 sections where the cross-section properties were evaluated. From these data a three-dimensional beam model was built with 30 constant-property elements. They are shown in Figure 14 by dashed lines vertical while triangles specify the beam portion to which constant properties were applied. In Figure 15 the geometry of some cross-sections, identified with A and B in Figure 14, are sketched; the
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli Table 8
Comparison of beam-end displacements under thermal loads Axial End Disp.
Bimetallic Composite
~
Present
3D
A%
Present
3D
A%
Present
3D
A%
3.080 x 10 -3 1.143 x 10 .3
0.27 3.12
13.43 x 10 .3 16.59 x 10 .3
13.48 )< 10 .3 16.52 X 10 3
0.37 0.41
2.984 x 10 3 3.686 × 10 3
3.014 x 10 -3 3.766 x 10 .3
0.98 2.12
I
f
I
I
I
I
I
I
I
i
i
~ i
I
I
i
i
r
i
Figure 14
~
End slope
3.072 x 10 -3 1.107 × 10 3
I
~
Transverse End Disp.
~
.
~
I t
~
J t
~
I i
I ~
'
I t
X
I I
I ~
I i
~
i
I
F
]
I
I
I
:
i
r
~
~
~
~
:
I t
i ~
~
ModeIDiscretisation Influence
Length
Analysed
Sections
Helicopter rotor blade
Table 9 Comparison of frequencies for a helicopter rotor blade Mode no.
Section A
Section B
Figure 15 Section geometry and mesh examples
finite element meshes used in these cases are also depicted: solid elements were adopted to model the main spar and the honeycomb while thin panels were used to model the skin and the anti-abrasive steel strip. Glue and paint are included in the finite element models of the cross-sections as nonstructural masses. Elastic and mass properties, such as stiffness, centre of gravity, shear and neutral axis locations, are obviously not constant along the blade axis; nevertheless, a linear model is used. In fact, the geometry of each cross-section is defined by using a reference frame placed on the axis of the leading edge. All the properties are correctly defined in this frame and no more transformations or manipulations are needed. In Table 9 the frequencies of the free blade are compared with experimental data, and show a satisfactory agreement.
1 2 3 4 5 6 7 8
Mode type
Exper.
Present
Error (%)
1 Flap. 2 Flap. 3 Flap. 1 Chord 1 Tors. 4 Flap. 5 Flap. 2 Chord
7.42 21.90 42.55 43.36 65.70 70.81 105.78 121.97
7.49 21.93 43.17 44.45 66.35 73.34 110.20 125.34
0.97 0.15 1.46 2.51 0.99 3.58 4.18 2.76
structure and capable of taking into account the overall deflection due to thermal and elastic anisotropies, is finally revealed. Simple frame models were built to demonstrate the validity of the approach for a large number of cases including both thin and thick, isotropic and anisotropic beams, even in the presence of thermal strain. Finally, the method has been used to evaluate the natural modes and frequencies of a helicopter rotor blade, good agreement with experimental results being obtained.
ACKNOWLEDGEMENTS The author is thankful to Prof. V. Giavotto due to many invaluable suggestions in the development of this work.
CONCLUDING REMARKS Following a short review of the theory of anisotropic beams, this paper described the formulation of the space stiffness matrix of a two-node beam element, exploiting the characterisation performed with a finite element, twodimensional analysis of the cross-section. This work is the natural completion of previous papers where the approach to the analysis of composite straight, untapered and untwisted beams was described. A way to compute an array of thermal equivalent loads, suitable to be applied at the ends of a beam in a 3D frame
REFERENCES 1. 2.
3. 4. 5.
Hodges, D. H., Review of composite rotor blades modeling. A1AA J., 1990, 2(3), 561-565. Giavotto, V., Borri, H., Puccinelli, L., Caramaschi, V. and Mussi, F.,. Evaluation of section properties for hollow composite beams. Presented at 'V European Rotorcraft and Powered Lift Aircraft Forum', Paper no. 35, 1979. Bauchau, O.A., A beam theory for anisotropic material. J. Appl. Mech., 1985, 52, 416-422. Bauchau, O.A. and Hong, C.H., Finite element approach to rotor blade modeling. J. Am. Helicopter Soc., 1987, 32(1), 60-67. Stample, A.D. and Lee, S.W. Finite element model for composite
625
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
6. 7.
8. 9.
10.
626
beams with arbitrary cross sectional warping. AIAA J. 1988, 26(12), 1512-1520. Cofie, E. and Bank, L.C., Analysis of thin-walled anisotropic frames by direct stiffness method. Int. J. Solids Struct., 1995, 32(2), 235249. Smith, E.C. and Chopra, I. Formulation and evaluation of an analytical model for composite box beams. Presented at '31st AIAA/ AHS/ASME/ASCE/ASC Structure, Structural Dynamics and Material Conference', Long Beach, CA, 2 - 4 April 1990. Barbero, E.J., Lopez-Anido, R. and Davalos, J.F., On the mechanics of thin-walled laminated composite beams. J. Compos. Mater., 1993, 27(8), 806-829. Giavotto, V., Borri, M., Mantegazza, P., Ghiringhelli, G. L., Caramaschi, V., Maffioli, G. C. and Mussi, F., Anisotropic beam theory and applications. Computers and Structures, 1983, 16, 303413. Kosmatka, J.B. 'Structural dynamic modeling of advanced composite propellers by the finite element method', Ph.D. Dissertation, University of California, Los Angeles, 1986.
11.
12. 13. 14. 15. 16. 17.
Hodges, D.H., Atilgan, A. R., Cesnik, C. E. S. and Fulton, M. V., On a simplified strain energy function for geometrically nonlinear behaviour of anisotropic beams. Compos. Eng., 1992, 2(5-7), 513 -526. Borri, M., Ghiringhelli, G.L. and Merlini, T., Linear analysis of naturally curved and twisted anisotropic beams. Compos. Eng., 1992, 2(5-7), 433-456. Ghiringhelli, G.L. On the thermal problem for composite beams using a finite element semi-discretisation. Composites: Part B, 1997, 28B, 483-495. Giavotto, V. Teoria elastica della trave anisotropa e non omogenea. In 'Scritti per Lucio Lazzarino', Pacini, Pisa, 1986, pp. 131-164. Ghiringhelli, G.L and Mantegazza, P., Linear straight and untwisted anisotropic beam section properties from solid finite elements. Compos. Eng., 1994, 4(12), 1225-1239. Borri, M. and Merlini, T., A large displacement formulation for anisotropic beam analysis. Meccanica, 1986, 21, 30-37. Bicos, A.S. and Springer, G.S., Design of composite boxbeam. J. Compos. Mater., 1986, 20(1), 86-109.
N;" I
PII: S1359-8368(96)00075-3
ELSEVIER
© 1997 Elsevier Science Limited Printed in Great Britain. All rights reserved 1359-8368/97/$17.00
On the linear three-dimensional behaviour of composite beams
G.L. Ghiringhelli Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, via C. Golgi, 40, 20133 Milano, Italy (Received 11 September 1995; accepted 20 November 1996) This paper describes a method of evaluating the linear stiffness matrix of a three-dimensional beam which accounts for all the possible couplings involved in the use of composite laminations appropriate to satisfy design requirements. The starting point of the procedure is represented by a complete evaluation of the beam crosssection properties by use of a two-dimensional finite element analysis. A mass matrix of the beam is also evaluated as well as the external forces, to account for overall thermal deflections. Some examples are presented and the results are compared either with their theoretical counterparts or with numerical results obtained by means of a full three-dimensional finite element analysis. © 1997 Elsevier Science Limited. (Keywords: B. mechanical properties; B. numerical modelling; composite beams)
INTRODUCTION Notwithstanding the availability of tools of ever-increasing power and reliability, the beam continues to be a model used widely in many branches of structural analysis, such as in the design of helicopter or aeolian generator blades, or the study of aircraft dynamics. Indeed, the advantages of an approach based on this model are well known. The development of tools able to predict the complex behaviour of composite beams has been driven by the increasing usage of these materials; there are several reasons for this, such as the significant improvement in fatigue life and damage tolerance, and the emphasis set on the possibilities offered from modem techniques of design, e.g. aeroelastic tailoring. In this context particular relevance is given to correct prediction of the elastic couplings between different load conditions. They can be introduced intentionally but can also represent an unwanted consequence of the design of composite laminations. Yet another problem is associated with both thermal load conditions and the anisotropic thermal characteristics of the materials, which can be responsible for global beam deflections. For example, in extreme environmental situations, owing to non-balanced thermal deformation, an aeroelastic tailored helicopter blade could show overall deflections that interfere with the designed structural couplings; the effort of the designer can then be nullified, possibly leading to a serious deterioration in performance. Since closed-form solutions are not general enough, the structural designer/analyst has to resort to numerical approximations which, when carried out on three-dimensional
finite element models, can be too costly. This is evident during the initial design phase, when the structural lay-out is not yet finalised and many analyses have to be carried out. The importance of the subject is confirmed by the available literature. In 1990 Hodges ~ carried out an exhaustive survey of the state-of-the-art on composite rotor blade modelling, citing almost 30 works. He pointed out that two approaches are mainly adopted, the first one being an analytical development of the problem while the second is based on some kind of finite element formulation. With regard to modelling generic, non-homogeneous and/or anisotropic beams showing any kind of thermoelastic coupling, the first category does not seem to be general enough in face of the design needs for composite beams and in particular for rotor blades. A further classification of finite-element-based methods can be made in terms of the in-plane and out-of-plane deformation of the beam cross-section. The availability of both these strain modes is strictly connected to a correct analysis of a solid beam. Several works allow out-of-plane cross-section warpings only facing the problem of thin-walled anisotropic beams; e.g. Giavotto et al. 2, Banchau 3 and Bauchau and Hong4 all adopted a two-dimensional (2D) approach to characterise the cross-section, while Stample and Lee 5 preferred a threedimensional (3D) formulation with warping nodes placed over the cross-section to describe bending, torsional and warping deformation. Similarly, Cofie and Bank 6 developed a 3D two-node beam element capable of accounting for warping and anisotropic coupling in thin-walled composite beams. These are suitable for a direct analysis of the stiffness matrix. Also papers by Smith and Chopra 7, Barbero et al. 8
613
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli and Bicos and Springer 18 suffer from some kind of limitation in blade topology. Other works finally remove this limitation: e.g. Giavotto et al. 9 and Kosmatka 1° both employ a similar finite element approach, and Hodges et al. 11 who instead exploit a variational-asymptotical beam sectional analysis. The work of Kosmatka is devoted to curved blades while the method of Hodges et al. is suitable for the analysis of straight, slender beams. Giavotto et al. consider straight beams that not are necessarily slender, taking into account a full cross-section characterisation. In fact, this work presents a completely general formulation for the evaluation of the cross-section elastic and mass properties and the stress distribution in a non-homogeneous and/or anisotropic beam. Irrespective of the geometrical shape of the crosssection (i.e. solid or hollow), thick or thin walls, single or multiple cells and type of materials adopted, the de Saint Venant approach is extended to straight, untapered and untwisted composite beams. The problem is solved by using a 2D finite element model of the cross-section without any undue simplification for geometry and material properties; the elastic characterisation of the cross-section is complete and the shear stiffnesses and related couplings are also furnished by the method. In this way great computational efficiency is gained during all analysis phases. Following the same basic idea of a semi-discretisation, the formulation was rewritten to extend the approach to curved and twisted beams by Borri et al. 12. A final extension was presented by the author t3, the basic formulation being completed to take into account strains due to the thermal anisotropy of the material and/or a complex temperature distribution on the beam cross-section. The present paper is a natural completion of these activities: it resumes previous works and describes a procedure to evaluate the stiffness and mass matrices of a two-node, three-dimensional beaml4; this element, which is suitable for use within the framework of a linear analysis, shows all the possible couplings introduced by the use of appropriate composite laminations. The starting point of the procedure is the characterisation of the beam cross-section performed with the cited 2D finite element analysis 9"15, the main points of which are briefly recalled. A theoretical section is then devoted to evaluation of the thermal equivalent loads that can be applied to a beam element to obtain a given thermal strain. The results obtained with the described procedures are compared either with their theoretical counterparts or with numerical results obtained by a full 3D finite element analysis. A final example compares the experimental frequencies of a helicopter rotor blade with results from an eigenvalue analysis on a beam model.
Tx2~_~~ M2
7j"
X3
M 3~
T1
XI
Figure 1 Beam reference frame and section load convention
to Giavotto et al. 9 and Ghiringhelli and Mantegazza 15 for an exhaustive discussion of the formulation. The classical theory of elasticity defines a beam as a prismatic, slender and straight body having both geometric and elastic properties that are constant along each generatrix. A section is obtained by cutting the beam along a plane normal to its axis in the undeformed configuration. According to the de Saint Venant hypotheses, we assume that loads applied along the axis produce negligible effects with respect to those due to extremity loads. Furthermore, the effects of their distribution at the beam extremities vanish far from the end at a distance that is small with respect to the beam length. To define the deformation of each cross-section we assume small displacements and strains. A linear behaviour of material mechanical properties is also assumed. Arbitrarily distributed anisotropy or non-homogeneity can characterise the cross-section but, according to the previous hypothesis, they must be constant along the beam axis. If these hypotheses are matched the elastic problem for a beam is greatly simplified, being reduced to a monodimensional problem. So the stress distribution on each cross-section can be analysed independently from each other, once the distribution of internal forces is known. It must be noted that linearity assumptions, in both displacements-strain and stress-strain laws, do not prevent a nonlinear analysis of a body exhibiting large displacements but small strains; in fact, the linearity assumptions refer only to the local behaviour of a section of the beam t6. The beam geometry, depicted in Figure 1, is defined in an orthogonal Cartesian reference frame, {x} [XlX2X3] T, with the axis x3 parallel to the beam generatrix. Let us define the strain and stress columns in the following form: =
{/3} ~- [/311 /322 2/321 2/331 2/332 /333] T
{ o / = I°,1 °22 o21 °3, ANISOTROPIC BEAM THEORY In this section some basic elements of the formulation of the elastic problem for a non-homogeneous/anisotropic beam are presented. They are fundamental for the arguments that will be presented in the following text. The reader is referred
614
so that the stress components acting on a section, i.e. tY3i, i ---1,3, can be easily isolated by formally extracting the second half of the stress vector with an appropriate matrix operator: {p} = [s]{~}
where [S] T = [0 I3], [I3] being a third-order identity matrix.
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
x2l TI,2
M T
Ml:~
~ X3
a
1,2
~
x,
Tl,i
T2.1
Figure 2 Three-dimensional elemental frame
In a beam model the displacement of a generic point is described by means of rigid displacements and rotations of the origin of the reference frame, while the stress state is reduced to resultants of internal forces and moments. If {T} and {M} are the arrays of resultant forces and moments, respectively, the conventions adopted being depicted in Figure 2, we can define the cross-section load vector as {0} T = [{T}T,{M}T]. It can be evaluated from: {0} =
~A[Z]T{p} dA
(1)
i.e. by integrating the elementary forces and moments due to stresses acting on the cross-section. The matrix [Z] is defined in such a way that the cross-product matrix related to the point with coordinate {x} is included:
[z] = [I (XA)T] The displacements, {X}, and rotations, {~o}, of a reference point belonging to the corresponding cross-section can be organised into a single array, {r}: {r} = [{x}TI{so}T] T When significant warping is present, i.e. in composite beams, the displacement of a generic point of the crosssection is profitably described as the sum of a term depending on the rigid movement of the cross-section and a warping term straining it, {g}: {s} = {g} + [Z]{r}
(2)
The sixth-order matrix [O] has only two non-null terms according to eqn (4):
[O]=
,
0
0
0
{g(xl,x2,x3) } = [Ni(xl,x2)] {u(xil,xi2,x3) }
(4)
(5)
eqn (5) states that nodal warpings are discrete parameters on the cross-section, still depending on the position along the beam axis; i.e. a semi-discretisation of warping is used. The following corresponding approximations for {s } and {e} are obtained: {s} = [Nl{u} + [Z]{r} {el = [B]{u} +
[S][NI{uI/3[S][ZI{d/}
(6) (7)
where [B] = [©][N] and
(3)
It takes into account the strain due to differential rigid displacement of adjacent cross-sections. The subscript /3 indicates differentiation with respect to the x3 direction; this notation will be used in the following developments.
1
The first two elements of {~b}, i.e. ~b, = Xl/3 - ~o2 and ~2 = X2/3 + ~°1, represent the transverse differential displacement, i.e. shear strain of the cross-section; the third term, i.e. ~b3 = X3z3, is the axis elongation; while the three last elements are associated with bending, ffl = ¢t/3 and ~b2 = ~o2/3, and torsion, ~b3 = ~o3/3, of the reference line, respectively. After discretising the beam section with appropriate twodimensional finite elements (Figure 2), we approximate the warping function {g} by interpolating the nodal values {u} with usual shape functions [N] with:
By writing strains as differentiation of displacements allows us to introduce the cross-section strain array, {~b}: {f} = {r}/3 q- [O]{r}
[t]=
0
"/1
0
0-
0
/2
0
12
/1
0
0
0
/1
0
0
/2
0
0
0
[D] =
615
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli The virtual work principle per unit length of the beam is written as
fA(~{s}T{p})/3dA=
IA~{8}T{0r} dA
(8)
By substituting eqns (6) and (7) into eqn (8) and carrying out the explicit differentiation of the term appearing in the lefthand side, we can write, recalling eqn (1):
{u/3/TI il{u/3/ LTCcT u { U/3}T/i 3 + {rrlT({o}13 -- [oIT{0})
(9)
where the following definitions are introduced: [M] =
JA[N]T[s]T[D][S][N]dA
{0} = ([I] +X3[O]T}{00}
(12)
where {00} is the resultant array at the origin. We then try to find a particular solution for {u}, {P} and {~b} having an equivalent linear form, i.e.
~A[B]T[D][B]dA
[El =
[C] =
(10)
depends on the way the beam ends are constrained and is generally negligible after a length of the order of the largest linear dimension of the section itself. Beyond this length at both the beam extremities, i.e. in its central portion, only the particular integral, which will be seen to depend only upon {0}, is significant. This is in accordance with Saint Venant's principle and consequently the general and particular integral were called by Giavotto et al.9 the extremity and central solution, respectively. The central solution is thus related to the determination of the usual global elastic properties and stresses of a beam section. The reader can refer to 9,a4for a comprehensive explanation of the complete solution of equations (11). The central solutions must satisfy eqn (1 ld) for generic extremity loads, so we directly look for linear internal resultants {0} of the type:
{u}= ([.o] +x3[ulJ){0o}
~A[B]T[D][S][N]dA
(13)
{P}= ([P0] +x3[P1]){0o} [R] =
~A[B]T[D][S][Z]dA
[L] =
~A[N]T[s]T[D][S][Z]dA
[A] =
fA[Z]T[s]T[D][S][Z] dA
{P} =
{P}/3 =
=
where [Uo], [ud, [Po], [P1], [~bo] and [~bl] are solution matrices; each of their six columns represents the solution for the corresponding unitary term of the related constant and linear part of {0 } respectively. After substituting the above equations into equations (11) and equating the terms having the same power of x3, we end up with the following system of linear equations with unknowns [u0], [ui], [Po], [P,], [¢'o] and [~bl]:
fA[N]T{p} dA fA[N]T{p}/3 dA
These matrices are evaluated successfully by using isoparametric elements and then assembled by the usual procedure in a finite element code. For arbitrarily varying {u}, {U/3}, {~} and {r}, eqn (9) leads to the following linear system of algebraic equations in {~b} and ordinary, first-order differential equations in {u }, {P} and {0}:
(14)
[pl] = tc]T[u,] + u0
-- C
(15) (16)
{P} = [M]{u}/3 + [c]T{u} -{- [L]{~}
(1 la)
[Po] = [M] [ul] + [c]T[uo] + [L] [~bo]
[C]{u}/3 -t- [E]{u} + [R] {~b}
(1 lb)
[L]T{u}/3-q-[R]T{u} -4- [A]{~b}
(1 lc)
The solution of the above sets of equations is naturally carried out by solving eqns (14), (15), (16) and (17). The right-hand side of eqn (14) reduces to only two columns because of the structure of [®] [see eqn (4)]; it must be noted that eqns (14) and (16) have the same coefficient matrix, thus a single factorisation is required for their solution. This is carded out after removing, in any convenient way, the six singularities due to the unknown warping {u } which, because of the properties of the finite element shape functions, already contains the rigid motions explicitly accounted for by {r}. Because the warping {g} contains
{P}/3 = {0} =
{0}/3 = [®]T{0}
(lld)
It can be seen that the equilibrium equations of a slice of infinitesimal length, i.e. Eqn (1 ld), are obtained. The solution of the above system is composed of a general and a particular integral. For a sufficiently long and slender beam, the general integral is related to the diffusion of stresses from the terminal sections of the beam itself. It
616
+x3
(17)
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli arbitrary rigid contributions owing to the applied constraints, a procedure to remove them is needed, as explained previously 14,16.
u/3/aNlS/j
additional term to eqn (9), the expression of the virtual work principle:
/
CROSS-SECTION STIFFNESS/COMPLIANCE MATRIX To recover the elastic properties of the beam we note that the central solution can be written by using eqn (13), so that the external virtual work of the beam slice, i.e. the first term on the right-hand side of eqn (9), can be written in a very concise form:
~{00}T([~0] T q-[u0]T [PI] "1=[t/lIT [P0]) {00} =6{00}T[F]{00}
(18)
The analysis of this equation allows us to define the compliance matrix of the section in the following form:
[F]=[J/o]TW[uo]T[P,]+[ul]T[po]
(19)
This can be shown to be symmetric by taking into account the first of eqn (16), eqn (17) and the expression for the identity matrix [/] arising from the second of eqn (16). Furthermore, it is independent of any rigid motion of the section, as demonstrated by Ghiringhelli and Mantegazza 15. In particular, it does not depend on the way in which the partition of eqn (2) is performed; in fact, constraints applied to nodal warpings to remove singularities from eqns (14) and (16) affect only displacements while stresses and the compliance matrix are not influenced. This agrees with the definition of that matrix as an intrinsic characterisation of the overall elastic properties of the beam. Conventional section parameters, used in standard engineering beam theory, can be obtained from the compliance matrix of the section. If the [F] matrix is diagonal or if couplings are restricted to axial/bending and/ or shear/torsion load conditions, engineering parameters, such as stiffness, shear centre and normal stress centre locations, can be evaluated so that the stiffness matrix of a beam is built with known techniques. But when more complex couplings are present, e.g. in the case of composite beams, the matrix could be fully populated. Some terms can be physically interpreted: by means of coefficients (4,6) and (5,6) the orientation of the torque moment that does not couple with bending moments can be defined as:
=
~
fA[ZIT[s]T[D]{¢T}dA
The temperature behaviour along the beam axis can be modelled with a polynomial expansion across the analysed cross-section, but higher order contributions are not of interest because, according to the Saint Venant principle, local effects due to distributed loads should be negligible. So a linear development is assumed, i.e. using nodal temperatures {t} and nodal temperature gradient {t'} vectors, and the thermal strain is expressed by: {eT} = {~}[N]T({t} +x3{t' })
(21)
where {o~} is the vector form of the thermal expansion tensor, i.e.
331T
=
In eqn (21) the distribution on the cross-section of the temperature and its axial gradient are obtained by the solution of the thermal problem 13. By using this equation the following matrices of nodal thermal influence coefficients can be defined:
IHT] fa [N]T [S]T[D]{~ }[N]T dA =
~A[B]T[DI{ot}[N]TdA
IVy] =
=
(22)
fA[Z]T[s]T[D]
{or}[U] T dA
and eqn (20) becomes:
~£Therm=
aU ~¢
[gl({t} q-x3{/'}) [W] ({1} +x3{t'})
(23)
Coherently with eqn (21), we still look for a linear {u}, {~b} and {P}:
° t x l = t a n - l ( -F46~F44j' ° t x 2 = t a n - l ( -F56~F55j A general procedure is needed to obtain a model showing the same behaviour of the original beam; this is the goal of a subsequent section.
THERMAL STRAIN The presence of thermal strain can determine the global deflection of a beam and this problem has great importance because when some constraints are applied, a severe tension state can be introduced. Thermal strain results in an
{u} = {u 0} + x 3 {u, }
(24a)
{~b} = {~b0} + x3{ ~bl }
(24b)
{P} = {Po} +x3{P, }
(24c)
In the absence of external loads, applying the same procedure as used previously leading to eqns (14)-(17), we obtain the following form of these equations:
Ul
617
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
{P1 } = [c]W{ul } + [L] {¢~ } -- [H] {t'}
U0
-- C
(25b)
the forces acting at the two ends of the beam, i.e. {0~} and {02}, are related by the following equation:
u1
{0, } = [H]{02 }
(28)
+ [:]{t}
(25C)
In addition, remembering that {0} is relative to cross-sections having the outgoing normal in the direction of x3, the nodal forces of the three-dimensional element can be written as:
{Po}=[M]{u,} +[c]T{u0} +[L]{¢o}-[H]{t}
(25d)
{S}~'{ -{01} }{02}
(29)
{s} = [L] {o2}
(30)
where the unknowns {u0}, {ul}, {P0}, {P1}, {¢0} and {¢1} are related to a specific temperature distribution given by {t} and {t' }. While eqns (14)-(17) are solved only once for unit external loads, in combining these solutions to obtain actual loads, the analysis of several thermal load conditions needs forwards/backwards substitution processes to be repeated for each set of temperature distributions. A restart procedure can be used successfully in this case. It can be appreciated that the assumption of linear behaviour along the axis enables us to keep the way of solving eqns (14)-(17) with the only modification being evaluation of the right-hand term. Nonetheless, higher terms could be accounted in a very simple way leading to further computational steps of the same kind as in eqns (25c) and (25d).
THE EQUIVALENT THREE-DIMENSIONAL BEAM
or, by eqn (28):
wherein the following matrix has been introduced: [LI=
From previous considerations, the internal forces in a genetic point, i.e. ~', can be interpolated from the end values, and from eqn (29) we obtain: {0} = [M]{S}
(26)
where load and displacement arrays, {S} and {W} respectively, are defined with the usual assumptions These terms are defined in the elemental frame so that the equilibrium equations expressed by the stiffness matrix must be rewritten in the basic frame by means of an appropriate transformation. To develop the formulation let us consider a generic beam, of length l, whose cross-section is characterised by the compliance matrix [F]. A natural local frame is conveniently assumed, the origin of which is placed at the middle of the beam. A normalised abscissa ~"is used for convenience, it ranges from - 1 to 1 when spanning from the first to the second node. According to beam theory, the equilibrium equations of a slice of infinitesimal length are
(Figure2).
{0}/3=[o]T{0} {0}/33=0
[M]=I~(~-l)[16] 1(~+1)[16]J Remembering equations (25), the virtual variation of the internal complementary work of the beam can be written as:
Thus the internal forces vary linearly along the beam and their end values are correlated; in fact, by integrating eqn (27) and defining a matrix [H] in the following form: [H]=
618
l[t]
[I1
6Li= ~lfl _ 1{60}y{F]{O}d~
(32)
while the term related to external forces is given by: 6Le = {rs}T{ W}
(33)
By introducing eqn (31) into eqn (32) and integrating we obtain the following form of the internal complementary work: ~Li = {c5S}T[p] {S }
(34)
where: I I2[F] [F]=6[-[F]
-[F]] 2[F] J
The complementary work equation is written by using eqn (30) to obtain an equation showing independent virtual variations only, i.e. Because of the arbitrariness of these virtual variations we obtain:
{602}.
{02} = ([L]T[pl[L]) -I[L]T{w} (27)
(31)
with
The aim of this section is to formalise a procedure to evaluate the 12th-order stiffness matrix of a simple, linear, 3D beam element with two nodes, in the form: {S} = [K]{ W}
[I] J
(35)
Finally, premultiplying by [L] and remembering eqn (30), we obtain eqn (26) where the stiffness matrix of the beam is defined as: [K] = [L] ([L]r [p] [L]) - 1[L]T
(36)
It must be stressed that this matrix accounts for all the global effects of the cross-section strain; so there is no need to
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli append other generalised degrees of freedom to nodal displacements accounting for cross-section warpings. During the evaluation of elastic matrices, i.e. Eqn (10), the mass matrix of the cross-section [m] can also be evaluated by a simple summation of inertia products, until the second order, integrated over each element. As far as a beam slice is concerned, once the terms due to the second time derivatives of warpings are neglected, the inertial forces per unit length related to the motion of the cross-section can be approximated by the inertial forces associated with rigid motion of the cross-section only: {FI} = [m]{}:}
with prescribed boundary conditions on the lateral surface. The effects due to the actual boundary conditions on the extremity cross-sections are neglected. This approach furnishes a temperature distribution throughout the beam volume: the behaviour is linear along the beam generatrix while the distribution on the cross-section depends on the thermal properties. This temperature distribution defines the initial thermal strain which leads to the thermal load defined by eqn (23) and the elastic problem is solved by equations (25). The displacement solution is still obtained in linear terms of strain cross-section parameters and nodal warpings. In particular, in order to describe the global deflection of the beam, we are interested in the first element expressed in the form of eqn (24). Once the section strains given by the solution of equations (25) are known, eqn (3) can be easily integrated to obtain the displacement components of the reference line at any point; i.e. in case of constant cross-section strain parameters it can be verified that:
(37)
In order to allow dynamic analyses these distributed loads must be concentrated at the beam nodes: a consistent inertial load array can be evaluated by introducing the interpolation of the rigid displacement of the reference frame at any position along the beam axis. Appropriate shape functions [N] can be applied to nodal displacements and rotations:
{ r} = [U]{ W}
(38) {r}=
According to a possible organisation of nodal unknowns and forces arrays, the [IV] operator can assume the form:
['o
{rl} + x 3 I
['o
{¢1}
(41)
where the subscript 1 indicates the boundary conditions at
"N 1
0
0
0
N2
0
N3
0
N1
0
N2
0
0
0
0
0
N5
0
0
0
0
0
Nl/x3
0
N2/x3
0
0
NI/x3
0
0
0
0
0
0
0
0
0
0
N4
0
0
N4
0
0
0
N6
0
0
0
0
N3/x3
0
N4/x3
0
0
N3/x3
0
0
0
N4/x3
0
0
0
0
0
0
N6
~
[N] =
N2/x3 0 0
A well-known linear polynomial can be applied to axial and torsional displacement, i.e. N5 and N6, while Hermitian cubic functions and their derivatives are needed for remaining terms, i.e. N1 to N 4. Then, the equivalence of the virtual work exploited by the inertial forces can be forced by: 6{W}'r[91//]{ff "} = IL(6{r}T[m]{~'}) dL
(39)
where the 12 × 12 mass matrix of the two node beam is introduced: f [Yff] = JL[N]T[m][N] dL
(40)
In eqns (26) and (37) the nodal displacement and force arrays are defined in the local reference frame of the beam, and generally a transformation should be used to write equilibrium equations in the basic reference frame.
Thermal equivalent loads A finite element approach to the thermal problem for a beam element 13 leads to the particular solution associated
N5
the first beam end. The application of eqn (41) is trivial, but a more general and practical procedure can be useful in a three-dimensional frame analysis. If the elastic and thermal problems are uncoupled, the presence of thermal strain simply leads to an additional term depending on cross-section strain parameters. The following thermal contribution to the complementary virtual work is identified:
l I' _ i {~0}T{~T} d~
6Li, Therm = ~
(42)
The stress distribution resulting from thermal anisotropy or non-homogeneity is self balanced. So, the term of the internal complementary work related to the cross-section force resultants is not affected by their presence and the contribution of eqn (42) is simply added to the complementary work equation. Remembering eqn (31) we obtain:
6Ci,zhorm={rS} T
l[Mlr ¢0}+5(~+1){¢, } d (43)
619
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli x A
EXAMPLES
V 10t
l
S.C.
Table 1
~" xI
N.C.
In this final section some examples are presented covering a large range of applications and proving the capability and efficacy of the two-dimensional approach to the analysis of composite beams.
b
A=200mm 2 t=-lmm E=70 GPa
I;:>
<:~-q>
<3
Figure 3
V t
A
Hollow isotropic two-cells beam
G=2300 GPa
The first example refers to an hollow, isotropic, two-cells beam. Theoretical solutions under the semi monocoque hypothesis are available. The geometry is presented in Figure 3: it can be considered unreasonable but it outlines
2b Hollow isotropic test
Stiffness properties of the cross-section ANBA
Axial stiffness, EA (N) Bending stiffness, E1 (N nlrfl 2) Torsional stiffness, G J 0 (N m m 2) Shear stress coordinates (mm) Normal stress coordinates ([ram)
Table 2
y
x
y
0.84 × 107 0.56 X 10 u 0.1877 X 10 II - 71.43 0
0.21 x l 0 II -0 0
0.84 × 107 0.56 × l 0 II 0.1877 × 1011 - 71.43 0
0.21 × l 0 II -0 0
Comparison of end displacements
Hollow beam
Shear stiffness
Theory
Included No Included No Included
3D FEM Present
End displacement ([mm) 1.669 1.614 1.669 1.614 1.669
The structure of this equation allows the identification of a displacement array which, when integrated, assumes the form:
{Or}=[Ho]{~bo}+[Ht]{~bl}
(44)
where:
[Ho]=21I- [I6]
-[•6]
{Or} is a nodal displacement array corresponding to thermal cross-section strain parameters. Thus the loads able to determine them are simply evaluated by premultiplying with the stiffness matrix of the element:
{Sr}=[KI([Ho]{~o}+[H1]{~bl})
(45)
However, this load array can be derived in a straightforward fashion by writing the complementary virtual work equation, which leads to
{S}=[K]({W}-{OT}) and then to eqn (45). These self-balanced nodal loads are 'thermal equivalent' in the sense that they determine a beam deflection characterised by {~bT} cross-section strain parameters and can be merged with the external forces in static problems.
620
Theory
x
× x x X x
10 -2 10 -2 10 -2 10 2 10 -2
End slope (rad)
End twist (rad)
3.805 X 10 -6 -3.805 x 10 -6 __ 3.805 x 10 -6
2.381 x 10 -5 2.381 x 10 -5 2.381 x 10 -s
the effects of the non-alignment of shear and normal stress centres. Two thin elements, working in plane stress state, were employed to model each panel, while stiffeners were modelled as concentrated axial elements. In Table 1 the section properties obtained with the computer code ANBA, developed by Giavotto et al.9, are compared with theoretical counterparts in terms of stiffness (axial, bending and torsional) and coordinates of the normal and shear stress centres. A clamped beam having this cross-section is loaded with a unit transverse force at its tip. In Table 2 the displacements and rotations at the free end are reported. The results obtained with the present formulation assume that only out-of-plane cross-section warpings are present. The values from the theoretical approach are equal to those obtained by using a BEAM element in the MSC/NASTRAN code once the gap between neutral and elastic axis is furnished. It can be seen that the present results automatically take into account all the elastic properties and also the shear stiffness contribution.
Thin-walled composite beam This example matched the experiment of Smith and Chopra7: it refers to thin-walled composite beams with rectangular cross-section, the geometric dimensions of
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli Table 3
Physical geometry of the box-beam
Dimensions (in) L d c Ply thickness Wall thickness (six plies)
Table 4
Lid = 56
L/d = 29
30 0.537 0.953 0.005 0.030
30 1.025 2.060 0.005 0.030
l/
3D FEM Model L/d=29 Specimen
/"
/
Properties of ply material (unidirectional graphite/epoxy)
Ell (msi) E22 (msi)
20.59 1.42 0.87 0.42
G 12.23,12 ( m s i )
u12
Y Z ,[.... X Table 5
Summary of reference cases
Specimen #
Lay-up
L/d
Top and bottom
Sides
Antisymmetric Antisymmetric Antisymmetric Symmetric Symmetric Symmetric Symmetric Symmetric Symmetric
56 56 56 29 29 29 56 56 56
(15)0 (0/30)3 (0/45)3 (15)6 (30)6 (45)6 (15)6 (30)6 (45)6
(15)6 (0/30)3 (0/45)3 (15/ -- 15)3 (30/-- 30)3 (45/--45)3 ( 1 5 / - 15)3 (30/-30)3 (45/-45)3
Symmetric layup
Figure 5
Anti-Symmetric layup
J d
Elastic Couplings: Bending/Torsion Extensio~'Shear
Elastic Couplings: Extension/Torsion Bending/Shear
Figure 4 Composite lay-ups and elastic couplings
which are summarised in Table 3. The structure is clamped at one end and loaded with transverse force, torque and axial load at the free end. Material properties are shown in Table 4; the units adopted in the reference paper are maintained in this comparison. The test allowed Smith and Chopra to investigate various methods to describe experimentally verified elastic couplings. In fact, several stacking sequences, summarised in Table 5, were analysed showing different coupling behaviour according to the scheme depicted in Figure 4. For the sake of brevity only the results relating to the coupled cases are compared. Specimens are classified into symmetric and antisymmetric configurations (specimens # 1 - 3 and #4-9, respectively). Two aspect ratios were analysed in the symmetric case, L/d = 56 (#4-6) and L/d = 29 (#7-9), while only L/d = 56 was analysed for antisymmetric specimens (#1-3). The present results are obtained by using a 2D mesh
3D FEM model for L/d = 29 specimen
of the cross-section, realised with thin laminate panels, to evaluate the stiffness matrix: isoparametric, three-noded elements were used, upper and lower walls were modelled with four elements while lateral walls were modelled with two elements. Then, a frame analysis was performed by means of a single beam element built using eqn (36). Because of some discrepancies found during the comparisons of these results, 3D meshes were prepared to confirm the results of the present approach. The commercial code I-DEAS from SDRC was used. In Figure 5 the mesh related to the L/d = 29 specimen is shown: 600 isoparametric, eight-noded laminate elements were adopted. The model relating to the slender case (L/d = 56) has the same topology and is obtained by simply scaling the node coordinates. The boundary conditions at the restrained end of the beam consist of perfect clamps on all the nodes. In this way, according to the Saint Venant principle, a stress state arises owing to the inconsistency between central and extremity solutions but, as the beam is slender enough, these effects are negligible. The results obtained by the author are referred to as 'Present' and '3D FEM' in figure legends; they are compared with data gathered from Smith and Chopra 7 and related to tests (Experimental), analytical tools (Chopra) and a finite element approach (Beam FEM) using the formulation presented by Stample and Lee 5. According to the reference paper, results are compared in terms of rotations because of their sensitivity to the precision of the method. Figures 6-8 show some comparisons in terms of the twist and bending slope at the middle of the beam due to tip torque, and twist due to axial load. The results obtained with the present formulation and the 3D FEM analysis overlap and agree very well with the cited Beam FEM. The comparison is also generally good with respect to experimental and analytical counterparts, but exceptions were found in some cases. Significant
621
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli 2.0
(a)
[] Present • 3D FEM
1.8 -- ' 7 " ~
• Experiment [] Chopra • Beam FEM
1.6 "
1.4
O
[] [] • [] •
1.4 - (b) 1.2
g
-[
1.0
Present 3D FEM Experiment Chopra Beam FEM
1.2
2 08
1.0
o
~
0.8 .~ 0.6 [0.4
0.6
"~ 0.4 0.2
0.2 Spec #1
Spec #2
0 ....d.__
Spec #3
Spec #2
Spec #1
Spec #3
Twist for specimen with L/d = 56 (asymmetrical lay-up): (a) unit tip torque and (b) unit axial load
Figure 6
1.6
[] • • • •
1.4 1.2 O
1.0
0.8
Present 3D FEM Experiment Chopra Beam FEM
g 0.6
0.8
o
[] Present • Chopra • 3D FEM • Beam FEM [] Experiment
(b)
eL O
=~ 0.6
0.4
~0
.=. 0.2
",~ 0.4 0.2 0
Spec #7
Spec #8
Spec #9
'
0
Spec #7
Spec #8
Spec #9
Twist (a) and bending slope (b) for unit tip torque in specimen with L/d = 56 (symmetrical lay-up)
Figure 7
1.0-
2.0(
[] [] • • •
(a)
1.75
~1.50 o
1.25
0)
Present 3D FEM Experiment Chopra Beam FEM
0 . 9 - (b) --, 0 . 8 -
[] Present [] Chopra [ ] 3D FEM • Beam FEM • Experiment
0.7 0.6 0.5-
1.00
O
0.4
0.75
08
0.3 "~ 0.50
¢~ 0.2
0.25
0.1 Spec #4
Figure 8
Spec #5
Spec #6
Spec #4
Spec #5
Spec #6
Twist (a) and bending slope (b) for unit tip torque in specimen with L/d = 29 (symmetrical lay-up)
discrepancies are present with respect to Experimental and Smith & Chopra's results for specimen #7 and #6, Figure 6(a) and Figure 7(b) respectively. So, for symmetric cases, the coupling parameters were evaluated by stepping the stacking sequence every 5 °. These results are depicted in Figures 9 and 10 from which it can be seen that the mismatched points
622
0
are quite far from a regular curve describing the phenomenon. It must be noted that the shape of these curves is characteristic of the change in properties due to a progressive orientation of fibre. Usually, a significant variation first occurs with lamination angle between 0 and 2 0 - 3 0 °, while little changes arise from 30 to 45 ° .
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
i (a) i 0 0 0 0
I
....
'
~x 0,6
&
I
~
;
0,6
!
underio.ueJ
........................................................................................
I!Symm. Config. I/d=56
0 0
4
:
,
Q-
~
..................... -----U-,I-+7-~7
(a)
I
I
0,8
0,8
+ :i: ~ ~ : ; i _ Symm. Conflg. L/d=29 under torque I
0,4 +
- !
!
•
,
m ~ 0,2
0,4
03
Present
¢-
• - Present
Exper
~5 c 0,2
•
Exper
Chopra
El
•
Chopra
o
Beam F E M
Beam FEM 0
/
..............................
0
5
10
15
20
25
0
-
.
30
35
40
0
45
5
10
i
2,5 (b)
I i
_
. . . . . . Symm. Conflg. I/d=29 under torque I i
i
2
20
25
30
35
40
!
i
__
~___
45
]
Symm. Config. I/d=56 under torque J :
I
~
15
Fibre orientation (deg)
Fibre orientation (deg)
I
1
-
g O O
o
o°° i
p. . . . .~. . . . ._. . . . .~. . . ~
~.
+l
'- ~
.
•
°
-~|
~- 0,5
r~ Chopra
I
+
-
~
-
" BeamF EMI 0
5
10
0
15
20 25 30 F i b r e odentation (deg)
35
40
-
, • tJ
•
I
0
;. :
• Exper 0,5
1
I 5
45
Present ~ l Exper
1 --
Chopra
/
'
o
+
+
•
i
Beam FEM I 10
15
20
25
30
35
40
45
Fibre o r i e n t a t i o n (deg)
Figure 10 Bending slope (a) and twist (b) for unit tip torque versus fibre orientation for symmetric lay-up (x/L = 0.5, L/d = 56)
Figure 9 Bending slope (a) and twist (b) for unit tip torque versus fibre orientation for symmetric lay-up (x/L = 0.5, L/d = 29)
Finally, Figure l l shows a comparison of the axial behaviour of bending slope in the case of symmetric lay-up for specimens #1 and #2. These results demonstrate that the present beam model automatically includes the shear flexibility properties of the beam and that no specialised
Table 6 Mechanical properties of isotropic materials E (kg mm 2) v (mm mm ~ °C ~)
Steel
Aluminium
21 000 0.31 12 × 10 -6
7000 0.3 24 × 10-6
Table 7 Mechanical properties of graphite fibre E u (kg mm -2) 14074
E22,33
(kg mm 2)
v12,31
923
1478
0.209
0.33
formulations are needed to account for this aspect of the problem. Thermal deflection o f a solid beam
This example refers to a beam 18 mm long, with solid rectangular cross-section, made by two layers with different mechanical properties (Figure 12) subject to uniform heating. Two situations were examined: a bimetallic realisation, made of steel and aluminium alloy, and a composite set-up using unidirectional graphite laminae in a cross-ply configuration, i.e. (0°/90°). The mechanical properties assumed for the simulation are
GI2,23,12 (kg mm 2) c~z (ram mm -I °C-l) 598
0.0
Ol2, 3
(mm mm I oC ]) 20 × 10 -6
summarised in Table 6 and Table 7. In Table 7 a thermal expansion coefficient typical for an epoxy resin is assumed along the transverse direction. No corrections were applied to take into account the actual graphite expansion coefficient to outline the effect of this kind of anisotropy. The specimen section was 6 mm width and 1 mm thick and was subjected to a +20°C uniform heating. The dimensions employed are the minimum required to avoid any interference from free-edge effects at the middle of the specimen. The results presented are related to a meshes having 12 two-dimensional, isoparametric, parabolic elements along the specimen width and eight elements along the thickness
623
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
~
(a) 16 /' [] - .....
~, 12 --
Present Smith & Chopra experiment Smith & Chopra analitycal without shear effect Smith & Chopra with shear effect []
X2
[]
zX []
A X3
Figure 13
The structural analysis with thermal strain gave the following results in terms of cross-section strain parameters for the bimetallic and composite specimen:
SS ~ s J ssSS SSSSSS
4
Bimetallic Composite
s SSS
~
(b) 1 2 -
I I 0.4 0.6 Normalised abscissa
0.2
t,
[] - 10.....
I 0.8
I 1.0
Present Smith & Chopra experiment Smith & Chopra analitycal without shear effect Smith & Chopra with shear effect [] []
/x
~ 6 ,~
[]
s SS
ssS
4
ssss
I
0.2
Solid specimen: two-dimensional mesh
I
I
0.4 0.6 Normalised abscissa
I
0.8
I
1.0
Figure 11 Effect of the stiffness for L/d = 56: (a) specimen #1; (b) specimen #2
{~br} = [000.3413 × 10-30.3316 X 10-300] T {~br} = [00 0.1230 × 10 -3 - 0 . 4 0 9 6 x 10-30 0] T
Corresponding to this thermal strain are a uniform axial force (bimetallic, 25.1876 kg; composite, 9.6201 kg) and a bending moment (bimetallic, - 1.2624 kg m; composite, -2.7578 kg m), these values being determined from eqn (42). As far as displacements of the free end are concerned, the results shown in Table 8 were obtained by a static analysis performed using the stiffness matrix of eqn (36) and the external loads given previously. Because of the symmetry of the problem, the overall deflection of the equivalent beam was evaluated by solving a frame model with a single element, 9 mm long, clamped at one end. These results are compared with the outputs of a threedimensional reference analysis in Table 8. It has been carried out by using a mesh composed of 216 isoparametric, parabolic, brick elements leading to a total of 3675 degrees of freedom. It represents a quarter of the complete beam (Figure 12), appropriate boundary conditions being applied. The 3D displacements were evaluated by computing the least-squares fit of axial and transverse displacements at the middle of the specimen end. The results well agree in the bimetallic case. Some discrepancies occur in the composite test but they are small enough, a few per cent, and are probably a consequence of specimen dimensions. In fact, a test with a shorter beam length, i.e. 3 mm only, gave worse results.
Modal frequencies of a helicopter rotor blade
Figure 12
The solid specimen and three-dimensional mesh
(Figure 13). Mesh statistics refers to 96 elements and 993 degrees of freedom. Because of the very small mesh dimensions, the whole specimen was modelled without paying any attention to structural or load symmetry.
624
The last example refers to evaluation of the frequencies and modes of a helicopter rotor blade, data for which are available by courtesy of Agusta. The blade, depicted in Figure 14, is constructed from an aluminium alloy spar wrapped in an extemal sheet also made of an aluminium alloy. A rear spar closes the cross-section, while the back of the profile is honeycombed. In this figure solid lines indicate the 14 sections where the cross-section properties were evaluated. From these data a three-dimensional beam model was built with 30 constant-property elements. They are shown in Figure 14 by dashed lines vertical while triangles specify the beam portion to which constant properties were applied. In Figure 15 the geometry of some cross-sections, identified with A and B in Figure 14, are sketched; the
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli Table 8
Comparison of beam-end displacements under thermal loads Axial End Disp.
Bimetallic Composite
~
Present
3D
A%
Present
3D
A%
Present
3D
A%
3.080 x 10 -3 1.143 x 10 .3
0.27 3.12
13.43 x 10 .3 16.59 x 10 .3
13.48 )< 10 .3 16.52 X 10 3
0.37 0.41
2.984 x 10 3 3.686 × 10 3
3.014 x 10 -3 3.766 x 10 .3
0.98 2.12
I
f
I
I
I
I
I
I
I
i
i
~ i
I
I
i
i
r
i
Figure 14
~
End slope
3.072 x 10 -3 1.107 × 10 3
I
~
Transverse End Disp.
~
.
~
I t
~
J t
~
I i
I ~
'
I t
X
I I
I ~
I i
~
i
I
F
]
I
I
I
:
i
r
~
~
~
~
:
I t
i ~
~
ModeIDiscretisation Influence
Length
Analysed
Sections
Helicopter rotor blade
Table 9 Comparison of frequencies for a helicopter rotor blade Mode no.
Section A
Section B
Figure 15 Section geometry and mesh examples
finite element meshes used in these cases are also depicted: solid elements were adopted to model the main spar and the honeycomb while thin panels were used to model the skin and the anti-abrasive steel strip. Glue and paint are included in the finite element models of the cross-sections as nonstructural masses. Elastic and mass properties, such as stiffness, centre of gravity, shear and neutral axis locations, are obviously not constant along the blade axis; nevertheless, a linear model is used. In fact, the geometry of each cross-section is defined by using a reference frame placed on the axis of the leading edge. All the properties are correctly defined in this frame and no more transformations or manipulations are needed. In Table 9 the frequencies of the free blade are compared with experimental data, and show a satisfactory agreement.
1 2 3 4 5 6 7 8
Mode type
Exper.
Present
Error (%)
1 Flap. 2 Flap. 3 Flap. 1 Chord 1 Tors. 4 Flap. 5 Flap. 2 Chord
7.42 21.90 42.55 43.36 65.70 70.81 105.78 121.97
7.49 21.93 43.17 44.45 66.35 73.34 110.20 125.34
0.97 0.15 1.46 2.51 0.99 3.58 4.18 2.76
structure and capable of taking into account the overall deflection due to thermal and elastic anisotropies, is finally revealed. Simple frame models were built to demonstrate the validity of the approach for a large number of cases including both thin and thick, isotropic and anisotropic beams, even in the presence of thermal strain. Finally, the method has been used to evaluate the natural modes and frequencies of a helicopter rotor blade, good agreement with experimental results being obtained.
ACKNOWLEDGEMENTS The author is thankful to Prof. V. Giavotto due to many invaluable suggestions in the development of this work.
CONCLUDING REMARKS Following a short review of the theory of anisotropic beams, this paper described the formulation of the space stiffness matrix of a two-node beam element, exploiting the characterisation performed with a finite element, twodimensional analysis of the cross-section. This work is the natural completion of previous papers where the approach to the analysis of composite straight, untapered and untwisted beams was described. A way to compute an array of thermal equivalent loads, suitable to be applied at the ends of a beam in a 3D frame
REFERENCES 1. 2.
3. 4. 5.
Hodges, D. H., Review of composite rotor blades modeling. A1AA J., 1990, 2(3), 561-565. Giavotto, V., Borri, H., Puccinelli, L., Caramaschi, V. and Mussi, F.,. Evaluation of section properties for hollow composite beams. Presented at 'V European Rotorcraft and Powered Lift Aircraft Forum', Paper no. 35, 1979. Bauchau, O.A., A beam theory for anisotropic material. J. Appl. Mech., 1985, 52, 416-422. Bauchau, O.A. and Hong, C.H., Finite element approach to rotor blade modeling. J. Am. Helicopter Soc., 1987, 32(1), 60-67. Stample, A.D. and Lee, S.W. Finite element model for composite
625
On the linear three-dimensional behaviour of composite beams: G. L. Ghiringhelli
6. 7.
8. 9.
10.
626
beams with arbitrary cross sectional warping. AIAA J. 1988, 26(12), 1512-1520. Cofie, E. and Bank, L.C., Analysis of thin-walled anisotropic frames by direct stiffness method. Int. J. Solids Struct., 1995, 32(2), 235249. Smith, E.C. and Chopra, I. Formulation and evaluation of an analytical model for composite box beams. Presented at '31st AIAA/ AHS/ASME/ASCE/ASC Structure, Structural Dynamics and Material Conference', Long Beach, CA, 2 - 4 April 1990. Barbero, E.J., Lopez-Anido, R. and Davalos, J.F., On the mechanics of thin-walled laminated composite beams. J. Compos. Mater., 1993, 27(8), 806-829. Giavotto, V., Borri, M., Mantegazza, P., Ghiringhelli, G. L., Caramaschi, V., Maffioli, G. C. and Mussi, F., Anisotropic beam theory and applications. Computers and Structures, 1983, 16, 303413. Kosmatka, J.B. 'Structural dynamic modeling of advanced composite propellers by the finite element method', Ph.D. Dissertation, University of California, Los Angeles, 1986.
11.
12. 13. 14. 15. 16. 17.
Hodges, D.H., Atilgan, A. R., Cesnik, C. E. S. and Fulton, M. V., On a simplified strain energy function for geometrically nonlinear behaviour of anisotropic beams. Compos. Eng., 1992, 2(5-7), 513 -526. Borri, M., Ghiringhelli, G.L. and Merlini, T., Linear analysis of naturally curved and twisted anisotropic beams. Compos. Eng., 1992, 2(5-7), 433-456. Ghiringhelli, G.L. On the thermal problem for composite beams using a finite element semi-discretisation. Composites: Part B, 1997, 28B, 483-495. Giavotto, V. Teoria elastica della trave anisotropa e non omogenea. In 'Scritti per Lucio Lazzarino', Pacini, Pisa, 1986, pp. 131-164. Ghiringhelli, G.L and Mantegazza, P., Linear straight and untwisted anisotropic beam section properties from solid finite elements. Compos. Eng., 1994, 4(12), 1225-1239. Borri, M. and Merlini, T., A large displacement formulation for anisotropic beam analysis. Meccanica, 1986, 21, 30-37. Bicos, A.S. and Springer, G.S., Design of composite boxbeam. J. Compos. Mater., 1986, 20(1), 86-109.