C O M P U T A T I O N A L M O D E L I N G O F S T R A N D - B A S E D W O O D C O M P O S I T E S I N B E N D I N G by P E G G I L Y N N C L O U S T O N B .A .Sc , The University of British Columbia, 1989 M.A.Sc., The University of British Columbia, 1995 A THESIS S U B M I T T E D I N P A R T I A L F U L F I L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S • (Facu l ty -o f Forestry) We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH COLUMBIA October 2001 © Peggi Lynn Clouston, 2001 UBC Special Collections - Thesis Authorisation Form http://www.l ibrary.ubc.ca/spcoll / thesauth.html In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission for extensive copying of t h i s thesis for s c h o l a r l y purposes may be granted by the head of my department or by his or her representatives. It i s understood that copying or p u b l i c a t i o n of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. The University of B r i t i s h Columbia Vancouver, Canada Date 1 of 1 07-Oct-01 15:40 A B S T R A C T A stochastic finite element approach is presented herein for simulating the nonlinear behaviour of strand-based wood composites with strands of varying grain angle. The approach is based on the constitutive properties of the individual strands making it practical and versatile in application. It can be used to gauge the effects of varying strand characteristics in product manufacturing or it may be useful in the design of wood composite structures. Both 2 and 3 dimensional, stochastic, materially nonlinear finite element codes are developed. The nonlinear constitutive behaviour of the wood strands is characterized within the framework of rate-independent theory of orthotropic elasto-plasticity. The constitutive model employs the Tsai-Wu yield criterion with the associated flow rule of plasticity. Failure is marked by an upper bound surface whereupon either perfect plasticity or an abrupt loss of strength and stiffness ensues. This constitutive model is implemented into the finite element codes. The programs are based on the conventional displacement formulation using linear isoparametric elements. The nonlinearities are resolved by a modified Newton-Raphson procedure. The programs are further formulated in a probabilistic manner using random variables as input for principal material strengths and stiffnesses. To generate entire data samples for comparison with experimental samples, the programs were written with extended capacity to perform Monte Carlo simulations. The mechanical properties of the strands are derived through both experiment and analysis. An experimental database of principal material strengths and stiffness parameters is acquired for Douglas fir heartwood strands. Statistical parameters for shear strength and stiffness as well as the interaction parameter of the Tsai-Wu criterion are estimated, however, through a least square minirriization of error between simulated and experimental compression strength of [ +15]sangle-ply laminates. Weibull weakest-link theory is employed to adjust experimental tensile strength values for size effect. The general performance of the programs is verified through comparison of results for several analyses solved using analytical techniques or alternate programs. Following this, the ability of the models to reproduce experimental findings for angle-ply laminates in tension, compression and 3 point bending is validated. A preliminary investigation is conducted to compare numerical simulations with experimental data for Parallam* in tension and 3 point bending . The favourable comparisons of the model to experimental results attest to the effectiveness of the proposed technique. /// T A B L E O F C O N T E N T S Abstract • ii Table of Contents iv List of Figures viii List of Tables xiii Notation xv Acknowledgments xix Chapter 1 - Introduction 1 1.1 Introduction 1 1.2 Model Considerations 3 1.3 Objective and Scope 5 Chapter 2 - Background 7 2.1 Introduction 7 2.2 Failure Theories 8 2.2.1 Tsai-Wu Criterion 10 -2.2.1.1 Significance of the Interaction Term, F 1 2 12 2.2.1.2 Evaluation of the Interaction Term, F 1 2 13 2.3 Literature Review of Wood Composite Modeling 17 Chapter 3 - Constitutive Model 21 3.1 Introduction 21 3.2 Constitutive Model 22 3.2.1 Linear Elastic Regime 22 3.2.1.1 Constitutive relations for off-axis coordinates 25 3.2.2 Elasto-plastic Regime 29 3.2.2.1 Initial Yield Surface 29 3.2.2.2 Hardening Rule 31 3.2.2.3 Flow Rule 34 iv 3.2.3 Post Failure Regime 40 3.3 Commentary on Plane Stress Formulation 40 Chapter 4 - Finite Element Formulation 42 4.1 Introduction 42 4.2 Two-Dimensional Linear Finite Element Model ; 43 4.2.1 Displacement-Based Element Representation 43 4.2.2 Stiffness Matrix Formulation 44 4.2.2.1 Compatibility 44 4.2.2.2 Equilibrium 46 4.2.2.3 Constitutive Relations 46 4.2.3 Two-Dimensional Model Limitations 49 4.3 Three-Dimensional Linear Finite Element Model 50 4.3.1 Displacement-Based Element Representation 50 4.3.2 Stiffness Matrix Formulation 51 4.3.2.1 Compatibility 51 4.3.2.2 Equilibrium 52 4.3.2.3 Constitutive Relations 52 4.4 Nonlinear Formulation 54 4.4.1 Modified Newton Raphson Method 54 4.4.2 Elastic-Plastic Formulation 57 4.4.3 Post-failure Regime 63 4.5 Program Verification 64 4.5.1 Isotropic Cantilever Analysis 65 4.5.2 5-Point Continuous Beam Analysis 68 4.5.3 AnglePly Laminate Analysis 71 Chapter 5 - Strand Property Database 75 5.1 Introduction 75 5.2 Material 75 5.3 Experimentally Derived Properties 76 v 5.3.1 Compression Tests 76 5.3.1.1 Material Preparation 76 5.3.1.2 Test Method 78 5.3.1.3 Experimental Results 79 5.3.2 Tension Tests 89 5.3.2.1 Material Preparation 89 5.3.2.2 Test Method 89 5.3.2.3 Experimental Results 90 5.3.3 Strength Modification due to size effect 95 5.3.3.1 Weibull Weakest-Link Theory 96 5.4 Analytically Derived Properties 98 5.4.1 Minimization Technique 99 5.4.2 Experimental Tests 101 5.4.2.1 Material Preparation and Test Method 101 5.4.2.2 Experimental Results 102 5.4.3 Minimization Results 104 5.4.3.1 Minimization using 2-Dimensional Model 104 5.4.3.2 Minimization using 3 Dimensional Model 108 Chapter 6 - Model Verification 112 6.1 Introduction 112 6.2 Compression Verification 113 6.2.1 Numerical Analyses 113 6.2.2 Experimental Tests 118 6.2.3 Results 118 6.3 Tension Verification 124 6.3.1 Numerical Analyses 125 6.3.2 Experimental Tests 127 6.3.3 Results 129 6.4 Bending Verification 135 6.4.1 Numerical Analyses 135 vi 6.4.2 Experimental Tests 144 6.4.3 Results • 144 6.5 Parallam* Investigation 153 6.5.1 Simulation of Parallam* Strand Lay-up 153 6.5.2 Parallam* Tension Comparison 155 6.5.2.1 Numerical Analyses 155 6.5.2.2 Experimental Tests 156 6.5.2.3 Results 157 6.5.3 Parallam* Bending Comparison 159 6.5.3.1 Numerical Analyses 159 6.5.3.2 Experimental Tests 161 6.5.3.3 Results 161 Chapter 7 - Conclusion 164 7.1 Summary 164 7.2 Concluding remarks 166 7.3 Future related research 167 References 169 Appendices Appendix A - Hardening Model 175 Appendix B - Parametric Study 178 Appendix C - Program Layout 181 vii L I S T O F F I G U R E S Chapter 1 ^ Figure 1.1 - Parallam*, PSL Manufacturing Process 3 Figure 1.2 - Principal Material Coordinates of PSL 4 Chapter 2 Figure 2.1 - Lamina and Laminate Coordinate System 7 Figure 2.2 - Theoretical strength envelopes for wood with different F 1 2 values 13 Figure 2.3 - Six original tests proposed by Tsai and Wu (1971) 14 Chapter 3 Figure 3.1 - Constitutive model for wood laminae 22 Figure 3.2 - Principal and Off-axis Coordinate Systems of a Lamina 26 Figure 3.3 - Subsequent Yield Surfaces for Hardening Material 34 Figure 3.4 - Illustration of Normality of Plastic Flow Vector for Associated Plasticity . . . . 35 Chapter 4 Figure 4.1 - Two-dimensional isoparametric element: (a) rectangular element; (b) distorted element 43 Figure 4.2 - Geometry of N - Layered Laminate Qones, 1975) 47 Figure 4.3 - Linear Solid Isoparametric Element 50 Figure 4.4 - Initial Stiffness Solution Algorithm for a Single Variable Situation (Owen and Hinton, 1980) 56 Figure 4.5 - Incremental Stress Changes at Initial Yield 60 Figure 4.6 - Stepped Process to Reduce Stress Point to Yield Surface 62 Figure 4.7 - Cantilever Beam Example 65 Figure 4.8 - Load versus Displacement for Isotropic Cantilever Beam Example 67 Figure 4.9 - Finite Element Mesh for 2 and 3 Dimensional Analysis of Isotropic Cantilever Beam 68 Figure 4.10 - 5 Point Continuous Beam 68 viii Figure 4.11 - Finite Element Mesh for Continuous Beam Example: a) 2 dimensions b) 3 dimensions 69 Figure 4.12- 3 Dimensional ANSYS* Solution to Continuous Beam Example 71 Figure 4.13- Exploded View of Symmetric Angle-ply Laminate 72 Figure 4.14 - Finite Element Mesh for Angle Ply Laminates a) 2 dimensions b) 3 dimensions 73 Chapter 5 Figure 5.1- Compression Test Set-up 78 Figure 5.2 - Typical Stress-Strain Curve for Compression Parallel-to-Grain Specimens . . 79 Figure 5.3 - Typical Stress-Strain Curve for Compression Perpendicular-to-Grain Specimens 80 Figure 5.4 - Tri-linear Approximation of Compression Stress-Strain Curve 81 Figure 5.5 - Compression Parallel-to-Grain Elastic Modulus Data 85 Figure 5.6 - Compression Parallel-to-Grain Tangent Modulus Data 86 Figure 5.7 - Compression Parallel-to-Grain Yield Stress Data 86 Figure 5.8 - Compression Parallel-to-Grain Ultimate Stress Data 87 Figure 5.9 - Compression Perpendicular-to-Grain Elastic Modulus Data 87 Figure 5.10 - Compression Perpendicular-to-Grain Tangent Modulus Data 88 Figure 5.11 - Compression Perpendicular-to-Grain Yield Stress Data 88 Figure 5.12 - Compression Perpendicular-to-Grain Ultimate Stress Data 89 Figure 5.13 - Tension Test Set-up 90 Figure 5.14 - Typical Parallel-to-Grain Tensile Behaviour 92 Figure 5.15 - Typical Perpendicular-to-Grain Tensile Behaviour 92 Figure 5.16 - Lognormal fit of Parallel-to-Grain Tensile Strength Data 93 Figure 5.17 - Lognormal Fit of Parallel-to-Grain Tensile Modulus Data 93 Figure 5.18 - Lognormal Fit of Perpendicular-to-Grain Tensile Strength Data 94 Figure 5.19 - Lognormal Fit to Perpendicular-to-Grain Tensile Modulus Data 94 Figure 5.20 - Select Sample of Stress-Strain Diagrams for [± 15]s Angle-ply Laminates . . . 103 Figure 5.21 - [ + 15]s Angle-ply Laminate Compression Specimen 103 ix Figure 5.22 - Cumulative Probability Distribution of [+ 15]s Laminate in Compression (using 2-D Finite Element Model) 107 Figure 5.23 - Stress-Strain Curves of [± 151 Angle-ply Laminate in Compression (using 2-D Finite Element Model) 108 Figure 5.24 - Cumulative Probability Distribution of [+ 15]s Laminate in Compression (using 3-D Finite Element Model) 110 Figure 5.25 - Stress-Strain Curves of [+ 15]s Angle-ply Laminate in Compression (using 3-D Finite Element Model) I l l Chapter 6 Figure 6.1 - Finite Element Mesh for [±30] s Angle Ply Laminates in Compression: a) 2 dimensions b) 3 dimensions 113 Figure 6.2 - Explanatory Depiction of 2 Dimensional Element 114 Figure 6.3 - Input Parameters for [±301 Angle Ply Laminate in Compression (2 Dimensional Program) 116 Figure 6.4 - Cumulative Probability Distribution of [ + 30]s Laminate in Compression (2-d Model) 119 Figure 6.5 - Cumulative Probability Distribution of [±30] s Laminate in Compression (3-d Model) 120 Figure 6.6 - Location of Sample Points Used to Show Stress Redistribution: a) 2 Dimensions b) 3 Dimensions 121 Figure 6.7 - Stress - Strain Diagram Depicting Stress Redistribution Upon Stress Point Failure (2 Dimensions) 122 Figure 6.8 - Stress - Strain Diagram Depicting Stress Redistribution Upon Stress Point Failure (3 Dimensions) 122 Figure 6.9 - Stress-Strain Curves of [±30] s Angle-ply Laminate in Compression (using 2 - Dimensional Model) 123 Figure 6.10 - Stress-Strain Curves of [±30] s Angle-ply Laminate in Compression (using 3 -Dimensional Model) 124 Figure 6.11- Tension Test Set-up 129 Figure 6.12- Tensile Failure of [ ± 30]s Angle Ply Laminate Specimen 129 x Figure 6.13 - Cumulative Probability Distribution of [+ 15]s Laminate in Tension (2 - Dimensional Model) 132 Figure 6.14 - Cumulative Probability Distribution of [± 15]s Laminate in Tension (3 - Dimensional Model) 132 Figure 6.15 - Cumulative Probability Distribution of [±30] s Laminate in Tension (2 - Dimensional Model) 133 Figure 6.16 - Cumulative Probability Distribution of [±30] s Laminate in Tension (3 - Dimensional Model) 133 Figure 6.17 - Stress-Strain Curve Comparison of [± 151 and [±30] s Laminates in Tension . 134 Figure 6.18 - Finite Element Mesh for Bending Configuration a) 2 Dimensions b) 3 Dimensions 136 Figure 6.19 - Comparison of Coarse and Fine Mesh for Bending Analysis 137 Figure 6.20 - Exploded View of Origin of Simulated Load Displacement Curves for [ +151 Angle-Ply Laminates 143 Figure 6.21 - 3 Point Bending Test Setup 144 Figure 6.22 - Cumulative Probability Distribution of Academic Laminates in 3 Point Bending (2 - Dimensional Model) 148 Figure 6.23 - Cumulative Probability Distribution of Academic Laminates in 3 Point Bending (3 - Dimensional Model) 148 Figure 6.24 - Load-Displacement Curves of [0°] Laminate in Bending (using 2 - Dimensional Model) 149 Figure 6.25 - Load-Displacement Curves of [0°] Laminate in Bending (using 3 - Dimensional Model) 149 Figure 6.26 - Load-Displacement Curves of [± 15]s Angle Ply Laminate in Bending (using 2 - Dimensional Model) 150 Figure 6.27 - Load-Displacement Curves of [± 15]s Angle Ply Laminate in Bending (using 3 - Dimensional Model) 150 Figure 6.28 - Load-Displacement Curves of [±30] s Angle Ply Laminate in Bending (using 2 - Dimensional Model) 151 Figure 6.29 - Load-Displacement Curves of [±30] s Angle Ply Laminate in Bending (using 3 - Dimensional Model) 151 xi Figure 6.30 - Load-Displacement Curves of [90°] Laminate in Bending (using 2 - Dimensional Model) 152 Figure 6.31 - Load-Displacement Curves of [90°] Laminate in Bending (using 3 - Dimensional Model) 152 Figure 6.32 - Cutting and Grain Angle Mapping Template for Parallam* 154 Figure 6.33 - Frequency Distribution of Measured Parallam* Grain Angle 155 Figure 6.34 - Finite Element Mesh for Parallam* in Tension 156 Figure 6.35 - ParallarrT Tensile Specimen Dimensions 157 Figure 6.36 - Tensile Failure of Parallam* Specimen 157 Figure 6.37 - Cumulative Probability Distribution of Parallam* in Tension (3 - Dimensional Model) 159 Figure 6.38 - Finite Element Mesh for Parallam* in 3 Point Bending 160 Figure 6.39 - 3 Point Bending Test Setup for Parallam* 161 Figure 6.40 - Cumulative Probability Distribution of Parallam* in 3 Point Bending (3 - Dimensional Model) 163 Figure 6.41 - Load-Displacement Curves of Parallam* in 3 Point Bending (using 3 - Dimensional Model) 163 Figure A. 1- Bilinear Stress - Plastic Strain Curve 176 Figure B . l - Various Mesh Grades Investigated for Uniaxial Analyses: a) 3 x 3 grid, b) 4 x 4 grid, c) 5 x 5 grid 178 Figure B.2 - Stress - Strain Diagrams of 3 Investigated Grids 180 xii L I S T O F T A B L E S Table 4.1 - Criteria for Tension vs. Compression Dominance 59 Table 4.2 - Brittle Failure Stress Reduction Ratios 64 Table 4.3 - Linear Elastic Orthotropic Comparison of COMAP with ANSYS* 70 Table 4.4 - Laminate Stiffness and First Ply Failure Comparison of P C L A M and COMAP. 74 Table 5.1 - Compression Properties of Strands 81 Table 5.2 - Analysis of Variance Results for Heartwood and Sapwood for Parallel-to-grain Compression 82 Table 5.3 - Analysis of Variance Results for Heartwood and Sapwood for Perpendicular-to-grain Compression 83 Table 5.4 - Correlation Matrix for Compression Parallel-to-Grain Specimens 84 Table 5.5 - Lower Triangle of Correlation Matrix for Compression Parallel-to-Grain Specimens 84 Table 5.6 - Correlation Matrix for Compression Perpendicular-to-Grain Specimens . . . 84 Table 5.7 - Lower Triangle of Correlation Matrix for Compression Perpendicular-to-Grain Specimens 85 Table 5.8 - Tension Properties of Strands 91 Table 5.9 - Compression Properties of [± 15]s Angle-ply Laminate 102 Table 5.10- Result of Nonlinear Minimization using 2-dimensional F.E. Model 104 Table 5.11- Result of Nonlinear Minimization using 3-dimensional F.E. Model 109 Table 6.1- Categories of Model Verification Tests 112 Table 6.2 - Finite Element Mesh Properties of Angle Ply Laminate Uniaxial Analyses . . 114 Table 6.3 - Input Properties of 3 Dimensional Program which differ from that of 2 Dimensional Program 117 Table 6.4 - Experimental and Simulated Data for [ + 30]s Angle Ply Laminate in Compression 119 Table 6.5 - Input Properties for [+ 15]s and [±30] s Angle Ply Laminates in Tension . . . 128 Table 6.6 - Experimental and Simulated Data for [+ 15]s and [±301 Angle-Ply Laminates in Tension 130 xiii Table 6.7 - Finite Element Mesh Properties of Angle Ply Laminate Bending Analyses . . 136 Table 6.8 - Tensile Strength Properties Used for 3 Point Bending Analyses 142 Table 6.9 - Experimental and Simulated Data for [0°],[± 15]s, [±30] s and [90°] Laminates in 3 Point Bending 145 Table 6.10- Experimental and Simulated Data for Parallam* in Tension 158 Table 6.11- Tensile Input Properties for Parallam* Bending Analysis 160 Table 6.12 - Experimental and Simulated Data for Parallam* in 3 Point Bending 162 Table B . l - Input Properties for Parametric Study 179 xiv N O T A T I O N The following symbols are used in this study: [A] = direction cosine matrix; a; = defined in Equation 3.41; a;j = direction cosine components; [B] = strain displacement matrix; b = width of specimen (mm); b; = randomly generated bivariate standard normal parameters; [C] = 3 dimensional material stiffness matrix in principal material coordinates; [C] = 3 dimensional material stiffness matrix in principal material coordinates; C;j = components of 3 dimensional material stiffness tensor (MPa) ; d = depth of specimen (mm); E ; = elastic moduli (MPa); E;' = tangent moduli (MPa); F = loadfN); {F} = external force vector; F|j = components of strength tensor (N); f = yield function; G = in-plane shear modulus (MPa); g = plastic potential function; ET = hardening modulus (MPa) ; I = moment of Inertia (mm4); [ J] = Jacobian matrix; [K] = global stiffness matrix of structure; [k]e = element stiffness matrices; k = threshold stress (MPa); L = length (mm); L;J = components of lower triangle of correlation matrix ; M = bending moment (N-mm); n m N: = components of strength tensor ; scale parameter of Weibull distribution; shape functions; total number of layers; P = number of probability levels; {P} = internal resisting force vector; [Q] = 2 dimensional material stiffness matrix in laminate coordinate system [Q'] = 2 dimensional material stiffness matrix in principal material coordinate system {R} = resultant in-plane force vector; R = ratio to proportion stress value; S = in-plane shear strength (MPa); SD = standard deviation; [T] = transformation matrix; t = thickness of lamina (mm); {u} = vector of nodal displacements; u, v, w = displacement in x, y, z directions (mm); V = volume (mm3); W ; = weighting factors for Gaussian integration; W i n t , W e x t = internal and external work (MPa) ; W p = plastic work(MPa) ; X = parallel-to-grain strength (MPa); Xj = parallel-to-grain tensile strength (MPa); X c , X c u = parallel-to-grain compressive yield and ultimate strength, respectively (MPa) ; Y = perpendicular-to-grain strength (MPa) ; Y t = perpendicular-to-grain tensile strength (MPa) ; Y c , Y C U = perpendicular-to-grain compressive yield and ultimate strength, respectively; (MPa) z ; = standard normal random variable; <X; = parameters which define the offset of the yield surface; P = shape parameter of Weibull distribution; Y; = engineering shear strain components; = components of 2 dimensional material stiffness tensor (MPa); xvi A = incremental value; vertical displacement (mm); Ei = strain components; e;P = plastic strain components; i " p = effective plastic strain; C = curvilinear coordinate of isoparametric element; •n = curvilinear coordinate of isoparametric element; 0 = grain angle (deg); A = vector of unknowns in minimization procedure; X = plastic multiplier; M = mean value; u ; j = Poisson's ratio; 5 = curvilinear coordinate of isoparametric element; n = defined in Equation 3.53; p; = portion of effective stress (MPa) ; O; = stress components (MPa) ; O = effective stress (MPa) ; T = non-specific material strength; $ = minimization function; X = hardening parameter; {¥} = residual force vector; Subscripts and Superscripts x, y, z = laminate coordinate directions; 1,2,3 = principal material directions; c = designates center value; e = designates elastic quantity; ep = designates elastic-plastic quantity; k = designates layer number; max = designates maximum value; o = designates original value; designates plastic quantity; iteration number; designates ultimate value; designates value at incipient yield. A C K N O W L E D G M E N T S I would like to extend a warm thank you to my supervisory committee, Dr. Frank Lam, Dr. Helmut Prion and Dr. Ricardo Foschi for their gentle guidance and support throughout the project. Thanks also go to members of my examining committee, Dr. J. David Barrett, Dr. Gary Schajer, and Dr. Robert Leichti for reviewing my work and providing helpful commentary. Acknowledgment is made to my friends and colleagues in the department of Wood Science, UBC, as well as our timber engineering lab manager, Bob Myronuk, who helped to make my graduate school days pleasantly memorable and rewarding. My gratitude goes to my Mother and Stepfather for providing support and encouragement when most needed and for keeping me focused and true to my goal throughout these years. Finally, I wish to express my deepest gratitude to my best friend and husband, Alexander Schreyer, for his unwavering support, ever-insightful comments and always sage advice - both technical and personal. xix ere order in variety we see And where, though all things differ, all agree Alexander Pope 1 - Introduction 1.1 I N T R O D U C T I O N A new division of structural building materials is rapidly gaining recognition and acceptance in today's construction industry. Generically known as Structural Composite Lumber (SCL), this group consists of Laminated Veneer Lumber (LVL), Parallel Strand Lumber(PSL), Laminated Strand Lumber (LSL) and thick oriented strand board/rimboard. These new engineered products have important advantages over conventional lumber owing to their sophisticated manufacturing processes. They are typically made with thin wood veneers, strands or flakes which are arranged and bonded together using adhesive under controlled heat and pressure. Reconstitution of the wood in this manner disperses the natural defects throughout the material resulting in more consistent mechanical properties than that of conventional lumber. This high consistency leads to a more efficient utilization of the wood fiber resource. Consequently, structural composite lumber is replacing conventional lumber in traditional wood applications (ie. beams and columns) and, moreover, are starting to be used in applications typically dominated by steel or concrete (ie. long span commercial roof trusses and shell structures). In contrast to the current growth of structural composite lumber, development of these products was relatively slow, depending predominantly on expensive empirical-based research initiatives. Typically, the products were refined gradually through experimental manipulation of the material constituents. The process was costly and time consuming. In light of this, an alternative approach for the development of future wood composite products is highly desired, if not outright necessary. A logical alternative is the implementation of a computational model. Computer models can be used to estimate the effects of varying raw material characteristics on the final product's mechanical properties thereby reducing fabrication and testing costs for a new product. In the same manner, the model could be used to optimize or customize current products. For example, with an accurate working model, companies could gauge the effect of strategically placing lower quality or less expensive species in the board. Moreover, an accurate computer model could be an invaluable design aid for analyzing (or optimizing for) complex applications, such as connection details. The question then becomes one of developing an accurate working model. Unfortunately, little work has been done in this area for wood-based composites. Only very recently have some studies attempted to analytically model relatively simple tensile specimens of strand or veneer-based composites. However, a vast amount of research has been done on the mechanical behaviour of advanced composite materials, (ie. those made from high modulus fibers such as graphite, silicon carbide, or boron, and used in aerospace applications). A multitude of modeling techniques and concepts have been established to characterize their stress-strain and failure behaviour, as discussed in any standard text on mechanics of composites (Jones, 1975; Hull, 1981; Gibson, 1993). One can fortunately draw ideas from this resource, while incorporating established wood mechanics principles. Such is the intent of this thesis. Specifically, the focus of this thesis is the development of a viable model for predicting the mechanical response of a parallel-aligned wood strand composite under bending. This material is necessarily a simplified version, with convenient and controllable geometry, of the more complex material parallel strand lumber (PSL). Although beyond the scope 2 of this thesis, the ultimate aim would be to provide a tool for commercial wood composites. As such, the current characteristics of the commercial product Parallam*, PSL have been used as a general reference from which to determine the model's primary considerations. 1.2 M O D E L C O N S I D E R A T I O N S Parallam*, PSL is a complex wood composite material fabricated through a technologically advanced processing technique. A visual depiction of the manufacturing process, taken from Sharp (1996), is shown in Figure 1.1. In forming a billet, the strands, ranging in length from approximately 0.6m to 2.4m, are delivered to a moving trough where they are randomly arranged through the depth and width and roughly aligned along the length. This deposition system results in a product with a highly complex geometric layout. From a macroscopic standpoint, the product could be considered as homogeneous and product behaviour characterized by averaged mechanical properties. I. Veneers peeled 2. Clipped.into strands .3. Adhesive! applied 4. Aligned & pressed 3. Parallam* PSL Figure 1.1 - Parallam*, PSL Manufacturing Process However, this approach has no regard for the constituent materials (ie. the strands and glue). It lacks the ability to account for varying properties or geometric arrangement of the strands in order to customize composite performance. Ergo, the present model's methodology is based on the mechanical properties of the individual strands and the geometric lay up is implicit in the analysis. The strand properties are assumed to be homogeneous within each strand and are established 3 empirically. It is noted that although gross longitudinal misalignment of strands is rare in PSL, gross grain deviation is considered in this thesis because individual strands may contain them, and future product development may investigate variation of fibre orientation to improve transverse strengths and stiffnesses. The strand deposition system also results in a product with high strength and stiffness in the longitudinal direction and comparatively poor properties in other directions. Macroscopically, wood is assumed to be orthotropic whereby properties differ in three orthogonal directions only, (Bodig and Jayne, 1982; Perkins, 1967) which greatly simplifies the stress-strain relationships for the material. Similarly, it is assumed that PSL is orthotropic with principal material coordinates as shown in Figure 1.2. The mechanical properties of the strands are inherently variable. This variability can significantly influence the accuracy of the strength prediction model and should be included integrally. In this thesis, parameter variation is handled through a probabilistic approach whereby measured quantities Transverse (y) axis Through-thickness (z) axis Figure 1.2 - Principal Material Coordinates of PSL 4 are expressed in statistical terms, known as 'random variables'. Probabilistic solutions are obtained through Monte Carlo simulations. When using data obtained from small specimens to predict the strength of large members, one must consider the variation of strength with size. This 'size effect' is most prominent with brittle modes of failure and is dependent on the strength variability of the material. Size effect has been acknowledged and addressed in many wood and timber studies such as, Barrett et al. (1975), Madsen and Buchanan (1986), and Barrett et al. (1995), as well as studies on structural composite lumber by Sharp and Suddarth (1991), and Clouston et al. (1998). Each of these studies espoused the use of Weibull weakest link theory to quantify size effect. As such, this theory is adopted for use in the current study. l .3 O B J E C T I V E A N D S C O P E The fundamental objective of the thesis is to develop a computational model for simulating the bending behaviour of parallel-aligned wood strand composites. An orthotropic finite element code will be developed for the general nonlinear analysis of wood laminates with layers of varying fiber orientation. The code will be developed within a probabilistic framework using random variables as input enabling stochastic-based Monte Carlo simulation analyses. It will be shown that orthotropic plasticity theory is effective in describing the nonlinear behaviour up to and beyond the point of laminate failure. Chapter 2 will provide a general review from the literature on the background of the yield criterion chosen for use in the constitutive model, as well as an outline of published studies on modeling of 5 wood composites. Chapter 3 aims at a complete description of the elastic-plastic-failure constitutive model employed in the model. Chapter 4 provides a review of the finite element formulation as well as a thorough program performance verification. Chapter 5 outlines the acquisition of the strand database while chapter 6 details application of the model to a series of problems. The numerical simulations are compared to experimental data for validation. Finally, in chapter 7, a summary is provided as well as a discussion on further areas of research. 6 2 - Background 2.1 I N T R O D U C T I O N Despite its importance, only a small number of studies have focused on strength prediction of structural wood composites. These few studies, all the same, employed conventional strength analysis techniques which have been proven successful for laminated composite materials. The basic building block in the analysis is the mechanical behavior of the individual layer, or lamina, (Figure 2.1) which can be studied from either a micro-mechanical or macro-mechanical perspective. Ay y Figure 2.1 - Lamina and Laminate Coordinate System Micro-mechanical theories determine the mechanical properties of a lamina by addressing constitutive information about the fiber and the matrix phase. The interaction of the phases is examined in detail through mechanics of materials principles - for example, the rule of mixtures -to determine strength and stiffness properties of the heterogeneous material. Macro-mechanical theories, on the other hand, involve the use of gross effective lamina properties wherein the material is presumed homogeneous. The effects of the constituent materials as well as 7 inherent cracks, notches or other discontinuities which lead to failure are not directly modeled; rather, they are incorporated indirectly in an averaged or "smeared" manner. A macro-mechanical approach has been adopted for this thesis. Lamina properties are established empirically and used in a semiempirical failure criterion. This chapter presents first an introduction to orthotropic failure criteria, which are integral in a macro-mechanical analysis, followed by a review of the specific criterion used for the present model, and finally a survey of wood composite publications focused on material behaviour modeling. 2.2 F A I L U R E T H E O R I E S Failure theories are, in general, mathematical models which may be determined theoretically by rational modeling of the material's physical characteristics or empirically by simply representing experimental observations. In the latter case, the theories are termed phenomenological 'treating the heterogeneous material as a continuum' and making no attempt to explain the mechanisms which lead to material failure (Wu, 1974). Early classical theories were based on limiting a physical variable such as stress, strain, or strain energy assuming a homogeneous and isotropic material. For example, one of the first theories, the Maximum Normal Stress theory, maintains that material failure occurs when the maximum normal stress reaches a critical value independent of other stresses at that point. The critical value can be determined from a uniaxial tensile test. Three simple, independent equations result: CTi = ° c ; CT2 = < * c ; ° 3 = ° c " (2-i) 8 where; O ; (i = 1, 2, 3) denote principal stresses and o c is the critical stress. These equations may be plotted with each normal stress component constituting an axis forming a triad in the stress space. This is known as a failure envelope (or surface). Combinations of stresses contained inside the failure envelope signify survival and on or beyond the envelope indicate material failure. The concept is common for both isotropic and orthotropic criteria, however, for the latter the plotted stresses usually correspond to those along the principal material axes. Many orthotropic criteria are based on isotropic criteria. The first orthotropic theory, devised by Hill in 1948, extended the well known isotropic von Mises theory for use with specially orthotropic materials (ie. materials for which the direction of the applied stresses coincide with the principal material directions). Hi l l formulated the criterion such that if the orthotropy vanished, the theory reduced to von Mises theory. The von Mises theory was also used by C.B. Norris in 1950 to develop a strength criterion for orthotropic materials.1 His theory applied the isotropic failure criterion to a simplified geometrical model. By his own admission, Norris's approach was not 'rigorously correct', representing an orthotropic material by an isotropic material with regularly spaced voids. Nonetheless, his equation produced good results when compared to data for 3 ply plywood and fiberglass laminate for which most of the glass fibers were aligned in one direction (Norris, 1962). Norris's theory has been recommended for design of glued-laminated beams in the 1974 Timber Construction Manual. However, it was also criticized for being overly conservative for many beam and arch situations in This theory succeeded his so called 'interaction formula' suggested in 1945 for use with plywood plates. 9 practice (Kobetz and Krueger, 1976). There are numerous orthotropic yield theories available in the literature as compiled and compared in surveys by researchers Sandhu (1972), Rowlands (1985) and Nahas (1986). Each criterion has strengths and limitations and no one criterion is suitable for all materials. However, one criterion, made popular by Tsai and Wu in 1971, received widespread attention due to its simplicity and generality. The two-dimensional (i.e. plane stress) form of this theory will be used in this study. 2.2.1 Tsai-Wu criterion Tsai and Wu proposed a simplified version of a tensor polynomial. Their theory predicts that failure will occur when the following inequality is satisfied: F -o - i+fV 7 .^ 1 ' (2-2) where; i , j = 1, 2,... 6 (repeated indices imply summation), and F^nd F ; j are second and fourth rank strength tensors, respectively. Under plane stress, o3, o4, and o 5 are assumed negligible and the F 6, F 1 6 and F 2 6 terms are zero due to special orthotropy. In expanded form Eq. (2.2) becomes: F ^ + F 2 a 2 + Fnax2 + F 2 2 a 2 2 + 2Fna1a2 + F 6 6 a 6 2 > 1 (2.3) The coefficients F i through F 6 6 , with the exception of F 1 2 , are described in terms of the strengths in the principal material directions. For a lamina in plane stress, these are: longitudinal strength in both tension and compression ( X T , X j , transverse strength in both tension and compression (Y x, and in-plane shear strength (S). Considering a uniaxial tension load on a specimen in the 1 10 direction, the above equation at failure becomes: F i X T + F n X , - 2 = l (2.4) and for compression is: F . X c + F ^ X c ^ l (2.5) By solving the equations 2.4 and 2.5 simultaneously and regarding the compression strength as negative, the expression for the strength parameters F t and F n are found to be: 1 1 1 F i = ; Fu = (2.6) Through similar mathematical manipulations, it can be shown that: _ J \_ 1 _ J _ A drawback in the application of the tensor polynomial theory is that there is no consensus among researchers for a method of determining the value of the combined stress parameter, F 1 2 which accounts for the interaction between normal stresses, ol and o 2. The only certainty is that in order for the failure surface to be a closed ellipsoid2, the coefficient F 1 2 must be bounded by the stability condition: F n F 2 2 - F 1 2 2 >0 (2.8) 2. Plotting equation 2.3 with c^ and o2 as coordinates of a point will produce an ellipsoidal failure surface, providing F 1 2 is within the prescribed range. Else, the surface would become open-ended, meaning no matter how large the stresses became, failure would never ensue, which is physically impossible. 11 The 'correct' determination of F 1 2 has been the topic of discussion for many researchers for many years. 2.2.1.1 Significance of the Interaction Term, F ] 2 The interaction term, F 1 2 characterizes the interaction of the normal stresses, o± and o2. Because this term involves both normal stresses, its evaluation must occur under a biaxial loading situation unlike the other strength tensors whose values are determined from uniaxial loading. The Tsai-Wu failure envelope takes the form of an ellipse in the Oj ,o2plane as shown in Figure 2.2. The strength tensors F 1 ( F 2 , F n and F 2 2 establish the axes intercepts whereas the component F 1 2 determines the rotation of the ellipsoid with respect to the principal material axes. Figure 2.2 shows the influence of varying F 1 2 in equation 2.2 when no shear stress is present. It has been argued that the value of F 1 2 "determines the effectiveness of tensorial-type failure criteria." (Suhling et al., 1984) 12 -30-1-Figure 2.2 - Theoretical Strength Envelopes for Wood with Different F 1 2 Values 2.2A.2 Evaluation of interaction term, F 1 2 Tsai and Wu (1971) reported that a separate combined stress test was necessary for determining F 1 2 . They provided six viable test options in their original paper: a simple biaxial tension or compression test with the applied load equal in each direction; a 45 degree off-axis tension or compression test; or a 45 degree off-axis positive or negative shear test as illustrated in Figure 2.3. 13 A y Biaxial tension Biaxial compression Off-axis tension Off-axis compression Positive shear Negative shear Figure 2.3 - Six Original Tests Proposed by Tsai and Wu (1971) Strict care must be taken in determining F 1 2 due to its high sensitivity to experimental variation. Tsai and Wu showed that for a graphite epoxy composite, slight inaccuracies in measurements of strength for most of the above tests would result in large inaccuracies in the calculated value of F 1 2 . This problem is accentuated by the fact that in order for F 1 2 to be physically admissible, it must satisfy the prescribed stability bounds (Eqn. 2.8). Many researchers conducted similar studies on a variety of materials in an attempt to determine F 1 2 experimentally. Pipes and Cole (1973) performed a limited number of off-axis tensile tests on boron epoxy composites to determine F 1 2 . They tested two coupons at 60 degrees, three at 45 degrees, four at 30 degrees, and three at 15 degrees. Using the mean of each set they found that only the value obtained for the 15 degree off-axis tests satisfied the stability criterion. Consequently, they concluded that the off-axis tensile test was not an adequate method of determining F 1 2 for boron epoxy composites. Suhling, et al. (1984) conducted a comprehensive study to establish F 1 2 for paperboard employing « 14 all of the six test methods suggested by Tsai and Wu (1971). They found the optimum value of F 1 2 by fitting a linear regression through experimental results for all four quadrants and compared this value with that obtained from the individual tests. The tension biaxial tests produced the closest value to the optimum value for a zero shear situation. They found that F 1 2 = 0 was a satisfactory solution for four shear levels considered; o 6 = 0, 6.9, 10.3, and 15.9 MPa. They also concluded that due to the highly sensitive and unstable nature of F 1 2 when calculated using off-axis tests that this test was not a suitable method to determine the interaction parameter. Anticipating the problems that researchers might encounter when determining F 1 2 , Wu published a paper in 1974 demonstrating a method to find an optimal biaxial stress ratio for calculating F 1 2 experimentally. This method was not strongly espoused because it involved a series of experimental iterations further complicating the tensor polynomial theory. Difficulties in experimental determination of F 1 2 prompted several researchers to seek a theoretical solution to the problem. Narayanaswami and Adelman (1977) performed a numerical analysis to prove that, from a practical standpoint, the arbitrary assignment F 1 2 = 0 was acceptable for filamentary composites. For the six tests described above, they computed the percentage error when setting F 1 2 = 0 for 10 different composite materials. In all cases the error was found to be less than ten percent and they therefore concluded this to be an acceptable error for practical engineering applications. Cowin (1979) derived Hankinson's formula: X Y X s i n 2 9 + Y c o s 2 9 15 (an empirical equation which has had remarkable and repeated success in predicting wood strength, o e, at an angle to grain, 0) using only the linear term of the tensor polynomial function, F ; o;. By considering the quadratic terms as well, he derived a formula for F 1 2 : Fl2 = V F H F 2 2 V (2-10) 2S Z that produced a similar, yet more accurate, equation to Hankinson's criterion to represent bone strength with respect to angle to grain. Both van der Put (1982) and Liu (1984) derived an expression for F 1 2 that reduces the tensor polynomial theory into the Hankinson formula: F 1 2 -i f 1 + 1 _ j J 2 \X-j-Y(- X(-Y-p S / (2.11) van der Put approached the problem by transforming the applied stresses to the principal material coordinate system while Liu transformed the strength tensors. Despite all efforts, no standard method to determine F 1 2 has been established. As such, a probabilistic approach has been developed in this study and will be addressed in chapter 5. 16 2.3 L I T E R A T U R E R E V I E W O F W O O D C O M P O S I T E M O D E L I N G One of the earliest attempts to analytically predict wood composite behaviour was by Hunt and Suddarth (1974). They predicted tensile modulus of elasticity and shear modulus of rigidity of medium density homogeneous flakeboard. The method involved a linear elastic finite element analysis of a tensile specimen modeled as an assemblage of four noded plate elements (representing the flakes) embedded in a rectangular grid of rigidly connected beam elements (representing the resin). The grain orientation in each plate element was randomly assigned simulating a random flake deposition in the mat. Young's modulus was calculated using Hooke's law where the stresses are simply the applied stress over the area and strain is the analog's average axial strain, Poisson's ratio was calculated from the usual relation v = e1 / e2 and shear modulus was determined from the isotropic relation G = E / 2(l + v). A Monte Carlo technique was used to produce distributions of the desired values (E and G) and the average analytical values were compared with experimentally obtained averages for two separate species, Douglas-fir and aspen. The model underestimated the experimental tensile modulus by 8 percent for aspen and 6 percent for Douglas-fir whereas the shear modulus was overestimated by 10 percent and 13 percent for aspen and Douglas-fir respectively. More recently, Shaler and Blankenhorn (1990) evaluated the suitability of two composite theories -the rule of mixtures and the Halpin-Tsai equation - to predict flexural modulus of elasticity of oriented flakeboard. They considerd flake geometry and orientation, density, resin content and species. The modeling approach was a two-stage process. The first stage treated the resin as matrix and the wood as fiber and the second stage treated the resin and wood as matrix and air as fiber. The method underestimated measured moduli by an average of 25 percent. 17 In 1993, Triche and Hunt developed a linear elastic finite element model capable of predicting the tensile strength and stiffness of a PWSC with controlled geometry. The model was micromechanical in nature considering each strand to have three layers ({e. pure wood, resin and an interface between them) and used, as input, the properties of the individual constituents. For' efficiency, the concept of a 'superelement', (an element representing one entire strand system) was introduced. Superelements were built by forming a mesh of elements for the model of one strand and "statically condensing internal degrees of freedom so that the remaining model (contained) only degrees of freedom that (were) essential in building the board model." Composite boards were constructed from a group of superelements and corresponding board stiffness (using the stiffness matrices of the superelements) was computed in a conventional finite element procedure. A large experimental database on the longitudinal tensile stiffness (and to a lesser extent, strength) of yellow poplar strands was developed. Also, a smaller sample of laminated assemblies was developed to assess the thickness change due to pressing of the strands and to acquire the properties of the wood-resin interface. The model was enhanced to consider the effect of overlapping strands. This was handled by "building on the existing model until its predicted E matched the ..." mean value for a sample of lap sandwich specimens of the average E within the specimen. The veracity of the model was checked through comparison of predicted and experimental results of 14 distinct categories of small composite boards with controlled geometry. Maximum strength was predicted according to several failure criteria including maximum stress theory and Tsai-Wu theory and board elasticity was calculated based on the nodal displacements of the model. Excellent accuracy was reported for the predicted E (from 0.0 to 11.1% error); however, prediction of 18 maximum stress was inconsistent and, in at least one case, unacceptable (from 1.2 to 101.1% error). Cha and Pearson (1994), developed a two-dimensional finite-element model which could predict the elastic tensile properties of a 3 ply veneer laminate consisting of an off-axis core ply of varying angles with or without a crack and clear, straight-grained outer plys. The laminate was modeled as a single plate and the stress components were defined across the laminate thickness. Good agreement was obtained between predicted and experimental strains at maximum load (maximum difference of 14.3%) as well as predicted and experimental stresses (maximum difference of 7.7%). Further, the model was used to investigate the influence of cracks and grain angle of the middle ply. For this, stresses within the layers were estimated as the product of the effective strain from the model and the corresponding elastic modulus of the lamina, determined experimentally. Resulting transverse and shear stresses were charted for varying grain angles and specified crack lengths. Recently, a very progressive study was carried out by Wang and Lam (1998) in which a three dimensional, nonlinear, stochastic finite element model to estimate the probabilistic distribution of PWSC tensile strength was developed. The model is based fundamentally upon longitudinal tensile strength and stiffness data of individual strands. A single strength modification factor (a) was incorporated, defined as a function of individual strength factors (aL, aT, a i n t, aoth) as well as the number of plys in the composite. The individual strength factors accounted for, respectively, 1) Length and 2) Thickness effect (defined by Weibull weakest-link theory), 3) the presence of resin and its influence at the interface, and 4) the effect of other factors, estimated by way of a finite element analysis together with nonlinear, least square regression with experimental data for two and three ply assemblies. 19 The finite element model considered strictly longitudinally laminated strands employing a three dimensional layered element with 8 nodes. Load was applied incrementally and a simple, uniaxial failure criterion (a function of the strength modification factor) was used to establish individual element failure. Upon detection of failure, a nonlinear algorithm was implemented to calculate the subsequent displacements and hence strains and stresses depending on the modified structural stiffness matrix. The model was verified through comparison with experimental data for four and six ply laminates. Further, various finite element grids were used, (affecting the size factor) to ensure the robustness of the prediction model. It was reported that, in all cases, excellent agreement was found. 20 3 - Model Development 3.1 INTRODUCTION It is of interest to note that in the past, research and research guidelines (e.g. ASTM D143 and D198) for wood and timber bending strength have focussed to a large degree on linear elastic constitutive theory. This is logical in that, as first indicated by Perkins (1967), the bending strength of the conventional building material, timber, is generally governed by strength reducing flaws such as knots or slope of grain contributing to a brittle, tensile failure, which is accepted as being linear elastic (Madsen and Buchanan 1986; Barrett et al. 1995). However, the stress-strain behaviour of wood in compression, both parallel and perpendicular to grain, is known to be nonlinear (Goodman and Bodig, 1971; Maghsood et. al., 1973; Conners, 1988). Consequently, bending fracture in high grade timber (with few defects) is considered to be quasi-brittle, where the load-displacement behaviour may be nonlinear but ultimate failure is brought on by brittle fracture. It is speculated that due to the minimized and controlled dispersement of flaws in wood composites such as PSL, compressive stresses play a more prominent role in a bending analysis when compared to that for timber. Hence, restriction of a bending analysis for a wood composite to the elastic range would result in an inefficient design because it would not account for the material strength beyond the proportional limit. In light of this, the constitutive model for the wood laminae of this study is comprised of four basic behavioural domains: elastic, elasto-plastic, post-failure brittle, or post-failure ductile, as illustrated by a simple uniaxial stress-strain curve in Figure 3.1. 21 o elastoplastic elastic post-failure ~ ductile brittle e Figure 3.1 - Constitutive Model for Wood Laminae It is presumed therein, that damage which occurs beyond the proportional limit can be described within the framework of classical plasticity theory. Upon initial loading, the material is assumed to behave linearly-elastic. Beyond the elastic domain the material may either 1) fail in a brittle manner or 2) first strain-harden and then ultimately fail in either a ductile or a brittle mode. For ductile failure, the material is assumed to lose all stiffness but retains strength, whereas for brittle failure both stiffness and strength are lost. Strain-hardening is characterized by successive growth of the yield surface defined by the Tsai-Wu criterion. The objective of this chapter is to derive the constitutive relationships for each behavioural regime. 3.2 C O N S T I T U T I V E M O D E L 3.2.1 Linear elastic regime The fundamental constitutive equation to describe linear elastic material behaviour is the well known Hooke's law, written here in generalized tensorial form (Sokolnikoff, 1956), where repeated indices imply summation. This law defines a one-to-one mathematical relationship between the nine components of the stress tensor, a- and the nine components of the strain tensor, ijkiski ' i , j ,k , l = 1,2,3 (3.1) 22 €;J. The stiffness tensor, C;jkl, is of 4 t h order, consisting of 81 stiffness coefficients. From equilibrium requirements, both stresses and strains are symmetric (i.e., o- = a- and e;j = e ^ , so that there are truly only six independent components of each. This means that the stiffness tensor is also symmetric with respect to these indices (i.e. C i j k l = C j ; k l = C i j I k = C j i l k) resulting in 36 stiffness coefficients. For the convenience of general mathematical manipulation, the double indexed system is often replaced by a single indexed system having a range of i= 1,2,....6 as follows: = a, S , i = E l CT22 = ° 2 8 22 = s 2 e>33 = a 3 8 33 = s 3 o 1 2 z . c 1 2 = Yl2 °~13 = °~5 2 s 1 3 = Y l 3 = s 5 a 2 3 = CT6 2 e 2 3 = Y23 = £ 6 (3.2) where Y23) Yn> a n d Y12 denote standard engineering shear strains, which physically represent the distortional change in angle from 90°. Equation (3.1) can now be written in contracted notation: 0 ; = ^ ; i,j = l,2,...6 (3.3) It can be shown that the total possible number of independent stiffness coefficients may be reduced further from 36 to 21, independent of material symmetry. Given the strain energy density function, W, such that the stresses are defined as: 5W 23 we differentiate with respect to e:, to find d 2 W _ c de.de. 11 If we reverse the order of differentiation, we have a2w _ c de-de{ ]> (3.5) (3.6) As the order of differentiation on W has no influence on the result, Qj = Cj ; and the stiffness matrix is symmetric, reducing the total number of independent coefficients to 21. Wood, in general, is assumed to have three orthogonal planes of symmetry. When the stress-strain relationships are developed within this principal coordinate system, there is no interaction between normal stresses and shearing strains, nor between shearing stresses and normal strains or shearing strains. As such, the number of independent components of the stiffness tensor reduces from 21 to 9 as follows: °1 c n C12 Cl3 0 0 0 c 1 2 C22 ^ 2 3 0 0 0 ° 3 C23 C 3 3 0 0 0 0-4 0 0 0 0 0 CT5 0 0 0 0 2 C 5 5 0 ° 6 . 0 0 0 0 0 2C 66 Y12/2 Y13/2 Y23/2J (3.7) The stiffness coefficients, in terms of engineering constants (i.e. Young's moduli (E), shear moduli (G), and Poisson's ratio (v)) are: 24 _ 1 - u 2 3 o 3 2 ^ 1 E , E 3 A r _ 1 ~ ^ 13^31 E , E 3 A C 4 4 = G 1 2 r • — r = ^23 E 2 E , A E , E 3 A C 5 5 - G 1 3 C , 3 = r = ^33 _ U 31 + U 2 1 U 3 2 E , E 3 A l - u 1 2 u 2 1 (3.8) E , E , A ^66 = ^23 ^here A = 1 ^12^21 ^23^32 ^31^13 2 U 2 1 U 32^ '13 E j E 2 E 3 (3.9) Poisson's ratios are dependent and can be estimated through the reciprocal relationship 1,1=1,2,3 (3.10) Lamina analysis often assumes a two-dimensional stress-state (plane stress) where o 3, o 5, and o 6 = 0. Considering this, the lamina stiffness matrix, denoted Q ; j, contains only four independent components as follows: -1 1 - Ui,U 0 12 u 21 1 2 E 2 1- o 1 2 o 2 1 0 U 1 2 E 2 1 - u 1 2 u 2 1 0 0 0 2G 12 E l £ 2 ' = [Ql £ 2 J 12 /2 712/2 (3-11) 3.2.1.1 Constitutive Relations for Off-axis Coordinates If the stresses and strains are defined in some non-principal material directions, (or off-axis coordinates) such as x and y in Figure 2.1, as is often the case for a lamina within a laminated member, the lamina is termed 'generally' orthotropic. The stiffness matrix in this case is fully 25 populated (ie. no zero terms present) and can be derived in terms of Qj by considering the transformation laws for Cartesian tensors as follows: If we consider two rectangular coordinate Figure 3.2 - Principal and Off-axis systems with a common origin as shown in Coordinate Systems of a Lamina Figure 3.2, we can calculate the relative orientation of one axis with respect to the other through use of the direction cosines a:i = cos (x;, x'). In the case of the lamina in Figure 3.2, a l i = A = a l 1 a 1 2 a 1 3 a 2 1 a 2 2 a 2 3 a 3 1 a 3 2 a 3 3 cosG sinG 0 - s inG cosG 0 0 0 1 (3.12) Referencing Figure 2.1, we can now transform the stress (or strain) tensor from the laminate axes (x,y,z) to the principal axes (1,2,3), using the direction cosine matrix in accordance with the transformation law (Mase, 1970), ° 1 1 ° 1 2 CT13 o 1 2 o" 2 2 a 23 _ ° 1 3 ° 2 3 ° 3 3 Simplifying the right hand side and using contracted notation, we have [A] 0"x 0" o „ 0" „ , xy y y z °* 2 o y 2 [A] 1 (3.13) 26 '<*1 r j 2 0 6 c2 s2 0 2sc 0 0 s2 c2 0 -2sc 0 0 0 0 1 0 0 0 -sc sc 0 c 2-s 2 0 0 0 0 0 0 c s 0 0 0 0 -s c xy yz (3.14) where c s cos 0 and s = sin 0 and [T] is termed the transformation matrix. The stresses in the x-y system can be written ° x ° 1 < • - [T]~ 1 . ° 3 > ° x y CTxz a 6 . (3.15) The strains are transformable in the same manner: S x £ 2 s y £ 3 Y12/2 ' = [ T } Y x y / 2 Y 1 3 / 2 Y x z / 2 J 2 3 / 2 ^ / 2 (3.16) Substituting Equation 3.16 into Equation 3.7 and then the result into Equation 3.15 we have: 27 ° x ay o\ s. z Y x z /2 YyZ/2 [c 'K«'} (3.17) The same is true for a 2 dimensional state of stress. Referencing Equation 3.11: ° " x . = [Tr[Q][T} • = [ Q ' K s ' } a*y, Y „ / 2 (3.18) 28 3.2.2 Elasto-Plastic Regime Beyond the elastic range, there is no longer a one-to-one relationship between stress and strain and the material becomes load path dependent. The classical incremental theory of plasticity accounts for this loading history dependent behaviour while considering the effects of stress interaction under multiaxial stress states. In essence, the incremental theory of plasticity relates the infinitesimal increment of plastic strain components (dep;j) to the stress state (o;j) and the stress increment (do;j). The constitutive relationships for the elasto-plastic regime are formulated through consideration of three fundamental concepts (Chen and Han, 1988): 1) the existence of an initial yield surface, 2) the hardening rule and 3) the flow rule. These concepts are explained in the following three sections. 3.2.2.1 Initial Yield Surface Similar to the concept of a failure surface, the initial yield surface defines the combination of stress components (o) at which elastic behaviour ends and plastic deformation begins. The quadratic function for yielding of an anisotropic plastic material is written in general form as (Shih and Lee, 1978) {(oi,ai,Mli,^) = 0 (3.19) In this study, the specific form of the yield function is taken to be f = a 2 ( a i , a , , M i j ) - k 2 = 0 (3.20) where: a is an effective stress or equivalent stress and k is a threshold stress which is equal to the size 29 of the yield surface. For the purposes of plastic formulation, the Tsai-Wu criterion is adapted to this general form as follows: The square of the effective stress is conveniently defined as a2 = M,(ai-ai)(p,-ai) (3.21) where the term M ; j describes the shape of the yield surface and a; = {at , a 2, a4}T = {a{ , a2, 0}T "describes the strength differential between the tensile and compressive strength or the 'offset' of the origin of the yield surface" (Shih and Lee, 1978), so that: f = M l i ( a i - a i ) ( a j - a . ) - k 2 =0 (i,j = 1,2,4) ie., f = M n a , 2 -2a}a^ +ax2 + M 22 a 22 - 2 a 2 a 2 + a 2 2 + 2M 1 2 [ a ] a 2 - a t a 2 -cs2ax + a 1 a 2 ] + M 4 4 - a 4 2 - k 2 = 0 Setting the equalities (3.22a) L j = 2 M i j C t j ; K = - M y a j a j + k 2 (3.23a ; 3.23b) one obtains: f S M i j o i o j - L i a i - K = 0 (3.24) Comparing Equation 3.24 with Equation 2.3, it is concluded that: 30 M . ^ K F ^ K ! F n F 1 2 0 F ] 2 F 2 2 0 O O F , 66 L: = - K F : = - K (3.25a; 3.25b) It is noted that the components of M ; j are not independent. For example, if M n = 1 then K = 1/F U and M 1 2 = F 1 2 / F u , etc... (Shih and Lee 1978). Furthermore, <x; can be found by substituting Equation 3.25a into 3.23a, then equating the result to 3.25b L , = 2 M l j a j = 2 K F i j a j = - K F , (3.26) and solving the resulting simultaneous equations: Fi = -2F l j a j (3.27) Finally, the square of the threshold stress is calculated from Equation 3.23b as k = K + M^oc ; a j (3.28) When referring to initial yielding, Equation 3.20 is written in terms of the initial parameters, a° , M::°, and k„ . ie.: f E a ^ a . ^ . M , 0 ) - ^ 2 =0 (3.29) 3.2.2.2 Hardening Rule Beyond the initial yield surface, material response is modeled according to inherent material behaviour. A material whose plastic deformation is assumed to occur under a constant yield stress level is termed perfectly plastic. For such a material, the yield surface remains that of Equation 3.29 31 throughout plastic flow. However, a material that exhibits stress dependency on plastic straining is said to undergo work or strain hardening and is modeled by specifying a new yield surface at each stage of plastic deformation known as subsequent yield surfaces. This post yield response, whereby the yield surface changes shape and/or position, is managed by way of a hardening parameter (x) such that each new yield surface satisfies: f ^a 2 (a 1 , a 1 ( x ) ,M l j ( x ) ) - k 2 ( 5 C ) = 0 (3.30) The hardening parameter relates the affected variables (a ;, M ; J , and k) to the degree of plastic deformation. In this s t u d y , w i l l be taken to be equivalent to a strain variable called the effective plastic strain e p (actually, JdsFp ). This variable is a scalar value function of the plastic strains. Referencing Equation 3.30, one can establish the state of loading for a material under a multiaxial stress state. The criterion for a "loading state" (ie. when additional plastic deformation will occur) is df i f f = 0 and - — d o ; > 0 , then d% * 0 (3.31) 9a: i and the stress point remains on an expanding yield surface. During "neutral loading" df i f f = 0 and -—da: = 0 , then dx - 0 (3.32) cc j which is plastic loading for a perfectly plastic material and the stress point remains on the constant yield surface and finally, for. a state of "unloading" we have i f f = 0 and - — d a : < 0 , then d% = 0 (3.33) da; 32 and the stress point returns inside the yield surface. Equation (3.30) represents a general anisotropic hardening model which specifies the rule for the evolution of the subsequent yield surfaces. Different hardening rules are formed by imposing different conditions on parameter variation with %. If parameters a ; and M ; j are assumed to remain constant through plastic flow and k varies with %, then the subsequent yield surfaces expand uniformly around the original yield surface, as in Figure 3.3, and the hardening model is said to be isotropic. If, on the other hand, both k and M ; j remain constant, while a : vary with %, the subsequent yield surfaces do not change shape but instead translate in stress space as a rigid body. This model is referred to as kinematic hardening. In the present study, a special type of general anisotropic hardening model is used whereby k varies proportionally with j(, « i remain constant while the nonlinear yield parameters of M ; j vary non-proportionally with plastic deformation. This model leads to an expansion and at the same time distortion of the yield surface upon plastic strain. A similar approach was taken by Vaziri et al. (1991). This study, however, did not account for the difference in strength between tensile and compressive modes (ie. at was ignored). Under the present formulation, Equation 3.30 becomes f = a 2 ( a i , a i , M l j ( x ) ) - k 2 ( x ) = 0 = M l j ( x ) ( a 1 - a 1 ) ( a j - a j ) - k 2 ( X ) = 0 The mathematical details of this model are presented in Appendix A . 33 Isotropic hardening Kinematic hardening Figure 3.3 - Subsequent Yield Surfaces for Hardening Material 3.2.2.3 Flow Rule In accordance with plasticity theory, under yielding, the total incremental strain vector is assumed to be composed of both incremental elastic and plastic components through simple superposition such that det = def + ds^ (3.35) The relationship between elastic strain increment and stress increment is defined through an infinitesimal form of Hooke's Law as dorj = C J : dE: (3.36) In what follows, the relationship between the total strain and stress increment beyond yield is derived. The flow rule enables a relationship to be made between the plastic strain component and 34 the stress increment by introducing the concept of a plastic potential function, g (o;, a ;, k). In particular, the flow rule is defined by classic plasticity theory as d s 1 p = d ? , — ( 3 . 3 7 ) do j where de ; p is the plastic strain increment vector, dX (termed the plastic multiplier) is the magnitude of the vector and dg/3o;, the stress gradient of the plastic potential function, defines the direction of the vector. For simplicity as well as computational efficiency, an associated flow rule will be assumed whereby the plastic potential function, g (o;, a ;, M ; J , k), is considered equal to the yield function, f (o;, a ;, M ; J , k) so that the direction of the plastic flow vector is normal to the yield surface as illustrated in Figure 3.4. Figure 3.4 - Illustration of Normality of Plastic Flow Vector for Associated Plasticity The hardening process of a material under a multiaxial stress state is characterized through one calibrated uniaxial stress-strain curve, using the effective stress (o) and the effective plastic strain (ep), which are scalar functions of the stresses and plastic strains respectively. The effective stress is 35 defined by Equation 3.21. To define the effective plastic strain we consider a simple combination of infinitesimal plastic strain increments, the simplest type with the correct 'dimension' is given by Chan and Han (1988) as ( d s ? ) 2 = C - d S i 1 5 - d S j p (3.38) where C is a constant. If we choose C = M^" 1 (the inverse of M ; I , providing | M ; J | *0), then we can simplify 3.38 using the flow rule to ( d B P ) 2 = M 1 - 1 . d ^ 2 - ^ ^ ' 9a j 9a j = M 1 ~ 1 -dX2 ^ M ^ a j - c t j ) - 2 M i k ( a k - a k ) (339) = dX2 • 4 M i j ( a i - a A ) • (aj - aj) - dX2 • 4Vk Thus, dX = — (3-40) 2k From expression 3.34, we find that di — da „ , 3a = 2 a - — = 2 k - — (3.41) da^ 3CTJ dcfl Designating = a j , and substituting Equation 3.40 and 3.41 into the flow rule we have 36 d S i P = d ^ - 2 k - a i d s p 2 k - a , (3-42) 2k = d £ p -a; A one to one relationship between effective stress and effective plastic strain is assumed as follows: G = H(S~ p ) (3.43) where, upon differentiating D C 7 = H ' or, - ^ - = H ' (3.44a), (3.44b) ds"p d s p H ' is termed the hardening (or plastic) modulus and is associated with the rate of expansion of the yield surface. In the present study, a trilinear stress-strain curve is assumed such that H ' is a constant. Thus, the yield surface expands in accordance with k = k o + H ' 8 p (3.45) Substituting expression 3.35 into 3.44 while considering uniaxial yielding (typically in one of the principal material directions), we can determine the function H ' as: H , = da da 1 ds p d s - d s p ds__ds^ (3-46) da da where, in the present study, we obtain the elastic stiffness, E l c = de/do and the tangent stiffness E l c ' = dep /do from the stress-strain curve of the compression parallel-to-grain test data. Thus, we have 37 H' = F ' (3.47) During plastic flow, each new state must satisfy the subsequent yield function so that f ( a i , a i , M i j , k ) + df = 0 (3.48) or, in this case A F J d( J A / R di da{ dM{- ' 5k df = - — d a ; + - d M i ; + — dk = 0 (3.49) which is known as the consistency condition. This condition assures that the subsequent stress and strain states remain on the subsequent yield surface. We consider each term, in turn. Using Equations 3.41 and substituting 3.35 into 3.36, then considering 3.42, we have di dd-. -da. = 2ka ; •dai = 2 k a i - C i j ( d s j - d s j p ) ( 3 , 5 ° ) = 2 k a i - C i j ( d s j - a j d s p ) Next, referencing expression 3.34, noting d% = de p and using Equation 3.44b 38 di dMij =(CTJ -a i )(o- j -ct jJdMij 5 M X J = ( a i - a i ) ( a j - a i ) - ^ — d x (rji -ajXCTj - a j ) — f dsP (CTi -ai)(aj -a : ) 3X 5Mij 5e p 5k dk C i V l : : ( ^ 1 - a 1 ) ( ^ j - a j ) - ^ - L H ' d 8 P (3.51) Referencing Equation 3.34 and 3.44b — dk = - 2 k H ' d s p 5k (3.52) The consistency Equation 3.49 becomes: df = i L d C T i + J L d M , . + ^ d k = 0 8G1 8Mi} " dk = 2 k a i C l j e ( d s j - a j d e p ) + II • H ' d e p - 2 k H ' d e p = 0 (3.53) w here IT = (o; - ai)(aj - o^ ) dk Solving for de p , we have ds~f H ' V 2k7 -de: + a„Cl„a m mn n (3.54) 39 Finally, substituting this into Equation 3.42 and the result into Hooke's law do; = Cife (def - d6jP) the relationship between the incremental stress and total incremental strain beyond yield is der = C?: Cjk ak aiCij H' 1-2k d s : = c 5 ds, = C^d8, (3.55) The plastic stiffness tensor, Qj P represents the degradation of the elastic stiffness tensor, due to plastic flow. A similar formulation for fiber-reinforced composite laminates is provided in Vaziri et al. (1991). 3.2.3 Post Failure Regime Onset of material failure is determined through implementation of a failure surface which serves as an upper bound to all other loading surfaces and remains unchanged throughout loading. It is similar in form to Equation 3.34, ie. f - a V i . o ^ M ^ - k 2 , = 0 ( 3 5 6 ) = M D ( a 1 - a 1 ) ( a j - a j ) - k 2 =0 where the suffix V denotes ultimate value. The values M^" and k u are material constants which come from ultimate strengths obtained from experimental tests in the principal material directions described in Chapter 5. 40 3.3 Commentary on Plane Stress Formulation It is noted here that the same two dimensional form of the Tsai-Wu criterion is used for both the two dimensional and three dimensional formulation. The rationale behind using this form of the yield/failure criterion for the 3 dimensional model is the following: a) The cases investigated in the current study (which reflect most applications for structural composite lumber) are assumed as being in a state of plane stress. Thus, a plane stress criterion was deemed adequate. In any case, the three dimensional model will produce estimates of the out-of-plane stresses, at which point, a judgement could be made as to whether or not the assumption of plane stress is correct. b) A thorough investigation of parameters required for the three dimensional form of the Tsai-Wu criterion is too time and cost prohibitive within the scope and realm of the current thesis. As the present study is the first investigation of the chosen modeling technique, the two-dimensional model was deemed to suffice. It was rationalized that, if found necessary from this thesis, a more extensive criterion could be used for future analysis. 41 4 - Finite Element Formulation 4.1 I N T R O D U C T I O N The finite element method is a well established numerical technique used for analyses of structures and continua (Cook et al. 1989; Zienkiewicz, 1971). Formulated in a computer program, the method is a powerful and versatile engineering tool with a wide range of applicability. For the present study, the finite element method provides a procedure for determining displacement and stress distribution throughout structural members considering nonlinear material behaviour. The constitutive model outlined in the previous sections has been implemented into a stochastic, materially nonlinear finite element based FORTRAN 77 program with extended capacity to perform Monte Carlo simulations entitled COMAP (COMposite Analysis Program). The program was adapted from a 2 dimensional nonlinear program provided by Owen and Hinton (1980). Both a 2 dimensional and 3 dimensional version of the program have been developed and investigated. This chapter systematically presents the general foundations of, as well as the underlying assumptions made in developing both versions of the program, the numerical implementation of the constitutive model of the previous chapter into the finite element analysis and finally, verification of both models through comparison with proven numerical or analytical solutions to given problems. 42 4.2 Two D I M E N S I O N A L L I N E A R F I N I T E E L E M E N T M O D E L Very often in the analysis of laminates the lamina is assumed to be in a simple 2 dimensional state of stress (or plane stress). As such, it was decided that a 2 dimensional approach to strength prediction would first be investigated. It was reasoned that if successful, the 2 dimensional approach would provide a cheaper and faster solution than that of a 3 dimensional model. In this section we summarize the general finite element technique for analysing composites using 2 d formulation. This formulation is described in more detail in Ochoa and Reddy (1992). 4.2.1 Displacement-Based Element Representation A structural member may be discretized and mathematically represented as a series of finite elements. The isoparametric formulation makes it possible to represent elements which are distorted as shown in Figure 4.1b. For the 2 dimensional model used in this thesis, plane bilinear isoparametric elements with 4 nodes and 2 degrees of freedom at each node are used. 1) \ p (1,1) ' ® ® ® (-1,-1) (1,-1) (a) (x4,y4) Owi) (X2,y 2) (b) Figure 4.1 - Plane Linear Isoparametric Element: (a) Rectangular element; (b) Distorted element While the four corners of the distorted element are defined by global Cartesian coordinates (x, y), displacements and coordinates are expressed in terms of curvilinear coordinates (£, r\) through 43 mapping. In keeping with isoparametric formulation, the displacements or the coordinates within the distorted element can be interpolated from nodal values as, u = ! N i ( 4,Ti ) - u i ; v = X N , ( ^ T i ) - v i (4.1) i i and x = l N i ( 4,Ti)-x i ; y = i N i ( ^ T , ) . y i (4.2) i i where, N;(£, r|) , known as shape, basis or interpolation functions, are chosen so as to preserve continuity of the displacements between elements. For this study, they are N i (4,Ti) = i ( l + ^ i ) ( l + T,Ti i) (4.3) 4 where i is the number of the element nodal point, and , r); are the nodal coordinates in the £ - r\ plane. 4.2.2 Stiffness Matrix Formulation 4.2.2.1 Compatibility With a view to solving for the nodal displacements, the stiffness matrix for each element must first be obtained. Letting the elemental nodal displacements be written in vector form as {u}e = {u;, v J T , the strains in the x , y coordinate system are defined using the strain-displacement relations as 44 Y x y u i du dx dx 3N, 5x 5x dx 0 0 0 0 U 2 U 3 cV dy 0 0 0 0 5N, ay ay aN 3 ay aN 4 ay U 4 V l du -dv 3N 2 5N 4 aN 2 3N 3 aN 4 V 2 dy dx [ 5y 5y dy 5y ax dx dx dx V 3 7 4 . (4.4) or, aN, "ax~ o aN, ay 0 aN, ay 5N-(4.4a) where [B] is designated as the strain displacement matrix. Here, the shape function derivatives are written with respect to the (x, y) coordinates, whereas the shape functions themselves are defined as functions of £ and r\. We can write by use of the general chain rule of differentiation aN ; L x_^_ '~~b\~b\^~b\"bx (4.5) 3N ; dN ; dx 5N, dy ~dx~~dj + 3N.. dx dy at; dNi dy dr\ dx dr\ dy dx\ • + such that the required derivatives 5N/dx and 5 N / 5 y are obtained by inversion as ' a N , ' ' a N / dx ' aN, 5 • aN ( U j an. (4.6) where [ J ] is called the Jacobian matrix defined by 45 [j] = dx dx an dj_ (4.7) 4.2.2.2 Equilibrium Equilibrium of the structure is assured through application of the principle of virtual work whereby the internal virtual work (8W i n t) is equal to the external virtual work (5W e x t), ie.: § W i n t = f v { 8 e f ^ o - ' } d V = 5W e x t = j ^ S u ^ - f F j d V (4.8) where {a'} = {ax, o y ) T ^ } 1 , {F} are external body forces and {8u} are the virtual displacements at the points of application of the forces. Considering Equation (4.4a) and simplifying, the equilibrium equation for an element is written in matrix form as U B ] T { a f d V ' = {F} e ( 4.9) where {F} e is the equivalent external force acting on the nodal points. 4.2.2.3 Constitutive Relations In order to resolve Equation 4.8, we must first consider the constitutive relationship between {o'}e and {e'}e. In this 2 dimensional analysis, we employ classical lamination theory (CLT) to analyse the properties of the structural laminate. Fundamental to the C L T , it is assumed that a state of plane stress exists in each layer; ie. the stresses in the through thickness (z) direction (reference Figure 2.1) are negligible. Also, the layers are 46 assumed to be bonded perfectly together, so that there is no slip between plies. Furthermore, reflecting the conditions of this study, the analysis presented here is restricted to that of a symmetric laminate (with respect to its mid-plane), subjected to solely in-plane forces such that no bending or twisting moments exist. The consequence of using a symmetric laminate is that no coupling is assumed between bending and extension (ie. in-plane loads will not create bending or twisting curvature). Referencing Figure 4.2, and using Equation 3.18, the stress-strain relationship for the k t h lamina is expressed in terms of the mid-surface strain {e °} as {-'}:=[Q '];KK=[Q']:K o } e ^ > \ A 2 1 z 2 Middle Surface > t z* ZN-1 k f Z N > 1 n 1 V ^ Layer Number Figure 4.2 - Geometry of N - Layered Laminate (Jones, 1975) The resultant laminate forces, {R} measured per unit width, are obtained by integrating the stress components for each lamina over the entire laminate, ie. 47 t/2 n CTx '= 1 • a y • dz = I" ! • CTy • dz > - t /2 k k=l CT*y. k . (4.11) where t = laminate thickness, n = total number of layers, and z k and zk.j are as defined in Figure 4.2. Substituting Equation 4.9 into 4.10 and noting that the stresses are assumed constant through each lamina, {N} = E[Q'Ltt{e'°} k=l (4.12) At this point, we can solve for the element stiffness matrix. Substituting Equation 4.10 into 4.9, we have U B ] T [ Q ' ] ; K d V - = { F } -Further substitution of Equation 4.4a yields or, (4.13) (4.14) Thus, considering Equation 4.11, the element stiffness matrix is calculated as M'=xaA.MT[Q'];[BH*dy k=l (4.15) where the summation is taken over all layers in the laminate. Converting to curvilinear coordinates, we use the relation dxdy = derjjjd^dr) (4.16) 48 where det[ J] is the determinant of the Jacobian, to get M* = Z *k J J[B]T[Q'K[B] det[j]d^ dTi (4.17) k=i -1-1 The integrand in this equation is solved numerically using Gaussian quadrature with a 2 x 2 sampling scheme as follows: M = 111 tk [ B ] T [ Q ' ] ; M det[j]W i W l (4.18) k=l j=l i=l where, in this case, the calculation is carried out at % = ± l/-J3 , r| = ± l/V3 and W ; and Wj are associated weighting factors equal to 1. The element stiffness matrices [k]e are added element-by-element, using standard assembly guidelines, to form the global stiffness matrix of the structure [K]. Knowing [K], the applied loads and appropriate boundary conditions, the unknown nodal displacements for the entire structure are found from the simultaneous solution of the equations [K]-{u} = {F}. 4.2.3 Two Dimensional Model Limitations As noted in section 4.2.2.3, it is implicitly assumed in using 2 dimensional planar elements that interlaminar stresses associated with the through-thickness (z) direction are negligible. This is a valid assumption if 1) the laminate is symmetric, 2) the applied loads on the laminate are statically equivalent to in-plane forces producing neither bending nor twisting and 3) "free-edge" effects are negligible. The present study complies with the first two criteria. The third criteria leads to a phenomenon known as delamination which has been extensively studied both theoretically and experimentally for many advanced composite materials (Pipes and Pagano 1971; Whitney and 49 Browning 1972; Gibson 1994). These studies show that within a "boundary region" (roughly one laminate thickness inward from the free edge of the laminate), a 3 dimensional stress state exists. Unfortunately, no such studies have been carried out on wood composites. It was decided that for further understanding of the material behaviour as well as providing a more versatile tool which could be used for more complex, non-symmetric composites, a 3 dimensional model would be investigated. 4.3 T H R E E D I M E N S I O N A L L I N E A R F I N I T E E L E M E N T M O D E L The 3 dimensional isoparametric procedure closely resembles the 2 dimensional development discussed in the previous section. For brevity, the following will concentrate primarily on the differences between the two models. This 3 dimensional finite element formulation is standard and can be found in any basic finite element text (Cook et al., 1989; Yang, 1986). 4.3.1 Displacement-Based Element Representation The 3 dimensional element geometry is that of a linear solid (eight-node-brick) element as shown in Figure 4.3. n Figure 4.3 - Linear Solid Isoparametric Element 50 The displacements and the coordinates of a point are now defined in 3 directions as u = E N i ( U O - u i ; v = 2N i (4 ,Ti ,0-v i ; w = 2 N i ( ^ , i l , 0 - w i (4.19) and i = E N i ( U O - i i ; y = ZN i(^,Ti ,0 .y i ; 2 = l N i ( 4 , T , , Q . 2 i (4.20) The corresponding linear shape functions are found to be N i (^n , g = i ( l + ^,)(l + nni)(l+C<;i) (4.21) where ^, n ;, and £; are the coordinate values of node i , with i = 1, 2, ... 8. 4.3.2 Stiffness Matrix Formulation 4.3.2.1 Compatibility The strain displacement relationship in 3 dimensions is written as dN, 0 0 ~dy~ aN, ~~dV 0 0 aN, ay 0 aN; dx 0 aN, "aT" o o a N , ~az~ o aN, ax aN, = [ B ] W (4.22) Because N ; are defined in terms of I;, n, and C, we use the inverse of the Jacobian matrix 51 dx dy dz 34 34 34 ax ay dz dr\ dx dy dz (4.23) to find the derivatives of the shape functions in terms of x, y and z, as 'aN;' aN,' dx aN, aN ; dy ar| aN, aN, I dz \ (4.24) 4.3.2.2 Equilibrium By applying the principle of virtual work, the equilibrium equation for an element is the same as for the 2 dimensional case, ^[BfjafdVMF} 6 (4.25) 4.3.2.3 Constitutive Relations From Equation 3.17, the stress-strain relationship for each element in the (x,y,z) laminate coordinate system (see Figure 2.1) is defined as {a ' } e =[C' ] e {s '} e (4.26) Substituting Equation 4.26 into 4.25, we have 52 By further substitution of the strain displacement relationship (Equation 4.22) into Equation 4.27, L [B] T [C 'r [B]dV' .{u} ' = {F} ' (4.28) we find the element stiffness matrix in the usual manner as [ k] = U B ] T [ C ' ] e [ B ] d x d y d 2 (4.29) Converting to curvilinear coordinates, we use the relation d x d y d z = det[j]dc;dridc; (4.30) to get M = 1 I 1[B]T[CT[B] det[j]d4dndC (4.31) Finally, using Gaussian quadrature with a 2 x 2 sampling scheme: = i Z i : [ B ] T [ C ' ] e [ B ] d e t [ j ] W i W j W k (4.32) k=i j=i i=i The calculation is carried out at % = ± l/V3 , r\ = ± l/V3 , C, = ± l/>/3 and W ; , W j and W k are the associated weighting factors equal to 1. As in the 2 dimensional model, the element stiffness matrices [k]e are assembled to form the global stiffness matrix of the structure [K]. This matrix relates the nodal displacements for the entire structure with the applied loads through the relationship, {F} = [K]{u} (4.33) 13 4.4 N O N L I N E A R F O R M U L A T I O N In a nonlinear analysis, the governing equation is 4.25. This equation represents the equilibrium of the externally applied force vector {F} with an internal resisting force vector {P} = J[B]T ja'(u)jdV which is a (materially) nonlinear function of the nodal displacements {u}. v The difference between the two vectors can be interpreted as a load imbalance. Typically, the finite element solution proceeds on an incremental basis and solution of these equations requires an incremental/iterative numerical algorithm which effectively minimizes, to an acceptable tolerance, the residual { ¥ ( u ) } = { F } - { P ( u ) } (4.34) There are numerous methods for solution of the nonlinear system of equations 4.34. The method used in this study is a modification of the classical Newton-Raphson method. 4.4.1 Modified Newton-Raphson Method As found in any standard text on nonlinear finite element formulation (Owen and Hinton, 1980; Cook et al. 1989), we first consider a truncated Taylor series expansion of {x¥} r about {u}r_j : 3{u] M (4.35) where the subscript r denotes the iteration number and {Au} = {u}r - {u}^. Referencing Equation 4.34, it is noted that a{u} 8{u} (4.36) where [K] r 4 is the global stiffness matrix evaluated at the beginning of the r t h iteration. Substituting 54 this into Equation 4.35 and considering that ultimately we want {XP}[ = 0, we have { * L . = [ K L - { M (4-37) We use this set of equations to compute {Au} whereupon we can update the displacement estimate to R = R-, + M (4-38) The new found displacements are used to calculate a new internal resisting force vector {P} whereby we re-evaluate the residual {¥} (equation 4.34). For the next iteration, we obtain a new structural (tangent) stiffness matrix [K] and repeat the process until the residual converges to the point of satisfying a prescribed tolerance. In the present study, we have adopted an initial stiffness routine (illustrated in Figure 4.4 for a single variable situation), in which the tangent stiffness corresponds always to the initial value to avoid zero or negative definite stiffness beyond the peak load. This also reduces the costs associated with reproducing the tangent stiffness matrix for each iteration. However, more iterative cycles are needed within each iteration to reach the set tolerance. 55 u 0 u 2 u 3 u Figure 4.4 - Initial Stiffness Solution Algorithm for a Single Variable Situation (Owen and Hinton, 1980) The convergence criterion used is one based on the residual force values (Owen and Hinton, 1980): )/fflr W < T O L E R A N C E / 4 _ 3 9 N WW) For this study, the tolerance was set to be 0.01. 56 4.4.2 Elastic-Plastic Formulation As alluded in the previous section, the equations of equilibrium are resolved in an incremental manner. During any iteration within an increment in load (or displacement), a part of an element may yield or fail. To determine this, all stress and strain quantities are monitored at each Gaussian point throughout the entire stress history of the member. The purpose of this section is to detail the procedure used in monitoring and adjusting the stresses and strains throughout the elastic-plastic regime. Yielding or failure of a point is established through use of the Tsai-Wu criterion. In the criterion, stresses are defined with respect to the principal material coordinate system. Thus, the stresses are transformed from the x-y coordinate system in accordance with the transformation relationships established in Chapter 3, before being used in the criterion. In the following discussion, stresses are denoted {o} = {at, o2, o 4}T or {o1 ; o2, o3, o 4 o 5 a 6 } T for the 2-dimensional model or 3-dimensional model, respectively. Similarly, strains in the principal material directions are denoted, {e}. We follow the perspective of one Gaussian point within the member during the r t h iteration. Al l stress and strain quantities associated with that point as well as the state of the yield surface are known from the (r-l) th iteration. 1. Using the residual forces from the previous iteration (ie. {¥},..!) which is the load imbalance, we calculate the resulting displacement (Equation 4.37) and from this, strain increment { Ae} (Equation 4.4a (2-D)or 4.22(3-D)). 2. Assuming elastic behaviour, we compute the elastic increment in stress {Aoe} from 57 Equation 3.7 (3-D) or 3.11 (2-D). 3. We calculate a trial stress, {oe} r which is the accumulated total stress for the point according to { a l = M . - , + M r (4-40) 4. We use our trial stress to calculate the effective stress o r = (M^^-a^o-a^2 and check this against the initial yield criterion (Equation 3.29); ie. we check if o r ^kD . If the answer is no, then the elastic assumption is valid and the trial stress is indeed correct. If the answer is yes, however, the point has either a) yielded and the point continues to deform ductilely (as a result of being compression - dominant) or b) failed in a brittle manner (as a result of being tension • dominant). Tension vs. Compression Dominance The decision between compression or tension dominance is made depending on the combination of stresses at the point of failure. We define conditions for tension dominance for the two and three dimensional models as shown in Table 4.1. To describe these conditions, we need to categorize the individual components of the effective stress (referencing Equation 3.22b) as: p2 = M 2 2 2 2 a 2 - 2a 9a ? + a 2 (4.41) p4 = M 4 4 - a 4 These components (p1; p2, and p4) are chosen to reflect the magnitude of the stresses parallel-58 to-grain, perpendicular-to-grain and in-plane shear, respectively. Table 4.1 - Conditions for Tension Dominance 2 D Model 3D Model o2^Y t o2^Y t l°6|.^s 1 o6|*S o,^ 0 and p! ^p2 o^O and Pi ^ 2 o2>0 and p2 > pj o2>0 and P2 * Pi P4 * Pi and p4 * P 2 S > or At an integration point, the stress state is deemed tension-dominant, and thereby brittle, when the failure criterion is violated and any one of the conditions shown occurs; otherwise, it is deemed compression-dominant. The rationale as to whether the point yields or fails is albeit subjective. It was chosen to provide good results when compared to experimental data. Tensile failure leads to the post-failure regime and will therefore be deferred for discussion in the next section. Stress combinations other than that defined by Table 4.1 result in ductile behaviour and follow plasticity formulation. We continue with elastic-plastic formulation. Given that the trial stress has been determined to be compression dominant and or ^ k0, then the stress vector has crossed the initial yield surface. This stress path is illustrated 59 vectorially in Figure 4.5 (Owen and Hinton, 1980). Aa; fl-RVAoV J ^ ^ ^ ^ ^ \ A ' " a. Initial yield surface 0" = k 0 i ^ I Susequent yield surface r o" = k. Figure 4.5 - Incremental Stress Changes at Initial Yield For yielded Gauss points, we wish to proportion the total stress appropriately into elastic and elasto-plastic portions. For this, we first compute a contact stress {oc} (entirely elastic) where the stress path just comes into contact with the initial yield surface, ie., { ° e } = { ° L . + f l - R ) { A c e } ( 4 - 4 2 ) where, R is determined through simple linear interpolation: a ! - k R = • a - a (4.43) r - l The remaining portion of the elastic stress increment, R-{Ao e} lays beyond the initial yield surface (which violates the yield criterion) and must be corrected to account for plastic 60 deformation. It is noted that if the Gaussian point had previously yielded, then R = 1 and the entire stress state is elastoplastic. Referring to Figure 4.5, we wish to find the total stresses {or} which satisfy elasto-plastic conditions. Substituting incremental changes for infinitesimal changes in Equation 3.55 we have {Aa r } = [ G - ] { A s r } = [ c e ] { A 8 r } - [ c " ] { A 8 r } (4.44) = { A < } - { A c ? } Thus, when incrementing from {o^}, we can calculate , , [clMM T[ci { ° . H a , _ , } + A < - I I I { A s , } (4.45) where all values are as defined in Section 3.2.2. The final stress point, corresponding to point A on the schematic, may not fall directly on the yield surface (ie. at A ' ) . The discrepancy is aggravated with larger increments however can be resolved by simply scaling the individual stress components of {or}. As the effective stress, o should coincide with the threshold stress, k, an appropriate scaling factor is found to be (4.46) If large increments in load (or displacement) is allowed in the preceding process, the final point A ' may contain some error. As explained in Owen and Hinton (1980), the inaccuracy may be reduced by relaxing the excess stress to the yield surface in several stages. 61 Referencing Figure 4.6, the excess stress is divided into m parts where m is the closest integer less than G—k •8 + 1 (4.47) Figure 4.6 - Stepped Process to Reduce Stress Point to Yield Surface Using a trilinear representation of the stress-strain curve, the threshold stress, k is updated after each iteration according to Equation 3.45: k r = k 0 +H'e p r . In accordance with the anisotropic hardening rule adopted in this study, the nonlinear strength variables, X c and Y c are also updated after each iteration, as described in Appendix A . 7. Finally, from the total stresses established for each Gaussian point, whether in an elastic state or an elastoplastic state, appropriately transformed back to the x - y coordinate system, we compute the equivalent nodal forces 62 (PL = J[ B] TM d V (4.48) v which are used in calculating the residual forces {¥} for the next iteration, if deemed necessary by the convergence criterion Equation 4.39. This iterative process repeats for each load increment until the stresses reach the ultimate strengths per Equation 3.56. At this point, failure occurs. 4.4.3 Post Failure Regime When failure has been established, post-failure behaviour ensues. Failure is described as being either brittle or ductile, depending on whether the combination of stresses at the point of failure result in a tension or compression dominant state (see section 4.4.2). Referencing Figure 3.1, when ductile failure occurs, the yield quantities, X c and Y c have reached their ultimate values and the stress level remains constant, reflecting that the material cannot support additional stress. Stresses are not permitted to move outside the ultimate yield surface and thus can only traverse the surface until both equilibrium and the constitutive relations are satisfied. When brittle failure occurs, however, the stress level is reduced gradually to facilitate convergence of the iterative procedure. The plane stresses {o1 ; o2, o4} are reduced proportionate to the same values in the previous load increment. The ratios used differ slightly for the 2-dimensional and 3-dimensional model and are summarized in Table 4.2. These values were established through comparison with test results to provide good post-peak results. 63 Table 4.2 - Brittle Failure Stress Reduction Ratios Stress 2 - d Model 3 - d Model o, + (tensile) 0.70 0.90 o{ (compressive) 0.98 0.98 o 2 + (tensile) 0.90 0.90 ° 4 0.95 0.95 In summary, when a local integration point is detected to have failed either in a brittle or ductile manner, the global material stiffness is modified until satisfaction of equilibrium is attained within the iterative process. Failure of a portion of one lamina is compensated for by an increase in the load carried by adjacent Gaussian points or laminae. The load is then further increased until final failure is reached. Final failure occurs when the load displacement curve begins to descend and the program is automatically halted when the load reaches 90 percent of the peak load. In rare cases, the program may also be halted due to nonconvergence beyond peak load. The 90 percent condition was chosen both in consideration of computer time as well as numerical instability. 4.5 P R O G R A M V E R I F I C A T I O N The theory presented in the foregoing chapters has been implemented into a finite element based computer program, COMAP. As with any new program, however, a thorough check is initially necessary before it can be used with confidence. It is the objective of this section to verify the general performance and results of COMAP. As no other software exists which has the same capabilities as our program, the verification process involves comparison with several programs -each targeting a specific aspect of COMAP. As well, a problem with an exact analytical solution is investigated to gauge the robustness of the program's prediction capabilities (within the bounds 64 of the theoretical assumptions). 4.5.1 Isotropic Cantilever Analysis In order to assess the plastic deformation capabilities of the program, the classic example of an inelastic cantilever beam, shown in Figure 4.7, is investigated. The example considers the case of a elastic - perfectly plastic material (such as structural steel), with rectangular cross section, limited to small deflections and neglecting the effects of shear on deflection. The problem has known analytical solutions for both elastic and elastic-plastic Figure 4.7 - Cantilever Beam behaviour (Timoshenko, 1972). For the elastic range, it is Example known that the yield load F y which first produces yielding of the beam is given by M v F Y = ^ (4.49) The deflection at incipient yield is F L 3 5 = i L (4.50) y 3EI In nondimensional form, we have a simple relationship between deflection and load in the entire elastic range, _8___F_ 5. ~ F. F F ^ 0< —<1 F (4.51) ' y V y In the region of elastic-plastic behaviour, the relationship is shown to be 65 f F v l 2 f F ! y 5 - 3+ — I F J 3 - ^ F y F ^ 1 < — < - (4.52) A graph of load F / F y versus deflection § / 8 y through the entire stress path is shown in Figure 4.8. This result forms the basis of comparison for the numerical analysis using the program C O M A P . The load displacement curve is reproduced using both the 2-dimensional and 3-dimensional version of the program assuming an isotropic material. The cantilever beam is discretized for both cases as shown in Figure 4.9. Linear 4-node, or 8-node elements are used as shown, together with a 2x2 (or 2x2x2) Gauss integration scheme. Load was applied incrementally at the free end of the beam lumped proportionally on the nodal points. The assumed isotropic material properties are summarized as follows: Elastic Modulus: E = 30 x 103 MPa Poisson's Ratio: v = 0.3 Uniaxial Yield Stress: o y = 30 MPa Hardening Modulus: H ' = 0 66 1.6 t 8 / 5 y Figure 4.8 - Load versus Displacement for Isotropic Cantilever Beam Example As the program is configured for a stochastic analysis with regards to the elastic moduli and yield stresses, the values aforementioned are considered to be the average values and standard deviations are zero. The load F was incremented until plastic collapse of the beam was attained. Collapse was deemed to have occurred when the iterative procedure diverged for an incremental load increase. The load increments are shown as crosses in Figure 4.8. Comparing the three curves, the program reproduces the analytical results reasonably well. The slight discrepancy in the elastic-plastic region, specifically for the 3-dimensional program, is attributed to the fact that the failure criterion does not take into consideration the stresses in the through-thickness direction. Also, the stresses were evaluated at the Gaussian integration points. Evaluating the conditions at the extreme fiber instead (where the longitudinal stresses would be approximately 12 percent higher) would lead to a lower yield load and thus a closer approximation of the analytical result. 67 - f -4@0.5 J.CH-A , --p - 3 - a - D 6@1 •4-4, 1 4 ® 0.5 1 4, 4@0.5 v v v 2-X (b) 6@1 4@0.5 Figure 4.9 - Finite Element Mesh for a) 2 and b) 3 Dimensional Analysis of Isotropic Cantilever Beam 4.5.2 5 Point Continuous Beam Analysis A further verification of the program entailed a comparison of the displacement and stress results obtained from a common commercial finite element program, ANSYS®, with the results obtained from COMAP. The purpose of this analysis is to confirm the accuracy of the linear elastic, orthotropic finite element formulation. The problem under consideration is that of a continuous beam illustrated in Figure 4.10. 1 I Y/ i V///A Figure 4.10-5 Point Continuous Beam As the beam is symmetric with respect to the longitudinal axis, only one half of the physical problem is modeled. For both programs, COMAP and ANSYS* in 2 and 3 dimensions, the geometrical parameters and finite element mesh are as shown in Figure 4.11. The applied load was 68 chosen such that the resultant internal stresses remained within the elastic range. A total load of approximately 360 N/mm width was applied, as shown. Linear 4-node, or 8-node elements were used, together with a 2x2 (or 2x2x2) Gauss integration scheme for both programs. Plane of symmetry £ ») 100 m: b) Point 'O' Figure 4.11 - Finite Element Mesh for Continuous Beam Example: a) 2 dimensions b) 3 dimensions The deterministic elastic properties for the material are: Longitudinal Elastic Modulus: E x = 11 x 103 MPa Transverse Elastic Modulus: E y = 0.4 x 103 MPa Through-thickness Elastic Modulus: E z = 0.62 x 103 MPa (for 3 d. analysis only) Poisson's Ratios: v^ = 0.32; v^ = 0.012; v x z = 0.29; v z x = 0.01; v y z = 0.20 and v 2 y = 0.20 (for 3 d. analysis only) Shear Moduli: G x v = 0.7 x 103 MPa G x z = 0.76 x 103 MPa; G y z = 0.08 x 103 MPa (for 3 d. analysis only) 69 Poisson's ratios were calculated based on v x y and the Elastic moduli in accordance with Equation 3.10. The general behaviour of the beam according to each program was observed. The longitudinal stress (oj contours, as well as deflection behaviour are demonstrated in Figure 4.12. For both the 2-dimensional and 3-dimensional cases, the COMAP results agreed favorably with that for ANSYS*. The results for the specific point 'O ' (depicted in Figure 4.11) is provided in Table 4.3. Any discrepancies are minor and are attributed primarily to extrapolation of stresses from Gaussian points to nodes. Table 4.3 - Linear Elastic Orthotropic Comparison of COMAP with ANSYS* Program Displacement at •O* (mm) Stress at 'O ' (MPa) 6X ° y 2 - Dimensional C O M A P -0.17 -1.8 - -14.2 -2.5 - 1.4 - -2 - Dimensional ANSYS* -0.18 -1.8 - -14.0 -2.7 - 1.3 - -3 - Dimensional C O M A P -0.18 -1.9 -0.07 -13.3 -2.9 -0.02 1.2 -0.01 0.47 3 - Dimensional ANSYS* -0.18 -2.0 -0.06 -14.5 -4.3 -0.4 1.3 -0.02 0.47 70 AKSYS 5.3 OCT 11 1997 U:29i09 PLOT NO. 1 MODAL SOLUTION STEP=1 SOB =10 TIMK.1 SZ (AVG) RSY8-0 IXX -2.057 SMK —15.887 SMX -14.675 1 HFQR RFOR -15.887 -12.491 -9 .096 - 5 . 7 - 2 . 3 0 * I. 092 4.487 7.883 I I . 279 14.875 Figure 4.12 - 3 Dimensional A N S Y S * Solution to Continuous Beam Example 4.5.3 Angle Ply Laminate Analyses The program C O M A P is formulated to analyze laminates with varying ply grain angles. To verify this aspect of our program, we use for comparison a program called PC-Laminate, version 1. This program was written by D.W. Radford from T U T E Technology Services. PC-Laminate employs simple mechanics with classical lamination theory (section 4.2.2.3) to calculate linear elastic mechanical responses of laminates with varying ply angles and also provides an estimate for the load to cause the initial failure of the laminate (first ply failure) under in-plane loads using the basic Tsai-W u criterion (Equation 2.3) with the interaction F 1 2 parameter equal to zero. 71 For the present analysis, [± 15]s and [±30] s angle-ply laminates were investigated. The stacking sequence of these laminates is illustrated in Figure 4.13 where, in our case, 0 = 15° or 30°. Figure 4.13 - Exploded View of Symmetric Angle-ply Laminate For both programs, the deterministic elastic and strength properties of the material were input as: Parallel-to-grain Elastic Modulus : Ej = 11 x 103 MPa Perpendicular-to-grain Elastic Modulus : E 2 = 0.4 x 103 MPa Through-thickness Elastic Modulus : E 3 = 0.62 x 103 MPa (for 3-d analysis only) Poisson's Ratios: v 1 2 = 0.32; v 2 1 = 0.012; (for 2-d analysis) v 1 3 = 0.29; v 3 1 = 0.01; v 2 3 = 0.20 and v 3 2 = 0.20 (for 3-d analysis only) Shear Moduli: G 1 2 = 0.7 x 103 MPa (for 2-d analysis) G I 3 = 0.76 x 103 MPa; G 2 3 = 0.08 x 103 MPa (for 3-d analysis only) Parallel-to-grain Tensile Strength: Xj = 80 MPa Perpendicular-to-grain Tensile Strength: Y t = 5 MPa Parallel-to-grain Compressive Strength: X c = 60 MPa 72 Perpendicular-to-grain Compressive Strength: In-Plane Shear Strength: Interaction Parameter: Y c = 15 MPa S = 6 MPa F 1 2 =0.0 MPa"2 Both angle ply configurations, [± 15]s and [±30] s ) were discretized in the same manner, illustrated in Figure 4.14. In the 2 dimensional model, the plane, 4 node elements were comprised of 4 layers each representing the plys. For each layer, a 2x2 Gauss quadrature rule was used resulting in 16 different stress points per element. The load was applied proportionately at the 5 end nodes, as shown. For the 3 dimensional model, 8 node brick elements were used. The width of each element coincided with the thickness of the ply. The load was spread over the cross section as shown. Figure 4.14 - Finite Element Mesh for Angle Ply Laminates a) 2 dimensions b) 3 dimensions Table 4.4 outlines the elastic stiffness and first ply failure results for all analyses. The values agree very well. The only difference worthy of mention is between the first ply failure load result for C O M A P (maximum 60.1 MPa) and PC-Laminate (55.9 MPa) with a maximum percentage error of 73 7%. This is due to the fact that, unlike the result by PC-Laminate, the ultimate stress computed C O M A P includes some stress redistribution prior to failure. Table 4.4 - Laminate Stiffness and First Ply Failure Comparison of P C L A M and COMAP Program [± 15]s Angle - Ply Laminate [ ± 3 0 ] s Angle - Ply Laminate E x Ultimate stress, ox E x Ultimate stress, ox (MPa) (MPa) (MPa) (MPa) P C L A M 9080 55.9 4870 23.1 2 - D C O M A P 9083 60.1 4866 23.6 3 - D C O M A P 9072 59.3 4855 23.2 5 - Strand Property Database 5.1 INTRODUCTION Evident from the discussion of the yield criterion in Chapter 3, the constitutive properties of the strands (in tension, compression and shear) are required as input to the model. Due to the unique failure characteristics of wood in response to these loadings and with respect to the principal material directions, different methodologies were considered to best estimate each property. Both experimental and analytical approaches were taken. The purpose of this chapter is to both describe fabrication and preparation of test specimens and also delineate the methods used in acquiring the constitutive properties of the strands. As mentioned in the introduction, one potential function of the model developed in this thesis would be to gauge the influence of changing the host material on the final product strength. For example, one could investigate the capacity of a PSL product made from a different wood species. Conceivably, one could perform simple constitutive tests, as outlined in this chapter, on any wood species and use the results as input to the program to obtain an estimate of the properties of the 'new' product. 5.2 M A T E R I A L A total of 100 - 3mm x 1220mm x 2440mm sheets of Coastal Douglas-fir veneers were generously supplied by Trus Joist MacMillan Ltd. The veneers were grouped in equal amounts of heartwood and sapwood and were representative of those used in the making of Parallam*, PSL. Veneers chosen for use were deliberately chosen free of large defects such as knots or localized slope of grain. 75 All specimens described in the following sections were prepared and fabricated by the author from these veneers. 5.3 EXPERIMENTAL MEASUREMENT OF STRAND PROPERTIES Standard test methods have been established for testing of wood in tension and compression by the American Society for Testing and Materials (ASTM). Although no specific tests are identified by these standards for testing of individual wood strands (specimens were nominally 19mm in width to emulate that used in PSL), ASTM standards were used as general guidelines. Specimens were conditioned to the ambient laboratory environment as opposed to being conditioned in a controlled humidity chamber due to a lack of this equipment at one of the test facilities. A l l compression specimens were tested in ambient environmental conditions of average 21.9°C ± 1.3°C and 30.2% ± 5.8% relative humidity. This correlates to an average of approximately 6.3% + 1% moisture content (Siau, 1984). For all other cases, specimen moisture content was monitored by use of the oven dry method (ASTM D2016) (mean = 7.7%, standard deviation = 0.9%). 5.3.1 Compression Tests 5.3.1.1 Material Preparation To overcome practical difficulties with testing equipment when dealing with the potential of buckling of thin specimens (such as working with an unpractically short gauge length to satisfy short column requirements), the compression tests were performed on 6-ply laminated specimens. These specimens (as well as specimens for other configurations discussed in the next sections) were 76 cut from boards which were fabricated in the laboratory as follows: 3 x 610 x 610 mm2 Douglas-fir heartwood veneers were first conditioned to a moisture content of approximately 4 percent to avoid blowouts. Then, a phenol-formaldehyde resin (PF 355H - a thermosetting resin used in making commercial wood composites) was combined in accordance with a standard procedure with pure dried cob, flour, soda ash and water, and then applied to the veneers using a roll spreader. The roll spreader provided for a consistent method of applying a prescribed amount of resin to each veneer. Six sheets of veneer oriented with the face grain in the longitudinal direction (such that the lathe checks for the bottom three veneers faced and opposed that for the top three veneers to avoid cupping) were laminated and pressed using a 1600 K N capacity Pathex 760 x 760 mm2 hot-press at 150°C and 1.38 MPa pressure for 6 minutes. The same procedure was done using Douglas-fir softwood veneers. Compression 'parallel-to-grain' specimens (nominally 17 x 19 x 60 mm3) and compression 'perpendicular-to-grain' specimens (17 x 19 x 50 mm3) were cut from 8 separate boards - 4 boards of heartwood and 4 of softwood. The specimens were cut using extreme care to ensure right angles to minimize eccentricities in loading. Furthermore, in the center of the cross-section of each specimen, a small divot (diameter less than 1 mm) was drilled to align the specimen with the center of the load platens which, in turn, had a matching small point. Parallel-to-grain specimen length was determined in accordance with ASTM D198 for short columns with no lateral support. This requirement states that the length to least radius of gyration not exceed 17. 77 5.3.1.2 Test Method Compression tests were conducted at the Laboratoire de Mechanique et Technologie (LMT), Ecole Normale Superieure (ENS), Cachan - University Paris, France. Tests were performed on a servo- hydraulic Material Testing System (MTS) machine with capacity 50KN under deflection control mode. The setup was equipped with a Linear Variable Displacement Transducer (LVDT) to measure cross head displacement. Swivel F i § u r e 5 A " Compression Test Set-up bearings, both top and bottom, were used to prevent eccentric loading on the specimen as shown i n F i g u r e 5 . 1 . A n M T S 6 3 2 . c 2 0 e x t e n s o m e t e r w h i c h is typically used with smooth surface advanced composites, was employed for some tests but gave inconsistent results due to poor contact with the wood surface. This approach was therefore abandoned. Instead, elastic moduli were determined using the cross-head displacement. To remove error caused by the inherent flexure in the testing machine, a compression test was performed with steel plates only. The resulting displacement, attributable to the machine, was thereby deducted from the overall displacement in the calculation of the material stiffness. Machine displacement was found to account for up to 30 percent of total displacement. Cross sectional dimensions were measured for each specimen using calipers with an accuracy of ±0.005 mm. A uniform loading rate of 0.12 mm/min. was used to produce failure between 5 and 10 minutes in compliance with A S T M D198. 78 5.3.1.3 Experimental Results Typical stress-strain curves for the parallel-to-grain and perpendicular-to-grain specimens are given in Figures 5.2 and 5.3, respectively. The curves have been zero-adjusted to remove nonlinearites at the curve origin. This was done by shifting the curve in the x direction only, such that the projection of the linear portion of the stress-strain curve passes through the origin. These nonlinearities occur due to the settling of the specimen within the rotating spherical bearing blocks. Beyond the initial linear portion of the curve, however, there is a notably nonlinear material response, as is expected for wood in compression (Goodman and Bodig 1971; Maghsood et al. 1973). 90 -80 -70 -CO a. S 5 0 -CO cn CD 55 40 -30 -20 -10 -/ Specimen H1-4L -0.002 0.002 0.006 0.01 0.014 Strain (mm/mm) Figure 5.2 - Typical Stress-Strain Curve for Compression Parallel-to-Grain Specimens 79 Strain (mm/mm) Figure 5.3 - Typical Stress-Strain Curve for Compression Perpendicular-to-Grain Specimens In order to describe this mechanical behaviour for use in the computer model, the curve is simplified to a tri-linear approximation with 4 defining variables: the elastic modulus prior to yielding, the yield stress, the elasto-plastic tangent modulus beyond yield and the ultimate stress (denoted for a perpendicular-to-grain specimen: E 2 c ) Y c ) E 2 c ' , and Y c u , respectively). These values were determined for each specimen by equating the area under the experimental stress-strain curve (i.e. the strain energy stored in the specimen) to the area of the fitted tri-linear stress-strain curve as illustrated in Figure 5.4. Perfect plasticity beyond ultimate load (as opposed to strain softening) has been assumed. This assumption not only simplifies the analysis greatly, but also agrees with experimental results up to an acceptable strain level. 80 Specimen S3 - 12T -0.02 0.04 Strain (mm/mm) 0.06 0.08 0.1 Figure 5.4 - Tri-linear Approximation of Compression Stress-Strain Curve A comprehensive summary of all measured compressive properties is given in Table Descriptive statistics are provided for each set of data. Table 5.1 - Compression Properties of Strands Configuration Count Property Elastic Yield Tan gent Ultimate Modulus Stress Modulus Stress Mean C O V . Mean C O V . Mean C O V . Mean C O V . (GPa) (%) (MPa) (%) (GPa) (%) (MPa) (%) Heartwood Parallel 53 10.09 19.13 67.32 18.04 1.93 33.11 76.46 7.06 Perp. 54 0.49 15.23 15.37 11.75 0.11 35.05 18.19 9.28 Sapwood Parallel 51 11.72 27.65 65.16 10.14 2.47 36.56 73.61 9.64 Perp. 56 0.37 13.54 14.75 20.68 0.10 33.13 16.71 18.09 A single factor Analysis of Variance (ANOVA) was performed for the heartwood and sapwood results for each property. The results of this analysis is outlined in Tables 5.2 and 5.3. The two data sets broken down by descriptive property elastic modulus, yield stress, etc., for both parallel-to-grain and perpendicular-to-grain samples, were found to be significantly different at a 0.05 percent level of significance. Consequently, the two sets of data could not be combined. It was decided arbitrarily that further development of the strand database and thus the model, would focus solely on the heartwood material. Table 5.2 - Analysis of Variance Results for Heartwood and Sapwood for Parallel-to-grain Compression Dependent Source of Variation Degree of Sum of Mean F To.05 Variable Freedom Squares Square Value Elastic M o d u l u s Between groups 1 14434405.5 14434405.5 2.1 3.9 W i t h i n groups 97 676023834.2 6969317.9 To ta l 98 690458239.6 Y i e l d Stress Between groups 1 116.1 116.1 3.2 3.9 W i t h i n groups 97 3534.4 36.4 Tota l 98 3651.5 Tangent M o d u l u s Between groups 1 7228897.3 7228897.3 11.91- 3.9 W i t h i n groups 97 58960353 607838.7 To ta l 98 66189250 Ul t imate Stress Between groups 1 184.3 184.3 4.6 1 ' 3.9 W i t h i n groups 97 3858.3 39.8 Tota l 98 4042.7 significantly different at a 0.05 percent level of significance 82 Table 5.3 - Analysis of Variance Results for Heartwood and Sapwood for Perpendicular-to-grain Compression Dependent Source of Variation Degree of Sum of Squares Mean F F0.05 Variable Freedom Square Value Elastic Modulus Between groups 1 330105.8 330105.8 79.51' 3.9 Within groups 100 415164.6 4151.6 Total 101 745270.4 Yield Stress Between groups 1 9.5 9.5 1.6 3.9 Within groups 100 607.5 6.1 Total 101 617.0 Tangent Modulus Between groups 1 19.1 19.1 0.02 3.9 Within groups 100 127398.8 1274.0 Total 101 127418.0 Ultimate Stress Between groups 1 51.0 51.0 8.81- 3.9 Within groups 101 577.4 5.8 Total 102 628.3 1 significantly different at a 0.05 percent level of significance The nonlinear compressive behaviour of the strands is randomly generated within the model based on the descriptive statistics of the four defining variables which are appropriately correlated. Each variable is simulated according to a bivariate standard normal distribution (Wang and Lam, 1998): bi = H b i + E L i j Z j (i = l,2,3,4) (5.1) i=i where: (for parallel to grain properties, for example) bx = E 1 C , b2 = X c , b3 = E 1 C ' , b4 = X c u , is the respective variable mean, L;j is the lower triangle of the correlation matrix of the variables, derived through Cholesky decomposition, and z ; is an independent standard normal random variable. The correlation matrix and corresponding lower triangle matrix are given for both parallel-to-grain and 83 perpendicular-to-grain properties in Tables 5.4 through 5.7. Table 5.4 - Correlation Matrix for Compression Parallel-to-Grain Specimens Property Elastic Modulus Yield Stress Tangent Modulus Ultimate Stress Elastic Modulus 1.00 0.48 0.20 0.55 Yield Stress 0.48 1.00 -0.21 0.84 Tangent Modulus 0.20 -0.21 1.00 0.14 Ultimate Stress 0.55 0.84 0.14 1.00 Table 5.5 - Lower Triangle of Correlation Matrix for Compression Parallel-to-Grain Specimens Property Elastic Modulus Yield Stress Tangent Modulus Ultimate Stress Elastic Modulus 1.00 0.00 0.00 0.00 Yield Stress 0.48 0.88 0.00 0.00 Tangent Modulus 0.20 -0.35 0.91 0.00 Ultimate Stress 0.55 0.65 0.28 0.43 Table 5.6 - Correlation Matrix for Compression Perpendicular-to-Grain Specimens Property Elastic Modulus Yield Stress Tangent Modulus Ultimate Stress Elastic Modulus 1.00 0.13 0.10 0.37 Yield Stress 0.13 1.00 -0.42 0.81 Tangent Modulus 0.10 -0.42 1.00 0.07 Ultimate Stress 0.37 0.81 0.07 1.00 84 Table 5.7 - Lower Triangle of Correlation Matrix for Compression Perpendicular-to-Grain Specimens Property Elastic Modulus Yield Stress Tangent Modulus Ultimate Stress Elastic Modulus 1.00 0.00 0.00 0.00 Yield Stress 0.13 0.99 0.00 0.00 Tangent Modulus 0.11 -0.43 0.90 0.00 Ultimate Stress 0.37 0.77 0.40 0.34 Figures 5.5 through 5.12 show the experimental data points in terms of cumulative probability distributions of the four defining variables together with fitted bivariate normal distributions. Through visual inspection, the experimental results are very well represented by the simulated data. 4000 6000 8000 10000 12000 14000 16000 Elastic Modulus (MPa) Figure 5.5 - Compression Parallel-to-Grain Elastic Modulus Data 85 • Experimental data Bivariate normal distribution -+- + 500 1000 1500 2000 2500 Tangent Modulus (MPa) 3000 3500 Figure 5.6 - Compression Parallel-to-Grain Tangent Modulus Data 0 A I 1 1 1 1 1 1 50 55 60 65 70 75 80 85 Yield Stress (MPa) Figure 5.7 - Compression Parallel-to-Grain Yield Stress Data Ultimate Stress (MPa) Figure 5.8 - Compression Parallel-to-Grain Ultimate Stress Data 200 300 400 500 600 700 Elastic Modulus (MPa) Figure 5.9 - Compression Perpendicular-to-Grain Elastic Modulus Data 87 Tangent Modulus (MPa) Figure 5.10 - Compression Perpendicular-to-Grain Tangent Modulus Data Ultimate Stress (MPa) Figure 5.12 - Compression Perpendicular-to-Grain Ultimate Stress Data 5.3.2 Tension tests 5.3.2.1 Material Preparation For parallel-to-grain tests, individual strands (3 x 19 x 50 mm3) were cut from 8 separate sheets of Douglas-fir heartwood veneer in the longitudinal orientation. For perpendicular to grain tests, in order to avoid breakage during handling, specimens were cut from 6-ply laminated veneer boards (fabricated as described in section 5.3.1). These specimens (17 x 19 x 152 mm3) were sawn from 4 separate boards using a fine carbide-tipped cut-off saw blade to minimize chipping. 5.3.2.2 Test Method Al l tension tests were performed at the University of British Columbia on a 250KN MTS machine together with MTS mechanical wedge action grips and an extensometer for calculation of elastic modulus, as depicted in Figure 5.13. 89 25.4mm gauge length_£^j MTS extensometer MTS wedge-action grip " tension specimen Figure 5.13 - Tension Test Set-up Cross sectional dimensions were measured using calipers. The loading rate was set at 0.64 mm/minute and to produce failure in no less than 5 minutes and no greater than 10 minutes in compliance with A S T M D198. Ultimate load and description of failure were recorded. Although rare, specimens which failed in the grip region were discarded. 5.3.2.3 Experimental Results A comprehensive summary of all measured mechanical properties is given in Table 5.8. Typical stress-strain curves for the parallel-to-grain andperpendicular-to-grain specimens are given in Figures 5.14 and 5.15, respectively. The strain in these graphs is calculated using crosshead movement divided by gauge length. It is noted, however, that the crosshead movement includes inherent grip movement. Referencing Figure 5.13, as the load increases, the specimen crushes around the grip teeth, which precipitates vertical displacement within the wedge action grips. This occurs prominently at initial load (leading to nonlinearities), and continues throughout loading to a lesser degree as the specimen becomes more dense. Consequently, the results have been zero-adjusted to remove these nonlinearities at the curve origin per section 5.3.1.3. Also, as a result, the graphs 90 appear more compliant and do not represent a true account of tensile stiffness. Tensile stiffness, instead, was determined through use of the extensometer. The extensometer measures only local displacement and must be removed prior to ultimate failure of the specimen to avoid instrument damage. As such, it was not an alternative to calculate strain in Figures 5.14 and 5.15. The graphs demonstrate the linear elastic behaviour of wood when loaded in tension. In view of this, two variables only were used to define and regenerate the tensile behaviour in each direction for the model: E l t , Xj, and E 2 t , Y t - elastic modulus and ultimate stress in the parallel and perpendicular direction, respectively. Elastic modulus was determined for each specimen through linear regression of the stress-strain curve where strain was measured using the 25.4 mm gauge extensometer. Figures 5.16 through 5.19 show the data points in terms of cumulative probability distributions of elastic modulus and ultimate stress together with fitted log-normal distributions. Table 5.8 - Tension Properties of Strands Configuration Count Property Elastic Modulus Ultimate Stress Mean (GPa) COV (%) Mean (MPa) COV (%) Parallel-to-grain 36 15.46 30.50 68.77 26.68 Perpendicular-to-grain 45 0.09 24.40 1.91 18.31 91 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 Strain (mm/mm) Figure 5.14 - Typical Parallel-to-Grain Tensile Behaviour 2 T -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 Strain (mm/mm) Figure 5.15 - Typical Perpendicular-to-Grain Tensile Behaviour 1 T Ultimate Stress (MPa) Figure 5.16 - Lognormal fit of Parallel-to-Grain Tensile Strength Data Ultimate Stress (MPa) Figure 5.18 - Lognormal Fit of Perpendicular-to-Grain Tensile Strength Data 1 T 20 40 60 80 100 120 140 160 Elastic Modulus (MPa) Figure 5.19 - Lognormal Fit of Perpendicular-to-Grain Tensile Modulus Data 94 The tensile strengths (both parallel and perpendicular-to-grain) reflect only the strength of the corresponding tested volume. Prior to implementing the values into the finite element code, adjustments must first be made for size. 5.3.3 Strength Modification due to Size Effect It has long been recognized for brittle materials that large members tend to display lower strengths than smaller members when subjected to the same environmental and loading conditions. The phenomenon has come to be known as 'size effect' and has implications on the tensile strength data discussed in the previous section. Numerous researchers have acknowledged and addressed size effect for many different materials (Barrett et al., 1975; Madsen, 1990; Sharp and Suddarth, 1991; Zweben, 1994). Each of these studies espoused the use of Weibull weakest-link theory to quantify size effect. As such, this theory was adopted for this study. 95 5.3.3.1 Weibull Weakest-Link Theory Weibull (1939) proposed that a variation in strength properties existed due to the statistical probability of a strength controlling defect occurring in a given volume. He used the weakest link concept to show how strength of a 'perfectly brittle' material could be described by a specific cumulative distribution function. Weibull's theory enabled the prediction of the probability of failure of a homogeneous isotropic material at a given volume according to the following: F(T) = 1 _ e ~ i l b na v M 1 ™ ! L dV where : F (x) = probability of failure x = material strength x m i n = minimum material strength (location parameter) m = scale parameter P = shape parameter (xm i n , m, and P are material constants and V Q is a reference volume) Commonly, T m ; n is accepted as being equal to zero, simplifying the problem, and the resulting formulation is called a two-parameter Weibull distribution. This distribution may be used to explain the strength difference observed in material bodies subjected to the same stress distribution but differing in volume as follows: 96 It is argued that a larger member has a higher probability of containing a larger flaw (or weaker zone) than does the smaller member and thus has a lower strength. In general, the strength of a volume of material at a given probability of failure is predicted given the strength and shape parameter of a common material when both are subjected to the same stress distribution. The solution is found by equating the probability of failure of one volume, V[ (corresponding to Tj), to that of another volume, V 2 (corresponding to T2), which yields: This simplifies to I v , T i P d V i = l v 2 T 2 P d V 2 (5-4) If it may then be assumed that the stresses in both volumes are uniform throughout the cross sections, then Equation 5.4 becomes: x , \ = T 2 P V 2 (5.5) Or, given that the cross section dimensions are equal: T /L , = x 2 P L 2 (5.6) As illustrated in Figures 5.14 and 5.15, wood displays brittle behaviour under tensile loading. Based on this, a size effect was assumed to exist and consequentially, the average values of the tensile strengths (both parallel and perpendicular-to-grain) were adjusted from the test volume to the representative size in the finite element analysis according to either Equation 5.5 or 5.6. 97 5.4 A N A L Y T I C A L L Y D E R I V E D P R O P E R T I E S Various approaches have been taken in the past for determining shear strength and stiffness of wood - the more popular being the ASTM shear block method for small clear specimens described in ASTM D143. This method, however, has limitations in that the test setup creates a complex distribution of stresses in the specimen and hence a questionable estimation of pure shear strength. Also, an added challenge lays in finding the pure shear characteristics of a wood strand (ie. from practical difficulties in testing the strand due to its geometry). At the same time there are complications, as outlined in Section 2.2.1.2, with establishing the interaction parameter of the Tsai-Wu criterion. An alternate method to obtain all three required variables (shear strength, S; modulus of rigidity, G; and the interaction parameter, F 1 2), based loosely on the unidirectional polymer matrix composite standard, ASTM D3518 / D3518M, has been developed for this study. This standard employs a [±45] s angle-ply laminate subjected to uniaxial tension. It can be shown that the lamina shear stress (o 1 2) is related to the uniaxial tensile stress on the laminate (o^ by (5.7) and hence an estimate of in-plane shear strength (S) can be calculated from the ultimate tensile stress. Although the standard advocates this method for shear strength estimation, it is noted that the lamina is still not in a state of pure shear, as normal stress components are present in the transformed coordinate system. As there is a known multiaxial state of stress present, it would be more logical to employ a strength criterion. Accordingly, the following approach has been developed to estimate the mean and standard deviation of the shear properties and the interaction parameter of wood strands. 98 5.4.1 Minimization Technique The statistical parameters of the variables were approximated simultaneously following a minimization technique originally developed for the Foschi-Yao damage accumulation model (Foschi.and Yao, 1986). The procedure entails a nonlinear least square minimization of error (using a quasi-Newton approach) between predicted and experimental compression strengths of a [ + 15]s angle-ply laminate. The compressive strengths of the laminates were simulated using the finite element model described in Chapter 4 together with Monte Carlo simulations. The predicted laminate strength is a function of the mean and standard deviation of the three unknown variables (S, G, F 1 2). Thus, a vector of six unknowns in total A = | ( J , S ,plG , \lP , S D S , S D G , S D F j is solved in the minimization process. An outline of the procedure follows: 1. Initial estimates for the mean and standard deviation of the three variables were provided for the model. A normal distribution was deemed appropriate for F 1 2 , as this enables either positive or negative values, reflecting this characteristic of the parameter. Lognormal distributions were chosen to represent the shear properties. 2. Strength was calculated for a sample of 300 replications. A larger sample could also have been used; however, the number of replications was needed to address computer time restrictions. For each replication the entire load displacement behaviour was modelled in order to establish ultimate stress. For the three dimensional model, each replication took between 20 to 30 minutes using a MS Windows * based, 233 MHz, Pentium II system. 3. These predicted strengths were ranked (ie. sorted in ascending order and given the associated probability of failure). 4. A residual function O was calculated as: 99 p ( T" 1 -2 exp (5.8) where: P denotes the number of probability levels for consideration and the superscripts pred and exp refer to the predicted and experimental laminate strengths, respectively. The predicted strength was determined for the same i t h probability level as the experimental strength. 5. As required by the quasi-Newton method, the gradient of the residual function with respect to f dO dO } ^ the unknowns V O = < ,... > was calculated. This was done numerically, based [DA, 3A6 J on a perturbation process as follows: each unknown A ; was perturbed by a positive increment 0.001A;. Then $ was recalculated, designated 0 + . The same was done for the same negative increment producing The gradient was then calculated as: do o + - o ~ T T - = , A A (1 = 1>-6) (5-9) oA- 2 -AA-6. New, adjusted, statistical parameters replaced the initial estimated values and the residual function was re-evaluated. This procedure was repeated until the difference between residual function values for subsequent iterations satisfied a set tolerance of l.OxlO"3. It is noted that the [+ 15]s angle-ply laminate was determined to be the most suitable configuration for determining both shear and the interaction parameter together. The reason stems from a finding in Clouston et al. (1998). Although the 45° configuration would produce the largest in-plane shear stresses (which would be beneficial in estimating shear strength), it was shown through a sensitivity analysis that the data from 15° off-axis tests had more tolerance for inaccuracies in experimental strength results and consequently were more reliable than for other angles (30°, 45°, or 60°) for 100 calculating the interaction parameter, Fu. 5.4.2 Experimental Tests 5.4.2.1 Material Preparation Inherent in the preceding minimization technique are the experimental compressive strengths of the [±15] s angle-ply laminates. These laminates were manufactured for testing. Four sheets of conditioned 3mm x 610mm x 610mm Douglas-fir heartwood veneer were laminated with a phenol-formaldehyde resin (PF 355H) based adhesive and pressed using a 1600 K N capacity Pathex 760 x 760 mm2 hot-press at 150°C and 1.38 MPa pressure for 6 minutes. The sheets were oriented in a symmetrical layout with the outside veneer angle at +15° to the longitudinal direction and the inside two veneers at an angle of -15° to the longitudinal direction. Naturally, the grain angle of each sheet was not uniform, so a visual average was used. The 15° angle was established using a protractor aligned with the average grain. The lathe checks of the bottom two veneers faced and opposed that of the top two veneers to avoid cupping. The specimens (nominally 11 x 19 x 40 mm3) were cut from 4 separate boards. They were prepared in the same manner as described in section 5.3.1.1. (ie. with right angles and a small divot at center cross section to minimize eccentricities). A length of 40 mm was established in accordance with ASTM D198 for short columns with no lateral support. The specimens were tested in compression using the same test method described in Section 5.3.1.2. 101 5.4.2.2 Experimental Results The mechanical properties of the [+ 15]s angle-ply laminates are summarized in Table 5.9. Figure 5.20 demonstrates the range of results for the stress-strain behaviour showing 14 of the total 31 curves including the strongest, weakest, stiffest, least stiff and the computed average. The results have been zero-adjusted to remove nonlinearites at the curve origin per section 5.3.1.3. The behaviour is predominantly ductile, although ultimate failure was generally abrupt. Visual inspection of the failed specimens (Figure 5.21) showed a combination of tension perpendicular-to-grain and/or in-plane shear failure. Table 5.9 - Compression Properties of [+ 15]s Angle-ply Laminate Statistics Count Elastic Modulus Ultimate Stress Average (MPa) C O V (%) 31 6505.27 26.28 51.65 15.39 102 70 t 0.025 Strain (mm/mm) Figure 5.20 - Select Sample of Stress-Strain Diagrams for [±15] s Angle-ply Laminates it f Figure 5.21 - [± 15], Angle-ply Laminate Compression Specimen 103 5.4.3 Minimization Results 5.4.3.1 Minimization using the 2 Dimensional Model The results of the minimization technique using the 2-dimensional model to calculate compression strength of the [±15] s angle-ply laminate, are summarized in Table 5.10. (Details of the 2 dimensional model were outlined in section 4.2 and a thorough description of the application of the program is deferred to Chapter 6). Table 5.10 - Results of Nonlinear Minimization using 2-dimensional Model Statistics [± 15]5 Compression S G Experimental (MPa) Simulated (MPa) (MPa) (MPa) (MPa'2) Mean 51.57 51.65 5.99 232.8 5.14X10"0 4 C O V (%) 14.42 15.39 11.7 17.76 71.73 The model predicts the average ultimate compression strength of a [± 151 angle-ply laminate to be 51.65 MPa. This is very close to the experimental result, 51.57 MPa, as would be expected since the parameters were fitted using this experimental data. The predicted variability is also very accurate, with predicted and experimental coefficients of variations of 15.39 percent and 14.42 percent, respectively. The average values and variability obtained for S, G and F 1 2 appear reasonable. Shear and shear modulus for Douglas-fir clear wood have been reported by Bodig and Jayne (1993) as 6.4 MPa and approximately (using a 14:1 ratio for E t:G) 800 MPa. The difference between the calculated and published results for the latter is attributed to the thin nature of strands as opposed to solid wood 104 specimens. The interaction parameter F 1 2 is more difficult to evaluate as there is no directly comparable data in the literature for wood strands. As such, one can evaluate it based on its conformity to a deterministic evaluation of the stability bounds (Equation 2.8) as follows: From Tables 5.1 and 5.8, the average strength values are: Xc=67.32 MPa ; Y c = 15.37 MPa; X, = 68.77 MPa ; Y t = 1.91 MPa The tensile strengths, X, and Y t , are first adjusted for size effect. For tension parallel-to-grain, the shape parameter of the 2-parameter Weibull distribution is calculated from maximum likelihood approach as P = 4.23 . Strength is adjusted from the length of the test specimen (LJ to the length of one strand in the prediction specimen (L^) using Equation 5.6 as follows: Given: X t l = 68.77 MPa = 50.8 mm L 2 =40.0 mm (3 = 4.23 For tension perpendicular-to-grain, k = 6.66 and strength is adjusted from the volume of the test specimen (V^ to the volume surrounding one integration point (V^) using Equation 5.5: then: = 68.77 40.07 50.8 H 2 3 = 72.8 MPa 105 Given: Y t l = 1.91 MPa V r = 19.09x17.35x152.0 = 50,344.2 mm 3 V 2 =2.375x5.0x2.55 = 30.28 mm 3 p = 6.66 then: X t 2 - X . , ~ V, P [ 50344.2 V - 6 6 = 1.91 =5.82 MPa 30.28 Using Equation (2.8), the stability bounds are found to be ± 1.51xlO"3MPa2. The result of 5. HxlO"4 MPa"2 falls within these bounds. As a visual reference, the cumulative probabilities of the 31 experimental ultimate strengths and 500 simulated strengths using the finalized statistical parameters for S, G and F 1 2 have been plotted in Figure 5.22. As expected, the distributions are very close. A log-normal distribution to the experimental data was also plotted for reference. Comparing the three curves, the experimental data tends to deviate slightly at the lower end of the curve, possibly due to manufacturing or test imperfections. Overall, however, the model predictions fit the experimental data extremely well. 106 Ultimate Stress (MPa) Figure 5.22 - Cumulative Probability Distribution of [+ 15]s Laminate in Compression (using 2-D Finite Element Model) 500 stress-strain curves for the [±15] s compression laminates were computer generated and compared with the experimental curves. The results are given in Figure 5.23. For clarity, only 5 curves of each (experimental and simulated) have been shown, which represent the average as well as the upper and lower bound for both stiffness and strength. The average curves were obtained by sequentially averaging stress along constant lines of strain. It is noted that the average curves serve only as a reference and is not indicative of material response, per se ; particularly within the higher strain range where the average values are based on fewer points because of varying total strain range. The simulated curves show the loading increments used in the model. 107 to 50 Q . I H ^ ' Experiment Average •Simulated Average - Experimental Simulated -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 -10 Strain (mm/mm) Figure 5.23 - Stress-Strain Curves of [± 15]s Angle-ply Laminate in Compression (using 2-d Finite Element Model) The simulated curves clearly lay within the experimental bounds. Furthermore, they capture the nonlinear behaviour of the [± 15]s laminates. 5.4.3.2 Minimization using 3-Dimensional Model The results of the minimization technique using the 3-dimensional model (described in section 4.3), are summarized in Table 5.11. 108 Table 5.11 - Result of Nonlinear Minimization using_3-dimensiona! Statistics [±15] s Compression S G F« Experimental (MPa) Simulated (MPa) (MPa) (MPa) (MPa2) Average 51.57 51.28 11.4 392.2 1.067xl0-°3 C O V (%) 14.42 13.84 21.8 21.7 8.78 F.E. Model Again, the average values and variability obtained for S, G and F 1 2 appear reasonable. The shear strength, in this case however, is twice that found using the 2 dimensional model. It should be understood that, although they compare well with published data, these values serve to calibrate the model and are expected to produce different results for different models. The interaction parameter is also slightly higher, although still within the deterministic bounds. The cumulative probabilities of the 31 experimental ultimate strengths and 500 simulated strengths using the 3-D model have been plotted in Figure 5.24. As in the 2-D model the experimental distribution is well predicted. 109 30.0 35.0 40.0 45.0 50.0 55.0 60.0 65.0 70.0 Ultimate Stress (MPa) Figure 5.24 - Cumulative Probability Distribution of [± 15]s Laminate in Compression (using 3-D Finite Element Model) Strain (mm/mm) Figure 5.25 - Stress-Strain Curves of [± 15]s Angle-ply Laminate in Compression (using 3-d Finite Element Model) Again, 500 stress-strain curves for the [± 15]s compression laminates were computer generated - this time using the 3 dimensional model, and compared with the experimental curves. The results are given in Figure 5.25. It is noted that the programs are adept at identifying mode of failure. Although precise initial mode of failure is difficult to ascertain from experiments due to microscopic cracks, it would appear from the stress - strain curves above that ductile behaviour is prevalent prior to ultimate tensile collapse as is predicted by the simulated results. In both cases, yielding occurs throughout the specimen due to predominant compressive stresses until a combination of tension perpendicular-to-grain and shear stresses cause abrupt brittle failure. Ill 6 - Model Verification 6.1 I N T R O D U C T I O N Having obtained all necessary strength parameters described in the previous chapter, the model is now in a complete form. Both the 2 dimensional and 3 dimensional program can be utilized to predict the nonlinear load - deflection behaviour of structural laminates up to and beyond failure. Within the objectives of this thesis, however, it remains to determine its accuracy. It is the purpose of this chapter to investigate the ability of the present model to reproduce experimental findings for a variety of combinations of geometrical and loading configurations on selected laminates. A total of eighteen categories were investigated as outlined in Table 6.1. Table 6.1 - Categories of Model Verification Tests Laminate Compression Tension Bending [0/0/0], - - 2 - D, 3 - D [90/90/90]s - - 2 - D , 3 -D [±15i 2 - D , 3 - D 2 - D , 3 -D 2 - D , 3 -D [±301 2 - D , 3 - D 2 - D , 3 -D 2 - D , 3 - D Parallam* - 3-D 3 - D This chapter is broken into three sections according to loading regime investigated. We begin by looking into simple uniaxial analyses; compression and tension, and progress to bending behaviour investigation. A preliminary investigation into modeling Parallam*, PSL in tension and bending is provided thereafter. 112 6.2 C O M P R E S S I O N V E R I F I C A T I O N The compressive behaviour of the [±15] s angle-ply laminate was utilised for the minimization process in the foregoing chapter. Although the good results obtained from this process do attest to the accuracy of the program, it would be redundant to discuss the findings here. The compression verification analysis will therefore focus only on the results for the [ ± 30]s angle-ply laminates. 6.2.1 Numerical Analysis Two dimensional and 3 dimensional analyses were performed to simulate the stress - strain response of the [±30] s, angle-ply laminates subjected to compressive loading. The geometrical properties, finite element mesh and boundary conditions are given in Figure 6.1. The finite element mesh properties are outlined in Table 6.2. Figure 6.1 - Finite Element Mesh for [ ± 3 0 ] S Angle Ply Laminates in Compression: a) 2 dimensions b) 3 dimensions 113 Table 6.2 - Finite Element Mesh Properties of Angle Ply Laminate Uniaxial Analyses Finite Element Mesh Properties 2 Dimensional Model 3 Dimensional Model Number of nodal points 25 125 Number of degrees of freedom 50 375 Number of elements 16 64 This mesh size was chosen for all uniaxial tests following a preliminary parametric study which is discussed in Appendix B. The 2 dimensional model utilized 4-node elements. Referencing Figure 6.2, each element was comprised of 4 layers where each layer represented an individual ply with a given ply angle. Using the assumptions of the classical lamination theory, outlined in Section 4.2.2.3, no slip was assumed between plies and strains are established in terms of the mid-surface strain. Therefore, using a 2x2 Gauss quadrature rule, there were 4 Gaussian points per element to establish the strains and 16 different stress points (4 points x 4 layers) per element or a total of 256 points per specimen monitored by the yield/failure criterion throughout the entire stress path. Stress points Figure 6.2 - Explanatory Depiction of 2 Dimensional Element 114 The 3 dimensional model employed 8-node brick elements as outlined in Section 4.3. Using a 2x2x2 quadrature rule, there were 8 integration points per element or a total of 512 integration points per specimen. As a result, the 3 dimensional program involved approximately 8 times more computer time. As the post-failure behaviour was sought, the laminate was loaded by prescribing an increasing value of displacement (as opposed to load). For each increment in displacement, the nonlinearities in the equilibrium equations were resolved using the modified Newton-Raphson iterative procedure. This displacement approach avoids problems with global stiffness matrix singularities beyond peak load. To assist in the making of input files for COMAP, an MS Excel® Visual Basic* macro was created. A property input dialog was used to manage the rather large number of statistical parameters required for COMAP. To demonstrate the use of this for the present analysis, the material stiffness and strength parameters for input into the 2 dimensional program are summarized in 'Windows' format in Figure 6.3. 115 PROPERTIES 3 Deterministic Propert ies-DK | Laminae thickness |2.55 mm Cancel | Poisson's ratio 0. _ s t ra in angle |30 Parallel to gram Tension -Perpendicular-to-grain Tension iBSiplIllllilllllllllllllllllii^ mean st dev. st. dev. Elastic Modulus (EXt] 154G3.0 4716 2 M P a \- Elastic Modulus (EYt) 91 2 22.3 M P a Strength (Xt) |72 8 |194 M P a ',. Strength (Yt) 5 83 |1.05 M P a -Parallel-to-grain Compres -ion mean st de^ • Perpendicular-to-grain Compression t dev Elastic Modulus (EXc) |10090 0 19D0 0 M P a Elastic Modulus (EYc] 490 0 |74 6 M P a Yield Strength (Xcl G7 3 1 3 0 M P a I Yield Strength (Yc) 1 5 4 1 8 M P a Tangent Modulus (EXc'j |1926 0 639 0 M P a Tangent Modulus (EYc'J 1 1 0 0 113 6 M P a Ult Strength (XcuJ |76 5 |5 40 M P a j Ult. Strength |YcuJ 1 8 2 |1 69 M P a In plane Shear mean st dev. Interaction Parameter IF12] 5.14E-04 MPa-2 Shear Modulus (GJ 2328 |41 3 MPa-jf :t dev 3.69E-04| MPa-". Shear Strength (S) |5 99 |0 70 M P a |; Figure 6.3 - Input Parameters for [±30] s Angle Ply Laminate in Compression (2 Dimensional Program) It is noted that the tensile strengths given in Figure 6.3 have been adjusted from those given in Table 5.6. The reason for this is that prior to implementing the values into the finite element code, adjustments must first be made for size effect as described in Section 5.3.3. The perpendicular-to-grain tension strength was adjusted from representing that of the tested volume to that of the tributary area for one Gauss point. Longitudinal strength was adjusted however, from the experimental gauge length to the model gauge length of 40 mm. This is exactly the same calculation using the same values as was outlined in Section 5.4.3.1. for the [± 15]s laminate in compression and will not be repeated here. 116 The calibrated factors, S, G and F 1 2 , as well as the tensile strength Y t differ for the 3 dimensional program from Figure 6.3. These input parameters are presented in Table 6.3. The mean value of Y t is different as a result of size effect (ie. due to the smaller volume surrounding each Gaussian point for the 3 dimensional model). The coefficient of variation of Y t is assumed to be the same for both models. Table 6.3 - Input Properties of 3 Dimensional Program Differing from those of 2 Dimensional Program Property Mean Standard Deviation Shear Strength, S (MPa) 11.4 2.6 Shear Modulus, G 1 2 (MPa) 390.0 89.7 Interaction Parameter, F 1 2 (MPa2) 1.067X10"03 9.368xl005 Perpendicular-to-grain Tension Strength, Y t (MPa) 6.46 1.14 The elastic moduli associated with the through-thickness direction, E 3 , G 1 3 and G 2 3 , have been estimated through general moduli relationships cited in Bodig and Jayne (1993) as E 3 : E 2 - 1.6 : 1 G n : G„ : G„ * 10 : 9.4 : 1 The major Poisson's ratios in the out-of-plane dimension are estimated by approximate relationships for Douglas-fir from Bodig and Jayne (1993): v 1 3 /E! = 2.6 x 10"05 MPa v 2 3 /E 2 = 4.9 x 10"04 MPa and the reciprocal relationships of Equation 3.10 are invoked for the reciprocal identities. These elastic moduli are 100 percent correlated with the generated random value; for example, for each 117 random value chosen for E 2 , a new value for E 3 is calculated based solely on the given ratio. The program is formulated for stochastic analyses. Al l strength and stiffness properties were assumed to vary between layers, with the exception of tension perpendicular to grain strength which was instead regenerated for each integration point. This strength was treated differently because in addition to yielding good results, material behaviour in this configuration is believed to be more in keeping with ideal brittle fracture theory. In this configuration, the within member variation is less likely to be correlated as it is, for example, for tension parallel to grain (Lam and Varoglu 1990). The failure mechanism is more likely 'perfectly brittle', failing completely when fracture occurs at the weakest point. The material flaws are assumed to be randomly distributed. 6.2.2 Experimental Tests For comparison to the analytical data, [±30]s angle ply specimens (nominally 19 x 10.5 x 40 mm3) were fabricated, prepared and tested under uniaxial compression in exactly the same manner as the [+ 15]s laminates described in Section 5.4.2. 6.2.3 Results The results of both the numerical simulations and experiments are summarized in Table 6.4 and Figures 6.4 and 6.5. 118 Table 6.4 - Experimental and Simulated Data for [ + 30]s Angle Ply Laminates in Compression Statistic Experiment Simulation (Count = 500) (Count = 39) 2d (% error) 3d (% error) Elastic Modulus Mean (MPa) 3052 2273 (25.5) . 2742 (10.2) Elastic Modulus COV (%) 27.5 13.3 (51.6) 8.3 (69.8) Strength Mean (MPa) 20.3 19.5 (4.0) 21.1 (3.8) Strength COV (%) 11.1 9.6 (13.9) 13.8 (23.9) Ultimate Stress (MPa) Figure 6.4 - Cumulative Probability Distribution of [±30] s Laminate in Compression (2 - Dimensional Model) 119 Ultimate Stress (MPa) Figure 6.5 - Cumulative Probability Distribution of [±30] s Laminate in Compression (3 - Dimensional Model) As evident from Table 6.4, the average value of the ultimate strengths is well predicted (maximum percent error of 4.0 percent). The average initial stiffness is underestimated for the 2 dimensional program yet quite accurate for the 3 dimensional program. The difference in the two programs is due to the lower shear stiffness calibrated for the 2 dimensional than for the 3 dimensional program (232.8 MPa vs. 390.0 MPa, respectively). Referencing Figures 6.4 and 6.5, the upper region of the cumulative probability distribution curve is slightly under-predicted by the 2 dimensional model, yet over-predicted by the 3 dimensional model. This is partly a result of the number of stress points for each model. The 3 dimensional model, having many more points for stress redistribution after localized failure, seemed to be better able to capture nonlinear behaviour and as a result produced higher laminate capacities. This may also account for the higher variability of the 3 dimensional model. 120 Stress Redistribution To illustrate how both programs manage stress redistribution upon integration point failure, as discussed in Section 4.4.3., the parallel-to-grain stress (ox) at prescribed points within a specimen (shown in Figure 6.6) have been tracked up to and beyond peak load. The stress-strain diagrams for each point, together with the overall longitudinal laminate behaviour, is given in Figures 6.7 and 6.8 for the 2 dimensional and 3 dimensional programs, respectively. Overall laminate stress is computed as the total applied force divided by the cross-sectional area of the specimen. a) Endview c — Figure 6.6 - Location of Sample Points used to Show Stress Redistribution: a) 2 Dimensions b) 3 Dimensions The location of the points which are plotted are marked as small crosses in cross section in Figure 6.6. For the 2 dimensional program, a total of 16 stress points are tracked. These are the 4 stress points per 4 layers of one element. For the 3 dimensional program, a total of 8 points are plotted which are simply the 8 Gaussian points of one element (one layer thick only). 121 30 T 0 0.002 0.004 0.006 0.008 0.01 0.012 Strain (mm/mm) Figure 6.7 - Stress - Strain Diagram Depicting Stress Redistribution Upon Stress Point Failure (2 Dimensions) Strain (mm/mm) Figure 6.8 - Stress - Strain Diagram Depicting Stress Redistribution Upon Stress Point Failure (3 Dimensions) Figures 6.7 and 6.8 show points which have failed in a brittle manner (gradual declining curves) and subsequent reaction of some surrounding points to compensate for this loss of capacity. For each increment in displacement shown, the yield criterion and equilibrium requirements are satisfied through the iterative procedure outlined in Section 4.4.2. Figure 6.9 and 6.10 attest to both program's capabilities to simulate the full range of stress - strain curves. Following the same format as for Section 5.4.3.1., the curves shown demonstrate the strongest, weakest, stiffest, least stiff, and calculated curve average for both experimental and simulated data. 28 T -0.005 0 0.005 0.01 0.015 0.02 Strain (mm/mm) Figure 6.9 - Stress-Strain Curves of [±30] s Angle-ply Laminate in Compression (using 2 - Dimensional Model) 123 32 T Strain (mm/mm) Figure 6 . 1 0 - Stress-Strain Curves of [ ± 3 0 ] S Angle-ply Laminate in Compression (using 3 - Dimensional Model) 6.3 T E N S I O N V E R I F I C A T I O N The general performance of the models to predict laminate behaviour under uniaxial tensile loading is investigated in this section using 2 different geometrical configurations: [±15]s and [±30]s angle ply laminates. 124 6.3.1 Numerical Analyses The 2 dimensional and 3 dimensional computation for constitutive tensile behaviour of the [ + 15]s and [±30] s angle ply laminates follows very closely to that conducted for compressive behaviour discussed in the previous section. The geometrical properties, finite element mesh, and boundary conditions are the same as described in Figure 6.1 save direction of applied displacement and specimen length. The [+ 15]s laminates were 120 mm in length and the [±30] s laminates were 60 mm in length. A longer length was used for the [± 151 laminates to ensure wood fibres were not continuous from one test grip to the other potentially increasing observed experimental strength. Given these specimen lengths, the tensile strengths were adjusted for size effect, per section 5.3.3, prior to being used in the programs. The calculation is as follows: For tension parallel-to-grain, the shape parameter of the 2 - parameter Weibull distribution of the experimental data was calculated from a maximum likelihood approach to be P=4.23 . Strength was adjusted from the length of the test specimen (Lj) to the length of one strand in the prediction specimen (Lj) using Equation 5.6: Given: P = 4.23 X,! = 68.77 MPa L t = 50.8 mm L 2 ([± 15]s laminate) = 120 mm ; L 2 ([±30]s laminate) = 60 mm then for the [± 15]s laminate 125 X t 2 - X t l P ( 50.8 V- 2 3 = 68.77 = 56.1 MPa 120.0 and for the [ + 30]s laminate X t 2 - X t l P ( 50.8^ 4.23 = 68.77 =66.1 MPa V60.0J The result is valid for both 2 and 3 dimensional analyses. For tension perpendicular-to-grain, P = 6.66 and strength is adjusted from the volume of the test specimen (VJ to the volume surrounding one integration point (V^ (which is different for both programs) using Equation 5.5: Given: p = 6.66 Y t l = 1.91 MPa V! = 19.09 x 17.35 x 152.0 = 50344.2 mm3 For the 2 dimensional program: V 2 ([±151 laminate) = 2.375 x 15.0 x 2.55 = 90.6 mm 3 ; Y t 2 = 1.91-^50.44-2! 6 6 6 = 4.9 MPa 90.! V 2 ([±30]s laminate) = 2.375 x 7.5 x 2.55 = 45.4 mm 3 ; Y t 2 = 1.91 ( ^ J ^ ) = 55 MPa and for the 3 dimensional program: ( 50344.2"! 6.66 V 2 ([±15], laminate) = 2.375 x 15.0 x 1.27 = 45.4 m m 3 ; Y t 2 = >-9l{ ^ j = 5.5 MPa 126 V 2 ([±301 laminate) = 2.375 x 7.5 x 1.27 = 22.7 mm 3 ; Y, 2 = 1.91-( 50344.2) 6.66 22.7 J = 6.1 MPa Coefficient of variation (COV) is assumed to remain the same as found in the experiment. Hence, the input parameters for the 2 dimensional and 3 dimensional programs are as outlined in Table 6.5. The elastic moduli associated with the through-thickness direction, E 3 , G 1 3 and G 2 3 , have been estimated through relationships cited in Bodig and Jayne (1993) as explained in Section 6.2.1. 6.3.2 Experimental Tests Specimens were cut from boards that were made as described in Section 5.4.2.1. The tension specimens were nominally 19 x 11 xl20 mm3 ([±15]s laminates) and 19 x 10.5 x 60 mm3 ([±30]s laminates) with minimal dimensional variation (1.4 % maximum in the through thickness dimension). 127 Table 6.5 - Input Properties for [± 15]s and [±30] s Angle Ply Laminates in Tension Property 2 Dimensional 3 Dimensional Mean St. Dev. Mean St. Dev. Parallel-to-grain Elastic Modulus (EXl) (MPa) 15463 4716.2 15463 4716.2 Tension Strength QQ [±15] 5 (MPa) 56.1 14.9 56.1 14.9 Strength (XJ [ + 30]s (MPa) 66.1 17.6 66.1 17.6 Perpendicular-to-grain Tension Elastic Modulus (EYt) (MPa) Strength QQ [±15] s (MPa) 91.2 4.9 22.3 0.87 91.2 5.5 22.3 0.96 Strength OQ [±301 (MPa) 5.5 0.96 6.1 1.07 Parallel-to-grain Elastic Modulus (Ex<) (MPa) 10090 1930 10090 1930 Compression Yield Strength (XJ (MPa) 67.3 13 67.3 13 Tang Modulus (EXc') (MPa) 1926 639 1926 639 Ultimate Strength (Xc u) (MPa) 76.5 5.4 76.5 5.4 Perpendicular-to-grain Elastic Modulus (EY() (MPa) 490 74.6 490 74.6 Compression Yield Strength (YJ (MPa) 15.4 1.8 15.4 1.8 Tangent Modulus (EYc') (MPa) 110 38.6 110 38.6 Ultimate Strength (Yc u) (MPa) 18.2 1.7 18.2 1.7 Interaction Parameter F 1 2 (MPa2) 5.1xl004 3.7xl0-°4 l.lxlO-03 9.4xl0-°5 In plane Shear Elastic Modulus (G) (MPa) 232.8 41.3 392.2 85.1 Strength (S) (MPa) 5.99 0.7 11.4 2.6 Poisson's Ratio Vl2 0.32 - 0.32 -128 The specimens were tested using a 250 k N MTS mechanical wedge action grips as described for the procurement of strand tensile data in Section 5.3.2.2. The loading rate was set at 1.8 mm/min. to produce failure in 6 minutes, on average. Tensile failure for both [ ± 15]s and [ ± 30]s were typical of that shown in Figure 6.12, where perpendicular-to-grain tension and in-plane shear stresses governed failure. Glue failures were not observed. universal test machine equipped with MTS Figure 6.11 - Tension Test Set-up Figure 6.12 - Tensile Failure of [±30] s Angle Ply Laminate Specimen 6.3.3 Results The results of the experimental tests and simulations are summarized in Table 6.6. The tension grips, depicted in Figure 6.11, were of a wedge type such that as the load increased, the measured displacement of the specimen included inherent displacement within the grips. As a result, the load-displacement curve exhibits, incorrectly, a lower stiffness. Unfortunately, laminate stiffness was not measured using an extensometer and so cannot be compared to the simulated data. Regardless, 129 descriptive statistics are included for the simulated elastic stiffness in Table 6.6 for completeness. The 2 dimensional model is less stiff than the 3 dimensional model due to the lower shear stiffness value used as input. Table 6.6 - Experimental and Simulated'. Data for [ + 15]s and [ ± 30]s Angle Ply Laminates in Tension Config. Statistic Experiment {Count} Simulation (Count =500) 2 d (% error) 3 d (% error) [±151 Elastic Modulus Mean (MPa) Elastic Modulus COV (%) Strength Mean (MPa) Strength COV (%) 39.4 {39} 16.7 6724.1 11.2 39.1 (0.6) 20.7 (23.9) 9128.6 12.8 36.4 (7.6) 17.5 (4.8) [±30] s Elastic Modulus Mean (MPa) Elastic Modulus COV (%) Strength Mean (MPa) Strength COV (%) 21.7 {41} 9.0 2273.0 13.3 21.3 (1.8) 12.1 (34.4) 3132.7 10.85 21.4 (1.4) 13.2 (46.7) Cumulative probability distributions for the [± 15]s and [±30] s angle ply laminates in tension are plotted in Figures 6.13 through 6.16. Ultimate stress is predicted relatively well by both models for both configurations (maximum percent error of 7.6%). The larger prediction error of the [± 15]s laminate strength for the 3 dimensional case is attributed to the fact that the 3 dimensional model has more Gaussian points than the 2 dimensional model for potential failure. We consider the stress state of a Gaussian point at incipient failure for the [+ 15]s laminate. For purposes of explanation, a likely stress state would be {o l s o 2 ) o 4}T = {39, -0.5, 4.1}T. From this, it is seen that tension parallel-to-grain governs. If the point is weak in tension perpendicular-to-grain and the layer is weak in tension parallel-to-grain, then the point will fail immediately in a brittle manner. Being a brittle 130 failure, as soon as one point fails, surrounding points (particularly in the weak layer) are heavily stressed and also fail, leading to catastrophic failure. Since the 3 dimensional model has more points, there is a higher probability of having a point which is weak in tension perpendicular-to-grain, and hence, of having overall lower laminate strengths. It is noted that the size adjustment to Y t for the 3 dimensional case does negate this somewhat; however, one could speculate that the Weibull approach may not be wholly appropriate for multi-axial stress states. A future study could investigate this further. The [±30] s laminate, on the other hand, fails numerically in a slightly more ductile manner (see Figure 6.17). A typical stress state of a Gaussian point at failure for the [±30] s laminate is {o,, o2, o 4}T = {22, -1.1, 6.2}T. In this case, the conditions of Table 4.1 may not be immediately met and the point may initially yield. Referencing Figure 3.1, when a point yields, it maintains strength and does not put the same burden on surrounding points as would a brittle failure. Thus, a weak point does not have the same impact on a [ ± 30]s laminate as for the [ ± 15]s laminate. For this reason, the 3 dimensional and 2 dimensional models are closer in prediction for the [ + 30]s laminate. The coefficient of variation of the [ + 301 laminate strength is somewhat high for both models. It is speculated that this occurs as a result of the variability for the fitted parameters, S (in-plane shear) and to a lesser degree, F 1 2 (interaction parameter). These parameters were established using the [± 15], laminate in compression (for reasons noted in section 5.4.1), which has a relatively low in-plane shear stress. Considering the stress state at failure given above for the [ + 30]s laminate, the strong presence of shear stress influences the ultimate laminate strength. If the shear variability is estimated as high, the laminate variability will be high. 131 Ultimate Stress (MPa) Figure 6.13 - Cumulative Probability Distribution of [± 15]s Laminate in Tension (2 - Dimensional Model) Ultimate Stress (MPa) Figure 6.14 - Cumulative Probability Distribution of [± 15]s Laminate in Tension (3 - Dimensional Model) 132 Ultimate Stress (MPa) Figure 6.15 - Cumulative Probability Distribution of [±30] s Laminate in Tension (2 - Dimensional Model) 133 45 T 0.012 Strain (mm/mm) Figure 6.17 - Stress-Strain Curve Comparison of [± 15]s and [±30] s Laminates in Tension 134 6.4 B E N D I N G V E R I F I C A T I O N In this section, we check the intended purpose of the model - to simulate the mechanical response of a wood strand composite in bending. In bending, all quadrants of the stress field are present simultaneously. Initial specimen failure may be a result of compressive stresses (present in the upper fiber of a positive bending specimen), while ultimate failure may result from predominantly tensile stresses (present in the lower fibers of said specimen). Hence, the bending analysis provides a rigorous check for the simulation model. This section reports the numerical and experimental methodology, as well as the comparative results for the response of [ ± 15]s and [ + 30]s angle ply laminates when subjected to 3 point bending. In addition, bending behaviour for laminated veneer specimens (with no grain angle variation - 0 and 90 degrees) is explored and discussed. 6.4.1 Numerical Analyses Figure 6.18 illustrates the geometric layout and corresponding finite element mesh of the 2 and 3 dimensional models for the [ ± 15]s and [ ± 30]s angle ply laminates under 3 point bending. The same mesh was used for the 0° and 90° laminates with the exception that these laminates were 6 plys in through thickness. The mesh properties for both programs are outlined in Table 6.7. As shown, the size of the 2 dimensional program is unaffected by the number of layers in the laminate. However, the 3 dimensional program becomes unmanageably large very quickly due to the number of elements, even with small specimens. 135 11 mm (4 ply); 17 mm (6 ply) a) 11 mm (4 ply); 17 mm (6 ply) b) Figure 6.18 - Finite Element Mesh for Bending Configuration a) 2 Dimensions b) 3 Dimensions Table 6.7 - Finite Element Mesh Properties of Angle Ply Laminate Bending Analyses Finite Element Mesh Properties 2 Dimensional Model 3 Dimensional Model 0 [±151 [±301 90 0 [±151 [±30] s 90 Number of nodal points 69 525 345 345 525 Number of degrees of freedom 138 1575 1035 1035 1575 Number of elements 44 288 176 176 288 The grading of the mesh, in this case, was partially dictated by computer capacity. Simply, the 3 dimensional model was not able to process the problem using the same element size as for the uniaxial tests (4.5 x 4.5 mm2 in the x - y coordinate system). Therefore, to weigh the effects of modeling error, a parametric study was carried out. Parametric study: In order to validate the element size chosen, two different element sizes were investigated using the 136 3 dimensional program for the above defined bending analysis. The results of an isotropic analysis (one element through the thickness to accommodate processor limitations) using the x - y coordinate mesh outlined in Figure 6.18 was compared to that of a mesh doubly refined in length and height. The input parameters of the analysis were the mean values in Table 6.5 (grain angle of 0°). The load - mid-span displacement curves for both results is shown in Figure 6.19. 8 0 0 . 0 T Displacement (mm) Figure 6.19 - Comparison of Coarse and Fine Mesh for Bending Analysis Although both curves are clearly distinct, it is important to note that both initial stiffness and ultimate load are approximately equal. However, the fine mesh produces a more ductile curve, due to the higher number of Gaussian points present. For this one layer analysis, this is a consequence of the stresses being more easily redistributed with a finer mesh. It remains to be seen if this feature is helpful in accurately predicting the experimental findings. For the time being, the coarse mesh 137 is accepted and its robustness will be borne out through comparison with experimental data. The input properties for the academic laminates differed from that outlined in Table 6.5 only in the values used for tensile capacities as a result of size effect. The tension perpendicular-to-grain values, varying with each Gaussian point, were adjusted using Equation 5.5 as described in Section 6.3.1.1. In so doing, it was assumed that at each Gaussian point the stresses were, on an infinitesimal scale, uniformly distributed about the point. The tension parallel-to-grain strength values were, on the other hand, adjusted for each strand. To do this, a load configuration effect was implemented. Based on similar principles as for size effect, load configuration effect is the result that brittle strength is dependant on the proportion of material that is highly stressed. For example, a pure tensile member, which has its entire volume highly stressed, has a higher probability of a critical flaw and hence a lower strength than does a bending member of the same size which has less of its volume highly stressed. The relationship between tensile and bending failure stress is derived as follows (Clouston, 1995): We equate the probability of failure in bending to the probability of failure in pure tension. Given Equation 5.3 Where we let: z1 = o x b = longitudinal tension stress in bending T 2 = oM = l o n g i t u d i n a l t e n s i o n s t r e s s i n p u r e t e n s i o n 138 Considering the 3 point bending specimen first: If omiX^ = maximum longitudinal tension stress in elastic bending M = bending moment along longitudinal axis M m a x = maximum bending moment then, _ M m a x _ F L max(b) F 6 ^ VbdV (6.1) S 4 Rearranging and simplifying, the expression for maximum elastic load, F is F = — a ,M (M max(b) v ' The longitudinal tension stress in bending at any point in the x - y plane is known from linear elastic beam theory as M(x)-y f F-xY 12 I V 2 Abd: o%b(x,y) = ^ ^ - = V1 flT -7 (M Substituting F from Equation 6.2 and simplifying, we have . 4 ° " max(b) ' X ' y / , 4 \ CTxb(x'y) = — i ~ d — The probability of failure, P^, in bending is now written as xyl dV P - 1 _ ~ v ° ^ ' - a m J L fb ~ ~ 1 e which resolves as follows 1 Cf 4 a p . -T\\T^-X1) D V (6.5) PA, = - e * 139 Finally, V o J,V Ldm xy dV = l - e b/2 0 L/2 m 4a Ldm xy dx-dydz -b/2-d/2 0 2 f 4c m „ V . V Ldm 1-e 0 L/2 xy) dx-dy 2bf40„„V(xy V V L d m ; (p+i)2 = 1 - e v ' v U d m J (p+l)2 = 1 - e v 7 2 b L d f a m „ V o 4 I m V f CT„ 2 V „ = l - e Considering the pure tension case, -V , a P f t = l -e V . I m Equating the two probabilities - V , V „ \ 0 max(t) V m J v u m a x ( b ) I m , 2-(P + l)2 whereupon we obtain the relationship max(t) fmax(b) , 2 V t (p + l ) 2 , (6.8) The tension parallel-to-grain strength for the 0° laminate (for example) is then adjusted as follows: Given: P = 4.2 ^ = om3x{t) = 68.77 MPa V t = 18.8 x 3.2 x 50.3 = 3026.0 mm3 V b (0°laminate) = 17.4/6 x 190 x 19 = 10,469 mm3 then X t 2 - X t 1 68.8 • 132.7 MPa 10469 ,2V t(p + l ) 2 , ^2-3026(4.2-l-l)2. The coefficient of variation is assumed constant so that the standard deviation becomes 132.7 MPa x 0.2667 = 35.4 MPa It is noted here that Equation 6.8 is valid for the current laminate under consideration (ie. with a depth equal to that of one strand). If laminates with several strands through the depth were considered, the formulation would need to be reassessed. For a laminate with 4 strands in depth, for example, the integration in Equation 6.5 would be carried out over the y coordinates of the strand in question. 141 This computation was carried out for all academic laminates and the final values used for input are tabulated in Table 6.8. The tension perpendicular-to-grain strengths were adjusted in the usual manner using Equation 5.5. Table 6.8 - Tensile Strength Properties used for 3 Point Bending Analyses Strength 2 Dimensional 3 Dimensional Mean St. Dev. Mean St. Dev. Parallel-to-grain [0°] (MPa) 132.7 35.4 132.7 35.4 Tension (XJ [±151 (MPa) 134.4 35.8 134.4 35.8 [±301 (MPa) 136.2 36.3 136.2 36.3 [90°] (MPa) 132.8 35.4 132.8 35.4 Perpendicular-to-grain [0°] (MPa) 5.18 0.95 5.75 1.05 Tension (Yj [±151 (MPa) 5.23 0.96 5.80 1.06 [±30] s (MPa) 5.27 0.96 5.85 1.07 [90°] (MPa) 5.19 0.95 5.75 1.05 The stresses are calculated and monitored at the Gaussian points for yielding or failure. However, failure naturally ensues at the extreme fiber in bending, where, in actuality, the longitudinal stresses are highest. To reflect this, the longitudinal stresses are extrapolated linearly from the Gaussian points to the outer edge of the beam prior to being checked in the failure / yield criterion. Referencing Table 6.5, the initial stiffness property of a strand is different depending on tensile or compressive stresses as well as material direction (parallel or perpendicular-to-grain). As previously noted, the bending analysis produces multiaxial states of stress with all combinations of tension, compression, parallel-to-grain, and perpendicular-to-grain stress states. The stress state of each point 142 is not known, however, until the first load (or displacement) is applied. Thus, it was necessary to devise a method to designate the appropriate initial stiffness corresponding to the unknown stress state of a Gaussian point. This was managed by first establishing the stress state for all Gaussian points (in the first iteration of the first increment) using default stiffnesses (tension assumed) with a small displacement increment. Figure 6.20 shows an extreme exploded view of the origin of the load displacement curves for [ ± 15]s angle ply laminates in bending. Being linear elastic in this small range, the problem converges on the first iteration. Thus, in the first iteration of the second increment, the stiffnesses are updated to reflect the stress sign (i.e. if negative, then compressive; if positive, then tensile). The program then proceeds with larger displacement increments based on those stiffness values. ~ D — Simulation Average ._+_ simulated 0 0.002 0.004 0.006 0.008 0.01 Displacement (mm) Figure 6.20 - Exploded View of Origin of Simulated Load Displacement Curves for [± 15]s Angle Ply Laminates 143 6.4.2 Experimental Tests Laboratory tests were conducted on the 4 categories of academic laminates: [0°], [± 15]s, [ + 30]5,and [90°]. Al l tests were performed using a 133 k N capacity Sintech model 30/D universal testing machine. The 3 point bending setup was arranged as illustrated in Figure 6.21. In all cases the clear span was 190 mm. Load was applied under displacement control mode at a constant rate of 1.02 mm/min. to achieve failure in, on average, 5 minutes. Specimen cross sectional dimensions were nominally 19 in depth and 17 mm in width for the [0°] and [90°]laminates -11 mm in width for the [±151 a n d [ ± 30]s angle-ply laminates. Dimensions were measured at ends and mid-span using digital calipers and were found to be reasonably consistent with a maximum coefficient of variation of 1.8 percent for thickness and 0.4 percent for depth. Figure 6.21 - 3 Point Bending Test Setup 6.4.3 Results The simulated and experimental ultimate loads for the 3 point bending analyses are plotted in cumulative probability form in Figures 6.22 and 6.23 for the 2 dimensional and 3 dimensional setups, respectively. The descriptive statistics are outlined in Table 6.9. 144 Table 6.9 - Experimental and Simulated Data for [0°], [± 15]s, [±30] s and [90°] Laminates in 3 Point Bending Config. Statistic Experiment Simulation (Count = 500) {Count} 2d (% error) 3d (% error) [0°] Initial Stiffness Mean (MPa) 12735.7 {20} 8456.4 33.6 9999.4 21.5 Initial Stiffness COV (%) 6.9 7.5 8.7 8.2 18.8 Ultimate Load Mean (N) 2967.5 2204.6 25.7 2622.0 11.6 Ultimate Load COV (%) 8.6 5.4 37.2 6.7 22.1 [±151 Initial Stiffness Mean (MPa) 9023.3 {20} 9172.5 1.6 8473.0 6.1 Initial Stiffness COV (%) 17.2 12.6 26.7 12.0 30.2 Ultimate Load Mean (N) 825.6 757.6 8.2 849.8 2.9 Ultimate Load COV (%) 16.0 9.0 43.8 11.5 28.1 [±30] s Initial Stiffness Mean (MPa) 4292.2 {23} 4870.8 13.5 4925.9 14.8 Initial Stiffness COV (%) 5.7 11.1 94.7 11.6 103.5 Ultimate Load Mean (N) 444.8 301.8 32.1 401.4 9.8 Ultimate Load COV (%) 8.9 7.3 18.0 8.7 2.2 [90°] Initial Stiffness Mean (MPa) 839.5 {19} 143.5 82.9 187.9 89.9 Initial Stiffness COV (%) 15.1 16.1 6.6 6.5 63.8 Ultimate Load Mean (N) 98.8 159.0 60.9 187.6 89.9 Ultimate Load COV (%) 13.0 7.4 43.1 4.7 63.8 The 2 dimensional program is clearly less accurate in predicting the experimental findings than its 3 dimensional counterpart. With regards to ultimate load prediction, in each case, with the exception of the 90° laminate, the 3 dimensional program is more accurate in describing the strength (i.e. both mean and variability). It is speculated that the ability of the 3 dimensional model to account for interlaminar shear stresses leads to a more realistic and therefore more accurate analysis. Although the through thickness stresses are not considered directly in the yield / failure criterion, 145 the mere presence of these stresses in the calculations provide a more precise computation of the in-plane stresses. It is fair to say that the present stochastic bending analysis, with different stiffnesses for each ply, requires a 3 dimensional analysis to account for the complexities. This said, the 3 dimensional program requires approximately 8 times more computer space and time (reference Table 6.7). The 2 dimensional analysis, as an alternative, provides a rough approximation for the average ultimate load (maximum error of 32.1 percent) for the 0°, [±15] s and [±30] s laminates. With regards to predicting initial stiffness of the curves, both programs produce reasonable results. Again, with the exception of the 90° laminate, the worst case result is that of the 0° laminate with an error of the mean of 33.6 percent. It is speculated that this relatively high error may be a result of local crushing at the node points which was prevalent for the 0° laminates in the experiments. The dominant compressive stresses may have contributed to more Gaussian points with compressive elasticity resulting in an overall more compliant laminate stiffness. Ultimate load for the 0° laminate is under-predicted for both the 2 and 3 dimensional model (25.7 % error for 2 D. and 11.6 % error for 3 D). Due to factors of scale, this fact is emphasized in Figures 6.22 and 6.23. It is speculated that this inaccuracy is a result of the programs' inability to simulate local crushing at the load head to the extent seen in the experimental specimens. Referencing Figure 2.2, the interaction parameter, F 1 2 , was determined in the second quadrant of the strength envelope, where the stress state is compression parallel-to-grain and tension perpendicular-to-grain. For the 0° laminate bending case, the third quadrant (compression parallel-to-grain and compression 146 perpendicular-to-grain) governs initially. It may be that the Tsai Wu strength envelope is overly restrictive in this quadrant. A full investigation could be carried out in a future study to assess the accuracy of this theory in all quadrants. In the same respect, alternative strength criterions for wood strand composites could also be surveyed. The 90° laminates are poorly estimated on all accounts. The cause of this may be experimental in nature. The 90° laminates were very fragile and susceptible to cupping. It is possible that the experimental specimens (for tensile tests as well as bending) contained micro-cracks or dimensional irregularities which negatively impacted the results. The load displacement curves for the entire range of experimental and simulated results are given in Figures 6.24 - 6.31. A total of 500 replications were run for each configuration; however, as before, only the stiffest, least stiff, strongest and weakest as well as curve averages for both sets are shown for clarity. A l l simulations, save the 90° laminates, seem to replicate the experimental data well. The ductile behaviour for the 0° and [± 15]s laminates is obvious, as is the onset of ultimate failure by abrupt brittle behaviour. The program, in its present state, however, is unable to capture the experimental feature of recovery where, upon abrupt local failure, the stresses redistribute to surpass the previous peak. This is likely a consequence of the slow release of stresses upon brittle failure (per Table 4.2) to assure convergence of the system. The stresses do indeed redistribute numerically, but more gradually. 147 Ultimate Load (N) Figure 6.22 - Cumulative Probability Distribution of Academic Laminates in 3 Point Bending (2 - Dimensional Model) 0 500 1000 1500 2000 2500 3000 3500 Ultimate Load (N) Figure 6.23 - Cumulative Probability Distribution of Academic Laminates in 3 Point Bending (3 - Dimensional Model) CP Laminate Results: 3500 T D is placement (mm) Figure 6.24 - Load-Displacement Curves of [0°] Laminate in Bending (using 2 - Dimensional Model) 3500 T Displacment (mm) Figure 6.25 - Load-Displacement Curves of [0°] Laminate in Bending (using 3 - Dimensional Model) 149 [±15]s Angle Ply Laminate 1200 T Displacement (mm) Figure 6.26 - Load-Displacement Curves of [+ 15]s Angle Ply Laminate in Bending (using 2 - Dimensional Model) Displacement (mm) Figure 6.27 - Load-Displacement Curves of [± 15]s Angle Ply Laminate in Bending (using 3 - Dimensional Model) 150 f±30]s Angle Ply Laminate 500 T Displacement (mm) Figure 6.28 - Load-Displacement Curves of [±30] s Angle Ply Laminate in Bending (using 2 - Dimensional Model) 600 T Displacement (mm) Figure 6.29 - Load-Displacement Curves of [±30] s Angle Ply Laminate in Bending (using 3 - Dimensional Model) 151 [±90]s Angle Ply Laminate 250 T Displacement (mm) Figure 6.30 - Load-Displacement Curves of [90°] Laminate in Bending (using 2 - Dimensional Model) 250 T tu o I 1 1 1 1 1 1 1 -2 0 2 4 6 8 10 12 14 Displacement (mm) Figure 6.31 - Load-Displacement Curves of [90°] Laminate in Bending (using 3 - Dimensional Model) 6.5 Parallam* Investigation As a preliminary assessment of the model, a comparison of simulated and experimental data for Parallam* was conducted. It should be understood that at this stage in the model development, complete accuracy is not expected. The model is not yet formulated to address all the complexities of the spatial relationship between individual wood elements in wood composites. It was felt that a general comparison at this stage would simply be useful in providing insight into ultimately simulating commercial products. This section addresses the current method involved in modeling grain angle deviations in Parallam", as well as a comparison between simulated and experimental findings for tensile and bending configurations. The 2 dimensional program, as a consequence of the assumptions made for the classical lamination theory, is restricted to symmetric laminate analysis. The 3 dimensional program, on the other hand, is configured to analyze laminates with any stacking sequence. Because the lay-up of strands - and therefore grain angle - is known to be random for Parallam*, the following investigation focuses solely on the 3 dimensional model. 6.5.1 Simulation of Parallam* Strand Lay-up Special preparations were required for the 3 dimensional model to accommodate the complex strand lay up of Parallam*. Considering the strands are randomly aligned with the longitudinal axis, it was decided that the grain angle could be randomly generated based on a uniform distribution of approximate observed measurements. These grain angle characteristics were physically measured as outlined below: 153 Ten individual boards of Coastal Douglas-fir 2.0E (ie. 2.0xl06 psi or 13.8 GPa modulus of elasticity) Parallam® were used for both tensile specimens as well as grain angle distribution measurements. The 3 dimensional mapping region lay directly adjacent to where the specimens were cut from the board, as shown in Figure 6.32. Visual measurements of maximum grain angle were estimated over 25 x 16 mm2 elements along the length of the board. Maximum grain angle was used because it was assumed that for the majority of cases, the maximum deviation of grain would govern strength (for example, a knot within an element would be the weakest point and would constitute a 90° grain angle). A layer of 3mm was then planed from the board and measurements were taken from the new surface. This was done for 4 successive layers to produce a 12 mm through-thickness profile. Figure 6.32 - Cutting and Grain Angle Mapping Template for Parallam® The data was then interpreted into a frequency distribution as seen graphically in Figure 6.33. From this data, a uniform distribution was created with equal probability of having the same angle as positive or negative. This uniform distribution formed the basis of random angle generation within the program. 154 1400 1200 1000 - H c 800 f 600 4H 400 200 0 n •+- + n 10 20 30 40 50 60 70 80 90 Angle to Grain Figure 6.33 - Frequency Distribution of Measured Parallam® Grain Angle 6.5.2 Parallam* Tension Comparison 6.5.2.1 Numerical Analysis Having completed this addition to the program, the input file was then formulated for the tensile specimen. The finite element mesh assumes a uniform stress field between the grips and is illustrated in Figure 6.34. This mesh considers the strand lay up to be substantially simplified in comparison to that of Parallam*. The model assumes a total of 4 strands through the thickness per specimen. This is not wholly realistic as the strands in Parallam* are not planar or continuous and tend to meander in three directions. Furthermore, voids are prevalent throughout the material. It is possible that future versions of the program could incorporate this geometric variability. 155 Figure 6.34 - Finite Element Mesh for Parallam® in Tension To account to some degree for the geometric variability, the tension parallel and perpendicular-to-grain strengths were assumed to vary from integration point to integration point. However, because the dimension in the through-thickness (z) direction was very narrow (3.1 mm ), the adjacent points in the z direction were given the same tension parallel-to-grain properties. These values were adjusted in the usual manner (described in Section 6.3.1.1) for size effect. The strand properties are that of the heartwood Douglas-Fir database established in Chapter 5. It is noted, however, that Parallam® may contain both of heartwood and sapwood. 6.5.2.2 Experimental Tests The rectangular specimens were shaped into a 'dogbone' profile to initiate failure within the gauge length. Without this special profile, failure would consistently occur in or near the grips due to stress concentrations from grip crushing. The dimensions of the specimens are given in Figure 6.35. 156 30 mrxJ Radius — 47.5 mm 4 16 mm ^ Jt -+ 407 mm +-590 mm -sir Figure 6.35 - Parallam® Tensile Specimen Dimensions A total of 56 specimens were tested using a 222 k N MTS machine with a mechanical wedge-type grip system. The loading rate was 1.27mm/min. An extensometer with a 25.4 mm gauge length was employed to determine tensile MOE at the mid-span of the specimen. Typical failure occurred near mid-span in a splintered manner as depicted in Figure 6.36. 6.5.2.3 Results The simulated data is compared with the experimental results in Table 6.10 and Figure 6.37. In this preliminary investigation, the elastic modulus was under-predicted which may be explained by the fact that each strand was allocated an angle to grain without variation. In actuality, the grain angle changes more frequently throughout the specimen. A more accurate depiction could be to vary the grain angle with each element. 157 The mean strength is substantially over-predicted. This is explained by the fact that the model, in its current state, is not formulated to manage all of the complexities of Parallam®. As previously noted, voids were not accounted for, which could lead to a significant over-estimation of strength. Moreover, the input properties were solely for heartwood Douglas-fir, although Parallam® is comprised of both heartwood and sapwood (the latter appearing slightly weaker in experiment - see Table 5.1) . Table 6.10 - Experimental and Simulated Data for Parallam® in Tension Statistic Experiment Simulation Percent Error (count = 56) (count = 500) (%) Elastic Modulus Mean (MPa) 14,222.9 11,945.1 16.0 Elastic Modulus C O V (%) 26.8 25.2 6.0 Strength Mean (MPa) 33.6 62.6 86.3 Strength C O V (%) 25.2 20.4 19.1 158 Ultimate Stress (MPa) Figure 6.37 - Cumulative Probability Distribution of Parallam® in Tension (3 - Dimensional Model) 6.5.3 Parallam® Bending Comparison A bending comparison was performed between simulated and experimental data for Parallam®, as was done for tension. It is again noted that this investigation is simply intended as a primary check. Modifications to the model to consider void volume, for example, should be further investigated. 6.5.3.1 Numerical Analyses Figure 6.38 depicts the specimen geometry and finite element mesh for the bending analysis. The material is assumed to be comprised of 4 separate strands through the thickness each with different stiffness and strength properties. As before, for Parallam® subjected to tensile loading, the tensile strengths (both parallel and perpendicular-to-grain) were assumed to vary between Gaussian points to better approximate the random nature of the strands. Treating the stress distribution about each 159 Gaussian point as uniformly distributed, the tensile strengths are adjusted for size effect using Equation 5.5. A summary of adjusted tensile strengths is provided in Table 6.11. Al l other input properties are taken from the acquired database for the Douglas-fir heartwood strands and are outlined in Table 6.5. 20 @ 9.5 = 190 mm Figure 6.38 - Finite Element Mesh for Parallam* in 3 Point Bending Table 6.11 - Tensile Input Properties for Paral am* Bending Analysis Tensile Property Mean Standard Deviation Parallel-to-grain Strength (Xj (MPa) 171.27 45.7 Perpendicular-to-grain Strength (YJ (MPa) 4.62 0.85 The grain angle for each strand was randomly generated in accordance with the uniform distribution shown in Figure 6.33. 160 6.5.3.2 Experimental Tests Bending tests were conducted using a Sintech model 30/D universal testing machine with a capacity of 133 k N - the same as for the academic material in section 6.4.2. The setup is illustrated in Figure 6.39 for a typical Parallam® specimen. Load was applied with a uniform loading rate of 0.9mm/minute under displacement control mode. Specimen cross sectional dimensions were measured at ends and mid-span and were found to have high consistency with a maximum coefficient of variation of 1.7 percent for thickness and 0.7 percent for depth. Figure 6.39 - 3 Point Bending Test Setup for Parallam® 6.5.3.3 Results A statistical summary of experimental and simulated data is outlined in Table 6.12 and visually represented in Figure 6.40. The entire range of load-displacement curves are shown in Figure 6.41. Considering the model is not regarded as complete at this stage, the results are quite reasonable. Ultimate load average is slightly over-predicted which is understandable given that voids are not considered and that the database comprised only of heartwood without the lower grade sapwood (which is present in equal proportions in Parallam*). The initial stiffness average is slightly under 161 predicted, although referencing Figure 6.41, the range of results lay close to the experimental bounds. The mode of failure is captured particularly well for Parallam*. In experiment it was observed that in some cases, the specimens tended to form compression wrinkles near the location of applied load prior to a splitting tensile failure at the bottom edge. This sequence of events was readily seen in both the experimental and simulated load - displacement curves. The curves follow a nonlinear path during predominant compressive failure and ultimately fail in a brittle manner due to predominant tensile stress failure. The assumptions made in this analysis greatly simplify the geometric properties of Parallam*. A more rigorous study to ascertain and subsequently model the Parallam* strand layup should be the next step. From this preliminary analysis, however, it appears that the theory presented herein is appropriate and capable of ultimately predicting the mechanical behaviour of wood strand composites in general. Table 6.12 - Experimental and Simulated Data for Parallam* in 3 Point Bending Statistic Experiment Simulation (Count = 500) {Count} 3d (% error) Initial Stiffness Mean (MPa) 9239.4{41} 8050.8 12.9 Initial Stiffness C O V (%) 12.3 12.7 3.3 Ultimate Load Mean (N) 1232.8 1317.9 6.9 Ultimate Load C O V (%) 16.4 13.5 17.7 162 400 600 800 1000 1200 1400 1600 1800 2000 Ultimate Load (N) Figure 6.40 - Cumulative Probability Distribution of Parallam* in 3 Point Bending (3 - Dimensional Model) 2000 T Displacement (mm) Figure 6.41 - Load-Displacement Curves of Parallam* in 3 Point Bending (using 3 - Dimensional Model) 163 7 - Conclusions 7.1 S U M M A R Y Despite the advantages of having a computational strength prediction model for wood composite research and design, very few studies have been devoted to this subject in the literature. Work that has been carried out has concentrated to a large degree on linear elastic constitutive theory. The physical nature, however, of the wood strand composites utilized in today's lumber products industry necessitate a more rigorous approach. In the present study, a continuum mechanics approach, considering elastic and inelastic behaviour up to and beyond failure, has been investigated. The primary objective of the present study was to develop a computational model to simulate the mechanical response of wood strand composites subjected to bending. A finite element based program, COMAP (COmposite Analysis Program), was developed. Two versions were created. A 2 dimensional version was developed which utilizes classical lamination theory (CLT) to evaluate resultant forces for symmetric multilayered laminates. This version, although efficient in computational time and general usage, is restricted by the assumption of symmetry and the assumptions of the CLT. Notably, plane stress is assumed within each layer and thus stresses in the through-thickness direction are neglected. A 3 dimensional version was created to provide a means of assessment for nonsymmetrical laminates and also to gauge the influence of through - thickness stresses. An orthotropic elastic - plastic - failure constitutive model was implemented into both versions of the program. The nonlinear formulation (beyond the elastic regime) was managed using classic orthotropic incremental plasticity theory. The constitutive model utilized the Tsai-Wu yield / 164 failure criterion to describe the initial and subsequent yield surfaces assuming an associated flow rule. A l l stresses were monitored for each iteration of each increment of loading (or displacement) and were initially transformed into the principal material directions prior to use in the yield criterion. Ultimate failure was defined using ultimate principal strengths (experimentally determined) to establish an upper bound yield surface. Upon reaching the upper bound surface, an integration point either 1) failed in a brittle manner or 2) first strain - hardened and then eventually failed in either a ductile or brittle mode. The latter choice was decided based on whether the combination of stresses at the point at the time of failure were predominantly tensile or predominantly compressive. Stiffness and strength properties of the material were input as random variables for both programs. This allowed for many replications of one configuration through Monte Carlo simulations to generate entire sample output. Moreover, the programs were formulated for stochastic analyses with varying strength and stiffness values within each member. A strand property database of strength and stiffness for Douglas-fir heartwood was acquired which was necessary as input into the programs. Laminated boards were first fabricated from which specimens were cut and then prepared for experimental tests. For shear strength, modulus of ridgidity and the interaction parameter of the Tsai-Wu theory, a least square minimization of error between simulated and experimental compression strength of [±15] s angle ply laminates was conducted. This calibration procedure provided the mean and standard deviation of the 3 variables. 165 To verify the performance of the programs, the solutions to various problems were compared to known analytical solutions and that of other software. No one other software package has all of the capabilities of COMAP and consequently model verification was handled systematically to first address plastic deformation accuracy, then to check the linear elastic orthotropic finite element formulation, and finally to confirm the stress and strain transformation process required for multilayered laminates. The effectiveness of the programs to replicate experimental findings was demonstrated in the final chapter. The results of numerical simulations for angle ply laminates in tension, compression and bending compared favorably with experimental data. Furthermore, preliminary comparisons of tensile and bending results for Parallam" showed the proposed technique to have great potential to ultimately model commercial wood strand composites. 7.2 CONCLUDING REMARKS The favorable agreement between the numerical simulations and experimental data demonstrate the accuracy of the present technique. The analyses conducted herein point to the fact that orthotropic plasticity theory in conjunction with a stochastic finite element approach appear to offer an elegant solution to predicting wood composite behaviour. Indeed, this technique is general enough to be used for a variety of laminated composites provided the experimental database for the layer component of the composite is available. A full account of required data input for the present programs is provided in Appendix C along with sample outputs. With respect to this, it is cautioned that the present academic form of the programs requires further detailed investigation prior to use - particularly with commercial products. 166 7.3 FUTURE RELATED RESEARCH The subject offers many interesting areas for further research. The following discussion provides a few ideas for extension and improvement of the current model. An obvious improvement at this stage would be to upgrade the program's ability to describe the physical structure of commercial wood composites. One basic feature lacking from the present model is the ability to capture the complex 3 dimensional staggered stacking sequence of the strands which is the nature of wood composites although not common with advanced aero-spatial type composites. Several studies have been devoted to understanding the complex spatial relationship between individual wood elements in wood composites (Dai and Steiner, 1993; Ellis et al., 1994; Oudjehane et al., 1998; Lu, C. 1999). To incorporate these techniques into the present formulation may prove to be a difficult, although rewarding, task. It must entail consideration of density variations, potential null volumes or voids as well as varying 3 dimensional geometry of elements. Many smaller alterations to the program could be investigated. Given the stochastic nature of the model, one could make any presently deterministic parameter stochastic. For example, the grain angle, in reality, likely varies within, as well as between, strands. Given an accurate mathematical depiction of this variation, properly implemented into the program, the model's present accuracy to simulate the experimental data may be enhanced. Presently, the compression properties are correlated in generating the compressive stress - strain curves in the principal material directions. Correlation between all variables could be considered, for example, within a strand. Also, the technique is not limited to the Tsai-Wu criterion nor the assumptions made for the description of plastic flow. Other yield surfaces (possibly considering non-associated flow) may be more 167 appropriate for other materials. Further areas of interest would include extension of the model to study wood based composite panels or investigation of nonlinear viscoelastic behaviour of wood composites. There is much to be learned of creep and creep rupture performance of wood composites under long term loading and varying temperature and humidity conditions. Also of strong interest would be adaptation of the model to investigate the interaction between mechanical connections and wood composites. The technique could potentially be extended to optimize strand layup for prescribed connections. 168 R E F E R E N C E S ASTM Standard D143 1994. Standard Methods of Testing Small Clear Specimens of Timber. American Society for Testing and Materials, Philadelphia, Pa. ASTM Standard D198 1994. Standard Methods of Static Tests of Timbers in Structural Sizes. American Society for Testing and Materials, Philadelphia, Pa. ASTM D2016 1991. Standard Test Methods for Moisture Content of Wood. American Society for Testing and Materials, Philadelphia, Pa. ASTM Standard D3518 / D3518M 1991. Standard Practice for In-Plane Shear Stress-Strain Response of Unidirectional Polymer Matrix Composites. American Society for Testing and Materials, Philadelphia, Pa. Barrett, J.D.; Foschi, R.O.; Fox, S. P. 1975. Perpendicular-to-grain strength of Douglas-fir. Can. J. Civ. Eng., 2(1), pp. 50-57 Barrett, J.D.; Lam, F.; Lau, W. 1995. Size effects in visually graded softwood structural lumber. J. of Mat. in Civ. Eng., 7(1), pp. 19-30 Bodig, J.; Jayne, B A . 1993. Mechanics of Wood and Wood Composites. Krieger Publishing Co. Malabar, Florida Cha, J. K.; Pearson, R .G. 1994. Stress Analysis and Prediction in 3-Layer Laminated Veneer Lumber: Response to Crack and Grain Angle. Wood and Fiber Sc., 26(1), pp. 97-106 Chen, W.F.; Han, D.J. 1988. Plasticity for Structural Engineers. Springer-Verlag, New York Inc. New York 169 Clouston, P. 1995. The Tsai-Wu Strength Theory for Douglas-fir Laminated Veneer. Master of Applied Science Thesis, University of British Columbia, Vancouver, B.C. Canada Clouston, P.; Lam, F.; Barrett, J.D. 1998. The Interaction Term of the Tsai-Wu Theory for Laminated Veneer, Journal of Materials in Civil Engineering, 10(2), pp. 112-116 Clouston, P.; Lam, F; Barrett, J.D. 1998. Incorporating Size Effects in the Tsai-Wu Strength Theory for Douglas-fir Laminated Veneer. Wood Science and Technology, 32 (1), pp.215-226 Conners, T.E. 1989. Segmented Models for Stress-Strain Diagrams. Wood Sc. and Tech., 23(l),pp. 65-73 Cook, R.D.; Malkus, D.S.; Plesha, M.E. 1989. Concepts and Applications of Finite Element Analysis. John Wiley and Sons. New York Cowin, S.C. 1979. On the Strength Anisotropy of Bone and Wood. ASME Trans., Journal of Appl. Mech., 46(4), pp. 832-837. Dai, C ; Steiner, P. R. 1993. Spatial structure of wood composites in relation to processing and performance characteristics. Part 1. Rationale for model development. Wood Science & Technology, 28(1), pp. 45-51. Ellis, S.; Dubois, J.; Avramidis, S. 1994. Determination of Parallam Macroporosity by Two Optical Techniques. Wood and Fiber Science, 26(1), pp.0-77 Foschi, R.O.and Yao, Z.C. 1986b. Another Look at Three Duration of Load Models. Procedings, IUFRO Wood Engineering Group Meeting, September; Florence, Italy: Volume 2, 19-1-1 Gibson, R. F. 1994. Principles of Composite Material Mechanics. McGraw-Hill Companies. New York 170 Goodman, J. R.; Bodig, J. 1971. Orthotropic Strength of Wood in Compression. Wood Science 4(2), pp. 83-94. Hill, R. 1948. A Theory of the Yielding and Plastic Flow of Anisotropic Metals. Proceedings of the Royal Society, Series A, Vol. 193, pp. 281-297 Hull, D. 1981. An Introduction to Composite Materials. Cambridge University Press. Cambridge Hunt, M. O.; Suddarth, S.K. 1974. Prediction of Elastic Constants of Particleboard. Forest Prod. J., 24(5), pp. 52-57 Jones, R ,M. 1975. Mechanics of Composite Materials. McGraw-Hill, New York Kobetz, R.W.; Krueger, G. P. 1976 Ultimate Strength Design of Reinforced Timber-Biaxial Stress Failure Criteria. Wood Science, 8(4), pp. 252-261 Lam, F., and Varoglu, E. 1990. Effect of Length on the Tensile Strength of Lumber. Forest Products J., 40(5), pp. 37-42 Liu,J.Y. 1984. Evaluation of the Tensor Polynomial Strength Theory for Wood. J. Comp. Mat., 18, pp. 216-225 Lu, C. 1999. Organization of Wood Elements in Partially Oriented Flakeboard Mats. Ph.D. Thesis, University of British Columbia, Vancouver, B.C. Canada Madsen, B.; Buchanan, A . H . 1986. Size Effects in Timber Explained by a Modified Weakest-Link Theory. Can. J. Civ. Eng., 13(2), pp. 218-232 Madsen ,B. 1990. Size Effects in Defect Free Douglas-fir. Can. J. Civ. Eng., 17, pp. 238-242 171 Maghsood, J.; Scott, N.R.; Furry, R.B.; Sexsmith, R. 1973. Linear and Nonlinear Finite Element Analysis of Wood Structural Members. ASAE Trans. Paper No.70-922, pp. 490-496 Mase, G.E. 1970. Theory and Problems of Continuum Mechanics. Schaum's Outline Series, McGraw-Hill, Inc. Nahas, M . N . 1986. Survey of Failure and Post-Failure Theories of Laminated Fiber-Reinforced Composites. J. Comp. Tech. and Res., 8(4), pp. 138-153 Narayanaswami, R.; Adelman, H . M . 1977. Evaluation of the Tensor Polynomial and Hoffman Strength Theories for Composite Materials. J. Comp. Mat., 11, pp. 366-377 Norris, C.B. 1962. Strength of Orthotropic Materials Subjected to Combined Stress. U.S. Forest Products Lab. Rep. 1816, FPL, Madison, Wisconsin. Ochoa, O.O.; Reddy, J.N. 1992. Finite Element Analysis of Composite Laminates. Kluwer Academic Publishers, Dordrecht; Boston Oudjehane, A.; Lam, F.; Avramidis, S. 1998. Forming and Pressing Processes of Random and Oriented Wood Composite Mats. Composites Part B: Engrg., 29B, pp. 211-215 Owen, D.R.J.; Hinton, E. 1980. Finite Elements in Plasticity. McGraw-Hill, New York Perkins, R.W. 1967. Fundamental Concepts Concerning the Mechanics of Wood Deformation: Strength and Plastic Behaviour. Forest Prod. J., 17(4), pp. 57-68 Pagano, N.J.; Pipes, R .B. 1971. The influence of Stacking Sequence on Laminate Strength. J. Comp. Mat., 5, pp. 50-57 Pipes, R.B.; Cole, B.W. 1973. On the Off-Axis Strength Test for Anisotropic Materials. J. Comp. Mat., 7, pp. 246-256 172 Rowlands, R.E. 1985. Strength (Failure) Theories and Their Experimental Correlation. Handbook of Composites, Vol.3 - Failure Mechanics of Composites, pp. 71-125 Sandhu, R.S. 1972. A Survey of Failure Theories of Isotropic and Anisotropic Materials. U.S. Air Force Tech. Report No. AFFDL-TR-72-71, Wright Patterson AFB, O H Shaler, S.M.; Blankenhorn P.R., 1990. Composite Model Prediction of Elastic Moduli for Flakeboard. Wood and Fiber Science, 22(3), 1990, pp.246-261 Sharp, D.J.; Suddarth, S.K. 1991. Volumetric effects in structural composite lumber. In: Proceedings, Int. Timber Eng. Conf., London, England. 3, pp. 427-433 Sharp D.J. 1996. Manufactured Structural Composite Lumber. Inter. COST 508 Wood Mech. Conf., Stuttgart Germany, pp. 273-283 Shih, C.F.; Lee, D. 1978. Further Developments in Anisotropic Plasticity. Trans. ASME, J. Eng. Mat. and Tech., 100, pp. 294-302 Siau, J.F. 1984. Transport Processes in Wood. Springer Verlag, Berlin Sokolnikoff, I.S. 1956. Mathematical Theory of Elasticity. McGraw-Hill, New York Suhling, J .C; Rowlands, R.E.; Johnson, M.W.; Gunderson, D.E. 1984. Tensorial Strength Analysis of Paperboard. Experimental Mech., 25(1), pp. 75-84 Timoshenko, 1972. Mechanics of Materials. New York: Van Nostrand Reinhold, pp. 306-309 Triche M . H.; Hunt M.O. 1993. Modeling of Parallel-Aligned Wood Strand Composites. Forest Products J., 43(11/12), pp. 33-44 173 Tsai, S. W.; Wu, E. M . 1971. A General Theory of Strength for Anisotropic Materials. J. of Comp. Mat., 5, pp. 58-80. van der Put, T .A.C.M. 1982. A General Failure Criterion for Wood. IUFRO Timber Engineering Group Meeting, Paper 23 [Sweden], JJJFRO, Vienna Vaziri, R.; Olson, M.D. ; Anderson, D.L. 1991. A Plasticity-Based Constitutive Model for Fibre-Reinforced Composite Laminates. J. of Comp. Mat., 25, pp.512-535 Wang Y. T.; Lam, F. 1998. Computational Modeling of Material Failure for Parallel-Aligned Strand Based Wood Composites. Computational Material Science 11, pp. 157-165 Weibull, W. 1939. A Statistical Theory of the Strength of Materials. Proc, Royal Swedish Inst., No. 151, Stockholm, Sweden Whang, B. 1969 Elasto-Plastic Orthotropic Plates and Shells. Proc. Symposium on Application of Finite Element Method in Civil Engineering, Vanderbilt University Tennessee, pp. 481-515 Whitney, J .M. and Browning, C.E. 1972. Free-Edge Delamination of Tensile Coupons. J. of Comp. Mat., 6, pp.300-303. Wu, E. M . 1974. Phenomenological Anisotropic Failure Criterion Mechanics of Composite Materials (ed. G.P. Sendeckyj) New York: Academic Press Yang, T. Y. 1986. Finite Element Structural Analysis. Prentice-Hall, Inc., Englewood Cliffs, New Jersey Zweben, C. 1994. Size Effect in Composite Materials and Structures: Basic Concepts and Design Considerations. NASA Conference Publication (Vol/Iss:3271) pp. 197-217 Zienkiewicz, O.C. and Taylor, R.L. 1971. The Finite Element Method: Volume 1. Basic Formulation and Linear Problems. McGraw-Hill Book Company, London — Appendix A - Hardening Model "To determine the nature of the subsequent yield surfaces is one of the major problems in the work-hardening theory of plasticity" (Chen and Han, 1988). The manner in which the yield surfaces evolve beyond initial yield is dictated by the hardening rule. Many rules have been proposed, the most widely known being the simple extremes, isotropic and kinematic hardening as described in Chapter 3. In this study, we investigate an anisotropic rule whereby the subsequent yield surfaces satisfy f EE a2(ai,ai,Uii(X))-'k2(X) = 0 (A.l) such that k varies proportionally with effective plastic strain per Equation 3.45 and the nonlinear parameters of M ; j vary non-proportionally described as follows: Only two strength variables evolve plastically in our case, X c and Y c . The other strength variables are assumed to fail in an abrupt brittle manner. Following a procedure initially proposed by Whang (1969), X c and Y c are updated upon detection of yielding and for each subsequent iteration. The updated values are calculated by equating the plastic work done (W15) during plastic deformation in a uniaxial test to that produced by the effective stress and effective plastic strain for a multi-axial state of stress. Referencing Figure A. 1, the plastic work done in a uniaxial test which is approximated by a bilinear stress - plastic strain curve is 175 w p = — ( r , - - r 0 l 2 ) ^ri„ (A.2) where E p is the plastic modulus and Tt and r o i represent the initial and subsequent yield values for either X c or Y c . Figure A . l - Bilinear Stress - Plastic Strain Curve In the same way, the work done by the effective yield stress is 2H 7( k ~K) (A.3) It is noted that the effective stress - effective strain relationship is obtained from either one of the uniaxial stress-strain curves (ie. parallel-to-grain or perpendicular-to-grain). The analysis is independent of choice made. Equating Equation A.2 to A.3, we have iF(k2-k«) = ^i-(r-r«) (A.4) Rearranging, the updated values for X c or Y c is found to be 176 or specifically, r ^ ^ k ' - k ^ + r i (A.5) x L = ^ - ( k ? - k * ) + x * Y - = i r ( k ' - k o ) + Y - 2 (A.6a) (A.6b) 177 Appendix B - Parametric Study A preliminary study was carried out to establish a computer-time efficient yet numerically accurate finite element mesh for the uniaxial analyses. Three mesh sizes using simple rectangular elements were investigated. In one case, a rather coarse mesh was used dividing the x and y dimensions into thirds (3x3 grid). The second was refined to a 4 x 4 grid and the third was refined further to a 5x5 grid. Al l three are illustrated in Figure B . l . Figure B.l - Various Mesh Grades Investigated for Uniaxial Analyses: a) 3 x 3 grid, b) 4 x 4 grid, c) 5 x 5 grid 178 A [ + 30]s, angle-ply laminate subjected to compressive displacement was modeled using the 3 dimensional program. To ensure comparability between meshes, a deterministic analysis was performed. The same strength and stiffness variables, outlined in Table B. 1, were used for each case. Table B . l - Input Properties for Parametric Study P r o p e r t y M e a n St . D e v . Parallel-to-grain Elastic Modu lus (E X l ) (MPa) 15463 0 Tens ion Strength ( X J (MPa) 72.8 0 Perpendicular-to-grain Elastic Modu lus (E Y t ) (MPa) 91.2 0 Tens ion Strength ( Y J (MPa) 6.5 0 Parallel-to-grain Elastic M o d u l u s (E X ( ) (MPa) 10090 0 Compress ion Y i e l d Strength ( X J (MPa) 67.3 0 Tangent M o d u l u s ( E X c ' ) (MPa) 1926 0 Ul t imate Strength ( X c u) (MPa) 76.5 0 Perpendicular-to-grain Elastic M o d u l u s (E Y c ) (MPa) 490 0 Compress ion Y i e l d Strength ( Y J (MPa) 15.4 0 Tangent Modu lus ( E Y c ' ) (MPa) 110 0 Ul t imate Strength ( Y c u) (MPa) 18.2 0 Interaction Parameter F.2 1.067X10"03 0 In plane Shear Elastic M o d u l u s (G) (MPa) 392.2 0 Strength (S) (MPa) 11.4 0 Poisson's Rat io V ,2 0.32 -The results of the three analyses are given graphically in Figure B.2. The ultimate load for each case varies slightly with progressively lower predictions for the finer grids (26.8 MPa, 26.3 MPa and 25.5 MPa for the coarse, average and fine mesh, respectively). This was a result of stress concentrations 7/9 at the boundaries. The smaller elements produced higher stresses for equal displacement and thus failed earlier. Considering the results were relatively consistent, it was decided that the best element for both efficiency and accuracy would be the 4x4 grid. 30 T 0.012 Strain (mm/mm) Figure B.2 - Stress - Strain Diagrams of 3 Investigated Grids 180 Appendix C - Program Layout A stochastic, materially nonlinear finite element based F O R T R A N 77 program has been developed for this study. The program, entitled COMAP (COMposite Analysis Program), has extended capacity to perform Monte Carlo simulations. The program was adapted from a 2-dimensional nonlinear program by Owen and Hinton (1980). Both a 2-dimensional and 3-dimensional version of the program have been developed and investigated. This Appendix describes the most important aspects of the program. A flowchart highlighting the essential steps in the program is provided in Figure B . l . Sample input and output files for the 3 dimensional program have been included and described for completeness. 181 Start a, o o a a u o a CL, o a o Input data from input file 1 for random variable setup (See sample input file 1) I Generate arrays of random numbers for random strengths and stiffnesses I For each replication, enter finite element analysis subroutine (calculate load - displacement curve) No Input data from input file 2 for finite element data (See sample input file 2) Zero arrays for data which accumulates Increment loads (or displacments) according to factors supplied by input file 2 Calculate element stiffnesses for elastic (Equation 4.18 or 4.32) or elastic-plastic behaviour (reference Equation 3.55) using random variables I Solve system of simultaneous equations (Equation 4.37) Generate strength values Calculate equivalent nodal forces (Equation 4.48) Evaluate effective stress (Equation 3.21) Increment stress (Equation 4.40 or 4.45) Check to see if solution has converged (Equation 4.39) Yesl Output result for this increment (See sample output files) 31 Next replication End Figure B . l - Program Flowchart (Adapted from Owen and Hinton, 1980) 182 I N P U T F I L E l ( S A M P L E ) Line # 1.OOE+00 0.OOE+00 0.OOE+OO 0.OOE+OO 4.77E-01 8.79E-01 0 . OOE + OO 0.OOE+OO 2.00E-01 -3.52E-01 9 . 14E-01 0.OOE+OO 5.54E-01 6 . 54E-01 2 . 84E-01 4 . 30E-01 1.OOE+OO 0 . OOE + OO 0. 00E+OO 0 . OOE + OO 1. 33E-01 9.91E-01 0.OOE+OO 0. 00E+OO 1. 05E-01 -4.33E-01 8.95E-01 0.OOE+OO 3.69E-01 7.66E-01 4.04E-01 3.38E-01 Line 1-32: Values of lower triangle correlation matrix for compression properties (Tibles 5.5 and 5.7) 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 3 D _ 3 p n t l 5 . t x t 35 3D_3pnt3 0 . t x t 36 37 1408 50 Line 33: Number of input files for finite element analysis Lines 34-35: Name of input files for finite element analysis Line 36: Number of Gauss points per layer, number of layers Line 37: Number of replications per random seed pair 183 38 10 Line 38: Number of random seeds 39 1 1 1 2 2 2 Lines 39-48: Random seeds 40 2 2 2 3 3 3 41 3 3 3 4 4 4 42 4 4 4 5 5 5 43 5 5 5 6 6 6 44 6 6 6 777 45 111 888 46 888 9 9 9 4 7 1 2 3 4 5 6 48 4 5 6 7 8 9 184 w < < N W i - l VD 0> O CN i H H H CN CN CN CN CO LA VO CO Ifl OD Ol H CO Tjt L/l CO o i n H CN i n r~ CN CN CN H CN CN cn 01 C i n cd •* 2 H vo o\ CN i n M 3 =*fc CO o> O H CN "1 <* VO <*\ i n vo Cv CO 0\ o H CN cn m i n vo N Cv Cv OD 00 00 00 00 00 CN CM CN CN CN CN CN cn "1 cn cn cn i-H i-H t-H H H H rH i n i n i n i n i n i n i n i n i n i n i n i n 185 in H ro ro ro ro ro ro ro ro vo VD vo vo m n in co in o in o a in - H o T3 ' C H tu , CD XI ro CO CD CD L n v o r ^ c o c r i O r s r \ ) r o ^ i n V o N c o c ^ o r H r \ ) r o ^ i n v o t ^ o D a i O r H r s j r o r o r o r o r o r o ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ i n i n i n i n i n i n i o L o L o L o v o v o v o v o i n i n i n i n i n i n i n i n i n i n i o v n t n i n L n i n i o i n i o i n i n i n i n i n u i i n i n u ^ vo t-. in in oo oo 186 e o • »«* +-4 8 O tq >~i xs • to N 8 S t/j h) R tu s <u \ o R =tfc 8 S" "S ^ V J a N P\ s a, a O a ,o ^ ^ 1 « is -Si 5^ =tfc =tfc J S » a o 8 8 3 •« S If J *v r~ *U •>* s • 2 s 4-J s » S 3 ->* 4-i ^ v U". lq tq ^ 5 ~« tq O 8 8, e o •»-» 4-i • o 3 9 9 s R 6 R "v. tu R R rt a rt a , x W 3 a , 8 t- =tfc 0-3 tq < R R 3^> O 00 =*fc —~ O N R% 8" S ^ S J S I fc -tq lq 8 -4 i - l <u tu Lq k, >*• "A Oo ON CN "vh ^- " N ^\ 'SS v\ (Jj (U tu Cu 8 R R R Its tu <u 187 3 a, O o c o O a i-o a ! H O Pi «' o "a, <u l H <u a o a QJ B QJ "QJ QJ QJ o ' £ o a, QJ 1-c QJ _a -a QJ -O ' C O QJ - d 3 a, o QJ H o 3 u "a. QJ a B u . o QJ QJ N <U QJ QJ L n v O l D C O ^ M ^ H O V D II II M a o J a a II n II n n n n X U 10 IS 10 H Q H O < fa O < rt > f> 2 S O H a a a a a c fa-u oi H CO cn c rt) CO "S OS CO CO H < fa H CO H a H H a EH a fa- OS fa fa a H S fa K s s u o fa H O o IS fa a fa Q OS o u fa- < u D a ta a s fa EH o OS H o CO a oi a H fa a u H H Pi w EH o a o; fa M a CO CO w 2 CO w 2 o H CO s H W fa o b g w a fa- CO Q fa w CO a OS o J fa O H H fa o o H a w OS a Q fa CO j CO H OS fa fa fa fa fa 0 fa fa o o o o o fa fa o O o p fa OS OS Qi OS OS a OS OS te- u W H u OS a M ca fa rn m m m m u < CQ to q a p 5 5 S O o 5 a a a a a o EH a a a W O Pi cn o oj ro in H CN CN CN CN CN ro in vo co in vo co m H H •tf m > co o H M ^ in > H H CN CN rO H CN ro -tf in in ro ro ro -tf ro ro o H ro ro vo s> CN CN ro ro cn o CN ro ro ro CO CTi CN CN ro ro in vo CN CN ro ro >t o o co co CN CN CN CN w Q O 188 in LT1 LO o o o o o o 1 1 1 1 i 1 IS IE w U> 2E 2E 7E 3E Q UJ CN H rH o o o o o o o o o o o o o o o 1 in <T\ CN o o o o o o o o o o o o o o o N cn r> O U) o o o o o o o o o o o o o o o CN UJ un CN o H U> UJ O UJ O O C O O O O O O O O O C X v O O O O O O O O C N O O O O O O o o m o o o o o o o o t n o o o o o o o o m o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o r o r o r o r o r o o o o o o co W D Q W o X o H O o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o O O H O O O O O O O O O H O O H H r-l H iH rH rH H H H U) CTl CN in CO Q rH H CN CN CN O 2 tn a 00 H r> o U J cn CN in CO -rH UJ r- r- CO H H CN CN CN H rH rH rH H ro ro ro ro ro ben to B S3 BJ 2 W U a. co ro ro ro ro co ro rH rH H o o o o o o o o O i i i i i i I 1 1 0 ] H w W < w w w H O o VO o ^ CN CN UJ ro ro CO in CO Q CT\ cn o o H O "tf in cn CN rH CN UJ rH 1 00 UJ ro CN r- H I rH CN ro CO cn ^ r> H VO CO ro ro cn • • H ro CO H CO in CO CN CO H H o rH o CN o r> cn CN CN in U) 00 in CN CO r- t- r- r- r- UJ H rH H 00 in ro rH in ro CO in ro m o o O o o o i i I i i i CO M w w w W w H O in o r> vo a o o cn H CN i cn m rH CN r> 00 X U) ro in H f- H cn 00 ao •vf1 cn CN ro t- ro 00 w H CN ro in a o 2 CO fc o H En U <C W • o o U i i < w w ~ W cn ro r> in in ai r-- ^ cn --i cn ro C N o X ro ro oo cn O H \f oo ro C N C N i o o ' o o i o o ' o o ' o o • o o W v o c n c N i n c o c O r H ^ r ^ QrHrHCNCNCNVor>r>r> O rH H rH rH 189 CN O W O O O CO o o o o o o> o o O O O CN O O o o o m o o o o o in o o o o o ro o o i-l r-1 rH H o o O O ta w ro CO CN o cn r> ro H VO cn ro ro CO LD cn 00 VO ro UO ro CO VO cn H VO o cn Lfl VD H cn H o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o o vo cn (N LD 00 CO H H CN CN CN H ro ro ro ro ro OUTPUT FILE 2 (SAMPLE ) Following is a typical output file providing the load - displacement curve for many replications DISPLACEMENT LOAD -.300000E-02 -.114300E+01 -.11G700E+01 -.119100E+01 -.282300E+01 -.284700E+01 -.287100E+01 -.289500E+01 -.291900E+01 -.294300E+01 -.296700E+01 -.299100E+01 -.301500E+01 -.303900E+01 -.306300E+01 -.308700E+01 -.311100E+01 -.313500E+01 -.315900E+01 DISPLACEMENT . 130979E+01 .413121E+03 .421791E+03 .430460E+03 .963515E+03 .967138E+03 .969484E+03 .969407E+03 .970101E+03 .970393E+03 .969730E+03 .965100E+03 .958432E+03 .949766E+03 .935536E+03 .917108E+03 .899986E+03 .884894E+03 .870838E+03 LOAD -.300000E-02 -.114300E+01 -.116700E+01 -.119100E+01 -.248700E+01 -.251100E+01 -.253500E+01 -.255900E+01 -.258300E+01 -.260700E+01 -.263100E+01 -.265500E+01 -.267900E+01 -.270300E+01 -.272700E+01 -.275100E+01 .115447E+01 .416110E+03 .424846E+03 .433582E+03 .853558E+03 .855299E+03 .852577E+03 .846192E+03 .838516E+03 .831G85E+03 .824157E+03 .816414E+03 . 806546E + 03 .794212E+03 .780463E+03 .766277E+03
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Computational modeling of strand-based wood composites...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Computational modeling of strand-based wood composites in bending Clouston, Peggi Lynn 2001
pdf
Page Metadata
Item Metadata
Title | Computational modeling of strand-based wood composites in bending |
Creator |
Clouston, Peggi Lynn |
Date Issued | 2001 |
Description | A stochastic finite element approach is presented herein for simulating the nonlinear behaviour of strand-based wood composites with strands of varying grain angle. The approach is based on the constitutive properties of the individual strands making it practical and versatile in application. It can be used to gauge the effects of varying strand characteristics in product manufacturing or it may be useful in the design of wood composite structures. Both 2 and 3 dimensional, stochastic, materially nonlinear finite element codes are developed. The nonlinear constitutive behaviour of the wood strands is characterized within the framework of rate-independent theory of orthotropic elasto-plasticity. The constitutive model employs the Tsai-Wu yield criterion with the associated flow rule of plasticity. Failure is marked by an upper bound surface whereupon either perfect plasticity or an abrupt loss of strength and stiffness ensues. This constitutive model is implemented into the finite element codes. The programs are based on the conventional displacement formulation using linear isoparametric elements. The nonlinearities are resolved by a modified Newton-Raphson procedure. The programs are further formulated in a probabilistic manner using random variables as input for principal material strengths and stiffnesses. To generate entire data samples for comparison with experimental samples, the programs were written with extended capacity to perform Monte Carlo simulations. The mechanical properties of the strands are derived through both experiment and analysis. An experimental database of principal material strengths and stiffness parameters is acquired for Douglas fir heartwood strands. Statistical parameters for shear strength and stiffness as well as the interaction parameter of the Tsai-Wu criterion are estimated, however, through a least square minirriization of error between simulated and experimental compression strength of [+15]s angle-ply laminates. Weibull weakest-link theory is employed to adjust experimental tensile strength values for size effect. The general performance of the programs is verified through comparison of results for several analyses solved using analytical techniques or alternate programs. Following this, the ability of the models to reproduce experimental findings for angle-ply laminates in tension, compression and 3 point bending is validated. A preliminary investigation is conducted to compare numerical simulations with experimental data for Parallam* in tension and 3 point bending. The favourable comparisons of the model to experimental results attest to the effectiveness of the proposed technique. |
Extent | 7001482 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2018-05-08 |
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. |
IsShownAt | 10.14288/1.0075230 |
URI | http://hdl.handle.net/2429/13671 |
Degree |
Doctor of Philosophy - PhD |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2001-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2001-714500.pdf [ 6.68MB ]
- Metadata
- JSON: 831-1.0075230.json
- JSON-LD: 831-1.0075230-ld.json
- RDF/XML (Pretty): 831-1.0075230-rdf.xml
- RDF/JSON: 831-1.0075230-rdf.json
- Turtle: 831-1.0075230-turtle.txt
- N-Triples: 831-1.0075230-rdf-ntriples.txt
- Original Record: 831-1.0075230-source.json
- Full Text
- 831-1.0075230-fulltext.txt
- Citation
- 831-1.0075230.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-0075230/manifest