N O N L I N E A R FINITE E L E M E N T ANALYSIS OF M U L T I L A Y E R E D B E A M - C O L U M N S by Rajesh Chandra B . Tech. (Hons.), Indian Institute of Technology, Kharagpur, India A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF A P P L I E D SCIENCE in T H E FACULTY OF GRADUATE 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 in partial fulfilment of the requirements for an advanced degree at the University of Br i t i sh Columbia, I agree that the Library 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 depart-ment 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 Br i t i sh Columbia Vancouver, B . C . , C A N A D A V 6 T 1W5 October,1989 Abstract 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 vir tual 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. Table of Contents Abstract ii Tables v Figures vii Acknowledgement viii 1 Introduction and Literature Review 1 1.1 General Remarks 1 1.2 Objectives 5 1.3 Thesis Organisation 6 2 Formulation of Theory 7 2.1 General Assumptions 7 2.2 Kinemat ic Relationships 7 2.3 Stress-Strain Relationships 10 2.4 V i r t ua l Work Equation 18 3 Finite Element Formulation 21 3.1 Fini te Element Discretization 22 3.2 Problem Formulation 23 3.3 Interpolation Functions 25 3.4 V i r t ua l Work Equations 27 3.5 Method of Computation 30 i i i 3.6 Numerical Integration 32 3.7 Convergence Criterion 33 3.8 F E N M B C — t h e Computer Program 34 4 Verification and Numerical Results 35 4.1 F ixed Ended Beam Under Uniform Load • 35 4.1.1 Fini te Element Results 37 4.2 Other Test Cases 40 4.2.1 Simply Supported Multilayered Beam under B iax i a l Loading for Data-set I 41 4.2.2 Simply Supported Multi layered Beam under B iax ia l Loading for Data-set II 45 4.2.3 Eccentric biaxial loading on Simply Supported Beam for Data-set I 48 4.2.4 Simply Supported Beam-Column with axial loading having ec-centricity in y-direction for Data-set I 50 4.2.5 Simply Supported Beam-Column with axial loading having ec-centricity in y-direction for Data-set II 52 4.2.6 Simply Supported Multi layered Beam with a Nonlinear Com-pression Zone 54 5 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 Biax ia 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 Bazan (1980) 14 2.7 Stress-strain relationship used in the present study 16 2.8 Bil inear stress-strain relationships for various m as used by K o k a . . . 17 2.9 Exponential compressive stress-strain relationship used for F E N M B C 20 3.1 Fini te element discretization of the beam-column 23 3.2 Loca l coordinate system for the element and the window 24 3.3 Flowchart for iterative technique using Newton-Raphson Method . . . 31 3.4 Convergence criterion 34 4.1 Data 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 . . 39 4.3 Data for the test cases 40 4.4 Displacement results for biaxially loaded multilayered beams for data-set I 43 4.5 Stresses for biaxially loaded multilayered beams for data-set I . . . . 44 vi 4.6 Displacement results for biaxially loaded multilayered beams for data-set II 46 4.7 Stresses for biaxially loaded multilayered beams for data-set II . . . . 47 4.8 Displacement results for simply supported multilayered beams under eccentric biaxial loading for data-set I 49 4.9 Displacement results for simply supported multilayered beam-column with axial loading having eccentricity in y- direction for data-set I . . . 51 4.10 Displacement results for simply supported multilayered beam-column with axial loading having eccentricity in y-direction for data-set II . . . 53 4.11 Displacement results for simply supported multilayered beams with a nonlinear compression zone 55 4.12 Stress-strain curve for simply supported multilayered beam with a non-linear compression zone 56 vii Acknowledgement M y special thanks to my supervisor Dr . Ricardo 0 . Foschi, for his guidance and encouragement throughout the course of my research work and in the preparation of this thesis. I would also like to thank Felix Z. Yao, Bryan 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 Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Vlll To M y Parents K a m l a and Akhi lesh Chandra CHAPTER 1 Introduction and Literature Review 1.1 General Remarks The use of glued-laminated multilayered beams and beam-columns is quite common in timber buildings and could be used increasingly in 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. Also , reliability methods are being used in revising the timber codes in Canada and elsewhere. The 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. Although this kind of idealization has resulted in satisfactory designs, it may not allow for an optimum design. The main reason for this simplification has been the complexities involved in three dimensional analysis and the difficulties involved in their hand calculations. However, with computers being widely used in structural engineering, such analyses have become easier and more accurate. Beam-columns are often subjected to biaxial bending coupled wi th tension or compression. This 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). Although the two loading conditions 1 Chapter 1: Introduction and Literature Review 2 Figure 1.1: Biaxia 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 in 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 ig . 1.1(c). For example the bending moment in y-direction My may be applied first up to a certain value following path /, then the bending moment in z-direction M2 may be applied upto a certain value following path k, and then finally the axial loading P may be increased following path i unt i l the beam-column fails. The failure may also be caused following paths ghi, bed or any other path as shown in F ig . 1.1(c). The 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 wi l l vary depending on the load history. However, if P, Mx, My increase proportionally as radial loading, 4he two loading conditions are equivalent. In the present study, this proportional increase in loading has been used. Wood being an organic building material, the modulus of elasticity, tensile and compressive strengths have large variability, both within and across members. This variability in material properties of wood makes an accurate design procedure diffi-cult and the performance of such wood structures more uncertain. In recent years, this uncertainty has been addressed by using reliability based design procedures in wood structures. Timber codes are being revised to incorporate this change in design philosophy. The present study was undertaken to develop a general analysis tool for mul-tilayered beams and beam-columns, with 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. The model F E N M B C (Finite Element Nonlinear analysis of Multi layered Beam-Columns) developed in the present study was to overcome these limitations. It is based on a 3-dimensional analysis including axial deformation, biax-ial bending and torsion for multilayered beam columns. Both 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) in length. Each 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. The model uses a linear finite element method to obtain stresses in each cell, considering each one a five noded plate element. The model uses weakest link failure criterion to calculate the strength of the beam. The Foschi/Barrett model formed the basis for the development of the Karlsruhe model [7]. This finite element model takes into account the material nonlinearity of wood in 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 al . [3], Foschi [10], Taylor and Bender [28] and Col l ing [7] ). The simulation of material properties of built-up beam-columns is not the objective of the present study. The analysis assumes that these material properties of all laminations are given, or can be simulated along their length. The stress-strain relationship used herein are discussed in detail in 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. This was due to the lack of suitable connection details which prevented very high stresses from being developed in tension members. W i t h the availability of more effective connections, there was a renewed interest in tensile strength of timber. The phenomenon of size effect was observed by Kunesh and Johnson [18] (1974). They tested commercial sizes of clear Douglas-fir and Hem-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 con-tributed to the development of appropriate theory dealing wi th defects in timber. The 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 Kunesh [16] (1975). Attempts 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. This analysis would include nonlinear material and geometric behavior of the beam-column so that the ultimate load carrying capacity can be determined. The finite element model would take into account the continuous variability of material properties along the length. The loading would be axial , transverse and lateral as in a three dimensional beam-column. 2. Implementation of this analysis into a computer program which would allow flexibility in placing the layers in different configurations, various support condi-tions and loadings. It would also be user friendly and would require a minimum of data input. 3. Validation of this computer program by comparing its predictions with available analytical or experimental results. Chapter 1: Introduction and Literature Review 6 1.3 Thesis Organisation Chapter 2 describes the formulation of the theory used in the development of the mul-tilayered beam-column model. The general assumptions, kinematic and constitutive relationships are given in detail. Chapter 3 provides the outline and detail of the finite element formulation of the problem and the computer implementation. The vir tual 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. The 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 fur-ther research work in this area. CHAPTER 2 Formulation of Theory 2.1 General Assumptions The 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 discon-tinuity exists between layers. 4. Materia] properties may have continuous variation along the span of a layer, however, material properties are assumed to be constant within a small division of an element (called a window as explained later in Section 3.1). 5. The stress-strain relationship has been assumed to be linear elastic with brittle failure in tension, and nonlinear having a falling branch with plastic failure in compression as shown in F ig . 2.7. 2.2 Kinematic Relationships The orientation of the beam-column with respect to a Cartesian coordinate system is shown in F ig . 2.1. The 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 The 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: dv dw d6 1 dx v-zO uA - u-y-WA = w -f y6 dx ^ ^ dx (2.1) (2.2) (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 term in Eq.( 2.1) by replacing w and z wi th v and y respectively. The last term in Eq.( 2.1) has been included to account for warping of the cross-section due to the rotation 9. The nonlinear strain-displacement relationship used in the present study is now presented. From the geometry of the F i g . 2.3, after deformation, the change in contour length ds can be expressed as j 2 ^ J 2u, . 9UA.2 (dvA 2 dwA 2 ds Kdx [(1 + -^-) + ( ^ ) (2.4) Expanding it binomially and neglecting the higher order terms, ds can be written as (2.5) j , h , duA l,dvA2 1 dwA2 d s = d x [ l + _ _ + _ ( _ _ ) + - ( — ) ] Thus the axial strain e can be expressed as e = ds — dx duA 1. dvA 1 ,dw As2 dx ^ ^ *~> ^ 2{ dx ' dx 2 dx (2.6) Chapter 2: Formulation of Theory 10 Figure 2.3: Large deformation relationships Substituting uAlvA and wA from Eqs.( 2.1), (2.2) and (2.3), the strain displacement for a general point in the beam-column is given by du d 2v d 26 d 2w dx dx 2 dx 2 dx 2 1 dv d6 x , 1 dw d6 ^ , 2.3 Stress-Strain Relationships While wood is assumed to behave elastically in tension as discussed earlier in Sec-tion 1.1, many uniaxial stress-strain relationships have been proposed by various re-searchers to model the nonlinear behaviour of wood in compression. These are mostly empirical in nature, attempting to theoretically represent the stress-strain relation-ships 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 in F i g . 2.4 has been used quite often to calculate the ultimate load carrying capacity of wood mem-bers. Compression Tension (a) stress-strain relationship (b) stress distribution Figure 2.4: Bilinear elasto-plastic stress-strain relationship Yl inen [29] (1956) proposed a simple bilinear relationship, which was used by Malhot ra and Mazur [20] (1971) in their research work. This relationship can be given as e=±[co--(l-c)fcln(\-jJ] (2.8) where e is the normal strain, o is the normal stress, fc is maximum compressive stress, E is the ini t ia l modulus of elasticity and c is a parameter depending on the shape of the curve. As 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 12 Ylinen(1956) fc tan0 E=tana O'HalloranCWS) Qlos(1978) (a) (b) (c) Figure 2.5: Various stress-strain relationships Malhot ra and Mazur [20] have used this relationship extensively to predict the buckling strength of solid timber columns and described the curve as a good approx-imation 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. Goodman and Bodig [15] (1971) carried out extensive tests on clear dry wood in compression parallel to the grain and reported the test results. O'Hal loran [24] (1973) has used these data to propose a mathematical equation for the stress-strain curve at various grain angles and grain orientations. This equation has been given in 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—- (2.10) r — 1 a = m 4 ^ < 2-»> where fc is the maximum stress. A typical plot of the fitted curve is shown in F i g . 2.5(b). The equation cannot be used beyond the maximum compressive stress / c , because it drops rapidly to negative 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 maximum 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 stress-strain 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 in 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 maxi-mum 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). On the basis of extensive experimental testing, Glos proposed a curve of the shape shown in F ig . 2.5 which is similar to that reported for clear wood by Bechtel and Norris [4] (1952), Moe [21] (1961) and others. The curve is characterized by four parameters as shown in F ig . 2.7. The 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 Bazan (1980) where G4 = Is 6 / 5 ( 1 - £ ) G 2 = E Is where / is the stress, e is strain, E is modulus of elasticity, fc is maximum compressive stress, fs is the asymptotic compression stress for large strain and ex is the strain at maximum 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. Mul t ip le curvilinear regression techniques have been used to ex-press 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 Bazan [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 = E0e - [E0e + | / c | ( l + m) + m £ 0 e ] ( l - A ( e + |e c |)) (2.13) where a is stress, EQ is elasticity modulus, e is strain, fc is maximum compressive stress, m is the slope of the falling branch on compression side, ec is the maximum strain corresponding to the maximum compressive stress fc and A ( e + |e c |) is the step function defined as follows e > - | e c | =• A ( e + | e e | ) = 1 (2.14) e < -\ee\ A ( e + |e e |) = 0 (2.15) The bilinear relationship for different values of m has been shown in F i g . 2.8. A further extension of this relationship was made to obtain a more realistic ap-proximation of the stress-strain behaviour of wood in 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 maximum tensile stress Ft, upon which a brittle failure occurs. a = Ee - [Ee - F{e)]{l - A{e)) The function F(e) for e < 0 has the form F(e) = [-\Po\-m1Ee](l-eft\) (2.16) (2.17) Chapter 2: Formulation of Theory 16 l«cl Plastic failure Brittle failure IPfJ Figure 2.7: Stress-strain relationship used in 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 in 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(£) (2.18) 1, i f £ > 0 0, if £ < 0 In Eq.( 2.16) only the three parameters E, fc and mj need be specified. The re-maining parameters of the equation can be calculated in the following manner. From Eq.( 2.16), for e = e c the following expression can be readily obtained [-\p0\-m1Eec)(l-e%\) = -\FJ (2.19) Differentiating Eq.( 2.16) with respect to e and substituting e — £ c , the following Chapter 2: Formulation of Theory 17 c j 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 - mxE(\ - e - l ^ k c | ) + (|po| - mxE\ec\)^-e-^^ bo I = 0 (2.20) From Eq.( 2.19), Fc bo | - m\E\ec (2.21) Defininig x = e (2.22) and taking the logarithm of both sides, Inx = E l l - b o T | £ c l (2.23) Substituting x into Eqs.( 2.20) and (2.21), and rearranging terms, —miE(\ — x) + Fc E — — — x - 0 1 - x \p0\ (2.24) Since from Eq.( 2.24), |p 0 | c a n he written as bol = \FC X (2.25) mx(l — x) 2 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 maximum strain |e c | corresponding to the maximum stress Fc can be written as | £ c | = _ % i / n z (2.27) hi 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 + 2mi + milnxi A typical stress-strain plot used in F E N M B C has been shown in F i g . 2.9, correspond-ing to E = 1 4 , 0 0 0 M P a , m ! = 0.25, Fc = UQMPa. The stress-strain relationship used in this study is a better representation of exper-imental results, as it can have a continuous curvilinear and smooth peak as opposed to the abrupt change of curvature in 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 expo-nential stress-strain relation, as given by O'Hal loran [24], has been taken care of by the variable slope m i . 2.4 Virtual Work Equation The Principle of Vi r tua l Work has been used to determine the governing equilibrium equations of the beam-column. If we apply a vir tual displacement 6a to the system, Chapter 2: Formulation of Theory 19 compatible with boundary conditions, the resulting internal work (6U) and external work (8W) done by the stresses and the external forces respectively are given as SU = jv8e TodV (2.30) 8W = Js8a TFdS (2.31) 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. The principle of vir tual work states that internal and external work must be equal: J 8e TadV = Jv 8a TFdV (2.32) Chapter 2: Formulation of Theory 20 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) Strains (m/m) (0.005) Figure 2.9: Exponential compressive stress-strain relationship used for F E N M B C CHAPTER 3 Finite Element Formulation The basic concept underlying the finite element method is that a structure can be modelled analytically by subdividing it into a finite number of elements. W i t h i n each finite element, a set of functions are assumed to approximate the stresses or displacements in that region. The set of approximating functions contain unknown parameters and are chosen to ensure continuity throughout the system. Appl icat ion of the principle of vir tual work results in a system of algebraic equations for the parameters in the approximating functions. Although finite element formulation can be based either on stress fields or displace-ment 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. The 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 Formulation 22 6. Calculation of strains and stresses based on nodal displacements. In case of nonlinear finite element analysis, an incremental/iterative technique is re-quired. The 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 wi l l be discussed in detail later in this chapter. 3.1 Finite Element Discretization The beam-column arrangement and the finite element discretization is shown in F i g . 3.1. The span is subdivided into elements of length A . The layers for each element are further subdivided into windows of length Ato.The numbering of layers has been done from top to bottom. The elements are numbered from left to right and so are the windows inside every element. The nodes are also numbered from left to right. The 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. The same could be achieved by employing a larger number of elements without windows, but then a larger stiffness matrix would be obtained. Thus while keeping the size of the stiffness matrix wi thin reasonable limits, the concept of windows permits consideration of the continuous variation of the material properties along the span. This results in 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 nonlineari-ties, 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 y.v 1 2 3 ,;*w, — « N| * • • : 4 SC I— ' / / / r 'z,w a typical node A layers x,u Figure 3.1: Fini te element discretization of the beam-column 3.2 Problem Formulation A beam element of length A with two end nodes, as shown in F ig . 3.2, is used in this formulation. A local coordinate <f (— 1 <<f < l ) i s used for each element along the x-axis. The x-coordinate of any point inside the element can be expressed as A , (3.1) where xc is the x-coordinate of the center of the element length. The element has been divided into a number of windows (NW) of equal length, hence A Aw -NW (3.2) A second local coordinate 7;( — 1 < 77 < 1) is used within a window. A transforma-tion 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. Chapter 3: Finite Element Formulation 24 4 -+i node 1 <> node 2 4W r - 1 • + 1 Figure 3.2: Local coordinate system for the element and the window Thus, for the i th window, F ig . 3.2, A A xi = xc - — + (i - l ) A t u = z c + — 6 %2 = %c ~ + *Aio = a:c + y f 2 A Aw A, XQ = x c - — + IAW — = xc + —£o From the above we can calculate £i,Cj2 and £o as = | ( . t - l ) A « ; - l 6 Also —iAw — 1 A 2 A A W 1 — tAw 1 A A fc - 6 Au> 2 A Hence we can write £ in terms of 77 as 2 A u ; A w A Atw - 1 ) + ^ (3.3) (3.4) (3.5) (3.6) Chapter 3: Finite Element Formulation 25 The 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 = {uuu[,vuv[,wuw[,91J[,u2,u2,v2,v2,w2,w'2,92,92} (3.7) 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 in the strain-displacement relationship. Also the displacement function must be continuous within the element and the displacement must be compatible between adjacent elements. The 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 sec-ond derivative in lateral displacements and twist and the first derivative in 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. This 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 approx-imation 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 du u = N L U I + N2{—)i + N3UJ + NA{—)j (3.8) dv dv v = N^ + ^i—^ + NzVj + N^ — h (3.9) dx dx w = + N2(^){ + N3Wj + N*^), (3.10) 9 = A ^ , . + + JV 3 0 j + J V 4 ( £ ) > (3.11) where subscripts i and j refer to the first and second nodes respectively of each element. Also N i = ( 3 - 1 2 ) ^3 = ( | + f f - ^ 3 ) (3-14) N4 = | ( - i - £ + ^ + f ) (3.15) Alternatively we can write these as u = P T ( 0 a « ' = P i T ( 0 a (3.16) t, = Q T ( £ ) a i ^ Q ^ a t>" = Q 2 T ( c > (3.17) u; = R r ( 0 a w' = R i T ( e ) a = R 2 T ( 0 a (3.18) 6 = S r ( £ ) a 0' = S i T ( O a 0" = S 2 r ( O a (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 = {NuN2l0,0,0,0,0,0,^3,^4,0,0,0,0,0,0} (3.20) P x 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) Chapter 3: Finite Element Formulation 27 Q i T = { 0 , 0 , ^ , ^ , 0 , 0 , 0 , 0 , 0 , 0 , ^ , ^ , 0 , 0 , 0 , 0 } (3.23) Q 2 T = {0,0, ^ " , ^ ' , 0 , 0 , 0 , 0 , 0 , 0 , ^ ' , ^ ' , 0 ,0,0,0} (3.24) ... etc. where Nx, N2, N3, N4 are as given earlier and N{, N^, N!$, N'4 are the first derivatives and N", N?, N£, N% are the second derivatives with respect to x. 3.4 Virtual Work Equations From Eq.( 2.7), the strain can now be written in the following form e = [ P i r -yQ2 T + y z S 2 T - z R 2 T ] a + ^ a T [ ( Q i - * S i ) ( Q i r - 2 S X T ) + ( R i + y S i ) ( R i T + y S i r ) ] a (3.25) = B o T a + i a T B i a (3.26) where B 0 = P i - 2/Q2 + yzS2 - 2R2 (3.27) B1 = ( Q i - 2 S 1 ) ( Q i T - z S 1 T ) - r ( R i + j / S 1 ) ( R 1 T - r j / S 1 T ) (3.28) B o corresponds to the linear part of e and B i to the nonlinear. App ly ing a virtual displacement 6a, as a —* a + A<5a, the following can be obtained e(A) = B 0 T ( a - f A<5a) + ^ ( a T + A<5a T )B 1 (a + A6a) = B 0 T a + A B 0 T 6 a + ^ a T B i a + ^ ( 5 a T B i a Zj Li + ^ a r B i * a + y<5a r Bi<5a Thus the vir tual strain can be written as 6e = (^)x-.0 = B o ^ a + ^ B i a + i a ^ a = B 0 r <5a + i a T B i T ( 5 a + ^a T Bi<$a Chapter 3: Finite Element Formulation 28 = ( B 0 r + ^ a r B x + ^ a r B 1 r ) 6 a = (B0T + aT(^±^))8a (3.29) and since B i is symmetrical (refer E q . 3.28) i.e. B i = B i r , Eq.( 3.29) can be written as 8e = ( B 0 T + a T B 1 )<5a (3.30) The stress a can be written as ( from E q . 2.16 ) a = £ ( B 0 T + ^ a T B i ) a + G(£) (3.31) where G(e) = -[Ee - {-\Po\ - m^e}^ - e W £ ] ( l - A(e)) (3.32) da can be written as da = E{B0T + a T Bi)<5a + H{e)(B0 T + a TB1)8a (3.33) where dG{e) H(e) de = ( - £ - m , £ ( l -eMic~Co)) - ( - | p o | - mxE(e - e0))^e^ (e- £o) (3.34) IPol Hence the vir tual work equation can be written as Jv 8e TadV - J 8a TFdS = 0 or f 6a T{B0 + &xa)c-dV - [ 8a TFdS = 0 (3.35) Jv Js Since this must hold for any 6a, j ^ ( B 0 + B i a ) < 7 < n / - JsFdS = 0 (3.36) Chapter 3: Finite Element Formulation 29 Thus we can write a function 0(a) as 0(a) - Jv{B0 + Bla)adV - JsFdS (3.37) The solution vector a must satisfy 0(a) = 0. To find this solution an iterative proce-dure is used. The 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 (3.38) da or A a = 0(a + A a ) - 0 ( a ) [ ^ ^ ] - 1 (3.39) aa Setting 0(a + A a ) = 0, A a = a 1 + 1 - a, = (3.40) or a , + 1 = a, - [K r]-V(a,) (3.41) Since 60(a) = J BxSaLardV + J (B0+ B1&)6adV = JvaB16a +Jv[(B0 + B1a)(E + H(e)}(B0 T+ a TB1)6a (3.42) Hence the tangent stiffness matrix can be written as [ K T ] = Jv[BiEBo Ta+^B1Ea TB1a + B0EBo T+ B0Ea TB1 + B i a £ B o T + B i a £ a r B i + B 0 H ( e ) B 0 T + B o # ( e ) a r B ! + B 1 a / / ( e ) B 0 r + B 1 a J r Y ( e ) a T B 1 + BxG(e)}dV (3.43) Chapter 3: Finite Element Formulation 30 This 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. Each value in the incremental solution vector A a is compared against an accept-able 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 The 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 matrix is utilised and only the band in 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 31 Q Start ^ Read geometry, material properties and load characteristics I Initialise [K^ = 0 ; {R} • 0 {u } = 0; iter = 0 r = r+1; increment the load vector {R} = {R} + {R} r r-1 r 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 Check convergence \ converged^ rice^-If there is no further load increment f Stop no Solution displacement vector {X} Calculate stresses {o} Figure 3.3: Flowchart for iterative technique using Newton-Raphson Method Chapter 3: Finite Element Formulation 32 3.6 Numerical Integration Recourse must be made to numerical integration methods to obtain the matrix [ 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 polynomial of degree (2n — 1) could be constructed and exactly integrated. The 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 (3-45) fc=ij=ii=i where n l , n2, n3 are number of integration points in x, y and z directions respectively. In the present case the integral I is of the form 1 = Jo J-zJo f( x>yi z) dxdy dz ( 3- 46) L / - i L i v ' (WWt (3-47) AwBT /•! 2 2 2 AwBT nl n2 n3 , 9 9 E E E W J ( ( , ^ a ) (3-48) z z z k=\j=i i=i 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 polynomial 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 in 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. This has been done because the program can handle multilayered beam columns and a number of windows within an element, to take care of the continuous variation of material properties along the span. For example, if there is a beam column with 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 in 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 in various directions. F ina l ly 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 ar is the solution vector at the end of r th iteration and a r + 1 is the solution at the end of (r + l) tk iteration in F i g . 3.4, then A x = |a r+i — ar \ represents the difference between the two displacement vectors. | a r | and |a r +i | can be given as follows NDOF K l 2 = E «?(«') i=i NDOF :=1 where ar(i) and aT+i(i) are components of ar and a r + i respectively. Chapter 3: Finite Element Formulation 34 a r+1"a r ^ ^ ^ ^ a r + 1 0 Figure 3.4: Convergence criterion Then A a can be written as A a = E y/(ar+i{i)-ar(i)y i=l (3.49) The convergence criterion for this norm is given as A a — T T < specified tolerance \ar\ (3.50) 3.8 FENMBC—the Computer Program A computer program called F E N M B C (Finite Element Nonlinear analysis of Multilayered Beam Column) 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, experi-mental or numerical results published in the literature, whenever possible. The cases presented below involve geometric nonlinearity, material nonlinearity, biaxial bend-ing, 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 Premium ( I B M P C - A T compat-ible) at the University of Br i t i sh 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. Each window can have a maximum of 5 Gauss points in the three principal axes directions. The 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. The 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: Data for Timoshenko's fixed ended beam bending problem with 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 = (a2 ~ x2){b2 - y2)x{b00 + b02y2 + b20x2 + b22x2y2) (4.1) v = (a2-x2)(b2 -y2)y{co0 +Co2y2 + c20x2 + c22x2y2) (4.2) w = (a2-x2)2(b2-y2)2{a00 + a02y2+ a20x2) (4.3) The values of all the parameters were estimated for various intensities of the load q and for different aspect ratios | , assuming v = 0.3. The results reported for an infinitely long plate (for | = 0) have been used to compare with the results of F E N M B C . Chapter 4: Verification and Numerical Results 37 4.1.1 Finite Element Results The program F E N M B C was run using the data shown above with some additional data required for the program: number of layers nlayer = 1, number of windows nwindo — 3, number of Gauss points along span in x-direction ngauss = 5, number of Gauss points along the width in y-direction ngausy = 2, Number of Gauss points along the depth in z-direction ngausz — 5, slope of the compression curve m\ = 0.2, and high values for maximum allowable tensile and compressive stresses of Fc = 6.0 x l0 15N/m 2 and Ft = 6.0 x 1 0 1 5 i V / m 2 to ensure full elastic behavior. The distributed load in the z-direction was given as qz = 10.7315./V/m and the number of load steps was specified as 16 for a smooth load-displacement relationship. The results have been tabulated in Table 4.1 and subsequently shown graphically in 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, q is the uniformly distributed load, / is the span, D — EI, 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 wmax is the maximum displacement in z- direction. From 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 Timoshenko FENMBC FENMBC FENMBC theoretical nonlinear nonlinear nonlinear nonlinear 6 elements 10 elements 20 elements Dh ^ Wrnnr j 0.0 0.0 0.0 0.0 0.0 0.0 10.0 0.4166667 0.40 0.4081488 0.4082008 0.4082064 20.0 0.8333333 0.62 0.6828256 0.6832096 0.6832456 30.0 1.2500000 0.84 0.8791920 0.8798320 0.8799200 40.0 1.6666667 0.98 1.0311120 1.0321920 1.0323440 50.0 2.0833333 1.10 1.1571440 1.1586880 1.1589200 60.0 2.5000000 1.21 1.2654880 1.2675280 1.2678400 70.0 2.1966667 1.31 1.3610000 1.3635520 1.3639520 80.0 3.3333333 1.39 1.4467600 1.4498240 1.4503200 90.0 3.7500000 1.47 1.5248480 1.5284240 1.5290160 100.0 4.1666667 1.54 1.5967040 1.6007920 1.6014800 110.0 4.5833333 1.62 1.6634000 1.6680080 1.6688000 120.0 5.0000000 1.67 1.7257520 1.7308720 1.7317680 130.0 5.4166667 1.73 1.7861440 1.7917760 1.7927760 140.0 5.8333333 1.79 1.8413200 1.8474560 1,8485600 150.0 6.2500000 1.84 1.8937040 1.9003360 1.9015520 160.0 6.6666667 1.90 1.9436160 1.9507440 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 40 4.2 Other Test Cases Other cases were run and the results have been reported below. These have been compared with the linear beam bending and beam column theories for lack of available nonlinear results. The same two data sets were used throughout the investigation unless otherwise specified. The first data set contained 10 layers wi th their elasticity modul i placed symmetrically wi th respect to the y-axis on the beam cross-section as shown i n F i g . 4.3. The second data set contained 10 layers wi th their elasticity modul i placed asymmetrically with respect to the y-axis of the beam cross-section as shown in F i g . 4.3. M.O.E.(kN/m*) M.O.E.(kN/m») 0.14e+08 0.066+08 0.12e+08 0.066+08 E 0.106+08 E 0.086+08 E 9 © 0.086+08 ^ 0.066+08 0.086+08 0.106+08 in k_ Q> 0.066+08 |2 0.106+08 >. 0.086+08 J2 0.126+08 o 0.106+08 0.126+08 0.126+08 0.146+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/m2 Maximum allowable comp. stress FQ =0.4e+08 N/m2 Slope of the comp. curve m == .25 Figure 4.3: Data for the test cases Chapter 4: Verification and Numerical Results 41 4.2.1 Simply Supported Multilayered Beam under Biaxial Loading for Data-set I The case of a multilayered beam with data-set I as shown in F i g . 4.3 with concen-trated loads in y- and z-directions (refer F i g . 4.4) was run. The additional data re-quired for running F E N M B C were: nelem = 10, nwindo = 3,ngauss = 3,ngausy = b,ngausz = 5. The maximum tensile and compressive stress values was taken as 4 . 0 x l 0 l o / V / m 2 for al l the layers.The concentrated loads were Py = l . O x 10 6 A^andP 2 = 1.0 x 1 0 6 N . The displacement results have been reported in F i g . 4.4. The total tensile 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 in F i g . 4.5. The linear deflections v and w were calculated from the linear beam bending the-ory using the method of transformed cross-section for comparison with the nonlinear results obtained from F E N M B C . The nonlinear displacements v and w are the max-imum displacements (at node no. 6) and are the same for al l the layers. As both the ends are hinged, the displacements converge to a lower value from the linear dis-placement values due to stiffening of the beam resulting from the large deformation terms. The rotation 0, as obtained in the result shown in F ig . 4.4, needs some explana-tion. 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. The loadings in transverse and lateral directions are applied incrementally. The load in y-direction, APy, produces a displacement Av. As the first load increment APZ is applied, there exists a moment A M . = APZ x Av, which causes a rotation in the beam cross-section along x-axis. Also the force APZ produces the displacement Aw in z-direction. As 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 S2, these being the derivatives of the shape functions for 6. The odd terms in y and z after integration cancel out due to the symmetry of layers, however, the even terms add up and produce the desired rotational deforma-tion. The total tensile stresses are seen to be linear as indeed they should be, since they are in the tensile zone below the maximum 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 45 4.2.2 Simply Supported Multilayered Beam under Biaxial Loading for Data-set II The other multilayered beam with data-set II as shown in F i g . 4.3 with concen-trated loads in y- and z-directions (refer F ig . 4.6) was run. The additional data re-quired for running F E N M B C were: nelem = 10, nwindo — 3,ngauss = 3,ngausy = 5,ngausz — 5. The maximum tensile and compressive stress values was taken as 4.0 x 1 0 l o / V / m 2 for all the layers.The concentrated loads were Py = 1.0 x 10 6./VandP* = 1.0 x 1067V. The displacement results have been reported in F i g . 4.6. The total tensile 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 in F ig . 4.7. The linear deflections v and w were again calculated from the linear beam bending theory using equivalent cross-section. The transverse and lateral deflections resulting from the rotation 6 were ignored in this calculation due to their small magnitude. The nonlinear displacements v and w again converge to a. lower value from the linear displacement values as discussed earlier. The rotation 6 in 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 in the tensile zone below the maximum tensile stress Fj. Chapter 4: Verification and Numerical Results 46 Figure 4.6: Displacement results for biaxially loaded multilayered beams for data-set II Chapter 4: Verification and Numerical Results 47 Stresses (MPa) Strains (m/m) Figure 4.7: Stresses for biaxially loaded multilayered beams for data-set II Chapter 4: Verification and Numerical Results 48 4.2.3 Eccentric biaxial loading on Simply Supported Beam for Data-set I The multilayered beam with data-set I as shown in F i g . 4.3 wi th eccentrically located concentrated loads in y- and z-directions (refer F i g . 4.8) was run. The additional data required for running F E N M B C were the same as in Section 4.2.1 with eccentricities ey = 20mm and ez = 20mm. The displacement results have been reported in F i g . 4.8. The nonlinear displacements v and w again converge to a lower value from the linear displacement values as expected. A t the beginning the rotation 0 is zero as are the loads applied. However as the first load increment is applied, 0 takes on a negative value. As can be seen from the loading in 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. As the load is increased, the rotation also gets affected by the displacements v and w as explained in Section 4.2.1. The interaction of these moments produces the rotational deformations as shown in F i g . 4.8. The stress-strain diagram has not been included here as the tensile and compressive stresses are sti l l 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 50 4.2.4 Simply Supported Beam-Column with axial loading having eccentricity in y-direction for Data-set I The 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 mod-ulus E. The Euler buckling load thus calculated was PE = 173.49&./V. A n axial load was applied with a very small eccentricity in the y-direction, and increased in steps of 8.6754kN. As 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. The results have been reported in F ig . 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 52 4.2.5 Simply Supported Beam-Column with axial loading having eccentricity in y-direction for Data-set II The program F E N M B C was further tested with data-set II for buckling and the results were compared with 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 in the y-direction. The results have been reported in F i g . 4.10. Chapter 4: Verification and Numerical Results 53 2.5 v (nonlinear) — B — w (nonlinear) ---A---© (nonlinear) G 50,000 100,000 150,000 P 200,000 Axial Compressive Load P (N) Figure 4.10: Displacement results for simply supported multilayered beam-column with axial loading having eccentricity in y-direction for data-set II Chapter 4: Verification and Numerical Results 54 4.2.6 Simply Supported Multilayered Beam with a Nonlin-ear Compression Zone The data-set I was used to test the response of the multilayered beam in the nonlinear compression zone. In the data-set I, the maximum allowable compressive stress values were divided by 100 to force the compressive stresses to fall in the nonlinear range. The 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. This case was run to test the response of the multilayered beam in the nonlinear range. The 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 with a nonlinear compression zone Chapter 4: Verification and Numerical Results 56 (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 non-linear compression zone CHAPTER 5 Summary, Conclusions and Scope for Future Research 5.1 Summary and Conclusions A direct vir tual work formulation of multilayered beam-columns is developed incor-porating geometric as well as material nonlinearities. F in i te element total and incre-mental equilibrium equations are derived based on the assumed displacement model. This method has been used to examine the strength and stability of glued-laminated beam-columns. The computer program developed has been verified and tested for several beam-column problems. The generality of the model makes it a powerful and versatile tool for solving a large variety of problems. The most important question for the user of the finite element method is whether the method yields sufficiently accurate results for his purpose. As 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 in the program. As in any other approximate numerical method, engineering judgement should be exercised in interpreting the results obtained by the finite element method since ac-curacy of the results depends on the assumptions employed in 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 The present model does not include shear deformations in 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 in 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] Abramowitz , 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, Mad 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.R1910. [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 Br i t i sh Columbia, 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 be-havior and design, M c G r a w - H i l l Book Company, New York. [7] Col l ing, 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 Wood Engineering Group Meeting, Turku, Finland, June 1988. 59 Bibliography 60 [8] Dawe, P.S., 1964, The Effect of Knot Size on the Tensile Strength of European Redwood, Wood 29( l l ) :pp 49-51. [9] Foschi, Ricardo 0 . and Barrett , David J . , 1980, Glued-Laminated Beam Strength: A Model Journal of the Structural Divis ion, American Society of C i v i l Engineers, V o l . 106, No. ST8, August, 1980, pp. 1735-1754. [10] Foschi, Ricardo 0 . , 1987, A procedure for the determination of localized modulus of elasticity, Holz als Roh- und Werkstoff 45 (1987), pp. 257-260. [11] Foschi, Ricardo O. and Barrett David J . , 1976, Longitudinal shear strength of Douglas-fir, Canadian Journal of C i v i l Engineering, V o l . 3, No. 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, At le , 1981, The Theory of Thin Walled Bars, John Wi ley & 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 Struc-tures: Determination of Compression Strength Behaviour of Glulam Components from Interaction of Material Properties), Heft 34/1978, Laboratorium fur den Konstructiven Ingenieurbau, Technische Universitat Miinchen. [15] Goodman, J .R . and J . Bodig , 1970, Orthotropic Elastic Properties of Wood, Proc. 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 , Exaud 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 Br i t i sh Columbia, 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 32-36. [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] Malhotra , S .K. and S.J .Mazur, 1971, Buckling Strength of Solid Timber Columns, Trans. Eng. Inst, of Canada, published in Eng. Journal 13(A-4):I-VII. [21] Moe , J . , 1961, The Mechanism of Failure of Wood in Bending, Publ icat ion of International Association for Bridge and Structural Engineering, 21:163-178. [22] Nemeth,L.J . , 1965, Correlation between Tensile Strength and Modulus of Elas-ticity for Dimension Lumber, Proc. 2nd Symposium on Non-Destructive Testing of Wood, 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 No. 1816. [24] O'Hal loran, M . R . , 1973, Curvillinear Stress-Strain Relationship for Wood in Compression, Ph . 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 Agricul tura 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 Inelas-tic 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062878/manifest