UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 MESHLESS LOCAL PETROVGALERKIN MICROMECHANICAL ANALYSIS OF STRUCTURAL COMPOSITES By THI D. DANG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2007 PAGE 2 2 2007 Thi D. Dang PAGE 3 3 To my parents, my wife and son for all the hard work, love and sacrifices they made toward my education PAGE 4 4 ACKNOWLEDGMENTS This dissertation marks the end of a journey that started in 2003 at the Department of Mechanical and Aerospace Engineering of the Univ ersity of Florida. I would like to thank the following people for their s upport in this endeavor. First, I would like to express my deep thanks to Dr. Bhavani V. Sankar for his exceptional guidance and support over the years. I will foreve r be indebted to him. He took me under his wing, took care of me and gave me opportunities to grow and mature taught me what I needed to know, and helped me make it. Here, I have b een extremely fortunate to have him as my academic advisor, he is not only an expert at the very forefront of the stateoftheart in composites science and technology, but also he is si mply the best advisor. I would like to thank him so much for passing knowledge on to me, a nd teaching me how to do research and solve new problems with perseverance. I learned much fr om him. I thank him very much for all what he did for me during my time here at the University of Florida. I am thankful to Dr. NamHo Kim, a co mmittee member of mine, for his time and for having taught the interesting course of Advanced Finite Element, I am very grateful to him for many valuable comments and useful discussions. I would like to express my grateful thanks to Dr. Ashok V. Kumar and Dr. William W. Hage r for their time and effort serving on my committee. Their valuable advice has improved the quality of this dissert ation. I would like to extend my special thanks to Dr. Wei Shyy, for his he lp and advice in the difficult time of mine at the University of Florida. The financial support of the Vietnam Edu cation Foundation during my Ph. D study is gratefully acknowledged. In particular, I am extr emely thankful to Mr. Kien Pham, the former executive director of Vietnam Education Founda tion, for his help and encouragement. I will never forget his help and kindness. PAGE 5 5 Special thanks go out to my friends and collea gues, Dr. HsuKuang Ching for the valuable comments and enlightening discussions on mesh less methods, Dr. Ryan Karkkainen and Dr. Satish Bapanapalli for their help on ABAQUS, and others at the Center for Advanced Composites of the University of Fl orida for their help in various wa ys. I also wish to thank all of my friends both near and far w ho have always wished me well and encouraged me through out. I would like to express my d eepest appreciation to my pare nts, wife, and soontobeborn child for all the hard work, love, and sacrifices they made toward my education. This dissertation is dedicated to them. PAGE 6 6 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........8 LIST OF FIGURES................................................................................................................ .........9 ABSTRACT....................................................................................................................... ............12 CHAPTER 1 INTRODUCTION..................................................................................................................14 Composite Materials and Computational Methods................................................................14 Objective and Outline of Study..............................................................................................19 2 MESHLESS LOCAL PETROVGALERKIN METHOD.....................................................23 Summary........................................................................................................................ .........23 Background..................................................................................................................... ........23 Moving Least Squares Approximation...................................................................................25 Weight Function................................................................................................................ .....28 Meshless Local PetrovGalerkin Weak Formulation............................................................29 Conclusions.................................................................................................................... .........33 3 MESHLESS LOCAL PETROVGALERKIN MICROMECHANICAL ANALYSIS..........37 Summary........................................................................................................................ .........37 Introduction................................................................................................................... ..........37 Unit Cell Analysis for Three Dimensional Elastic Constants................................................42 Meshless Local PetrovGalerkin Formulation of the Problem...............................................44 Moving Least Squares Approximation Scheme..............................................................44 Meshless Local PetrovGalerkin Weak Formulation and Discretization........................46 InPlane Compression, Tension and Shear Cases...........................................................52 Generalized Plane Strain Problem...................................................................................52 Out of Plane Shear Test: Special Generalized Plane Strain Problem..............................55 Treatment of Material Discontinuity...............................................................................57 Treatment of Periodic Boundary Conditions...................................................................61 Conclusions.................................................................................................................... .........63 4 STIFFNESS AND SHEAR CORRECTION FACTOR PREDICTION OF A CLASS of COMPOSITE BEAMS...........................................................................................................67 Summary........................................................................................................................ .........67 Background..................................................................................................................... ........67 PAGE 7 7 Exact Solution for the Shear Correction Factor of Laminated Beams...................................71 Shear Correction Factor of Isotropic Beams...................................................................71 Shear Correction Factor of Laminated Beams................................................................71 Micromechanical Analysis of Thin Walled Textile Composite Beams.................................79 Geometric Model of Textile Composite Unit Cell..........................................................79 Constitutive Relation of Textile Composite Unit Cell....................................................80 Unit Cell Boundary Conditions.......................................................................................81 Periodic boundary condition for unit curvature.......................................................82 Periodic boundary condition for un it transverse shear strain...................................84 Conclusions.................................................................................................................... .........86 5 RESULTS AND DISCUSSION.............................................................................................92 Introduction................................................................................................................... ..........92 Timoshenko Beam Problem...................................................................................................92 TwoPhase Verification Problem...........................................................................................93 Prediction of Stiffness Properties...........................................................................................95 Example 1...................................................................................................................... ..95 Example 2...................................................................................................................... ..97 Example 3...................................................................................................................... ..99 Stiffness and Shear Correction Factor Prediction of Composite Beams................................99 Conclusions.................................................................................................................... .......102 6 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK........................124 Summary........................................................................................................................ .......124 Major Conclusions and Contributions..................................................................................125 Recommendations for Further Work....................................................................................126 APPENDIX IMPORTANT NOTES FOR USI NG HALPINTSAI EQUATIONS...............128 LIST OF REFERENCES.............................................................................................................130 BIOGRAPHICAL SKETCH.......................................................................................................137 PAGE 8 8 LIST OF TABLES Table page 31 Periodic boundary conditions for the MLPG method........................................................64 41 Geometric parameters of plai nweave textile beam unit cell.............................................88 42 Periodic displacement boundary conditions for the beam unit cell...................................88 51 Material properties of Eglass/PP composite...................................................................103 52 Comparison of elastic constants between the current method and HalpinTsai equations...................................................................................................................... ....103 53 Material properties of Eglass/epoxy composite..............................................................103 54 Comparison of elastic constants betw een the current method and HalpinTsai equations...................................................................................................................... ....103 55 Material properties of Eglass/PP composite...................................................................103 56 Comparison of elastic constants.......................................................................................104 57 Comparison of shear correction fact ors for various types of laminates...........................104 58 Geometric parameters of plai nweave textile beam unit cell...........................................104 59 Constituent material proper ties for beam examples.........................................................104 510 Comparison of beam stiffness coefficients......................................................................105 PAGE 9 9 LIST OF FIGURES Figure page 11 Schematics of woven composites......................................................................................21 12 Illustration of possible unit cell ge ometries for the model material..................................21 13 Alternativ e unit cells..................................................................................................... .....22 14 Finite element model of the uni t cell for structural composites.........................................22 15 Mesh of yarn reinforcements in a unit cell for textile composites.....................................22 21 Illustration of do main of influence....................................................................................34 22 Distinction between andiiuu ............................................................................................34 23 Shape of weight function...................................................................................................35 24 Schematic representation of a body in MLPG method......................................................35 25 Schematics of MLPG.........................................................................................................36 31 Unit cell of a unidirectional composite..............................................................................64 32 Illustrative example of an inhomogeneous body including multipoint constraints...........65 33 Longitudinal shear problem...............................................................................................65 34 Illustration of an inhomogeneous body..............................................................................66 35 Unit cell and periodic boundary conditions for the deformation121M...........................66 41 Coordinate system for the laminated beam........................................................................88 42 Textile composite b eam under shear force........................................................................89 43 Unit cell of thin walled textile composite beam................................................................90 44 Half unit cell model of thin walled textile composite beam..............................................90 45 Cantilever beam............................................................................................................ .....91 51 Timoshenko beam............................................................................................................105 52 Nodal scheme with 72 nodes...........................................................................................105 53 Comparison of displacements U & V at the line y = 0.5...................................................106 PAGE 10 10 54 Comparison of normal stresses x at the line y= 0.5........................................................106 55 Comparison of normal stresses x at the line x = 4...........................................................107 56 Comparison of shear stresses xy at the line x = 4.............................................................107 57 Two phase problem..........................................................................................................108 58 Comparison of radial disp lacements at the interface.......................................................109 59 Comparison of horizontal displacements along line y =0................................................109 510 Comparison of radial displacements along line x =1.......................................................110 511 Radial stresses in the two materials at the interface by MLPGbased micromech anics...............................................................................................................110 512 Hoop stresses in the two materials at the interface by MLPGbased micromechanics...111 513 Tangential (shear) stresses in the two materials at the interface by MLPGbased micromech anics...............................................................................................................111 514 Comparison of radial stresses at the interface..................................................................112 515 Comparison of tangential stresses at the interface...........................................................112 516 Distribution of x ............................................................................................................113 517 Distribution of y ...........................................................................................................113 518 Distribution of z ............................................................................................................113 519 Distribution of x y ...........................................................................................................113 520 Distribution of interf acial stresses for the case11 =1.......................................................114 521 Distribution of x ............................................................................................................114 522 Distribution of y ...........................................................................................................114 523 Distribution of x y ...........................................................................................................115 524 Distribution of z ............................................................................................................115 525 Distribution of interf acial stresses for the case331 .....................................................115 PAGE 11 11 526 Displacement along zdirection over domain for the case131 ....................................116 527 Unit cell for shear test (121M ) in the x1x2 plane.........................................................116 528 Deformed boundary of unit cell for case of121M .........................................................117 529 Deformed unit cell for case of121M and shear stress distribution by FEM.................118 530 Periodic deformation of unit cell for shear test (121M ) by the current method............118 531 Tension tests in x1 and x2 directions (111M and221M )...............................................119 532 Comparison of transverse elastic modulus ET of the composite......................................119 533 Comparison of shear modulus GLT of the composite.......................................................120 534 Unit cell of thin walled textile composite beam..............................................................120 535 Half unit cell model of thin walled textile composite beam............................................121 536 Deformed unit cell on shear problems.............................................................................121 537 Deformed unit cell........................................................................................................ ...122 538 Textile composite beam under bending...........................................................................123 PAGE 12 12 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MESHLESS LOCAL PETROVGALERKIN MICROMECHANICAL ANALYSIS OF STRUCTURAL COMPOSITES By Thi D. Dang August 2007 Chair: Bhavani V. Sankar Major: Mechanical Engineering This study aims at the development of a micr omechanical model for structural composites using the meshless local PetrovGa lerkin method (MLPG) to predic t the stiffness properties, and the shear correction factors of structural compos ites from the analysis of the representative volume element (RVE). Micromechanical analysis of structural composites is possible due to the presence of the RVE also called the unit cell. On the microscale, comparable to the scale of unit cell, the composite is heterogeneous due to the pr esence of the reinforcing yarn and the matrix. However, on the macroscale, compared to the structural scale, the composite is assumed to be homogeneous and orthotropic. Th e homogeneous composite proper ties are then predicted from the properties of the constituent materials and their distribution. The highlight of the study is the treatment of some challenging problems in the meshless local PetrovGalerkin me thodbased micromechanical analysis of structural composites. The unit cell is discretized and periodi c boundary conditions are set up and imposed between opposite endfaces of the unit cell, essential boundary co nditions are enforced by the penalty method. The treatment of material discontinu ities at the interfaces between different phases of the composite is presented by means of the direct imposition of interface boundary conditions. In addition, an algorithm for handling of periodic boundary conditions in the MLPG method using the PAGE 13 13 multipoint constraint technique is also presente d. The MLPG formulation is presented for the six linearly independent deformations of the unit cell. From the forces acting on the unit cell for each of the six deformations, the elastic cons tants of the composite can be computed. In this study, a micromechanical model usi ng both the finite element method and the MLPG method is suggested to predict the flexur al stiffness properties and the shear correction factors for a class of composite beam problems including textile composite beams. Examples are presented to illustrate the effec tiveness of the current method, and it is validated by comparing the results with available anal ytical and numerical solutions. The current method shows great potential for applications in micromechan ics, especially for textile composites. PAGE 14 14 CHAPTER 1 INTRODUCTION Composite Materials and Computational Methods Composites can be generally defined a select ed combination of dissimilar materials with a specific internal structure. The unique combination of materi al components leads to singular mechanical properties and superi or performance characteristics not possible with any of the components alone. Additionally composite material s are often overwhelmingl y superior to other materials on a strengthtoweight or stiffnesstoweight basis. With this advantage, the range of applications for composite mate rial appears to be limitless. Due to their high specific stiffness and stre ngth, advanced fiber reinforced composite materials have long been used in the aerospace in dustry, and with the increasing focus on lightweight vehicle manufacture due to environmental legislation, automotive applications are becoming more widespread. Other notable engineer ing applications includ e pressure vessels and waste water pipes and fittings. Th e need to utilize the mechanical performance of materials and to avoid unreasonable overspecification for ae rospace applications was highlighted by Boeing, who estimated that it cost $10,000 per pound to la unch a satellite into orbit (Boeing Press Release, 1996). These composite materials are not only used in the aerospace and automobile industry but also are widely used in marine and biomedical applications. Recent developments in textile manufacturi ng processes show many superior properties, for example the strength and stiffness of the comp osite material can be oriented in required directions, and no material we ight is wasted in providing reinforcement in unnecessary directions. Textile composites pr esenting a class of advanced composites can be defined as the combination of resin system with a textile fi ber, yarn or fabric system (Figure 11). PAGE 15 15 Textile structural composites have been the focus on many aircraft and automotive manufacturers for more than two decades. They offer a good alternative to develop composite primary structures for commercial transport airplanes with costs that are competitive with those of current metallic airplanes. Textile com posites were considered for many components to improve structural performance and to reduce costs. Boeing and LockheedMartin evaluated textile composites for fuselage frames, wi ndow belt reinforcements, and various keel components of the fuselage. Northrop Grumman ev aluated textile concepts for making stiffened skins using 3D textile composites, and McD onnell Douglas evaluated knitted, braided, and stitched textile fabrics for a wing box. A large part of the early investigations on textile structural composite application were c onducted within th e scope of the textile working group within NASAs Advanced Composites Technology (A CT) program (Poe and Harris, 1995). The heterogeneity and anisotr opy of structural composites ma ke the application of the standard analysis methods unsuitable. As a resu lt of this increasing in terest in structural composites, a number of micromechanical models to predict effective strength and stiffness of such composites have been developed (Is hikawa and Chou, 1982; Chou and Ko, 1989; Whitcomb, 1991; Naik, 1994; Marrey and Sankar, 1995; Poe et al, 1997; Cox and Flanagan, 1997). Most of these models, even though they ar e based on rather very simplifying assumptions, provide a fairly good approximation of effectiv e mechanical properties. The conventional micromechanical models for textile structural co mposites assume that the state of stress is uniform over a distance comparable to the dime nsions of the represen tative volume element (RVE) or the unit cell, this woul d require the prediction of the e ffective (macroscopic) properties of the composites from the constituent material (microscopic) characteristics such as yarn and matrix properties, yarn matrix interface character istics and the yarn architecture (Figure 12). PAGE 16 16 Most of all, programming of these models provide engineers with powerful tools to compute effective properties for use in FEM software and a costeffective alternative to extensive testing. This is possible if we assume that there is a representative volume element or a unit cell that repeats its self throughout the volume of the composite. The unit cell can be considered as the smallest possible building block for the text ile composite such that the composite can be constructed from spatially transl ated copies of it, w ithout the use of rota tions or reflections (Figure 13), and the unit cell is not unique. The prediction of the effective macroscopic properties from the constituent mate rial characteristics is one of the aspects of the science known as micromechanics. The effective properties incl ude thermomechanical properties like stiffness, strength and coefficients of the thermal expansion as well as thermal conductivities, electromagnetic and othe r transport properties. Most of the early work to determine the properties of heterogeneous materials was restricted to particular (spherica l) inclusions, which were assumed to be isotropic. For example, Hashin (1962) derived expressions for the bounds for the elastic m oduli of heterogeneous materials using variational theory A review of analytical met hods for predicting the effective properties of particulate composites was pres ented by Christensen (1990). However, the complex geometry of the textile performs make s such precise analytic al modeling difficult. Traditionally, the micromechanical analysis fo r textile composites can be classified into three categories: mechanics of material models energy method, and finite element of the unit cell. All of the above models use a unit cell in th e composite material, and they attempt to model the material as homogeneous orth otropic material. In the mechanics of material models, the yarns are approximated as simple structural elemen ts, such as beams, plates, laminates, etc., and their deformation behavior is assumed to be governed by the corresponding structural PAGE 17 17 constitutive relations. The kinematics is also si mplified to a great extent, and a relation between the overall deformation of the unit cell and the aver age forces are derived. The energy approach is similar to the previous one, except that stra in energy in the unit cell is evaluated based on some assumed displacement field, which is usua lly an oversimplification of exact displacements. The elastic constants are derived by equating the strain energy in the approximate model to that in an idealized homogeneous composite. Most en ergy based approaches provide bounds for the homogeneous properties; can be used as a ch eck for experimental ob servations or other theoretical models. The th ird method is the rigorous micromech anical analysis of the unit cell, which often requires the use of numerical methods like FEM (Fi gure 14), and uses the exact threedimensional constitutive relations for the yarn and the matrix material. The finite element method (FEM) has been the dominant technique in computational mechanics in the past decades; it has made signi ficant contributions to the developments in engineering and science. The finite element method has been successfully applied to many problems in mechanics of composites material s. It is a robust and thoroughly developed technique, but it is not without shortcomings. The reliance of the method on a mesh leads to complications for certain classes of problems; fo r example, crack proble ms with arbitrary and complex paths, and simulation of phase transforma tions is also difficult. In the development of advanced composite materials, especially textile composites, one of the major technical barriers in modeling textile composites such as braided an d woven composites is the finite element mesh generation (Marrey and Sankar., 1995, 1997). Actu ally, for composite materials with complex yarn architectures, the meshing of individual yarns in the unitcell is quite simple. However, the meshing of interfacial region be tween matrix phases and individua l yarns is much more difficult (Kim and Swan, 2003). The region is multiply connected, and hence mesh in different phases PAGE 18 18 may not be compatible (Figure 14B and Figure 15), furthermore it is difficult to get a suitable mesh on which opposite faces of th e unit cell have identical nodes so that periodic boundary conditions can be implemented using multipoint c onstraints, and it is impossible to completely overcome those meshrelated difficulties by a mesh based method. The available finite element based methods are satisfactory for stiffness pred iction because stiffness properties are based on volume averaging of stresses and strains in the representative volume element of the composite, and the approximation involved in the FE meshi ng does not affect the results significantly. However modeling the damage, especially progre ssive damage, requires accurate description of the stress field in different phases and requires a very fine mesh (Marrey and Sankar, 1997; Zhu et al., 1998). The FEM based micromechanical models have been successfully employed in predicting thermoelastic constants of fiber reinforced composites materials; their use for strength prediction under multiaxial loading co nditions is not practical (Zhu et al, 1998, Karkkainen and Sankar, 2006). We expect that the tediousness and inaccura cies involved in mesh generation, and hence inaccuracies in the results can be avoided us ing the new meshless techniques such as the meshless local PetrovGalerkin (MLPG) met hod. The MLPG approach proposed by Atluri (1998, 1999, 2000) is one of the several meshless schemes. The main advantage of this method compared to other meshless methods is that no ba ckground mesh is used to evaluate the various integrals appearing in th e local weak formulation of the pr oblem. Therefore, this method is a truly meshless approach in terms of both interpolation of variable s and integration of energy. The meshless methods have been demonstrated to be efficient in solving di fferent problems (Atluri and Zhu, 1998; Atluri, 2004; Raju and Chen, 2001; Raju and Phillips, 2003; Ching and Batra, PAGE 19 19 2001; Batra and Ching, 2002; Batra et al., 2004). In this study, the MLPG me thod is applied to micromechanics problems of composite materials. Objective and Outline of Study The objective of our study is to present the meshless local Petrov Galerkin (MLPG) formulation for problems of composite microm echanics, called the meshless local Petrov Galerkin methodbased micromechanical model fr om the analysis of uni t cells. The theoretical methodologies, computer implementations and pract ical application are ca rried out. The MLPG micromechanical model is proposed, firstly to predic t the elastic constants/s tiffness properties of structural composites, and sec ondly to predict the flexural s tiffness properties and the shear correction factor for a class of composite beam problems including textile composite beams. For the above objective, in this study periodic boundary conditions of the unit cell under various loadings are set up. Especially we suggest a technique to treat mate rial discontinuity at interface, while periodic boundary conditions are handled us ing the multipoint constraint technique. The MLPG method is formulated for generalized plane strain problems. Examples are presented to illustrate the effectiveness of the current micromechanical model, and it is validated by comparing the results with available analytical and numerical solutions. The general outline of this study is as follows: In Chapter 2, the literature on meshless methods is reviewed, with emphasis on the meshless local Petrov Galerkin method which is us ed through our study. The description of the moving least squares approximation, the choice of the weight func tions and the strong and weak forms, and its discretization ar e also given in this chapter Chapter 3 gives a theoretical methodology to com pute the macroscale elastic constants, in this chapter the meshless local Petrov Galerkin formulation is presented for each of the six marcrostrains imposed the unit cell. Periodic boundary conditions for each case of macrostrains PAGE 20 20 are set up. Problems arising in micromechanics ar e treated such as the treatment of material discontinuity at the interface and periodic boundar y conditions are handled using the multipoint constraint technique. In Chapter 4, based on our study in Chapter 3, a theoretical framew ork and procedure are given to predict stiffness prope rties and shear correction factor s for composite beams including textile composite beams. The suggested model is straightforward and a truly micromechanical model. In addition, an exact solution for the shea r correction factor of general laminated beams is also derived in this chapter. In Chapter 5, Computer implementations and numerical examples are treated. The examples are presented to illustrate the effectiveness of the current method, discussion on results for each example is given, and it is validated by comparing the results obtained from the current micromechanical model with available analyti cal and numerical soluti ons. Finally Chapter 6 summarizes the work done in our study, in whic h conclusions and contribu tions are pointed out, and some recommendations for further possible deve lopments are also mentioned briefly in this chapter. PAGE 21 21 Figure 11. Schematics of woven co mposites (without matrix pockets). Figure 12. Illustration of possi ble unit cell geometries for the model material. A) Homogeneous and orthotropic composite (macrosca le). B) Composite (microscale). y2 y1 x2 x1 A B PAGE 22 22 Figure 13. Alternative unit cells (short and l ong dash lines at upper le ft) in a plane woven laminate. The smaller, solidlined rectangl e (lower right) shows a reduced cell that takes advantage of reflection symmetry a bout a vertical plane. [Reprinted with permission from Cox, 1997. (Page 57, Figure 53)]. A B Figure 14. Finite element model of the unit cel l for structural composites. A) Periodically deformed unit cell of fiber composites. B) Unit cell of textile composites. Figure 15. Mesh of yarn reinforcements in a un it cell for textile composites. A) Voxelbased mesh. B) Tetrahedral mesh. [Reprinted with permission from Kim and Swan, 2003. (Page 1003, Figure 19)]. PAGE 23 23 CHAPTER 2 MESHLESS LOCAL PETROVGALERKIN METHOD Summary In this chapter, we give the literature review on meshless methods. The meshless local PetrovGalerkin (MLPG) method is used to solve micromechanics pr oblems of composite materials. The MLPG formula tion including the movi ng least squares approximation, the choice of the weight functions and the strong and weak forms are also introduced. Background The initial idea of meshless methods can be traced back to 1977 when Lucy (1977) and Gingold and Monaghan (1977) proposed a smooth particle hydrodynamics (SPH) method that was used for modeling astrophysi cal phenomena without boundaries, such as exploding stars and dust clouds. Extensive developments have been ma de in several varietie s since then and with many different names: SPH (Monaghan, 1982, 1988,1992), Generalized Finite Difference Method (Liszka and Orkisz, 1980), Diffuse Element Method (Nayrole s et al., 1992), Particle in Cell Method (Sulsky et al ., 1992), Wavelet Galerkin Method (Qian and Weiss, 1993), Reproducing Kernel Particle Met hod (RKPM) (Liu et al., 1995a,b ) Element Free Galerkin (EFG) (Belytschko et al., 1994), Parti tion of Unity (PU) (Babuska and Melenk, 1995, 1996), Hp Cloud Method (Duarte and Oden, 1995,1996) Finite Point Method (Onate et al., 1996a,b), Natural Element Method (NEM) by Sukumar, Moran and Belytschko (1998), Meshless Galerkin Methods using Radial Basis Functions (RBF) (W endland, 1999). The major differences in these socalled meshless methods come only from th e techniques used for interpolating the trial function. However, even though no mesh in required in these methods for th e interpolation of the trial and test functions for the solution variables, the use of background meshes is inevitable in PAGE 24 24 these methods, for the integration of the weak form, or of the en ergy. Therefore, these methods are not truly meshless. Recently, Atluri et al have developed two truly meshless methods, those are the meshless local boundary integral equation (LBIE) (Zhu et al., 1998a, b), a nd the meshless local PetrovGalerkin (MLPG) method (Atluri and Zhu, 1998a, b) Both these methods are truly meshless, as no domain or boundary meshes are required in th e two methods, either for the purposes of interpolation of the trial and test functions fo r the solution variables, or for the purpose of integration of the weak form. All integrals can be easily evaluated over overlapping, regularly shaped, domains (in general, spheres in th reedimensional problems) and their boundaries. The local weak forms are the basis of the meshless local PetrovGa lerkin (MLPG) method (Atluri, 2004). Methods that base d on the local unsymmetric weak formulation are the methods MLPG2, MLPG3, and MLPG4 (LBI E). MLPG2 uses the collocati on Diracs Delta function as the test function in each subdom ain, while MLPG3 employs the error function as the test function in each subdomain, MLPG4 (LBIE) is developed by employing the modified fundamental solution as the test function. Methods that based on the local symmetric weak form (LSWF) are the methods MLPG1, MLPG5 and MLPG6. MLPG1 is developed by employing the weight function as the test function in each loca l subdomain. In MLPG5, a constant function for second order PDE is taken to be the test functi on in each subdomain. In MLPG6, the trial and test function come from the same space, this MLPG6 reduces to the Galerkin method and leads to a symmetric stiffness matrix. It should be noted that in all the MLPG methods, the local weak forms are constructed over overlapping subdomains of arbitrary shapes, such as spheres, tetrahedral, hexahedra, etc. in 3D. In as much as the local subdomains are overlapping, spheres are convenient to be used as local subdomai ns because of their directional isotropy. PAGE 25 25 In this study, we use the ML PG1 method to solve micromechan ics problems of composite materials. For simplicity, this method will be ge nerally termed as the meshless local PetrovGalerkin (MLPG) method in this dissertation. The MLPG formulation including the moving least squares approximation, the choice of the weight function, strong form, and the local symmetric weak form (LSWF) are presented. A computer code in MATLAB is developed and a classical example is given to validat e the code and the method. Moving Least Squares Approximation In the MLPG method, the field variable u(x) is approximated by the moving least squares (MLS) technique. This approximation is based on three components: a weight function of compact support associated with each node, polynomia l basis functions, and a set of coefficients that depend on the position x of the point. First we consider a subdomain x called the domain of de finition of the MLS approximation for the trial function at the point x, which is located in the problem domain The unknown trial approximation uh(x) of the function u(x) is defined by 1()()(),m h jjx jupxaxTxp(x)a(x)x (21) where 12()(),(),...,()T mppp pxxxx is a vector of the comple te monomial basis of order m, and a(x) is a vector containing unkn own coefficients(),1,2,...,jajm x Examples of ()T x p in the 2D problems are 12[1,,]forlinearbasis,3xxmTP(x) (22) 22 121122()[1,,,,,]forquadraticbasis,6xxxxxxmTPx (23) PAGE 26 26 The m unknown parameters ()jax can be determined by minimi zing the weighted discrete 2L norm, defined as 2 1 ()n iii iJwu Tx(x)p(x)a(x)(24) where wi(x) is the weight function as sociated with the node i with wi(x) >0 for all x in the support of wi( x ) xi denotes the value of x at node i n is the number of point s in the neighborhood of x ( x) for which the weight functions wi( x )>0 (Figure 21), and iu refers the fictitious nodal value and not the nodal valu e of the function uh at the pointi x Finding the minimum of J ( x ) in Equation 24 with respect to a ( x ) leads to the following linear relation ()()() Axax=Bxu (25) where 1()()()()n iii iwTAxxpxpx (26) 1122()()(),()(),...,()()nnwpwpwp Bxxxxxxx (27) 12 ,,...,T nuuu u (28) Equation 25 can be solved for a ( x ) and yields 1 ()()() ax=AxBxu (29) As can be seen from Equation 29, the unknown coefficients a ( x ) can be obtained only if A ( x ) defined by Equation 25 is nonsingular, so a necessary condition for a welldefined MLS approximation is that at least m weight functions are nonzero (i.e. n m ) for each sample point x. PAGE 27 27 1 ()(),(),n hh iiiiix iuuuxuu xxx (210) where 1()()()()m ij ji jp 1xxAxBx (211) ()i x is usually called the shape function of the MLS approxi mation corresponding to node i The spatial derivatives of the shape function ()i x are 111 ,,,, 1()()m ikjkjijkkji jpp ABABAB (212) where,/i i x The derivative of A1 in Equation 212 can be computed by taking the derivative of A1A = I Thus 111 ,, kk AAAA (213) It should be noted that the shape functions()ij x derived from the MLS approximation do not satisfy the Kronecker delta criterion()ijij x Therefore, they are not nodal interpolants and the name approximation is used, i.e. ()h iiuu x (Figure 22 for a simple one dimensional case for the distinction between andiiuu ). This property makes the satisfaction of essential boundary conditions more difficult than that in the finite element method. Several techniques have been developed to enfor ce essential boundary conditions, incl uding the method of Lagrange multipliers (Belytschko et al., 1994), modified variational principl es (Lu et al., 1994), and the penalty method (Atluri and Zhu, 1998). In the current study, we use the penalty method to enforce the essential boundary conditions. PAGE 28 28 Weight Function An important ingredient in the meshless met hod is the weight function used in Equations from 24 to 213, the weight function shoul d be nonzero only over a small neighborhood of xi, called the domain of influence of node i, in order to generate a set of sparse discrete equations. This is equivalent to starting that the weight functions have compact support. The precise character of wi( x ) seems to be unimportant, although it is almost mandatory that it be positive and increase monotonically as i xxdecrease. Furthermore, it is desirable that wi( x ) be smooth, if wi( x ) is C1 continuous, then for polynomial basis, the shape functions ()i x are also C1 continuous (Dolbow and Belytschko, 1998). In the current study, we consider the two we ight functions which are commonly used in meshless methods: one is the Gaussian weight function and the other is the spline weight function. They have compact support and are con tinuously differentiable f unctions (Figure 23). The Gaussian weight function and the spline we ight function are given, respectively, in Equation 214 and Equation 215 below. 22 2exp(/)exp(/) 0 () 1exp(/) 0kk iiii ii k i ii iidcrc dr w rc dr x (214) 23416830 () 0iii ii i iii iiddd dr w rrr dr x (215) where iidxx is distance from the sampling point x to the node xi, and ri is the radius of the domain of influence of the weight function wi( x ). The parameters ci and k in Equation 214 control the shape of the Gaussian weight function wi( x ). The parameter k = 1can be taken in Equation 214. So far, there is no theory to determine an optimal value of the parameter ci and is PAGE 29 29 chosen empirically. Lu et al (1994) reco mmended a method to choose the constant ci. We follow their suggestion and define ci as the distance from the node xi to the third nearest neighboring node. The domain of influence ri can be chosen as 3.5iirc so that the weight function wi( x ) covers sufficient number of nodes to ensure the n onsingularity of A in Equation 211. It has been reported in the literature th at the Gaussian weight function usually performs better than the spline weight function. Therefore, unless otherwise stated, the Gau ssian weight function is used in the current study. Meshless Local PetrovGalerkin Weak Formulation We consider the following twodimens ional elastostatic problem on the domain bounded the boundary (Figure 24). Consider a 2D elasticity problem in the domain bounded by The equilibrium equations are ,0inijjib (216) where ij is the Cauchy stress tensor, and bi is the body force. The boundary conditions are as follows onijjitnt (217) oniiuuu (218) where it is defined as the prescrib ed traction on a surface, and iu is the prescribed displacement field, and nj is the unit outward normal to the boundary u and t are complementary subsets of Unlike the element free Galerkin method (B elytschko et al., 1994) wh ich is based on the global Galerkin formulation, the present local Pe trovGalerkin formulation is constructed over a PAGE 30 30 local subdomain s which is located inside the global domain This local subdomain s is taken to be either a circle or a part of a circle in a 2D problem. A generalized local weak form of Equati ons from 216 to 218 over a local subdomain s can be written as follows ,()()0ssuijjiiiiibvduuvd (219) where ui and vi are the trial and the test functions, respectively, and su is the part of the boundary s of s over which essential boundary cond itions are specified. In general, s ssL with s being the part of the local boundary located on the global boundary and Ls being the other part of the local boundary ove r which no boundary conditions are specified, i.e., ss with s ssL (Figure 25). In Equation 219, is a penalty parameter ( >>Youngs modulus/ Length) with is used to im pose the essential boundary conditions. In our study, the value of791010 gives good results. Using ,,,()ijjiijijijijvvv and the divergence theorem in Equation 219 leads to ,()()0sssuijjiijijiiiiinvdvbvduuvd (220) where ni is a unit outward normal to the boundary s Unlike in the conventional Galerkin method in which the trial and test func tions are chosen from the same space, the PetrovGalerkin method uses the trial functions and the test functions from different spaces. In particular, the test functions need not vani sh on the boundary where the essential boundary conditions are specified. Imposing the natural boundary condition stoniijjitnt in Equation 220 we obtain ,()()0ssustssuiiiiiiijijiiiii Ltvdtvdtvdvbvduuvd (221) PAGE 31 31 in which st is part of s over which the natural boundary condition, iitt is specified. For a subdomain located entirely within the global domain, there is no intersection between s and s sL and the integrals over su and st vanish. In order to simplify the above equation, we deliberately select the test functions vi are chosen such that they vanish on Ls, the circle (for an internal node ) or the circular arc (for a node on the global boundary ). This can be accomplished by using the weight function wi(x) in the MLS approximation as also the test function vi (Atluri and Zhu, 1998b), but the radius ri of the support of the weight function is replaced by the radius ro of the local domain s, such that the test function vanishes on Ls. Using these test functions and rearranging Equation 221, we obtain the following local symmetric weak form in linear elasticity, as ,ssusustssuijijiiiiiiiiiivduvdtvdtvdbvduvd (222) For 2D problems, two independent sets of test functions should be applied in Equation 222, which gives ,ssusustssuijkijikiikiikiikiikivduvdtvdtvdbvduvd (223) where vki denotes the ith component of the test function in the kth set. For simplicity, Equation 223 can be written form as ssusustsusdadddaddv +vuvt=vt+vu+vb (224) In Equation 224, v denotes the strain matrix deri ved from the test functions, and is the stress vector derived from the trial functions. That is 11 (1)(1)(1) 112212 22 (2)(2)(2) 112212 12,v (225) PAGE 32 32 where the superscript i denotes the ith test function. Functions v, u, t, and b are defined as follows: 1112111 2122222,,, vvutb vvutb vutb (226) The two sets of test functions v in Equation 226 should be linearly independent. The simplest choice for v is or ijijvvv vI where ij is the Kronecker delta and I is the identity matrix. As long as the union of all local subdomains covers the global domain, Equations from 216 to 218 will be satisfied in the global domain and on its boundary respectively. Substituting the MLS approximation Equation 210 into Equation 224, and summing over all nodes leads to the following discre tized system of linear equations 111 ()()() ()()()ssusu stsusnnn vijjijjijj jjj iiidadd dadd x,xDBu+vx,xSuvx,xNDSBu =vx,xt+vx,xSu+vx,xb (227) where v (x, xi) is the value at x of the test functi on corresponding to node i and 12 21 ,1 ,2 ,2,1 20 0 0 0 10 10 1 00(1)/2j jj jjnn nn E N B D (228) 2forplanestressforplanestress ; forplanestrainforplanestrain (1)(1) E E E (229) PAGE 33 33 and 1 21if is prescribed 0 ; 0 0if is not prescribed on i i isuu S S S u S (230) Equation 228 can be simplified into the follo wing system of linear algebraic equations in ju 1 ,1,2,...,N ijji jiNKuf (231) where N is the total number of nodes and u is the generalized fictitious displacement vector. The socalled stiffness matrix K and the load vector f are defined by ()()()ssusuijvijijijdaddK x,xDB+vx,xSvx,xNDBS (232) ()()()stsusiiiidddfvx,xtvx,xSuvx,xb (233) As can be seen from Equation 232 the system stiffness matrix is banded but unsymmetric. The locations of the nonzero entire s in the stiffness matrix depend upon the nodes located inside the domain of influence of the nodes. Conclusions A brief review of the literature on meshless me thods has been given in this chapter. The meshless local PetrovGalerkin (MLPG) method (e xactly the MLPG1 method) is used to solve miromechanics problems of composite material s. The MLPG method including the moving least squares approximation, the choice of the weight function, the local symmetric weak form (LSWF), and the discretization of the weak form are presented. A computer code is written in MATLAB to solve problems in elastoelastics. Th is code is validated by solving the classical problem of Timoshenko beam. Results and discussion are given in Chapter 5. PAGE 34 34 Figure 21. Illustration of domain of influence. Figure 22. Distinction between andiiuu x Boundary node x1 x2 xi u iu iu ()hux x1 x3 x2 x4 x5 x Sampling point Domain of influence Domain of integration Node PAGE 35 35 Figure 23. Shape of weight function. Figure 24. Schematic representa tion of a body in MLPG method. t NodeI I te u I te 1 x 2 x ()iw xx PAGE 36 36 Figure 25. Schematics of MLPG. [Reprinted with permission from Atluri et al., 1999. (Page 353, Figure 5)]. PAGE 37 37 CHAPTER 3 MESHLESS LOCAL PETROVGALERKIN MICROMECHANICAL ANALYSIS Summary In this chapter a complete set of the meshless local PetrovGalerkin (MLPG) formulation is presented for the periodic composites containing mate rial discontinuities from the analysis of the unit cell. The MLPG formulation is also presente d for the generalized plane strain problem. The treatment of material discontinuity at the inte rface between the two phase s of the composite is presented by means of direct imposition of interf ace boundary conditions. The unit cell is discretized and periodic boundary conditions ar e set up and imposed between opposite faces of the unit cell. The treatment of periodic boundary conditions in the MLPG method is handled by using the multipoint constraint technique. The MLPG method is used in the micromechanical model for predicting the elastic constants/s tiffness properties of the composite. To our knowledge, this is the first study in which the MLPG method is formulated for the socalled MLPG method based micromechanical analysis. Examples are presented to illustrate the effectiveness of the current me thod, and it is validated by compar ing the results with available analytical and numerical solutions (Chapter 5). The current method has the potential for use in micromechanics, especially for textile composite s, where the meshing of the unit cell has been quite difficult. Introduction The traditional methods of computation such as finite element (FEM) or finite difference methods have successfully been applied to many problems in mechanics of composite materials. They are robust and thoroughly developed technique s, but they also have many drawbacks, and are impractical in a class of problems in compos ites. The reliance of the method on a mesh or a grid leads to complications. For example, in the modeling of large deformation processes, PAGE 38 38 considerable loss in accuracy arises when the elem ents in the mesh become severely distorted. In simulation of failure process, we need to m odel the propagation of cr acks with arbitrary and complex paths, and the simulation of phase tran sformation is difficult as reported by Belytschko at al (1996). The use of a mesh in micromechan ics creates difficulties in the treatment of discontinuities, which do not coinci de with the original mesh lines. In the development of advanced composite mate rials, especially textile composites, one of the major technical barriers in modeling te xtile composites such as braided and woven composites is the finite element mesh generati on. For composite materials with complex yarn architectures, actually, the meshing of individual yarns in the unitcell is quite simple; however the meshing of interfacial region between matrix phases and individual yarns is much more difficult. The region is multiply connected, and mesh in different phases may not be compatible (Kim and Swan, 2003). Recently, efforts have been devoted to develo ping methods for the pr ediction of effective material properties (Marrey and Sankar, 1995; Sankar and Marrey, 1997; Zhu et al, 1998; Karkkainen and Sankar 2006; Peng and Cao, 2000). The methods developed include the homogenization method, the finite element method, analytical mode l and experimental approach. Due to the immense variety of available compos ite materials and possible fabric construction geometry, it is impractical and very time c onsuming to characterize them by experimental approach. Analytical models, on the other hand, cannot deal with complex microstructures. Though the homogenization method is effective in pr edicting material proper ties (Yi at al, 1998; Takano at al, 1999) the huge computational cost li mits its application in simulating the forming of complex microstructures as textile composites. The finite element based micromechanical models are satisfactory for stiffness prediction because stiffness properties are based on volume PAGE 39 39 averaging of stresses and strains in the represen tative volume element of the composite. Further, the approximation involved in the FE meshing doe s not affect the result s significantly. However modeling the damage, especially progressive dama ge, requires accurate description of the stress field in different phases and re quires a very fine mesh as re ported by Marrey et al (1995, 1997). The FEM based micromechanical models have been successfully employed in predicting thermoelastic constants of ma ny composites materials; their us e for strength prediction under multiaxial loading conditions is not practical (Zhu et al, 1998). We expect that the tediousness and inaccura cies involved in mesh generation and hence inaccuracies in the results can be avoided by using a new class of meshless techniques that do not require a mesh to discretize the problem. On e such technique is the meshless local PetrovGalerkin method (MLPG) recently developed by Atluri et al (1998, 1999, 2000). The MLPG approach proposed by Atluri (1998) is one of th e several meshless schemes. The main advantage of this method compared to other meshless meth ods is that no background mesh is used to evaluate various integrals appearing in the local weak formulation of the problem. Therefore, this method is a truly meshless approach in terms of bo th interpolation of vari ables and integration of energy. The meshless methods have been demonstrat ed to be efficient in solving a variety of problems (Atluri, 2004; Raju and Chen, 2001; Raju and Phillips, 2003; Ching and Batra, 2001; Batra and Ching 2002). In this study the MLPG method is applie d to micromechanics of composites. A major drawback in applying meshless methods to problems of inhomogeneous materials is the treatment of material discontinuity at the interface. The highorder continuity of the moving least squares approximations (MLS), which is at least C1, allows for continuity of displacements and stresses throughout the subdomain. However, th e high order continuity imposes a difficulty PAGE 40 40 when considering discontinuitie s of derivatives at the inte rface of inhomogeneous bodies, because the shape functions from the MLS a pproximations do not have the Delta function properties. For the analysis of linear elastostatic problems by the meshless methodthe element free Galerkin method (EFG), Cordes and Moran ( 1996) used the method of Lagrange multipliers; Krongauz and Belytschko (1998) employed a special jump function at the line or the surface of discontinuity with parameters governing the st rength of discontinuity; and Cai and Zhu (2004) used the direct imposition of essential boundary and interface conditions. Whereas Cordes and Moran studied a twodimensional elastostatic problem, Krongauz and Belytschko as well as Cai and Zhu analyzed a onedimensional elastosta tic problem, all based on EFG method. Recently, Batra et al. (2004) also used the MLPG met hod to analyze heat conduction in which the continuity of the normal component of heat flux at the interface betw een two materials is satisfied either by the method of Langrange multipliers or by using a jump function. Li et al (2003) introduced a method for the treatment of material discontinuity by combining the MLPG 5 method and the MLPG 2 method. The MLPG 5 me thod is used inside the homogenous domain, the MLPG 2 method is used only on the interface The weak form of the MLPG 2 method uses the Dirac delta function on the nodal location. A nother main drawback in applying the MLPG method to the micromechanical analysis of comp osite materials is the treatment of periodic boundary conditions. Unlike in the finite elem ent method in which we can use a penalty formulation or Langrange multipliers to enforce the periodic boundary conditions, but it seems to be not able to use the penalty approach to en force the periodic boundary conditions in the MLPG method due to its nature of local weak forms over overlapping subdomains. PAGE 41 41 In the current study, a comple te set of the meshless loca l PetrovGalerkin (MLPG) formulation has been presented and including shear lo adings in directions of out of and inplane, thus permitting a more complete micromechanical analysis of a unidirectional composite material subjected to combined loading states. Fo r the treatment of material discontinuity at the interface between different phases of the composite in this study, we used the technique of direct imposition as shown by the authors (Da ng and Sankar., 2005, 2006), in which the inhomogeneous medium is considered separate ly as several homogeneous bodies, the MLS approximation is used separately within each of the homogeneous domains, the local domains of nodes that are close to but not ex actly on the material interface ar e adjusted not to cross over the interface boundary. The actual displa cements are to be computed at the nodes on the material interface, and after that, conditions of continuity of displacements at th e interface can also be directly enforced as in th e FE method. Especially for th e treatment of periodic boundary conditions, in this study, we initially use a tr ansformation method, and then use the multipoint constraint technique (Shepha rd, 1984; Mechnik, 1991) to handle the periodic boundary conditions. This chapter is organized as follows. Firstly, the unit cell analysis, constitutive relations and periodic boundary conditions are presente d. Secondly, a complete set of the MLPG formulation for modeling the unit cell of a unidi rectional fiber composite corresponding to six linearly independent deformations imposed on the unit cell. In this section we give a brief description of the MLS approxima tion, weak form and discretizati on for each case of the applied macrostrains. Especially, the treatment of materi al discontinuity and the treatment of periodic boundary conditions are also presented. Finally, conclusions are summarized. PAGE 42 42 Unit Cell Analysis for Three Dimensional Elastic Constants The macroscale properties of the composite ar e determined at a scale much larger than the dimensions of the unit cell, but comparable to the dimensions of the structural component. The average stresses at a point at the structural scale will be called the marcroscale stresses or macrostresses. The actual stresses at a point at the continuum level will be called the microscale stresses or microstresses. To distinguish the macroscale deformations and stresses from their microscale counterpart sa superscript M will be used to denote the macroscale deformation and stresses. The micromechanical analysis of a stru ctural composite is performed by actually analyzing the unit cell of the composite usi ng the MLPG method. We assume that uniform macrostresses exist th rough the composite. For a unidirectiona l composite, it is assumed that the fibers are circular in cross section packed in a square array. Thus the unit cell or the representative volume element is a square. The un it cell is shown in Figur e 31, the edges of the unit cell are assumed to be parallel to the coordinate axes x1, x2 and x3, with unit cells repeating in all three directions. The length of the unit cell in the xi direction is defined as 2 a The unit cell analysis assumes that the composite is under a unif orm state of strain at the macroscopic scale. However the actual stresses in the fiber and the matrix within the unit cell will have spatial variation. These stresses are ca lled microstresses. The macrostresses are average stresses M ij required to create a given state of macrodeform ations, and they can be computed from the microstresses ij obtained from the MLPG method as 1M ijij VdV V (31) PAGE 43 43 where V is volume of the unit cell, on the macroscale the composite is assumed to be homogeneous and orthotropic. The composite behavior is characterized by the following constitutive relation 111213141516 1111 2223242526 2222 33343536 3333 444546 2323 5556 3131 66 1212. orMM MM MM MM MM MM MMcccccc ccccc cccc ccc Symmcc c C (32) We can calculate all stiffness constants cij corresponding to the six cases of the macrostrains which are imposed on the boundaries of the unit cell. Therefore the elastic constants can be computed by the relation below MMS (33) where 31 21 123 32 12 123 1323 1 123 23 31 121 000 1 000 1 000 1 00000 1 00000 1 00000 EEE EEE EEE SC G G G (34) PAGE 44 44 In the micromechanical analysis, the unit ce ll is subjected to six linearly independent macroscopic deformations. In each deformation case one of the six macrostrains is assumed to be nonzero and the rest of the macrostrains are set equal to zero. The six cases are: case 1:111M; case 2:221M; case 3:331M ; case 4: 121M ; case 5: 231M ; case 6: 311M. The periodic boundary conditions for these six cases ar e shown in Table 31. Details of deriving the periodic boundary conditions can be found in (Marrey and Sankar, 1995). Except the case 4 (121M) which is not shown in Figure 31, for ot her cases, because of load and geometry symmetry of the unit cell, only one quadrant of th e unit cell need be considered to describe the behavior of the unit cell, and the entire unit ce ll (Adams and Crane, 1984) for these cases. Note that Case 3 of331Mis to be solved as a generalized plane strain problem, and the Cases 5 and 6 are longitudinal shear loading problems and also to be solved as special generalized plane strain problems (Adams and Doner, 1967; Gibson 1994).The procedure will be to solve the problem for specified uniform displacements which are derived from the six linearly independent macroscopic deformations. Meshless Local PetrovGalerkin Formulation of the Problem Moving Least Squares Approximation Scheme In the following we provide a brief descrip tion of the MLS approximation and also the MLPG formulation for the sake of completion a nd also to introduce the various notations and definitions (Atluri, 2004). In the ML PG method, the shape functions ()i x of the unknown trial function are found by the MLS approximation (L ancaster and Salkauskas, 1981). First, we consider a subdomain x, called the domain of definition of the MLS approximation for the trial function at point x, which is located in the problem domain The unknown trial approximation uh(x) of the function u (x) is defined by PAGE 45 45 1()()()m h jj juxpxaxTp(x)a(x) (35) Where 12()(),(),...,()T mppp pxxxx is a vector of the comp lete monomial basis of order m Examples of ()Tpxin the 2D problems are 12(1, )for linear basis, 3Tpxxm (36) 22 121122(1, ,,,)for linear basis, 6Tpxxxxxxm (37) The m unknown parameters ()jax can be determined by minimi zing the weighted discrete 2L norm, defined as 2 1 ()()()()n iii iJxw Txxpxaxu (38) Where n is the number of points in the neighborhood of x for which the weight functions w(x xi)>0, and iu refers to the nodal parameter of the function uh at the pointix. We choose the weight function to have the Gaussian distribution as 22 2exp(/)exp(/) 0 () 1exp(/) 0kk iiii ii k i ii iidcrc dr wxx rc dr (39) where iid xx is distance from the sampling point x to the node xi, and ri is the radius of the domain of influence of the weight function w(xxi) Finding the extremum of J(x) in Equation 38 with respect to a(x) leads to the following system of linear equation for the determination of a(x) ()() () AxaxBx (310) where 1()()()()n T iii iwAxxxpxpx (311) PAGE 46 46 1122()()(),()(),...,()()nnwww Bxxxpxxxpxxxpx (312) Solving a(x) from Equation 311 and substituting it into Equation 35, the MLS approximation can be defined as 1 ()()n h ii iuuxx (313) where the shape function()i xis defined by 1()()()m ijji jp1xxA(x)B(x) (314) ()i xis usually called the shape function of the MLS approxi mation corresponding to node i Note that()ij xdoes not satisfy the Kron ecker delta criterion()ijij x. Therefore, they are not interpolants, and the name approximation is used, i.e. ()h iiuxu (Figure 22 for a simple one dimensional case for the distinction between and iiuu ). For the matrix A to be invertible, the number of points, n must at least equal m ()nm In our study, we choose m = 3 and k = 1 in Equation 313, and take 4iirc (315) where ci is the distance from node i to its third nearest neighboring node. Meshless Local PetrovGalerkin Weak Formulation and Discretization Consider the general case of an i nhomogeneous body along with boundary conditions including periodic boundary conditions presenting to the cases of six macrostrains imposed. In this study, periodic boundary conditions are imposed using multipoint constraints. For simplicity, the problem has two phases se parated by a single interface int f in the domain (1)(2) PAGE 47 47 which is bounded by (1)(1)(1) uput(Figure 32). The interface is defined by(2) jn, the unit outward normal of(2) along the material boundary. The equilibrium equations in the two phases are (1)(1)(1) ,0inijjib (316) (2)(2)(2) ,0inijjib (317) where (1)(2)andijij are the Cauchy stress tensors, (1)(2)andiibb are the body forces in the two media. The boundary conditions are as follows (1)(1)(1)onijjitnt (318) (1)(1)oniiuuu (319) where it is the prescribed traction on a surface(1)t and iu is the prescribed displacement field(1) u and (1) jn is the unit outward normal to the boundary (1) u and (1) t are complementary subsets of The multipoint constraints for the faces pu & pu of the body are (1)on ; and 1...pupuiiipuim uuc (320) where m is number of nodes on the face pu (orpu ), pu and pu are complementary subsets of(1) pu, and ci is constant. Continuity of tractions a nd displacements on the interface int f are as (1)(2) inton iifuu (321) (1)(1)(1)(2)(2)(2) inton iijjijjiftnnt (322) PAGE 48 48 where superscripts 1 and 2 re fer to variables belonging to (1)(2)and respectively. In our study, we use the transformation approach and the multipoint constraint technique to handle the periodic boundary condi tions (Equation 320); therefore we dont need to enforce the periodic boundary conditions at variational level. The trea tment of the periodic boundary conditions using the multipoint constraint t echnique and the transformation approach is presented in detail at the end of this chapter. Also in our study, we consider the inhomogeneous medium as two separate homogeneous bodies, the weak form and its discretization are presented for each of homogeneous bodies. And then we apply interface displacem ent conditions (Equation 321) to reconnect the bodies via common nodes, and this work presented is similar to the collocation method in FEM. The treatment of material discontinuity at the interface int f using the transformation approach and then collocation method can be found later in this chapter. In addition, we use the penalty method to enforce the essential boundary conditions (Equation 319) for all cases of macrostrains imposed. The e ssential boundary conditions here may be fixed displacements at supports or may be marostra ins given through deformations imposed. As shown in the above section, the periodic boundary conditions are applied only when we use the full unit cell for the case of inplane shear121M to model the problem. Now, we derive the weak form and its disc retization for each of homogeneous mediums. For simplicity, we use general notations to presen t the weak form and its discretization below for a homogeneous medium. A generalized local weak form of Equations from 316 to 319 over a local subdomain s can be written as follows ,()()0suijjiiiiibvduuvd (323) PAGE 49 49 where ui and vi are the trial and the test functions, respectively, and u is the part of the boundary s of s over which essential boundary conditions are specified. In general, s ssL with s being the part of th e local boundary located on the global boundary and Ls being the other part of the local boundary ove r which no boundary conditions are specified, i.e., ss with s ssL In Equation 323, is a penalty parameter ( >>Youngs modulus/Length), which is used to impose the essential boundary conditions. Our experience has shown that the value of691010 gives good results. Also, the test functions vi are chosen such that they vanish on Ls, and this can be accomplished by using the weight function wi in the MLS approximation as also the test function vi, but the radius ri of the support of the weight function is replaced by the radius ro of the local domain s. Using integration by parts and the divergen ce theorem in Equation 323, after some algebraic operations, finally yields the expression in the matrix form as ssusustsusdadddaddv +vuvt=vt+vu+vb (324) In Equation 324, v denotes the strain matrix derive d from the test functions, and is the stress vector derived from the trial functions. That is 11 (1)(1)(1) 112212 22 (2)(2)(2) 112212 12,v (325) where the superscript i denotes the ith test function. Functions v, u, t, and b are defined as follows 1112111 2122222,,, vvutb vvutb vutb (326) PAGE 50 50 The two sets of test functions v in Equation 326 should be linearly independent. The simplest choice for v as proposed by Atluri and Zhu (2000) is or ijijvvv vI where ij is the Kronecker delta and I is the identity matrix. As long as the union of all lo cal subdomains covers the gl obal domain, Equations from 316 to 319 will be satisfied in the global domain and on its boundary respectively. Substituting the MLS approximation Equation 313 into Equation 324, and summing over all nodes leads to the following discreti zed system of linear equations 111 (,)(,)(,) (,)(,)(,)ssusu stsusnnn vijjijjijj jjj iiidadd dadd xxDBuvxxSuvxxNDSBu =vxxt+vxxSu+vxxb (327) where v ( x ,xi) is the value at x of the test functi on corresponding to node i and 12 210 0 nn nn N (328) ,1 ,2 ,2,10 0j jj jj B (329) 210 10 1 00(1)/2E D (330) 2forplanestressforplanestress ; forplanestrainforplanestrain (1)(1)E E E (331) and 1 2 u1if is prescribed 0 ; 0 0if is not prescribed on i i iu S SS S u (332) PAGE 51 51 Equation 327 can be simplified into the followi ng system of linear algebraic equations in ju. 1 1,2,...,N ijjj jKufiN (333) Ku=f (334) where N is the total number of nodes and u is the generalized fictitious displacement vector. The socalled stiffness matrix K and the load vector f are defined by (,)(,)(,)ssusuijvijijijKdadd xxDBvxxSvxxNDBS (335) (,)(,)(,stsusiiii f dddvxxtvxxSuvxx)b (336) As can be seen from Equation 335 the system stiffness matrix is banded but unsymmetric. And the subscript i stands for the test function of the ith node, the subscript j denotes stands for the trial function of the jth node. The entry Kij in the global stiffness matrix is generalized by integrating over the inte rsections of the subdomain of node j on which the trial function is nonzero, with the subdomain centered at node i on which the test function is nonzero. Integration over subdomains in the MLPG method at the same area can happen more than one time. Unlike in the finite element method, overlapping subdomains in the MLPG method does not increase the global stiffness of the st ructure. This can be explained as follows We view the MLPG method as a weighted residual method and FEM as a minimum potential energy method. In the weighted residua l method we require the residual be equal to zero in some average sense over the domain. In that sense it does not matter we overlap the integrals. In fact that emphasi zes certain overlapping regions more than the other. There is PAGE 52 52 nothing wrong in it. In the FEM we integrate over an element to calculate the strain energy, which contributes to K If we overlap, then we will be over estimating the stiffness by a large factor. If we treat both the MLPG method and FEM as a weighted residual method and if we do not use integration by parts then there is no pr oblem with overlapping. However, when we use integration by parts to reduce th e continuity requirements, the bounda ry terms come into picture. In the MLPG method the boundary terms vanish, as the test function is equal to zero along the local domain boundary. However in the FEM the boundary terms over an element boundary exists, but they cancel when you consider adjacent element. T hus all the boundary terms cancel except the global domain boundary terms, which is then accurate. When one overlaps or adds on elements one over the other, the boundary terms do not cancel out a nd contribute to the residual. InPlane Compression, Tension and Shear Cases For the cases of compression/tension/shear, the macrostrains111M ,221M and 121M are given respectively, we can obtai n the fictitious displacements from Equation 334 to Equation 336, and the actual displacements from Equa tion 313. After that, the microstrains and microstresses can be computed. Finally the macr ostresses and elastic co nstants can be obtained from Equation 31 and Equation 32. The treatmen t of the periodic bounda ry conditions for the cases 121M is required and will be described in this chapter. Generalized Plane Strain Problem Consider a unidirectional fiber reinforced com posite in Figure 31. Introduce the coordinate system as follows. The x3axis is aligned with the direction of fibers (assuming fibers are all straight and laid parallel to each other), and the x1and x2axes are lie in a crosssection of the material perpendicular to the fibers. The material is homogeneous in the x3 direction. When PAGE 53 53 investigating material properties, it is reasonable to assume that composite is subjected to a uniform deformation in the x3direction, i.e. all the strains are independent of x3. The uniformity of the deformation and the material homogeneity in the x3direction enables one to simplify the problem of the microscopic deformation of th e material from a general threedimensional problem in the x1x2x3 space to a special twodimensional one in the x1x2 plane, a generalized strain problem as shown by Adams and Crane (1984) in which all the strains are functions of x1 and x2 only. It implies that the stresses are not functions of x3 from the constitutive relations for linearly elastic and homogeneous materials. Th en the displacements can be written as: 12 12 123(,) (,) (,)ouuxx vvxx wwxxx (337) where o is the constant macroscopic direct strain in the z direction. 33o (338) Expressions for the strains can now be simplified according to Equation 337. 1112 121 2213 213; ; uuv x xx vwu x xx 3323 323; wwy x xx (339) where represents normal strain and represents engineering shear strain. The constitutive equations in terms of stressst rain relations are given as PAGE 54 54 11 11 22 22 33 33 12 1210 10 10 (1)(12) 12 000 2 E (340) The relations can be expressed in a reduced form as 11 22 12(1)(12) 10 10 (1)(12)(1)(12) 00(12)/2 0oE EE (341) or oDF (342) and 33(1) 0 (1)(12)(1)(12)(1)(12)oEEE (343) or 33(1) (1)(12)T oE F (344) where the matrix D is defined as shown in Equation 330 for the conventional plane strain problem and F* (1)(12) (1)(12) 0 E E *F (345) PAGE 55 55 Reactions t at u can be expressed as oNNDNF t (346) where the matrix N is defined in Equation 328. Since 33 is constant and prescribed we can choose the components 330v in the strain matrix in Equation 325 which comes from the test function. Therefore stress 33 will not contribute to Equation 324, when Equation 342 and Equation 346 are substituted into the first and third term, respectively, of the right hand side of Equation 324. After some algebraic operations, finally two terms are added to the lo ad vector in Equation 336 as shown below, while the stiffness matrix K in Equation 335 remains unchanged. (,)(,)(,) (,)(,)stsus ssuiiii oo vifddd dd ** ivxxtvxxSuvxxb xxFvxxNSF (347) Out of Plane Shear Test: Special Generalized Plane Strain Problem The case of outofplane shear test is a l ongitudinal shear problem (Adams and Doner, 1967; Gibson, 1994). The macroscopic deformations 13231or1MM is applied to the unit cell, and for reasons mentioned above, only the quadrant of the unit cell need be considered as shown in Figure 33 with the shear loading. The problem of longitudinal shear loading is defined by a displacement filed of the form 120,(,) uvwwxx (348) For such a system of displacement, the only nonvanishing stress components are 31133223 12, ww GG x x (349) PAGE 56 56 where G13 and G23 are the shear moduli of the constituent materials. Thus we implicitly assume that the principal material axes of the fiber a nd matrix are parallel to the coordinate axes. The assumed displacement field automatically satisfies the equilibrium equations in the x1 and x2 directions, while equilibrium in the x3 direction requires that 3132 12 22 22 120 or 0 xx ww G xx (350) Using procedures similar to that given in the above section of the weak formulation and discretization, we obtain the weak form below 3132 3 12()0subvdwwvd xx (351) where w and b3 are the prescribed displacement field and the body force, respectively in the x3 direction, and the di scretized equations 1 ;1,2,...,N ijji jKwfiN (352) where N is the total number of nodes and w is the fictitious displacement in the x3 direction, the stiffne ss and load vector are (,)(,)(,)ssusuijvijijijKdvdvd=xxGB+xx xxNGB (353) 33(,)(,)(,)stsusiiii f vtdvwdvbdxxxxxx (354) where (,)vi xx is the value of the test function corresponding to node i, evaluated at the point x. 3t is defined as the prescrib ed traction on a surface in the x3 direction, and PAGE 57 57 ,1 13 132312 ,2 230 ,,,and 0j vj jG nn G NBG (355) Treatment of Material Discontinuity In general, the MLS approximation in Equation 313 does not pass through the nodal data, which are approximate values at nodes. This leads to some difficulties in imposing essential boundary conditions and treating material discont inuity. As mentioned above, a penalty method is used to enforce the essen tial boundary conditions in the current study. Treatment of material discontinuity is as follows. The current method involves considering th e inhomogeneous medium as two separate homogeneous bodies and then applying interf ace conditions to rec onnect the bodies. For example let us consider a twophase problem as shown in Figure 34. The inhomogeneous body is separated into two homogeneous parts. In each part, we have th e weak form of the problem as described above. We will determine the actual di splacements at the interface, and after that, conditions of continuity of displacements at the in terface can also be direc tly enforced as in the FE method. In this method, common nodes are used at the interface that belongs to both materials. In Figure 34, nodes (1012) are define d at the interface of materials 1 and 2. We use the following example to illustrate the proposed method. For material 1: Assume that neighbors of point 10 include 7, 8, 10 and 11; neighbors of point 11 include 7, 8, 9, 10, 11 and 12; ne ighbors of 12 include 8, 9, 11 and 12. Consider the following MLS approximation 1 ()()n ii iuxxu (356) Such that, at points xj, we have the actual displacements as PAGE 58 58 1 ()n jiji iuxu (357) For nodes at the interface, we have 1071078108101010111011 1171178118101110111111121112 1281289129111211121212 ()()()() ()()()()() ()()()()uxuxuxuxu uxuxuxuxuxu uxuxuxuxu (358) In matrix form 7 8 1071081010101110 9 11711811911101111111211 10 1281291211121212 11 12 ()()0()()0 ()()()()()() 0()()0()() u u uxxxx u uxxxxxx u uxxxx u u (359) The variables can be partitioned into two subsets as follows r brbb bu uu u (360) where b = 10, 11, 12 denote the nodes on the interface, and r = 7, 8, 9 denote the nodes which do not lie on the interface. Expanding the matrix bto explicitly show the x and y directions, we have 12 12..................... ()0()0...()0 0()0()...0() .....................bibiNbi b bibiNbixxx x xx (361) where bi x bi = b 1 b 2,, bnb are the coordinates of nodes at interface. bis a 2nbN matrix, r is a 2nbnr matrix and b is a 2nbnb matrix; ru is a 2nr 1 vector, and bu is a 2nb 1 vector, all are known matrices. PAGE 59 59 Furthermore nb: number of nodes on interface, nr: number of nodes which do not lie on the interface, nb + nr = N : total number of nodes in body. According to Equation 360, we can compute the fictitious displacements at interface as: 1bbrrbuuu (362) We can rewrite Equation 334 in the form as r rb bu KKf u (363) where bKis a known 2Nnb matrix, and rKis a known matrix of size 2Nnr Substituting Equation 362 into Equation 363 an d after some algebraic operations we obtain the following system of linear al gebraic equations for material 1. (1)(1)(1)Kuf (364) 11 (1)(1)(1)(1)(1)(1)(1) (1) (1) (1) (1) (1) (1)rbbrbb r b r bKKKK u u u f f f (365) bu is the actual displa cement at interface. We use the same procedure to obtain the fo llowing system of linea r algebraic equations for material 2 (2)(2)(2)Kuf (366) PAGE 60 60 11 (2)(2)(2)(2)(2)(2)(2) (2) (2) (2) (2) (2) (2)bbrbbr b r b rKKKK u u u f f f (367) The conditions of continuity of the interface displacements can be directly enforced to reconnect the two separate homogeneous bodies through the assembly of the global stiffness matrix of material 1 and 2 to obtain the global stiffness matrix of the composite (inhomogeneous body) through common nodes at the interface. ()()() (1) () (2) (1) () (2) 0ccc r c b r r c rKuf u uu u f f f (368) where the superscripts 1, 2 and c stand for material 1, 2 a nd composite, respectively. Therefore, the method presented here is si milar to the collocation method in the FE method. It means that we can use co mmon nodes at the interface. Note that bis a very sparse matrix, and one should choose as many points as possible on the interface. In fact, we can use this technique to enfor ce the essential boundary conditions, instead of using the penalty method. However, we use the tec hnique only to treat the material discontinuity and not for essential B.Cs in order to avoid the computational burden of d ealing with large size matrices and also to avoid reshaping the global stiffness matrix many times. PAGE 61 61 Treatment of Periodic Boundary Conditions In the current study, the treatme nt of periodic boundary conditi ons is required for the case of shear test in the x1x2 plane corresponding to th e imposed macrostrain 121M on the unit cell. The quadrant of the unit ce ll in Figure 31B is not suitable for this analysis. Therefore we use the full model of the unit cell as show n in Figure 35 to analyze the problem. The periodic boundary conditions for the unit cell in this case (121M ) are given in Table 31. We need to only treat the periodic boundary conditions on displacements. We rewrite the periodic boundary condition for this case as 0 2 0 0RL RL MN MNuu vva uu vv (369) Below is an algorithm in general for the treatment of periodic boundary conditions on displacements in the MLPG method using the mu ltipoint constraint technique (Shephard, 1984; Mechnik, 1991). Assume the periodic boundary conditions for the opposite faces L & R of the unit cell (Figure 35) is: 1...RL iiiim uuc (370) where m is number of nodes on the face R (or L ), and ci are constants. The periodic boundary conditions that relate to the complete degree of freedom of the unit cell u can be written in the form A uC (371) PAGE 62 62 where the matrices A and C contain constants. There are more degrees of freedom in u than the constraint equations, so the matrix C has more columns than rows. u is partitioned as: where,andrLrrL rLRrL RLRuuu AACuu uuu = (372) where A is a known 2mN matrix, u is an unknown 2N vector, N is total number of nodes in the unit cell. ruis a (2 N 4 m ) 1 vector of fictitious displace ments of nodes which do not lie on the boundary R & L (superscript r stands for remained displacements). Luis a 2 m 1 vector of actual displacements of nodes on the boundary L Ruis a 2 m 1 vector of actual displacements of nodes on the boundary R rLuis a (2 N 2 m ) 1 vector. ArL is a known 2 m (2 N2m ) matrix, and AR is a 2 m m matrix. Actually, columns of the matrix ArL that are corresponding to the added fictitious displacement vector ru in Equation 372 are zero. Because there are as many degrees of freedom Ruas there are indepe ndent equations of the constraint equations, the matrix AR is square and nonsingular. From Equation 372, we obtain RRrLrLuACAu (373) The complete array of degrees of freedom can be written as 0rL rL R RrL Ru I u AC AA u (374) or rL ouTuQ (375) where I is a (2 N 2 m ) (2 N 2 m ) unit matrix, and the matrix T is a 2 N (2 N 2 m ) matrix, while Qo is a 2 N 1 vector. PAGE 63 63 The global stiffness matrix can be transformed into the form below by using the previously developed method to treat the material discontinuity at the interface with the displacement vector partitioned as shown in Equation 372. KuF (376) Premultiplying the global equations by TT and substituting for u from Equation 375, we obtain the reduced equation set as rLrLrLKuF (377) where rLT rLT oKTKT FTFKQ (378) After urL is computed from Equation 377 and Equation 373 yields Ru. Conclusions Meshless local PetrovGalerkin micromechanical analysis has been presented for modeling the unit cell of the composites. A complete set of the MLPG formulation is formulated for six cases of the applied macrostrains, for the fi rst time, the treatmen t of periodic boundary conditions using multipoint constraints, and the tr eatment of the material discontinuity using the technique of direct imposition of interface boundar y conditions in the MLPG method have been presented, the generalized stra in plane problem has been form ulated to model the unit cell by using the MLPG method for the first time. The current method shows prom ise in application to the micromechanics of textile composites where the complexities and inaccuracies involved in the FE mesh can be avoided. Numerical exampl es, results and discussi on to illustrate the effectiveness of the current method ar e presented in the next chapter. PAGE 64 64 Table 31. Periodic boundary c onditions for the MLPG method. Case Constraints between left and right faces Constraints between top and bottom faces Out of plane strains 111 u1(2a,x2) u1(0,x2) = 2a u2(2a,x2) u2(0,x2) = 0 ui(x1,2a) ui(x1,0) = 0 i = 1,2 330 ,310 ,230 221 ui(2a,x2) ui(0,x2) = 0 i = 1,2 u1(x1,2a) u1(x1,0) = 0 u2(x1,2a) u2(x1,0) = 2a 330 ,310 ,230 331 ui(2a,x2) ui(0,x2) = 0 i = 1,2 ui(x1,2a) ui(x1,0) = 0 i = 1,2 331 ,310 ,230 121 u1(2a,x2) u1(0,x2) = 0 u2(2a,x2) u2(0,x2) = 2a ui(x1,2a) ui(x1,0) = 0 i = 1,2 330 ,310 ,230 231 u3(2a,x2) u3(0,x2) = 0 u3(x1,2a) u3(x1,0) = 2a 330 ,310 311 u3(2a,x2) u3(0,x2) = 2a u3(x1,2a) u3(x1,0) = 0 330 ,230 Figure 31. Unit cell of a unidirect ional composite. A) Full unit cell. B) Quadrant to be analyzed due to symmetry. X2 X2 2a 2a a a X1 X3 X3 X1 1 13 23 2 3 1 2 13 23 3 A B Fiber Matrix PAGE 65 65 Figure 32. Illustrative example of an inhom ogeneous body including multipoint constraints. Figure 33. Longitudinal shear problem. X1 Matrix Fiber X3 13 X2 u pu pu t n ( 1 ) n ( 1 ) n ( 2 ) x2 x1 x3 (2) intf (2) (1) (1)(1)(1) (1)(1)(1)(1) intuput uputf (1)(2) PAGE 66 66 Figure 34. Illustration of an inhomogeneous body. Figure 35. Unit cell and periodic boundary conditions for the deformation121M. N M L R x1 x2O2a 2a PAGE 67 67 CHAPTER 4 STIFFNESS AND SHEAR CORRECTION FACTOR PREDICTION OF A CLASS OF COMPOSITE BEAMS Summary The current theories to determine shear correc tion factors have been the subject of much previous research. However, th e computation of shear correcti on factors for textile composite beams is still a challenging problem. In this chapter we suggest a mi cromechanical model to predict stiffness properties of a class of compos ite beams with specific emphasis on prediction of shear correction factors. Three linearly indepe ndent deformations, namely, pure extension, pure bending and shear are applied to the unit cell. The top and bottom surfaces of the beam are assumed to be traction free, periodic boundary conditions are set up and are enforced on the lateral boundaries of the unit ce ll. The suggested model based on both the finite element method and the meshless local PetrovGalerkin (MLPG) method is straightforward and is a truly micromechanical model. In addition, an exact so lution for the shear correction factor of general laminated beams is also derived in this chap ter. Examples are presented to illustrate the effectiveness of the present micromechanical mo del. It is first validated for isotropic and laminated beams by comparing the nu merical results with available analytical solutions (Chapter 5). Prediction of the flexural s tiffness properties and the shear co rrection factor of the isotropic and laminated beams using both the finite el ement method and the meshless local PetrovGalerkin (MLPG) method are in very good agreemen t with the analytical solutions. Finally, the model is illustrated for a thin walled textile composite beam. Background The shear correction factor k was introduced by Timoshenko (1921) to compensate for the error resulting from the assumption of a constant shear strain in the first order shear deformation beam theory (FSDT), as opposed to the classi cal parabolic distributi on. The aforementioned PAGE 68 68 commercial programs usually ignore shear deforma tions due to the fact that they are very difficult to compute shear co rrection factors (Sapountzakis and Mokos, 2006). Though these deformations are quite small in many engineerin g applications, they may be dominant in some situations, where bending moment s are small compared to shear forces acting on the member. This is normally true in short span be ams, curved box girder bridges, etc. The shear correction factors of beams, plat es, and shells have been studied by some researchers during past decades. Cowper (1966) derived expressi ons for shear coefficient of beams with different crosssecti on. Using Saint Venant's flexur e problem, the shear coefficient was derived by assuming transverse shear stress distributions due to a linearly varying shear force. Chow (1971) presented a theory on the pr opagation of flexural waves in an orthotropic symmetric laminated plate. The dynamic equations of orthotropic laminate d plates were derived using Timoshenko beam theory to include the effect s of transverse shear a nd rotatory inertia, and a one dimensional procedure was adopted to obta in the shear correction factor by equating the shear strain energy from Timoshenko beam theo ry to the shear strain energy obtained from equilibrium conditions. Whitney (1973) extended the procedure of Chow (1971) to derive an expression for the shear correction factor of a cro ssply laminated plate, a nd the accuracy of this method was illustrated by comparing the transverse deflections of a rectangular laminated plate obtained from FSDT against elasticity solutions for symmetric and unsymmetric layers. Also, Whitney (1972) developed a procedure for accurate ly calculating the mechanical behavior of thick laminated composite and sandwich plates including the effects of transverse shear deformation in the analysis. In this work, an expression for the shear correction factor was derived by considering cylindrical bending about the length and breadth of the plate, and discrete values of shear correction factors were presente d for symmetric angleply, antisymmetric angle PAGE 69 69 ply, and symmetric crossply la minated plates and sandwich pl ates. Whitney (1972) discussed the variation of shear correction factor due to vari ation in number of layers of the laminate, and it was shown that the shear correction factor does not approach th e classical value for homogeneous plates as the number of layers is in creased. A simplified analysis of static shear correction factor for composite beams was presented by Bert (1973), and an expression was derived for the shear correction f actor of crossply laminated beams. The results obtained from this analysis for the shear corr ection factor of beams were comp ared with the results presented by Whitney (1973) for plate geometries, and the nu merical values matched up to four significant figures. From this work, it can be observed that the shear correction f actors for symmetric and asymmetric crossply laminates are identical fo r beam and plate geometries. Dharmarajan and McCutchen (1973) extended Cowper's formulatio n (1966) for orthotropi c beams and presented expressions of shear coefficient for rectangular and tubular crosssections. A predictorcorrector approach to obtain shear corre ction factors was proposed by Noor and Peters (1989), who obtained the shear correction fact ors for laminated cylinders. In the predictor phase, FSDT was used to predict the gross respons e characteristics of th e structure by selecti ng initial values for the shear correction factor. The corrector phase involved the computation of transverse stresses from the predicted solution using three dimens ional equilibrium equatio ns. Subsequently, the shear correction factors were computed by matching the integral of the tr ansverse shear strain energy in the thickness direction obtained from 3dimensional el asticity with that of FSDT. Recently, the thinwalled composite beams have been studied by Barbero et al (1993) and Kim et al (1996). Hutchinson (2001) derived a new formula for k of dynamic beams, and Sapountzakis and Mokos (2006) used a pure BEM approach to calculate k for 3D dynamic beams. PAGE 70 70 The first effort made for the thinwalled text ile composite beams actually is the previous work of the authors (Sankar a nd Marrey, 1993). In the work, Sanka r et al suggested a two unitcell method to predict the shear stiffness by equating the shear st rain energy obtained from the constitutive relations of FSDT to that obtained using the actu al shear stress distribution calculated. From the computed shear stiffness on a cross section, we can easily calculate the shear correction factor k. Most of the previous work has placed empha sis on deriving shear corr ection factors of the some specific structures such as homogeneous an d anisotropic beams, plates, shells, sandwich, and laminated beams. Derivation of k by most of the authors was based on the comparison of the shear strain energy for the actua l beam to the counterpart for an equivalent Timoshenko beam. Most of the existing expressions for the shear correction factor are qu ite complex involving integral expressions and are not explicit enough to be used in engineering practice. In the current study, we suggest a truly mi cromechanical model using both FEM and the meshless local PetrovGalerkin method (Dang and Sankar, 2007) to predict k for a class of structural composites with special emphasis on the thinwalled textile composite beams, the suggested model is straightforward and is able to capture the effect of more complicated fiber architectures. The model is first verified by applying it to isotropic and laminated beams for which the results are known, and then illustrated fo r the thinwalled textile composite beams. In addition, a method to derive the shear correction factor of general laminated beams is presented in this chapter. This chapter is organized as follows. Firstly, we derive an exact expression of the shear correction factor for laminated beams. Secondly presents the micromechanical analysis of composite beams including a description of the unit cell, constitutive relations, and periodic PAGE 71 71 boundary conditions of the unit cell is set up for each case of applied loadings. Finally, conclusions are summarized. Exact Solution for the Shear Correction Factor of Laminated Beams Shear Correction Factor of Isotropic Beams Consider a homogeneous beam with r ectangular cross section, with width b and height h. The actual shear stress distribution through the thickness of the beam, from a course on mechanics of materials is given by 232 1, 222c xzQzhh z Lhh (41) where Q is the transverse shear force. The transverse shear stress in the firstorder theory is a constant as f xzQ Lh (42) The strain energies due to transverse shear stresses in the two theories are 2 2 1313 2 2 131313 25 1 22cc szx A ff szx AQ UdA GGLh Q UdA GGLh (43) The shear correction factor is the ratio ofto f c s sUU, which gives k=5/6. This value is considered as the standard value and is wide ly used for homogeneous beams. Shear Correction Factor of Laminated Beams In this section we reintroduce an exact ex pression of the shear correction factor for general laminated beams, which has been deri ved by Sankar (1994). Consider a laminated beam in the xzplane. We assume a state of plain strain parallel to the xzplane. The case of plane PAGE 72 72 stress can be handled in an analogous manner. The layers of the beam are assumed to be orthotropic and satisfy the constitu tive relations of the following form 1113 1333 x xxx zzzzcc cc (44) 55 x zxzc (45) The strainstress will be expressed by 1113 1333 x xxx zzzzss ss (46) The coordinate system is shown in Figure 41. Unlike in the traditi onal beam theories, the xaxis is assumed to lie at the bottom surface of the beam. We introduce two integral operators Iz and Ih which are defined as follows 0()()z z I fzfd (47) 0()()h h I fzfd (48) The choice of location of the x axis and the above integral operators make it very convenient to write the many integrals in the current work in a compact form. The displacement field in the beam can be assumed to be of the form 01(,)()(,) wxzwxwxz (49) 0 01(,)()(,) w uxzuxzuxz x (410) where w0 and u0 are displacements of points on the x axis, and w1( x z ) can be considered as corrections to the classical beam theory approx imation. In the following we will show that the functions u1 and w1 can be expressed in the terms of the stress x x PAGE 73 73 We start with strain displacement relations zzw z (411) xzuw zx (412) The above relations can be integrated to obt ain formal expressions for displacement field in the beam as follows 0(,)(,)zzzwxzwxzI (413) 0 0(,)(,)zz zzzzzw uxzuxzzIII x x (414) When multiple operators are used as in E quation 414, we will assume that they are evaluated from right to left. This will avoid the use of parenthe ses to separate the integral operators. We will use the stressstrain relations Equation 45 and Equation 46 to replace zz and x z in Equation 413 and Equation 414 by terms containing xz,andxxzz to obtain the following expressions for displacement field. 01333(,)(,)()zxxzzwxzwxzIss (415) 0 01333 55(,)(,)xzxx zz zzzw uxzuxzzIIIss x cxx (416) The stress equilibrium equations (neglecting body forces) are 0xxxzxz (417) 0xz zzxz (418) By integrating the above equations andzxzz can be expressed in terms of x x as follows PAGE 74 74 x x zxzI x (419) 2 2 x x zzzzII x (420) In deriving Equation 419 and Equation 420 we have assumed that no tractions act on the bottom surface of the beam, z =0. Thus the external forces can only be applied to the top surface of the beam, z= h The case of a beam loaded at the bottom can be treated in an analogous manner, or by simple shifting the coordinate system. The case where loads act on both top and bottom surfaces can be dealt by superposition. Substituting from Equation 419 and Equation 420 into Equation 415 and Equation 416, the expressions for displacements take the form. 2 01333 2(,)() x x zxxzzzwxzwxIsIsII x (421) 3 0 01333 3 551 (,)() x xxxxx zzzzzzzzw uxzuxzIIIIsIIsII x cxxx (422) By comparing Equation 421 and Equation 422 with Equation 49 and Equation 410 respectively, one can identify the corrections terms11(,)and(,) uxzwxz Equation 421 and Equation 422 can be thought of as a set of assumed displacement fields in terms of 00()and() uxwx and (,)xx x z where x x is an arbitrary function continuous in x and piecewise continuous in z Now we will use the standard proce dures employed in beam theories to determine these functions. Next we derive expressions for normal stra ins from Equation 421 and Equation 422. By differentiating both sides of Equation 421 and Equation 422 with respect to z and x respectively, we obtain PAGE 75 75 2 131 x x zzxxsL x (423) 224 00 23 24 x xxx xxuw zLL xxxx (424) where the operators L1, L2, and L3 are defined as follows 133(.)(.)zzLsII (425) 213 551 (.)(.)zzzzLIIIIs c (426) 333(.)(.)zzzzLIIsII (427) From Equation 44 the normal stress x x can be written as 1113 x xxxzzcc (428) Substituting for the normal strains from Equation 423 and Equation 424 into Equation 428 we obtain 2242 00 11231213131 242 x xxxxx xx xxuw czLLcscL x xxxx (429) Bringing the x x terms to one side of the above equation, we obtain 2242 00 1123131 242 x xxxxx xxuw czLLcL x xxxx (430) where 11 11 13131 c c cs (431) 13 13 13131 c c cs (432) The axial force and bending mome nt resultants are defined as PAGE 76 76 x hxxNI (433) x hxxMIz (434) I must be noted that the bending moment Mx is referred to the bottom surface of the beam. If the axial force resultant is zero, then this bending moment will be the same as that referred to the middle surface of the beam. Submitting for x x from Equation 430 into Equation 433 and Equation 434, we obtain 224 00 213 224 x xxx xuw NABMMM x xxx (435) 224 00 213 224 x xxx xuw MBDKKK x xxx (436) where the operators Mi and Ki are defined as follows 1131(.)(.)h M IcL (437) 2,3112,3(.)(.)hMIcL (438) 1131(.)(.)hKIzcL (439) 2,3112,3(.)(.)hKIzcL (440) The beam axial, coupling and flexural stiffness terms A B and D respectively are given by 2 11 0,,1,,hABDczzdz (441) The above definition of the laminate stiffn ess coefficients is different from that of traditional composite literature where the midplan e of the beam is considered as the reference plane. Thus the coupling coefficient B does not vanish for beams wh ich are symmetric about the midplane. Equation 435 and Equation 436 can be writte n conveniently as PAGE 77 77 0 2 0 2 x xc x xcu NN AB x MM BD w x (442) The terms that Nxc and Mxc denote can be obtained by comparing Equation 442 with Equation 435 and Equation 436. 24 213 24 x xxx xcNMMM x x (443) 24 213 24 x xxx xcMKKK x x (444) Equation 442 can be solved for the reference plane strain and curvature to obtain 0 2 0 2 xxc x xcu NN AB x MM BD w x (445) where 1AB AB B D BD (446) The exact solution for a laminated beam can be compared with the Timoshenko beam equation to obtain an expression for shear correction factor for laminated beams. The procedure described here is similar to that in Bolley and Tolins (1956). The governing equation for shear deformable laminated beam is (Whitney, 1987): 22 0 22 551 0x xwM DM xkAx (447) where 5555 zAIc (448) PAGE 78 78 and k is the shear correction factor. It is assu med that there is no net axial force acting on the bean, i.e., Nx=0. We will derive a similar equation us ing the present method. From Equation 445 we obtain 2 0 20xxcxcw DMBNDM x (449) The expression for Nxc and Mxc can be found in Equation 443 and Equation 444. We will assume that iv x M and its higher derivatives are equal to zero, and use the beam theory expression for x x i.e., 11 x xxcBzDM in Equation 443 and Equati on 444 and substitute into Equation 449 to obtain 2 2 0 212111 2 20x xwM DMBMMDKKcBzD xx (450) Comparing Equation 450 and Equation 447, we can readily obtain an expression for the shear correction factor for a laminated beam as 1 55212111kABMMDKKcBzD (451) It should be remembered that K1, K2, M1 and M2 are all integral operators acting on 11cBzD For a homogeneous beam the integrati ons are much simpler, and we obtain 1 133355 5513 111.21csc kcs c (452) For an isotropic beam under pl ane stress parallel to the xz plane the above expression reduces to k = 0.8333(1+ ), and for the case of plain strain we obtain k = 0.8333/(1). If the Poissons ratio is neglected ( c13 = s13 = 0), then we obtain the standard result k = 0.8333. PAGE 79 79 Micromechanical Analysis of Thin Walled Textile Composite Beams Geometric Model of Textile Composite Unit Cell The thin walled textile composite beam is assumed to be in the xzplane with unit cells repeating in the x direction. For an example, the textile composite beam with 9 unit cells is prescribed in Figure 42B. A state of plane strain parallel to the xz plane is assumed. The idea behind using a beam model to illustrate the propos ed method and its scope are similar to the one dimensional fabric strip models presented by Ishikawa and Chou (1983), Yoshino and Ohtsuka (1982). In our study, a unit cell as shown in Figur e 43 is adopted for the thin walled textile composite beams, the left end of the unit cell x = 0 has minimum support constraints to prevent rigid body translation and rota tion. Geometric parameters are given in Table 41. The elliptic curves of transverse yarns, z1, z2, and z3 are defined, respectively, as 22 11 22 12/2 1 xaze ee (453) 2 2 21 22 121 0 ze x ee x (454) 22 31 22 121 xaze ee xa (455) The lines of longitudinal yarns, z4 and z5 are defined, respectively, as 42 cos1 22 hxd z a (456) 52 cos1 22 hxd z a (457) PAGE 80 80 Because of the symmetry of unit cell about the axis at the middle of the unit cell, we can use a half unit cell with boundary conditions as shown in Figure 44 to model the problem. Based on our study results, this unit cell model on the shear problem can capture material microstructure and give a good pr ediction. Therefore it is used to predict the shear stiffness coefficient ( K33 or k A55) and the shear correction factor k. Constitutive Relation of Te xtile Composite Unit Cell For a thin walled textile composite beam, the whole composite structure can be represented by repeating the unit cell as shown in Figure 42. A state of plane strain parallel to the xzplane is assumed. On the macroscale it is assumed that beam is homogeneous and its behavior can be characterized by the following beam constitutive relations 1112130 122223 1323330PKKK MKKK VKKK (458) where [ K ] is the symmetric matrix of beam stiffness coefficients; 00,,and are the midplane axial strain, curvature and shear strain, respectively; P M and V are the axial force, bending moment and shear force resultants, resp ectively, on cross sec tion in the homogeneous beam, and can be determined as: 11n ee xx iPV a (459) 11n ee xz iVV a (460) 11n ee xx i M zV a (461) PAGE 81 81 where,ee x xxz are the bending and shear stress at the center of the ith element, a is the unit cells length, Ve is the volume of the ith element. The midplane deformations are relate d to the midplane axial displacement u0, transverse displacement w and rotation as 0 00,, u w x xx (462) Actually, K11, K12, K22 and K33 are similar to the laminate stiffness coefficients A11, B12, D22 and kA55, respectively. There is no equivalent for K13 and K23 in the laminate theory, because the layers are assumed to be orthotr opic, and they are rotated about the z axis only. However, such a coupling between inplane deformations a nd transverse shear deformations may exist in textile composites as the fibers are inclined to the xyplane unlike in the laminates. Further the present model cannot be used to pred ict stiffness prope rties involving the y direction, such as A12, A22, D12 and D22, etc. However, their effects are not be ing neglected because a state of plane strain parallel to the xzplane has been assumed, and hence ,yyxy etc. are all taken as equal to zero. The beam constitutive relation in Equation 462 can also be expressed in terms of compliance coefficients. 0111213 122223 0132333SSSP SSSM SSSV (463) Unit Cell Boundary Conditions Consider a cantilever beam as loaded in Fi gure 45A. Using the superposition principle, we have the cases of the beam loaded as shown in Figures 45(B, C, and D). In order to determine stiffness coefficients Kij in Equation 458, now we need to set up periodic boundary conditions applied to the unit cell in Figure 43 or Figure 44 according to each scheme of loading in Figure PAGE 82 82 45(B, C, D), actually applyi ng three linearly independent de formations to the unit cell separately. It means that the differences in displacements between corresponding points on the left and right faces of the unit cell will be equal to that in a homogeneous beam subjected to the same type of deformation. The three macrosc opic deformations applied to the unit cell are ( i ) a unit axial strain; ( ii ) a unit curvature accompanied by a transverse deflection such that the shear strain x z vanishes; and ( iii ) a unit transverse strain. No w we derive periodic boundary conditions for the cases of unit curvature and unit transverse strain, the derivation of periodic boundary conditions for the case of unit axial strain can be found in (Zhu et al, 1998) as seen in Table 42. Periodic boundary condition for unit curvature Consider the unit cell as shown in Figure 43, the periodic displacement boundary condition corresponding to unit curvature along the x axis (1 ) was derived (Karkkaimen and Sankar, 2006). All other curvature will be zeros. Curvatures in this case are defined as follows 222 221,0,0xyxzwww xzxz (464) The periodic displacement boundary condition will be derived from the definition of curvature along the xaxis. Integrating once with respect to x yields ()w x fz x (465) where() f zis an arbitrary function of z. Differentiation of this expression with respect to z, together with the requirement the x zis set as zero, indicates that 2 '()0xzw fz xz (466) PAGE 83 83 Due to the above expression,() f zmust therefore be a constant, since0xz Furthermore, this constant must be zero due to the specification that slope at the origin is zero, i.e., 00xw x Next, Equation 465 is integrated with respect to x, giving rise to the foll owing expression and an arbitrary function ()hz. 2() 2x whz (467) Note coordinate values for opposite faces of the unit cell (Figure 43) may be substitute in Equation 467. 2(,)() 2 (0,)() a wazhz wzhz (468) These are then subtracted from each other, eliminating the unknown h(z), and effectively prescribing the relative displacem ent on opposite faces that can be used to apply unit curvature. 2(,)(0,) 2a wazwz (469) However, a further boundary condition is requir ed to remove the transverse shear forces that will be present due to the application of this displacement. In this way, the mechanical effect of curvature is isolated. Thus th e requirement is that transverse shear strain is zero, which is defined as 0zxuw zx (470) This can be rearranged as below, and the value of slope is known from Equation 465 above (where the arbitrary function has been shown to be zero). uw x zx (471) PAGE 84 84 Similar to the above procedure, this expression can then be integrated with respect to z, evaluated with coordinate values of opposite faces, which are then subtracted from each other to specify the relative displacemen t that must be prescribed. uzxc (472) (,)(0,)uazuzaz (473) Thus, the periodic boundary condi tions as shown in Equation 469 and Equation 473 must be simultaneously applied to isolat e the effects of an applied unit cu rvature, as seen in Table 42. Periodic boundary condition for unit transverse shear strain We introduce the following coordinate system. The xcoordinate is taken along the length of the beam, zcoordinate along the height of the beam and the ycoordinate is taken along the width of the beam (Figure 45A). In a general beam theory, all applied loads and geometry are such that the displacement (u,v and w) along the coordinates (x,y and z) are only function of x and z coordinates (Figure 45A). Here we fu rther assume that the displacement v is identically zero. There are a number of beam theories that used to present the kinematics of deformation. The EulerBernoulli beam theory neglects both tran sverse shear and transverse normal effects i.e., deformation is due entirely to bending and in plane stretching. Th e next theory in the hierarchy of beam theories is the Timoshenko beam theory (Timoshenko, 1921), which is based on the displacement field. (,)() (,)()ouxzzx wxzwx (474) where () x denotes rotation of the cross section, wo is the transverse deflection of the point (x,0) on the midplane (i.e. z=0) of the beam. The Timoshe nko beam theory relaxes the normality assumption of the EulerBernoulli beam theory and includes a constant state of transverse shear strain with re spect to the thickness coordina te. The Timoshenko beam theory PAGE 85 85 requires shear correction factors, which de pend not only on the material and geometric parameters, but also on the loading and boundary conditions. In this section, based on the Timoshenko beam theory, we need to set up periodic boundary conditions applied to the unit cell in Figure 44 according to the scheme of loading in Figure 45D. It means that the differences in displacements between corresponding points on the left and right faces of the unit cell (Figure 44) will be equal to that in a homogeneous beam subjected to the same type of deformation. Th e slope of the displacem ent is made up of two effects, the rotation of the cross section plus the additional slope caused by the shear. The assumed constant shear strain component is expr essed in terms of displacement and rotation as shown in Equation 462. Invoking the Timoshenkos beam theory (Timoshenko, 1921), the periodic boundary conditions of the unit cell (Fi gure 44) can be described as (/2,)(0,) (/2,)(0,) 2 uazuzz a wazwz (475) where and are constants, because of the unit shear strain condition, we take 1 is be determined from the condition of the axia l force resultant on cross section, set as zero 0 N (476) The determination of can be obtained from the two computational models as below Model 1 The periodic boundary conditions below are applied to the unit cell (Figure 44) (/2,)(0,)0 (/2,)(0,) 2uazuz a wazwz (477) From this model we can compute the axial force resultant 1P by Equation 459. PAGE 86 86 Model 2 The periodic boundary conditions below are applied to the unit cell (Figure 44) (/2,)(0,) (/2,)(0,)0uazuzz wazwz (478) Similarly, from this model we can compute the axial force resultant 2P by Equation 459. Therefore, the function in Equation 475 can be obtai ned by the superposition method as: 1 12 20P PP P (479) Finally, the periodic boundary conditions of the unit cell for the three cases can be summarized in Table 42. Conclusions Most practical analysis of laminated com posite beams are based on first order shear deformation theories in the form of Timoshenko beam, which require shea r correction factors to correctly account for shear stiffn ess and transverse shear stresses. In contrast to the research effort that has been invested in the developmen t of shear correction factors for laminated plates, the derivation of shear correction factors for lami nated beams has received less attention, and the available solutions are limited to specific la yup configurations. In the current study, we presented a general expression to evaluate shear correction factors for general laminated beams. In this chapter, a truly micromechanical model has been presented to predict stiffness coefficients and shear correcti on factors for a class of compos ite beams including the textile composite beams which have been never discussed before. In our method, the unit cell of a beam was analyzed, and its periodic boundary conditions have been set up for each case of applied loadings. The current model has been verified by applying to isotropic, laminated, and PAGE 87 87 biomaterial beams. Comparing the results with th e beam theory, the derived analytical solution and the lamination theory are given in Chapte r 5. The current micromechanical model is straightforward and can be used to predict stiffness coefficients and shear correction factors for other complex composites. PAGE 88 88 Table 41. Geometric parameters of plainweave textile beam unit cell. Table 42. Periodic displacement boundary conditions for the beam unit cell. u(L,z)u(0,z) w(L, z)w(0, z) (i). Unit axial strain 01 a 0 (ii). Unit curvature 1 za 2/2a (iii). Unit shear strain 01 z a/2 Note that: L is the unit cell length, L = a for cases (i) and (ii), L = a/2 for the case (iii). Figure 41. Coordinate syst em for the laminated beam. a (mm) b (mm) d (mm) e1(mm) e2(mm) h(mm) 3.6 1.8 0.65 0.325 0.375 0.75 z x h PAGE 89 89 A B C Figure 42. Textile composite beam under shear force. A) with 4 unit cells. B) with 9 unit cells. C) with 45 unit cells. V L z, w x, u PAGE 90 90 Figure 43. Unit cell of thin walled textile composite beam. Figure 44. Half unit cell model of thin walled textile composite beam. PAGE 91 91 A B C D Figure 45. Cantilever beam. A) under ax ial force.B) under bending. C) shear force. P M V x, u M V L z, w y z h b P PAGE 92 92 CHAPTER 5 RESULTS AND DISCUSSION Introduction Based the theories presented in Chapters 24. In this chapter, we ha ve developed computer codes in MATLAB and FORTRAN to solve micromechanics problems. Examples are presented to illustrate the effectiveness of the current model, and it is validated by comparing the results with available analytical and numerical solu tions. The algorithm for ha ndling periodic boundary conditions and the treatment of ma terial discontinuity are verifi ed through out these examples. Numerical results of prediction of elastic cons tants of fiber composites, stiffness and shear correction factor predicti on of a class of composite beams are given along with discussion in this chapter. Timoshenko Beam Problem Based on the MLPG method as described in Ch apter 3, we have developed a computer code in MATLAB to solve problems in elastostatics. This code is validated by solving the classical problem of Timoshenko beam. Consider a beam of the length L subjected to a parabolic trac tion at the free end (Figure 51). The beam has characteristic height D and is considered to be of unit depth and is assumed to be in a state of plane stress. The exact so lution to this problem is given (Timoshenko and Goodier, 1970). 2 2(63)(2)() 64xpyD uLxxy EI (51) 2 223()(45)(3) 64ypDx uyLxLxx EI (52) Where the moment of inertia I of the beam is given by PAGE 93 93 312D I (53) The stresses corresponding to the displacem ents in Equation 51 and Equation 52 are () (,)xPLxy xy I (54) (,)0yxy (55) 2 2(,) 24xyPD x yy I (56) In consistent set of units, we take P= 1, E = 1000, D= 1, L= 8, = 0.25 and the penalty parameter = 108. A regular nodal mesh of 72 nodes (Fi gure 52) with 18 nodes along the xaxis and 4 nodes along the yaxis is used. The radius ro of the local subdomain is taken to be the distance between two neighboring nodes in yaxis. The Gaussian weight function is to be adopted with ri = 4ci, k = 1; the parameter ci in the Gaussian weight function is chosen as the distance to the third near est neighboring node in the ydirection. In the co mputation, 8x8 Gauss points are used in each localdomain s and 8 Gauss points are used along each section of s. Figures from 53 to 56 show the comparis on of displacements, stresses obtained from the MLPG method and those of the analytical solutions. We can s ee that the MLPG solutions are in excellent agreement w ith the analytical ones. TwoPhase Verification Problem Based on the theory which has been devel oped in Chapter 4, we have developed a computer code in FORTRAN to solve micromech anics problems. The second example is a twophase problem for investigati ng the material discontinuity. The current method results are compared to the FEM numerical solution (ABAQUSTM) for the twophase problem (Figure 57) In consistent set of units, the material properties used for PAGE 94 94 this study are 70,0.2;3.5,0.35ffmmEE and the fiber volume fraction Vf = 0.503, the boundary is represented as a unit sq uare, and the radius of fiber R = 0.8. The MLPG scheme with 183 nodes (Figure 57B) is used and the uniform load P = 1 is applied at its ri ght side as shown in Figure 57A state of plane stress normal to the x1x2 plane is assumed. Figures from 58 to 510 show the compar ison of radial and horizontal displacements between the current method and the FE method at the fibermatrix inte rface, and along lines y = 0 and x = 1. We realize that the current methods solu tions are very close to those of the FEM, which are computed using the ABAQUS software. Figures from 511 to 513 show the distribution of interfacial stresses at interface, we can see that the radial and the tangential stresses in the two materials are iden tical at the interface. The hoop stresses in the two materials along the interface are not equal as expected. The comparison of interfacial stresses with those from FEM is shown in Figure 514 and Figure 515. They show that the radial stresse s and the tangential stress for the two materials obtained by the MLPGmicromechanics are in good agreement with th e FEM results. The accuracy of stress computation obtained from the current method may be a good aspect in modeling the damage, especially progressive dama ge, which requires accurate description of the stress field in different phases. Figures from 516 to 520 show the distri bution of stresses in the unit cell for the case111 Similarly, Figures from 521 to 525 show the distribution of stresses in the unit cell for the case331 Both the cases show that the radial and tangential stresse s are continuous at the material interface, but there is a jump at th e interface for hoop stresses as expected (Figure 520) and (Figure 525). PAGE 95 95 Figure 526 shows the distribu tion of displacement in the z direction for the case of out of plane shear test (131 ) Figure 527 shows the comparison of the defo rmed shape of period ic boundary between the FEM and the current method. For the case of shear test in the x1x2 plane (121M), the full unit cell as shown in Figure 527A is used, and the unit cell is discreti zed with 577 nodes (Figure 527B), the periodic bounda ry conditions for the case 121M as shown in Table 31 are imposed on the unit cell. Figures from 528 to 530 show the deformed shape of the unit cell. We can see that the deformation shape of the unit cell obtained by using the algorithm of the periodic boundary condition treatment as shown in Chapter 3 are in good agreement with th at obtained from the FEM. Prediction of Stiffness Properties Example 1 The unidirectional fiber composite was assumed to have circular fibers packed in a square array. The fiber and matrix materials were assu med isotropic, and their properties for this study are shown in Table 51. MLPG mi cromechanical model is developed to model the unit cell of the unidirectional fiber composite. The Because of symmetry, only onequarter of the unit cell is modeled, the boundary is represented as a square w ith a side dimension of 1, and the radius of fiber R= 0.8. The displacements applied on the boundaries correspond ing to the cases 1122331,1,and 1 are shown in Table 31. The MLPG method with 183 nodes is used (Figure 57B). Figure 531 pres ents the tension tests in the x1 and x2 directions11221and1 by imposing a specified uniform displacement 1 u on the unit cell. The composite is assumed to be in a state of plane strain normal to the x1x2 plane. The microstresses are computed by the PAGE 96 96 MLPG method, and then the average stresses ar e computed using Equation 31, and finally the stiffness properties Cij can be also computed by Equation 32. Finally allows us to compute elastic constants as shown in Table 52 below. For the same convention as shown in the composite literatures, the x3, x1 and x2 directions in Figure 32 are changed to the x1, x2, and x3 ones, respectively. The results of the current method are compared with those from the rule of mixture and HalpinTsais equation (1969). For the analytic al solutions, the longitudinal Youngs modulus E1 given by the rule of mixture as follows 1 f fmmEEVEV (57) where Vf is the fiber volume fraction, Vm is the matrix volume fraction. The major Poissons ratio 12 is 12 f fmmVV (58) The HalpinTsai equations for the longitudina l Youngs modulus and major Poissons ratio are the same as those in the rule of mi xture. While the transverse Youngs modulus E2 by the HalpinTsai is given as 21 1 f mfV E EV (59) Where is the reinforcing factor depending on fibe r geometry, packing geometry and loading conditions. For a circular fiber in square array packing geometry, has the value of two. is given as 1f m f mE E E E (510) PAGE 97 97 The major Poissons ratio 12 is 12 f fmmVV (511) The HalpinTsai for the inplane shear modulus G12 is 121 1 f mfV G GV (512) 1f m f mG G G G (513) Appropriate values of for each of specific cases can be chosen as shown in Appendix A. Table 52 shows the comparison of elastic consta nts. The elastic constants computed using the current method are very close to those obtained by the anal ytical solution (HalpinTsai equations). The maximum error is 2.45% for ET, 1.16 % forTL 1.82% forLT The minimum error obtained is 0.14% for ET. Note that the HalpinTsai fo rmulas are not available forTT Table 4.4 shows the comparison of elastic cons tants between the HalpinTsai equations, the FEM solutions, and the current method. Compari ng with the FEM solutions, we can see that most of the elastic constants predicted from th e current method are very close with those from the HalpinTsai equations. Example 2 The analysis of the unit cell for predicting th e full three dimensional elastic constants including shear moduli of a unidirectional Eglass/ epoxy composite is presented in this example. The fiber and matrix materials of the composite we re assumed isotropic. The elastic constants for the glass fiber and epoxy matr ix are given in Table 53. PAGE 98 98 Table 52 shows the comparison of elastic c onstants between the current method and the HalpinTsai equations. The elas tic constants such as transv erse and longitudinal modulus, Poissons ratio computed using the current me thod are very close to those obtained by the HalpinTsai equations (Halpin and Tsai 1969). The maximum error is 3.51% for LT, 2.35 % for TL, 2.27% for ET, and the error obtained is 0.59% for EL. Note that the HalpinTsai formula is not available for TT. The current method also give s a reasonable prediction for GLT and GTT. Figure 532 and Figure 533 show the compar ison of the computed elastic constants ET and GLT under various fiber volume fraction Vf with those from the mechan ics of material approach (Kaw, 1997) and HalpinTsais eq uation (Halpin and Tsai, 1969) as well as experimental data (Tsai, 1964; Noyes and Jones, 1968). As shown in Figure 532, the experimental data for the transverse elastic modulus (Tsai, 1964) are quite higher compared to those obtaine d from the mechanics ap proach, and are very slightly lower than those predicted by the meshless local PetrovGalerkin micromechanical analysis. The HalpinTsais equation gives a good approximation due to the fact that the HalpinTsais equation is a semianalytical empirical function. The mechanics approach provides a low bound for the transverse elastic modulus. The slig ht difference between the experimental data and the current method can be due to the impe rfect adhesion between fiber and matrix, the possible separation of fiber in one bundle unde r tension, and some possible voids in the composite. Figure 533 shows that shear modulus GLT of the composite as a function of fiber volume fraction. The mechanics approach underestimat es the shear modulus, while the current method gives more reasonable results than the HalpinTs ais equation compared with the experimental PAGE 99 99 data (Noyes and Jones, 1968). No experimental data of other elas tic constants is available for comparison. Example 3 Table 56 shows the comparison of elastic cons tants obtained from the current method with those from FEM (Marrey and Sankar, 1995) and HalpinTsai equations for the Eglass/epoxy composite with constituent materi als as shown in Table 55. We realize that the current method seems to give a better prediction than FEM for the elastic constants, especially for EL, we can see that the current method and FEM gi ve more reasonable result for GLT than the HalpinTsais equation compared with the experimental data (Noyes and Jones, 1968). Stiffness and Shear Correction Factor Prediction of Composite Beams In Table 57, the shear correction factors fo r various laminated beam are computed by the micromechanical model using the unit cell in Figu re 535 and compared to analytical solutions given by Equation 451 in Chapter 4 and other investigators. Firs t, the shear correction factor obtained from the current model is verified by ap plying them to an isotropic beam, we can see that the shear correction factor from the micromechanical model is in excellent agreement with the standard value 5/6; the e rror only is 0.54%. Next, the shea r correction factor obtained by Equation 451 for a three layer sa ndwich beam is compared with the result obtained from the current micromechanical model, the good agreement is obtained. The finite element model of the unit cell of the three layer lami nated beam and its deformed shape under the periodic boundary conditions according to the case (iii) in Table 42 of Chapter 4 is shown in Figure 536A. In this example, the shear correction factor for the sandwich beam computed using Equation 551 is found to be k = 0.36721, thus the transverse shear stiffness of the bean is given by kA55 = 41.128 x 106N/m. From the sandwich beam theory (Whitney 1987) the transverse stiffness will be equal to the product of the core shear modulus a nd the core thickness whic h is equal to 32 x 106N/m. PAGE 100 100 The analytical solution by Equation 451 is hi gher because it takes into account the shear stiffness contribution of face sheets also. Finall y, the shear correction factors obtained from the current model are compared to the analytical solution (Whitney, 1973) for a bimaterial beam, and to the solution for reported by Bert (1973) for a two unsymmetr ic crossply laminated beam, as expected the shear correction factors have e rrors, the there are minimal and acceptable. Figure 536B shows the finite element mo del of the unit cell of the bimaterial beam and its deformed shape under the periodic boundary co nditions according to the case (iii) in Table 42 of Chapter 4. In general, the current micromechanical mode l can accurately predict shear correction factors for various types of laminated rectangular beam s as illustrated by the good agreement with the analytical solution and th e existing values reported by previous studies. The micromechanical model described above is used to predict stiffness coefficients for an isotropic beam, a biomaterial beam, and a thin walled textile composite beam. The dimensions of the unit cell and the yarn architect ure are described in Table 58 and in Figure 534. The unit cell dimensions (1.8mm x 1.8mm) as shown in Figure 535 are used for the isotropic and bimaterial cases, the properties of materials fo r all the cases are given in Table 59. In the case of the textile composite the matrix material is modeled as an isotropic material, while the yarn is modeled as a transverse isotropic materi al, with the 23 plane as the plane of isotropy. The yarn direction is assumed to be parallel to the 1 axis, and note that the 2 axis and y axis are parallel to each other. The results for all the cases are presented in Table 510, the results for the isotropic and biomaterial beams are compared with exact beam theory solutions. It may be seen from Table 510 that the current model is able to predict th e stiffness coefficients of the isotropic and PAGE 101 101 biomaterial beams, the axial and bending stiffn ess coefficients are predicted accurately. As shown above the shear stiffness predictions have errors but they are minimal and acceptable. Figure 537 shows the finite element model of the unit cells of the thin walled textile composite beam and their deformed shape under various independent loading conditions. For the case of the textile composite beam, at present no reliable experimental data and analytical solution are available to compare with the curr ent model. When shear tractions are allowed on the top and bottom surfaces of the unit cell (Figure 537C), us ing the micromechanics method (Dang and Sankar, 2007), which was presented in th e previous chapter, we can obtain the shear modulus Gxz and is found to be 2.503GPa, this yields A55 = Gxzh = 4.505 x 106Nm1, from Table 510 the shear stiffness coefficient K33 = 6.68 x 106Nm1, we have the shear correction factor of the textile composite beam k = K33/A55 = 1.483. From K11 of the textile beam one may extract the Youngs modulus Ex as K11/h, which will yield Ex = 12.89 GPa. If this value is used to predict the flexural stiffness of a homogeneous beam as D11= Exh3/12, one can obtain D11 = 6.26Nm as a upper bound value, whereas the actual fl exural stiffness is equal to 4.567Nm. Once can estimate a lower bound value of D11 by considering the textile beam u nder bending as shown in Figure 538. The textile beam is created by a large number of the unit cells enough to be considered as a homogeneous beam, we can obtain D11 from the beam theory2/(2)tipEIMLV for an example the textile beam with 45 unit cells (Figure 538B), D11 is estimated to be 4.329Nm. We can see that the predicted value of K11 in Table 5.6 is reasonable, and the current model gives a good estimate. This further illustrates the importan ce of the current analysis in predicting beam stiffness properties directly. PAGE 102 102Conclusions The technique of direct impos ition of interface boundary conditi ons is used in the current model for the treatment of the material discontin uity and the algorithm for treatment of periodic boundary conditions in the MLPG method using multipoint constraints are validated by comparing the results with available analytical and numerical solutions. The current method shows good agreement with the FEM solution, and achieves significantly higher degree of agreement of the radial and the tangential stresse s in the two materials at the interface with smaller number of degrees of fr eedom compared to the FEM. The elastic constants obtained by the MLPG micromechanical analysis match well with HalpinTsais and experimental data. The curre nt method seems to give a better prediction for the elastic constants than FEM. In general the current method and FEM give better estimates for GLT than HalpinTsais equations. The micromechanical model suggested in Ch apter 4 to predict stiffness and shear correction factors for composite beams is straig htforward and is a truly micromechanical one. Prediction of the flexural stiffness properties and the shear correction factor of the isotropic and laminated beams using both the finite element me thod and the meshless lo cal PetrovGalerkin (MLPG) method are in very good agreement with th e analytical solutions. Results and discussion are also given for a thin walled textile com posite beam using the current micromechanical model. PAGE 103 103 Table 51. Material properties of Eglass/PP composite. Property Unit EGlass Epoxy Tensile modulus GPa 70 3.5 Poissons ratio 0.2 0.35 Axial shear modulus Gpa 29.167 1.296 Volume fraction 0.503 0.497 Table 52. Comparison of elastic constants between the current method and HalpinTsai equations. Moduli are in GPa. Method E1 E2 E3 12 21 13 31 23 32 Current method 36.90 11.848 11.848 0.280 0.0850.280 0.085 0.3150.315 HalpinTsai 36.95 11.565 11.565 0.275 0.0860.275 0.086 Table 53. Material properties of Eglass/epoxy composite. Property EGlass Epoxy Axial modulus (GPa) 73.1 3.45 Transverse modulus (GPa) 73.1 3.45 Axial Poissons ratio 0.22 0.35 Transverse Possons ratio 0.22 0.35 Axial shear modulus (GPa) 30.19 1.83 Volume fraction 0.503 0.497 Table 54. Comparison of elastic constants between the current method and HalpinTsai equations. Elastic Constants HalpinTsa i equations Current Method EL(GPa) 38.484 38.256 ET(GPa) 11.514 11.253 LT 0.285 0.275 TL 0.085 0.083 TT 0.322 GLT (GPa) 4.771 5.322 GTT (GPa) 3.06 2.897 Table 55. Material properties of Eglass/PP composite. Property EGlass Epoxy Tensile modulus (GPa) 70 3.5 Poissons ratio 0.2 0.35 Axial shear modulus (GPa) 29.167 1.296 Volume fraction 0.6 0.4 PAGE 104 104 Table 56. Comparison of elastic constants. Elastic Constants HalpinTsai FEM Current method EL(GPa) 43.40 43.12 43.658 ET(GPa) 14.792 18.15 13.718 LT 0.260 0.242 0.268 TL 0.252 0.222 0.247 TT 0.102 0.083 GLT (GPa) 4.451 5.590 5.390 GTT (GPa) 3.860 3.92 3.796 Table 57. Comparison of shear correction factors for various t ypes of laminates. Shear correction factor Current method Laminate type and ref. FEM MLPG Reference value % difference Isotropic beam E= 10GPa, =0.3 0.8378 0.8345 0.8333 0.54 3Layer laminated beam Es=100GPa, Gc= 10GPa, = 0.25, s= 1mm, c= 8mm 0.3617 0.3634 0.3672 1.5 Bimaterial beam E1= 100GPa, 1= 0.3, 1= 0.9mm E2= 10GPa, 2= 0.3, 2= 0.9mm 0.4906 0.506 0.510 3.8 2Layer unsym. crossply beam EL= 25x106 psi, ET= 106 psi, LT= 0.25 GLT= 0.5x106 psi, GTT= 0.2x106 psi 0.7873 0.8068 0.8212 4.13 Table 58. Geometric parameters of plainweave textile beam unit cell. Table 59. Constituent material properties for beam examples. Material 1E(GPa) 2E(GPa) 12G(GPa) 12 23 Isotropic beam 10 10 3.85 0.3 0.3 Material 1 100 100 38.46 0.3 0.3 Bimaterial beam Material 2 10 10 3.85 0.3 0.3 yarn 159 10.9 6.4 0.38 0.38 Textile beam matrix 3.5 3.5 1.3 0.35 0.35 a (mm) b (mm) d (mm) e1(mm) e2(mm) h(mm) 3.6 1.8 0.65 0.325 0.375 0.75 PAGE 105 105 Table 510. Comparison of beam s tiffness coefficients (SI units). Beam Model K11 K12 K22 K33 FEM 19.78x 1060 5.35 5.81x106 Current model MLPG 19.78x106 0 5.33 5.79x106 Isotropic beam Beam theory 19.78x 106 0 5.34 5.77x106 FEM 10.88x107 40.05x103 29.37 18.68x106 Current model MLPG 10.89x107 40.05x103 29.38 19.27x106 Bimaterial beam Beam theory 10.88x107 40.05x103 29.37 19.42x106 FEM 23.2x106 0 4.567 6.68x106 Textile beam Current model MLPG 24.36x106 0 4.403 6.59x106 Figure 51. Timoshenko beam. Figure 52. Nodal scheme with 72 nodes. x u y D= 1 L=8 P=1 PAGE 106 106 Figure 53. Comparis on of displacements U & V at the line y= 0.5. Figure 54. Comparison of normal stresses x at the line y= 0.5. PAGE 107 107 Figure 55. Comparison of normal stresses x at the line x= 4. Figure 56. Comparison of shear stresses xy at the line x= 4. PAGE 108 108 A B C Figure 57. Two phase problem. A) An illustration of tension test in x1 (111M ). B) Mesh scheme of 183 nodes. C) Radial, tangential and hoop stress at interface. Material 1 Material 2 P I PAGE 109 109 0.00E+00 1.00E04 2.00E04 3.00E04 4.00E04 5.00E04 6.00E04 7.00E04 8.00E04 9.00E04 020406080100Angle (in degree)Radial displacement Current method FEM Figure 58. Comparison of radial displacements at the interface. 2.00E04 0.00E+00 2.00E04 4.00E04 6.00E04 8.00E04 1.00E03 1.20E03 1.40E03 0.00E+002.00E014.00E016.00E018.00E011.00E+001.20E+00XHorizontal displacement Current method FEM Figure 59. Comparison of horiz ontal displacements along line y =0. PAGE 110 110 0 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 0.0014 0.0016 0.0018 0.00E+002.00E014.00E016.00E018.00E011.00E+001.20E+00YRadial displacement Current method FEM Figure 510. Comparison of radi al displacements along line x =1. 0 0.2 0.4 0.6 0.8 1 1.2 020406080100Angle (in degree)Radial stress Material 1 Material 2 Figure 511. Radial stresses in the two mate rials at the interface by MLPGbased micromechanics. PAGE 111 111 0 0.2 0.4 0.6 0.8 1 1.2 1.4 020406080100Angle (in degree)Hoop stress Material 1 Material 2 Figure 512. Hoop stresses in the two materials at th e interface by MLPGbased micromechanics. 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0.05 020406080100Angle (in degree)Tangential stress Material 1 Material 2 Figure 513. Tangential (shear) st resses in the two materials at the interface by MLPGbased micromechanics. PAGE 112 112 0 0.2 0.4 0.6 0.8 1 1.2 020406080100Angle (in degree)Radial stress Mat 1: Current method Mat 2: Current method Mat 1: FEM Mat 2: FEM Figure 514. Comparison of radi al stresses at the interface. 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0.05 020406080100Angle (in degree)Tangential stress Mat 1: Current method Mat 2: Current method Mat 1: FEM Mat 2: FEM Figure 515. Comparison of tangent ial stresses at the interface. Distribution of stresses over domain in the case 111M (conventional plane strain problem) PAGE 113 113 Figure 516. Distribution of x ( x is scaled by 1/50) Figure 517. Distribution of y (y is scaled by 1/10) Figure 518. Distribution of z (z is scaled by 1/10) Figure 519. Distribution of x y ( x y is scaled by 1/3) PAGE 114 114 Figure 520. Distribution of in terfacial stresses for the case11 =1 (r is scaled by 1/50, r and r are scaled by 1/30). Distribution of stresses over domain in the case 331M (generalized plane strain problem) Figure 521. Distribution of x ( x is scaled by 1/10) Figure 522. Distribution of y (y is scaled by 1/10) PAGE 115 115 Figure 523. Distribution of x y ( x y is scaled by 1/5) Figure 524. Distribution of z (z is scaled by 1/100) Figure 525. Distribution of in terfacial stresses for the case331 (r is scaled by 1/7,r andr are scaled by 1/10). PAGE 116 116 Figure 526. Displacement along zdi rection over domain for the case131 A B Figure 527. Unit cell for shear test (121M ) in the x1x2 plane. A) Unit cell and its support. B) Nodal mesh (577 nodes). R Fiber Matrix 2a 2a PAGE 117 117 0.2 0 0.2 0.4 0.6 0.8 1 0.200.20.40.60.811.21.41.61.8XY Current method FEM Figure 528. Deformed boundary of unit cell for case of121M (not to scale). A PAGE 118 118 B Figure 529. Deformed unit cell for case of121M and shear stress distribution by FEM. A) Not to scale. B) Scale 1:1. 1.50 1.00 0.50 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 1.501.000.500.000.501.001.50 XY FEM Current method Figure 530. Periodic deformation of unit cell for shear test (121M ) by the current method. PAGE 119 119 Figure 531. Tension tests in x1 and x2 directions (111M and221M ). 0 5 10 15 20 25 30 0.20.30.40.50.60.7Fiber Volume FractionTransverse modulus (GPa) Mechanics approach Current Method HT Equation Experimental Figure 532. Comparison of transverse elastic modulus ET of the composite. u = Fiber Matrix Rigid R a a Fiber Matrix Rigid a R a u = PAGE 120 120 0 1 2 3 4 5 6 7 8 9 10 0.40.450.50.550.60.650.7 Fiber volume fractionShear modulus (GPa) Experimental Mechanics approach HT Equation Current method Figure 533. Comparison of shear modulus GLT of the composite. Figure 534. Unit cell of thin walled textile composite beam. PAGE 121 121 Figure 535. Half unit cell model of thin walled textile composite beam. B A Figure 536. Deformed unit cell on shear problems. A) 3layer sy mmetric laminated beam. B) bimaterial beam. PAGE 122 122 A B C D Figure 537. Deformed unit cell. A) unit axial strain. B) unit curvature. C) shear problem. D) periodic BCs allowed on top and bottom surfaces. PAGE 123 123 A B Figure 538. Textile composite beam under bending. A) 9 unit cells. B) 45 unit cells. M M PAGE 124 124 CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK The purpose of this chapter is to summarize th e work done and to draw the findings of the work presented in this dissertation as a whole, to highlight the significan t contributions and to make recommendations for further studies in th e field, based on percei ved limitations of the technique suggested. Summary The meshless local Petrov Galerkin (MLPG) mi cromechanical analysis has been presented in this study for problems of composite micromec hanics; the micromechanical model is used to predict the elastic constants of structural composites from the analysis of the representative volume element (RVE). Micromechanical analysis of structural composites is possible due to the presence of RVE (unit cell). The unit cell is a ssumed to span the material continuously in all three dimensions. On the microscale, comparable to the scale of unit cell, the composite is heterogeneous due to the presence of the rein forcing yarn and the matrix. However, on the macroscale, compared to the structural scale, the composite is assumed to be homogeneous and orthotropic. The homogeneous composite propert ies are then predicted from the constituent material properties and the yarn geometry. The unit cell is discretized and periodic boundary conditions are imposed between opposite endfaces of the unit cell, essential boundary co nditions are enforced by the penalty method. The treatment of material discontinuity at the interf ace between different phase s of the composite is presented by means of direct imposition of interface boundary c onditions, and an algorithm for handling of periodic boundary conditions in the MLPG method using multipoint constraint elements is also presented. The MLPG micromech anical model is presented for the six linearly PAGE 125 125 independent deformations of the unit cell. From the forces acting on the unit cell for each of the six cases, the elastic constants of the composite can be computed. Also in this study, we introduced a theoretica l framework to predict the flexural stiffness properties and the shear correction factor for a class of co mposite beam problems with emphasis on textile composite beams. The suggested model is straightforward and a truly micromechanical model. Major Conclusions and Contributions A number of conclusions were drawn at th e end of each chapter. The most important conclusions from each chapter are listed: The technique of direct impos ition of interface boundary conditions is used for the first time for the treatment of the material discont inuity at the interface in the MLPG method. Current method shows good agreement with the FEM solution, and achieves significantly higher degree of agreement between the radial and tangential stresses in the two materials at the interface with smalle r number of degrees of freedom compared to the FEM. The MLPG method is formulated for the gene ralized plane strain problems and out of plane shear problems for the first time The treatment of periodic boundary conditions in the MLPG method using multipoint constraint technique is presented. The elastic constants obtained by the meshle ss local PetrovGalerkin micromechanical model match very well with available results. The MLPG micromechanical model gives a bette r prediction for the elastic constants than FEM. In general the current model and FEM give better estimates for GLT than HalpinTsais equations. The theoretical framework of the truly micromechanical model based on both the finite element method and the meshless local PetrovGalerkin method to predict the flexural stiffness properties and the shear correction factor for composite beams is suggested. Prediction of the flexural stiffness properties a nd the shear correction f actor of the isotropic and laminated beams obtained from the FE M and MLPG methodbased micromechanical model are in very good agreement with th e analytical soluti ons. The current micromechanical model is straightforward a nd can be easily used to predict the shear correction factor for textile structural co mposites and other complex composites which have been never discussed before. PAGE 126 126 A computer code is developed in FO RTRAN and MATLAB, based on the MLPG formulation for problems in composite micromechanics The MLPG methodbased micromechanical mo del is a truly meshless method, wherein no elements or background cells are involved, either in the interpolation or in the integration. Recommendations for Further Work While the major conclusions and achievements obtained as shown in the preceding section, a number of unsolved areas have become appa rent where further investigation would be beneficial. These are summarized below. One of the major factors influencing the succ ess of a methodology is the cost vs. accuracy tradeoff. Thus effective implementations of meshless methods are key to their success especially in composite micromechanic s as unit cell calcula tions are usually computationally intensive (Cox and Flanagan, 1997), the output of unit cell calculations includes both the spatially averag ed response of the cell and de tails of stress distributions in yarns and matrix. Comparison of computa tional cost between a meshless method and an FE solution with same number of unknowns in loworder finite elem ents was carried out by Belytschko et al (1996), and A tluri and Shen (2002). It has b een seen that the FE results are in general less expensive. However, comparing the costs based on a given level of accuracy or if highorder finite elements ar e compared, then the results can be quite different. The effort of research ers in improving the meshless met hods in this aspect is still on going. The latest effort devoted to impr oving the effectiveness of the MLPG method can be found in Atluri and Shen (2002) and A tluri (2004) in which the authors introduce MLPG5 using the Heaviside functio n as the test function in e ach domain. The test function is defined as s s1atx () 0atxvx (514) We start with the local symmetric weak fo rm (Equation 221). It is rewritten as ,()()0ssustssuiiiiiiijijiiiii Ltvdtvdtvdvbvduuvd (515) Because this test function does not vanish over the local boundary Ls, rearranging the terms, we obtain the local symmetric weak form as ,sssusustsusijijiiiiiiiiiiii Lvdtvdtvduvdtvduvdbvd (516) With the choice of the test function as Equation 61, this leads to the derivatives of the test functions vanished on s and Equation 63 is simplified as: PAGE 127 127 ssusustsusiiiiiiiiiiii Ltvdtvduvdtvduvdbvd (517) It is seen that in Equation 64, there is no domain integrati on involved in the left hand side, which leads to the stiffness matrix, after discretization. If the assumption of zero body force is made, then domain integration is tota lly eliminated. Thus, the domain integral on each domain is altogether avoided in compu ting the stiffness matrix. It involves only boundary integrals over each circle (domain), which will greatly improve the effectiveness of the MLPG5 method and will make the solution stable, fast and accurate. It is reported that the MLPG5 provides a simple and effici ent alternative to the finite element and boundary element methods (Atluri 2004). We st rongly recommend using the type of this MLPG method with the theoreti cal framework as shown thro ugh out the dissertation to problems in composite micromechanics, especi ally in 3D problems. Studies in these directions will have considerable impact. The MLPG methodbased micromechanical mode l in general is a truly meshless method, wherein no elements or background cells are involved, either in the interpolati on or in the integration. However as computing the m acrostresses (Equation 31), the axial force (Equation 459), the bending moment (Equati on 460), and the shear force resultants (Equation 461), we still use connectivity in our computer code to calculate these integrals, this may lead to lightly increase the comput ational cost. We realized that this is a shortcoming in our work, overcom e that would be beneficial. The current micromechanical model could also be employed to perform a parametric study of the effect of textile architecture, geomet ry and material properties on shear correction factors. For example types of textile compos ite beam with varying spacing and different dimensions, material properties ratio, and ya rn orientation angle, etc. could provide a useful insight into the effect of these para meters have on the stiffness coefficients and shear correction factors. This could also potentially lead to the ability to optimize the microarchitecture to a specific application. We hope that our work will be of some usefulness to those who want to do research in the fascinating field of co mposite micromechanics. PAGE 128 128 APPENDIX IMPORTANT NOTES FOR USIN G HALPINTSAI EQUATIONS The longitudinal Youngs modulus E1 given by the mechanics approach is 1 f fmmEEVEV (A1) where Vf is the fiber volume fraction, Vm is the matrix volume fraction. The major Poissons ratio 12 is 12 f fmmVV (A2) 1 1 f mfV M M V (A3) where (/)1 (/)fm fmMM MM (A4) in which M = composite material modulus E2, G12, or 23 Mf = corresponding fiber modulus E2, G12, or 23 Mm= corresponding matrix modulus E2, G12, or 23 is a measure of fiber reinforcement of the composite material and is dependent on fiber geometry, packing geometry, and loading conditions. Obtaining accurate values for is the most difficult part of using the HT equations. Approximations for Circular fibers in square array (excellent results for Vf~0.55) = 2 for E2 PAGE 129 129 = 1 for G12 Square fibers in square array (good results for Vf <=0.9) = 2 for E2 = 1 for G12 Rectangular cross section in diamond array (good results for Vf <=0.9) = 2a/b for E2 = 1.73log(a/b) for G12 where a= width of fiber, and b= thickness of fiber For the limiting cases When = 0 m m f fM V M V M 1 (Series Model) When M= VfMf + VmMm (parallel model) PAGE 130 130 LIST OF REFERENCES Adams, D.F., Doner, D.R., 1967. Longitudinal sh ear loading of a unidi rectional composites. Journal of Composite Materials 1, 4. Adams, D.F., Crane, D.A., 1984. Finite element micromechanical analysis of a unidirectional composite including longitudinal shear load ing. Computers & St ructures 18, 1153. Atluri, S.N., Zhu, T., 1998a. A new meshless lo cal PetrovGalerkin (MLPG) approach to nonlinear problems in comput ational modeling and simula tion. Computer Modeling & Simulation in Engineering 3, 187. Atluri, S.N., Zhu, T., 1998b. A new meshless lo cal PetrovGalerkin (MLPG) approach in computational mechanics. Computational Mechanics 22(2), 117. Atluri, S.N., Kim, H.G., Cho J.Y., 1999. A critical assessment of the truly meshless local PetrovGalerkin (MLPG), and local boundary integral equation (LBIE) me thods. Computational Mechanics 24(5), 349. Atluri, S.N., Zhu T.L., 2000. The meshless local Pe trovGalerkin (MLPG) approach for solving problems in elastostatics. Com putational Mechanics 25(23), 169. Atluri, S.N., Shen, S., 2002. The meshless local PetrovGalerkin (MLPG) method: A simple & lesscostly alternatives to the finite element and boundary element methods. CMES: Computer Modeling in Engin eering & Sciences 3, 11. Atluri, S.N., 2004. The meshless method (MLPG) for domain & BIE discretizations. Tech Science Press. Babuska, I., Melenk, J.M., 1995. The partition of unity finite element method. Technical Report BN1185, Institute for Physics, Science and Technology, University of Maryland, Maryland. Babuska, I., Melenk, J.M., 1996. The partition of un ity finite element method: Basic theory and application. Computer Me thods in Applied Mechanics and Engineeing 139, 3. Barbero, E.J., LopezAnido, R., Davalos, J.F., 1993. On the mechanics of thinwalled laminated composite beams. Journal of Composite Materials 27, 806. Batra, R.C., Ching, H.K., 2002. Analysis of el astodynamic deformations near a crack/notch tip by the meshless local PetrovGalerkin (MLP G) method. CMES: Computer Modeling in Engineering & Sciences 3(6), 717. Batra, R.C., Porfiri, M., Spinello, D., 2004. Treat ment of material discont inuity in two meshless local PetrovGalerkin (MLPG) formulations of axisymmetr ic transient heat conduction. International Journal for Numerical Me thods in Engineering 61(14), 2461. PAGE 131 131 Belytschko, T., Lu, Y.Y., Gu, L., 1994. Elementfree Galerkin methods. International Journal for Numerical Methods in Engineering 37, 229. Belytschko, T., Krongauz, Y., Organ, D., Flemi ng, M., Krysl, P., 1996. Meshless Methods: An Overview and Recent Developments. Comput er Methods in Applied Mechanics and Engineering 139, 3. Bert, C.W., 1973. Simplified analysis of stat ic shear factors for beams of nonhomogeneous crosssection. Journal of Co mposite Materials 7, 525. Boeing Press Release, 2nd Oct ober 1996. New Boeing composite st ructure leads the ways to lighter, cheaper satellites. Boley, B.A., Tolins, I.S., 1956. On the stresses an d deflections of rectangular beams. Journal of Applied Mechanics 23, 393. Cai, Y.C., Zhu, H.H., 2004. Di rect imposition of essential bounda ry conditions and treatment of material discontinuities in the EFG me thod. Computational Mechanics 34, 330. Ching, H.K., Batra, R.C., 2001. Determination of crack tip fields in line ar elastostatics by the meshless local PetrovGalerkin (MLPG) method. CMES: Computer Modeling in Engineering & Sciences 2(2), 273. Chou, T.W., Ishikawa, T., 1983. One dimensional micromechanical analysis of woven fabric composites. AIAA Journal 21, 1714. Chou, T.W., Ko, F.K., 1989. Textile structural co mposites. Composite Materials Series, Vol. 3, Elsevier, New York. Chow, T.S., 1971. On the propagation of flexural waves in an ort hotropic laminated plate and its response to an impulsive load. Journa l of Composite Materials 5, 306. Christensen, R., 1990. A critical evaluation for a class of micromechanics models. Journal of the Mechanics and Physics of Solids 38(3), 379. Cordes, L.W., Moran, B., 1996. Treatment of materi al discontinuity in the elementfree Galerkin method. Computers Methods in Applie d Mechanics and Engineering 139, 75. Cowper, G.R., 1966. The shear coefficient in Ti moshenko's beam theory. Journal of Applied Mechanics 33, 335. Cox, B.N., Flanagan, G., 1997. Handbook of analy tical methods for textile composites. NASA CR4750. PAGE 132 132 Dang, T., Sankar, B.V., 2005. A Meshless Loca l PetrovGalerkin (MLPG) Micromechanical Model for Fiber Composites. In: Proceeding of 20th Annual Technical Conference of the American Society for Composites, Dr exel University, Philadelphia, PA. Dang, T., Sankar, B.V., 2006. Meshless Local Pe trovGalerkin Micromechanical Analysis of Fiber Composites Including Axial Shear Load ing. 7th World Congress on Computational Mechanics (WCCM VII), Los Angeles, California, Paper Number 489. Dang, T., Sankar, B.V., 2007. Meshless Local Pe trovGalerkin Formulation for Problems in Composite Micromechanics. AIAA Journal 45(4), 912. Dharmarajan, S., McCutchen, H., 1973. Shear coe fficients for orthotropic beams. Journal of Composite Materials 7, 530. Dolbow, J., Belytschko, T., 1998. An introduc tion to programming the meshless element free Galerkin method. Archives in Co mputational Mechanics 5(3), 207. Duarte, C.A.M., Oden, J.T., 1995. Hp cloudsMeshless method to solve boundary value problems. TICAM Report, 95. Duarte, C.A.M., Oden, J.T., 1996. An hp adapti ve method using clouds. Computer Methods in Applied Mechanics and Engineeing 139, 237. Gibson, R.F., 1994. Principles of composite mate rial mechanics. New York, McGrawHill, Inc. Gingold, R.A., Monaghan, J.J., 1977. Smoothed pa rticle hydrodynamics: Theory and application to nonspherical stars. Royal Astronomi cal Society, Monthly Notices 181, 375. Halpin, J.C., Tsai, S.W., 1969. Effect of environm ental factors on composite materials. Air Force Technical Report AFMLTR. Hashin, Z., 1962. The elastic moduli of heterogene ous materials. Journal of Applied Mechanics 29, 143. Hutchinson, J.R., 2001. Shear coefficients for Timoshenko beam theory. Journal of Applied Mechanics 68, 87. Ishikawa, T., Chou, T.W., 1982. Stiffness and st rength behavior of woven fabric composites. Journal of Materials Science 17, 3211. Ishikawa, T., Chou, T.W., 1982. Elastic behavior of woven hybrid composites. Journal of Composite Materials 16, 2. Karkkainen, R.L., Sankar, B.V., 2006. A direct micromechanics method for analysis of failure initiation of plain weave textile composite s. Composites Science and Technology, 66(1), 137. PAGE 133 133 Kaw, A.K., 1997. Mechanics of compos ite materials. New York, CRC Press. Kim, S.J., Yoon, K.W., Jung, S.N., 1996. Shear correction factors for thin walled composite boxbeam considering nonclassical behavior. Journal of composite materials 30, 1132 1149. Kim, H.J., Swan, C.C., 2003. Voxelbased meshi ng and unitcell analysis of textile composites. International Journal for Numerical Methods in Engineering 56, 977. Krongauz, Y., Belytschko, T., 1998. EFG approxi mation with discontinuous derivatives. International Journal for Numerical Methods in Engineering 41, 1215. Lancaster, P., Salkauskas, K., 1981. Surface generated by moving least squares methods. Mathematics of Computation 37(155), 141. Li, Q., Shen, S., Han, Z.D, Atluri, S.N., 2003. Application of Meshless Local PetrovGalerkin (MLPG) to Problems with Singularities, and Ma terial Discontinuities, in 3D Elasticity. CMES: Computer Modeling in Engi neering & Sciences 4, 571. Liu, W.K., Chen, Y., Uras, R.A., 1995a. Enrich ment of the finite element method with the reproducing kernel particle method. Current Topics in Computational Mechanics, ed. Cory, J.F., Gordon, J.L., PVP305, ASME, New York, 253. Liu, W.K., Jun, S., Zhang, Y.F., 1995b. Repr oducing kernel particle methods. International Journal for Numerical Met hods in Fluid 20, 1081. Listzka, T., Orkisz, J., 1980. The finite difference method for arbitrary meshes. Computers & Structures 5, 45. Lu, Y.Y., Belytschko, T., Gu, L., 1994. A new implementation of the element free Galerkin method. Computer Methods in Applie d Mechanics and Engineering 113, 397. Lucy, L.B., 1977. A numerical approach to the testing of the fission hyp othesis. Astronomical Journal 8(12), 1013. Marrey, R.V., and Sankar, B.V., 1995. Microm echanical models for textile structural composites. NASA CR198229. Marrey, R.V., Sankar, B.V., 1997. A micromechan ical model for textile composite plates. Journal of Composites Ma terials 31(12), 1187. Mechnik, R.P., 1991. Considerati on of constraints within the fi nite element method by means of matrix operators. Internati onal Journal for Numerical Methods in Engineering 31, 909 926. PAGE 134 134 Monaghan, J.J., 1982. Why particle method works. SIAM Journal on Scientific and Statistical Computing 3(4), 422. Monaghan, J.J., 1988. An introduction to SPH Computer Physics Communications 48, 89. Monaghan, J.J., 1992. Smooth particle hydrodyn amics. Annual Review of Astronomy and Astrophysics 30, 543. Naik, R.A., 1994. Analysis of woven and brai ded fabric reinforced composites. NASA CR194930. Nayroles, B., Touzot, G., Willon, P., 1992. Gene ralizing the finite element method: Diffuse approximation and diffuse elements. Computational Mechanics 10, 307. Noor, A.K., Peters, J.M., 1989. A posteriori estimates for shear correction factors in multilayered composite cylinders. Journal of Engineering Mechanics 115 (6), 1225. Noyes, J.V., Jones, B.H., 1968. Analytical desi gn procedures for the strength and elastic properties of multilayer fiber composites. Proceedings of AIAA/ASME 9th Structure Dynamics and Material Conference, 68. Onate, E., Idelsohn, S., Zienkiewicz, O.C., Taylor, R.L., 1996a. Finite point method in computational mechanics. Application to conve ctive transport and fluid flow. International Journal for Numerical Methods in Engineering 39, 3839. Onate, E., Idelsohn, S., Zienkiewicz, O.C., Tayl or, R.L., Sacco, C., 1996b. Stabilized finite point method for analysis of fluid mechanics problems. Computer Methods in Applied Mechanics and Engineeing 139, 315. Peng, X.Q., Cao, J., 2000. Numerical determinati on of mechanical elasti c constants of textile composites. 15th Annual Technical Conference of the American Society for Composites, College Station, Texas, September, 25. Poe, C.C., and Harris, C.E., 1995. Mechanic s of textile composites conference. NASA Conference Publication 3311, Part 1 &2. Qian, S., Weiss, J., 1993. Wavelet and the numer ical solution of partial differential equations. Journal of Computational Physics 106(1), 155. Raju, I.S., Chen, T., 2001. Meshless PetrovGalerk in method applied to axisymmetric problems. 42nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit, Seattle, WA, Apr. 16. Raju, I.S., Phillips, D.R., 2003. Local coordinate approach in meshless local PetrovGalerkin method for beam problems. AIAA Journal 41(5), 975. PAGE 135 135 Sankar, B.V., 1994. An exact solution for stresses in laminated beams. In ternational Journal for Engineering Analysis and Design 1, 293. Sankar, B.V., Marrey, R.V., 1993. A Unit Ce ll Model of Textile Composite Beams for Predicting Stiffness Properties. Comp osites Science and Technology 49(1), 61. Sankar, B.V., Marrey, R.V., 1997. Analytical me thod for micromechanics of textile composites. Composites Science and Technology 57(6), 703. Sapountzakis, E. J., Mokos, V.G., 2006. Dynamic analysis of 3D beam elements including warping and shear deformation effects. Intern ational Journal of Solid s and Structures 43, 6707. Shephard, M.S., 1984. Linear multipoint constraints applied via transformation as part of a direct stiffness assembly process. International J ournal for Numerical Me thods in Engineering 20, 2107 Sukumar, N., Moran, B., Belytschko, T., 1998. The natural element method in solid mechanics. International Journal for Numerical Methods in Engineering 43, 839. Sulsky, D., Chen, Z., Schhreyer, H.L., 1992. The application of a materi alspatial numerical method to penetration. New Methods in Transien t Analysis, eds. Smolin ski, P., Liu, W.K., Hulbert, G., Tamma, K., AMD143/PVP246, ASME, New York, 91. Takano, N., Uetsuji, Y., Kashiwagi, Y., Zako, M., 1999. Hierarchical modeling of textile composite material and structures by the homogenization method. Modelling and Simulation in Materials Scie nce and Engineering 7, 207. Timoshenko, S.P., 1921. On the correction for shea r of the differential equation for transverse vibration of bars of prismatic ba rs. Philosophical Magazine 41, 744. Timoshenko, S.P., Goodier J.N., 1970. Theory of Elasticity. 3rd ed. New York, McGraw Hill. Tsai, S.W., 1964. Structural behavior of composite materials. NASA CR71. Wendland, H., 1999. Meshless Galerkin methods us ing radial basis function. Mathematics of Computation 68(228). 1521. Whitcomb, J.D., 1991. Threedimensional stress an alysis of plain weave composites. Composite Materials Fatigue and Fracture 3, ASTM STP 1110, 417. Whitney, J.M., 1972. Stress analysis of thick la minated composite and sa ndwich plates. Journal of Composite Materials 6, 426. Whitney, J.M., 1973. Shear correction factors for or thotropic laminates under static load. Journal of Applied Mechanics 40, 302. PAGE 136 136 Whitney, J.M., 1987. Structural an alysis of laminated anisotropi c plates, Technomic Publishing Co., Lancaster, Pennsylvania, USA. Yi,Y.M., Park, S.H., Youn, S.K., 1998. Asympto tic homogenization of viscoelastic composites with periodic microstructures. International Jour nal of Solids and Structures 35, 2039 2055. Yoshino, T., Ohtsuka, T., 1982. Inner stress anal ysis of plane woven fiber reinforced plastic laminates. Bulletin of the JSME 25(202), 485. Zhu, H., Sankar, B.V., Marrey, R.V., 1998. Eval uation of failure criteria for fiber composites using finite element micromechanics. Jour nal of Composites Materials 32(8), 766. Zhu, T., Zhang, J.D., Atluri, S.N., 1998a. A lo cal boundary integral equa tion (LBIE) method in computational mechanics, and a meshless discretization approach. Computational Mechanics 21, 223. Zhu, T., Zhang, J.D., Atluri, S.N., 1998b. A meshless local boundary in tegral equation (LBIE) method for solving nonlinear problems. Computational Mechanics 22, 174. PAGE 137 137 BIOGRAPHICAL SKETCH Thi Dang was born and raised in Hai Duong, Viet Nam. He received his bachelors degree in mechanical engin eering from Hanoi University of Technology in 1993. After that, he joined PETROVIETNAM, where he spen t 4 years as a Weldi ng Engineer at Steel Structures Enterprise. He also got his masters degree in 1999 in computational mechanics from the University of Liege in Be lgium, and has been a lecturer at Hanoi University of Technology since 2001. His research interests focus on com putational mechanics with emphasis on micromechanical modeling of structural com posites, meshless methods and finite element models. After completion of the current body of work at the Un iversity of Florida in the area of micromechanics of composite structur es, Thi Dang will start working at General Motors Research and Development Center Michigan as a materials researcher. 