UFDC Home  UF Theses & Dissertations   Help 
Material Information
Thesis/Dissertation Information
Subjects
Notes
Record Information

Full Text 
PAGE 1 1 HOMOGENIZATION AND UNCERTAINTY ANALYSIS FOR FIBER REINFORCED COMPOSITE S By JINUK KIM A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 201 1 PAGE 2 2 201 1 Jinuk Kim PAGE 3 3 To my Lord, parents, brother and friends PAGE 4 4 ACKNOWLEDGMENTS I would like to sincerely thank Professor Nam Ho Kim, for his impeccable guidance, understanding and patience throughout this work. His mentorship and encouragement have provided me with great e xcitement through out the entire research life at the U niversity of Florida It was his direction, insight and criticism that lead me to the better way of research and study. I express my gratitude to the members of my supervisory committee, Professor Susa n Sinnott, Professor Anthony Brennan, Professor Simon Phillpot, Professor Elliot Douglas and Professor Youping Chen for their willingness to review my Ph.D. research and provide constructive comments to help me c omplete this dissertation. Special thanks to Professor Susan Sinnott for her kind and thoughtful support as an advisor in materials science and engineering department. I would also like to thank all my lab colleagues and friends in the materials sci ence and engineering department and the MDO lab in MAE department. Most importantly, I express my deep gratitude and thanks to my loving family. My father and mother have supported me throughout the school year with devotion and love. My brother has continuously pray ed for me and showed me the right way to live by his every life Finally, and most of all, I thank my Father in heaven for supporting my dreams and ambitions and forgiving me wrongdoings Thank you and love you. You gave me all. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ................................ ................................ ................................ .. 4 LIST OF TABLES ................................ ................................ ................................ ............ 7 LIST OF FIGURES ................................ ................................ ................................ .......... 9 ABSTRACT ................................ ................................ ................................ ................... 11 CHAPTER 1 INTRODUCTION ................................ ................................ ................................ .... 13 Composite Analysis ................................ ................................ ................................ 13 Motivation ................................ ................................ ................................ ............... 17 Objective ................................ ................................ ................................ ................. 18 2 LITERATURE REVIEW ................................ ................................ .......................... 20 Homogenization ................................ ................................ ................................ ...... 20 Reuss and Voigt Methods ................................ ................................ ................ 2 1 Eshelby Method ................................ ................................ ................................ 23 Mori Tanaka Method ................................ ................................ ........................ 26 Representative Volume Element (RVE) ................................ ........................... 27 Homogenization for Inelastic Behavior ................................ ............................. 28 Uncertainty Analysis ................................ ................................ ............................... 29 3 MATERIAL MODEL ................................ ................................ ................................ 32 Introduction ................................ ................................ ................................ ............. 32 Fiber P hase ................................ ................................ ................................ ...... 32 Matrix P hase ................................ ................................ ................................ .... 33 Elastic Material Model for Fiber ................................ ................................ .............. 35 Elastoplastic Material Model for Matrix ................................ ................................ ... 37 4 HOMOGENIZATION USING RVE ................................ ................................ .......... 43 Introduction ................................ ................................ ................................ ............. 43 Representative Volume Element (RVE) ................................ ................................ .. 46 Boundary Condition ................................ ................................ ................................ 47 Periodic Boundary Condition ................................ ................................ ............ 48 Application of UVE in Structure ................................ ................................ ........ 51 Homogenization for Elastic Behavior ................................ ................................ ...... 53 Effect of the Volume Fraction ................................ ................................ ........... 59 Effect of the Fiber Crack ................................ ................................ ................... 60 PAGE 6 6 Homogenization for Inelastic Behavior ................................ ................................ .... 62 Three dimension Ramberg Osgood model ................................ ...................... 65 Parameters Calculation ................................ ................................ .................... 68 Combined Loading Case ................................ ................................ .................. 72 5 UNCERTAINTY ANALYSIS ................................ ................................ .................... 76 Introduction ................................ ................................ ................................ ............. 76 Stochastic Response Surface Method ................................ ................................ .... 77 Collocation Points ................................ ................................ ............................. 80 Parametric Design ................................ ................................ ............................ 82 Uncertainty P ropagation in E lastic B ehavior ................................ ........................... 83 Input Variables for FEM Analysis ................................ ................................ ..... 84 FE Analysis of the RVE ................................ ................................ .................... 85 Monte Carlo S imulation ................................ ................................ .................... 87 Uncertainty P ropagation in I nelastic B ehavior ................................ ........................ 89 Input Variables for FEM Analysis ................................ ................................ ..... 89 FE analysis of the RVE in Elastoplastic Behavior ................................ ............. 91 Monte Carlo Simulation ................................ ................................ .................... 92 6 SUMMARY AND CONCLUSIONS ................................ ................................ .......... 94 APPENDIX: PYTHON CODE FOR PARAMETRIC ANALYSIS ................................ .... 97 LIST OF REFERENCES ................................ ................................ ............................. 100 BIO GRAPHICAL SKETCH ................................ ................................ .......................... 108 PAGE 7 7 LIST OF TABLES Table page 3 1 Characteristics of several fiber reinforced materials ................................ ........... 34 3 2 Typical matrix property ................................ ................................ ....................... 34 4 1 Mechanical property of constituents ................................ ................................ ... 53 4 2 Effective elastic properties ................................ ................................ .................. 57 4 3 Effective stiffness matrix from Eshelby method ................................ .................. 58 4 4 Effective Stiffness matrix from Mori Tanaka method ................................ .......... 58 4 5 Young's modulus compared to upper bound (vf=50%) ................................ ....... 58 4 6 Young's modulus compared to lower bound ................................ ....................... 58 4 7 Effective Stiffness from ABAQUS (Vf= 5%) ................................ ........................ 59 4 8 Effective stiffness from Eshelby method (Vf= 5%) ................................ .............. 59 4 9 Effective Stiffness from Mori Tanaka method (Vf= 5%) ................................ ...... 59 4 10 Young's modulus compared to upper bound ................................ ...................... 60 4 11 Young's modulus compared to lower bound ................................ ....................... 60 4 12 Degradation by transverse direction fiber crack having 50% crack area ............ 61 4 13 Effective stress in UVE according to the position ................................ ............... 63 4 14 Optimized Ramberg Osgood parameters for given matrix materials .................. 69 5 1 The third order H ermite polynomial basis for five input variabl es ( P represents the order of polynomial ) ................................ ................................ .... 82 5 2 Input values of parameters for FE anlysis ................................ .......................... 86 5 3 Coefficients of polynomial expansion for stiffness matrix components ............... 87 5 4 Output statistics of polynomial expansion ................................ ........................... 88 5 5 Input values of parameters for elastoplastic behavior ................................ ......... 91 5 6 The coefficients of polynomial expansion for strain energy in x direction ........... 92 PAGE 8 8 5 7 The propagation of the uncertainty fo r elastoplastic behavior ............................ 93 PAGE 9 9 LIST OF FIGURES Figure page 1 1 Intermediate material forms for thermoplastic composites. ................................ 14 1 2 Stamp forming of thermoplastic composites. ................................ ...................... 15 2 1 Simple uniaxial long fiber composite s ................................ ................................ 22 2 2 Sc hematic illustration of homogenization by Eshelby method ........................... 24 2 3 Different periodical elements for hexagonal packing ................................ ......... 27 2 4 Schematic of steps involved in MCS and SRSM ................................ ............... 30 3 1 Tensile stress strain curves of pitch based carbon fibers (HM50, HM70) and PAN based (3C, 5C) ................................ ................................ .......................... 35 3 2 Stress strain curves of epoxy L135i measured at different temperatures .......... 38 3 3 Elastic material for fiber and elastoplastic material for matrix. ............................ 38 3 4 Example of data input for material model in Abaqus ................................ ........ 39 4 1 Unidirectional fiber composite s composed of nine RVEs ................................ .. 43 4 2 RVE used in this research for uniaxial fiber composite s ................................ .... 44 4 3 von Mises stress A ) not constrained B ) constrained to be flat ......................... 47 4 4 Volume element deformed without periodic boundary condition ........................ 48 4 5 Periodic displacement boundary condition in x direction ................................ ... 49 4 6 X direction stress in the RVE with periodic boundary condition and UVE .......... 52 4 7 The unit volume element UVE in the composite structure ................................ 52 4 8 Transverse direction crack in fiber ................................ ................................ ..... 61 4 9 Heterogeneous structure composed of UVE piled up A ) Y direction B ) Z direction ................................ ................................ ................................ ............. 63 4 1 0 Structure to calculate inelastic properties in UVE A ) with a symmetric boundary condition. B ) Mises stress distribution in the structure under x direction tension. ................................ ................................ ................................ 64 PAGE 10 10 4 11 X direction average stress in homogenized material and in heterogeneous material A) for the first yield stress of 39.75MPa. B) for 50.25MPa. .................. 70 4 12 XY shear direction average stress in homo genized material and in heterogeneous material varying the first yield stress ................................ ......... 70 4 13 XY shear direction average stress in homogenized material and in heterogeneous material varying second yield stress ................................ ......... 71 4 14 Z direction average stress in homogenized material and in heterogeneous material ................................ ................................ ................................ .............. 72 4 15 Deformed shape and x direct ion stress under shear and tensile combined loading A) Heterogeneous materials B) Homogeneous materials .................... 73 4 16 Inelastic behavior comparison for heterogeneous and homogenized model A ) X direc tion average stress under combined load B ) XY shear average stress ................................ ................................ ................................ ................. 73 4 17 Contour of x direction stress under biaxial combined loading A) Heterogeneous materials. B) Homogeneous materials ................................ ..... 74 4 18 Average stress comparison under biaxial loading. A ) x direction B ) y direction. ................................ ................................ ................................ ............. 75 5 1 Examples of RVE generated by parameteriz ed design ................................ ..... 83 5 2 Histogram of stiffness matrix component A) C11. B)C12. C)C13. D)C33. E)C44. F)C55. ................................ ................................ ................................ .... 88 5 3 Elastoplastic behavior of modulus of fiber and matrix and volume fraction of fiber A) Ef1=197.85GPa Em1=3.09GPa B) Ef5= 250.15GPa, Em1=3.09GPa C) Ef 3 = 224.00 GPa, Em 1 =3. 50 GPa D) Ef1=197.85GPa, Em 5 =3.91GPa ................................ ......... 90 5 4 Z direction behavior of heterogeneous RVE A) Ef1=197.85GPa, Em1=3.09GPa B) Ef5= 250.15GPa, Em1=3.09GPa ................................ ........ 91 5 5 Histogram of strain energy for each loadin g condition in elastoplastic behavior A) x direction under x tension. B) xy shear direction under xy shear load. C) x direction under combined load. D) xy direction under combined load. ................................ ................................ ................................ .................... 93 PAGE 11 11 Abstract of Dissertation Present ed to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy HOMOGENIZATION AND UNCERTAINTY ANALYSIS FOR FIBER REINFORCED COMPOSITE S By Jinuk Kim August 201 1 Cha ir: Susan Sinnott Cochair: Nam Ho Kim Major: Materials Science and Engineering Because of the geometrical complexity and multiple material constituents, the behavior of fiber reinforced composite s is nonlinear and difficult to model These complex and non linear behavior makes the computation al cost much higher than that of general homogeneous materials. To make use of computational advantage h om ogenization method is applied to the fiber reinforced composite mode l to minimize the cost at the expense of det ail of local micro scale phenomena It is nearly impossible to make a homogenized model that behaves exactly the same as the composites in every aspect However, it would be worth while to use it when overall or specific macroscopic behavior is a major conc ern, such as overall heat flux in a given area, overall conductivity or overall stress level in a beam instead of microscopic accuracy In this research, the effective stiffness matrix is calculated using homogenization of elastic behavior In this case t he homogenization method was applied to the representative volume element (RVE) to represent elastic behavior. F or homogenization of inelastic material, the anisotropic Ramberg Osgood model is applied and its parameters are PAGE 12 12 calculated using the effective s tress and effective strain relation obtained from the heterogeneous material analysis. Many conventional structural analyse s ha ve been carried out on the basis of constant value s for mechanical properties, including h eat ca pacity and so on. It means uncertainty in the material parameters is ruled out and while it is widely accepted there is uncertaint y and it is inevitable in the procedure of obtaining the values and even during the design process This research applied the uncertainty analysis technique which make s use of a statistical approach such as stochastic response surface method (SRSM) to the behavior of the composite material The main purpose of appl ying the uncertainty analysis is to see how the uncertainty in m echanical propert ies propagate s to the macro behavior of the entire composites Uncertainty in some properties could vanish away to the macro behavior and others can result in amplification or decrease of the uncertainty in macro behavior. Utilizing this a nalysis, design parameters that are important can be identified which can help make an e ffective approach in development, design or manufacturing processes PAGE 13 13 CHAPTER 1 INTRODUCTION Composite Analysis Composite materials are macroscopically composed of more than two materials. A narrow definition of composites is restricted to combinations of materials that contain high strength fiber reinforcement and matrix that supports fiber [1] As science and technology advance, the demand of high performance mater ials has been increased and its s uperior properties over single component system s have made the application of fiber reinforced composite s popular in the industry. For structural applications, the greatest advantage of composite s may come from high stiffne ss and strength per weight. I nexpensive fiberglass reinforced plastic composites ha ve been used in various industrial and consumer products from automotive and aircraft since the 1950 s [ 2 ]. C ommercial airplane companies also provide a market for advanced composites application The Boeing 787 makes great use of composite materials in its airframe and primary structure than any previous commercial airplane s When the cost is not the main issue, such as military systems advanced composites have been used m ore extensively. Not only jet fighter s but also military vehicles feature integrated composite armored body. Application of advanced composite materials for structures has continued to increase, but one of the biggest factors that limit the widespread appl ication is their high cost in manufacturing process compared to conventional metallic alloy systems [3]. The high cost of composites manufacturing partially comes from the trial and error approach in design, process development and labor intensive manufact uring process. It is due in PAGE 14 14 part to the complex geometry and mechanical and thermal characteristics of composites. T here are many different manufacturing processes for composite s Because of its relatively high ductility, thermoplastic is one of the famou s materials used for the matrix part. Thermoplastic composites are supplied in a variety of ready to use intermediate forms which is generally called prepreg, as shown in Figure 1 1 T hey can easily be controlled and processed th rough the standard production techniques [4]. Figure 1 1 Intermediate material forms for thermoplastic composites. When the behavior of composite structure is expected to be in elastic in working condition such as automo tive or aircraft working condition plasticity could not be the main area of interest. But, it becomes a main issue in manufacturing process where composite materials experience plastic deformation. The product processes that are currently being used in in dustry are basically adapted from sheet metal forming processes. Pr oduction techniques for thermoplastic such as roll forming pultrusion, compression molding, diaphragm forming, and stamp forming are similar to general metal forming processes. Figure 1 2 shows a schematic stamping forming process for the intermediate composite form. Since the purpose of this process is to produce permanent deformation to the material, it is called plasticity in the classical mechanics. PAGE 15 15 Figure 1 2 Stamp forming of thermoplastic composites. Because of the high complexity of both material behavior and structural geometry, it is difficult to predict the production process. Therefore it is essential to utilize numerica l approach es to understand and predict the behavior of heterogeneous materials. The computational model and the prediction of the process require understanding s of stress strain relations of the heterogeneous materials However, if heterogeneous materials are modeled directly using a full scale simulation, they requires too much computational source because they are usually very small features compared to the size of structures The c lassical plasticity theory is the most popular approach to model inelasti c behavior and widely being used in many applications But it requires too many parameters and complexity T his classical plasticity theory is reasonable for homogeneous or isotropic materials but not for heterogeneous, anisotropic yielding, and anisotrop ic hardening materials at least It has to contain more parameters and complex analytical expressions and more iteration to solve nonlinear equations leading to a g reat amount of computational cost and still having questions on accuracy. PAGE 16 16 T he purpose of co mputational model is to predict the response of a structure accurately, which requires using proper model parameters However, choosing model parameters for composites application have relied In addition experiment s show ed that most structural parameter s ha ve random, stochastic characteristic s When the effect of these uncertainties in parameters is understood, the analysis and design process to develop advanced composite materials can have benefit from it by giving more accurate and reliable results. If c ritical parameters can be identified, detail ed or optimized design can be possible by paying more attention on th at critical parameter s effect The process of estimating the effect of model parameter uncertainty on the r esponse uncertainty is called uncertainty propagation. Accordingly many probabilistic approaches have been introduced to predict more accurate behavior of system under uncertainty, and they showed the random character play a very important role in the dec ision making process I t seems that the random characteristic of heterogeneous materials would be more important than that of homogeneous materials For statistical approach es to engineering problem s t he most popular method for uncertainty propagation an alysis is Monte Carlo Simulation (MCS) T h e MCS method is one of the sampling methods which generates random samples according to the probability distribution of random parameters, and using physical model to calculate the random samples of response. This method has an advantage of simplicity, which can be applicable to virtually all engineering applications, and the accuracy of the method is independent of the number of model parameters. However, in order to estimate the distribution of response accuratel y, a large number of samples are required, which is PAGE 17 17 equivalent to a numerous evaluation of the physical model. For modern advanced computational models, this can be a significant bottleneck as they requires significant amount of computational resources. Mo tivation For simulation of composite structures, using homogenization and a simpler material model than the classical plasticity theory can reduce computational cost significantly. I t is virtually impossible to make a homogeneous material that behaves exac tly same as heterogeneous materials in every aspect and scale However when a specific behavior is of interest, a homogenized material can be modeled so that it can have similar behavior under a limited condition For example, when the temperature on the surface of heat barrier that is made of composite structure is the main interest, a homogeneous material that ha s the same or similar thermal property can replace the heterogeneous part. The s ame analogy can be applied to the constitutive model that define s the relation between stress and strain. Although homogenization of elastic properties ha s been studied intensively by many researchers, not much research has been reported for homogenization of inelastic properties. Most of inelastic homogenization rese arch is based on the classical plasticity theory which requires many parameters to be defined and consequently leaving assumptions that cannot be used for heterogeneous materials, e.g., isotropic hardening. Therefore more effective homogenization model than the classical one needs to be developed to make a good use of computational tools. Although MCS is widely used it becomes impractical for computationally expensive models because it requires a large number of samples and each sample means evaluatin g the expensive computational model Performing 50,000 finite element PAGE 18 18 analyses using input values generated by MCS for composite material is far beyond current computational ability especially when nonlinear behavior is considered. Although t here are many approximation methods to improve computational cost, all of them introduce some kind of approximation, such as linearization. Therefore, it would be desirable to make the MCS sampling computational ly inexpensive. T he complexity in modeling inelastic behavi or using plasticity theory prevents not only the homogenization process, but also the statistical study because the statistical application still requires a significant amount of simulations In addition when geometric parameters are included, such as the size of fiber, it would be difficult to take into account their effect using the classical plasticity model. Because of these reasons, uncertainty analysis on inelastic behavior has difficulty in application for heterogeneous materials. Objective The fir st objective of this research is to develop a homogeniz ation methodology for unidirectional fiber reinforced composites by using numerical method so that the homogenized material model can provide the same effective stress strain relation in macroscopic sc ale as heterogeneous materials. For the elastic behavior, effective are to be calculated. For inelastic behavior, the composites is considered as elastoplastic materials and is to be homogenized through the Ramberg Osgood model with anisotropy tensor imbedded in, which reduces significantly the number of plasticity parameters and numerical iterations. Throughout this research, representative volume element (RVE) and unit volume element (UVE) are to be used to understand the behavior of composites and homogenize the heterogeneous materials more accurately. PAGE 19 19 The periodic boundary condition that assigns the periodicity to RVE model is also discussed for its validity. The secon d objective of this research is to perform uncertainty analysis using stochastic response surface method ( SRSM ) that uses more efficient sampling method than Monte Carlo Simulation. The u ncertainty analysis is performed not only for elastic properties but also for in elastic properties. The h omogenized material with a simpler material model enables the uncertainty analysis for inelastic behavior. From the uncertainty analysis, the focus is given on how much uncertainties that exist in mechanical or geometric al properties propagate in to the macro scopic behavior of the composites Consequently, it gives information on which parameters are critical in a given specific conditions so that composites development, design and manufacturing process can have benefits f rom the analysis. PAGE 20 20 CHAPTER 2 LITERATURE REVIEW Homogenization Almost all structural material s are heterogeneous in micro and /or macro scale I t is generally known that the macro scale behavior of a structure is caused by the micro scal e behavior. T he micro scale behavior is often more complicated than that of the macro scale. In fiber reinforced composite s for example, the i nhomogeneous local displacement field can be develop ed even under uniform global displacement. In such a case, it would be necessary to have a micro scale model in order to describe this local behavior. However, b ecause of limited computational resources and modeling complexity it would be impractical to perform all in one analysis for composite structure from the m icro scale to the macro scale [ 7 ]. Instead of full scale analysis when effective properties of the composite materials are used by averaging the local behavior of individual fibers and matrices [8] the composite structure can be modeled as a homogeneous material in the macro scale considering only effective properties. Among numerous methods to predict the effective properties of composite materials, a variety of theories have been developed for homogenization such as effective medium model s of Eshelby [ 9 ], Hashin [ 10 ], and Mori and Tanaka [ 11 ] T he simple bounds on the effective moduli can be determined by the approaches of Voigt [ 12 ] and Reuss [ 13 ]. From a different perspective, Hill [14] and Christensen [15] proposed a self consistent method. Hill [ 1 6 ] also showed that Voigt and Reuss approach es provide rigorous upper and lower bound s, respectively M athematical or analytical homogenization methods had been pioneered by Bensoussan [ 17 ] and Sanchez Palencia [ 18 ]. The active c omputational aspects of ho mogenization have been PAGE 21 21 initiated by Guedes and Kikuchi [19 ] Over the past decade major contributions have been made to extending the theory of computational homogenization to nonlinear regime [ 20 23 ] and to improving fidelity and computational efficiency of numerical simulations [ 24 34 ]. The fundamental theory that helps to understand the homogenization concept will be shown in C hapter 4 Reuss and Voigt Method s The simplest, but not necessarily the best, bounds of homogenization are the Voigt (1 889 ) and Reuss (1929) bounds. The Voigt upper bound on the effective elastic modulus, of a mixture of N material phases is ( 1 ) where is the volume fraction of the i th constituen t and t he elastic modulus of the i th constituent. There is no way that nature can put together a mixture of constituents ( e.g. a rock) that is elastically stiffer than the simple arithmetic average of the individual constituent moduli given by the Voigt bound. The Voigt bound is sometimes called the iso strain average, because it gives the ratio of average stress to average strain when all constituents are assumed to have the same strain. The Reuss lower bound of the effective elastic modulus, M R is ( 2 ) There is no way that nature can put together a mixture of constituents that is elastically softer tha n this harmonic average of moduli given by the Reuss bound. The PAGE 22 22 Reuss bound is sometimes called the iso stress average, because it gives the ratio of average stress to average strain when all constituents are assumed to have the same stress. Mathematically the M in the Voigt and Reuss formulas can represent any modulus, the bulk modulus K the shear modulus The simplest example is the elastic behavior of aligned long fiber composites as shown in Figure 2 1 Figure 2 1 Simple uniaxial long fiber composite s The simplest way to homogenize the elastic behavior is to consider it as a parallel slabs bonded together. Both constituents are under the same strain and this condition is valid for loading along the fiber axis. ( 3 ) For composites having the reinforcement is subject to much higher stress than the matrix and there is a redistri bution of the load. The overall stress can be expressed in terms of the two contributions. Here f represents the volume fraction of fiber. ( 4 ) The elastic modulus of the composites can be written PAGE 23 23 ( 5 ) This can be simplified using equation ( 3 ) ( 6 ) This famous Rule of Mixture shows that the global stiffness is simply a weighted average between moduli of the two constituents, depending only on the volume fraction of fibers. The equal strain model is often called as Voigt model while the method that uses an equal stress condition is called Reuss model [ 35 ] act ing on the matrix when the transverse loading is applied. It is shown in equation ( 7 ) ( 7 ) Next, the net strain is the sum of the contributions from the matrix and the fiber as ( 8 ) Then, elastic modulus of the composites can be written as equation ( 9 ) from equation ( 7 ) and ( 8 ) ( 9 ) It is also called inverse rule of mixture Eshelby Method fibers) with one made of the matrix material. This is called the equivalent homogeneous PAGE 24 24 inclusion [9]. This equivalent inclusion is assumed to have an appropri ate strain, called the equivalent transform ation strain such that the stress field is same with the actual inclusion. The following is a summary of the steps followed in the homogenization procedure according to the Eshelby method which is shown in Figure 2 2 Figure 2 2 Schematic illustration of homogenization by Eshelby method Consider an initially unstressed elastic homogeneous material Imagine cutting an ellipsoidal region (i. e. inclusion) from this material. Imagine also that the inclusion undergoes a shape change free from the constraining matrix by subjecting it to a transformation strain ( B ). Since the inclusion has now changed in shape, it cannot be replaced directly into the hole in the matrix material. Imagine applying surface tractions to the inclusion to return to its original shape, then imagine returning it back to the matri x material ( C ). Imagine welding the inclusion and matrix material together then removing the surface tractions. The matrix and inclusion will then reach an equilibrium state when the inclusion has a constraining strain relative to the initial shape before it was removed ( D ). The stress in the inclusion can now be calculated as follows assuming the strain is uniform within the inclusion: PAGE 25 25 ( 10 ) w here are the components of the elasticity tensor of the matrix material. Eshelby has shown that the constraining strain can be described in terms of the transformation stra in using the equations: ( 11 ) The Eshelby tensor S is a fourth atio of This Eshelby tensor or the concept of the Eshelby tensor is widely being used in studies on heterogeneous materials [36]. Finally, the stress in the inclusion is determined by substituting eq uation ( 11 ) into equation ( 10 ) and simplifying to obtain: ( 12 ) w here I is the components of the fourth rank identity tensor. Th e above equation can be rewritten in matr ix forms as follows: ( 13 ) The braces are used to indicate a vector, while the brackets are use d to indicate a matrix. Next, the expressions of the Eshelby tensor, S, are presented for the case of long infinite cylindrical fibers. In this case, the values of the Eshelby tensor depend on of the fibers and are determined as follows: PAGE 26 26 ( 14 ) unbounded, which is infinite, space. This also means that that this is based on the assumption that the constitue nts are not influenced by each other. So, this method is only applicable to very low volume fraction of heterogeneous materials [ 37 ] ; i.e., the inclusion is very small compared to the size of matrix and the location of inclusion would not make any differen ce concentration factors and self consistent theory [ 38 ] Mori Tanaka Method The main difference be tween Mori Tanaka and Voigt method s comes from the assumption on the Poisson ratio of fiber and matrix are different, while the latter assumes they are the same. Mori Tanaka method is almost same as Eshelby met hod but slightly different in taking account of effective strain. Mori Tanaka method also uses Eshelby tensor. Thus, this method also cannot consider the size and position of fiber. However, Mori Tanaka method can be applied to higher concentration of fi ber while Eshelby method only can be applied to dilute concentration (V f << 1). In addition to that, when the inclusions is considered as rigid particles, voids (V f =0) and matrix (V f =1), Mori Tanaka method works better than Eshelby method in those limit ca se. PAGE 27 27 Representative Volume Element (RVE) In most micromechanical analysis of fiber reinforced composite materials, representative unit cell or representative volume element, RVE, is the first step into the analysis. The reason why RVE concept is widely use d is the periodic characteristics of fiber composites [ 39 ] In addition to that, unit cell approach can save computational cost for the study on heterogeneous materials. Based on this advantage, it can be applied to characteriz ing heterogeneous materials w ith macroscopically and statistically homogeneous structure [ 40, 41]. For the shape of the RVE, square packing and hexagonal packing are most popular in RVE application for unidirectional fiber composites [ 42 49 ]. Several shapes of the unit cell that is possible for the unidirectional fiber composites as shown in Figure 2 4 [ 50 ] Figure 2 3 Different periodical elements for hexagonal packing To be verified as a n appropriate RVE, the volume to calculate average behavior need to meet two criteria [51]. First, it must be small enough with respect to the dimensions of the macro scale so that it can be considered as a material PAGE 28 28 point in the equivalent homogeneous continuum. Second, it must be large enough respect to the scale of the inclusion phase to have elastic properties independent of the loading condition [52]. Although the RVE concept is popular in micromechanical analysis, it seems that the effort to find a proper bou ndary condition was not given much in researches compared to its popularity It is not difficult to find a research that only considered the displacement boundary condition for composite analysis [53 55]. To represent the behavior of RVE proper ly it is kn own that not only displacement boundary condition but also traction boundary condition should be applied so that periodicity condition can be satisfied [56]. For this reason, the mixed periodic boundary condition has been the acceptable boundary condition for RVE research [57]. Homogenization for I nelastic B e havior R esearch on homogenization or effective properties for nonlinear behavior is still in progress Most research in this field still consider s isotropic hardening and even isotropic behavior [58 60 ]. It can be easily found that most of inelastic modeling is based on the classical plasticity theory [61 63]. However, it is obvious that the material responses are anisotropic due to the prescribed direction of fiber. Griffin and Kamat [64] used rthotropic yield criteria and flow theory along with unidirectional Ramberg Osgood plastic material model Similar to the fiber composites the homogenized materials were modeled as orthotropic elastoplastic one Although the material behavior in the princ ipal directions is well described by Ramberg Osgood relation [65] the isotropic plasticity equivalent (IPE) material concepts applied was not well matched to the plastic behavior of composite material, which is the focus on the proposed research. O ther re search es make use of multiaxial Ramberg Osgood PAGE 29 29 model for anisotropic material to describe the local inelastic behavior [66 67]. However, th ose model s are introduced for anisotropic homogeneous material s Uncertainty Analysis To utilize statistical approa ch, sufficient number of samples within the distribution of the input variables should be prepared. One of the m ost popular method s to generate random samples is Monte Carlo Simulation (MCS) I t is a numerical method based on random sampling and statistica l estimation [ 68 ] MCS has been applied to various engineering application of composites ; e.g., fatigue and failure modeling [ 69 ] and modeling of random properties [ 70 ]. Although it is widely used, the classical method s for uncertainty analysis such as st andard MCS require a large number of samples in order to estimate the distribution of model output accurately [ 71 ]. This is not appropriate for computationally expensive models because it takes too much time to repeat all the samples created by MCS. There have been several approaches to reduce the number of samples for MCS [ 72]. How ever these approaches require special knowledge on the response. Because of this, other RSM methods based on series approximation have been developed to reduce the required samp le numbers. One of the methods is Stochastic Response Surface Method (SRSM) [ 73 74 ] which is a customized response surface technique for random inputs A schematic of the steps involved in the uncertainty analysis methods using MCS and SRSM is compared i n Figure 2 4 [ 75 ] PAGE 30 30 Figure 2 4 Schematic of step s involved in MCS an d SRSM Another method to overcome the disadvantage of MCS is Latin Hypercube Sampling (LHS). LHS samples based on equally prob able intervals while MCS samples random ly The less independence of LHS on sampling number is the advantage over MCS. There are other types of uncertainty analysis methods, which are based on perturbation techniques [76, 77]. Perturbation based method is e ffective for a linear problem or a small variation of a stochastic variable but it may not effective for a nonlinear problem or a large variation of variables [78 79 ]. A variational method for the derivation of lower and upper bounds for the effective ela stic moduli has been studied [ 80 ]. There have been several studies on uncertainty of elastic mechanical property T he influence of uncertainty in microscopic properties on homogenized elastic property using the perturbation based analysis [81 83] and polyn omial chaos expansion was used for prediction of complex physical system that has several subsystems [84]. A pplication of uncertainty analysis on elastic behavior can be easily found for aerospace PAGE 31 31 application. It has complex phenomena and the uncertainty propagation into aerothermoelastic or aeroelastic behavior has been studied [85 8 7]. Microstructure effect on effective constitutive property was studied using statistical volume element [ 8 8 ] but it does not consider the anisotropic properties. The study o n the uncertainty propagation in multiaxial nonlinear materials behavior is very difficult to find. Some has been done in one dimensional linear hardening or isotropic materials [ 89 9 1 ] but it is not easy to find one associated with anisotropic composite m aterials PAGE 32 32 CHAPTER 3 MATERIAL MODEL Introduction Design goals of fiber reinforce d composites include high specific strength or stiffness on weight bases. Low density fiber and matrix materials enable high specific strength and moduli composites. The overall mechanical property of composites depends not only on the properties of consti tuents, but also on the interfacial characteristics. It also depends on the geometrical characteristics of fiber itself. For example, it depends on if fiber is continuous or discontinuous, and if fibers are aligned or randomly oriented There are too many things to consider every aspects of the composites In this research, for a clear understanding of the validity on approach and result for homogenization and uncertainty analysis, the parameters to be considered are narrowed down to simp le but most critical parameters. T he bonding region is not considered in modeling. It is one of the most important aspects in composite materials but its effect is too wide to keep track of it It is important but too much is unknown to be used in homogeni zation or uncertainty analysis in this research. Thus, material models for fiber and matrix are the only constituents considered in this research for heterogeneous materials. Fiber P hase Because of strong bonding in elements of low atomic number, such as C B, A l and Si, they can form stiff and low density materials. These materials can be made from the elements themselves or from their compounds or oxides such as Al 2 O 3 and SiO 2 which can be made to ceramic fibers. T he key point is that the flaw exists in these materials, especially the one open to surface, leading to fracture failure. Thus, only the PAGE 33 33 form of fiber with small radius can enables very high strength applications and this feature is advantage of fiber reinforcement. It is widely known that the smaller the fiber diameter and the shorter the length, the higher strength, but the greater variability described by Weibulll statistics [ 9 2 ]. G lass fibers are most famous composite materials because of relatively high strength at low cost [ 9 3 ]. High stre ngth glass fibers have been used in structural applications such as pressure vessel since the 1960s that does not require specific stiffness. The other two most popular types are carbon and aramid fibers, while polymer is usually an epoxy, vinylester or po lyester thermosetting plastic. Carbon fibers are widely used for aerospace applications but the drawback of the mechanical properties of carbon fibers is its low ductility compared to glass, SiO 2 and Kevlar fibers applications [ 9 4 ]. Table 3 1 [95] shows the several fiber materials generally used. Matrix P hase The matrix phase for fiber reinforced composites can be a metal, polymer or ceramic. Generally, matrix works as binding materials that supports and protects the fibers and metals and polymers are used because of ductility. As a supporting medium, externally applied stress is transmitted and distributed to the fibers and only a small portion is loaded in matrix phase. It also transfers loading when the fiber is br oken. For protection roles, matrix protects fibers from surface damage by chemical reaction or impacts, which can induce surface cracks inducing fracture. Table 3 2 [9 3 ] shows properties of some matrix materials PAGE 34 34 Table 3 1 Characteristics of several fiber reinforced materials Material Specific G ravity Tensile Strength (GPa) Specific S trength ( GPa) Modulus of Elasticity (GPa) Specific Modulus (GPa) Whiskers Graphite 2.2 20 .0 9.1 700 318 Silicon nitride 3.2 5 .0 7 .0 1. 6 2.2 350 380 109 118 Aluminum oxide 4.0 10 .0 20 .0 2.5 5.0 700 1500 175 375 Silicon carbide 3.2 20 .0 6. 3 480 150 Fibers Aluminum oxide 3.9 1. 4 0. 4 379 96 Aramid (Kevlar 49) 1.4 3.6 4.1 2.5 2. 9 131 91 Carbon 1. 8 2. 2 1.5 4.8 0.7 2.7 228 724 106 407 E glass 2. 6 3.5 1.3 72 28 Boron 2. 6 3.6 11.4 400 156 Silicon carbide 3.0 3.9 1.3 400 133 UHMWPE (Spectra900) 1.0 2.6 2. 7 117 121 Metallic Wires High streng th steel 7.9 2. 4 0.3 210 26 Molybdenum 10.2 2.2 0.2 324 31 Tungsten 19.3 2.9 0. 2 407 21 Table 3 2 Typical matrix property Material Density(Kg/m 3 ) Et (GPa) t (MPa) (10 6 / C) Polyester 1200 1400 2.5 4.0 45 90 0.37 0.40 100 200 Epoxy 1100 1350 3.0 5.5 40 100 0.38 0.40 45 65 PVC 1400 2.8 58 50 Nylon 1140 2.8 70 100 Polyethylene 960 1.2 32 120 Epoxy, which is widely used as a matrix material, is a thermosetting polymer [96 ] and has wide range of applications including fiber reinforced composites and general purpose adhesives. Thermosetting is one of the subdivisions of p olymer classification according to the behavior when temperature rises. Thermoplasts are the other division. Thermosetting polymers become permanently hard when heat is applied and do not soften when heat is applied subsequently [ 9 5]. Thermosets are basically brittle materials, PAGE 35 35 while thermoplastics can undergo more plas tic deformation. In the same t hermosets category, epoxy is tougher and more expensive than vinyl e s ters, ha s better resistance to moisture and hea t distortion and shrink s less than polyesters when curing [ 9 7 ] Because of its better mechanical property, it has been extensively used for aerospace applications. Elastic Material Model for Fiber Major type of commercial carbon fiber has tensile modulus in range from 150GPa to 380GPa [9 8 ] a nd most of them have linear elastic behavior as shown in Figure 3 1 [ 99 ]. Figure 3 1 Tensile stress strain curves of pitch based carbon fibers (HM50, HM70) and PAN based (3C, 5C) Therefore the fiber is assumed to behave linear in this research and the Young s modu lus for fiber is chosen as 224GPa and Poisson ratio of 0.3 based on the data in published literatures. A linear elastic material model is valid for small elastic strains, normally less than 5% [10 0 ] and experimental data in Figure 3 1 shows that the strain range that fiber is in application stays less than 5%. Therefore it is reasonable to assume fiber is linear elastic material. PAGE 36 36 stress and total el astic strain as ( 1 ) The represents t he total true stress tensor, is the fourth order elasticity tensor, and is t he total elastic strain tensor. The simplest form of elasticity is isotropic material, which is considered as the material for fiber in this research and its stress strain relation is given matrix form in equation ( 2 ) ( 2 ) In equation (2), G is the shear and These input parameters are a function of temperature but the temperature effect is not considered in this research. One thing to be cautious i s that Abaqus use engineering strain for shear strain. However, the unidirectional fiber composites is an orthotropic material and the elastic compliance is defined as PAGE 37 37 ( 3 ) The physical in terpretation of is that it is the transverse strain in the j direction when the i is generally not same as When , and then it is called transversely isotropic material, which is a subclass of orthotropic material. The RVE used in this research is orthotropic material although the form of stiffness matrix seems alike transversely isotropic because 12 plane th at is perpendicular to 3 direction is not actually an isotropic plane that has no directionality on that plane. Elastoplastic Material Model for Matrix M atrix is generally made of polymer especially in aircraft application s Although the characteristic of the polymer is viscoelastic or viscoplastic behavior affected by time, temperature and moisture the isotropic homogeneous elastoplastic material behavior is considered in this research for the simplicity. Figure 3 2 [10 1 ] shows a stress strain curve for an epoxy material at different temperature. Although the temperature effect is not considered in this research, the matrix behavior is based on the data at 25 o C. PAGE 38 38 Figure 3 2 Stress strain curves of e poxy L135i measured at different temperatures Based on these data, the elastoplastic material behavior is implemented into Abaqus finite element ( FE ) analysis. The schematic drawing of the constit utive behavior of these two material models are shown in Figure 3 3 Figure 3 3 Elastic material for fiber and e lasto p lastic material for matrix. T he elastoplastic material is modeled with two yield point. One is for the end point of linear elastic regi on, which is 45MPa and the other is at 75MPa that begins perfec t plastic deformation. The plastic strain when the perfect plastic region begins is given as 0.015 based on the stress strain curve in Figure 3 2 For the plasticity m odel implemented in the library in Abaqus FE program, it approximates the smooth stress strain behavior of the material with piecewise linear PAGE 39 39 curves joining the given data points. Any number of points can be used to approximate the actual material behavior which means it can be used as a very close approximation of the actual material behavior as shown Figure 3 4 Figure 3 4 Example of data input for material model in Abaqus When the stress s train relation is given by strain dependent properties instead of constants, it behaves nonlinear ly Computational challenges come from the fact that the equilibrium equations should be written using material properties that depends on the strains, which a re not known in advance. In order to describe the plastic stresses strains relation, conventional plasticity theory uses yield criterion that define the condition of the onset of inelastic behavior. This yield condition is generally defined by a yield func tion described by stress, and hardening equivalent stress, R. When R=0, yield occurs when equivalent stress value reaches the initial yield stress, This hardening stress can be a constant for linear strain hardening or a function of plastic strain for isotropic hardening and a function of both stress and strain for kinematic hardening. So, the yield occurs when the equation ( 4 ) is sati sfied. ( 4 ) Equation ( 4 ) is called yield surface and it is worth to note that is a scalar equation. PAGE 40 40 The material is in elastic region when f < 0 and in inelastic region otherwise. However, f > 0 cannot be occur physically and f = 0 whenever it is in elastic region. It is c alled consistency condition. Since the stress must be on the yield surface and the size of the yield surface is related to the magnitude of the accumulated plastic strain, the magnitude of the plastic strain increment also must be related to the stress inc rement. This relation produced plastic flow rule defined as equation ( 5 ) [10 2 ]. ( 5 ) These yield function and flow rule are the fundamental of classical plasticity theory. Different yield criteria have their own definition for the terms used in equation ( 4 ) For example, isotropic von Mises yield criteria [10 3 ] the equivalent stress is defined as equation ( 6 ) and denotes deviatoric stress. ( 6 ) Beyond yield point, stresses are related to strains by incremental constitutive plastic strain, the constit utive relation becomes as ( 7 ) where i s the elastoplastic tangent operator that relate the stress and total strain. From equation ( 5 ) and ( 7 ) we can see that tangent operator is a function of which is called plastic consistency parameter or plastic multiplier. PAGE 41 41 Since the yield surface should be updated whenever the loading is applied beyond yield point, the yield surface also should be updated using ac cumulated plastic strain at that point. T he updated yield function can be written as ( 8 ) Without expanding the equation ( 8 ) in detail, we can see that it is a nonlinear equation in terms of This equation can be solved to compute using the local Newton Raphson method. If isotropic/kinematic hardening i s a linear function of or the effective plastic strain, then only one iteration is required to calculate the every components updated points. After is found, stress, effective strain and hardening parameters can be obtained. As has been shown a bove, the material constituent of elastoplastic material can be calculated algorithmically. E quation ( 7 ) relates infinitesimal stress increments with correspo nding infinitesimal strain increment at a given state. However, the iteration for equilibrium is carried out in a finite magnitude rather than infinitesimal one Therefore Abaqus carr ies out the global iteration to find structural equilibrium, while the u ser subroutine (UMAT) that models material behavior requires local iteration to calculate the incremental constitutive model. These are the core of elastoplastic materials FE analysis and also are the main reasons it cost a lot of time especially when it h as more than two materials. Another issue when approached using conventional plasticity theory represented by equation ( 4 ) is that the yield point is decided by equivalent stress which is a scalar even if three dimensional stress and strain components are involved. To compensate for this issue when it is applied to anisotropic materials, Hill introduced anisotropic yield PAGE 42 42 criteria [10 4 ] and modified versions h ave introduced since then [10 5 10 7 ]. modified versions of anisotropy yield criteria are to define anisotropic yield point which means that the direction dependency is considered However these yielding criteria have nothing to do with harden ing which describe the behavior of material after yielding Consequently, anisotropic hardening should be considered which is not easy as long as equivalent stress is used. The homogenization model in this research can reduce the local iterations by appl ying tensor form of Ramberg Osgood model and application of anisotropy tensor solves the anisotropy issues. Ramberg and Osgood introduced simple formula to describe stress strain curve in terms of three parameters in 1943 [ 10 8 ] as equation ( 9 ) ( 9 ) where s is stress, e is strain and K and n are constants. In contrast to the classical plasticity model, Ramberg Osgood model describe entire materials behavior in an equation without defining any yield point and it is adequate to represent stress strain curve that does not ha ve sharp yield point. It was originally introduced for uniaxial loading condition and the modified form for the multiaxial loading condition is used in this research. More d etails of multiaxial Ramberg Osgood model are discussed in Chapter 4. PAGE 43 43 CHAPTER 4 HOMOGENIZATION USING RVE Introduction In a small scale, all materials are heterogeneous. To understand and be convinced of all the mechanism at a high degree of accuracy, one should investigate all phenomena that occur at atom ic or molecular scale. When en gineering materials were to be designed at this level of accuracy, the required amount of computational re source s would be out of practical ity Figure 4 1 shows that only a portion of the unidirectional fiber composites and its analysis to calculate stress and strain components under uniaxial static tension loading took more than 20 hours based on 64bit Pentium 3.4GHz processor. Figure 4 1 Unidirectional fiber composite s composed of nine RVE s In this chapter homogenization of composite material for elastic behavior is performed on a representative volume element, RVE. For a brief comparison, it took 17 minutes for a simple tension loading. Heterogeneous materials can be geometrically represe nted by the concept of periodicity as shown in Figure 4 2 and t he uniaxial fiber PAGE 44 44 composites can be regarded as the periodic materials composed of blocks of representative volume element RVE. Figure 4 2 RVE used in this research for uniaxial fiber composite s Performing an analysis on the RVE, not on entire engineering structure, one can reduce the computational time to understand the mechanical behavior of composite material B ut, it requires app ropriate boundary condition so that the RVE can be analyzed as a component in the middle of the structure, not as an independent separate one. This is why the boundary condition for the RVE is important. The different results from different boundary condi tion are shown later. In this research, t he periodic boundary condition contains both periodic displacement boundary condition and periodic traction boundary condition, which also called as mixed periodic boundary condition. T he important issues on this pe riodic boundary condition are also discussed later A d ifferent approach from the RVE concept used for elastic application will be applied for the homogenization for inelastic behavior. However, when the simulation of entire structure is required, not micr o mechanical study, even this heterogeneous RVE cannot give any cost saving because it ha s to be assembled again to compose entire structure as heterogeneous. So, heterogeneous RVE can save cost for the micro scale PAGE 45 45 analysis alone but not for the macro scal e analysis or simulation such as manufacturing process One of the hypotheses to overcome this difficulty that arises for larger scale analysis is that the structure of the material is considered as a continuum This means that there exist measures associa ted with properties that govern the behavior of the media and t he properties of material at a point can be computed using an averaging scheme T hese properties are actually similar to the averages of very complicated interactions and phenomena in the atom ic scale. Like wise h omogenization is similar to the continuum concept in terms of the homogenized medium that has properties governing the behavior of the heterogeneous media. I n this research, the mean field theory is applied which is also known as the average field theory. Using this theory, macro field is defined as the volume average of corresponding micro fields. T his average scheme was introduced for analytical method s in earlier researches ; i.e. Reuss and Voigt method as discussed in the C hapter 2 The main issue of the analytical method is that it is valid only under a specific condition. For example, whenever the shape or the position of fiber change s the analytical expression for averaged propert ies should be derived again. Furthermore, it is v ery difficult to formulate all complex geometry effects and nonlinear behaviors. On the contrary, numerical method does not have th ese limitation s As the computational power continues to grow as technology advances, it enables calculation of more complica ted geometry and reduces the number of assumptions in the model For this reason, the numerical method is used in this research using a commercial program Abaqus. Python script is written in order t o calculate effective properties as a post process from n umerical model and to generate PAGE 46 46 parameterized heterogeneous RVEs. As the Abaqus has strength in capability of user customized material, anisotropic hardening Ramberg Osgood model was implemented using a FORTAN code for the user subroutine (UMAT) At this point, it needs to be mentioned that the RVE approach gives approximate macro scale understanding of the composites because there is no exact periodicity in real random media. It could get closer to the more accurate analysis results when the real com posites has higher periodicity. Representative Volume Element (RVE) The heterogeneous RVE used in this research has been built using a commercial FE program Abaqus I n this research, three dimension square unit cell is considered as a RVE because it is re latively easy to impose symmetry and periodic boundary conditions compared to the hexagonal one as discussed in Chapter 2 briefly. RVE is a unit size cub e having dimensions of a=1, b=1 and c=1 as shown in Figure 4 2 A total of 1 1781 C3D8 elements are used to model the RVE T he stress and strains are calculated at the eight integration points in each element and the output data from entire RVE are post processed to calculate required values for statistical analysis using Python a nd Matlab. The RVE part has mapped mesh on every outer surface so that constraint equations between opposite facing nodes can be constructed. How the constraint equation and periodic boundary condition are applied is discussed la ter The RVE has two mater ials. One is an isotropic homogeneous uniaxial fiber parallel to z axis and the other is an isotropic homogeneous matrix. T he materials for fiber and matrix are described detail in Chapter 3. The radius and the center of the fiber are parameterized as vari ables for statistical analysis. A perfect bonding between the PAGE 47 47 fiber and matrix is assumed for simplicity because the delamination or boundary effect is not covered in this research. Boundary Condition As the RVE is not distinguishable from the next in the periodic structure it can be said that the response of the entire composite structure under uniform macroscopic loading is same as the response of the RVE under the same loading condition. To apply the RVE concept to the periodic composite materials, the appropriate periodic boundary condition is an important part of modeling The Figure 4 3 shows that how the distributions of von Mises stress is different in the heterogeneous RVE according to different boundary conditions on the surface W hen the RVE does not have any constraints and can deform freely, the deformed shape has curved shape under normal tension and under shear tension A B Figure 4 3 von Mises stress A ) not constrained B ) constrained to be flat PAGE 48 48 W hen the boundary of the RVE has constrained not to move in transverse direction, it results in local or overall higher stress When the shear is applied on the RVE whose transverse direction surface to be flat, it has high stress induced inside. However, most stud ies on the effect of the boundary condition w ere performed for the elastic behavior or elastic mechanical property. In the process of the homogenization study in this research, the effect of boundary condition in inelast ic behavior is also discussed. Figure 4 4 shows that the deformed shape of the RVE that is under uniaxial x direction tension without constraining boundary condition. As both top and bottom surface has a bulged out shape due to ha rd fiber in center, it cannot compose a periodic structure and consequently cannot be used as a representative volume element. Thus, proper boundary condition is required so that the behavior of the RVE can represent entire periodic structure. T he RVE re quires constraints that related the nodes on the opposite side of RVE so that the opposite sides deform in the same shape, which make the geometric periodic. Figure 4 4 Volume element deformed without periodic boundary cond ition Periodic Boundary Condition First, it would be better to know about the displacement boundary condition that sometimes applied to the RVE as a periodic boundary condition, which is not PAGE 49 49 appropriate It is w hen a RVE is subject to displacement field on the boundary in the form: ( 10 ) T he s denotes each boundary surface and is the size of the RVE in j direction When is a constant strain the average or effective strain is same as the applied constant strain, i.e., assuming fiber and matrix is perfectly bonded. Application of this homogeneous displacemen t boundary condition to a RVE results in the flat surface remain flat after deformation This condition is in appropriate because the RVE is heterogeneous and has a fiber that is generally harder than the matrix in the center ; the deformation of the surface cannot remain flat after deformation If boundary surfaces are forced to be flat when it has to deform to have wavy surface, it will be over constrained and the result will be different. On the contrary, t he periodic displacement boundary condition constr ains the boundary in pair facing opposite each other, i.e. the plane at the coordinate x=0 and the plane at the coordinate x=1 when the dimension of the RVE is unit (d x =1). The periodic displacement boundary condition constrains the boundary to keep the re lative displacement constant according to the strain on that boundary as shown in Figure 4 5 Figure 4 5 Periodic displacement boundary condition in x direction PAGE 50 50 The periodic displacement boundary conditions can be assigned considering the deformation of a RVE relative to a fixed global coordinate. The displacement u(x 0 ) of a point x 0 in a n x coordinates in RVE need to be considered first. Then the characteristic distance d that does not dependent on a location of x 0 defines the size of the RVE in x direction. From this concept, the periodic boundary condition can be expressed as follows: ( 11 ) where is the average strain component. Y and Z direction also has perio dic displacement boundary condition same way. As this boundary condition is only available for the RVE that has opposite face, the periodic boundary condition cannot be applied when the RVE that has triangular shape or when the number of the nodes on the f acing planes are different. Besides to the periodic displacement boundary condition, the traction boundary condition should be satisfied as well. The traction boundary condition can be written as ( 12 ) w here t denotes traction on the boundary. Y and Z direction also has periodic traction boundary condi tion same way. Equation ( 11 ) and ( 12 ) define the periodic bo undary condition that hold for arbitrary microstructure in the RVE. T his periodic boundary condition is applied to both elastic and inelastic behavior analysis in this research. In practical implementation in Abaqus, the periodic boundary conditions are i mposed using multi point constraints (MPC). MPC can impose a relationship between PAGE 51 51 different degree of freedoms in the numerical model. In general, multiple dependent degrees of freedom can be related to a single independent degree of freedom. Application o f UVE in Structure The periodic boundary condition is applied to calculat ing elastic stiffness matrix from the RVE in this research and the detail is presented when homogenization is discussed The periodic boundary condition is inevitable to calculate eff ective stiffness matrix when macro sc ale strains are applied to the RVE. However, as the periodic boundary condition is not a perfect boundary condition, there ha ve been stud ies on its effect on homogenized material properties. Drago and Pindera presented the periodic displacement boundary condition as the upper bound of effective properties and the periodic traction boundary condition as the lower bound [57]. One disadvantage of periodic boundary condition is that it is limited to a symmetric shape of the RVE I t is well known that it has unrealistic ally stiff response on the boundary [1 09 ]. Since obtaining accurate effective mechanical properties is the key in homogenization it is necessary to verify the accuracy of periodic boundary condition first How ever, it should be noted that the periodic boundary condition is not a requirement for the homogenization process but for the application of the RVE concept to keep its periodic characteristics. It means that the periodic boundary condition itself has noth ing to do with the process of effective propert y calculation if the properties are available directly from the entire composite structure. F or the verification, the behavior of the RVE that has periodic boundary condition and the behavior of a volume regi on that is called UVE in the heterogeneous structure are compared in Figure 4 6 The UVE is a unit volume, as shown in Figure 4 7 which is geometrically same as RVE, in the PAGE 52 52 composite structure in which t he stress and strain is calculated The UVE is named to distinguish from the RVE. Figure 4 6 X direction stress in the RVE with periodic boundary condition and UVE Figure 4 7 The unit volu me element UVE in the composite structure It shows that the behavior of the RVE and the UVE is almost same in elastic region but not in the plastic region. Based on the assumption that the UVE represent more close to real behavior, the applic ation of the RVE with periodic boundary condition to the elastic analysis is allowable but not in plastic region. Thus, the properties obtained from the RVE are used to homogenize elastic behavior and UVE is PAGE 53 53 used to inelastic behavior. The size of composit e structure that surrounds UVE is discussed later for homogenization for inelastic behavior. Homogenization for E lastic B ehavior As the fiber is modeled elastic material, RVE will behave elastic ally when the matrix behaves in elastic region. E lastic behavi or of materials can be described by elastic mechanical properties such as elastic Young s modulus Poisson ratio and so on. Since individual constituents are elastic, the combination of them will also be elastic. These elastic mechanical properties can be narrowed down to the relation between stress and strain, represented by the stiffness matrix. Table 4 1 shows the elastic mechanical properties for the constituents used in the RVE. Table 4 1 Mechanic al property of constituents Ef (modulus of fiber) 224GPa Em (modulus of matrix) 3.5GPa f (Poisson of fiber) 0.3 m (Poisson of matrix) 0.26 Vf (fiber volume fraction) 50% The constitutive law for the RVE can be determined based on th ese mechanical properties and t o have effective stiffness that can represent overall or global stiffness T he averag e scheme is used to calculate the macro stress, which generally defined as volume average in RVE as [ 11 0 ] ( 13 ) where is the true stresses in the RVE or micro stresses. PAGE 54 54 As RVE itself is composed of small elements, to perform vo lume integral, the Gaussian Quadrature Integration method is applied for each element ( 14 ) To have all the stress values in every meshed element in RVE, a commercial finite element program is used. E ach stress component in every meshed element is calculated and then u sing vector operation command, all the stresses having same index in element are summed up. As there are six stress indices for isotropic material the resultant effective stress that represents whole RVE stress state will be a 6 by 1 vector The effective constituents can be expressed as ( 15 ) As the strain s are given (in fact, it is applied as a boundary condition) and the stress es in each element are calculated, the component s of the stiffness matrix , can be obtained. U sing simple algebr aic manipulation, the effective stiffness component is obtained using FEM program and averaging scheme. For example, when is applied and all other macro scale strains are zeros, the averaged stresses become the first column of the effective stiffness matrix, i.e.: PAGE 55 55 ( 16 ) As we have effective elastic stiffness matrix, we can obtain effective global stress directly from applied str ain using constitutive equation without modeling both fiber and matrix together, which is computationally expensive The RVE is an orthotropic material but it can be t ransversely i sotropic materials when the uniaxial fiber is at the center of the RVE. The t ransversely i sotropic materials ( five independent constants) stiffness matrix is obtained by Abaqus analysis for RVE. ( 17 ) Furthermore, from following relation between compliance matrix, S, and the elastic property, we can calculate effective elastic properties from the stiffness matrix of the RVE. ( 18 ) PAGE 56 56 ( 19 ) The and are the extensional moduli of elasticity along the 1, 2, and 3 directions, r e spectively. Also, ij ( i, j = 1 2 , and , and are the shear moduli. It is clear that there are only five independent material constants ( ) for a transversely isot ropic material is dependent of and can be calculated by ( 20 ) A materia l is called isotropic when its property is the same in all three directions. In that case, There are only two independent material constants (E, ) for an isotropi c material. However, anisotropic materials have nonzero entries at the upper right and lower left portions of their compliance and stiffness matrices. The components of the compliance matrix ha ve following relation with stiffness matrix components. ( 21 ) ( 22 ) PAGE 57 57 ( 23 ) ( 24 ) ( 25 ) ( 26 ) Table 4 2 shows th e calculated effective elastic mechanical properties calculated by the stiffness matrix components and compliance matrix when the fiber volume fraction is 50% and the material properties in Table 4 1 are used As expected, the stiffness in the fiber direct ion (E 33 ) is an order of magnitude larger than that of other two directions. Table 4 2 Effective elastic properties Effective e lastic property Value [GPa] 10.92 113.90 0.22 0.28 2.97 4.09 The stiffness matrix obtained from RVE model can be compared to the stiffness mat rix obtained from Eshelby method as shown in Table 4 3 and Mori Tanaka method in Table 4 4 PAGE 58 58 Table 4 3 Effective stiffness matrix f rom Eshelby method 7.31 2.43 4.42 0 0 0 2.43 7.31 4.42 0 0 0 2.69 2.69 110.96 0 0 0 0 0 0 3.44 0 0 0 0 0 0 4 0 0 0 0 0 0 4 Table 4 4 Effective Stiffness matrix from Mori Tanaka method 10. 40 3. 62 3.9 8 0 0 0 3. 62 10. 40 3. 98 0 0 0 3.9 8 3.9 8 116.01 0 0 0 0 0 0 3.3 9 0 0 0 0 0 0 4 0 0 0 0 0 0 4 The Voigt and Reuss method which is known to give upper and lower bounds of the stiffness components are calculated The results in Tables 4 5 and 4 6 show that the calculated elastic properties are within the bounds. The effective elastic properties are compared in Table 4 5 and Table 4 6 for volume fraction 50% In fiber volume fraction of 50% axial direction estimates , are close to each other but in transverse direction , the analytical methods give big discrepancy. The Table 4 5 Young' s modulus compared to upper bound ( V f=50%) Abaqus RVE Eshelby Mori Tanaka 113.90 106.95 113.75 114.35 Difference 0% 6.0% 0.1 % + 0. 4% Table 4 6 Young's modulus compared to lower bound ( V f=50%) Abaqus RVE Eshelby Mori Tanaka 10.93 6.42 9.08 6.93 Difference 0% 41. 3 % 16.9 % 36.6 % PAGE 59 59 Effect of the Volume Fra ction The fiber volume fraction is decreased from 50% to 5% and the stiffness matrix and the elastic properties calculated based on the RVE, Eshelby and Mori Tanaka method are compared. As can be seen in Tables 4 7, 4 8, and 4 9, all three methods yield si milar effective elastic properties. This is because for small volume fraction, the interaction between fibers can be negligible. Table 4 7 Effective Stiffness from ABAQUS (Vf= 5%) 4.62 1.61 1.64 0 0 0 1.61 4.62 1.64 0 0 0 1 .64 1.64 15.32 0 0 0 0 0 0 1.49 0 0 0 0 0 0 1.53 0 0 0 0 0 0 1.53 Table 4 8 Effective stiffness from Eshelby method (Vf= 5%) 4.59 1.60 0.91 0 0 0 1.60 4.59 0.91 0 0 0 1.62 1.62 14.95 0 0 0 0 0 0 1.59 0 0 0 0 0 0 1.65 0 0 0 0 0 0 1.65 Table 4 9 Effective Stiffness from Mori Tanaka method (Vf= 5%) 4.61 1.62 1.64 0 0 0 1.62 4.61 1.64 0 0 0 1.64 1.64 15.39 0 0 0 0 0 0 1.50 0 0 0 0 0 0 1.53 0 0 0 0 0 0 1.53 The lower the fiber volume fr action, the closer the Eshelby method to the Abaqus result. However, the Mori Tanaka method is always closer than that of the Eshelby method Table 4 10 PAGE 60 60 and Table 4 11 Th erefore it can be concluded that for a low fiber volume fraction, analytical methods can predict the effective elastic properties accurately, but the Mori Tanaka method is closer than that of Eshe l by method. Table 4 10 Young's modulus compared to upper bound (Vf = 5%) Abaqus RVE Eshelby Mori Tanaka 14.45 14.47 14.53 14.53 Difference 0% + 0.1 4 % + 0.5 5 % + 0.5 5 % Table 4 11 Young's modulus compared to lower bound (Vf = 5%) Abaqus RVE Eshelby Mori Tanaka 3.98 3.99 3.97 3.68 Difference 0% + 0.2 5 % 0.2 5 % 7. 5 % As expected, the discrepancy in transv erse direction Young s modulus, reduced when fiber volume fraction is reduced compared to Table 4 6 However, the difference of axial direction Young s modulus for Mori Tanaka is slightl y increased when fiber volume fraction is reduced. Effect of the Fiber Crack So far, the effective elastic properties for perfect fiber matrix composites are considered. However, it is often possible that some of fibers are damaged, and as a result, the ef fective elastic modulus can significantly different from the perfect one. The effect of the fiber crack on elastic stiffness is investigated in this section The crack, which is a discontinuity, is modeled as enriched feature. The extended finite element m ethod (XFEM) uses the enriched feature and is an extension of general FEM allowing the discontinuities to exist in an element by enriching degrees of freedom with special displacement functions. This method is efficient especially when the crack propagatio n is to be investigated as the conventional FEM have to re mesh and update continuously PAGE 61 61 through the entire structure. In this research, since the effect of the crack size in fiber is the focus the propagation of the crack is not considered. Stationary cra ck can be defined using a n enrichment command and assigning crack domain in Abaqus. When the elements are intersected by the defined stationary crack domain, the elastic strength of that element is regarded as zero, which can be regarded as discontinuous. Figure 4 8 Transverse direction crack in fiber Figure 4 8 shows the crack in fiber covering the fiber cross section area of 50%. d the shape is assumed as half circular. The front of real crack in propagation direction would not have a straight line but it is assumed straight in this research without considering the effect of the crack surface shape for simplicity. The effective sti ffness is calculated numerically in the same way as it was calculated for the RVE that does not have crack. Table 4 12 shows how much percentage the stiffness degrades due to the crack when the cracked surface covers 50% of the cr oss section of the fiber. Table 4 12 Degradation by transverse direction fiber crack having 50% crack area +0.56 % 0.27 % 24.52 % 0 0 0 0.22 % +0.54 % 24.80 % 0 0 0 24.17 % 25.75 % 36.58 % 0 0 0 0 0 0 +0.77% 0.11 % 0 0 0 0 0 1.22% 0 0.01 % 0.01 % 0.49 % 0 0 +0.44% PAGE 62 62 Note that a significant reduction in stiffness is observed in while and are relatively unchanged. The stiffness reduction in Table 4 12 is the case when the crack is appeared periodically, which is not a common expectation of crack in a fiber. Normally, a crack exists only at a specific l ocation, not periodically. Therefore, if the effect of crack is going to be modeled in macro scale structure, it is necessary to use the reduced stiffness in Table 4 12 only for the location that has a crack. All other locations should use the nominal effe ctive stiffness. In addition, it is important to note that the stiffness reduction in Table 4 12 is proportional to the depth of crack in a fiber. Homogenization for Inelastic B ehavior When a material shows a nonlinear behavior, its behavior cannot be des cribed using effective stiffness. Instead, a physical model, which comes with model parameters, that describes the behavior of the material is required. Therefore, homogenization of inelastic behavior means to identify the model parameters that make the ho mogenized material show a similar behavior with the heterogeneous material. In this section, the anisotropic Ramberg Osgood model is used as a model to describe the inelastic behavior of fiber matrix composites When the fiber is an isotropic elastic mater ial and the matrix is an isotropic hardening elastoplastic material, the fiber reinforced composites behaves like an anisotropic hardening elastoplastic material. Unlike to the elastic behavior, the constitutive relation for elastoplatic material cannot be expressed using one parameter such as stiffness matrix and even becomes more complex when the hardening behavior is anisotropic. Therefore homogenization is far more complicated than the elastic behavior. As discussed earlier, the UVE not RVE, is the vo lume of interest to obtain effective properties for inelastic homogenization. To obtain the properties from UVE, the PAGE 63 63 appropriate size of the heterogeneous structure should be decided first because the material behavior would be different from the volume el ement that is placed at the most outer area to the one placed in the middle of the structure A B Figure 4 9 Heterogeneous structure composed of UVE piled up A ) Y direction B ) Z direction Figure 4 9 shows that the deformed shape of the structures composed of several fibers in y direction and z direction under x direction displacement loading respectively From F igure 4 9 and Table 4 13 because of symmetry, three of UVE s in y direction pile up seems to be enough for the size of the heterogeneous structure base model and also three of UVEs in z direction would be enough. Table 4 13 Effective stress in UVE according to the ir position in each direction UVE position Y direction Z direction The 3rd UVE (red square) 10.54 MPa 10.32 MPa The 2nd UVE 10.54 MPa 10.32 MPa The1st UVE 10.32 MPa 10.21 MPa The difference in stress values between y direction and z direction come from different geometric constraints. However, the point is three UVE, when symmetry PAGE 64 64 boundary condition is applied, is enough for the size of the structure. Figure 4 10 shows the heterogeneous structure to be analyzed to obtain effective inelas tic properties. A B Figure 4 10 Struct ure to calculate inelastic properties in UVE A ) with a symmetric boundary condition B ) Mises stress distribution in the structure under x direction tension. As is discussed in C hapter 2 the classical plasticity approach has weakness in terms of the computational cost for homogenization and it is more complicated to model when statistical analysis should be performed. On the contrary, the Ramberg Osgood model does not require exp ensive iteration process es to obtain the plastic multiplier which is essential in the classical plasticity model and even the yield surface is not required. In addition, anisotropic hardening application to anisotropic Hill s criteria in the classical plas ticity theory is not simple The Ramberg Osgood model, which was originally introduced for one dimension nonlinear behavior, has been modified for three dimension application s But, most of them were for homogeneous material and mostly for the isotropic ha rdening model In this research, three dimensional Ramberg Osgood model is associated with anisotropic hardening heterogeneous fiber composites so that it can represent the anisotropic behavior of homogeniz ed composite material T he key PAGE 65 65 to homogenization o f multiaxial anisotropic hardening composite material is to obtain Ramberg Osgood parameters and anisotropy tensor from heterogeneous composite material behavior Nonlinear least square s method is used to find optimum values for its parameters so that it h as minimum discrepancy with the behavior of heterogeneous material Anisotropy tensor is origin ated from the Hill s anisotropy yield criteri on. However, t he Hills coefficients are found that they are limited in represent ing anisotropic behavior of RVE, esp ecially when one direction shows strong elastic behavior. Three dimension Ramberg Osgood model The three dimensional Ramberg Osgood model is implemented to Abaqus FE program through the user subroutine coded by FORTRAN. Equation ( 27 ) shows the tensor form of anisotropic Ramberg Osgood model used in this research The difference from the original Ramberg Osgood model is that it is written in the tensor form and has an anisotropy tensor, M. ( 27 ) : Effective strain and effective stress : Effective c ompliance tensor : Anisotropy tensor to be calculated : Effective d eviatoric stress vector : Anisotropic equivalent stress A,TY and B: constants to be calculated The compliance matrix , wa s obtained from the inverse of the effective stiffness matrix obtaine d from the homogenization procedure for elastic behavior of the RVE. PAGE 66 66 To calculate equivalent stress which is defined in equation ( 28 ) anisotropy tensor shou ld be calculated. ( 28 ) The equivalent hydrostatic stress The stress deviator For orthotropic materials, the components of anisotropy tensor, M, can be [10 4 ] parameters are known as material properties and they are independent on the loading state. The anisotropy tensor in the equation ( 27 ) restricted to orthogonal materials but our model has fiber in the center of matrix and the x direction and y direction symmetry made it identical to transversely i sotropic materials. ( 29 ) Components in anisotropy tensors can be found from six ma terial tests that are carried out on the heterogeneous UV E FE model. Three uniaxial tensile tests and three uniaxial shear tests in the principal directions are required but symmetry reduced the number of tests. Yield stress ratio is defined by the measured yield stress in particular direction compared to the reference yield stress, which is in this research. PAGE 67 67 ( 30 ) From the uniaxial tests, the components are obtained. ( 31 ) As the reference yield stress is x direction, symmetry in x direction and y direction makes G=0, and the elastic behavior in z direction made F=0. The shear part L, M and N is defined by the relation in Mises criteria, but this does not apply to the heterogeneous materials behavior. Most of all, Hill s anisotropy parameters are for anisotropic yield point and they do not have any information on anisotropic hardening behavior. Therefore anisotropic hardening should be included by Ramberg Osgood parameters and the shearing component s in anisotropy tensor can be obtained together when the parameters are calculated through optimization procedure. PAGE 68 68 Parameters C alculation First, the heterogeneous composite UVE is modeled using Abaqus, a nd the stress strain relation is calculated by gradually increasing strain in each direction. Then, the goal is to find the Ramberg Osgood model parameters of a homogeneous material that describe a similar behavior with the Abaqus model. The anisotropic mu ltiaxial Ramberg Osgood model s parameters and anisotropy tensor are obtained using nonlinear least square fit (NLSQ) function in Matlab. The data that Ramberg Osgood parameters are fitted is obtained from Abaqus FE analysis of UVE in heterogeneous structu re in x direction and x y shear direction The NLSQ finds parameters that minimize the error which is defined as the square sum of the difference between the Abaqus calculated stress and the stress from the Ramberg Osgood model at the same strain values (se e equation ( 32 ) ). ( 32 ) In addition to the Ramberg Osgood parameter s the anisotropy tensor component, M 44 that is calculated at equation ( 31 ) A ll parameters in the model are calculated as optimized values for the Ramberg Osgood model so that it can represent the behavior of the heterogeneous material model. Table 4 14 shows th e optimized Ramberg Osgood parameters obtained from this procedure for given mechanical pr operty of matrix. SY1 is the first yield point of the matrix material s stress strain curve and SY2 is the second yield point that matrix material begins perfect plastic deformation. PAGE 69 69 Table 4 14 Optimized Ramberg Osgood paramet ers for given matrix materials properties # SY1 SY2 A TY B M44 Error 1 39.75 61.83 7.67E 08 30.06 6.28 3.90 1.85E 01 2 67.40 3.91E 08 27.27 5.95 3.78 1.83E 01 3 70 .00 7.48E 08 31.84 5.98 3.81 1.82E 01 4 72.60 7.14E 08 32.29 5.90 3.81 1.81E 01 5 78.17 6.78E 08 32.34 5.55 3.81 1.78E 01 6 43.33 61.83 1.00E 09 20.40 7.50 3.84 1.83E 01 7 67.40 4.16E 08 31.02 6.51 3.82 1.80E 01 8 70 .00 7.37E 08 31.85 5.98 3.79 1.78E 01 9 72.60 6.81E 08 32. 30 5.90 3.79 1.77E 01 10 78.17 7.76E 08 32.65 5.52 3.74 1.75E 01 11 45 .00 1.83 1.00E 09 20.41 7.50 3.83 1.82E 01 12 67.40 1.00E 09 21.00 7.23 3.79 1.79E 01 13 70 .00 3.35E 08 31.11 6.49 3.79 1.77E 01 14 72.60 7.03E 08 32.31 5.90 3.79 1.76E 01 15 78.17 7.55E 08 32.65 5.52 3.75 1.74E 01 16 46.67 61.83 1.00E 09 20.88 7.60 3.80 1.82E 01 17 67.40 1.00E 09 21.30 7.29 3.79 1.78E 01 18 70 .00 1.00E 09 21.27 7.12 3.76 1.76E 01 19 72.60 7.23E 08 37.55 6.68 3.73 1.75E 01 20 78.17 7.76E 08 32.75 5.53 3.74 1.73E 01 21 50.25 61.83 1.00E 09 22.56 7.98 3.86 1.83E 01 22 67.40 8.02E 08 42.89 8.00 3.77 1.76E 01 23 70 .00 8.58E 08 44.42 8.00 3.67 1.75E 01 24 72.60 1.00E 09 22.45 7.19 3.72 1.73E 01 25 78.17 2.76E 08 32.54 6.23 3.73 1.70E 01 To see gra phically how close the homogenized model is to the heterogeneous RVE model, the effective stress strain curves for both cases are compared in Figure 4 11 Each figure has both heterogeneous materials data that is marked by crosses and homogenized Ramberg Osgood model data marked by circles. Figure 4 11 (a) shows the effective stress in x direction when the x direction displacement is applied to the heterogeneous materials and homogenized materials that has of 39.75MPa (SY1) and (b) of 50.25MPa(SY5) while keeping the second yield stress of 78.17MPa. It corresponds to the comparison between #5 and #25 in Table 4 14 PAGE 70 70 A B Figure 4 11 X direction average stress in homogenized material and in heterogeneous material A) for the first yield stress of 39.75MPa B) for 50.25MPa. The elastic region matches well for either case but as the strain increases, the discrepancy also inc reases. As the Ramberg Osgood parameters are from best curve fit by NLSQ to normal tension data and shear data at the same time, the shear behavior also needs to be verified. Figure 4 12 shows when xy shear displacement is applied for the same model. A B Figure 4 12 XY shear direction average stress in homogenized material and in heterogeneous material varying the first yield stress PAGE 71 71 This shows the optimized Ramberg Osgood model well matches for shear stress strain relation but a little off in normal loading state. One thing to note is that A ) and B ) curves are almost identical in both Figure 4 11 and Figure 4 12 It means the variation in the f irst yield point does not affect the elastoplastic behavior of the composite material. However, Figure 4 13 shows that the variation of second yield points affects the behavior of materials noticeably A B Figure 4 13 XY shear direction average stress in homogenized material and in heterogeneous material varying second yield stress As it is proven at this point that the first yield stress of matrix is not critical parameter for macroscopic beh avior, it is not considered in uncertainty analysis as an input parameter while the second yield stress is considered. In contrast to the x and y direction, the composite material behaves almost linear elastic in z direction which is axial direction of fib er. The comparison between heterogeneous materials and homogenized model is shown in Figure 4 14 It shows that the Ramberg Osgood model can represent linear elastic behavior in z direction accurately and even identical to the com posites behavior. PAGE 72 72 Figure 4 14 Z direction average stress in homogenized material and in heterogeneous material Combined L oading C ase The general loading cases in working environment is the combination of the two loading types, normal stress and shear stress. So, the parameters that were obtained under unidirectional tension and unidirectional shear loading should satisfy arbitrary complex loading condition as well when it is homogenized. To verify this combined loading c onditions are applied and compared their behaviors One is the mixture of shear and tensile displacement loading and the other is biaxial normal loading. For shear and tensile loading mixture, the same amount of displacement is applied to both heterogeneou s model MPa and second yield stress is 70 MPa and corresponding homogenized model. The shear strain of 0.04 and t he tensile strain in x direction of 0.04 are applied at the same time. Figure 4 15 shows the x directions stress contour under the combined loading case compose of shear and tensile loading. The calculation time used for ( A ) is 12 minutes and for ( B ) 4 seconds which is a significant difference considering the size of the analysi s. PAGE 73 73 A B Figure 4 15 Deformed shape and x direction stress under shear and tensile combined loading A) Heterogeneous materials B) Homogeneous materials However, as the homogenized model has only one element, it cannot make similar deformed shape as heterogeneous materials. However, its effect is important and need to be studied in the future. Figure 4 16 s hows that the homogenized model works well in shear and tensile combined loading conditio n. T he elastic behavior of the homogenized material is almost identical to heterogeneous one and the plastic region also shows very similar behavior. The stresses in y direction and z direction are not plotted because they show zero stresses under this loa ding condition. A B Figure 4 16 Inelastic behavior comparison for heterogeneous and homogenized model A ) X direction average stress under combined load B ) XY shear average stress PAGE 74 74 For the verification, another combine d loading condition is applied. X direction normal strain of 0.04 and Y direction compressive normal strain are applied to the heterogeneous RVE and homogenized model as well Figure 4 17 shows the x directions stress contour unde r biaxial combined loading condition. A B Figure 4 17 Contour of x direction stress under biaxial combined loading A) Heterogeneous materials. B) Homogeneous materials. Figure 4 18 s hows that the homogenized model works well in biaxial combined loading condition too. From th ese verification s it can be said that the homogenized materials using multiaxial anisotropic Ramberg Osgood model behaves very similar to heterogeneous RVE material not on ly in unidirectional but also in combined loading conditions PAGE 75 75 A B Figure 4 18 Average stress comparison under biaxial loading. A ) x direction B ) y direction Therefore the Ramberg Osgood parameters that are obtained from heterogeneous material data satisfy not only the uniaxial loading but also the elastoplastic stress strain relation under combined loading states. As we have a materials ana lysis. PAGE 76 76 CHAPTER 5 UNCERTAINTY ANALYSIS Introduction E xperimental researches in engineering have shown that model parameter s or material propert ies show a random, or stochastic character no matter what the material s are Heterogeneous material s such as a composite material ha ve a complex structure and, unfortunately, b ecause of high d emand in performance, the micro or macro structures become complicated This cause s the uncertainty in parameters relating to the composites to increase Even in a specific interest of scale, many parameters such as voids, crack, particles and grain bou ndar ies are involved but it is not easy to tell how much they are involved in failure or how these parameters may be correlated to one another The analysis based on exact values that do not take uncertaint y into account is called deterministic analysis. On the other hand, the uncertainty analysis considers uncertain ty in input parameters and their effects on results. For the uncertainty analysis, material parameters having statistical distribution are used as input s instead of deterministic values For un cert aint y modeling, the uncertainty in the model inputs should be characterized first. This can be represented by standardized normal random variables ( srvs ) with zero mean and variance of one. W hen other types of random variables are going to be used, a transformation is required An assumption on the input variables is that they are independent so that each variable can be expressed in terms of srvs through a transformation. When the input variables are ready, uncertainty propagation should be carried ou t with generated samples of input variables. Monte Carlo Simulation (MCS) is the most common method to generate random samples because of its simplicity in implementing on a computer PAGE 77 77 program and it does not require specific knowledge of properties However it has weakness that it requires large number of samples, on the order of thousands or even millions, to have acceptable precision in the response. This is not appropriate for computationally expensive models because it takes too long to consider all the samples created by MCS especially for FE analysis Stochastic Response Surface Method R esponse surface method is one of a regression method in statistics. Regression analysis is a statistical method that utilizes the relation between two or more variabl es in which dependent variables can be estimated by independent variables. General regression model is shown in equation ( 1 ) ( 1 ) As the equation shows the response is modeled with a combination of functions of variables and error. (i=1, ,p) are the coefficients, (i=1, p) are the terms or bases and is error. The indicates a vector of input variables and so on. When =1 and =x, then it is linear (first order) model. When =1 and =x and then it is a polynomial model. Both r esponse s urface m ethod (RSM) [ 11 1 ] and s tochastic r esponse s urface m ethod (SRSM) [ 73 ] has been used to approximate the output. Stochastic is synonymous with random and often used as c ounterpart of the word deterministic and consequently t he SRSM is an extension of RSM that is customized to the random inputs instead of deterministic values By transforming all input variables to the standard normal random variables and by using Hermi te polynomial s which are orthogonal to the Gaussian distribution, as the PAGE 78 78 bases functions, the SRSM provides many effective features in identifying probabilistic characteristics of random outputs. It provides robust estimates of the coefficients compared t o the one obtained using non orthogonal polynomials [11 2 ] and convergence is optimal for Gaussian processes [11 3 ]. The biggest difference of SRSM is the type of input. SRSM uses random variables that have forms of probability density function while RSM use s deterministic values. T hrough the SRSM, the uncertainty in input parameters propagates to the uncertainty in the response, which is explained in the following step. The first step for SRSM is to decide which parameters will be input variables and stati st i cally transform them to have standard random variables. As SRSM uses Gaussian distribution, the distribution of the random variable, let say X, can be represented by the transformation [ 11 4 ]. ( 2 ) where is the mean of X, is the standard deviation of X and is the first order orthogonal polynomial of the standard normal distribution which has a mean of zero and a standard deviation of one. The advantage of this standard normalization is that the same set of orthogonal polynomials can be used for any normal distribution instead of deriving a distribution specific set for each normally distributed parameter. The se cond step is to define output approximation. A functional approximation that takes into account all uncertainties in input should be chosen. Polynomial chaos expansion will be used and the collocation point method to select points will be used t o obtain th e coefficients of the response surface The third is to investigate the output of the approximation. As the approximation for the output represented by polynomial expansion is obtained, the statistical PAGE 79 79 characteristics of the output will be investigated th rough MCS. The statistical moments, probability density function, and correlation between input and output will be discussed. Since the evaluation of function using SRSM is very inexpensive, millions of samples can be generate to calculate accurate statist ical properties of the output. T he SRSM utilizes polynomial chaos expansion to represent the output and the coefficients of the polynomial expansion are calculated using a collocation approach [ 11 5 ]. The procedure is to evaluate the function at several po ints (and possibly the derivatives of the function) and then fit a polynomial to these known points by minimizing the error between function values and approximations The coefficients that minimize this error are then found and the approximation is consi dered to be a good estimate to the true function [ 83 ] T he polynomial expansion in SRSM uses Hermite polynomial bases and provides closed form solution of model output requiring significant ly lower number of simulation s than MCS. T he expression for the 3 rd order polynomials in n dimension is shown in equation ( 3 ) ( 3 ) is the model output, the represent coe represents standard random variables. The number of unknown coefficients is determined by dimension of the design space n which is the number of random PAGE 80 80 variables. For 3 rd order expansion, the number of unknown can be ca lculated using equation ( 4 ) ( 4 ) A lthough, the procedure using SRSM seems attractive, the challenge is how accurately the polynomial expansion can be established to replace computational model This includes developing proper method to generate and updating best representative points for surrogate to minimize the error. Collocation Points T h e collocation concept comes from t he idea of Gaussian Quadrature in numerical integration [ 11 6 ] The c ollocation points are selected based on a modification of the standard orthogonal collocation method. T his method is called the orthogonal collocation method [ 11 7 ]. The points are selected so that standard normal random variables, U i in equation ( 3 ) takes the values of either zero or one of the roots of the higher order Hermite polynomial. Thos e collocations points are in high probability region of the input distribution and it is the feature that makes this kind of effective collocation method cost effective method. The theory behind collocation method is based on the concepts of Gaussian Quadr ature It is a numerical integration technique and uses orthogonal polynomials for the selection of points. Equation ( 5 ) is the typical form of integrals in G aussian Quadrature. ( 5 ) PAGE 81 81 where g(x) is an orthogonal polynomial and f(x) is a non negative weighting function defined in the connected space F. Gaussian Quadrature seeks to obtain the best numerical estimate for the integral by choosing x values, evaluating g(x) at the points and computing the integral. The x values are the roots of the orthogonal polynomials. T he result of Gaus sian Quadrature integration is shown in equation ( 6 ) ( 6 ) The coefficients depend on the weighing function and g(x) is evaluated at abscissa values that are the roots of th e n th orthogonal polynomial calculated with respect to the weighting function f(x). The integral becomes exact when g(x) is a polynomial of degree of 2n 1. It means integral can be estimated using only n samples. Table 5 1 shows H ermite Polynomial when the number of unknown value is five and polynomial is 3 rd order. The collocation points determine the samples of input variables. Using the selected input variables which are deterministic values computational model is used to eval uate the deterministic values of the output. Then, the coefficient s of the polynomial expansion can be obtained using the least square method. This is the strength of regression method that presents robust means to estimate coefficients because it uses mor e points than collocation points. When we have the output represented by a polynomial expansion, not by time consuming computational model we can use the MCS method using the SRSM to analyze uncertainty propagation. C onsequently, we replaced the time con suming computational procedure with polynomial equation to apply MCS that requires huge amount of iterations. PAGE 82 82 Table 5 1 The third order H ermite polynomial basis for five input variables ( P represents the order of polynomial ) Order P=0 1 P=1 P=2 P=3 Parametric Design P arametric model creation of the RVE for the FE analysis is an important technique for uncertainty analysis because building a response surface requires a set of deterministic output data and it is impractical to make RVEs one by one for every change of variable manually. For example, when the position and size of the fiber in matrix are variables, they can be parameterized in a code that makes a loop to generate several RVE models automaticall y. The corresponding code in Abaqus is the Python script A Python script to create RVEs can be written from the scratch but utilizing files PAGE 83 83 that Abaqus GUI environment generates would save a lot of time for whom the scripts are not familiar. Abaqus GUI by which user makes parts and does analysis clicking menus generate the objective language Python automatically. T he files with extension of rec and jnl are the files written in Python code for temporary and permanently. For example when elastic modulus of 100GPa and Poisson ratio of 0.3 are given in Abaqus GUI menu, then Abaqus automatically write in a file the following code: mdb.models [ myModel 1 ] .materials [ myMatrl ] .Elastic (table=(( 1 0 0 000.0, 0.3), )) R eplacing the values to variable and making a loop, each model can be created automatically at ever y iteration Figure 5 1 shows examples of RVEs in various geometries of the fibers generated by the parameterized code. The Python code is also used to access and make use of th e output data Abaqus calculated during the analysis. Figure 5 1 Examples of RVE generated by parameterized design Uncertainty P ropagation in E lastic B ehavior Although there are many parameters in the com posite material, three input variables are chosen for calculating the components of the effective stiffness matrix for elastic behavior First two, Ef and Em, are respectively elastic moduli of fiber and matrix, and Vf is the volume fraction of the fiber. These three variables are the most fundamental ones in the effective properties of composites PAGE 84 84 As we have three input variables and applying to 3 rd order polynomial, 20 coefficients need to be calculated. T he roots of 4 th order polynomial, which are and zero are used as collocation points to fit data to polynomial. This method is called ECM (Efficient Collocation Method) [54]. For collocation points, the three combinations of these roots are used because we need U 1 U 2 and U 3 w hich represents three random input variables. This result in 5 3 (=125) combinations. Generally, more than two times of the unknown parameter is required for the sample number. This means about 40 samples are enough to find 20 unknown coefficients but the number of combination that five collocation points can make is 125. So, we have far more points than we need F inding efficient way to choose proper 40 points out of 125 could be another study topic but in this research all of 125 poin ts are used to calculate the coefficients using least square fit Input Variables for FEM Analysis In general, the input variables may have different distribution types but for the purpose of SRSM, all input variables are transformed to the standard norm al random variables These input variables are used in SRSM but it should be used in numerical analysis to have physical output values. F or the FEM analysis, these input variables should be transformed to have physical meaning, i.e. 100MPa or volume fracti on of 0.5. The transformation is performed by a transformation equation that transform s collocation points to the values that has physical meaning. ( 7 ) PAGE 85 85 where a standard normal random variable, while physical input variable In the above equation is is the standard deviation of the random variable X For the followin g numerical example, the s tandard deviation of input variable is assumed to be 5% coefficient of variance ( COV = / ) The procedure is to find the value of physical variable X for a given collocation point U and statistical parameters the input va riable. Then, using the physical variable X, the finite element analysis is performed with the RVE and obtained corresponding data output. As we input the radius of fiber not the volume fraction of fiber into FEM, we need to have specific volume fraction. Equation ( 8 ) shows the input variables for the radius. ( 8 ) FE A nalysis of the RVE The stiffness matrix will be calcu lated putting the deterministic values of the transformed input variables to have physical meaning into the heterogeneous Abaqus RVE model T his was enabled by parameterized design of the RVE. Instead of constructing RVEs for every size of the fiber one by one manually, the radius of the fiber is parameterized so that only changing the value of the radius can construct the corresponding RVE automatically. The mechanical parameter such as fiber s Young s modulus and matrix s Young s modulus are parameterized too. T he geometric parameterization is only possible by the Python code that Abaqus supports. Table 5 2 shows the input values of parameters that will be put into Abaqus FE analysis to have deterministic output, which is the comp onent of the stiffness matrix. As PAGE 86 86 each parameter has five values, the combination of all parameters will make total number of 125 FE analysis data. S o, each stiffness matrix component has 125 data and it is saved as a vector. Table 5 2 Input values of parameters for FE anlysis #1 #2 #3 #4 # 5 E f 197854 .00 215690 .00 224000 .00 232310 .00 250145 .00 E m 3091. 5 0 3370. 2 0 3500 .00 3629.8 0 3908.5 0 Vf 0.44 0.48 0.5 0 0.5 2 0.5 6 Equation ( 3 ) can be rephrased to have simpler expression in terms of coefficients and three input variables as C= ( 9 ) As we have {c} obtained from FE analysis and input variables {U i } in equation ( 9 ) we can calculate its coefficients by matrix algebra and we can have stiffness matrix components expressed in form of polynomial expansion As long as we have this polynomial expansion for stiffness matrix components Monte Carlo Simulation can be applied to th e expansion for the statistical investigation for the propagation of uncertainty Table 5 3 shows the twenty coefficients for six stiffness components calculated by least square fit method. PAGE 87 87 Table 5 3 Coefficients of polynomial expansion for stiffness matrix components Bas e s C11 C12 C13 C33 C44 C55 1 13.517 4.116 4.5 00 115.615 2.833 3.901 0.03 0 0.004 0.008 5.581 0.003 0.008 0.646 0.202 0.217 0.2 00 0.139 0.187 0.943 0.076 0.2 00 5.565 0.149 0.275 0.001 0 0 0 0 0 0.003 0 0.001 0 0 0.001 0.006 0 0.001 0.28 0 0 0.001 1 0.001 0 0 0 0 0 0.042 0.004 0.009 0.001 0.007 0.013 1 0.06 0 0 0.012 0.004 0.009 0.017 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.001 0 0 0 0 0 0.001 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.002 0 0 0 0 0.001 0.004 0 0.001 0.002 0.001 0.001 Monte Carlo Simulation To see the propagation of uncertainty 50,000 samples from the normally distributed random variables are generated in Matlab and used to calculate the output of stiffness matrix component values from polynomial. Then, using 50,000 values of ca lculated components, the mean and COV of the output are calculated. Table 5 4 shows the comparison between input COV and output COV. The COV of input variables (Ef, Em and Vf) was equally 0.05, which means the degree of dispersion of those three input variables were same. But, the output shows the degree of dispersion is different. The stiffness matrix component C33 shows the least variation while C11 and C55 show the most variation, which means uncertainty in input variables prop agated most. PAGE 88 88 Table 5 4 Output statistics of p olynomial expansion Method Statistics C11 C12 C13 C33 C44 C55 SRSM Input COV 0.05 MCS (Ef, Em, Vf) Output COV 0.09 0.05 0.07 0.04 0.07 0.09 Thus, we can pay more attention on the C55 direction shearing stiffness when we design uniaxial fiber composites considering modulus of constituents and volume fraction of fiber. Figure 5 2 is the histogram stiffness matrix components and shows how each component i s distributed in range A B C D E F Figure 5 2 H istogram of stiffness matrix component A) C11 B)C12 C)C13 D)C33 E)C44 F)C55. A lthough some COV values are increased from 5% to 9%, the uncertainty in the m icro scale does not amplified much in the macro scale response This can be understood if the scaling up process is considered as an integration process and t he PAGE 89 89 fluctuation in the local material properties are averaged out after integration. Thus, the COV of the macro scale properties remains in the same level with the micro scale properties. However, this is only true when both materials are linear elastic. Uncertainty P ropagation in I nelastic B ehavior F or elastic behavior of the heterogeneous RVE, it was described by an effective stiffness matrix and each component was expanded into polynomial expression for statistical investigation. However, the elastoplastic behavior of the heterogeneous RVE which is modeled using multiaxial anisotropic Ramberg Osgood model, has four parameters, which were A, TY, B and M 44 These parameters are dependent each other that they could not be described as one values such as an effective stiffness matrix component in elastic condition. For this reason, each parameter cannot b e expanded into polynomial expansion independently Therefore instead of analysis on the Ramberg Osgood parameters, the strain energy stored in loading direction is considered in uncertainty analysis Input Variables for FEM Analysis The same input varia bles as the previous elastic analysis are considered first. The Young s modulus of fiber, the Young s modulus of matrix and the volume fraction of the fiber are verified first. As it has same number of input variables and degree of polynomial, same colloca tions points are transformed to physical values and put into the Abaqus heterogeneous RVE model. Before calculat ing deterministic values of strain deformation energy, the elastoplastic behaviors of the heterogeneous RVE were investigated. Figure 5 3 shows that how the heterogeneous RVE behaves in inelastic fraction change. First, the only difference between A ) with C PAGE 90 90 fiber, but two curves are identical. The RVE in A ) and D modulus of matrix and the curves are not significantly different. This means the into the un certainty of deformation energy. On the contrary it clearly shows that the volume fraction of fiber makes change in the macro scale response. Figure 5 4 shows the Young s modulus of fiber affects the behavior of the heterogeneous RVE in z direction which shows a similar trend with Figure 5 3. A B C D Figure 5 3 Elastoplastic behavior of heterogeneous material modulus of fiber and matrix and volume fraction of fiber A) Ef1= 197.85GPa Em1= 3.09GPa B) Ef 5 = 250.1 4 GPa, Em 1 =3. 09 GPa C) Ef 3 = 224.00 GPa Em 1 =3. 50 GPa D) Ef1=197.85GPa, Em 5 =3. 91 GPa PAGE 91 91 inelastic behavior of homogenized material, whereas t he volume fraction shows the largest contribution to the output variation A B Figure 5 4 Z direction behavior of heterogeneous RVE A) Ef1= 197.8 5 GPa, Em1= 3.09GPa B) Ef 5 = 250.15 GPa, Em 1 =3. 09 GPa Thus, two input variab les are chosen for uncertainty analysis in elastoplastic behavior : fiber volume fraction, Vf, and second yield point of the matrix, SY2. First, for the volume fraction, the collocation points are transformed to physical values using the transformation equation used before and then so does the second yield point of the matrix. Table 5 5 shows the physical input var iables that will be put into Abaqus heterogeneous RVE model. Table 5 5 Input values of parameters for elastoplastic behavior #1 #2 #3 #4 # 5 Vf 0.44 0.48 0.5 0.5 2 0.5 6 SY2 61.83 67.40 70 .0 72.60 78.17 FE analysis of the RVE in E lastoplastic B ehavior X direction strain energy under x tension displacement loading (3 rd order polynomial with two variables) can be expanded into polynomial as PAGE 92 92 ( 10 ) W e have 25 points of design of experiment for 10 unknown coefficients of polynomial. The least square method is used to calculate coefficients In a similar way, three kinds of loading conditions ar e applied, which are x direction tension displacement xy shear displacement and xy shear tension combined displacement condition. T he x direction strain energy is calculated for x direction loading and xy shear direction strain energy is calculated for xy shear loading and so on. When the coefficients are obtained, MCS is performed to estimate the distribution of the output response which is the strain energy in each direction Using deterministic data from Abaqus FE analysis the coefficients polynomial r esponse surface for each loading conditions are calculate d and shown in Table 5 6 Table 5 6 The coefficients of polynomial expansion for strain energy in x direction X tension 3.707 0.131 0.222 0.002 0.007 0.013 0 0 0.001 0.001 XY shear 1.243 0.039 0.021 0.001 0.002 0 0 0 0 0 Combined ( tension, shear ) 3.506 0.127 0.214 0.002 0.007 0. 014 0 0 0 0 0.858 0.038 0.001 0 0.002 0.001 0.001 0 0 0 Monte Carlo Simulation To investigate the propagation of uncertainty 50,000 samples from the normally distributed random variables are generated in Matlab and used to calculate the output of s train energy in each direction. PAGE 93 93 Table 5 7 The propagation of the uncertainty for elastoplastic behavior Input COV (Vf, SY2) Output COV (Energy) X tension 0.05 0.07 XY shear 0.05 0.04 Combined ( tension, shear ) 0.05 0.07 ( x normal direction) 0.04 (xy shear direction) Table 5 7 and Figure 5 5 shows that the uncertainty in volume fraction and second yield point of the matrix propagated increasing the x direction strain energy but not significantly. It seems the uncertainty in volume fraction and the second yield point propagates to the strain energy without being changed. A B C D Figure 5 5 H istogram of strain energy for each l oading condition in elastoplastic behavior A) x direction under x tension B) xy shear direction under xy shear load C) x direction under combined load D) xy direction under combined load. PAGE 94 94 CHAPTER 6 SUMMARY AND CONCLUSI ONS A homogenization method is pro posed t o reduce the computation cost for fiber reinforced composite material simulation. Not only for elastic behavior but also elastoplastic behavior was homogenized so that both have same average stress strain constitutive relations. T he elastic behavior of the fiber reinforced composites was homogenized through the RVE concept and the effective elastic mechanical properties were obtained by numerical approach using the user subroutine in Abaqus FE program. The application of RVE concept and homogenizatio n in this research has significant meaning in that it significantly reduce s computation al time especially for elastoplastic behavior analysis. The Ramberg Osgood model parameters are successfully calculated from heterogeneous material s FE data using leas t square fit to find optimal values. The successful implementation of the multi axial anisotropic hardening Ramberg Osgood model into Abaqus FE program made it possible the fast and effective homogenization for inelastic behavior Through numerical example s, it is shown that the proposed model works well with combined loading conditions not to mention simple uniaxial tension and shear loadings. However, there was a little discrepancy in normal direction behavior while shear behavior matched identical to het erogeneous materials I t seems some coupling between normal and shear direction occurs for combined loading for heterogeneous. However it is shown that the boundary condition is important in RVE analysis for periodic heterogeneous materials. The results shows strong dependency on the boundary condition and it needs to be studied more in further researches. In this research, UVE model is used instead of RVE with periodic boundary condition. It could PAGE 95 95 be said that the boundary condition must be studied deep l y when the research is on multi scale analysis. As the homogenized model in this research is simple compared to complexity of real composite materials, it is possible that important factors and properties can be blurred and buried under homogenization. Al though it is true consideration of applying homogenization should not be disregarded because full scale analysis is still impractical and design and manufacturing requires computational tools More than that, it is shown that it is practical and worth to studied more. For uncertainty analysis, s modulus of fiber and matrix, volume fraction of fiber, and yield points of matrix were investigated to see how much of the uncertainty in them will propagate to macro behavior of homogenized composites F or elastic behavior, the uncertainty in those parameters propagated to the axial direction stiffness increasing the uncertainty while it did not in other directions. For elastoplastic behavior, among the input variables that considered in this research the f irst yield point of matrix part and Young s modulus of both fiber and matrix did not give much effect on the strain energy variation during inelastic deformation. But the uncertainty in volume fraction and second yield point of matrix part propagated to th e uncertainty in strain energy when it is under normal direction tension loading. Although the amplification of the uncertainty to macro behavior was not significant, it was enough to make it possible to distinguish which parameters give more effect on the behavior of the material in specific condition. The homogenization which is oriented to a specific purpose can save computational cost or time significantly and enables bigger scale analysis. Along with homogenization, the uncertainty analysis can reduce the time spent on trial and error PAGE 96 96 and the time spent on finding which parameter should have more attention in design or manufacturing process. If experimental data was obtainable for this research, it would provide better information A lthough the experim ent data is essential, it is worth to note that the detail investigation for micro scale using SEM or TEM in experiments only reveals the response, not providing much about mechanism. T his means it is possible to choose wrong parameters or miss critical pa rameters when micromechanics are applied. When more reliable and accurate uncertainty analysis is possible, it could present more accurate information than the one based on experimental data or experiences of individual engineers. PAGE 97 97 APPENDIX PYTHO N CODE FOR PARAMETRI C ANALYSIS # # To run the Python when you start Abaqus/CAE # type the following command in command line in Abaqus/CAE # >>abaqus cae script=***.py # or to run a script file from the Abaqus command line interface # >>execf ile('***.py') # ========================================== # Parametric study WITH periodic displacement BC (using *equation) # ========================================== from part import from material import from section import from assembly import from step import from interaction import from load import from mesh import from job import from sketch import from visualization import from connectorBehavior import from math import from Numeric import # ***************************** ************************** # Define mechanical properties and variables # ******************************************************* fYoung=[179200, 224000, 268800] nuf=0.3 mYoung=[2800, 3500, 4200] num=0.26 # Width of RVE w=1.0 # Fiber center(cx,cy) cxr=[0.0 0.02, 0.05] cyr=[0.0, 0.02, 0.05] #Volume fraction(=Vf/V, V=1 and Vf=4*pi*R^2*1) Vfr=[0.4,0.5,0.6] # Cosine 45 (radian) PAGE 98 98 c45=cos(45*pi/180) # # Define a fu n ction to create Uniax ial Fiber Reinforced Composites def createPartRVE(myModel, Em, Ef, Vf): # # Define a function to apply the constant in Periodic Displacement Boundary Conditions def CEqApply(flna me, direction, dof): # # Define a function to apply the Periodic Displacement Boundary Conditions def CEq(flname, direction): # #def ends # ******************************************************* # Analysis loop # ******************************************************* volume=0 iter=1 iter1=1 iter2=1 iter3=1 for Em in mYoung: for Ef in fYoung: for Vf in Vfr: for cx in cxr: for cy in cyr: ... ... # # calcualte average stress # def aveS(Nframe, filename, stepname): ... ... iter=iter+1 #end of cy loop # end of cx loop # no iteratio n for cx and cy #end of Vf for loop iter1=iter1+1 #end of Ef for loop PAGE 99 99 iter2=iter2+1 #end of Em for Loop iter3=iter3+1 del mdb.models['Model 1'] PAGE 100 100 LIST OF REFERENCES [1] S.T. Peters, Handbook of composites, second edition, Chapman & Hal l, New York, 1998. [2] S.T. Peters, Handbook of composites, Process research, Chapman & Hall, s econd ed. ,1998 [3] B. Ravi Deo, Strarnes H. James Jr., Holzwarth C. Richard Low cost composite materials and structures for aircraft applications, RTO Applied Vehicle Technology Panel (AVT) Norway, 7 11 May 2001 [4] H. Wittich K. Friedrich, Interlaminar f racture e nergy of l aminates m ade of t hermoplastic i mpregnated f iber b undles, Journal of Thermoplastic Composite Materials 1 ( 1988 ) 221 231. [5] Harris E. Charles Starnes H. James An a ssessment of the s tate of the a rt in the d esign and m anufacturing of l arge c omposite s tructures for a erospace v ehicles, NASA/TM 2001 210844 (2001) 20. [6] G. D. Wyss, K. H Jorgensen, A u ser's g uide to LHS: Sandia's l at in hypercube sampling software, Sandia National Laboratories, Albuquerque, New Mexico, 1998. [ 7 ] H. Berger U. Gabbert, H. Koppe, R. Rodrguez Ramos, J. Bravo Castillero, R. Guinovart Diaz, J.A. Otero, G.A. Maugin Fini te element and asymptotic homogenization methods applied to smart composite materials, Computational Mechanics 33 ( 2003 ) 61 67. [ 8 ] Jorg Hohe Wilfried Becker, Effective stress strain relations for two dimensional cellular sandwich cores: Homogenizat ions, material models and properties, Appl. Mech. Rev. 55 ( 2002 ) 61 87. [ 9 ] J.D. Eshelby, The determination of the field of an ellipsoidal inclusion and related problems, Proc. R. Soc. Lond. A 241 ( 1957 ) 376 396. [ 10 ] Z. Hashin, The elastic moduli of heterogeneous materials, J. Appl. Mech. 29 ( 1962 ) 143 150. [ 11 ] T. Mori K. Tanaka, Average stress in the matrix and average elastic energy of materials with misfitting inclusions, Acta Metall. 21 ( 1973 ) 571 574. [ 12 ] W Voigt, On the relation between the elasticity constants of isotropic bodies, Ann. Phys. Chem. 274 ( 1889 ) 573 587. [ 13 ] A Reuss, Determination of the yield point of Polycrystals based on the yield condition of single crystals, Z. Angew. Math. Mech. 9 ( 1929 ) 49 58. [ 14 ] R Hill, A s elf consistent mechanics of composite materials, J. Mech. Phys. Solids 13 ( 1965 ) 357 372. [15 ] R. M. Christensen K.H. Lo, Solution for effective shear properties in three phase sphere and cylinder models, J. Mech. Phys. Solids 27 ( 1979 ) 315 330. PAGE 101 101 [ 16 ] R. Hill, The elastic behaviour of a crystalline aggregate, Proc. R. Soc. London, Ser. A 65 ( 1952 ) 349 354. [17 ] A. Benssousan, J.L. Lions G. Papanicoulau, Asymptotic a nalysis for p eriodic s tructures, North Holland, 1978. [18 ] Sanchez Palencia, Non hom ogeneous media and vibration theory, Lecture notes in physics127, Springer Verlag, Berlin, 1980. [19 ] J. M. Guedes N. Kikuchi, Preprocessing and post processing for materials based on the homogenization method with adaptive finite element methods, Compu ter Methods in Applied Mechanics and Engineering 83 (1990) 143 198. [20 ] K. Terada, N. K ikuchi, Nonlinear homogenization method for practical applications, Computational Methods in Micromechanics, ASME, New York, AMD 212/MD 62 ( 1995 ) 1 16. [21 ] J. Fish K. Shek, M. Pandheeradi, M.S. Shephard, Computational p lasticity for c omposite s tructures b ased on m athematical homogenization : Theory and p ractice, Comp. Meth. Appl. Mech.Engng. 148 ( 1997 ) 53 73. [22 ] J. Fish and K. L. Shek, Finite d eformation p lastic ity of c omposite s tructures: Computational m odels and a daptive s trategies, Comp. Meth. Appl. Mech. Engng. 172 ( 1999 ) 145 174. [23 ] J. Fish Q.Yu, Multiscale d amage m odeling for c omposite m aterials: Theory and computational f ramework, International Journa l for Numerical Methods in Engineering 52 ( 1 2 ) ( 2001 ) 161 192. [24 ] Terada, K., Kikuchi, N., A class of general algorithms for multi scale analysis of heterogeneous media, Comput. Methods Appl. Mech. Eng. 190 ( 2001 ) 5427 5464. [25 ] K. Matsui, K. Terada K. Yuge, Two scale finite element analysis of heterogeneous solids with periodic microstructures, Computers and Structures 82 ( 2004 ) 593 606. [ 26 ] R.J.M. Smit, W.A.M. Brekelmans, H.E.H. Meijer, Prediction of the mechanical behavior of nonlinear heterog eneous systems by multilevel finite element modeling, Comput. Methods Appl. Mech. Eng. 155 ( 1998 ) 181 192. [ 27 ] C. Miehe, A. Koch, Computational micro to macro transition of discretized microstructures undergoing small strain, Arch. Appl. Mech. 72 ( 2002 ) 300 317. [ 28 ] V. Kouznetsova, W.A. Brekelmans, F.P.T. B aaijens, An approach to micro macro modeling of heterogeneous materials, Comput. Mech. 27 ( 2001 ) 37 48. [ 29 ] F. Feyel, J.L. Chaboche, FE2 multiscale approach for modeling the elasto viscoplastic behavior of long fiber sic/ti composite materials, Comput. Methods Appl. Mech. Eng. 1 83 ( 2000 ) 309 330. PAGE 102 102 [ 30 ] S. Ghosh, K. Lee, S. Moorthy, Multiple scale analysis of heterogeneous elastic structures using homogenization theory and V oronoi cell finite e lement method, Int. J. Solids Struct. 32 ( 1995 ) 27 62. [ 31 ] S. Ghosh, K. Lee, S. Moorthy, Two scale analysis of heterogeneous elasticplastic materials with asymptotic homogenization and Voronoi cell finite element model. Comput. Methods Appl. Mech. Engrg 132 ( 1996 ) 63 116. [ 32 ] J.C. Michel, H. Moulinec, P. Suquet, Effective properties of composite materials with periodic microstructure: a computational approach, Comput. Methods Appl. Mech. Eng. 172 ( 1999 ) 109 143. [ 33 ] M.G.D. Geers, V. Kouznetsova, W .A.M. Brekelmans, Gradient enhanced computational homogenization for the micro macro scale transition, J. Phys. IV 11 ( 2001 ) 145 152. [ 34 ] C. McVeigh, F. Vernerey, W. K. Liu, L. C. Brinson, Multiresolution analysis for material design, Comput. Methods Ap pl. Mech. Engrg. 195 ( 2006 ) 5053 5076. [ 35 ] Derek Hulk T.W, Clyne, An Introduction to composite materials, Cambridge University Press 1996. [36] P. Sharma S. Ganti Size dependent E shelby's tensor for embedded nano inclusions incorporating surface/interface energies Journal of Applied Mechanics 71(5) (2004) 663 671 [37] S.L i, R.A. Sauer, G. Wang, A circular inclusions in a finite domain I. The Dirichlet Eshelby problem. Acta Mech 179 ( 2005 ) 67 90 [ 38 ] R. Hill, Elastic p roperties of r einforced s olids: Some t he oretical p r inciples, J. Mech. Phys. Solids 2 ( 1963 ) 357 372. [39 ] A A. Gusev, Representative v olume e lement size for e lastic c omposites: A numerical study, Mech. Phys. Solids 45 ( 1997 ) 1449 1459. [40] S. Hazanov M. Amieur, On overall properties fo r e las t ic heterogeneous bodies smaller than the representative volu m e, Int. J. Engnr Sci. 33 (9) (1995) 1289 1301 [ 41 ] B.R. Kima H.K. Lee, An RVE based micromechanical analysis of fiber reinforced composites considering fiber size dependency, Composite Str uctures 90 (4) ( 2009 ) 418 427 [42] M.W. Hyer A.M. Waas Micromechanics of l inear e lastic c ontinuous f iber c omposites, Reference Works, Comprehensive Composite Materials, Elsevier 2000 [43] Z. Hashin BW. Rosen, The elastic moduli of fibre reinforce d materials, ASME, Journal of Applied Mechanics 31 (1964) 223 2 32. [44] G. Caligiana, F. Cesari I m ateriali c ompositi, Pitagora Editrice Bologna, 2002 [45] S. Li On the unit cell for micromechanical analysis of fibre reinforced composites Procee dings Royal Society London A 455 (1999) 815 38. PAGE 103 103 [46] D.F. Adams, D.A. Crane, Finite element micromechanical analysis of a unidirectional composite including longitudinal shear loading, Computers and Structures 18 (1984) 1153 65 [47] D.F. Adams, D.R. Doner, Longitudinal shear loading of a unidirectional composite, Journal of Compo s ite Materials 1 (1967) 4 17. [48] D.S. Li, M.R. Wisnom, Finite element micromechanical modeling of unidirectional fibre reinforced metal matrix composites Compos S ci Technol 51 (1994) 545 63 [49] M.R. Nedele, M.R. Wisnom, Finite element micromechanical modeling of a unidirectional composite subject to axial shear loading, Composites 25 (1994) 263 2 72 [50] S. Li, General unit cell for micromechanical analyse s of unidirectional composites, Composites Part A 32 (2000) 815 816 [51] Marek Jerzy Pindera, Jacob Aboudi, Steven M. Arnold, Limitations of the uncoupled, RVE based micromechanical approach in the analysis of functionally graded composites, Mechanics o f Materials 20 (1995) 77 94 [52] R. Hill, Elastic properties of reinforced solids: Some theoretical principles, J. Mech. Phys. Solids 11 (1963) 357 372 [53] Satish Bapananalli Ba Nghiep Nguyen, Prediction of elastic properties for curved fiber polym er composites, Polymer Composites ( 2008 ) 544 550 [54] Zihui Xia, Yunfa Zhang, Fernand Elly, A unified periodical boundary conditions for representative volume elements of composites and applications, International Journal of Solids and Structures 40 (20 03) 1907 1921 [55] G. Mathan, N. Siva Prasad, Evaluation of effective material properties of spiral wound gasket through homogenization, International Journal of Pressure Vessels and Piping 87 (2010) 704 713 [56] Yi Pan, Lucian Iorga Assimina A. Pel egri, Analysis of 3D random chopped fiber reinforced composites using FEM and random sequential adsorption, Computational Materials Science 43 (2008) 450 461 [57] Anthony Drago Marek Jerzy Pindera, Micro macro mechanical analysis of heterogeneous mater ials: Macroscopically homogeneous vs periodic microstructures, Composites Science and Technology 67 (2007) 1243 1263 [58] L. Capolungo C. Jochum Homogenization method for strength and inelastic behavior of nanocrystalline materials, International Jou rnal of Plasticity 21 (2005) 67 82 [59] P. Luciano Moreira Gerard Ferron, Finite element implementation of an orthotropic plasticity model for sheet metal forming simulations, Latin American Journal of Solids and Structures 4 (2007) 149 176 PAGE 104 104 [60] G. Mathan N. Siva Prasad, Evaluation of effective material properties of spiral wound gasket through homogenization, International Journal of Pressure Vessels and Piping 87 (2010) 704 713 [61] Brett A. Bednarcyk Marek Jerzy Pindera, Inelastic response of a woven carbon/copper composite PartII: Micromechanics m odel, Journal of composite materials 34 ( 2000 ) 299 331 [62] Zhenyu Xue Ashkan Vaziri, Non uniform hardening constitutive model for compressible orthotropic materials with application to sandwich plate cores, CMES 1 0 (1) (2005) 79 95 [63] Ji Hoon Kim Nyung Gyu Lee, Development of nonlinear constitutive laws for anisotropic and asymmetric fiber reinforced composites, Polymer composites ( 2008 ) 216 228 [64] O.H. Griffin Jr M.P. Kamat, Three D imensional Inelastic Finite Element Analysis of Laminated Composites, Journal of Composite Materials 15 ( 1981 ) 543 560. [65] Petri Makela, Soren Ostlund Orthotropic elastic plastic material model for paper materials, International Journal of Solids and Structures 40 (2003) 5599 5620 [66] Roland Mucke Otto Ernes Bernhardi, A constitutive model for anisotropic materials based 192 (2003) 4237 4255 [67 ] Ala Tabiei Jin Wu, Three dimensional nonlinear orthotropic finite element material model for wood, Composite Structures 50 (2000) 143 149 [ 68 ] Nicholas Metropolis S. Ulam, The M onte Carlo m ethod, Journal of the American Statistical Association 44 ( 1949 ) 335 341. [ 69 ] G. Amaniampong, C.J. Burgoyn e, Monte Carlo s imulations of the time dependent failure of bundles of parallel fibers, Eur. J. Mech. A/Solids 15(2) ( 1996 ) 243 266. [ 70 ] M.E. Cruz, A.T. Patera, A parallel Monte Carlo finite element procedure for the analysis of multicomponent media. I nt. J. Num Meth. Engrg. 38 ( 1995 ) 1087 1121. [71 ] R.L. Iman, W.J. Conover, Small sample sensitivity analysis techniques for computer models, with an application to risk assessment. Communications in Statistics, Part A. Theory and Methods ( 1980 ) 1749 1842 [72] R.E. Melchers, Structural Reliability Analysis and Prediction, John Wiley & Sons, New York 2001. [73] S.S. Isukapalli, A. Roy, P.G. Georgopoulos, Stochastic r esponse s urface m ethods(SRSMs) for u ncertainty p ropagation: Application to e nvironm ental and b iological s ystems, Risk Analysis 18(3) ( 1998 ) 351 363 [74] Shuping Huang, Sankaran Mahadevan, Ramesh Rebbaa Collocation based stochastic finite element analysis for random field problems, Probabilistic engineering mechanics 22 (2007) 194 205 PAGE 105 105 [ 75 ] S. S. Isukapalli, A. Roy, P. G. Georgopoulos, Efficient s ensitivity/ u ncertainty a nalysis u sing the c ombined s tochastic r esponse s urface m ethod and a utomated d ifferentiation: Application to e nvironmental and b iological s ystems, Risk Analysis 20 ( 5 ) ( 2000 ) 591 602. [ 76 ] A.A. Henriques, J.M.C. Veiga, J.A.C. Matos, J.M. Delgado, Uncertainty analysis of structural systems by perturbation techniques, Struct. Multidisc. Optim. 35 (2008) 201 212. [ 77 ] S. Sakata, F. Ashida M. Zako, Kriging based approxi mate stochastic homogenization analysis for composite materials, Computa. Methods Appl. Mech. Engrg. 197 (2008) 1953 1964. [78 ] Qing Sheng Yang Qing Hua Qin, Micro mechanical analysis of composite materials by BEM, Engineering Analysis with Boundary Ele ments 28 (2004) 919 926. [79] A.A. Henriques, Efficient analysis of structural uncertainty using perturbation techniques, Engineering Structures 30 (2007) 990 1001. [80] Z. Hashin S. Shtrikman, A variational approach to the theory of the elastic behav ior of polycrystals, J.Mech.Phys.Solids 10 ( 1962 ) 343 352 [81] Anthony Drago, Marek Jerzy Pindera, Micor macromechanical analysis of heterogeneous materials: Macroscopically homogeneous vs periodic microstructure, Composite science and technology 67 (20 07) 1243 1263 [8 2 ] S. Sakata, F. Ashida, T. Kojima, M. Zako Three dimensional stochastic analysis using a perturbationbased homogenization method for elastic properties of composite material considering microscopic uncertainty, International Journal of Solids and Structures 45 (2008) 894 907 [8 3 ] S. Sakata, F. Ashida M. Zako, Kriging based approximate stochastic homogenization analysis for composite materials, Comput. Methods Appl. Mech. Engrg. 197 (2008) 1953 1964 [84] Roger Ghanem John Red Hor se, Propagation of probabilistic uncertainty in complex physical systems using a stochastic finite element approach, Physica D 133 (1999) 137 144 [85] D. R. Millman, P. I. King Estimating the probability of failure of a nonlinear aeroelastic system, Jo urnal of Aircraft 43 (2006) 504 516. [86] Nicolas Lamorte, Bryan Glaz, Uncertainty p ropagation in h ypersonic a erothermoelastic a nalysis, 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural d ynamics, and Materials c onference,12 15 April 2010, Orlando, Flor ida [87] Hamed Haddad Khodaparast John E. Mottershead Kenneth J. Badcock Propagation of structural uncertainty to linear aeroelastic stability, Computers and Structures 88 (2010) 223 236 [88] Xiaolei Yin, Wei Chen, Albert To, Cahal McVeigh, Wing Ka m Liu, Statistical volume element method for predicting microstructure constitutive property relations, Comput. Methods Appl. Mech. Engrg. 197 (2008) 3516 3529 PAGE 106 106 [ 89 ] Boris Jeremic, Kallol Sett, Probabilistic elasto plasticity: formulation in 1D, Acta Geo technica 2 (2007) 197 210 [9 0 ] Maciej Anders Muneo Hori, Stochastic finite element method for elasto plastic body, Int. J. Numer. Meth. engng. 46 (1999) 1897 1916 [9 1 ] Maciej Anders Muneo Hori, Three dimensional stochastic fi nite element method for elasto plastic bodies, Int. J. Numer. Meth. Engng 51 ( 2001 ) 449 478 [9 2 ] R.B. McKee, Jr. George Sines, A statistical model for the tensile fracture of parallel fiber composites, J. ELASTOPLASTICS 1 ( 1969 ) 185 199 [9 3 ] Jack R. Vison Robert L. Sier akowski, The behavior of structures composed of composite materials, Kluwer academic publishers, NY, 2004 [9 4 ] Alan Baker Stuart Dutton and Donald Kelly, Composite materials for aircraft structures, second edition, AIAA, 2004 [95] William D. Callist er, Jr., Materials science and engineering an introduction, John Wiley & Sons, Inc.,2003. [9 6 ] Debdatta Ratna, epoxy composites: Impact r esistance and f lame r etardancy: Rapra Revie w Reports 185 16 (5), 2007. [9 7 ] Derek Hull T.W. Clyne An introductio n to composite materials, second edition, Cambridge University Press, 1996 [ 9 8 ] Alan Baker, Stuart Dutton Donald Kelly, Composite materials for aircraft structures, second edition, AIAA, 2004 [ 99 ] Deborah D.L. Chung, Carbon fiber composite, Butterwo rth Heinemann, 1994 [10 0 [10 1 ] Bodo Fiedler Thomas Hobbierunken, Influence of stress state and temperature on the strength of epoxy resins, 11 th International conference on fracture, Tur in, Italy, March 20 25, 2005 [10 2 ] Jacob Lubliner Plasticity t heory, Dover publication, 2008 [10 3 ] R. von Mises, Mechanik der festen Krper im plastisch deformablen Zustand. Gttin. Nachr. Math. Phys. 1 (1913) 582 592 [10 4 ] R. Hill A theory of the yie lding and plastic flow of anisotropic metals. Proc. Roy. Soc. London 193 (1948) 281 297 [10 5 ] R. Hill Theoretical plasticity of textured aggregates. Proceedings of the Cambridge Philosophical Society 75 (1979) 179 208 [10 6 ] R. Hill Constitutive modelin g of orthotropic plasticity in sheet metals. Journal of the Mechanics and Physics of Solids 38(3) (1990) 405 417 PAGE 107 107 [10 7 ] R. Hill Plastic anisotropy and the geometry of yield surfaces in stress space, Journal of the Mechanics and Physics of Solids 48 (2000) 1093 1106 [10 8 ] Walter Ramberg, Willam R. Osgood, Description of stress strain curves by three parameters, National Advisory Committee for Aeronautics, Technical Note 902,1943. [1 09 ] Sinisa Dj Mesarovic Jagan Padbidri, Minimal kinematic boundary conditi ons for simulations of disordered microstructures, Philosophical magazine 85 (1) January 2005, 65 78 [ 11 0 ] S.S. Isukapalli, Uncertainty a nalysis of t ransport t ransformation m odels. PhD Thesis, Rutgers University, Piscataway, NJ. 1999. [ 11 1 ] R. H. Myers D. C. Montgomery, Response s urface m ethodology Wiley New York 1995 [11 2 ] W. Gautschi, Orthogonal Polynomials: Applications and c omputations. Acta Numerica 5 (1996) 45 119. [11 3 ] D. Xiu, D. Lucor, C.H. Su, G.E. Karniadakis, Stochastic m odeling of f low s tr ucture i nteractions u sing g eneralized p olynomial c haos Journal of Fluids Engineering 124(1) (2002) 51 59. [11 4 ] Sastry S. Isukapulli, Panos G.Georgopoulos, Computationally efficient uncertainty propagation and reduction using the stochastic response surfa ce method, 43 rd IEEE Conference on Decision and Control, December 14 17, 2004 [11 5 ] A.A. Henriques, J.M.C. Veiga, J.A.C. Matos, J.M. Delgado, Uncertainty analysis of structural systems by perturbation techniques, Struct. Multidisc. Optim. 35 (2008) 201 21 2. [ 11 6 ] M.D. Webster, M.A. Tatang, G.J. McRae, Application of the p robabilistic c ollocation m ethod for an u ncertainty a nalysis of a s imple o cean m odel, MIT Joint Program on the Science and Policy of Global Change, 1996. [11 7 ] J. Villadsen, M.L. Michelsen, Solution of differential equation models by polynomial approximation, Prentice Hall, Englewood Cliffs, New Jersey, 1978. PAGE 108 108 BIOGRAPHICAL SKETCH Jinuk Kim was born in Pusan, South Korea (Republic of Korea) in 1972. He received a Bachelor of Science a nd Master of Science in m echanical e ngineering from Sungkyunkwan University in Korea in March 2000. He worked as an intern research er in SAFE research center in Sungkyunkwan University until 2002 studying on Tribology. At the University of Michigan he earned a Master of Science in mechanical engineering and a Master of Science in mat erials science and engineering i n July 2006. For further study, he went to the University of Florida to pursu e a Ph.D. degree in materials science and engineering where he earn ed his doctorate in August 2011 