N O N L I N E A R FINITE E L E M E N T ANALYSIS OF MULTILAYERED BEAM-COLUMNS by Rajesh C h a n d r a B . Tech. (Hons.), Indian Institute of Technology, Kharagpur, India A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E O F M A S T E R OF APPLIED SCIENCE in T H E FACULTY OF G R A D U A T E STUDIES CIVIL ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA October,1989 (c) Rajesh Chandra, 1989 In presenting this thesis i n partial fulfilment of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the L i b r a r y shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of C i v i l Engineering The University of B r i t i s h C o l u m b i a Vancouver, B . C . , C A N A D A V 6 T 1W5 October,1989 Abstract A finite element model has been developed i n this thesis for predicting strength and stiffness behavior of multilayered beam-columns. T h e analysis incorporates material and geometric nonlinearities i n order to determine the ultimate load carrying capacity. T h e finite element model takes into account the continuous variability of material properties along the length of layers so that multilayered wood beam-columns can be analyzed. Transverse as well as lateral bending in combination with axial tension or compression can be considered along with different layer configurations, various support and loading conditions. A computer program has been developed based on this formulation. C u b i c beam elements have been used. Numerical integration of the v i r t u a l work equations has been carried out using Gauss quadrature. T h e resulting set of nonlinear equations is solved by using the Newton-Raphson scheme. Numerical investigations have been carried out to verify the results and test the capabilities of the program. Table of Contents Abstract ii Tables v Figures vii Acknowledgement 1 2 3 viii Introduction and Literature Review 1 1.1 General Remarks 1 1.2 Objectives 5 1.3 Thesis Organisation 6 Formulation of T h e o r y 7 2.1 General Assumptions 7 2.2 K i n e m a t i c Relationships 7 2.3 Stress-Strain Relationships 10 2.4 V i r t u a l Work Equation 18 F i n i t e Element Formulation 21 3.1 F i n i t e Element Discretization 22 3.2 P r o b l e m Formulation 23 3.3 Interpolation Functions 25 3.4 V i r t u a l Work Equations 27 3.5 M e t h o d of Computation 30 iii 4 3.6 N u m e r i c a l Integration 32 3.7 Convergence Criterion 33 3.8 F E N M B C — t h e Computer Program 34 Verification and N u m e r i c a l Results 4.1 F i x e d Ended B e a m Under Uniform L o a d 4.1.1 4.2 35 • F i n i t e Element Results 37 Other Test Cases 4.2.1 40 Simply Supported Multilayered B e a m under B i a x i a l Loading for Data-set I 4.2.2 41 S i m p l y Supported Multilayered B e a m under B i a x i a l Loading for Data-set II 4.2.3 45 Eccentric biaxial loading on Simply Supported B e a m for Dataset I 4.2.4 48 Simply Supported B e a m - C o l u m n w i t h axial loading having eccentricity in y-direction 4.2.5 for Data-set I 52 Simply Supported Multilayered B e a m w i t h a Nonlinear C o m pression Zone 5 50 Simply Supported B e a m - C o l u m n w i t h axial loading having eccentricity in y-direction for Data-set II 4.2.6 35 54 Summary, Conclusions and Scope for Future Research 57 5.1 Summary and Conclusions 57 5.2 Scope for Future Research 58 Bibliography 59 iv Tables Comparison of Timoshenko and F E N M B C solution for large deflection of beams v Figures 1.1 B i a x i a l loading on beam-column 2 2.1 Coordinate system for the beam-column 8 2.2 A x i a l and lateral displacement relationship 9 2.3 Large deformation relationships 10 2.4 Bilinear elasto-plastic stress-strain relationship 11 2.5 Various stress-strain relationships 12 2.6 Bilinear stress-strain relationship proposed by B a z a n (1980) 14 2.7 Stress-strain relationship used in the present study 16 2.8 Bilinear stress-strain relationships for various m as used by K o k a . . . 17 2.9 E x p o n e n t i a l compressive stress-strain relationship used for F E N M B C 20 3.1 F i n i t e element discretization of the beam-column 23 3.2 L o c a l coordinate system for the element and the window 24 3.3 Flowchart for iterative technique using Newton-Raphson M e t h o d . . . 31 3.4 Convergence criterion 34 4.1 D a t a for Timoshenko's fixed ended beam bending problem 36 4.2 Comparison of Timoshenko's and F E N M B C ' s large beam-bending 4.3 D a t a for the test cases 4.4 Displacement results for biaxially loaded multilayered beams for data- . . 40 set I 4.5 39 43 Stresses for biaxially loaded multilayered beams for data-set I vi . . . . 44 4.6 Displacement results for biaxially loaded multilayered beams for dataset II 46 4.7 Stresses for biaxially loaded multilayered beams for data-set II . . . . 4.8 Displacement results for simply supported multilayered beams under eccentric biaxial loading for data-set I 4.9 47 49 Displacement results for simply supported multilayered beam-column w i t h axial loading having eccentricity i n y- direction for data-set I . . . 51 4.10 Displacement results for simply supported multilayered beam-column with axial loading having eccentricity i n y-direction for data-set II . . . 53 4.11 Displacement results for simply supported multilayered beams w i t h a nonlinear compression zone 55 4.12 Stress-strain curve for simply supported multilayered beam w i t h a nonlinear compression zone 56 vii Acknowledgement M y special thanks to m y supervisor D r . Ricardo 0 . Foschi, for his guidance and encouragement throughout the course of my research work and i n the preparation of this thesis. I would also like to thank Felix Z . Yao, B r y a n Folz, James D . Dolan and all my friends of the C i v i l Engineering Department, U B C , for their support and valuable suggestions. Financial support in the form of a Research Assistantship from the N a t u r a l Sciences and Engineering Research Council of C a n a d a is gratefully acknowledged. Vlll To M y Parents K a m l a and A k h i l e s h C h a n d r a CHAPTER 1 Introduction and Literature Review 1.1 General Remarks T h e use of glued-laminated multilayered beams and beam-columns is quite common in timber buildings and could be used increasingly i n larger commercial and industrial structures. Analysis and design methods for biaxial bending and axial loading are simplistic and do not represent the true behavior of the composite beam. W i t h the advent of computers, however, it is feasible to incorporate ultimate strength analysis and design. A l s o , reliability methods are being used in revising the timber codes in Canada and elsewhere. T h e present study attempts to provide a more rigorous analysis tool for beam-columns under biaxial bending and axial loading. It has been customary to treat a three-dimensional structure as a collection of two-dimensional, planar structures. A l t h o u g h this k i n d of idealization has resulted in satisfactory designs, it may not allow for an o p t i m u m design. T h e m a i n reason for this simplification has been the complexities involved in three dimensional analysis and the difficulties involved in their hand calculations. However, w i t h computers being widely used in structural engineering, such analyses have become easier and more accurate. Beam-columns are often subjected to biaxial bending coupled w i t h tension or compression. T h i s may result from the space action of the entire framing system or by an axial load eccentrically located with respect to the principal axes of a beam-column cross section, as shown in Figs. 1.1(a) and 1.1(b). A l t h o u g h the two loading conditions 1 Chapter 1: Introduction and Literature Review 2 Figure 1.1: B i a x i a l loading on beam-column shown above are statically equivalent and are considered identical, their nonlinear behavior may be quite different depending on the history of loading. However, in the present analysis, such differences i n loading history have been ignored largely because of the unavailability of data on such stress reversals and failure modes. Recent studies by Buchanan [5] and others have tried to address this issue. In the present analysis monotonically increasing loads have been considered and load reversals have not been allowed. It is also of interest to note that the failure in the beam-column may be caused following one of the several paths shown in F i g . 1.1(c). For example the bending moment in y-direction M y may be applied first up to a certain value following path /, then the bending moment i n z-direction M 2 may be applied upto a certain value following path k, and then finally the axial loading P may be increased following path i until the beam-column fails. The failure may also be caused following paths ghi, bed or any other path as shown i n F i g . 1.1(c). T h e two loading conditions Chapter 1: Introduction and Literature Review 3 as shown in Figs. 1.1(a) and 1.1(b) are different and the response of the beam-column w i l l vary depending on the load history. However, if P, M , M increase proportionally x y as radial loading, 4he two loading conditions are equivalent. In the present study, this proportional increase in loading has been used. W o o d being an organic building material, the modulus of elasticity, tensile and compressive strengths have large variability, both w i t h i n and across members. This variability i n material properties of wood makes an accurate design procedure difficult and the performance of such wood structures more uncertain. In recent years, this uncertainty has been addressed by using reliability based design procedures i n wood structures. T i m b e r codes are being revised to incorporate this change in design philosophy. T h e present study was undertaken to develop a general analysis tool for multilayered beams and beam-columns, w i t h a specific application to glued-laminated beam-columns. K o k a [17] has developed such a model for solid wood beam-columns. However it has limitations of being a planar model, i . e. bending in only one direction, lacking the capability of torsional analysis and can only handle a solid (i. e. single layered) beam-column. T h e model F E N M B C ( F i n i t e Element Nonlinear analysis of Multilayered Beam-Columns) developed i n the present study was to overcome these limitations. It is based on a 3-dimensional analysis including axial deformation, biaxial bending and torsion for multilayered beam columns. B o t h geometric and material nonlinearities were incorporated in order to be able to predict ultimate load carrying capacities. Foschi and Barrett [9] developed a model for predicting the strength and stiffness of glued-laminated beams. In this model, each lamination of thickness t is divided into cells of depth t and width W, and 6 inches (152mm) i n length. E a c h cell is randomly assigned a density and a knot diameter from knot frequency data for the specific lamination in consideration. A modulus of elasticity ( M O E ) and a tensile Chapter 1: Introduction and Literature Review 4 strength value is assigned to each cell. T h e model uses a linear finite element method to obtain stresses i n each cell, considering each one a five noded plate element. T h e model uses weakest link failure criterion to calculate the strength of the beam. T h e Foschi/Barrett model formed the basis for the development of the Karlsruhe model [7]. T h i s finite element model takes into account the material nonlinearity of wood i n compression, by using a linear elasto-plastic stress-strain relationship. It then considers the successive failure of cells in the tension zone. For the past few years research work has also been done in the area of simulation of material properties of lumber, particularly within board variations in modulus of elasticity and strength. ( Bender et a l . [3], Foschi [10], Taylor and Bender [28] and Colling [7] ). T h e simulation of material properties of built-up beam-columns is not the objective of the present study. T h e analysis assumes that these material properties of all laminations are given, or can be simulated along their length. T h e stress-strain relationship used herein are discussed i n detail i n Section 2.3 in Chapter 2, however, a brief survey of axial tension strength of timber has been presented here. It was not until the mid-seventies that the axial tensile strength of timber was given much attention. T h i s was due to the lack of suitable connection details which prevented very high stresses from being developed i n tension members. W i t h the availability of more effective connections, there was a renewed interest i n tensile strength of timber. T h e phenomenon of size effect was observed by K u n e s h and Johnson [18] (1974). They tested commercial sizes of clear Douglas-fir and H e m fir and noticed a significant decrease in tensile strength with increasing cross-sectional dimensions. Norris [23] (1955), Dawe [8] (1964), Nemeth [22] (1965) and others have contributed to the development of appropriate theory dealing w i t h defects i n timber. T h e effect of defects in timber was also ignored in the belief that modulus of rupture was a conservative estimate of tension strength. There was a renewed interest in Chapter 1: Introduction and Literature Review 5 this area after in-grade testing showed failures occuring at stresses much lower than the modulus of rupture. Since then there has been considerable research work done in this area. A general finding from these studies was that the large knots reduced the strength more than smaller knots, and edge knots more than center knots, a result confirmed by Kunesh and Johnson [19] (1972) and Johnson and K u n e s h [16] (1975). A t t e m p t s have been made to combine various characteristics such as knot size, flexural stiffness and slope of grain to predict tensile strength. 1.2 Objectives This study aimed at the following objectives: 1. Development of a finite element analysis model for predicting strength and stiffness behavior of multilayered beam-columns. T h i s analysis would include nonlinear material and geometric behavior of the beam-column so that the ultimate load carrying capacity can be determined. T h e finite element model would take into account the continuous variability of material properties along the length. T h e loading would be axial, transverse and lateral as i n a three dimensional beam-column. 2. Implementation of this analysis into a computer program which would allow flexibility in placing the layers i n different configurations, various support condi- tions and loadings. It would also be user friendly and would require a m i n i m u m of data input. 3. Validation of this computer program by comparing its predictions w i t h available analytical or experimental results. Chapter 1: Introduction and Literature Review 1.3 6 Thesis Organisation Chapter 2 describes the formulation of the theory used in the development of the multilayered beam-column model. T h e general assumptions, kinematic and constitutive relationships are given i n detail. Chapter 3 provides the outline and detail of the finite element formulation of the problem and the computer implementation. T h e v i r t u a l work equations have been set up and the Newton-Raphson iteration technique has been employed to solve the nonlinear equations. Gauss quadrature scheme has been employed for the numerical integration. T h e details of convergence criterion have been outlined. Chapter 4 compares the predictions of this program with available results wherever possible and presents cases including stability and strength analysis. Chapter 5 includes conclusion of the present study and recommendations for further research work in this area. CHAPTER 2 Formulation of Theory 2.1 General Assumptions T h e assumptions made in the analysis of multilayered beam-columns are as follows: 1. Plane section remains plane (with certain qualifications as discussed later in section 2.2). 2. Solid rectangular sections are assumed for each layer and thus for the entire cross-section. 3. A l l layers in the beam are in constant contact with each other, thus no discontinuity exists between layers. 4. Materia] properties may have continuous variation along the span of a layer, however, material properties are assumed to be constant w i t h i n a small division of an element (called a window as explained later in Section 3.1). 5. T h e stress-strain relationship has been assumed to be linear elastic w i t h brittle failure i n tension, and nonlinear having a falling branch w i t h plastic failure in compression as shown in F i g . 2.7. 2.2 Kinematic Relationships T h e orientation of the beam-column with respect to a Cartesian coordinate system is shown i n F i g . 2.1. T h e x-axis coincides with the longitudinal axis of the beam, the 7 Chapter 2: Formulation of Theory 8 y-axis is along the lateral direction of the beam and the z-axis is along transverse direction to the beam. Let u,v and w be the displacements of the geometric center 0 along the x,y and z axes respectively. Let 6 be the rotation of the cross-section about a:-axis. Using the right-hand rule, the rotation 6 is positive as shown in the F i g . 2.1. Figure 2.1: Coordinate system for the beam-column T h e displacements of a generic point A on the cross-section, located at (x,y, z) are related to those of O and approximated as follows: u A WA - = dv u-y-dx 1 dw d6 dx ^ ^ dx (2.1) v-zO (2.2) w -f y6 (2.3) In F i g . 2.2, only the transverse deformation has been shown. It can be readily deduced from this figure that the axial deformation of point A w i l l be reduced by Chapter 2: Formulation of Theory 9 Figure 2.2: A x i a l and lateral displacement relationship an amount z$j£ i f the cross-section rotation is assumed to be the slope (shear deformation neglected). A similar argument can be employed for the second t e r m i n E q . ( 2.1) by replacing w and z w i t h v and y respectively. T h e last t e r m i n Eq.( 2.1) has been included to account for warping of the cross-section due to the rotation 9. T h e nonlinear strain-displacement relationship used i n the present study is now presented. F r o m the geometry of the F i g . 2.3, after deformation, the change i n contour length ds can be expressed as j 2 ^ J 2u, . 9UA. ds Kdx [(1 + -^-) + 2 ( dv dw (^) A 2 A (2.4) 2 E x p a n d i n g it binomially and neglecting the higher order terms, ds can be written as j d s = , h d x [ l , du A + _ _ l,dv A2 + _ ( _ _ ) + 1 - dw A2 ( — ) ] (2.5) T h u s the axial strain e can be expressed as e = ds — dx dx 1 ,dw As du 1. dv *~> ^ dx ' dx ^ 2^ dx A A { 2 2 (2.6) Chapter 2: Formulation of Theory 10 Figure 2.3: Large deformation relationships Substituting u v Al A and w A from Eqs.( 2.1), (2.2) and (2.3), the strain displacement for a general point i n the beam-column is given by du dx 1 dv 2.3 dv dx 2 d6 dx 2 2 2 2 d6 , x dw dx 2 1 dw d6 ^ , Stress-Strain Relationships W h i l e wood is assumed to behave elastically i n tension as discussed earlier i n Section 1.1, many uniaxial stress-strain relationships have been proposed by various researchers to model the nonlinear behaviour of wood i n compression. These are mostly empirical in nature, attempting to theoretically represent the stress-strain relationships obtained from laboratory testing of different types of wood. In this section, the important stress-strain relationships in compression applicable to wood structures are Chapter 2: Formulation of Theory 11 reviewed first and then the one used in the present analysis is discussed. A simple elastic-perfectly plastic bilinear approximation as shown i n F i g . 2.4 has been used quite often to calculate the ultimate load carrying capacity of wood members. Compression Tension (a) stress-strain relationship (b) stress distribution Figure 2.4: Bilinear elasto-plastic stress-strain relationship Y l i n e n [29] (1956) proposed a simple bilinear relationship, which was used by M a l h o t r a and M a z u r [20] (1971) in their research work. T h i s relationship can be given as e=±[co--(l-c)f ln(\-jJ] (2.8) c where e is the normal strain, o is the normal stress, f is m a x i m u m compressive stress, c E is the initial modulus of elasticity and c is a parameter depending on the shape of the curve. A s can be readily found out, the curve described by Eq.( 2.8) and shown in F i g . 2.5a is tangent to the elastic modulus at the origin and tangent to the ultimate compressive stress for large strain. Chapter 2: Formulation of Theory Ylinen(1956) 12 O'HalloranCWS) Qlos(1978) fc tan0 E=tana (a) (b) (c) Figure 2.5: Various stress-strain relationships M a l h o t r a and M a z u r [20] have used this relationship extensively to predict the buckling strength of solid timber columns and described the curve as a good approximation to the results of laboratory testing of 144 clear eastern spruce wood having different moisture contents. However they have not considered the shape of the curve beyond the ultimate load. G o o d m a n and B o d i g [15] (1971) carried out extensive tests on clear dry wood i n compression parallel to the grain and reported the test results. O ' H a l l o r a n [24] (1973) has used these data to propose a mathematical equation for the stress-strain curve at various grain angles and grain orientations. T h i s equation has been given i n the form cr = Ee- Ae n (2.9) where a is stress, t is strain, E modulus of elasticity, A and n are equation constants determined by fitting the equation to a given set of experimental data. If the strain at peak stress is found to a certain ratio, r, of equivalent strain under elastic conditions, Chapter 2: Formulation of Theory 13 the parameter A and n can be found from n = — —r — 1 r a where f c = m 4 ^ (2.10) < -»> 2 is the m a x i m u m stress. A typical plot of the fitted curve is shown in F i g . 2.5(b). T h e equation cannot be used beyond the m a x i m u m compressive stress / , because it drops rapidly to negative c stress values. O'Halloran claims that this is not a serious problem failing to recognize that the shape of the falling branch of the curve is needed to predict the ultimate bending strength. It is true that the stress-strain relationship beyond m a x i m u m load cannot be qualified easily in an axial compression test, because it is largely a function of the test machine characteristics and the rate of loading, but a description of stressstrain behaviour beyond ultimate load is essential to the development of an ultimate bending strength theory. A simple bilinear proposal by Bazan [2] (1980) has been illustrated i n F i g . 2.6. Bazan assumed without any supporting argument that the slope of the falling branch is a variable which can be arbitrarily taken as that value which produces maxim u m bending moment for any neutral axis depth. A different assumption used by Buchanan [5] (1984) is that the slope of the falling branch is a material property which can be estimated as part of the calibration of the computer model to test results. A comprehensive study of the stress-strain relationship of timber with defects, in compression parallel to the grain, has been made by Glos [14] (1978). O n the basis of extensive experimental testing, Glos proposed a curve of the shape shown i n F i g . 2.5 which is similar to that reported for clear wood by Bechtel and Norris [4] (1952), Moe [21] (1961) and others. T h e curve is characterized by four parameters as shown in F i g . 2.7. T h e equation of the curve is given by f= a (2.12) Chapter 2: Formulation of Theory 14 stress / tension strain compression / (a) stress-strain relationship (b) stress distribution Figure 2.6: Bilinear stress-strain relationship proposed by B a z a n (1980) where Is 6/5(1-£) G 2 = G = 4 E Is where / is the stress, e is strain, E is modulus of elasticity, f is m a x i m u m compressive c stress, f s is the asymptotic compression stress for large strain and e is the strain at x m a x i m u m stress. Glos has estimated the four parameters to define the shape of the curve from four measurable wood properties: density, moisture content, knot ratio and percentage Chapter 2: Formulation of Theory 15 compression wood. M u l t i p l e curvilinear regression techniques have been used to express expected values of each parameter in terms of the four properties, using length regression equations. K o k a [17] (1987) has used the bilinear stress-strain relationship as proposed by B a z a n [2] and later modified by Buchanan [5] (1984), who pointed out that the slope of the falling branch of the stress-strain relationship is considered to be a material property. K o k a has used the following relationship o = E e - [E e + | / | ( l + m) + m £ e ] ( l - A ( e + |e |)) 0 0 c 0 (2.13) c where a is stress, EQ is elasticity modulus, e is strain, f c is m a x i m u m compressive stress, m is the slope of the falling branch on compression side, e is the m a x i m u m c strain corresponding to the m a x i m u m compressive stress f c and A ( e + |e |) is the step c function defined as follows e > - | e | =• A ( e + | e | ) = 1 (2.14) e < -\e \ (2.15) c e A ( e + |e |) = 0 e e T h e bilinear relationship for different values of m has been shown i n F i g . 2.8. A further extension of this relationship was made to obtain a more realistic approximation of the stress-strain behaviour of wood i n compression parallel to the grain by employing an exponential variation on the compression side. O n the tension side, the stress-strain varies linearly and elastically up to a m a x i m u m tensile stress F , upon which a brittle failure occurs. t a = Ee - [Ee - F{e)]{l - A{e)) (2.16) [-\ \-m Ee](l-eft\) (2.17) T h e function F(e) for e < 0 has the form F(e) = Po 1 Chapter 2: Formulation of Theory 16 Brittle failure l«cl Plastic failure IPfJ Figure 2.7: Stress-strain relationship used i n the present study In Eqs.( 2.16) and (2.17) above, o is the stress, E is the modulus of elasticity, e is the strain, |po| is the intercept on the stress axis resulting from the extrapolation of the tangent to the falling branch of the stress-strain relationship as shown i n F i g . 2.7 and mi is the slope of the falling branch. In Eq.( 2.16), A ( e ) is a step function which can be defined as follows A(£) 1, if£>0 0, if £ < 0 In Eq.( 2.16) only the three parameters E, f c (2.18) and mj need be specified. T h e re- maining parameters of the equation can be calculated in the following manner. From Eq.( 2.16), for e = e the following expression can be readily obtained c [-\p \-m Ee )(l-e%\) 0 1 c = -\FJ (2.19) Differentiating Eq.( 2.16) with respect to e and substituting e — £ , the following c Chapter 2: Formulation of Theory 17 cj stress / tension t strain compression m>0 ... / / . - m-0 m<0 / , m=-1 Figure 2.8: Bilinear stress-strain relationships for various m as used by K o k a equation can be obtained - m E(\ x - e-l^ k c | m E\e \)^-e-^^ = 0 bo I ) + (|po| - x c (2.20) F r o m Eq.( 2.19), bo| - F c m\E\e (2.21) c Defininig x =e (2.22) and taking the logarithm of both sides, E Inx = -boT l l (2.23) | £ c l Substituting x into Eqs.( 2.20) and (2.21), and rearranging terms, —miE(\ F E — x) + — — — x 1 - x \p \ c - 0 (2.24) 0 Since from Eq.( 2.24), | p | 0 c a n he written as bol = \F C X m x ( l — x) 2 (2.25) Chapter 2: Formulation of Theory 18 finally, from Eqs.( 2.20) and (2.25), x(l + m i ) + m-[xlnx — mj (2.26) For a given mi, x can be computed and from Eq.( 2.23) the m a x i m u m strain |e | c corresponding to the m a x i m u m stress F can be written as c | | = _%i/nz (2.27) hi £ c A simple iterative scheme employing Newton's method can be used to determine the unknown from Eq.( 2.26). Defining the function ip as V> = x(\ + mi) + mixlnx — mi (2.28) following the usual procedure in Newton's method, the following iterative equation can be obtained _ x , ( l + mi) + mixjnx, - mi 1 + 2 m i + milnxi A typical stress-strain plot used in F E N M B C has been shown i n F i g . 2.9, corresponding to E = 1 4 , 0 0 0 M P a , m ! = 0.25, F = c UQMPa. T h e stress-strain relationship used in this study is a better representation of experimental results, as it can have a continuous curvilinear and smooth peak as opposed to the abrupt change of curvature i n the simple bilinear stress-strain relationship.The falling branch can have different slopes and the failure in compression can be easily programmed. Also, the drawback of falling rapidly to negative values of the exponential stress-strain relation, as given by O ' H a l l o r a n [24], has been taken care of by the variable slope m i . 2.4 Virtual Work Equation T h e Principle of V i r t u a l Work has been used to determine the governing equilibrium equations of the beam-column. If we apply a virtual displacement 6a to the system, Chapter 2: Formulation of Theory 19 compatible w i t h boundary conditions, the resulting internal work (6U) and external work (8W) done by the stresses and the external forces respectively are given as SU = 8W = j 8e odV (2.30) J 8a FdS (2.31) T v T s where the internal work is integrated over the entire volume V and the external work is integrated over the surface area S. F is the external force acting per unit surface area. T h e principle of v i r t u a l work states that internal and external work must be equal: J 8e adV T = J v 8a FdV T (2.32) 20 Chapter 2: Formulation of Theory CO D_ CO <D CO to Q> k— "to CD > "(0 CO CD Q. E o o (100) (0.025) (0.02) (0.015) (0.01) (0.005) Strains (m/m) Figure 2.9: Exponential compressive stress-strain relationship used for F E N M B C CHAPTER 3 Finite Element Formulation T h e basic concept underlying the finite element method is that a structure can be modelled analytically by subdividing it into a finite number of elements. Within each finite element, a set of functions are assumed to approximate the stresses or displacements in that region. T h e set of approximating functions contain unknown parameters and are chosen to ensure continuity throughout the system. A p p l i c a t i o n of the principle of v i r t u a l work results in a system of algebraic equations for the parameters i n the approximating functions. A l t h o u g h finite element formulation can be based either on stress fields or displacement fields, most often a displacement based finite element formulation is applied. Such an approach is followed here. Correspondingly, the treatment for the complete structure is then accomplished by studying the stiffness matrices for the elements. T h e steps involved in a linear finite element analysis are 1. Discretization of the body into finite elements. 2. Evaluation of element stiffness matrices by deriving nodal force-displacement relationship. 3. Assemblage of the stiffness and force matrices for the system of elements and nodes. 4. Introduction of boundary conditions. 5. Solution of resulting equations for nodal displacements. 21 Chapter 3: Finite Element 22 Formulation 6. Calculation of strains and stresses based on nodal displacements. In case of nonlinear finite element analysis, an incremental/iterative technique is required. T h e stiffness matrix calculated above is either updated for each load increment or the same stiffness matrix is used until the solution displacement vector converges. These procedures will be discussed i n detail later i n this chapter. 3.1 Finite Element Discretization T h e beam-column arrangement and the finite element discretization is shown i n F i g . 3.1. T h e span is subdivided into elements of length A . T h e layers for each element are further subdivided into windows of length A t o . T h e numbering of layers has been done from top to bottom. T h e elements are numbered from left to right and so are the windows inside every element. T h e nodes are also numbered from left to right. T h e material properties within each window are assumed to be constant and thus, a finer division would allow closer representation of a continuous variation of material properties along the span. T h e same could be achieved by employing a larger number of elements without windows, but then a larger stiffness m a t r i x would be obtained. Thus while keeping the size of the stiffness matrix w i t h i n reasonable limits, the concept of windows permits consideration of the continuous variation of the material properties along the span. T h i s results i n increased efficiency of the program while solving for the displacements and stresses. Also, as this program is to be used for the ultimate load analysis involving material and geometry nonlinearities, the smaller the stiffness matrix, the greater the savings of computer time in the overall solution. Chapter 3: Finite Element Formulation 23 a typical window 1 2 3 4 / y.v ,;*w, I— — «SC N| *' • • : / x,u / r a typical node 'z,w A layers Figure 3.1: F i n i t e element discretization of the beam-column 3.2 Problem Formulation A beam element of length A with two end nodes, as shown in F i g . 3.2, is used in this formulation. A local coordinate <f (— 1 <<f < l ) i s used for each element along the x-axis. T h e x-coordinate of any point inside the element can be expressed as A, where x c (3.1) is the x-coordinate of the center of the element length. T h e element has been divided into a number of windows (NW) Aw - A NW of equal length, hence (3.2) A second local coordinate 7;( — 1 < 77 < 1) is used within a window. A transformation function relating n inside the window and £ is required (refer to F i g . 3.2), since numerical integration has to be performed inside every window and summed over the number of windows and the number of layers to obtain the element stiffness matrix. 24 Chapter 3: Finite Element Formulation 4 node 1 -+i <> node 2 4W •+1 - 1 r Figure 3.2: L o c a l coordinate system for the element and the window Thus, for the i th window, F i g . 3.2, A A x - — + (i - l ) A t u = z + — 6 (3.3) %2 = %c ~ (3.4) XQ A x - — + IAW xi = = c c + *Aio = a: + y f c Aw — c 2 A, x + —£o = (3.5) c From the above we can calculate £i,Cj2 and £o as = |(.t-l)A«;-l — 1 —iAw A 6 2 A A — tAw A W A 1 1 Also fc - 6 Au> 2 A Hence we can write £ i n terms of 77 as 2Au; A Atw Aw -1) + ^ (3.6) Chapter 3: Finite Element Formulation 25 T h e displacement vector consists of nodal degrees of freedom u,v, w, 8 and their first derivatives u',v',w' and 9' respectively. Thus there are eight degrees of freedom (dof) per node giving a displacement vector of sixteen dof for a beam element. The displacement vector is arranged according to the following form: a T = {u u[,v v[,w w[,9 J[,u ,u ,v ,v ,w ,w' ,9 ,9 } u u u 1 2 2 2 2 2 2 2 (3.7) 2 where the subscripts 1 and 2 refer to node 1 and node 2 respectively. 3.3 Interpolation Functions To satisfy the completeness criterion, the assumed displacement field must be at least a complete polynomial of order equal to the highest derivative occurring i n the straindisplacement relationship. Also the displacement function must be continuous w i t h i n the element and the displacement must be compatible between adjacent elements. T h e compatibility condition is ensured if the displacement field is continuous on the inter-element boundary upto the derivative of one order less than the highest derivative appearing in the strain-displacement relation. For the beam-column element the strain-displacement relationship contains second derivative in lateral displacements and twist and the first derivative i n axial displacement. Hence it is necessary to choose the displacement function such that u,u',v,v' ,w,w' ,9 and 9' are continuous at nodes. T h i s can be achieved by adopting linear displacement field for u and cubic displacements for other degrees of freedom. However, complete cubic interpolations are used to approximate u,v,w and 9 within each element. In an earlier work by K o k a [17], which employed this element, it was found that a cubic interpolation function for u gives a much improved approximation of the axial stresses. A complete cubic polynomial requires 4 parameters to define the function.The displacements and their first derivatives at the two nodes provide sufficient parameters to fully describe a cubic polynomial function. Thus we Chapter 3: Finite Element Formulation 26 can write u, u, w and 6 as du N {—)i u = v dv = N^ + ^i—^ dx w = 9 = N + L U I + 2 3 + du N {—)j (3.8) A dv + NzVj + N^ — h dx + N (^) 2 N UJ + N { 3Wj A^,.+ + N*^), + JV 0 + J V ( £ ) 3 (3.9) j 4 (3.10) (3.11) > where subscripts i and j refer to the first and second nodes respectively of each element. Also N i ( = ^3 = (| + f f - ^ N = | ( - i - £ + ^ + f ) 4 ) 3 3 - 1 2 ) (3-14) (3.15) Alternatively we can write these as u = P (0a «' = P i ( 0 a t, = Q (£)a i u; = R (0a w' = R i ( e ) a 6 = S (£)a 0' = S i ( O a T T r r (3.16) T ^ Q ^ a T T t>" = Q T 2 = R 0" = S r 2 (c> T 2 (0a (Oa (3.17) (3.18) (3.19) where single primes denote first derivatives with respect to x and double primes denote second derivatives with respect to x. P , Q , R , S, P i , Q i . . . etc. are vectors as given below P T = {N N 0,0,0,0,0,0,^3,^4,0,0,0,0,0,0} (3.20) Px r = {N[, ^ , 0 , 0 , 0 , 0 , 0 , 0 , ^ , ^ , 0 , 0 , 0 , 0 , 0 , 0 } (3.21) Q T = {0,0,^,^2,0,0,0,0,0,0,^,^4,0,0,0,0} (3.22) u 2l Chapter 3: Finite Element Formulation 27 Qi T = {0,0,^,^,0,0,0,0,0,0,^,^,0,0,0,0} (3.23) Q T = {0,0, ^ " , ^ ' , 0 , 0 , 0 , 0 , 0 , 0 , ^ ' , ^ ' , 0 , 0 , 0 , 0 } (3.24) 2 ... etc. where N , N , N , N x 2 3 4 are as given earlier and N{, N^, N!$, N' are the first derivatives 4 and N", N?, N£, N% are the second derivatives w i t h respect to x. 3.4 Virtual Work Equations F r o m Eq.( 2.7), the strain can now be written i n the following form e [ P i -yQ = r + yzS T 2 - T 2 + ^a [(Qi - *Si)(Qi T = Bo T a+ia B T i zR - 2S r T 2 T X ]a ) + (Ri + ySi)(Ri + ySi )]a T r (3.25) (3.26) a where B 0 B 1 = P i - 2/Q2 + yzS - 2R2 = (Qi - 2 S ) ( Q i - z S (3.27) 2 T 1 T 1 ) - r ( R i + j/S )(R 1 T 1 - j/S r T 1 ) (3.28) B o corresponds to the linear part of e and B i to the nonlinear. A p p l y i n g a virtual displacement 6a, as a —* a + A<5a, the following can be obtained e(A) = B 0 = B 0 T T ( a - f A<5a) + ^ ( a a + AB T 0 + A < 5 a ) B ( a + A6a) T T 1 6a + ^ a B i a + ^(5a Bia T T Zj Li + ^ a B i * a + y<5a Bi<5a r r Thus the v i r t u a l strain can be written as 6e = (^)x-. 0 = Bo^a + ^ B i a = B <5a + i a B i ( 5 a + ^a Bi<$a r 0 T T + i a ^ a T Chapter 3: Finite Element Formulation = (B = (B 28 + ^a Bx + ^a B r r 0 r )6a r 1 a (^±^))8a T (3.29) T 0 + and since B i is symmetrical (refer E q . 3.28) i.e. B i = B i , Eq.( 3.29) can be written r as 8e = ( B T 0 + a B )<5a (3.30) T 1 T h e stress a can be written as ( from E q . 2.16 ) a = £(B T 0 + ^ a B i ) a + G(£) (3.31) T where G(e) = -[Ee - {-\ Po \ - m^e}^ - e W ] ( l - A(e)) £ (3.32) da can be written as da = E{B + a Bi)<5a + H{e)(B T T 0 + a B )8a T T 0 1 (3.33) where dG{e) H(e) de = (-£-m,£(l -M ~ ) ic Co) e - ( - | p o | - m E(e - e ))^e^ IPol (e x £o) 0 (3.34) Hence the v i r t u a l work equation can be written as J 8e adV - J 8a FdS T v T = 0 or f 6a {B Jv + & a)c-dV T 0 x - [ 8a FdS Js T = 0 (3.35) Since this must hold for any 6a, j^(B + 0 Bia)<7<n/ - J FdS s = 0 (3.36) 29 Chapter 3: Finite Element Formulation Thus we can write a function 0(a) as J {B + B a)adV 0(a) - v 0 J FdS - l (3.37) s T h e solution vector a must satisfy 0(a) = 0. To find this solution an iterative procedure is used. T h e Newton-Raphson method is a commonly used iterative technique to solve non-linear equations. It uses a truncated Taylor series expansion of the function 0(a) 0 ( a + A a ) = 0(a) + [ ^ ] A a da (3.38) or A a = 0(a + A a ) - 0 ( a ) [ ^ ^ ] aa (3.39) 1 Setting 0(a + A a ) = 0, Aa = a - a, = 1 + 1 (3.40) or a, = a, - [K r ]-V(a,) + 1 (3.41) Since 60(a) = J BxSaLardV + J (B + = J aB 6a +J [(B + B a)(E 0 v v 1 0 B &)6adV 1 + H(e)}(B + T 1 0 a B )6a T 1 (3.42) Hence the tangent stiffness matrix can be written as [K ] T = J [BiEBo a+^B Ea B a T +B 1 i a £Bo T r 0 a T 0 £ a B i + B H(e)B r i B Ea B T 1 + B +B a//(e)B 1 + B EBo + T v 0 0 T 0 + Bo#(e)a B! + B a r Y ( e ) a B + BxG(e)}dV T 1 J 1 1 r (3.43) Chapter 3: Finite Element Formulation 30 T h i s tangent stiffness matrix is computed for each window and summed over the layers to find the element stiffness matrix. These are then assembled to obtain the global stiffness matrix. E a c h value in the incremental solution vector A a is compared against an acceptable tolerance specified by the user and it is determined whether further iteration is needed to obtain a sufficiently accurate solution vector a. 3.5 Method of Computation T h e Cholesky decomposition method is used in this program to invert the global stiffness matrix [K<sr] and to determine the solution displacement vector a. The symmetrical and the banded character of the tangent stiffness m a t r i x is utilised and only the band i n the lower triange of the matrix is computed. After entering the boundary conditions, the [Kj] matrix is decomposed and the incremental solution displacement vector A a is obtained using the residual load vector 0(a,). A flowchart for the computation is given in F i g . 3.3. Chapter 3: Finite Element Formulation Q Start 31 ^ Read geometry, material properties and load characteristics I Initialise [K^ = 0 ; {R} • 0 {u } = 0; iter = 0 If there is no further load increment r = r+1; increment the load vector f {R} = {R} + {R} r r-1 r Stop iteration n=0 I n = n+1 (— Calculate stiffness matrix IK] Apply boundary condition I Decompose and solve for {Aa} Calculate the incremented displacement vector rice^Check convergence \ no converged^ Solution displacement vector {X} Calculate stresses {o} F i g u r e 3.3: Flowchart for iterative technique using Newton-Raphson M e t h o d Chapter 3: Finite Element Formulation 3.6 32 Numerical Integration Recourse must be made to numerical integration methods to obtain the m a t r i x [ K j ] . Gauss quadrature allows the integral / to be computed as /= f f (£)% = £, Witt (- ) 3 44 where £, are fixed abscissae and Hi corresponding weights for the n chosen integration points. If a polynomial expression is to be integrated, it is easy to see that for n sampling points, there are 2n unknowns (Hi and £,) and hence a p o l y n o m i a l of degree (2n — 1) could be constructed and exactly integrated. T h e error thus is of the order of 0(A ). 2n In the present problem, since there are terms involving x,y and z, a Gaussian integration scheme is employed, which is of the form 1 LLLj^^OdZdndC = nl = n2 n3 E E E W * / ( f n 1 ; . ( 0 fc=ij=ii=i (3-45) where n l , n2, n3 are number of integration points i n x, y and z directions respectively. In the present case the integral I is of the form 1 = J-zJo Jo AwBT 2 2 2 z9 x z dxd /•! AwBT ,z f( >yi ) y L /-i L i nl z9 n2 v (- ) dz 3 46 ' (WWt - (3 47) n3 E E E W J ( ( , ^ a ) k=\j=i i=i (3-48) where Aw is the length, B the width and T the thickness of each window. As was pointed out earlier, n Gauss points integrate a p o l y n o m i a l of the order (2n — 1) exactly. In the stiffness expression in E q . 3.43, there is a polynomial of order 8 in the x—direction, of order 4 each in y— and z—directions. Hence a 5-point Gaussian integration in x-direction and a 3-point Gaussian integration in y- and z-directions are needed respectively. However it is important to note that i n the program F E N M B C Chapter 3: Finite Element Formulation 33 there is a flexibility of choosing the number of Gauss points from 1 to 5. T h i s has been done because the program can handle multilayered beam columns and a number of windows w i t h i n an element, to take care of the continuous variation of material properties along the span. For example, if there is a beam column w i t h 8 layers and 5 windows per element, one Gauss point in both x- and z-directions would give fairly accurate results as shown later. Also depending on the orientation of the loading the number of Gauss points i n y- and z-directions should be suitably selected. For example if there is a distributed load in the z-direction on a beam column with only one layer, three Gauss points in z-direction and one Gauss point in ^-direction would give accurate results. However if the applied load is acting in the y-direction then three Gauss points in y-direction and one Gauss point in z-direction would suffice. In case of biaxial loading, sufficient number of Gauss points should be chosen in both y- and z-directions. Thus one has to be careful in choosing the number of Gauss points i n various directions. F i n a l l y it should be mentioned that since numerical integration has been employed, stresses and strains can be determined only at the Gaussian points. If stresses and strains are needed at other points, a suitable interpolation may be employed. 3.7 Convergence Criterion A n Euclidean norm criterion has been used to check the convergence at every load step. If a is the solution vector at the end of r r the end of (r + l) tk th iteration and a r + 1 is the solution at iteration in F i g . 3.4, then A x = |a +i — a \ represents the difference r r between the two displacement vectors. | a | and | a + i | can be given as follows r Kl = 2 r NDOF E i=i «?(«') NDOF :=1 where a (i) r and a (i) T+i are components of a and a r r + i respectively. Chapter 3: Finite Element Formulation 34 a r+1"a r ^ ^ ^ ^ a r + 1 0 Figure 3.4: Convergence criterion T h e n A a can be written as A a = E i=l y/(a i{i)-a (i)y r+ r (3.49) The convergence criterion for this norm is given as Aa — T T < specified tolerance \a \ (3.50) r 3.8 F E N M B C — t h e Computer Program A computer program called F E N M B C (Finite Element Nonlinear analysis of Multilayered B e a m C o l u m n ) based on the above finite element formulation was developed. It is a. versatile finite element program capable of analyzing a 3-dimensional multilayered beam-columns with numerous possibilities in terms of geometric, material properties and loading conditions. Extensive numerical testing was done and some of the results have been reported in the next chapter. CHAPTER 4 Verification and Numerical Results Numerical results obtained from F E N M B C were compared with analytical, experimental or numerical results published in the literature, whenever possible. T h e cases presented below involve geometric nonlinearity, material nonlinearity, biaxial bending, buckling, lateral stability and torsion. These cases have been tested for both single as well as multilayered beam-columns and the results are reported here. The program F E N M B C has been implemented and tested on a mainframe, A m dahl 5840, a Sun 4/260 workstation and an A S T 286 P r e m i u m ( I B M P C - A T compatible) at the University of B r i t i s h Columbia. Double Precision (Real*8) arithmetic is used throughout the program to reduce the effect of round-off. The program has the capability of handling 20 layers, 20 elements and 10 windows per element. E a c h window can have a m a x i m u m of 5 Gauss points i n the three principal axes directions. T h e dimensions in the program can be easily modified to accomodate a greater number of layers, elements, windows and Gauss points. The large deflection results given by Timoshenko [25] for beam-bending were used to verify the program's numerical results. 4.1 Fixed Ended Beam Under Uniform Load The numerical results for a fixed ended beam were obtained. T h e data for the beam have been given below. Timoshenko used the energy method to obtain an approximate solution for a plate 35 Chapter 4: Verification and Numerical Results 36 MQht>0.125 Figure 4.1: D a t a for Timoshenko's fixed ended beam bending problem w i t h fixed edges. He obtained the total strain energy of the plate by adding the energy due to bending and the energy due to the strain of the middle surface and then used the principle of virtual work. He used the following displacement functions satisfying the boundary conditions imposed by the clamped edges for a rectangular plate with sides 2a and 2b. u = (a ~ x ){b - y )x{b + b y v = (a -x )(b -y )y{co +Co y w = (a -x ) (b -y ) {a 2 2 2 2 2 2 2 2 00 2 2 2 2 2 2 00 0 2 02 + b x 2 2 2 (4.1) + c xy) (4.2) 2 + c x 2 2 2 2 22 a x) 2 20 + b xy) 22 20 + a y+ 02 2 20 (4.3) T h e values of all the parameters were estimated for various intensities of the load q and for different aspect ratios | , assuming v = 0.3. T h e results reported for an infinitely long plate (for | = 0) have been used to compare w i t h the results of F E N M B C . Chapter 4: Verification and Numerical Results 4.1.1 37 Finite Element Results T h e program F E N M B C was run using the data shown above w i t h some additional data required for the program: number of layers nlayer = 1, number of windows nwindo — 3, number of Gauss points along span i n x-direction ngauss = 5, number of Gauss points along the width in y-direction ngausy = 2, N u m b e r of Gauss points along the depth in z-direction ngausz — 5, slope of the compression curve m\ = 0.2, and high values for m a x i m u m allowable tensile and compressive stresses of F = 6.0 x c l0 N/m 15 2 and F = 6.0 x 1 0 i V / m 1 5 t 2 to ensure full elastic behavior. T h e distributed load in the z-direction was given as q = 10.7315./V/m and the number of load steps z was specified as 16 for a smooth load-displacement relationship. T h e results have been tabulated in Table 4.1 and subsequently shown graphically i n F i g . 4.2. The program F E N M B C ' s results are in good agreement with Timoshenko's results. In the Table 4.1, the theoretical linear, Timoshenko's nonlinear and F E N M B C ' s nonlinear results have been presented, span, D — EI, q is the uniformly distributed load, / is the where E is modulus of elasticity and / is moment of inertia of the beam cross-section along y-axis, h is the depth of the beam and w max is the m a x i m u m displacement i n z- direction. F r o m Table 4.1, It should be noted that the order of accuracy is quite good even with 6 elements as compared to 20 elements. Chapter 4: Verification and Numerical Results Table 4.1: Comparison of Timoshenko and F E N M B C solution for large deflection beams Linear theoretical ^ Wrnnr Dh 0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 120.0 130.0 140.0 150.0 160.0 Timoshenko nonlinear 0.0 0.4166667 0.8333333 1.2500000 1.6666667 2.0833333 2.5000000 2.1966667 3.3333333 3.7500000 4.1666667 4.5833333 5.0000000 5.4166667 5.8333333 6.2500000 6.6666667 0.0 0.40 0.62 0.84 0.98 1.10 1.21 1.31 1.39 1.47 1.54 1.62 1.67 1.73 1.79 1.84 1.90 FENMBC nonlinear 6 elements FENMBC FENMBC nonlinear nonlinear 10 elements 20 elements j 0.0 0.4081488 0.6828256 0.8791920 1.0311120 1.1571440 1.2654880 1.3610000 1.4467600 1.5248480 1.5967040 1.6634000 1.7257520 1.7861440 1.8413200 1.8937040 1.9436160 0.0 0.4082008 0.6832096 0.8798320 1.0321920 1.1586880 1.2675280 1.3635520 1.4498240 1.5284240 1.6007920 1.6680080 1.7308720 1.7917760 1.8474560 1.9003360 1.9507440 0.0 0.4082064 0.6832456 0.8799200 1.0323440 1.1589200 1.2678400 1.3639520 1.4503200 1.5290160 1.6014800 1.6688000 1.7317680 1.7927760 1,8485600 1.9015520 1.9520720 Chapter 4: Verification and Numerical Results Figure 4.2: Comparison of Timoshenko's and F E N M B C ' s large beam-bending Chapter 4: Verification and Numerical Results 4.2 40 Other Test Cases Other cases were run and the results have been reported below. These have been compared w i t h the linear beam bending and beam column theories for lack of available nonlinear results. T h e same two data sets were used throughout the investigation unless otherwise specified. T h e first data set contained 10 layers w i t h their elasticity m o d u l i placed symmetrically w i t h respect to the y-axis on the beam cross-section as shown i n F i g . 4.3. T h e second data set contained 10 layers w i t h their elasticity m o d u l i placed asymmetrically with respect to the y-axis of the beam cross-section as shown i n F i g . 4.3. M.O.E.(kN/m*) 0.066+08 0.12e+08 0.066+08 0.106+08 E E 0.086+08 9 © in M.O.E.(kN/m») 0.14e+08 0.086+08 E 0.086+08 ^ 0.066+08 0.106+08 k_ Q> 0.066+08 |2 o 0.086+08 0.106+08 J2 0.106+08 >. 0.126+08 0.126+08 0.146+08 0.126+08 0.146+08 0.146+08 , 150 mm 150 mm (a) Data set i (b) Data set II Span L==6.0 m (unless otherwise specified) Maximum allowable tensile stress f\ =0.4e+08 N/m 2 Maximum allowable comp. stress F =0.4e+08 N/m Q Slope of the comp. curve m ==.25 Figure 4.3: D a t a for the test cases 2 Chapter 4: Verification 4.2.1 and Numerical 41 Results Simply Supported Multilayered Beam under Biaxial Loading for Data-set I T h e case of a multilayered beam with data-set I as shown i n F i g . 4.3 with concentrated loads in y- and z-directions (refer F i g . 4.4) was run. T h e additional data required for running F E N M B C were: nelem = 10, nwindo b,ngausz = 5. = 3,ngauss = 3,ngausy = T h e m a x i m u m tensile and compressive stress values was taken as 4 . 0 x l 0 / V / m for a l l the layers.The concentrated loads were P = l . O x 10 A^andP = l o 2 6 y 2 1.0 x 1 0 N . T h e displacement results have been reported i n F i g . 4.4. T h e total tensile 6 stresses comprising of biaxial bending stresses and axial stresses at the Gauss point located at nelem — 6,nlayer = 10,nwindo = l,ngauss = l,ngausy = 5,ngausz = 5 has been shown i n F i g . 4.5. T h e linear deflections v and w were calculated from the linear beam bending theory using the method of transformed cross-section for comparison with the nonlinear results obtained from F E N M B C . T h e nonlinear displacements v and w are the maxi m u m displacements (at node no. 6) and are the same for a l l the layers. A s both the ends are hinged, the displacements converge to a lower value from the linear displacement values due to stiffening of the beam resulting from the large deformation terms. T h e rotation 0, as obtained in the result shown in F i g . 4.4, needs some explanation. For symmetric loading in biaxial bending of beam finite element model, one would usually expect displacements only in transverse and lateral directions. In the present case of biaxial bending of a simply supported rectangular multilayered beam, the support condition may be thought of as a knife-edge support. T h e loadings in transverse and lateral directions are applied incrementally. T h e load in y-direction, AP , y produces a displacement Av. exists a moment A M . = AP Z As the first load increment AP Z is applied, there x Av, which causes a rotation in the beam cross-section along x-axis. Also the force AP Z produces the displacement Aw in z-direction. A s Chapter 4: Verification and Numerical Results 42 the loads are incremented gradually, the rotation also increses. In the mathematical formulation of the stiffness matrix given by the E q . 3.43, the rotation comes from the coupled terms involving S i and S , these being the derivatives of the shape functions 2 for 6. T h e odd terms i n y and z after integration cancel out due to the symmetry of layers, however, the even terms add up and produce the desired rotational deformation. T h e total tensile stresses are seen to be linear as indeed they should be, since they are in the tensile zone below the m a x i m u m tensile stress FT- Chapter 4: Verification and Numerical Results Figure 4.4: Displacement results for biaxially loaded multilayered beams for data Chapter 4: Verification and Numerical Results Stress (MPa) 500 | 400 h Strain (m/m) Figure 4.5: Stresses for biaxially loaded multilayered beams for data-set I Chapter 4: Verification and Numerical Results 4.2.2 45 Simply Supported Multilayered Beam under Biaxial Loading for Data-set II T h e other multilayered beam with data-set II as shown i n F i g . 4.3 with concentrated loads i n y- and z-directions (refer F i g . 4.6) was run. T h e additional data required for running F E N M B C were: nelem = 10, nwindo — 3,ngauss = 3,ngausy = 5,ngausz — 5. T h e m a x i m u m tensile and compressive stress values was taken as 4.0 x 1 0 / V / m for all the layers.The concentrated loads were P = 1.0 x 10 ./VandP* = l o 2 6 y 1.0 x 10 7V. T h e displacement results have been reported i n F i g . 4.6. T h e total tensile 6 stresses comprising of biaxial bending stresses and axial stresses at the Gauss point located at nelem = 6,nlayer = 10,nwindo — l,ngauss = l,ngausy = 5,ngausz = 5 has been shown i n F i g . 4.7. T h e linear deflections v and w were again calculated from the linear beam bending theory using equivalent cross-section. T h e transverse and lateral deflections resulting from the rotation 6 were ignored in this calculation due to their small magnitude. T h e nonlinear displacements v and w again converge to a. lower value from the linear displacement values as discussed earlier. The rotation 6 i n this case results from the torsional rigidity as explained in Section 4.2.1 as well as due to asymmetric placement of the beam layers in data-set II. The total tensile stresses are seen to be linear again, as they are i n the tensile zone below the m a x i m u m tensile stress Fj. Chapter 4: Figure 4.6: data-set II Verification and Numerical Results Displacement results for biaxially loaded multilayered beams 46 for Chapter 4: Verification and Numerical Results Stresses (MPa) Strains (m/m) Figure 4.7: Stresses for biaxially loaded multilayered beams for data-set II 47 Chapter 4: Verification and Numerical Results 4.2.3 48 Eccentric biaxial loading on Simply Supported Beam for Data-set I T h e multilayered beam with data-set I as shown i n F i g . 4.3 w i t h eccentrically located concentrated loads in y- and z-directions (refer F i g . 4.8) was run. T h e additional data required for running F E N M B C were the same as i n Section 4.2.1 w i t h eccentricities e y = 2 0 m m and e z = 20mm. T h e displacement results have been reported in F i g . 4.8. T h e nonlinear displacements v and w again converge to a lower value from the linear displacement values as expected. are the loads applied. a negative value. A t the beginning the rotation 0 is zero as However as the first load increment is applied, 0 takes on A s can be seen from the loading i n the F i g . 4.8, the moments produced at the centroidal axis of the beam due to the two eccentric loadings are in opposite directions. Depending on the load magnitude, eccentricity and torsional rigidity, a negative rotation occurs on applying the first load increment. A s the load is increased, the rotation also gets affected by the displacements v and w as explained i n Section 4.2.1. T h e interaction of these moments produces the rotational deformations as shown i n F i g . 4.8. T h e stress-strain diagram has not been included here as the tensile and compressive stresses are still in the linear elastic range. A n example of nonlinear stress-strain behavior under biaxial loading has been discussed later in this chapter. Chapter 4: Verification and Numerical Results 49 Figure 4.8: Displacement results for simply supported multilayered beams under eccentric biaxial loading for data-set I Chapter 4: Verification and Numerical Results 4.2.4 50 Simply Supported Beam-Column with axial loading having eccentricity in y-direction for Data-set I T h e program F E N M B C was tested for buckling and the results were compared with Euler buckling loads for an equivalent cross-sectional area and a mean elasticity modulus E. T h e Euler buckling load thus calculated was PE = 173.49&./V. A n axial load was applied w i t h a very small eccentricity in the y-direction, and increased i n steps of 8.6754kN. A s the load increases, the deflections become larger and finally the beam-column fails due to instability as it approaches the Euler buckling load. This demonstrates the program's capability to analyse beam-columns for buckling failures. T h e results have been reported in F i g . 4.9. Chapter 4: Verification and Numerical Results Figure 4.9: Displacement results for simply supported multilayered beam-col with axial loading having eccentricity in y-direction for data-set I Chapter 4: Verification and Numerical Results 4.2.5 52 Simply Supported Beam-Column with axial loading having eccentricity in y-direction for Data-set II T h e program F E N M B C was further tested w i t h data-set II for buckling and the results were compared w i t h Euler buckling loads for an equivalent cross-sectional area and a mean elasticity modulus E. It should be noted, however, that the centroidal axis and the shear center axis do not coincide and the theoretical buckling load can be calculated using the transformed cross-section method. However, ignoring this fact, the Euler buckling load was calculated just to have an idea of the magnitude and was found to be PE = l73A9kN and was applied with a very small eccentricity i n the y-direction. T h e results have been reported in F i g . 4.10. Chapter 4: Verification and Numerical Results 53 2.5 v (nonlinear) —B— 50,000 w (nonlinear) ---A--- 100,000 © (nonlinear) G 150,000 P 200,000 Axial Compressive Load P (N) Figure 4.10: Displacement results for simply supported multilayered beam-column w i t h axial loading having eccentricity in y-direction for data-set II Chapter 4: Verification and Numerical Results 4.2.6 54 Simply Supported Multilayered Beam with a Nonlinear Compression Zone T h e data-set I was used to test the response of the multilayered beam i n the nonlinear compression zone. In the data-set I, the m a x i m u m allowable compressive stress values were divided by 100 to force the compressive stresses to fall i n the nonlinear range. T h e loading and boundary conditions have been shown in F i g . 4.11 along with the deflections v and w. Once again the nonlinear deflections converge to lower values compared to the linear deflections due to the stiffening of the beam. T h i s case was run to test the response of the multilayered beam in the nonlinear range. T h e stress-strain plot has been shown in the F i g . 4.12. Defining failure criteria for the beam-column in tensile as well as compression zone, the program can be used to predict member strengths and stability. Chapter 4: Verification and Numerical Results 55 Figure 4.11: Displacement results for simply supported multilayered beams w i t h a nonlinear compression zone 56 Chapter 4: Verification and Numerical Results (0.025) (0.02) (0.015) (0.01) (0.005) 0 Strains(m/m) Figure 4.12: Stress-strain curve for simply supported multilayered beam with a nonlinear compression zone CHAPTER 5 Summary, Conclusions and Scope for Future Research 5.1 Summary and Conclusions A direct v i r t u a l work formulation of multilayered beam-columns is developed incorporating geometric as well as material nonlinearities. F i n i t e element total and incremental equilibrium equations are derived based on the assumed displacement model. T h i s method has been used to examine the strength and stability of glued-laminated beam-columns. T h e computer program developed has been verified and tested for several beam-column problems. T h e generality of the model makes it a powerful and versatile tool for solving a large variety of problems. T h e most important question for the user of the finite element method is whether the method yields sufficiently accurate results for his purpose. A s conforming ele- ments have been used for displacement functions satisfying completeness criteria, any desired degree of accuracy, with the limitations of finite element being an approximate method, could be obtained with suitably modifying the dimensions i n the program. A s in any other approximate numerical method, engineering judgement should be exercised i n interpreting the results obtained by the finite element method since accuracy of the results depends on the assumptions employed i n the formulation and the limitations of the material idealization. Since the finite element model developed in the present study is quite general 57 Chapter 5: Summary, Conclusions and Scope for Future Research 58 in its application to glued-laminated beam-columns, it is now possible to carry out parametric studies for the codification of design procedures to include instability and effects of plastic members for realistic loading and material response. 5.2 Scope for Future Research T h e present model does not include shear deformations i n its formulation and hence it may not give very accurate results for short and stocky beam-columns. A n inclusion of shear effects would enhance its capabilities. A load path dependent plastic response of the material would make the model more realistic i n terms of predicting long term behavior of the beam-column undergoing load reversals. A n inclusion of creep behavior of wooden members may be required. Bibliography [1] A b r a m o w i t z , M i l t o n and Irene A . Stegun, 1965, Handbook of Mathematical Func- tions, Dover Publications, Inc., New York. [2] Bazan, I . M . M . , 1980, Ultimate Bending Strength of Timber Beams, P h . D . thesis Dissertation, Nova Scotia Technical College, Halifax, N . S . [3] Bender, D . A . , F . E . Woeste, E . L . Schaffer and C M . M a r x , 1985, Reliability for- mulation for the strength and fire endurance of glued-laminated beams, U S D A Forest Service Research Paper F P L 460, U . S . Forest Products Laboratory, M a d i son, W I . [4] Bechtel, S . C . and C . B . Norris, 1952, Strength of Wood Beams and Rectangu- lar Cross-section as affected by Span-Depth Ratio U S D A Forest Service, Forest Product Laboratory Report, M . R 1 9 1 0 . [5] Buchanan, Andrew H . , 1984, Strength Model and Design Methods for Bending and Axial Load Interaction in Timber Members, thesis presented to the Unversity of B r i t i s h C o l u m b i a , at Vancouver, B . C . , in partial fulfilment of the requirements for the degree of Doctor of Philosophy. [6] Chen, W . F . and Atsuta, T . , 1976, Theory of Beam-Columns, Vol. 2, Space behavior and design, M c G r a w - H i l l Book Company, New York. [7] Colling, Francois., 1988, Estimation of the effect of different grading criterions on the bending strength of glulam beams using the "Karlsruhe calculation model", paper presented at I U F R O W o o d Engineering Group Meeting, T u r k u , F i n l a n d , June 1988. 59 Bibliography 60 [8] Dawe, P.S., 1964, The Effect of Knot Size on the Tensile Strength of European Redwood, W o o d 2 9 ( l l ) : p p 49-51. [9] Foschi, Ricardo 0 . and Barrett, D a v i d J . , 1980, Glued-Laminated Beam Strength: A Model Journal of the Structural Division, American Society of C i v i l Engineers, V o l . 106, N o . S T 8 , August, 1980, pp. 1735-1754. [10] Foschi, Ricardo 0 . , 1987, A procedure for the determination of localized modulus of elasticity, Holz als R o h - und Werkstoff 45 (1987), pp. 257-260. [11] Foschi, Ricardo O . and Barrett D a v i d J . , 1976, Longitudinal shear strength of Douglas-fir, Canadian Journal of C i v i l Engineering, V o l . 3, N o . 2, June 1976, pp. 198-208. [12] Foschi, Ricardo O . , 1977, Longitudinal shear strength in wood beams: a design method, Canadian Journal of C i v i l Engineering, V o l . 4, 1976, pp. 363-370. [13] Gjelsvik, Atle, 1981, The Theory of Thin Walled Bars, John W i l e y & sons, New York. [14] Glos,P., 1978, Berichte zur Zuverlassigkeitstheorie der Bauwerke: Zur Bestim- mung des Festigkeitsverhaltens von Brettschichtholz bei Druckbeanspruchung aus Werkstoff-und Einwirkungskenngrossen (Reliability Theory for Timber Structures: Determination of Compression Strength Behaviour of Glulam Components from Interaction of Material Properties), Heft 34/1978, L a b o r a t o r i u m fur den Konstructiven Ingenieurbau, Technische Universitat M i i n c h e n . [15] G o o d m a n , J . R . and J . B o d i g , 1970, Orthotropic Elastic Properties of Wood, P r o c . A S C E 9 6 ( S T l l ) : p p 2301-2319. [16] Johnson, J . W . and R . H . Kunesh, 1975, Tensile Strength of Special Douglas-fir and Hem-fir 2-inch Dimension Lumber, Wood and Fiber 6(4):pp 305-318. Bibliography 61 [17] K o k a , E x a u d N . , 1987, Laterally Loaded Wood Compression Members: Finite Element and Reliability Analysis, M . A . S c . thesis, October 1987, Department of C i v i l Engineering, The University of B r i t i s h C o l u m b i a , Vancouver, B . C . , Canada. [18] Kunesh, R . H . and J . W . Johnson, 1974, Effect of Size on Tensile Strength of Clear Douglas-fir and Hem-fir Dimension Lumber, Forest Product Journal 24(8):pp 3236. [19] Kunesh, R . H . and J . W . Johnson, 1972, Effect of Single Knots on Tensile Strength of2x8-inch Douglas-fir Dimension Lumber, Forest Product Journal 22(l):pp 32- 36. [20] M a l h o t r a , S . K . and S.J.Mazur, 1971, Buckling Strength of Solid Timber Columns, Trans. E n g . Inst, of Canada, published in E n g . Journal 13(A-4):I-VII. [21] M o e , J . , 1961, The Mechanism of Failure of Wood in Bending, Publication of International Association for Bridge and Structural Engineering, 21:163-178. [22] Nemeth,L.J., 1965, Correlation between Tensile Strength and Modulus of Elasticity for Dimension Lumber, P r o c . 2nd S y m p o s i u m on Non-Destructive Testing of W o o d , Spokane, Washington, U S A . [23] Norris, C . B . , 1955, Strength of Orthotropic Materials Subjected to Combined Stresses, U S D A Forest Sevices, F P L Report N o . 1816. [24] O'Halloran, M . R . , 1973, Curvillinear Stress-Strain Relationship for Wood in Compression, P h . D . thesis, Colorado State University. [25] Timoshenko S. and Woinowsky-Krieger S., 1959, Theory of Plates and Shells, M c G r a w - H i l l Book Company, New York. [26] Timoshenko, S.P.and Gere, J . M . , 1961, Theory of Elastic Stability, M c G r a w - H i l l Book Company, New York. Bibliography 62 [27] Timoshenko, S.P., 1953, History of Strength of Material, M c G r a w - H i l l Book Company, New York. [28] Taylor, S.E. and Bender D . A . , 1987, Modelling Localized Tensile Strength and MOE Properties in Lumber, paper presented at the 1987 International Winter Meeting of the American Society of A g r i c u l t u r a l Engineers, Paper No. 87-4519. [29] Y l i n e n , A . , 1956, A Method of Determining the Buckling Stress and Required Cross-sectional Area for centrally loaded Straight Columns in Elastic and Inelastic Range, I A B S A Publications, 16:529-549. [30] Zienckiewicz, O . C . , 1971, The Finite Element Method in Engineering Science, M c G r a w - H i l l Book Company, Londan, England.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Nonlinear finite element analysis of multilayered beam-columns
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Nonlinear finite element analysis of multilayered beam-columns Chandra, Rajesh 1989
pdf
Page Metadata
Item Metadata
Title | Nonlinear finite element analysis of multilayered beam-columns |
Creator |
Chandra, Rajesh |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | A finite element model has been developed in this thesis for predicting strength and stiffness behavior of multilayered beam-columns. The analysis incorporates material and geometric nonlinearities in order to determine the ultimate load carrying capacity. The finite element model takes into account the continuous variability of material properties along the length of layers so that multilayered wood beam-columns can be analyzed. Transverse as well as lateral bending in combination with axial tension or compression can be considered along with different layer configurations, various support and loading conditions. A computer program has been developed based on this formulation. Cubic beam elements have been used. Numerical integration of the virtual work equations has been carried out using Gauss quadrature. The resulting set of nonlinear equations is solved by using the Newton-Raphson scheme. Numerical investigations have been carried out to verify the results and test the capabilities of the program. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-30 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062878 |
URI | http://hdl.handle.net/2429/27917 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1989_A7 C43.pdf [ 2.81MB ]
- Metadata
- JSON: 831-1.0062878.json
- JSON-LD: 831-1.0062878-ld.json
- RDF/XML (Pretty): 831-1.0062878-rdf.xml
- RDF/JSON: 831-1.0062878-rdf.json
- Turtle: 831-1.0062878-turtle.txt
- N-Triples: 831-1.0062878-rdf-ntriples.txt
- Original Record: 831-1.0062878-source.json
- Full Text
- 831-1.0062878-fulltext.txt
- Citation
- 831-1.0062878.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062878/manifest