- Email: [email protected]
A simple and efficient curved beam element for the linear and non-linear analysis of laminated composite structures
A SIMPLE AND EFFICIENT CURVED BEAM ELEMENT FOR THE LINEAR AND NON-LINEAR ANALYSIS OF LAMINATED COMPOSITE STRUCTURES DMTRIOS KARAMANLIDIS University of Rhode Island, Department of Civil Engineering, Kingston, RI 02881, U.S.A. (Received 30 January 1987) Abstract-The paper reports on the development and application of a curved, two-noded beam finite element capable of analyzing the pre- and post-buckling behavior of laminated anisotropic structures undergoing large elastic deformations (displacements and rotations). The total number of the elemental degrees-of-freedom equals six. The formulation of the geometrically non-linear problem is performed along the lines of the ‘updated Lagrangian’ (UL) description of motion. The development of the pertinent element matrices is based on a modified version of the variational theorem due to Helhnger and Reissner which-unlike the ‘classical’ assumed displament fo~ulation-allows for stress resultants and displacements to be approximated independently from one another. In order to assess the performance of the element, a number of sample problems were investigated. Some of the numerical results obtained are presented and discussed in the final part of the paper.
1. INTRODUCIION Modern structural design has been enhanced by the introduction of advanced composite materials such as carbon-fiber-reinforced plastics. Typical application areas are aircraft and space vehicles, automobiles, pressure vessels, etc. While offering significant strength-to-wei~t advantages, composite materials have been shown to introduce mechanical couplings which require more careful analysis than is usually necessary for isotropic or orthotropic materials. For example, application of the so-called aeroelastic tailoring procedure, which consists in employing ‘unbalanced’ laminates (i.e. laminates which do not have a ply oriented at -6 for every one oriented at +6) in order to improve the aeroelastic efficiency, results in coupling between the bending and twisting behavior of the laminate, as well as between the in-plane stretching and shearing. Further, even more unique couplings may arise in using uns~me~c laminates since the stretching and bending behavior become interrelated. In fact, the use of unsymmetric laminates may result in fully anisotropic behavior such that bending, stretching, shearing and twisting are completely coupled. Another aspect that complicates the analysis of 1amiMted beams, plates or shells is due to the fact that the classical bending theory (CBT), which neglects transverse shear effects on the assumption that plane sections that were normal to a reference plane before deformation remain plane and normal after deformation, is no longer applicable. For a more realistic analysis, it is mandatory to include transverse shear effects, together with related, thickness-dependent effects such as that of rotary inertia in the vibration problem. The so-called shear deformation or first-order 623
bending theory (SDBT) removes the normalcy ass~ption by allo~ng sections initially normal to the reference plane to become non-normal (while remaining plane) after deformation. Closed-form analytical solutions for laminated structures are generally available only for relatively simple cases regarding lamination, geometry, boundary and loading conditions, etc. Thus, in general, one has to resort to numericat techniques in order to analyze practically relevant composite laminated structures. Apart from the early paper by Bufler [1], who utilized the transfer matrices method, the vast majority of the work published so far is based on the finite-element method (FIN). Surveys concerning currently available commercial FEM packages which allow for the analysis of composite laminated structures can be found, for example, in [2,3]. It appears that thus far a rather substantial effort has been devoted towards the development of finite-element models for plates (e.g. [4-71) and to a lesser extent for shells (e.g. @HO]). In contrast to this, it seems that the development of finite-element models for the analysis of free-form curved beams and arches has not been pursued as intensively as in the cases of plates and shells. To be more specific, a review of the pertinent open literature seems to indicate that the only laminated composite beam elements developed so far are those reported in [I 1-l 51. This is quite surprising if one considers that a good many composite material structural members (e.g. ring stiffeners for shells, plywood girders in roof structures, turbine and helicopter blades) can be modelled most adequately by means of beam-type elements. The aforementioned previous beam element developments are relatively limited in scope for they: (i) are applicable to certain geometries only (i.e. straight or
DIMITRIOSKARAMANLIDIS
624
circular); (ii) consider symmetric laminations only; and (iii) deal with infinitesimally small deformation (‘geometrically linear’) problems. In addition, in all these cases the classical assumed displacement method was employed On the other hand, previous work on plates (e.g. [ 16,171 and others) has demonstrated that when applied to the analysis of laminated composite structures, so-called mixed or hybrid stress formulations offer certain advantages over conventional assumed displacement elements. The research work reported herein was carried out with the intention of improving currently available modelling capabilities for free-form laminated composite curved beams and arches by proposing a novel methodology which is both simple to implement and computation~ly efficient. Precisely, a two-m&d, curved beam element (Fig. 1) with a total of six degrees-of-freedom has been developed on the basis of a modified version of the variational theorem due to Hellinger and Reissner (‘mixed-hybrid’ element). Transverse shear, non-symmetric stacking as well as large deformation effects have been inclued in the formulation.
described in terms of an orthotropic ~lationship between the Truesdell stress increments and the (linear parts of the) incremental Green strains, namely:
It is assumed sub~uently that the macromechanical properties of a single layer can be
0
0
G2
0
0
0
0
C&
0
0
0
0 c,,
G, =
4 1 - vf2. (EJE,);
E2
t2=1-v:**(E2/EI); vrtE2
c I2 =
1 - v&. (EJE,);
C, = G,z; C,, = G,,t
(2)
where E, and E2 are the Young’s moduli in principal directions 1 and 2, respectively; v,~ is the Poisson’s ratio for transverse strain in the Zdirection due to stress in the l-direction; and G,r and GL3are the shear moduli in the 1-2 and l-3 planes, respectively. In order to obtain the constitutive matrix for the laminated beam, it is necessary first to transform the constitutive equations for every layer, i.e. eqn (l), from the principal material axes (l-2-3) to the beam axes (x-y-t). Owing to the fact that the z and 3 axes
-h
T
1 (a)
(1)
wherein m is the number of layers in the laminated beam. It is, furthermore, noteworthy that the reduced stiffness coefficients C, defined w.r.t. the principal material axes (l-2-3) of a lamina (cf. Fig. 2) can be expressed in terms of the more familiar ‘engineering’ constants, i.e.
2. MATHEMATICAL FORMULATION
2.1. Macromechanical constitutive properties
C,*
G2
(P = I,. ..,m)
c
In the following, a two-~mensiona1 free-form frame is considered which undergoes arbitrarily large displacements and/or rotations (but small strains). In describing its geometrically non-linear behavior, it is presumed that: (i) the current configuration is used as reference configuration (‘updated Lagrangian’, UL, description); (ii) the finite-element mesh consists of curved beam elements which are shallow (cf. [18]) before and after deformation; (iii) the kinematics of the beam is described along the lines of SDBT (‘Timoshenko beam theory’).
C,,
(b)
Fig. I. Curved, laminated composite beam: (a) local coordinates; (b) ply coordinates in the thickness direction.
625
A simple and et&Sent curved beam element
point at the centroidal line. Then, by combining eqns (6) and (7) with each other the ‘integrated’ constitutive law for the laminated beam is obtained, i.e.
42
A22
0 0 41 42
0
0
-433
I[1 AK
(8a)
& Ae
or Au= EAiz. Fig. 2. Relation between structural axes (x-y) and material axes (l-2) for a single ply.
are parallel with each other, the following equations are obtained after the transformation
(W
In the foregoing, AN, AM and AQ are the beam stress resultants, namely: (AN; AM; AQ) =
(AU,,; Aa,, . c; Au,) dc (9)
whereas the cross-sectional properties A, are defined as follows:
@=l,...,m).
(3)
(10) The relations which allow for the calculation of the transformed reduced stiffnesses i$ from the coefficients C, can be found in any composite materials textbook (e.g. [19-211) and shall, therefore, not be presented here. Since the laminated member under consideration is a beam subject to stretching and bending in the z-x plane only
where k is the shear correction factor. 2.2. V~ri~~io~l theorem
Aayv = Aa,, = 0.
(4)
Then, from eqn (3), the following relations can be obtained by means of a standard ‘static condensation’ procedure: [%I=[?
Eq:,][iz]
(p=L...,m)
(5)
The finite-element model to be presented sequently is based on the variational equation
sub-
+ , Nf(Awd2 dx + [AN .du s
where +(AQ +AN .q).Aw
+AM,At?g
- (WEF) + (correction terms) = stationary
Et5 = E,,. In accordance with SDBT it is subsequently assumed that the normal strain As, varies linearly over the thickness, whereas the shear strain Ay,, is constant (‘Timoshenko beam theory’), namely: AsXX = Ae + 5 . AK;
Ay,, = Ay
(7)
where Ae and AK represent the stretching and bending strain increments, respectively, for a generic beam C.A.S. 29/4--o
(11)
which is a modified version of the well-known theorem due to Hellinger and Reissner (cf. [IS]). In the foregoing I denotes the distance between the two element nodes; Au = LAu, Aw, A0] T and Aa = LAN, AQ, AM J T are the incremental displacement and stress resultant vectors at a generic element point, respectively; (. . _):, denotes differentiation w.r.t. the local coordinate x; z = z(x) is a function describing the curved element centoidal line;
DIMITIUCJSK ARAMANLIDIS
626
D is the inverse material matrix of the laminated composite beam (i.e. D = E-l); and (WEF) stands for the work by externally applied forces. Inclusion of ‘correction terms’ in eqn (11) is necessary in order to prevent the accumulation of inherent linearization errors causing the obtained approximate solution to drift away from the actual one. The independent variables in eqn (1 1), namely Au and Aa, are subject to the following subsidiary constraints: (i) Au satisfies the geometric boundary and interelement continuity conditions, i.e. Au =AII (ii)
and
Au+ = Au-
(12a,b)
Au satisfies the homogeneous part of the incremental equilibrium equations, i.e. AN., = 0; AM., - AQ = 0; AQ., + (ANz.,)., = 0.
(13a,b,c)
In order to derive the pertinent element matrices from eqn (1 1), the independent variables Au and Au are approximated by means of ‘admissible’ trial functions
ment vector, 1 is a proportionality factor and r denotes the ‘fictitious load’ or ‘residual force vector’. The standard approach to solve eqn (16) consists in imposing in each solution step incremental load changes (represented by AL) and solving for Ai provided K, and r have been calculated from the data from previous solution step(s). In general, either a pure incremental or a selfcorrecting incremental or an incremental/iterative (Newton-Raphson-like) procedure can be adopted. However, all of these procedures cease to give reasonable results as soon as one approaches the vicinity of an instability point (i.e. buckling or bifurcation limit) in the load deflection curve. The reason for that is the fact that the augmented tangent stiffness matrix becomes singular (i.e. det K, = 0). In such a case it is necessary to employ a so-called continuation algorithm so that the post-critical behavior of the structure under consideration can be analyzed. The continuation algorithm utilized in this paper is the ‘energy incrementation scheme’ described in detail in a previous publication by this author [22]. In order to make the present paper reasonably selfcontained, the major features of the algorithm are outlined subsequently. First, the incremental energy by the externally applied forces is written down, namely AS = AfiT(AP + r).
Au= NAi
and
Au =P/I
(14a, b)
which are chosen so as to a priori fulfil eqns (12) and (13). In eqns (14), standard finite-element notation is used, i.e. N denotes the matrix of ‘shape functions’, Ai is the nodal displacement vector, P represents the matrix of stress trial functions and jl is the vector of generalized stress coefficients. Introduction of eqns (14) into (11) and subsequent element-by-element elimination of the j?s leads to elemental tangent stiffness matrix, namely: K, = GTH-‘G + K,.
(13
(17)
Assuming proportional load changes, the vectors AP and Ai can be scaled as shown below AP = ALP,,,;
Ai = A1 x + y
(18)
where K,x = P,,,;
K,y = r.
(19)
Then, by virtue of eqns (18) eqn (17) yields a quadratic equation in terms of A1, namely: a(AL)‘+ b(A1) + c = 0
The definition of the matrices that appear in eqns (14) and (15) is given in Appendix I.
(20)
where 3. INCREMENTAL
SOLUTION
STRATEGY
a = PL,x;
The differences between the finite-element model presented herein and the standard assumed displacement approach notwithstanding, the augmented matrix equation that needs to be solved for each incremental step has exactly the same form in both cases, i.e.
K,Aii = A1 Pier + r
(16)
where K, represents the augmented tangent stiffness matrix, P, is an arbitrarily chosen reference load vector, Ai stands for the incremental nodal displace-
b = rTx + P&y; c =rTy-AS.
(21)
It is clear from the above that by specifying in each solution step the energy increment AS (‘independent’ variable) the pertinent values for A1 (‘dependent’ variable) can be determined. Readers interested in more details about the algorithm itself or its performance in comparison with so-called ‘arc-length incrementation schemes are urged to consult [22].
A simple and efficient curved beam element 4. TEST STUDIES
The curved, laminated composite beam proposed herein has been implemented into the element-library library of the general-purpose, non-linear code MELINA [23]. Extensive numerical experiments for structures with a varying degree of complexity were carried out in order to validate the element. Quite recently, a similar development has been reported on by Whitcomb [24]. The element proposed in [24] is a Cnoded rectangle with eight degrees-offreedom altogether. The standard assumed displacement isoparametric element formulation has been utilized in conjunction with the so-called substitute shape function technique. The purpose of the latter is to avoid ‘shear locking’ usually exhibited by isoparametric elements when applied to bending problems. The performance of the element proposed in [24] was reportedly quite satisfactory for the problems considered. In a recent paper 1251,it has been shown by this author that Withcom~s element can be formulated in a mathematically consistent yet simple way by utilizing a modified version of the variational theorem due to Hellinger and Reissner 1181.Furthermore, by using the so-called degeneration concept the elemental eight degrees-of-freedom were reduced to the standard six (two translations plus one rotation per nodal point). The resulting shear-flexible beam element (subsequently called ‘hybrid strain’ element) has also been put into the MELINA library and used for comparison purposes. The examples considered herein include: (i) a simply supported beam subject to linearly dist~buted symmetric loading; (ii) a cantilever beam made from glass-reinforced plastic (GRP} material; and (iii) a hinged-hinged deep circular arch under a concentrated force at its apex. It should be mentioned that for none of these problems are reference solutions-either analytical or numerical-available from the open literature. As a matter of fact, to the best of this author’s knowledge
+-----I_-l-
Fig. 3. Simply supported cross-ply laminated beam under distributed loading.
DIMTRIOS KARAMANLIDIS
P 2d W
-Lh
=
2
L
=
loo
e., =
0.1
deflection
parameter
rotation
6 (=u/L)
79
Fig. 4. Symmetric cross-ply laminated cantilever beam: loaddeflection curves for material 1.
the development of a laminated, composite curved beam element undergoing arbitrarily large displacements and/or rotations has never been attempted before. 4.1. Simply supported beam The three-layer, cross-ply, simply supported beam subjected to a linearly distributed transverse loading (cf. Fig. 3) was considered as the first test problem. Taking advantage of the symmetry, only one half of the structure was modelled using both the hybrid stress element proposed herein and the hybrid strain element derived from Whitcomb’s formulation [24]. The obtained numerical results for different (L/h) values are shown in Table 1. (The converged solutions are printed in bold characters.) The comparison of these results clearly shows that with regard to convergence the hybrid stress element outperforms
the hybrid strain element. This result was typical for all numerical studies carried out within the course of the present investigation. 4.2. Cantilever beam In this example the large deflection behavior of a composite laminated beam subject to a concentrated force at its free end (Fig. 4) was investigated. Two Table 2. Material
(2) Material Material Material Material Material
1 2 3 4 5
27 39 46 57 64
properties of unidirectionally reinforced materials
(lcNzm2) 21.90 30.20 35.04 45.65 47.50
(kNzm’) 4.34 5.13 5.75 7.07 8.29
G
(kN/gm2) 1.61 1.90 2.13 2.63 3.09
glass-
vL7. 0.31 0.29 0.28 0.26 0.25
A simple and efficient curved beam element
629
deflection
parameter
15 (=u/L)
deflection
parameter
d (=u/L)
deflection
parameter
TJ (=r/L)
deflection
parameter
q (=w/L)
rotation
*
rotation
23
Fig. 5. Laminated composite cantilever beam: comparison of numerical results for different materials and laminations,
DIMITRIOSKARAMANLIDIS
630
I
2P
(=w/H)
deflection
parameter
7)
deflection
parameter
q (=w/H)
8 .:
deflection
parameter
77 (=w/H)
Fig. 6. Hinged-hinged laminated circular arch under a concentrated force at its apex: q-1 curves.
types of lamination were considered, i.e. a symmetric cross-ply (0°/900/900/Oo) laminate and an antisymmetric angle-ply (- 45”f +45”) laminate. The numerical values chosen for the material properties (cf. Table 2) are typical for GRP composite materials. Figure 4 shows the numerically obtained, nonlinear load-deflection curves for the structure. TWO and 10 hybrid stress and hybrid strain elements, respectively, were used to model the beam. For both elements, the obtained results agree very well with the converged solutions. The purpose of Fig. 5 is to illustrate how the non-linear structural response of the laminated composite cantilever beam is being affected by (i) the type of lamination, and (ii) the material properties of a single lamina. It is observed that the cross-ply symmetric laminate is noticeably stiffer than the angle-ply antisymmetric laminate. In contrast to
this, the influence of the different CRP materials considered appears to be relatively marginal. 4.3. Hinged-hinged cir~lar arch The third test problem studied herein was the pre- and post-buckling non-linear response of a hinged-hinged deep circular arch subject to a concentrated force at its apex (cf. Fig. 6). In order to gain some insight about the behavior of this structure, four different types of lamination for two different high-modulus, fiber-reinforced materials were considered (cf. Table 3). Due to the symmetry of the structure, only one half of the arch was modelled using four and 12 hybrid stress and hybrid strain elements, respectively. First, it is seen from Fig. 6 that the solutions obtained by these two different element models agree very well with each other. Second, in the case of an
A simple and efficient curved beam element
631
Table 3. Material properties for two typical fiber-reinforced materials &a)
(G?a)
(&
Graphite-epoxy (V,= 60%)
ls4 la4 1~8 la8
130 130 130 130
I 7 7 1
6 6 6 6
“LT 0.28 0.28 0.28 0.28
Boron-epoxy (V,= 50%)
2s4 2a4
200 200
19 19
6 6
0.23 0.23
arch made from graphite-epoxy laminae, the nonlinear response is virtually insensitive to variations in the number of layers or in the stacking sequence. The situation changes entirely in the case of an arch made from boron-epoxy laminae where the symmetric laminate appears to be considerably stiffer than the antisymmetric one. 5. CONCLUSIONS In the present paper a curved beam element has been proposed for the numerical analysis of the preand post-buckling non-linear response of laminated composite beam and frame structures. The formulation of the element is based on a modified version of the variational theorem due to Hellinger and Reissner wherein deflections and stress resultants represent independent variables (‘mixed model’). Inter-element continuous trial functions were used to approximate the displacement field, whereas the stress resultant inter-element continuity is fulfilled in an approximate sense only (‘hybrid model’). The performance of the proposed ‘hybrid stress’ element has been evaluated in comparison with a laminated composite beam element (‘hybrid strain’) derived from a two-dimensional element model proposed in a recent paper by Whitcomb [24]. Extensive test studies carried out on a number of carefully chosen examples clearly show that the proposed element: (i) represents a reliable and computationally inexpensive numerical tool for solving geometrically non-linear laminated composite beam problems; and (ii) is computationally much more efficient than the ‘hybrid strain’ element. REFERENCES
H. Bufler, Die Bestimmung des Spannungs- und Verschiebungszustandes eines geschichteten Kijrpers mit Hilfe von Ubertragungsmatrizen. 1ng Arch. 31,229-240 (1962). 0. H. Griffin, Jr, Evaluation of the finite-element software packages for stress analysis of laminated composites. Composites Technol. Rev. 4 (1982). J. Mackerle, Finite element software for stress analysis of laminated composites. Structl Anal. Syst. 2, 157-219 (1986). J. N. Reddy, Analysis of layered composite plates accounting for large deflections and transverse shear strains. Adv. Non-lin. Comput. Mech. 155-202 (1982).
Lamination +30”/-30”/-30”/+30 + 30”/- 30”/+ 30”/- 30 + 30”/- 30”/+ 30”/- 30”/- 30”/+ 30”/- 30”/+ 30” +30”/-30”/+30”/-30”/+30”/-30”/+30”/-30 + 30”/- 30”/- 30”/+ 30 +30”/-30”/+30”/-30
5. J. N. Reddy and W. C. Chao, Non-linear bending of thick rectangular, laminated composite plates. ht. J. Non-lin. Mech. 16, 291-301 (1982). 6. A. K. Noor and M. D. Mathers, Finite element analysis of anisotropic plates. Int. J. numer. Meth. Engng 11, 289-307 (1977). I. A. K. Noor ‘and J. M. Peters, Multiple-parameter reduced basis technique for bifurcation and postbuckling analyses of composite plates. Int. J. numer. Meth. Engng 19, 783-803 (1983). 8. T. Y. Chang and K. Sawamiphakdi, Large deformation analysis of laminated shells by finite element method. Comput. Struct. 13, 331-340 (1981). 9. F. Tabaddor, Large deformation of cord-reinforced multilavered’ shells. J. Pike Sci. Technol. 12. 253-267 (1979)._ 10. S. S. Murthy and R. H. Gallagher, Anisotropic cylindrical shell element based on discrete Kirchhotf theory. Int. J. numer. Meth. Engng 19, 1805-1823 (1983). 11. R. K. Gupta, A. Venkatesh and K. P. Rao, Finite element analysis of laminated anisotropic thin-walled open-section beams. Comput. Strut. 3, 19-31 (1985). 12. A. T. Chen and T. Y. Yang, Static and dynamic formulation of a symmetrically laminated beam finite element for a microcomputer. J. Composite Mater. 19, 459-475 (1985). 13. P. V. Ramana Murthy and K. P. Rao, Finite element analysis of laminated anisotropic beams of bimodulus materials. Comput. Struct. 18, 779787 (1984). 14. P. V. Ramana Murthy and K. P. Rao, Analysis of curved laminated beams of bimodulus composite materials. J. Composite Mater. 15. A. Venkatesh and K. P. Rao, A laminated anisotropic curved beam and shell stiffening finite element. Comput. Struct. 15, 197-201 (1982). 16. S. T. Mau and E. A. Witmer, Static, vibration, and thermal stress analysis1 of laminated plates and shells by the hybrid-stress finite-element method, ASRL TR 169-2, MIT (1972). 17. H. Dittrich, Berechnung von Sandwichplatten mittels gemischter finiter Elemente. Doctoral Dissertation, University of Stuttgart (1979). 18. K. Washizu. Variational Methods in Elasticitv and Plasticity, 3rd edn. Pergamon Press, Oxford (1982). 19. R. -M. Jones, Mechanics of Composite Materials. McGraw-Hill. New York (1975). of an Anisotropic 20. S. G. Lekhnitskii, Theory ofEla&ity Body. Mir, Moscow (1981). 21. L. R. Calcote, The Analysis of Laminated Composite Structures. Van Nostrand Reinhold, New York (1969). 22. D. Karamanlidis, A note on incrementation schemes in nonlinear structural analysis. Engng Anal. 3,225-234 (1986). element 23. D. Karamanlidis, MELINA-Multipurpose language for incremental nonlinear analysis (user’s manual), Internal Report, Civil Engineering Department, University of Rhode Island (1983). 24. J. D. Whitcomb, A simple rectangular element for two-dimensional analysis of laminated composites. Comput. Struct. 22, 387-393 (1986).
632
DIMITRIOS KARAMANLIDIS
25. D. Karamanlidis and R. Jasti, On the variational foundations of the substitute shape function technique. Cornput. Srrucf. 26, 1041-1042 (1987).
APPENDIX
l-5
[
1
P=
-Z.x -z
0
H=
1-r 00r1;<=x/l 1 0
0
PTDPdx
0
0 -1
0
r
r
0 0
‘0
0
1 0 0 0
0 N/I
00 0
1 I
0 0
00
-N/I
0 0
0
x1 01 0 0
0
001
I: MATRICES
K, =
0
fI
NII symmetric
1. INTRODUCIION Modern structural design has been enhanced by the introduction of advanced composite materials such as carbon-fiber-reinforced plastics. Typical application areas are aircraft and space vehicles, automobiles, pressure vessels, etc. While offering significant strength-to-wei~t advantages, composite materials have been shown to introduce mechanical couplings which require more careful analysis than is usually necessary for isotropic or orthotropic materials. For example, application of the so-called aeroelastic tailoring procedure, which consists in employing ‘unbalanced’ laminates (i.e. laminates which do not have a ply oriented at -6 for every one oriented at +6) in order to improve the aeroelastic efficiency, results in coupling between the bending and twisting behavior of the laminate, as well as between the in-plane stretching and shearing. Further, even more unique couplings may arise in using uns~me~c laminates since the stretching and bending behavior become interrelated. In fact, the use of unsymmetric laminates may result in fully anisotropic behavior such that bending, stretching, shearing and twisting are completely coupled. Another aspect that complicates the analysis of 1amiMted beams, plates or shells is due to the fact that the classical bending theory (CBT), which neglects transverse shear effects on the assumption that plane sections that were normal to a reference plane before deformation remain plane and normal after deformation, is no longer applicable. For a more realistic analysis, it is mandatory to include transverse shear effects, together with related, thickness-dependent effects such as that of rotary inertia in the vibration problem. The so-called shear deformation or first-order 623
bending theory (SDBT) removes the normalcy ass~ption by allo~ng sections initially normal to the reference plane to become non-normal (while remaining plane) after deformation. Closed-form analytical solutions for laminated structures are generally available only for relatively simple cases regarding lamination, geometry, boundary and loading conditions, etc. Thus, in general, one has to resort to numericat techniques in order to analyze practically relevant composite laminated structures. Apart from the early paper by Bufler [1], who utilized the transfer matrices method, the vast majority of the work published so far is based on the finite-element method (FIN). Surveys concerning currently available commercial FEM packages which allow for the analysis of composite laminated structures can be found, for example, in [2,3]. It appears that thus far a rather substantial effort has been devoted towards the development of finite-element models for plates (e.g. [4-71) and to a lesser extent for shells (e.g. @HO]). In contrast to this, it seems that the development of finite-element models for the analysis of free-form curved beams and arches has not been pursued as intensively as in the cases of plates and shells. To be more specific, a review of the pertinent open literature seems to indicate that the only laminated composite beam elements developed so far are those reported in [I 1-l 51. This is quite surprising if one considers that a good many composite material structural members (e.g. ring stiffeners for shells, plywood girders in roof structures, turbine and helicopter blades) can be modelled most adequately by means of beam-type elements. The aforementioned previous beam element developments are relatively limited in scope for they: (i) are applicable to certain geometries only (i.e. straight or
DIMITRIOSKARAMANLIDIS
624
circular); (ii) consider symmetric laminations only; and (iii) deal with infinitesimally small deformation (‘geometrically linear’) problems. In addition, in all these cases the classical assumed displacement method was employed On the other hand, previous work on plates (e.g. [ 16,171 and others) has demonstrated that when applied to the analysis of laminated composite structures, so-called mixed or hybrid stress formulations offer certain advantages over conventional assumed displacement elements. The research work reported herein was carried out with the intention of improving currently available modelling capabilities for free-form laminated composite curved beams and arches by proposing a novel methodology which is both simple to implement and computation~ly efficient. Precisely, a two-m&d, curved beam element (Fig. 1) with a total of six degrees-of-freedom has been developed on the basis of a modified version of the variational theorem due to Hellinger and Reissner (‘mixed-hybrid’ element). Transverse shear, non-symmetric stacking as well as large deformation effects have been inclued in the formulation.
described in terms of an orthotropic ~lationship between the Truesdell stress increments and the (linear parts of the) incremental Green strains, namely:
It is assumed sub~uently that the macromechanical properties of a single layer can be
0
0
G2
0
0
0
0
C&
0
0
0
0 c,,
G, =
4 1 - vf2. (EJE,);
E2
t2=1-v:**(E2/EI); vrtE2
c I2 =
1 - v&. (EJE,);
C, = G,z; C,, = G,,t
(2)
where E, and E2 are the Young’s moduli in principal directions 1 and 2, respectively; v,~ is the Poisson’s ratio for transverse strain in the Zdirection due to stress in the l-direction; and G,r and GL3are the shear moduli in the 1-2 and l-3 planes, respectively. In order to obtain the constitutive matrix for the laminated beam, it is necessary first to transform the constitutive equations for every layer, i.e. eqn (l), from the principal material axes (l-2-3) to the beam axes (x-y-t). Owing to the fact that the z and 3 axes
-h
T
1 (a)
(1)
wherein m is the number of layers in the laminated beam. It is, furthermore, noteworthy that the reduced stiffness coefficients C, defined w.r.t. the principal material axes (l-2-3) of a lamina (cf. Fig. 2) can be expressed in terms of the more familiar ‘engineering’ constants, i.e.
2. MATHEMATICAL FORMULATION
2.1. Macromechanical constitutive properties
C,*
G2
(P = I,. ..,m)
c
In the following, a two-~mensiona1 free-form frame is considered which undergoes arbitrarily large displacements and/or rotations (but small strains). In describing its geometrically non-linear behavior, it is presumed that: (i) the current configuration is used as reference configuration (‘updated Lagrangian’, UL, description); (ii) the finite-element mesh consists of curved beam elements which are shallow (cf. [18]) before and after deformation; (iii) the kinematics of the beam is described along the lines of SDBT (‘Timoshenko beam theory’).
C,,
(b)
Fig. I. Curved, laminated composite beam: (a) local coordinates; (b) ply coordinates in the thickness direction.
625
A simple and et&Sent curved beam element
point at the centroidal line. Then, by combining eqns (6) and (7) with each other the ‘integrated’ constitutive law for the laminated beam is obtained, i.e.
42
A22
0 0 41 42
0
0
-433
I[1 AK
(8a)
& Ae
or Au= EAiz. Fig. 2. Relation between structural axes (x-y) and material axes (l-2) for a single ply.
are parallel with each other, the following equations are obtained after the transformation
(W
In the foregoing, AN, AM and AQ are the beam stress resultants, namely: (AN; AM; AQ) =
(AU,,; Aa,, . c; Au,) dc (9)
whereas the cross-sectional properties A, are defined as follows:
@=l,...,m).
(3)
(10) The relations which allow for the calculation of the transformed reduced stiffnesses i$ from the coefficients C, can be found in any composite materials textbook (e.g. [19-211) and shall, therefore, not be presented here. Since the laminated member under consideration is a beam subject to stretching and bending in the z-x plane only
where k is the shear correction factor. 2.2. V~ri~~io~l theorem
Aayv = Aa,, = 0.
(4)
Then, from eqn (3), the following relations can be obtained by means of a standard ‘static condensation’ procedure: [%I=[?
Eq:,][iz]
(p=L...,m)
(5)
The finite-element model to be presented sequently is based on the variational equation
sub-
+ , Nf(Awd2 dx + [AN .du s
where +(AQ +AN .q).Aw
+AM,At?g
- (WEF) + (correction terms) = stationary
Et5 = E,,. In accordance with SDBT it is subsequently assumed that the normal strain As, varies linearly over the thickness, whereas the shear strain Ay,, is constant (‘Timoshenko beam theory’), namely: AsXX = Ae + 5 . AK;
Ay,, = Ay
(7)
where Ae and AK represent the stretching and bending strain increments, respectively, for a generic beam C.A.S. 29/4--o
(11)
which is a modified version of the well-known theorem due to Hellinger and Reissner (cf. [IS]). In the foregoing I denotes the distance between the two element nodes; Au = LAu, Aw, A0] T and Aa = LAN, AQ, AM J T are the incremental displacement and stress resultant vectors at a generic element point, respectively; (. . _):, denotes differentiation w.r.t. the local coordinate x; z = z(x) is a function describing the curved element centoidal line;
DIMITIUCJSK ARAMANLIDIS
626
D is the inverse material matrix of the laminated composite beam (i.e. D = E-l); and (WEF) stands for the work by externally applied forces. Inclusion of ‘correction terms’ in eqn (11) is necessary in order to prevent the accumulation of inherent linearization errors causing the obtained approximate solution to drift away from the actual one. The independent variables in eqn (1 1), namely Au and Aa, are subject to the following subsidiary constraints: (i) Au satisfies the geometric boundary and interelement continuity conditions, i.e. Au =AII (ii)
and
Au+ = Au-
(12a,b)
Au satisfies the homogeneous part of the incremental equilibrium equations, i.e. AN., = 0; AM., - AQ = 0; AQ., + (ANz.,)., = 0.
(13a,b,c)
In order to derive the pertinent element matrices from eqn (1 1), the independent variables Au and Au are approximated by means of ‘admissible’ trial functions
ment vector, 1 is a proportionality factor and r denotes the ‘fictitious load’ or ‘residual force vector’. The standard approach to solve eqn (16) consists in imposing in each solution step incremental load changes (represented by AL) and solving for Ai provided K, and r have been calculated from the data from previous solution step(s). In general, either a pure incremental or a selfcorrecting incremental or an incremental/iterative (Newton-Raphson-like) procedure can be adopted. However, all of these procedures cease to give reasonable results as soon as one approaches the vicinity of an instability point (i.e. buckling or bifurcation limit) in the load deflection curve. The reason for that is the fact that the augmented tangent stiffness matrix becomes singular (i.e. det K, = 0). In such a case it is necessary to employ a so-called continuation algorithm so that the post-critical behavior of the structure under consideration can be analyzed. The continuation algorithm utilized in this paper is the ‘energy incrementation scheme’ described in detail in a previous publication by this author [22]. In order to make the present paper reasonably selfcontained, the major features of the algorithm are outlined subsequently. First, the incremental energy by the externally applied forces is written down, namely AS = AfiT(AP + r).
Au= NAi
and
Au =P/I
(14a, b)
which are chosen so as to a priori fulfil eqns (12) and (13). In eqns (14), standard finite-element notation is used, i.e. N denotes the matrix of ‘shape functions’, Ai is the nodal displacement vector, P represents the matrix of stress trial functions and jl is the vector of generalized stress coefficients. Introduction of eqns (14) into (11) and subsequent element-by-element elimination of the j?s leads to elemental tangent stiffness matrix, namely: K, = GTH-‘G + K,.
(13
(17)
Assuming proportional load changes, the vectors AP and Ai can be scaled as shown below AP = ALP,,,;
Ai = A1 x + y
(18)
where K,x = P,,,;
K,y = r.
(19)
Then, by virtue of eqns (18) eqn (17) yields a quadratic equation in terms of A1, namely: a(AL)‘+ b(A1) + c = 0
The definition of the matrices that appear in eqns (14) and (15) is given in Appendix I.
(20)
where 3. INCREMENTAL
SOLUTION
STRATEGY
a = PL,x;
The differences between the finite-element model presented herein and the standard assumed displacement approach notwithstanding, the augmented matrix equation that needs to be solved for each incremental step has exactly the same form in both cases, i.e.
K,Aii = A1 Pier + r
(16)
where K, represents the augmented tangent stiffness matrix, P, is an arbitrarily chosen reference load vector, Ai stands for the incremental nodal displace-
b = rTx + P&y; c =rTy-AS.
(21)
It is clear from the above that by specifying in each solution step the energy increment AS (‘independent’ variable) the pertinent values for A1 (‘dependent’ variable) can be determined. Readers interested in more details about the algorithm itself or its performance in comparison with so-called ‘arc-length incrementation schemes are urged to consult [22].
A simple and efficient curved beam element 4. TEST STUDIES
The curved, laminated composite beam proposed herein has been implemented into the element-library library of the general-purpose, non-linear code MELINA [23]. Extensive numerical experiments for structures with a varying degree of complexity were carried out in order to validate the element. Quite recently, a similar development has been reported on by Whitcomb [24]. The element proposed in [24] is a Cnoded rectangle with eight degrees-offreedom altogether. The standard assumed displacement isoparametric element formulation has been utilized in conjunction with the so-called substitute shape function technique. The purpose of the latter is to avoid ‘shear locking’ usually exhibited by isoparametric elements when applied to bending problems. The performance of the element proposed in [24] was reportedly quite satisfactory for the problems considered. In a recent paper 1251,it has been shown by this author that Withcom~s element can be formulated in a mathematically consistent yet simple way by utilizing a modified version of the variational theorem due to Hellinger and Reissner 1181.Furthermore, by using the so-called degeneration concept the elemental eight degrees-of-freedom were reduced to the standard six (two translations plus one rotation per nodal point). The resulting shear-flexible beam element (subsequently called ‘hybrid strain’ element) has also been put into the MELINA library and used for comparison purposes. The examples considered herein include: (i) a simply supported beam subject to linearly dist~buted symmetric loading; (ii) a cantilever beam made from glass-reinforced plastic (GRP} material; and (iii) a hinged-hinged deep circular arch under a concentrated force at its apex. It should be mentioned that for none of these problems are reference solutions-either analytical or numerical-available from the open literature. As a matter of fact, to the best of this author’s knowledge
+-----I_-l-
Fig. 3. Simply supported cross-ply laminated beam under distributed loading.
DIMTRIOS KARAMANLIDIS
P 2d W
-Lh
=
2
L
=
loo
e., =
0.1
deflection
parameter
rotation
6 (=u/L)
79
Fig. 4. Symmetric cross-ply laminated cantilever beam: loaddeflection curves for material 1.
the development of a laminated, composite curved beam element undergoing arbitrarily large displacements and/or rotations has never been attempted before. 4.1. Simply supported beam The three-layer, cross-ply, simply supported beam subjected to a linearly distributed transverse loading (cf. Fig. 3) was considered as the first test problem. Taking advantage of the symmetry, only one half of the structure was modelled using both the hybrid stress element proposed herein and the hybrid strain element derived from Whitcomb’s formulation [24]. The obtained numerical results for different (L/h) values are shown in Table 1. (The converged solutions are printed in bold characters.) The comparison of these results clearly shows that with regard to convergence the hybrid stress element outperforms
the hybrid strain element. This result was typical for all numerical studies carried out within the course of the present investigation. 4.2. Cantilever beam In this example the large deflection behavior of a composite laminated beam subject to a concentrated force at its free end (Fig. 4) was investigated. Two Table 2. Material
(2) Material Material Material Material Material
1 2 3 4 5
27 39 46 57 64
properties of unidirectionally reinforced materials
(lcNzm2) 21.90 30.20 35.04 45.65 47.50
(kNzm’) 4.34 5.13 5.75 7.07 8.29
G
(kN/gm2) 1.61 1.90 2.13 2.63 3.09
glass-
vL7. 0.31 0.29 0.28 0.26 0.25
A simple and efficient curved beam element
629
deflection
parameter
15 (=u/L)
deflection
parameter
d (=u/L)
deflection
parameter
TJ (=r/L)
deflection
parameter
q (=w/L)
rotation
*
rotation
23
Fig. 5. Laminated composite cantilever beam: comparison of numerical results for different materials and laminations,
DIMITRIOSKARAMANLIDIS
630
I
2P
(=w/H)
deflection
parameter
7)
deflection
parameter
q (=w/H)
8 .:
deflection
parameter
77 (=w/H)
Fig. 6. Hinged-hinged laminated circular arch under a concentrated force at its apex: q-1 curves.
types of lamination were considered, i.e. a symmetric cross-ply (0°/900/900/Oo) laminate and an antisymmetric angle-ply (- 45”f +45”) laminate. The numerical values chosen for the material properties (cf. Table 2) are typical for GRP composite materials. Figure 4 shows the numerically obtained, nonlinear load-deflection curves for the structure. TWO and 10 hybrid stress and hybrid strain elements, respectively, were used to model the beam. For both elements, the obtained results agree very well with the converged solutions. The purpose of Fig. 5 is to illustrate how the non-linear structural response of the laminated composite cantilever beam is being affected by (i) the type of lamination, and (ii) the material properties of a single lamina. It is observed that the cross-ply symmetric laminate is noticeably stiffer than the angle-ply antisymmetric laminate. In contrast to
this, the influence of the different CRP materials considered appears to be relatively marginal. 4.3. Hinged-hinged cir~lar arch The third test problem studied herein was the pre- and post-buckling non-linear response of a hinged-hinged deep circular arch subject to a concentrated force at its apex (cf. Fig. 6). In order to gain some insight about the behavior of this structure, four different types of lamination for two different high-modulus, fiber-reinforced materials were considered (cf. Table 3). Due to the symmetry of the structure, only one half of the arch was modelled using four and 12 hybrid stress and hybrid strain elements, respectively. First, it is seen from Fig. 6 that the solutions obtained by these two different element models agree very well with each other. Second, in the case of an
A simple and efficient curved beam element
631
Table 3. Material properties for two typical fiber-reinforced materials &a)
(G?a)
(&
Graphite-epoxy (V,= 60%)
ls4 la4 1~8 la8
130 130 130 130
I 7 7 1
6 6 6 6
“LT 0.28 0.28 0.28 0.28
Boron-epoxy (V,= 50%)
2s4 2a4
200 200
19 19
6 6
0.23 0.23
arch made from graphite-epoxy laminae, the nonlinear response is virtually insensitive to variations in the number of layers or in the stacking sequence. The situation changes entirely in the case of an arch made from boron-epoxy laminae where the symmetric laminate appears to be considerably stiffer than the antisymmetric one. 5. CONCLUSIONS In the present paper a curved beam element has been proposed for the numerical analysis of the preand post-buckling non-linear response of laminated composite beam and frame structures. The formulation of the element is based on a modified version of the variational theorem due to Hellinger and Reissner wherein deflections and stress resultants represent independent variables (‘mixed model’). Inter-element continuous trial functions were used to approximate the displacement field, whereas the stress resultant inter-element continuity is fulfilled in an approximate sense only (‘hybrid model’). The performance of the proposed ‘hybrid stress’ element has been evaluated in comparison with a laminated composite beam element (‘hybrid strain’) derived from a two-dimensional element model proposed in a recent paper by Whitcomb [24]. Extensive test studies carried out on a number of carefully chosen examples clearly show that the proposed element: (i) represents a reliable and computationally inexpensive numerical tool for solving geometrically non-linear laminated composite beam problems; and (ii) is computationally much more efficient than the ‘hybrid strain’ element. REFERENCES
H. Bufler, Die Bestimmung des Spannungs- und Verschiebungszustandes eines geschichteten Kijrpers mit Hilfe von Ubertragungsmatrizen. 1ng Arch. 31,229-240 (1962). 0. H. Griffin, Jr, Evaluation of the finite-element software packages for stress analysis of laminated composites. Composites Technol. Rev. 4 (1982). J. Mackerle, Finite element software for stress analysis of laminated composites. Structl Anal. Syst. 2, 157-219 (1986). J. N. Reddy, Analysis of layered composite plates accounting for large deflections and transverse shear strains. Adv. Non-lin. Comput. Mech. 155-202 (1982).
Lamination +30”/-30”/-30”/+30 + 30”/- 30”/+ 30”/- 30 + 30”/- 30”/+ 30”/- 30”/- 30”/+ 30”/- 30”/+ 30” +30”/-30”/+30”/-30”/+30”/-30”/+30”/-30 + 30”/- 30”/- 30”/+ 30 +30”/-30”/+30”/-30
5. J. N. Reddy and W. C. Chao, Non-linear bending of thick rectangular, laminated composite plates. ht. J. Non-lin. Mech. 16, 291-301 (1982). 6. A. K. Noor and M. D. Mathers, Finite element analysis of anisotropic plates. Int. J. numer. Meth. Engng 11, 289-307 (1977). I. A. K. Noor ‘and J. M. Peters, Multiple-parameter reduced basis technique for bifurcation and postbuckling analyses of composite plates. Int. J. numer. Meth. Engng 19, 783-803 (1983). 8. T. Y. Chang and K. Sawamiphakdi, Large deformation analysis of laminated shells by finite element method. Comput. Struct. 13, 331-340 (1981). 9. F. Tabaddor, Large deformation of cord-reinforced multilavered’ shells. J. Pike Sci. Technol. 12. 253-267 (1979)._ 10. S. S. Murthy and R. H. Gallagher, Anisotropic cylindrical shell element based on discrete Kirchhotf theory. Int. J. numer. Meth. Engng 19, 1805-1823 (1983). 11. R. K. Gupta, A. Venkatesh and K. P. Rao, Finite element analysis of laminated anisotropic thin-walled open-section beams. Comput. Strut. 3, 19-31 (1985). 12. A. T. Chen and T. Y. Yang, Static and dynamic formulation of a symmetrically laminated beam finite element for a microcomputer. J. Composite Mater. 19, 459-475 (1985). 13. P. V. Ramana Murthy and K. P. Rao, Finite element analysis of laminated anisotropic beams of bimodulus materials. Comput. Struct. 18, 779787 (1984). 14. P. V. Ramana Murthy and K. P. Rao, Analysis of curved laminated beams of bimodulus composite materials. J. Composite Mater. 15. A. Venkatesh and K. P. Rao, A laminated anisotropic curved beam and shell stiffening finite element. Comput. Struct. 15, 197-201 (1982). 16. S. T. Mau and E. A. Witmer, Static, vibration, and thermal stress analysis1 of laminated plates and shells by the hybrid-stress finite-element method, ASRL TR 169-2, MIT (1972). 17. H. Dittrich, Berechnung von Sandwichplatten mittels gemischter finiter Elemente. Doctoral Dissertation, University of Stuttgart (1979). 18. K. Washizu. Variational Methods in Elasticitv and Plasticity, 3rd edn. Pergamon Press, Oxford (1982). 19. R. -M. Jones, Mechanics of Composite Materials. McGraw-Hill. New York (1975). of an Anisotropic 20. S. G. Lekhnitskii, Theory ofEla&ity Body. Mir, Moscow (1981). 21. L. R. Calcote, The Analysis of Laminated Composite Structures. Van Nostrand Reinhold, New York (1969). 22. D. Karamanlidis, A note on incrementation schemes in nonlinear structural analysis. Engng Anal. 3,225-234 (1986). element 23. D. Karamanlidis, MELINA-Multipurpose language for incremental nonlinear analysis (user’s manual), Internal Report, Civil Engineering Department, University of Rhode Island (1983). 24. J. D. Whitcomb, A simple rectangular element for two-dimensional analysis of laminated composites. Comput. Struct. 22, 387-393 (1986).
632
DIMITRIOS KARAMANLIDIS
25. D. Karamanlidis and R. Jasti, On the variational foundations of the substitute shape function technique. Cornput. Srrucf. 26, 1041-1042 (1987).
APPENDIX
l-5
[
1
P=
-Z.x -z
0
H=
1-r 00r1;<=x/l 1 0
0
PTDPdx
0
0 -1
0
r
r
0 0
‘0
0
1 0 0 0
0 N/I
00 0
1 I
0 0
00
-N/I
0 0
0
x1 01 0 0
0
001
I: MATRICES
K, =
0
fI
NII symmetric