Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 2, 39, 2001 NUMERICAL ANALYSIS OF THICK, MULTI-LAYERED COMPOSITE GIRDERS USING HYBRID-STRESS FINITE ELEMENTS Jerzy Gołaś Faculty of Civil andEnvironment Engineering, University of Technology andAgriculture, Bydgoszcz e-mail: tomaszj@mail.atr.bydgoszcz.pl This paper presents an analysis of numerical solutions for multi-layer composite girders under static loading. In the algorithms of presented solutions the hybrid-stress model of the finite element method based on Reissner’s modified variational functional was used. Two versions of numerical programs were developed for an N-layer finite element. The programs were tested on the example of a three-layer laminated beam of carbon fibre-reinforced epoxy composite, as well as on the example of numerical solution for a cantilevered plate. In addition, numerical examples concerning laminated glued timber beams and a retainingwall of reinforced soil were analysed. Key words: multi-layer composite, transverse inhomogeneity, hybrid- stress finite element, numerical analysis 1. Introduction Multi-layer composite girders are characteristic by large transverse defor- mability (warping) of cross section due to variable mechanical properties of individual layers. It is, therefore, necessary to take into account the effect of transverse shear in the strength analysis of these girders and, often, the influ- ence of elongation of normal elements as well, see Lo et al. (1977), Mau et al. (1972), Pagano (1969), Pian and Chen (1982, 1983), Rikards (1988), Spilker et al. (1977), Spilker (1980). This is particularly necessary in cases of deep girders of the depth/span ratio h/L > 0.25. This can, however, also apply to more slender girders, in which the ratios of the corresponding moduli of elasticity are large, i.e. Eix/E i+1 x > 100 (see Gołaś, 1995, 1997). Thus correct 284 J.Gołaś determination of stresses anddisplacements inmulti-layer girders becomes po- ssible only when it is based on the strength analysis, which takes into account the effect of the transverse shear and elongation of the normal elements. What is presented below is the algorithm, on the basis of which a statical computation program of multi-layer cylindrically bent girders was developed. This is followedbyexamples of actual solutions.Thehybrid stressmodel of the finite elementmethod, based on themodifiedReissner’s variational functional was used to formulate algebraic equations. Such a model was proposed and analysed by several authors, namely:Mau et al. (1972), Pian andChen (1982, 1983), Rikards (1988), Spilker et al. (1977), Spilker (1980). 2. Formulation of the problem The subject of the following analysis are thick, multi-layer composite gir- ders composed of N layers of diversified stiffness, bent to a cylindrical surface in plane x,z. The stiffness of individual layers across their thickness is con- stant, but can vary along the girder span – in the direction of the x axis. The material of each layer is ortothropic and linear-elastic. The girder is subjec- ted to known static loading, applied to the external surfaces z = ±h/2 and mass loading inside the volume. In addition, suitable geometrical boundary conditions are imposed on the parts of girder edges. The girder is subdivided into an arbitrary number ne of laminar elements by cross-sections perpendicular to the longitudinal axis x. Typical laminar finite element is shown in Fig.1. At the edges of the elements contact at x=0, l thedisplacementboundary conditions u = ũ are satisfied. The modified Reissner variational functional assumes the following form (see Spilker et al., 1977; Spilker, 1980; Rikards, 1988) ΠmR(u,σ)= ∑ ne { N∑ i=1 [ − 1 2 ∫ Vni σ i⊤ S i σ i dV + (2.1) + ∫ Vni σ i⊤ e i(ui) dV − ∫ Vni f̃ ⊤ i u i dV − ∫ Sσ ni t̃ ⊤ u i dS ]} where σi,ei,ui – stress, strain and displacement vectors, respectively Numerical analysis of thick, multi-layered composite... 285 Fig. 1. Multi-layer finite element. Geometry and numeration of layers S i – flexibility coefficientmatrix of the ortothropicmaterial in the ith layer Vni – volume of the ith layer in nth element Sσni – surface on which the external loadings t̃ are imposed f̃i – mass force vector. The stress field σi ⊤ = [σix,σ i z,σ i xz], displacement field u i⊤ = [ui,wi] and strain field ei ⊤ = [εix,ε i z,γ i xz] are approximated by following functions (Spilker, 1980) σ i =Pi   β i β i β i+1   =Piβi ui =Nidi ei =Bidi (2.2) in which Pi,Ni,Bi are appropriate matrices of the approximation functions, di is the nodal displacements vector and βi – vector containing nodal stress parameters βi assigned to the inside of a given layer, aswell as theparameters β i and β i+1 connectedwith its bottomand top surfaces.The stress conditions of continuity at the contact surfaces of the ith and (i+1)th layers σiz(zi+1)=σ i+1 z (zi+1) σ i xz(zi+1)=σ i+1 xz (zi+1) (2.3) are carried out by the equality of the corresponding parameters β i and β i+1 286 J.Gołaś of the surfaces in contact. If, at the external surfaces of the extreme layers 1 and N, the static conditions are imposed in the form σ1xz(z1)=σ 1 z(z1)= 0 σ N xz(zN+1)= 0 (2.4) then β 1 =0 β N+1 σxz =0 (2.5) where the lower index in (2.5)2 denotes the choice of these parameters β N+1 only, which are connected with the shear stress components. Substituting relationships (2.2) into functional (2.1), after suitable sum- mation across all layers i=1,2, ...,N, we obtain Πmc(de,β)= ∑ ne Πemc = ∑ ne ( − 1 2 β ⊤ Hβ+β⊤Qde−d ⊤ eFe ) (2.6) where the matrices H, Q and Fe correspond to the given layered element and are appropriately composed of the following sub-matrices H i = ∫ Vni P i⊤ S i P i dV Q i = ∫ Vni P i⊤ B i dV (2.7) F i = ∫ Vni N i⊤ f̃ i dV + ∫ Sσ ni N i⊤ t̃ i dS and the vectors β and de the components of which correspond to parameters of individual layers of the element. The nodes for the description of displace- ments are located only at the contacting edges of the elements x= const. From the condition of stationariness of functional (2.6) with respect to the mutually independent parameters β and de, we obtain β=H−1Qde (2.8) and ∑ ne (kde−Fe)=0 or Kd=F (2.9) where k=Q⊤H−1Q (2.10) denotes the stiffness matrix of the layered element. Numerical analysis of thick, multi-layered composite... 287 The displacements are determined from displacement algebraic equations of equilibrium (2.9) referring to the global system, following by stresses com- puted in individual layers of the system using relationships (2.8) and (2.2)1. Twoversionsof thecomputerprogramfor thefiniteelementof N layers are presented in the paper. The first having the number of stress approximation parameters nβ = 14N − 5 and the number of displacement approximation parameters nu =8N+4, and the second 14N−5 and 4N+4, respectively. In both versions of the program the two-dimensional stress field in each layer of the element is approximated by the relationship σ i(x,z) ︸ ︷︷ ︸ (3×1) =Pi(x,z) ︸ ︷︷ ︸ (3×19) βi︸︷︷︸ (19×1) (2.11) where the polynomial functions matrix Pi(x,z) contains – relevantly for the components: σix – a third degree approximation in the two directions of the x and z coordinates; σiz – a linear approximation in the x direction and a fifth degree approximation in the z direction; σixz – a quadratic approximation in the x direction and a fourth degree approximation in the z direction. For example, P i1,14,P i 2,14 and P i 3,14 elements of the P i(x,z) matrix are as follows P i1,14 =Fi(z)x 3 P i2,14 = 3 10 Ki(z)x P i 3,14 = 3 20 Ri(z)x 2 (2.12) where Fi(z)= 1 5 (z3i +4z 2 izi+1+4ziz 2 i+1+z 3 i+1)− 3 10 (3z2i +4zizi+1+3z 2 i+1)z+z 3 Ki(z)= 2(z 3 iz 2 i+1+z 2 iz 3 i+1)−z(4z 3 izi+1+7z 2 iz 2 i+1+4ziz 3 i+1)+ +2z2(z3i +4z 2 izi+1+4ziz 2 i+1+z 3 i+1)−z 3(3z2i +4zizi+1+3z 2 i+1)+z 5 (2.13) Ri(z)= (4z 3 izi+1+7z 2 iz 2 i+1+4ziz 3 i+1)−4z(z 3 i +4z 2 izi+1+4ziz 2 i+1+z 3 i+1)+ +3z2(3z2i +4zizi+1+3z 2 i+1)−5z 4 On the other hand the parameter vector βi has the form βi = [β i 1,β i 2, ...,β i 5;β i 1,β i 2, ...,β i 9;β i+1 1 ,β i+1 2 , ...,β i+1 5 ] ⊤ In addition, the adopted description of stresses (2.11) satisfies homogeneous equilibrium conditions in the area of the layer and static conditions on the external surfaces of the extreme layers. 288 J.Gołaś In order to describe the displacement field ui(x,z) = [ui,wi]⊤ in a typical layer element, three nodes were assumed in the first version of the program along the edges x=0, l and twonodes in the secondversion.At each node two translational degrees of freedom were assumed. Thus, in each of the versions we have nβ =9, nu =12 and nβ =9, nu =8 independent stress parameters and displacement degrees of freedom per layer, respectively. 3. Numerical analysis and concluding remarks The developed program was tested on the example of a three-layer plate strip considered in papers by Lo et al. (1977), Spilker (1980), Rikards (1988), bent to a cylindrical surface, for which exact analytical solutions, obtained on the basis of the theory of elasticity, are known (see Pagano, 1969). For a strip of h/L ratio equal 0.25, divided into 12 finite elements in both versions of the developed program, the relative error of the deflections reduced todimensionless values w(L/2,0)didnot exceed 1.3%.Ontheotherhand, the greatest deviations of the reducedvalues of the stress field in the cross-sections σx(L/2,z), σz(L/2,z), σxz(0,z) were: for the first version of the program, correspondingly 4.8, 2.1 and 8.2% – with subdivision into 12 elements, and 3.9, 1.6 and 5.6% – with subdivision into 24 elements. For the second version of the program thedeviationswere 5 up to 8%greater than in the first version. In view of higher precision of the solutions, further numerical analysis shall be carried out using the first version of the program. Fig. 2. Cantilever plate. Statical scheme and grid of elements The program was also verified on the example of a FEM displacement model of cantilever plate, presented byKleiber (1989). In this case the homo- geneous cantilever plate of dimensions 2.0m×10.5m, was supported on the Numerical analysis of thick, multi-layered composite... 289 edge x= 0 and loaded by a concentrated force P at the point (10.5,0) on the free edge x=L (Fig.2). The values of the vertical displacements w(L,0) at the point of application of the force P with subdivision into 10 and 14 four-layer elements differed from the solutions in Kleiber (1989) by 6.0 and 5.9%, respectively. The greatest differences in the stresses σx and σxz in the cross-section x = 0.375m (in order to compare the solutions of both FEM models) did not exceed 5.8%. Full picture of the obtained stress field u(x,z) and σx stress component are shown in Fig.3. The displacements presented there are magnified ten-fold with respect to the geometric dimensions of the plate. Intensity of shading of the surface corresponds to the increase of the stress values. wmax =0.14192m; umax =0.01968m; σxmax =7634.591kPa Fig. 3. Displacement field u(x,z) and normal stress component σ x (x,z) for isotropic cantilever plate loaded by force P =500kN (Fig.2). Material constants: E=2.1 ·106kPa, ν =0.25 Using the developed program, bending of a glued timber laminated beam was analysed, wheremechanical properties of individual layers were subjected to degradation during operation. Results of the numerical computations of the simply supported and fixed-end beam, both composed of seven layers, the properties of which changed in time from the state I to state III (according to Table 1) are presented below.Thebeam spanwas L=12m, depth h=1.0m and theuniformlydistributed load qz =0.100MN/mwas applied to theupper surface z =0.5m. The beam was subdivided along the span into 40 laminar elements. The results of computations for beams: (a) simply supported and (b) fixed at both ends are presented in Fig.4 and Fig.5, respectively. They present the deflection of the beam axis and distribution of the normal stresses σx for each of the examined states I, II, and III. It was assumed that the 290 J.Gołaś (a) State I (b) State II (c) State III Fig. 4. Simply supported inhomogeneous beam of glued laminated wood, composed of seven layers, with stress constants changing over the operation period, according to Table 1. Deflection line of neutral axis of the beam and displacement field σ x (x,z) in: (a) state I, (b) state II, (c) state III – over the beam segment 3.0¬x¬ 3.6m. For the remaining length – according to state II. Note: relation between the horizontal and vertical scales is x : z=1 : 5 Numerical analysis of thick, multi-layered composite... 291 (a) State I (b) State II (c) State III Fig. 5. Double-sided encastered beam of glued laminated wood. Description as for Fig.4, but in case (c) state III occurs over the beam segment 0¬x¬ 0.6m 292 J.Gołaś state III occurs locally, i.e. between x = 3.0m and 3.6m in beam (a) and between x=0 and 0.6m in beam (b). The intensity of shading of the beam surfaces in these figures corresponds to the increase of stress. In both cases an increase of deflection and redistribution of the stress field can be noticed. Table 1 Material constants State I State II State III Number Thickness of layer of layer Ex Gxz Ex Gxz Ex Gxz [m] [MPa] [MPa] [MPa] [MPa] [MPa] [MPa] 1 0.096 11500 700 9200 560 920 56 2 0.048 9500 600 7600 520 7600 520 3 0.156 9000 550 7200 440 7200 440 4 0.400 8000 500 6400 400 6400 400 5 0.156 9000 550 7200 440 7200 440 6 0.048 9500 600 7600 520 7600 520 7 0.096 11500 700 9200 560 920 56 Fig. 6. Diagram of reinforced soil retaining wall structure with backfill reaction An engineering application of the program to a statical analysis of a re- inforced earth structure is given below. Theoretical bases and applications of structures of this type are presented among others by Sawicki and Kulczy- kowski (1986), Sawicki and Leśniewska (1993). A retaining wall of reinforced ground, loaded at its upper surface by surcharge qz = 60kN/m and by side Numerical analysis of thick, multi-layered composite... 293 pressure of the backfill is shown in Fig.6. Ownweight is also taken into acco- unt in the calculations. The structural reinforcement in the formof aluminium 0.6×2.8×800cm flat bars is placed unidirectionally in eight layers. The di- stance between the reinforcement in a layer is 0.5m, the same as the distance between the layers. The following properties were assumed in the calculations: – for the ground γg =18 kN/m 3 Eg =1.2 ·105 kPa νg =0.20 Gg =0.5 ·105 kPa – for the reinforcement Er =0.72 ·108 kPa νr =0.34 µr =0.00067 where γα,Eα, να,Gα, µα, denote correspondingly the bulk density, Young’s modulus, Poisson’s ratio, shear modulus, volumetric proportion of the com- ponent α. The technical moduli of elasticity of a single composite layer are assumed to be the following (Sawicki and Kulczykowski, 1986; Sawicki and Leśniewska, 1993) Ex =1.682 ·10 5 kPa Ez =1.215 ·10 5 kPa νzx =0.145 νxz =0.201 Gxz =0.68 ·10 5 kPa Itwas assumed in the considerations that in a block of soil,measuring L×h= 12.00× 5.00m zero displacement state u = 0, w = 0 occurs in the plane z=−h/2. The computations were carried out for two alternatives: Alternative I: Block of soil is homogeneous and without reinforcement, lo- aded by surcharge qz and own weight. Alternative II: Retaining wall is reinforced in the above-described way, lo- aded by surcharge qz, own weight and triangular backfill side load ac- cording to the formula (see Sawicki and Kulczkowski, 1986) qx = 1 2 ξ tan2 ( 45◦− Φ 2 ) γh where Φ = 37◦ is the angle of internal friction of the ground, γ = 18kN/m3, and ξ =0.618 is the coefficient of interaction of the backfill and retaining wall. Thecomputeddisplacement state forbothalternatives is presented inFig.7 and Fig.8. 294 J.Gołaś wmax =0.00477m; umax =0.00128m; σzmax =148.061kPa Fig. 7. Displacement field u(x,z) and distribution of stress component σ z (x,z) for the first loading alternative. Displacements in the diagram are 250 times enlarged with respect to dimensions of the structure wmax =0.00462m; umax =0.00117m; σzmax =145.841kPa Fig. 8. Displacement field u(x,z) and distribution of stress component σ z (x,z) for the second loading alternative Numerical analysis of thick, multi-layered composite... 295 On the grounds of the numerical analyses of the above-mentioned exam- ples, both versions of the developed program can be recommended for deter- mination of the displacement and stress states in thick,multi-layer ortothropic composite girdershaving the h/L ratio less than 1/2.Theaboveprogramscan also be used for computation of miscellaneous multi-layer engineering struc- tures operating in a state of plane strain. References 1. Gołaś J., 1995, On Limits of Application of Kirchhoff’s Hypothesis in the Theory ofViscoelastic FibrousComposite Plates,Engng Trans., 43, 4, 603-626 2. Gołaś J., 1997, On Necessity of Making Allowance for Shear Strain in Cylin- drical Bending of Fibre Composite Viscoelastic Plates,Arch. Civil Engng, 43, 2, 121-147 3. Kleiber M., 1989, Introduction to Finite Element Method (in Polish), PWN, Warszawa-Poznań 4. LoK.H.,ChristensenR.M.,WuE.M., 1977,AHigh-OrderTheoryofPlate Deformation. Part 2: Laminated Plates, J. of Applied Mechanics, 44, 669-676 5. Mau S.T., Tong P., Pian T.H.H., 1972, Finite Element Solution for Lami- nated Thick Plates, J. Comp. Materials, 6, 304-311 6. Pagano N.J., 1969, Exact Solutions for Composite Laminates in Cylindrical Bending, J. Comp. Materials, 3, 398-411 7. Pian T.H.H., Chen D., 1982, Alternative Ways for Formulation of Hybrid Stress Elements, Intern. J. for Numer. Meth. in Engng, 18, 1679-1684 8. PianT.H.H.,ChenD., 1983,On theSuppressionofZeroEnergyDeformation Modes, Intern. J. for Numer. Meth. in Engng, 19, 1741-1752 9. RikardsR.B., 1988,The Finite ElementsMethod in Plates and Shells Theory (in Russian), Zinatne, Riga 10. SawickiA.,KulczykowskiM., 1986,ElasticAnalysis ofReinforced –Earth Structures,Arch. Hydrotechniki, 33, 3, 299-311 11. Sawicki A., Leśniewska D., 1993,Reinforced Soil. Theory and Applications (in Polish), PWN,Warszawa 12. Spilker R.L, Chou S.C., Orringer O., 1977, Alternate Hybrid-Stress Ele- ments for Analysis of Multilayer Composite Plates, J. Comp. Materials, 11, 51-70 13. Spilker R.L., 1980, A Hybrid-Stress Finite Element Formulation for Thick Multilayer Laminates,Computer and Structures, 11, 507-514 296 J.Gołaś Analiza numeryczna grubych wielowarstwowych dźwigarów kompozytowych z wykorzystaniem hybrydowo-naprężeniowych elementów skończonych Streszczenie W pracy przedstawiono analizę rozwiązań numerycznych dla grubych wielowar- stwowychdźwigarówkompozytowychobciążonych statycznie.Walgorytmachrozwią- zań wykorzystano hybrydowy model naprężeniowy metody elementów skończonych bazujący na zmodyfikowanym funkcjonale wariacyjnymReissnera. Opracowano dwie wersje programu numerycznego dla elementu skończonego o N warstwach. Progra- my przetestowano na przykładzie trójwarstwowej belki z laminatu kompozytowego (z epoksydu wzmocnionego włóknami węglowymi) oraz na przykładzie rozwiązania numerycznego dla tarczy wspornikowej. Analizowano ponadto konkretne przykłady liczbowe dotyczące belek z drewna klejonego warstwowo oraz konstrukcji muru opo- rowego z gruntu zbrojonego. Manuscript received February 24, 2000; accepted for print May 22, 2000