MECHANICAL BEHAVIOR OF WOVEN FABRIC COMPOSITES UNDER MESO-LEVEL UNCERTAINTIES: MODELING AND SENSITIVITY ANALYSIS by Mojtaba Komeili B.Sc. Iran University of Science and Technology, Iran, 2008 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The College of Graduate Studies (Mechanical Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Okanagan) August 2010 © Mojtaba Komeili, 2010 ii Abstract Unit cell modeling of the woven fabric composites is a strong tool for studying fabric behavior at meso-level. Among the suggested methods the three-dimensional finite element analysis is found to be promising and of greater interest to researchers in the field. However, due to the multi-scale nature and particular behavior of fibrous yarns, numerical procedures applicable to woven fabrics differ from conventional finite element routines. Moreover, from a practical point of view, most of the models in the literature focus on the ideal unit cells and excludes the meso-level uncertainty factors. On the other hand, experimental measurements indicate non-repeatability of test data which is a result of embedded inherent uncertainties. In this study, a finite element model capable of defining the representative volume element of a plain weave under homogenous loading has been created. Because of particular behavior of dry fabrics, a special constitutive material behaviour is defined via user-defined subroutines, which links to the solver of a finite element package (Abaqus). Material properties of yarns for dry glass fabrics are extracted by fitting their numerical responses to results of experiments. Subsequently, geometrical and material meso-level uncertainty factors are studied separately using two-level factorial designs. The output of runs are inserted into a commercial design package (Design Expert) for producing results in terms of probability plots, effects and percentage contribution of each factor and their interactions. Result shows that depending on the loading type, there are factors that show significant contribution toward the final response, whereas others are quite negligible. More elaboration on the design results provides informative conclusions regarding modeling of unit cells, as well as their behavior during different loading steps. These conclusions can be used in more comprehensive unit-cell homogenization formulations, as well as for defining appropriate tolerances for meso-level defects in woven fabric composites. iii Table of Contents Abstract ....................................................................................................................................................... ii Table of Contents ....................................................................................................................................... iii List of Tables .............................................................................................................................................. vi List of Figures ........................................................................................................................................... vii Acknowledgements ......................................................................................................................................x Dedication ................................................................................................................................................... xi Chapter 1 : Introduction ............................................................................................................................1 1.1 Textile composites and woven fabrics .................................................................................................1 1.2 Multi-scale modeling of woven fabrics ...............................................................................................5 1.3 Meso-level modeling ...........................................................................................................................7 1.4 Motivation and objectives of the work .............................................................................................. 17 1.5 Thesis outline .................................................................................................................................... 18 Chapter 2 : FEM modeling of fabric unit cell at meso-level ................................................................. 20 2.1 Introduction ....................................................................................................................................... 20 2.2 Standard procedures for 3D FE modeling of unit cells ..................................................................... 21 2.3 Geometry of unit cell ......................................................................................................................... 22 2.4 Material modeling ............................................................................................................................. 25 2.4.1 Micro-scale homogenization ...................................................................................................... 25 2.4.2 Constitutive model based on physical observations ................................................................... 30 2.5 Periodicity boundary conditions ........................................................................................................ 30 2.6 Loading boundary condition .............................................................................................................. 32 2.7 Nonlinearity in analysis ..................................................................................................................... 34 2.8 Finite element solver ......................................................................................................................... 36 2.9 Constitutive model for dry fabrics ..................................................................................................... 38 2.10 Summary ......................................................................................................................................... 40 Chapter 3 : Creating FE model in Abaqus ............................................................................................ 41 3.1 Reference unit cell ............................................................................................................................. 41 3.2 CAD modeling of yarns ..................................................................................................................... 42 3.3 Defining model in Abaqus FE software ............................................................................................ 44 3.3.1 Kinematic conditions for axial loading ...................................................................................... 45 3.3.2 Kinematic conditions for shear loading ...................................................................................... 47 3.4 Material modeling ............................................................................................................................. 50 iv 3.4.1 Homogenized material properties ............................................................................................... 51 3.4.2 Postulating a suitable constitutive model in meso-level ............................................................. 53 3.5 Implementing the material constitutive model in Abaqus ................................................................. 55 3.5.1 Explicit integrator ....................................................................................................................... 55 3.5.2 Explicit subroutine...................................................................................................................... 56 3.5.3 Implicit integrator ....................................................................................................................... 59 3.5.4 Implicit subroutine...................................................................................................................... 60 3.6 Postulating material model constants ................................................................................................ 64 3.6.1 Material model ........................................................................................................................... 64 3.6.2 Longitudinal yarn behavior ........................................................................................................ 65 3.6.3 Transversal yarn behavior .......................................................................................................... 65 3.7 Inverse identification ......................................................................................................................... 66 3.8 Summary ........................................................................................................................................... 71 Chapter 4 : Uncertainty factors and their modeling ............................................................................. 72 4.1 Introduction ....................................................................................................................................... 72 4.2 Uncertainty of test results .................................................................................................................. 72 4.3 Geometrical uncertainty factors......................................................................................................... 74 4.4 Material uncertainty factors ............................................................................................................... 77 4.5 Summary ........................................................................................................................................... 79 Chapter 5 : Sensitivity analysis on uncertainty factors ......................................................................... 80 5.1 Base model and its validation ............................................................................................................ 80 5.2 Factorial design: Fundamentals ......................................................................................................... 83 5.3 Planning factorial design ................................................................................................................... 85 5.3.1 Choosing response function ....................................................................................................... 85 5.4 Geometrical sensitivity analysis ........................................................................................................ 87 5.4.1 Axial loading .............................................................................................................................. 88 5.4.2 Shear loading .............................................................................................................................. 93 5.5 Material sensitivity analysis .............................................................................................................. 97 5.6 Non-repeatability of test results ....................................................................................................... 103 5.7 Some more discussions on the results ............................................................................................. 105 5.7.1 Important effect of misalignment in axial loading ................................................................... 105 5.7.2 Effect of fricion in axial loading .............................................................................................. 106 5.8 Summary ......................................................................................................................................... 108 v Chapter 6 : Conclusions and remarks .................................................................................................. 109 6.1 Summary of methodology ............................................................................................................... 109 6.2 Summary of inverse method and sensitivity analysis results .......................................................... 110 6.2.1 Geometrical sensitivity ............................................................................................................. 110 6.2.2 Material sensitivity ................................................................................................................... 112 6.3 Recommended future work ............................................................................................................. 112 Bibliography ............................................................................................................................................. 114 Appendix A: Transformation of second order tensors in vector form ............................................... 123 Appendix B: Matrix form of square root .............................................................................................. 125 vi List of Tables Table 3-1: Material properties of E-Glass/PP yarns predicted by Peng and Cao (2002) ................... 51 Table 4-1: Geometrical uncertainty factors ............................................................................................ 76 Table 4-2: Material uncertainty factors .................................................................................................. 78 Table 5-1: Symbols for factors used in geometrical sensitivity analysis ............................................... 87 Table 5-2: The value of each factor in geometrical sensitivity analysis ................................................ 88 Table 5-3: The values of response indicators in each run relative to those of the base unit cell (axial loading) ....................................................................................................................................................... 90 Table 5-4: Results of the 2-level material factorial design for shear tests; A (yarn spacing), B (yarn width), C (yarn height) and D (misalignment) ........................................................................................ 91 Table 5-5: The difference between values of different indicators in each run using the basic unit cell (shear loading) ........................................................................................................................................... 94 Table 5-6: Results of the 2-level geometrical factorial design for shear tests; A (yarn spacing), B (yarn width), C (yarn height) and D (misalignment). ............................................................................ 95 Table 5-7: Symbols for factors used in material sensitivity factorial analysis ..................................... 97 Table 5-8: The value of each factor for material sensitivity analysis .................................................... 97 Table 5-9: The values of response indicators in each run relative to those of the base unit cell (shear loading) ....................................................................................................................................................... 99 Table 5-10: Results of the 2-level factorial design for geometrical shear tests; A (longitudinal stiffness), B (transverse stiffness) and C (friction coefficient) ............................................................. 101 vii List of Figures Figure 1-1: Four common types of textiles (a) woven (b) braided (c) knitted (d) stitched fabric (Potluri and Manan 2007) ...........................................................................................................................2 Figure 1-2: Examples of woven fabrics (a) plain ng=2 (b) twill ng=3 (c) 4 harness satin ng=4. ............3 Figure 1-3: Multi-scale nature of woven fabrics .......................................................................................6 Figure 1-4: Simplified geometry used by Kawabata et al. (1973 a-c) for study of woven fabrics ........7 Figure 1-5: Geometrical assumption of woven fabrics used in Ishikawa methods (a) mosaic model (b) fibre undulation model (c) bridge model (Ishikawa and Chou 1982) ...............................................8 Figure 1-6: 3D finite element models used by Boisse et al. (2001) ......................................................... 12 Figure 1-7: Discrete Element Method; A: warp, B: weft, C: transverse yarn compression, D: Rotational Resistance, mi: concentrated mass. ....................................................................................... 14 Figure 1-8: Geometry and mesh of the model for a plain weave fabric including all fibres in each yarn and their interactions (Durville 2005) ............................................................................................. 15 Figure 1-9: Interpenetration between the yarns (Hivet and Boisse 2008) ............................................ 16 Figure 2-1: Balanced plain weave fabric (Badel et al. 2008a) ............................................................... 20 Figure 2-2: The unit cell geometry used by Gasser (2000); (a) an elementary pattern (b) the quarter of elementary pattern. ............................................................................................................................... 23 Figure 2-3: Two unit cell types used by Badel et al. (2007) for modeling of shear loading ................ 23 Figure 2-4: Unit cells used by Chen et al. (2001) a) Basic unit cell b) Extended unit cell. .................. 24 Figure 2-5: Unit cell for a micro-structure (Peng and Cao, 2002) ........................................................ 26 Figure 2-6: Initial and deformed meso-structure (Boisse et al. 2007) ................................................... 31 Figure 2-7: Boundary condition on a unit cell (Boisse et al. 2007) ........................................................ 32 Figure 2-8: Displacement boundary conditions for different types of loading (a) uni-axial tension (b) bi-axial tension (c) shear ........................................................................................................................... 33 Figure 2-9: A sample picture frame test setup (Cao et al. 2008) ........................................................... 34 Figure 2-10: Undulation in the yarn (Boisse et al. 1997) ........................................................................ 35 Figure 2-11: Drifting error in explicit finite element (Felippa 2001) .................................................... 37 Figure 2-12: Implicit finite element procedure, prediction followed by correction (Felippa 2001). .. 38 Figure 3-1: Geometrical parameters for defining the basic yarn. ......................................................... 42 Figure 3-2: Dialog box for defining the geometrical parameters of the yarn ...................................... 43 Figure 3-3: Defining the parametric cross section. ................................................................................. 43 Figure 3-4: Creating solid yarn using two end cross sections and the connecting line. ...................... 44 viii Figure 3-5: Node sets created for defining periodic boundary condition on one of the unit cell surfaces. ...................................................................................................................................................... 46 Figure 3-6: Defining the constraint for surrounding surfaces in Abaqus/CAE ................................... 47 Figure 3-7: Rigid frame and the coordinate system attached to each link before and after deformation. ............................................................................................................................................... 48 Figure 3-8: (a) Pinned nodes at the corner of the unit cell. (b) Node sets and reference points for keeping the corner lines perpendicular to the deformation plane. ....................................................... 49 Figure 3-9: Coupling of yarn (a) side surfce (b) end surfaces. .............................................................. 50 Figure 3-10: Corners of the rigid frame constrained to the movement of a reference point. ............. 50 Figure 3-11: Validation of the refernece finite element model .............................................................. 52 Figure 3-12: Comparison of Fibre Frame and other possible frame of rotation in 2D; (a) before deformation frames coincide, (b) the separation of frames after deformation. ................................... 53 Figure 3-13: Flowchart for programming the user-defined subroutine in explicit solver .................. 57 Figure 3-14: Flow chart for programming the user-defined subroutine in implicit solver ............... 63 Figure 3-15: Photo from sample tested by Buet-Gautier and Boisse (2001) ........................................ 66 Figure 3-16: Local coordinate system for defining material orientation along yarn .......................... 67 Figure 3-17: Response force of fitted model for uni-axial test data ...................................................... 68 Figure 3-18: Response force of fitted model for equi-biaxial test data ................................................. 69 Figure 3-19: Response force of fitted model for unbalanced bi-axial test data.................................... 69 Figure 3-20: Results of response curves from fitted model vs. two experiments from (Cao et al. 2001) ............................................................................................................................................................ 70 Figure 4-1: A tomographic image of a woven fabric and its section (Badel et al. 2009) ..................... 73 Figure 4-2: Results of shear tests on similar woven fabric samples by different laboratories (Cao et al. 2008)....................................................................................................................................................... 74 Figure 4-3: Misalignment in a typical woven fabric (Milani et al. 2007) .............................................. 75 Figure 4-4: The plain woven fabric used by Cao et al. (2008) ............................................................... 76 Figure 4-5: Positive misalignment for axial loading ............................................................................... 77 Figure 4-6 : Negative (a) and positive (b) misalignments in shear test ................................................. 77 Figure 5-1: Sensitivity of basic model to mesh size and number of sections under axial loading ...... 81 Figure 5-2: Diagonal length of the shear frame used for calculating stretch ratio (before and after deformation) ............................................................................................................................................... 82 Figure 5-3: Sensitivity of basic model to mesh size and number of sections under shear loading ..... 82 ix Figure 5-4: Imaginary frame defined for shear frame test under misalignment, (a) positive misalignment and (b) negative misalignment ......................................................................................... 86 Figure 5-5: Schematic of two indicators used for the assessment of case studies ................................ 87 Figure 5-6: Result of factorial design runs for axial mode (the corresponding values of indicators for each run are given in Table 5-3) ......................................................................................................... 89 Figure 5-7: Half normal probability plots for geometrical factorial design (axial loading) ............... 91 Figure 5-8: Effects of main factors in their normalized range .............................................................. 92 Figure 5-9: Effect of AD interaction in normalized range of A (uni-axial mode) ................................ 93 Figure 5-10: Result of factorial design runs for shear mode (the corresponding values of indicators for each run are given in Table 5-5) ......................................................................................................... 94 Figure 5-11: Half normal probability plots for geometrical factorial design (shear loading) ............ 95 Figure 5-12: Effects of important main factors in their normalized range .......................................... 96 Figure 5-13: Effect of BD interaction in the normalized range of B ..................................................... 96 Figure 5-14: Result of factorial design runs (the corresponding values of indicators for each run are given in Table 5-9) ..................................................................................................................................... 98 Figure 5-15: Half-Normal probability plots of material factorial design ........................................... 100 Figure 5-16: Effect of main factors in material factorial design ......................................................... 102 Figure 5-17: Non-repeatability of test data and comparing it to the factorial design results. Red line: nominal unit cell; Grey lines: factorial designs; Green lines: experiment ................................. 104 Figure 5-18: Effect of shear stiffness on the importance of misalignment in axial loading .............. 106 Figure 5-19: Superposition of deformed and un-deformed unit cell in bi-axial loading. Solid shape is the deformed and shadow shape is the un-deformed unit cell............................................................. 107 Figure 5-20: Effect of friction and longitudinal stiffness in larger strain values; A (longitudinal stiffness) and B (friction coefficient) ...................................................................................................... 108 x Acknowledgements First of all I want to show my gratitude to almighty God for blessing me with the strength and opportunity to study as a graduate student at this level and everything that I have in my life. This research was done under supervision of Dr. Abbas S. Milani. I wish to express my great appreciation for the patient help, encouragement and excellent inspiration I have received from him. He was not only a supervisor, but also a very supportive friend of mine during my study in Canada. I also want to thank all my friends for their kindness and keeping me entertained during my stay in Kelowna. The financial support of the Natural Sciences and Engineering Research Council of Canada (NSERC) is greatly acknowledged. I also want to thank my thesis examination committee members (Drs. A.S. Milani, A. Rteil, S. Tesfamariam, H. Najjaran from UBC Okanagan and Dr. X. Peng from Shanghai Jiao Tong University) and also Prof. J. Cao from Northwestern University for their valuable time and advice. At last but not least I want to express my gratefulness to my family and, especially, my parents for their unconditional love at all times; for this I am eternally indebted. I also want to thank the Baird family who didn’t let me feel lonely away from home. xi Dedication To my parents for their unlimited love and care! 1 Chapter 1: Introduction In recent years, there has been an increase in the use of composite materials, especially in the aerospace and transportation industries. In comparison to traditional materials, composites offer higher strength to weight ratios, non-corrosive properties, dimensional stability and good conformability. Their high strength to weight ratio and anisotropic properties are particularly used in structural design and optimization. Recently, percentage weight of composites in aerospace applications has increased from 3% in Airbus A320 to over 20% in A380. Similarly, Boeing’s 787 uses over 50% of composites (Potluri and Manan 2007). In naval architecture, Fibre-Reinforced Polymer (FRP) composites have been increasingly used in the production of hulls, topsides, masts, sails and other parts of the vessels. Furthermore, composite materials have played key roles in reducing the magnetic, acoustic, hydrodynamic, radar, thermal signatures, as well as increasing payload, top speed, and operation range in marine vessels (Lua 2007). Additionally, the Comanche helicopter, the Composite Armored Vehicle, the Resin Fusion Infusion (RFI) NASA (Chung and Tamma 1999), the centre wing box of A380, and wings and fuselage of Boeing and Airbus’ future aircrafts (Boisse et al. 2007) are among the numerous examples of composite material usage in modern industries. 1.1 Textile composites and woven fabrics Textile composites are produced by mechanically bonding two or more fibre bundles (also known as yarns) together in a specific architecture to meet the desired mechanical, thermal, etc. properties in a final product. The four most commonly used textile composites are woven, braided, knitted and stitched fabrics. The difference between these types of woven fabrics arises from the way that yarns are placed and held together. Figure 1-1 shows a schematic for these four different textile types. The most common type of textiles is the woven fabric, which is made through the interlacing of yarns. Braided fabrics are produced by interweaving three or more strands of yarn, whereas knitted fabrics are fabricated by inter-looping yarns in a horizontal or vertical direction. The former is known as weft knit and latter as warp knit, respectively. In stitched fabrics, also called non-crimp fabrics, yarns are stitched to each other. 2 Figure 1-1: Four common types of textiles (a) woven (b) braided (c) knitted (d) stitched fabric (Potluri and Manan 2007) Essentially, all textiles composites are made of two-dimensional or three-dimensional repetitions of a pattern, known as Representative Unit Cell (RUC) or shortly Unit Cell. This pattern can be used for the classification of any type of fabric. For woven fabrics (which are the type of textiles used in this thesis), there are two geometrical parameters that must be defined in order to identify fabric (Ishikawa and Chou 1982) nfg and nwg, where the former denotes a warp thread is interlaced with every nfg-th fill and the later denotes a fill thread interlaced with every nwg-th warp. Figure 1-2 shows three different types of woven fabrics most commonly used in industry. The yellow box shows the common unit cell for each fabric type. Here only the balanced fabrics were considered i.e. the case of nfg= nwg= ng. 3 Figure 1-2: Examples of woven fabrics (a) plain ng=2 (b) twill ng=3 (c) 4 harness satin ng=4. Over the past decade, considerable attention from composite manufacturing sector has been devoted to textile composites and especially woven fabrics. These materials have many advantages over unidirectional fibre reinforce composites, such as enhanced dimensional stability over a large range of temperatures (Tabiei and Yi 2002), more balanced properties in the fabric plane, and better impact resistance. Generally, unidirectional fibre reinforced and laminated composite plates are used where the in-plane mechanical properties of the final part is as important as an overall design factor. However, they have proven to show weaknesses in thickness direction and they are prone to interlaminar delamination (Aitharaju and Averill 1999). To overcome this limitation of composites, textile composite laminates are utilized to provide tri- directional reinforcement. In comparison to unidirectional composites, impact toughness in woven textile composites is higher because of interlacing fibre bundles that prevent the growth of damage (Peng and Cao 2000). In addition, the reorientation of the fibres can cause significant change in the fabric stiffness with high degree of anisotropy, which affects the material behavior in both crash and impact loading (Tabiei and Ivanov 2003). These features make woven fabrics popular for use in structures that are subjected to transverse impact loading, such as airbags, human body armors and protective jackets in airplane jet engines (Ivanov and Tabiei 2002). In addition, the out of plane strength of woven fabrics, which is a result of interlacing yarns, can be useful in carrying secondary loads due to load-path eccentricities, local buckling, etc. (Seng et al. 1997). One of the other important advantages of woven fabrics for composite manufacturers is their ease of formability and low production cost. Woven fabrics are easy to handle in dry or pre- impregnated pre-form and offer good drapeability. Many different forming techniques, such as (a) (c)(b) 4 rolling, molding, compression and diaphragm forming, as well as machining, have been utilized for the manufacturing of woven fabric composites (Lim and Ramakrishna 2002). Reports have also shown that woven composites are particularly suited for manufacturing of doubly curved components (Page and Wang 2002). Furthermore, in relation to increasing concerns regarding environmental aspects of industrial products, thermoplastic woven composites are seen as superior materials for mass production, because of their recycling potential compared to those of thermoset composites (Lebrun et al. 2003). When considering the huge amount of raw material consumption in modern industries, especially in transportation, the use of thermoplastic textile composites can be an efficient way of reducing environmental waste. There are also several new research areas that have brought the capabilities of other technologies, such as flexible electronics, microfluidics and actuated materials, into the woven fabrics field and vice versa. This helps in design and manufacturing of highly advanced products like body armor with embedded medical sensors or communications equipment, apparel with microfluidic cooling or heating capabilities (King et al. 2005). As previously stated, one of the primary advantages in woven fabric composites is their flexibility for making 3D complex surfaces. For example, in their manufacturing stage, dry fabrics are first formed into the 3D shape by punching, manual laying, etc., and then impregnated by resin to give them a solid shape. Deformations and loadings during these processes can affect the mechanical behavior of the final product, as well as the behavior of the material in additional manufacturing processes (if there is more than one process involved). Prediction of material properties in the final product can, however, introduce some difficulties due to the highly nonlinear behavior of woven fabrics. This is a result of the extremely low resistance of fabrics to in-plane shear deformation, which allows for yarns to orient easily and gives fabrics an enormously anisotropic behavior (Xue et al. 2003). There can be the possibility of defects, such as wrinkles, folds and tearing during the forming process of woven fabrics. Moreover, the deformation of fabric changes the permeability of parts that can significantly affect the RTM (Resin Transfer Molding) efficiency. Consequently, there is a substantial need for detailed models that can predict the mechanical behavior of woven fabrics during the manufacturing process, as well as their mechanical properties after the deformation in the final product. The above discussion implies that the first step in utilizing composite materials, especially the woven fabric composites in engineering platform and their spectacular benefits in design, is a 5 reliable prediction of their material properties and mechanical behavior. The development of woven fabric textiles hinges on comprehensive and convenient models that can be used in real life applications. The mechanical properties of woven fabrics, including stiffness, strength and thermal coefficients, are important in optimizing manufacturing of composites. It should also be noted that the established material models must be neither too complicated to make their usage practically feasible, nor too simple to underestimate the accurate properties and characteristics of the composite materials. Additionally, they should be flexible enough to be incorporated into common modeling procedures, such as finite element methods. 1.2 Multi-scale modeling of woven fabrics The mechanical behavior of woven fabrics is difficult to predict due to complex interactions of yarns in the fabric and the interaction of fibres in each yarn. Indeed, the associated problem of characterizing the multiple scales is the greatest obstacle to unrestricted implementation of woven fabrics. While the traditional short-fibre and unidirectional long fibre composites are often treated at micro- and macro-scales (Potluri et al. 2006), textile composites and specifically woven fabrics have an additional intermediate level known as meso-level. Essentially, woven fabrics can be considered as structured, hierarchical materials, having three structural levels (Lomov et al. 2007). These three levels are called macro-level, meso-level and micro-level, where the first has the largest and the latter has the smallest order of magnitude (Figure 1-3). The order of magnitude for each of these scales and the definition of material behavior is as follows: 6 Figure 1-3: Multi-scale nature of woven fabrics Macro-level: Can be defined as the level of geometry of the final composite product. In this scale, which is usually in the order of 10−1- 100 meters, the geometry of mechanical parts including curvature, fibre volume fraction and shear angle are most frequently dealt with. The macroscopic behavior is, however, extremely dependent on the interactions of yarns at the meso- scale. Meso-level: Is characterized by the yarn dimensions, the interactions of individual yarns and fabric textures. Meso-scale defines both the internal structure of the yarns and the variations of the yarn direction. The order of magnitude here is 10−3- 10-2 meters, depending upon the yarn size and pattern of fabric. Indeed, this level greatly influences the mechanical behavior of the woven fabrics and, consequently, is the most necessary to analyze (Guagliano and Riva 2001). Nevertheless, each yarn itself is made of bundles of fibres, which makes the next hierarchical scale. The analysis at meso-scale includes the geometry of the woven fabric and interlacing of the yarns. Here, the yarns are considered as continuous domains. This assumption can, however, be disputed considering the micro-level behavior. Micro-level: Is defined by the arrangement of fibres in each impregnated yarn or fibrous ply. The order of magnitudes in this scale is close to the fibre size (typically around 5–20 μ meters), which is relatively small in comparison to the final part or even fibre yarns. Nonetheless, fibre interactions and behaviors have a significant influence on the properties at the meso-level, which eventually has an effect on the macro-level behavior. Micro‐level fibers ~10‐6 m Macro‐level part ~100 m Meso‐level yarns ~10‐3 m 7 1.3 Meso-level modeling As previously discussed, it is important to study the woven fabrics and their behaviour in the meso-level. Response at this scale can be used in a homogenization process to find the equivalent material properties for models at macro-level, as well as to investigate local deformation mechanisms that are often necessary in damage analysis. Details of the homogenization of material properties for macro level can be found in the study conducted by Takano et al. (1999), Carvelli and Poggi (2001) and Peng and Cao (2002). To study the behavior of the parts made from woven fabrics, it is necessary to use their homogenized material properties at macro-level, which requires a detailed study of fabric behavior at meso-level. Kawabata et al. (1973a-c) presented general form theories for uni-axial, bi-axial and shear deformation properties of plain weave fabrics based on a simplified model representing the structure of the fabric unit cells. The basic geometry used in this theory is based on one- dimensional stiffness elements (bars) representing yarns. Lateral compression between the yarns is also modeled by stiffness elements that are located at intersection of the yarns (Figure 1-4). Mathematical formulation includes some geometrical and mechanical constants, which can be determined using real specimens and experimental data (usually yarn tensile test and fabric bi- axial tension test). They also showed that compressive properties of the yarn, which is a result of lateral compressive force, have a profound influence on the tensile properties of the fabrics. Figure 1-4: Simplified geometry used by Kawabata et al. (1973 a-c) for study of woven fabrics Ishikawa and Chou (1982) compared three different analytical models for the investigation of stiffness and strength of an impregnated unit cell woven fabric. These three methods were 8 already developed by the same authors and included “mosaic model,” “fabric undulation” and “bridge model.” In the mosaic model, the resin impregnated woven fabric was idealized as an assemblage of asymmetrical cross ply laminates and yarn waviness is neglected. By applying iso-stress and iso-strain assumption on this model, lower and upper bounds of elastic modulus can be determined (Ishikawa and Chou 1982). Fabric undulation method (also known as crimp method) was designed for assessing the validity of the mosaic model and for studying the elastic and knee behavior of woven fabrics (Ishikawa and Chou 1983). The bridge model idea was based on the success of the “threadwise” idealization, which was suggested in the undulation model. This method is effective in the study of satin composites; further details about the method are presented in Ishikawa and Chou (1982). It should be noted that all of these methods are based on the assumption of one-dimensional behavior of woven fabrics and classical laminate theory. Furthermore, the last two models only consider the undulation of yarns in the given loading direction. These three methods are based upon assumptions for simplifying the yarns geometry (Figure 1-5) and/or mechanical behavior of the unit cell, which leads to loosing accuracy. However, they can still effectively predict the equivalent material properties of impregnated woven fabrics for specific loading modes. Results from experimental studies conducted on woven fabrics have revealed that for 8H satin composites, the methods proposed by Ishikawa and Chou (1982) can, in many cases, predict the material properties of woven fabrics. Generally, however, there are still discrepancies between the results of the experiment and the values determined by the prediction (Ishikawa et al. 1985). (a) (b) (c) Figure 1-5: Geometrical assumption of woven fabrics used in Ishikawa methods (a) mosaic model (b) fibre undulation model (c) bridge model (Ishikawa and Chou 1982) 9 Naik and Shembekar (1992) proposed a two-dimensional model based on Ishikawa and Chou (1982) models that considered the effect of undulation on both warp and weft yarns for more detailed analysis of woven fabrics. This model, however, still relies on geometrical simplifications and may not be used for finding the deformation details in all possible modes. In particular, it is not effective in including the geometrical changes that are observed during the loading of woven fabrics (especially dry fabrics), which can be a significant source of nonlinear behavior. Taking advantage of numerical methods, Whitcomb (1989) conducted studies on the numerical analysis of woven fabrics using a 3D finite element method. Here, he emphasized the aspects of the FEA that are different from the classical methods used for unidirectional laminated composites. As a continuation to this idea, Woo and Whitcomb (1993) and Whitcomb et al. (1994) developed a new element type for the finite element method in order to account for fabric microstructure. They called this new element type “macro-elements.” Cox et al. (1994) and Xu et al. (1995) introduced binary models as a computationally effective tool for analyzing different types of textile composites. This model analyzes the textiles on the tow gauge. Binary model formulation represents the tow architecture by a set of strings of one-dimensional line elements, which identifies the location of the axis of one of the tows in the three-dimensional space. The one-dimensional tow elements constitute one finite element mesh. The solid geometry of the component is also meshed by the conventional three-dimensional solid finite element called “effective medium element.” The effective medium elements constitute a second finite element mesh. In real terms, this means that the one-dimensional elements represent the axial properties of tows (axial stiffness), while the effective medium elements represent the matrix dominated properties of the composite (transverse and shear stiffnesses, and Poisson’s effect). In the next step, the dual element system of the two meshes and the corresponding constitutive laws are combined into a single problem by coupling the meshes using multi-point constraints, which can be done by customized or commercial finite element software. In the end, the combined problem is solved and the stress and strain values are calculated using finite element software (McGlockton et al. 2003). Vandeurzen et al. (1996a-b) developed a general library of building blocks for textile composites, which includes both geometrical models of fabrics and the prediction of their engineering constants for laminate-based models on macro-cells and micro-cells. This method 10 was named the Combi-Cell model by the authors. While their model is very accurate in predicting the fabric tensile behavior, it is, however, not able to predict the shear modulus correctly. They also implemented their model into a comprehensible software package called TEXCOMP, which is based on Microsoft Excel®. It has been suggested that more accurate results can be achieved using the finite element method, although this was impractical because of computational limitations at that time. More recently, the authors developed an analytical model of stress distribution in a general stress model for woven fabrics (Vandeurzen et al. 1997). In this model, a multilevel-multistep approach was applied for the first time in woven fabrics. Indeed, the model proposed by Chen and Cheng (1993) for unidirectional composites was extended for woven fabrics in each step. Bigaud and Hamelin (1997) developed a computer program called TIS3D for the geometrical modeling of textile composites along with a numerical procedure for the estimation of elastic and failure properties of textiles in a polymer matrix. The mechanical analysis of stiffness behavior is done in multi-scale: micro, meso and macro. The equivalent properties at each scale are determined from the lower scale and are based on an energy method, Lagrange multipliers (Chen and Cheng 1993), and averaging of strain and stress. A similar model was developed by Tan et al. (1997) to evaluate the three-dimensional elastic properties of woven fabrics using micro/macro elements. This model can be used for the evaluation of different weaving types. Comparisons between the results from theoretical data and the finite element method show a good contingency in the linear deformation range. Boisse et al. (1997) proposed a finite element simulation method for the 3D shaping process of fabrics. The constitutive model for the construction of finite elements was based upon the bar model (Kawabata 1973a-c). First, a series of experimental bi-axial tests were designed to investigate the fabric behavior and warp/weft interactions. Next, the finite element model was used for the simulation of stamp forming using a square box punch, as well as an ellipsoidal punch and die. Tabiei and Jiang (1999) developed a micromechanical composite material model for woven fabrics including the nonlinear stress-strain behavior. They coded a user defined subroutine in the Abaqus finite element package to define the constitutive behavior of woven fabrics including the nonlinearity. Using iso-stress and iso-strain assumptions and meso-mechanical analysis, the constitutive equations were averaged in each representative unit cell. Then effective stress and 11 strains in each unit cell were passed to the finite element code for a global analysis of the material structure. Comparing results from case studies using this model and linear analysis showed that, in axial loading, there is a slight difference between the results from nonlinear analysis and linear analysis. However, under shear deformations, the effect of nonlinear material behavior was noticeable. Indeed, experimental results from picture frame and trellising tests highlighted that shear nonlinearity exists in the majority of polymer composites. Tabiei and Yi (2002) compared several methods previously presented by Tabiei and Jiang (1999), Tabiei et al. (2000), Jiang et al. (2000) , Tanov and Tabiei (2001) and arrived at a new simplified method for woven fabrics. They discovered that the simplified formulation could make their method more cost effective for numerical calculation of multi-scale problems. Kuhn et al. (1999) developed a semi-analytical and computationally effective method for calculating the elastic micro-fields in a woven fabric under in-plane loadings including uni-axial, bi-axial and shear. They used the Rayleigh-Reitz method to estimate the displacement field in the laminate. While most the models previously discussed in the above section were designed for systems with matrix surrounding the yarns, in many other applications the designers have to, instead, deal with dry fabrics. Yarns in dry fabrics are made of bundles of fibres without any consolidation to hold them together. This gives dry yarns negligible shear stiffness values in their constitutive models, as well as highly nonlinear and strain dependant transverse stiffness. As a result, the study of mechanical properties in these types of materials requires modified techniques for accuracy and constitutive models that are close to the real system. Gasser et al. (2000) studied the mechanical behavior of dry fabrics by focusing on the local (meso-scale) behavior for better understanding of their influence on the macro-level behavior. They developed a procedure for nonlinear finite element analysis specifically for studying mechanical behavior of dry fabrics based on their constituents. Their results were compared to biaxial tests and showed an agreeable accuracy of this method for predicting the material behavior in the macro-level. Additionally, they were able to establish a general foundation to simulate the drawing of woven fabrics in the manufacturing processes (Boisse et al. 2001; Buet-Gautier and Boisse 2001; Hivet et al. 2002; Boisse et al. 2005). Initially, experimental bi-axial tests were utilized to study the effect of nonlinearity in meso-scale on the macroscopic behavior. Afterwards, the results from bi-axial tests were incorporated into models based on bars and yarns (Kawabata 1989; Boisse et al. 12 2005), which were then used in a special finite element formulation for fabrics. Boisse et al. (2001) and Hivet et al. (2002) also included a meso-level 3D finite element analysis (Figure 1-6) for better understanding of the local phenomena in the fabric meso-level, especially yarn crushing, which is one of the main sources of nonlinearity. It was shown that nonlinear behavior is mainly due to geometrical changes rather than contact nonlinearities. The authors also used virtual test data from 3D modeling of unit cells under biaxial and shear tests in meso-level in order to develop a new three-node finite element model, which can be used effectively in the simulation of draping in woven fabrics (Hamila and Boisse 2007). It should be further emphasized that because of the special mechanical behavior of yarns at meso-level, classical methods of finite element are not suitable. Most routine methods and finite element packages use rate constitutive equations (or hypo-elastic constitutive equations) for finite strain analyses. Hamila and Boisse (2007) used the standard Green Naghdi or Jaumann objective derivatives, which are not suitable for the analysis of finite strain in fibrous materials where there exists an immense drop in stiffness in the transverse direction, as compared to the longitudinal direction (Badel et al. 2007). The constitutive equations in dry fabric materials should be handled using fibre rotation frame, which requires modifications to the classical equations used for finite rotations (Badel et al. 2008a; Badel et al. 2009). Figure 1-6: 3D finite element models used by Boisse et al. (2001) Peng and Cao (2000 and 2002) applied a dual homogenization technique to predict the mechanical behavior of woven fabrics from a sequence of models for fibres. In the first step, a unit cell was made for the yarn section to estimate the effective elastic constants of fibre yarn. Then, another unit cell in meso-level was constructed to represent the periodic structure of woven fabric. Three basic numerical tests (uni-axial, bi-axial and shear) were performed using 13 this unit cell. Finally, by correlating the force versus displacement of this unit cell and a four node shell element with the same size, the effective mechanical properties of the fabric were evaluated. A comparison between the results from the numerical shear test and experimental results (Chen et al. 2001) demonstrated the effectiveness of this procedure. However, their analysis was based on the classical finite element procedures and particular mechanical behavior of dry yarns, which was considered by Gasser et al. (2000), was neglected. While this makes the method easier to apply, it also reduces the accuracy of predicted results to some degree. Carvelli and Poggi (2001) used a similar method for estimating the stiffness and strength in woven fabrics, the implementation of the failure criteria into dual homogenization techniques and 3D finite element analysis of unit cells in meso-level. Lomov et al. (2001) suggested the use of an integral modeling and design tool, which can be effectively used for the meso-level modeling of textile structures. It is also efficient in stiffness analysis for finding homogenized material properties, studying local phenomena, and permeability analysis. The suggested method could be used as a foundation to start the analysis of textiles. This method was also shown to be effective in constructing models for complex shapes and textile architecture. The authors incorporated some of their techniques into an easy to use software package called WiseTex, which can be used as a powerful tool in the manipulation and visualization of fabrics. Roy and Sihn (2001) derived a mixed three-dimensional variation model for stress analysis of the Representative Volume Element (RVE) of woven fabric composites based on the Reissner variational principle (Reissner 1950; Pagano 1978; Harrison and Johnson 1996). Their formulation yields a set of second order differential equations that can be solved numerically by the finite element or finite difference methods (Sihn and Roy 2001). However, the numerical results using polynomial shape functions revealed singular stress behavior in the interface region of the material junction, which could not be handled even by using higher order polynomials or finer mesh (Sihn and Roy 2001). Sihn et al (2004a) modified their proposed method by using B- spline displacement approximation, which could provide continuity of stress and displacement within each yarn and matrix sub-region. The new method could reduce singularity effects in the presence of perfect bonding between the yarns and the matrix with small crimp. However, for increased crimp angle it was still sensitive to singularities with a small variation in the stiffness of matrix materials (Sihn et al. 2004b). 14 Page and Wang (2002) developed a finite element model for the prediction of shear reaction force in fabrics. They used a simplified geometry for yarn cross section and generated the yarn waviness by applying the appropriate boundary condition and loading, which simulates the weaving process in yarns. Transverse stiffness modulus was estimated by a series of fabric compression tests using ASTM 1388. Results showed a highly nonlinear response for the transverse stiffness, which led to an iterative solution method for determining its approximate value during the analysis. This procedure was very similar to the method used by Gasser et al. (2000), which was based on an inverse method for finding the stress-strain relationship for transverse stiffness. Eventually, using these results and a previously derived analytical model, shear force was calculated as a function of shear angle, which agreed well with the experimental observations. Ballhause et al. (2005 and 2006) used a Discrete Element Model (DEM) for dynamic analysis of woven fabrics and fabric reinforced composite materials. Here, the unit cell was described by an assembly of discrete mass points and rheological elements. The system of elements were assembled in order to represent relevant deformation mechanisms like crimp interchange, trellising or locking. This model was similar to the model proposed by Kawabata et al. (1973a-c), however, it also includes lumped mass elements, which can make it suitable for dynamic modeling (Figure 1-7). Figure 1-7: Discrete Element Method; A: warp, B: weft, C: transverse yarn compression, D: Rotational Resistance, mi: concentrated mass. Durville (2005) proposed a new approach for including all fibres of a fabric and their interactions in the finite element simulation of fabrics under different loading conditions. Geometrical representation of their model is shown in Figure 1-8. The originality of this approach was the model used for contact-friction interactions between fibres, and the solution of the nonlinear problem that followed by robust algorithms. The advantage of this method is that it 15 provides a better simulation of the phenomena taking place at the level of fibres (micro-level). This can be fascinating for a small sample of fibres with a low amount of fibre in each yarn. However, for cases with enormous number of fibres in each yarn and huge samples, it can be computationally very expensive. Figure 1-8: Geometry and mesh of the model for a plain weave fabric including all fibres in each yarn and their interactions (Durville 2005) King et al. (2005) proposed a new approach for continuum modeling of fabrics. In this approach, the main idea was to select a geometric model for the fabric weave coupled with constitutive models for the yarns behavior. An energy minimization method was used to relate the fabrics structural configuration to the macroscopic deformation. The geometrical representation of the fabric in their models is very similar to Kawabata (1973a-c), although, the continuum models developed in their approach could not directly model individual yarns. Instead, they contain information about the behavior of the fabric meso-structure. The main advantage of their continuum models was to account for meso-level phenomena, such as crimp interchange, locking and resistance to relative yarn rotation caused by yarn interactions. Xue et al. (2005) developed an integrated meso-macro model for the prediction of mechanical properties of woven fabrics under large deformation. They had previously suggested a non-orthogonal constitutive model for the prediction of mechanical properties in macro scale (Xue et al. 2003) where they used experimental measurements for finding the constants of their model. Their integrated model begins with a geometrical description of the yarn and the unit cell including geometrical changes under trellising shear deformation (Mcbride and Chen 1997). 16 Consequently, mechanical analysis on a unit cell is conducted to determine the equivalent shear properties. Boubaker et al. (2007a 2007b) developed a meso-level discrete model in which an element of fabric was modeled by a set of grid nodes endowed with a mass and connected with flexional and stretching springs. Although the concept of using a discrete element was a method that had been previously developed, taking into account a better model of yarn-yarn interactions and friction was an original approach in their model. Here, mechanical behavior was analyzed using an energy analysis, taking into account the compression strain energy of the yarns and the work of the reaction forces exerted between yarns. Next, a discrete variational principle was used to account for the presence of the non-conservative forces arising from friction. Hivet and Boisse (2008) developed a simplified meso-level model using an axial element for the mechanical behavior of woven fabrics. The new model was based on the 3D geometrical description of the unit cell and ensures 3D consistency, meaning there are no voids or interpenetration between the yarns (Figure 1-9) due to geometrical assumptions. This new model can be considered as a modification of Kawabata et al. (1973a-c) model for more consistency, which was one of the problems in the original model. Figure 1-9: Interpenetration between the yarns (Hivet and Boisse 2008) Recently, Lin et al. (2009) used a finite element model for the prediction of shear behavior in woven fabrics. Here, they attempted to incorporate the details of picture frame kinematics in their model. Based on their results, the kinematic boundary conditions in a finite element at a 17 small shear angle and a large shear angle regime should be different. This is a result of the different mechanisms that are dominant in each regime. 1.4 Motivation and objectives of the work There is a significant amount of work done on the meso-level modeling of woven fabrics. On the other hand, most of these works consider the fabric at its perfect condition with no material flaws or imperfections. However, in a real case, what occurs to a fabric during its production or handling process can induce some flaws into the material. Furthermore, it is known that in composites there are always effects caused by intrinsic uncertainty in the constituent material properties and fibre geometries, which can lead to non-repeatable behavior of final products (Skordos and Sutcliffe 2008). Conversely, there are a number of studies on the identification and evaluation of these flaws using NDT (Non-Destructive Testing) and image processing techniques (Creighton et al. 2001; Bu et al. 2009). Even in laboratory samples that are used for experimental measurements of mechanical properties there are reports of discrepancy in the fabric mechanical behavior, which is a result of uncertainty and imperfections (Milani and Nemes 2004; Lomov and Verpoest 2006). Elleithy (2000) described some of the most common imperfections in the woven fabric reinforced composites. He concluded that imperfections are common in all scales of study from lamina plies and yarns to single fibres. Actually, he reported that at any level, all imperfections contribute to the mechanical response of the fabric structure. Thus, developments of precise mathematical models, which can take into account these imperfections and uncertainties, are necessary. Milani and Nemes (2004) used an inverse method for characterizing material behavior in textile composites, including the effect of noises in their response. Following from this, Milani et al. (2007) developed an inverse methodology for the approximation of the non-uniform state of fibre misalignment. Their results can be used for the enhancement of constitutive models and numerical procedures. Zeman and Sejnoha (2004), followed by Sejnoha and Zeman (2008), noted the importance of including the geometrical imperfections in the finite element modeling of a unit cell. They derived a mathematical framework for determination of a reliable representative volume element (unit cell) under imperfections and used morphology of imperfect material systems by an appropriate set of statistical descriptors. Descriptions were then 18 introduced into an optimization environment that provided a tool for a unit cell to statistically resemble the actual system. Eventually, the finite element model of a woven composite was formulated based on geometrical parameters and used to predict the overall response of the composite by a numerical homogenization method. However, in their method they considered all imperfections as part of the given model and did not concentrate on the effect of each factor and at a higher level or the possible interactions between them. The aim of this study is to develop a comprehensive analysis on the effect of imperfections and uncertainty factors on the mechanical response of woven fabrics at meso-scale. The imperfect models are analyzed using 3D finite element modeling in the Abaqus software and defining user defined subroutines for simulating the specific behaviors that exist in woven fabrics. The imperfections include material and geometrical uncertainties. Komeili and Milani (2009) presented a basic study via a one-factor-at-a-time sensitivity analysis on the unit cell of a balanced plain weave fabric that was also used by Peng and Cao (2002). As a continuation of their study, this thesis will be a comprehensive full factorial design analysis on the effect of a larger set of uncertainty factors and their interactions. In addition, a different constitutive model for yarns, which can represent physical properties of dry fabrics more accurately (Badel et al. 2008a), has been implemented in both an explicit and a new implicit subroutine. 1.5 Thesis outline This thesis is organized into six chapters. The first chapter, presented above, was a brief introduction on the general applications of composite materials and, more specifically, woven fabric textile composites. This chapter also included a general literature review on the modeling aspects of woven fabrics. The second chapter provides an overview about finite element modeling of woven fabric unit cells in meso-level. This discussion is further extended in Chapter 3, which explains details of the finite element model employed in this thesis. Here the steps to produce the unit cell model in Abaqus is outlined, including the creation of the CAD model and the writing of the user defined FORTRAN subroutines for specific material constitutive models of yarns. Chapter 4 examines the uncertainty factors that exist in plain weave fabrics and their range of variation. These uncertainties are separated into two categories, which will be studied in the subsequent sensitivity analysis. Chapter 5 presents the results and the discussions on the 19 sensitivity analysis via full factorial designs. The important factors and interactions are identified using half-normal probability plots, factorial effects, and percentage contributions. Finally, Chapter 6 concludes the thesis with a summary of the procedures used and the results obtained in this study, as well as recommendations for future work. 20 Chapter 2: FEM modeling of fabric unit cell at meso-level In this chapter, the general procedure for finite element modeling of a fabric unit cell is presented. Entire steps from geometrical modeling to choosing the appropriate solver are covered. The main goal is to provide a perspective on building meso-level unit cell models for 3D finite element analysis of dry woven fabrics. 2.1 Introduction The fabric considered in this section is a balanced plain woven fabric (Figure 2-1). This fabric is comprised of the most basic pattern, because of its simplicity in term of the manufacturing process and, as a result, it has a very wide range of applications in a variety of industries. Thus, the study of its properties using 3D finite element (FE) modeling is of great interest and there are several studies that have been conducted on this fabric (Gasser et al. 2000; Peng and Cao, 2002; Boisse et al. 2007). Using 3D finite element modeling in meso-level can not only be helpful in finding the material behaviour of unit cell, but also provides insight into the deformation phenomenon and interactions that occur in this scale. 3D finite element modeling is a very suitable methodology to provide results for factorial designs used during the sensitivity analysis, especially in this study that will address uncertainties and imperfections at meso-level. Figure 2-1: Balanced plain weave fabric (Badel et al. 2008a) 21 2.2 Standard procedures for 3D FE modeling of unit cells Procedures used for 3D FE modeling of unit cells can be broken into seven steps: 1- Geometrical modeling and meshing: First, the geometry of unit cell should be described clearly and CAD models should be prepared to be meshed in a finite element package. Deciding on the basic representative volume element (unit cell), which must be neither too complicated to increase the computational cost of the model, nor too simple to neglect details of the analysis or results in the conclusions is crucial. The geometry of each yarn and assumptions on the parameters for defining the yarns arrangement is also important. Attention should be paid in defining the geometry of each yarn, as well as the yarns position in a unit cell in order to ensure that geometrical consistency in the CAD model is satisfied (Hivet and Boisse 2005). 2- Material properties: As previously discussed, due to the multi-scale nature of the study, the assumption of solid homogenous material for yarns is merely a simplification for modeling to be computationally feasible. There are some approaches that use special material constitutive models for yarns to take into account effects from fibres in micro- level (Gasser et al. 2000; Sherburn 2007). On the other hand, there are others that utilize classical elastic material properties and finite element procedures (Peng and Cao 2000). Material properties can be determined using experimental measurements on yarns, an inverse characterization method on experimental results (Gasser et al. 2000), or results of homogenization from material properties in micro-level (Takano et al. 1999; Peng and Cao 2002). 3- Periodicity boundary condition: The unit cell is modeled as a system, which is cut off from a larger system (whole fabric) and is assumed to represent the entire fabric. Hence, its mechanical response cannot be considered merely by using an isolated structure. There are a set of equations and constraints on the boundaries of a unit cell, which make it act like a representative volume element in a periodic structure. These sets of boundary conditions should be applied on any unit cell regardless of loading conditions. 4- Loading boundary conditions: Depending upon the type of loading (axial, shear, combined etc) that is applied on the macro-level model, boundary conditions should be applied on the unit cell to simulate the loading of the fabric. For homogenization 22 purposes, usually the stress-strain curves are of importance so that a load control analysis can be conducted on the model to find the reaction force at boundaries. 5- Nonlinearity sources: Nonlinear material properties and contact interaction between the yarns can be identified as the sources for nonlinear behaviour. Moreover, because of the undulation of yarns under axial loading and large displacement under shear loading, the finite strain theory should be used in finite element analysis to obtain accurate results. 6- Finite element solver: A sensible choice of a finite element solver is one of the most challenging steps to obtain reliable results and to reduce computational costs. In addition to the option of writing a customized finite element code, there are a variety of commercial packages, such as ANSYS, Abaqus etc. The available options for each package differ based upon the packages capabilities, selected method and numerical solvers. The Abaqus 6.8-1 package is used for finite element modeling of unit cell in this project. It includes two solvers (integrators), which are appropriate for the current study: (1) “Dynamic/Explicit” solver, which uses an explicit non-iterative time dependant solver and (2) “Static/General,” which uses an implicit static iterative solver. Both of these integrators can use features 1-5 previously mentioned, which are essential for fabric modeling. 7- Post processing: Finite element packages can report many different variables as output of an analysis. These variables can be selected from built-in variable lists or from customized variables. The reaction force versus displacement (stress vs. strain) for displacement control analysis can be the first choice of post processing output to assess the response of unit cells. 2.3 Geometry of unit cell The preliminary step to start the analysis is defining the appropriate unit cell (representative volume element). The unit cell should be as small as possible to reduce modeling and computation costs, however, at the same time, it should not be too small so that it will not reflect the actual behaviour of a fabric. There is a variety of suggested unit cells for dry woven fabrics in the literature. The difference between these unit cells is based upon the assumptions that are made on the periodicity, type of loading and geometry. 23 In order to analyze fabrics under axial loading Gasser et al. (2000) used an elementary pattern depicted in Figure 2-2(a) for general plain woven fabrics, which can be reduced to a quarter of elementary pattern for a balanced fabric (Figure 2-2(b)). Cross sectional geometry of the yarns are assumed to be lenticular, symmetric and truncated at the extremities (Hivet and Boisse 2005). This geometry was determined using different optical observations. (a) (b) Figure 2-2: The unit cell geometry used by Gasser (2000); (a) an elementary pattern (b) the quarter of elementary pattern. Continuing the same modeling approach Badel et al. (2007) used two more generalized representative volume elements for modeling of shear loading. These representative volumes needed two different sets of boundary conditions to satisfy periodicity. However, the final results from the finite element analysis proved that both models lead to similar results and experimental measurements. Figure 2-3 shows these two representative volume elements. Figure 2-3: Two unit cell types used by Badel et al. (2007) for modeling of shear loading Peng and Cao (2000) used a different definition for unit cell in the analysis of plain woven fabrics under general loading conditions (Figure 2-4a). To check the accuracy of this unit cell they compared results obtained from experimental data and an expanded unit cell composed of 24 four basic unit cells as shown in Figure 2-4b (Chen et al. 2001). In this model, four sinusoidal curves define the geometry of the cross section for each yarn and its axis; this geometry is proposed based on the measurement by Mcbride and Chen (1997). Due to the flexibility of this unit cell definition for different types of loading modes, it has been selected for the current study during the finite element analysis. Figure 2-4: Unit cells used by Chen et al. (2001) a) Basic unit cell b) Extended unit cell. According to McBride and Chen (1997), a unit cell can be defined using three constants, which include yarn width w, yarn spacing s and yarn height h (Figure 2-4a). The sinusoidal equations for defining the geometry of unreformed yarn are: ݕଵሺݔሻ ൌ ଶ ቂܿݏ గ ௫ ௦ 1ቃ ሺ0 ൏ ݔ ൏ ݏሻ (2-1) ݕଶሺݔሻ ൌ ଶ ቂܿݏ గ ௫ ௦ െ 1ቃ ሺ0 ൏ ݔ ൏ ݏሻ (2-2) ݕଷሺݔሻ ൌ െ݄ ܿݏ గ ௫ ఉ ቀ0 ൏ ݔ ൏ ௪ ଶ ቁ (2-3) ݕସሺݔሻ ൌ െ݄ ܿݏ ቂ గ ൫௫ିሺ௦ିఉሻ൯ ఉ ቃ ቀݏ െ ௪ ଶ ൏ ݔ ൏ ݏቁ (2-4) ߚ ൌ గ ௪ ଶ ௦ቂ௦మቀഏ ೢ ర ೞ ቁቃ By applying these formulas, yarn length (L) and the cross sectional area of yarn in the unit cell (A) can be derived: 25 ܮ ൌ ට1 ቀ ௗ ௬భ ௗ௫ ቁ ଶ௦ ݀ݔ ൎ 1 ଵ ଶ ቀௗ ௬భ ௗ௫ ቁ ଶ ൨௦ ݀ݔ ൌ ݏ మ గమ ଵ ௦ (2-5) ܣ ൌ ሺݕଷ െ ݕଶሻ݀ݔ ೢ మ ൌ ቀ మା௪మ ସ ቁ ଶ ܽݎܿݏ݅݊ ቀ ଶ ௪ మା௪మ ቁ െ ௪ ଶ ௪ మିమ ସ (2-6) 2.4 Material modeling The multi-scale nature of unit cell modeling implies that all of the mechanical properties that are used for meso-level model should consider the mechanical behaviour of fibres in micro-level. There are different approaches for finding the constituent mechanical properties of yarns. Some of these methods are based on the measurements that are done on single yarns or fabrics. They can use inverse methods to find constants and parameters that are defined in their meso-level model. On the other hand, a homogenization (Bakhvalov and Panasenko 1989) can be used to find the properties of a yarn based on the mechanical properties of fibres in micro-level. Nonetheless, when using any of these approaches the material behaviour can be defined as linear or nonlinear, elastic or plastic or temperature/rate dependant. 2.4.1 Micro-scale homogenization In order to find the effective material properties of a yarn from fibres the basic unit cell of a fibre must be defined. This means that the yarn itself is a heterogeneous periodic medium, which is made of fibre unit cells. Peng and Cao (2002) defined a micro-structure unit cell as shown in Figure 2-5. Using this unit cell, they presented the homogenization formulation that can be used in a numerical package (e. g. Abaqus) to find the averaged properties of yarn material. 26 Figure 2-5: Unit cell for a micro-structure (Peng and Cao, 2002) Global and local coordinates in this micro-structure model is defined by x and y, respectively. The relationship between the two coordinate systems can be defined using ߳ as a very small positive number that is the ratio between the dimension of a unit cell and whole structure. ൌ /߳ (2-7) The asymptotic expansion of displacement field is the essential part in homogenization. Hence, the displacement field can be written as: ݑ ൌ ݑ ሺ, ሻ ߳ ݑ ଵሺ, ሻ ߳ଶ ݑ ଶሺ, ሻ … (2-8) And keeping in mind that for a general function, ߶ : డ థሺ,ሻ డ௫ ൌ డ థ ሺ,ሻ డ௫ ଵ ఢ డ థሺ,ሻ డ௬ (2-9) For strain tensor ߝ: ߝ ൌ ଵ ଶ ൬డ ௨ డ௫ೕ డ ௨ೕ డ௫ ൰ ൌ ଵ ఢ ߝ ିଵሺ, ሻ ߝ ሺ, ሻ ߳ ߝ ଵ ሺ, ሻ … (2-10) where ߝିଵሺ, ሻ ൌ ߝ௬ሺሻ (2-11) ߝ ሺ, ሻ ൌ ߝ௫ሺሻ ߝ௬ሺଵሻ (2-12) ߝ ଵ ሺ, ሻ ൌ ߝ௫ሺଵሻ ߝ௬ሺଶሻ (2-13) ߝ௫ሺݒሻ ൌ ଵ ଶ ൬డ ௩ డ௫ೕ డ ௩ೕ డ௫ ൰, ߝ௬ሺݒሻ ൌ ଵ ଶ ൬డ ௩ డ௬ೕ డ ௩ೕ డ௬ ൰ (2-14) 27 Through the assumption of periodic media over the entire structure, elastic coefficients can also be defined as a periodic function of ߳ in x (global) coordinate system. ܥ ఢ ൌ ܥሺ ߳⁄ ሻ ؠ ܥሺሻ (2-15) Stress function can be expressed as: ߪ ఢ ൌ ଵ ఢ ߪ ିଵሺ, ሻ ߪ ሺ, ሻ ߳ ߪ ଵ ሺ, ሻ … (2-16) ߪ ሺ, ሻ ൌ ܥ ఢ ߳ ሺ, ሻ (2-17) Therefore, the basic differential equations on this periodic media are presented as: ߪ, ఢ ݂ ൌ 0 ݅݊ Ω (2-18) ߪ ఢ ݊ ൌ ܶ ݊ Γଵ (2-19) ݑ ఢ ൌ 0 ݊ Γଶ (2-20) where Ω is the domain of the problem and Γଵ are Γଶ the force and displacement boundary conditions respectively. By substituting (2-17) in (2-18) and equalling the same powers of ߳ we obtain: డ ఙೕ షభ డ௬ೕ ൌ 0 (2-21) డ ఙೕ బ డ௬ೕ ൌ 0 (2-22) డ ఙೕ బ డ௫ೕ డ ఙೕ భ డ௬ೕ ݂ ൌ 0 (2-23) Equation (2-21) can be written in the variational form of: డ ఙೕ షభ డ௬ೕΩച ߜݑ ݀Ω ൌ ቀܥఢ డ ௨ೖ బ డ௬ ቁ ,Ω ച ߜݑ ݀Ω ൌ 0 (2-24) In the next step, a mean operator is defined that works on a Y-periodic function Φሺݕሻ as: lim ఢืశ Φ ቀ ఢ ቁ ݀ΩΩച ൌ ଵ || Φሺሻ ܻ݀ ݀Ω Ω (2-25) where |ܻ| represents the volume of the unit cell (micro-level unit cell). In this problem, the homogenization method consists of finding the limit solutions of (2-18)-(2- 20) as ߳ ՜ 0. Using (2-25) for equation (2-24) there is: lim ఢืశ ቀܥ ఢ డ ௨ೖ బ డ௬ ቁ ,Ω ച ߜݑ ݀Ω ൌ ଵ || ቀܥ డ ௨ೖ బ డ௬ ቁ , ߜݑ ܻ݀ ݀Ω Ω ൌ 0 (2-26) By applying the divergence theorem on (2-26): ଵ || ቀܥ డ ௨ೖ బ డ௬ ቁ , ߜݑ ܻ݀ ݀Ω Ω ൌ ଵ || ׯ ܥ డ ௨ೖ బ డ௬௦Ω ݊ ߜݑ ݀ݏ ݀Ω (2-27) 28 that yields డ ௨ೖ బ డ௬ೕ ൌ 0 (2-28) In other words ݑ is a function of x only. Therefore, (2-8) can be simplified as follows: ݑ ൌ ݑ ሺሻ ߳ ݑ ଵሺ, ሻ ߳ଶ ݑ ଶሺ, ሻ … (2-29) This formula implies that physically, ݑሺሻ is the macroscopic displacement while ݑ ଵሺ, ሻ, ݑ ଶሺ, ሻ, … are perturbing displacements relevant to micro-structure. In the next step, equations (2-12) and (2-17) are written as: ߪ ൌ ܥ ఢ ቀడ ௨ೖ బ డ௫ డ ௨ೖ భ డ௬ ቁ (2-30) Using (2-30) in the variational form for (2-22): ቂܥ ఢ ቀడ ௨ೖ బ డ௫ డ ௨ೖ భ డ௬ ቁቃ ,Ω ച ߜݑ ଵ ݀Ω ൌ 0 (2-31) lim ఢืశ ቂܥ ఢ ቀడ ௨ೖ బ డ௫ డ ௨ೖ భ డ௬ ቁቃ ,Ω ച ߜݑ ଵ ݀Ω ൌ ଵ || ቂܥ ቀ డ ௨ೖ బ డ௫ డ ௨ೖ భ డ௬ ቁቃ , ߜݑ ଵ ܻ݀ ݀Ω Ω ൌ 0 (2-32) By applying the integration by part and noting that we have displacement boundary conditions here (it implies ߜݑଵ ൌ 0 at unit cell boundaries) and keeping in mind that ݑ is only a function of x, we obtain: ܥ డ ௨ೖ భ డ௬ డ ఋ௨ భ డ௬ೕ ܻ݀ ݀Ω Ω డ ௨ೖ బ డ௫Ω ൬ ܥ డ ఋ௨ భ డ௬ೕ ܻ݀ ൰ ݀Ω ൌ 0 (2-33) To solve this formula Ψ is introduced which satisfies the following equation: ܥ డ Ψೖ డ௬ డ ఋ௨ భ డ௬ೕ ܻ݀ଢ଼ ൌ ܥ డ ఋ௨ భ డ௬ೕ ܻ݀ (2-34) The substitution of (2-34) into (2-33) yields: ܥ డ ௨ೖ భ డ௬ డ ఋ௨ భ డ௬ೕ ܻ݀ ݀Ω Ω డ ௨ೖ బ డ௫Ω ൬ ܥ డ Ψೖ డ௬ డ ఋ௨ భ డ௬ೕ ܻ݀ ൰ ݀Ω ൌ 0 (2-35) that can be written in the following form by applying the divergence theorem: ׯ ܥ ݑଵ ݑ డ ఋ௨ భ డ௬ೕ௦Ω ݀ݏ ݀Ω ׯ ܥ Ψ ݊ డ ௨ೖ బ డ௫ డ ఋ௨ భ డ௬ೕ ݀ݏ ݀Ω௦Ω ൌ 0 (2-36) Hence: 29 ݑ ଵ ൌ െΨሺ, ሻ డ ௨ೖ బ డ௫ (2-37) and substituting it into (2-30) is obtained: ߪ ൌ ൬ܥ െ ܥ డ Ψೖ డ௬ ൰ డ ௨ೖ బ డ௫ (2-38) By integration over the unit cell volume the stress-strain relations for the homogenized elastic body is presented as: ߪ ൌ ܥ ு డ ௨ೖ బ డ௫ (2-39) where ߪ is the average of ߪ over the unit cell volume and ܥு is the homogenized elastic constants. They can be written in the form of: ߪ ൌ ଵ || ߪ ሺ, ሻ ܻ݀ (2-40) ܥ ு ൌ ଵ || ൬ܥሺሻ െ ܥሺሻ డ Ψೖ డ௬ ൰ ܻ݀ (2-41) The key to finding the mechanical properties of the micro-scale unit cell is solving (2-41), which can be easily done using a finite element tool. To make it suitable for finite element solvers, (2-41) can be written as: ቀ ் ܻ݀ ቁ Ψ ൌ ் ܥ ܻ݀ (2-42) in which is the stress-strain matrix of the constitutive materials, is the displacement-strain matrix and ܥ is a vector of column ij of the . The finite element solver should be used to find the characteristic displacement matrix which is Ψ. Eventually, homogenized elastic constants of a micro-level unit cell can be found by following averaging integration: ሾுሿ ൌ ଵ || ሺܫ െ શሻ ܻ݀ (2-43) having શ ൌ ሼΨଵଵ,Ψଶଶ,Ψଷଷ,Ψଷଶ,Ψଷଵ,Ψଵଶሽ (2-44) This formulation was used by Peng and Cao (2002) in Abaqus/Standard solver to find the material properties of yarns in an E-glass/PP yarn with 70% E-glass volume fraction. Results showed a reasonable match for some constants compared to conventional analytical methods (Kaw 1997; Halpin and Tsai 1969) and experimental results in the available cases. 30 2.4.2 Constitutive model based on physical observations Another alternative approach for finding the equivalent material properties of yarns in meso- level is using some physical interpretations from their characteristics along with experimental measurements. Moreover, sometimes applying inverse identification methods to find material properties indirectly can be very helpful (Gasser et al. 2000). For example, from the fibrous nature of fabrics it is apparent that there are some zero or quasi zero stiffness values in the constitutive model of yarn materials. For instance, the potential of yarns in bearing longitudinal compression is negligible. In addition, because there is no resin consolidated in the space between the dry fibres, they can easily slide on each other, which means a very small magnitude of shear modulus. For the same reason, the crushing law is developed to model the transverse behaviour of fibbers. According to the crushing law, the transverse modulus of yarns is a function of contact conditions between yarns. This law is developed based on the observation that the more compression applied on the yarn, the stiffer it becomes. Kawabata et al. (1973a-c) attempted to write the transverse modulus as a function of contact force. However, if the yarn is under longitudinal tension, it is more difficult to crush. As a result, Gasser et al. (2000) developed transverse modulus as a function of both longitudinal and transverse strain using the following formula, assuming direction 1 is along the yarn axis: ܧ ൌ ܧ |ߝ | ߝଵଵ ݅ ൌ ሼ2,3ሽ (2-45) where ܧ, n and m are constants that can be determined using inverse identification method. There are also some other crushing functions available in the literature (Badel et al. 2008b-c). 2.5 Periodicity boundary conditions As previously mentioned, the entire fabric system is simplified into a Representative Unit Volume (or unit cell) for the sake of simplicity in numerical analysis. However, to keep the nature of the fabric yarns, periodicity boundary conditions should be implied on the model in order to consider the effects from adjacent unit cells. Applying simple symmetry boundary conditions may not work properly in this case, because, if symmetry conditions are applied on four sides of a unit cell, it is confined in a rectangular box, which prohibits it from deforming. Instead, the periodicity can be simply obtained by forcing nodes in opposite corners to have the same transverse displacement and the same rotation over global opposite edges (Gasser et al. 31 2000). To express this condition in a mathematical form (Badel et al. 2007), assume that the whole fabric is made of repetitions of one unit cell using translational vectors , which is made of two orthogonal vectors and located in the plane of fabric (Figure 2-6). ൌ ∑ ݉ఈ ఈଶఈୀଵ ݉ఈ א Գ (2-46) Figure 2-6: Initial and deformed meso-structure (Boisse et al. 2007) Assuming that boundary values for a property in a unit cell is defined by ߲ܸ, it can be split into four separate boundary conditions on four surfaces or two pairs ሼ߲ ఈܸି, ߲ ఈܸାሽఈୀଵ,ଶ as shown in Figure 2-7. Paired points are presumed to be two points of boundary that are images of each other by an elementary translation ఈ for ߙ ൌ 1, 2. However, because of the periodicity in the structure these two points are superimposed with the images of each other. Hence ఈି א ߲ ఈܸି; ఈା א ߲ ఈܸା ఈା െ ఈି ൌ ఈ (2-47) Deformed structure can be assumed to remain periodic during deformation (Figure 2-6). Defining ߮ to be the mapping function between the initial configuration () and the current configuration (), transformation of the structure can be decomposed into macroscopic (߮) and periodic (ݓ) terms: ߮ሺሻ ൌ ߮ሺሻ ݓሺሻ (2-48) In a homogenous loading case (such as uni-axial tension, bi-axial tension and pure shear), which is considered in this project, the macroscopic loading function is known and the aim of meso-level simulation is finding the periodic behaviour that is related to the local stress distribution in the unit cell. Periodicity of ݓሺሻ implies that: 32 ݓሺࢻିሻ ൌ ݓሺࢻାሻ (2-49) for paired points ࢻേ א ߲ ఈܸേ, ߙ ൌ 1, 2. In other words, every point in two opposing boundaries should have identical values for the periodic term of the transformation. Moreover, for every point that belongs to two boundaries (corner points) identical values are expected. To avoid rigid body motion in the unit cell ݓ ൌ 0 is applied at four corner points around the unit cell and the motion of these four corners is used to describe transformation in the system. Figure 2-7: Boundary condition on a unit cell (Boisse et al. 2007) Repetition of the unit cell shown in Figure 2-4a along with equation (2-4) does not represent the entire fabric. However, considering the symmetry of shape for balanced plain weave (Badel et al. 2007) it can reasonably predict the behaviour of the material. In other words, Figure 2-4a shows a unit cell that is simplified based on symmetry in geometry of the fabric. The periodic boundary condition for homogenous loading cases in this weave can be attained by allowing the side surfaces of yarns to remain plane and normal to the unit cell mid surfaces during the deformation (Peng and Cao 2002). By comparing the results using this unit cell and those from an extended unit cell that could satisfy equation (2-46), proves the efficiency of this unit cell (Chen et al. 2001), keeping in mind the noticeably lower amount of computational time. 2.6 Loading boundary condition Depending on the type of exerted load on the unit cell, different boundary conditions can be used for the simulation of fabric loading. In this study, three different types of loading modes are of interest: uni-axial, bi-axial and shear. The numerical experiments are designed to simulate the 33 behaviour of the unit cell in a woven fabric when tested in a universal tension machine (for tension modes) and a picture frame test (for shear mode). Therefore, the boundary conditions in the numerical methods are designed in a way to simulate these tests. For the uni-axial test, one edge of the unit cell is fixed and the opposite edge is loaded by a specified total displacement. To avoid the rigid body motion, one of the free sides is constrained in a direction normal to loading (Figure 2-8a). Bi-axial load is implied in a similar way, but the difference is that the displacement boundary condition is on all four edges of the unit cell (Figure 2-8b). For the shear test, the idea is based on what occurs in a real experiment using a shear diaphragm (Figure 2-9). In a shear diaphragm (frame), the rectangular specimen is fixed in a frame that is pinned in the corners so that it can be deformed to a diamond shape. This theoretically applies pure shear on the specimen. During the experiment, one corner of the frame is fixed while the opposite corner is moving on a slider. Reaction force in any of these corners is the output of the experiment. Similarly, in the finite element model, one corner should be kept fixed while the opposite corner is moving in the diagonal direction of the unit cell. The length of each edge should also be kept fixed to simulate the rigid frame of the test (Figure 2-8c). (a) (b) (c) Figure 2-8: Displacement boundary conditions for different types of loading (a) uni-axial tension (b) bi-axial tension (c) shear 34 Figure 2-9: A sample picture frame test setup (Cao et al. 2008) 2.7 Nonlinearity in analysis There are a number of reasons why using a linear finite element method for unit cell modeling is impractical. Some of the sources of nonlinearity are intrinsic in woven fabrics, while others can be considered or ignored during modeling based on the compromise that can be made between the analysis cost and accuracy of the model. The first source of nonlinearity is the contact between yarns. It is evident that the main purpose of weaving yarns together is to establish a contact between them to stand the applied loads via transverse pressure between the yarns and carry it through yarns tension. The contact condition between the yarn elements and possibility of a decrease or increase in the number of nodes that are engaged in the contact zone is a major source of nonlinearity. In addition, there is usually some tangential (frictional) behaviour that allows for the material system to be non- conservative and contributes to the nonlinearity. Contact is considered to be an intrinsic character of woven fabrics and cannot be neglected during the meso-level analysis. Undulation in the yarns in axial loading is another source of nonlinearity in the system (Boisse et al. 1997). In dry fabrics, because there is no consolidated material system to control the lateral movement of fibres, they can easily lose their shape and become straightened or wavier. This means that if the dry fabrics are under unbalanced tension, yarns that are located in 35 a direction parallel to the higher load tend to become straightened, while yarns in the other direction (normal) tend to keep their waviness or even increase it, especially under uni-axial loading (Figure 2-10). To consider this phenomenon, the numerical analysis should be carried out under finite strain conditions and large deformation. However, a small deformation option can still be used to find the results for very small loading conditions. Figure 2-10: Undulation in the yarn (Boisse et al. 1997) Material constitutive models can be another source of nonlinearity. As mentioned in the material properties section, it has been proven that transverse young modulus can be a nonlinear function of both longitudinal and transverse strain. The longitudinal modulus can also be nonlinear for particular type of fibres. There have, however, been studies in the literature that consider material constitutive linearity and have obtained acceptable results for a specific range of deformation (Peng and Cao 2002). 36 2.8 Finite element solver Generally, in a finite element solver an appropriate integration method must be selected. The selected integrator is usually chosen based on the nature of the problem, simplicity level that is expect in formulation convergence criteria and computational cost tradeoffs. Essentially, a finite element analysis can be broken down into increments, in which each increment shows the configuration of the system at a portion of loading. The total load is exerted on the system step by step, starting from zero (or any reference condition) and ending with the final value of the load. The configurations of the system from the start of the analysis to the last step are known as increments. Splitting the answer into increments is usually done due to a variety of reasons; For example, having a better insight about the behaviour of structure as a function of applied load, helping the convergence and avoiding unusual answers that frequently arise for nonlinear differential equations if solved in one increment, etc. Normally, the integrators for the finite element method can be separated into two major groups: Explicit and Implicit. The difference between these two groups is their incremental approach for finding the solution in the next step based on the configuration of the current step. Finite element packages usually offer the option of conducting analysis using either of these integrators. The first approach is known as “Purely Incremental,” “Predictive” or “Explicit” method whereby the configuration of the system in the next step is extrapolated based on the current data. In this method, increments are usually chosen to be very small in order to reduce the amount of error, therefore, a larger number of increments should be applied because there are no iteration increments. In Explicit method there is no strong mechanism to check the results and, therefore, there is usually an amount of drifting error, which means the system is not at equilibrium at the end of analysis (Figure 2-11). 37 Figure 2-11: Drifting error in explicit finite element (Felippa 2001) The second approach is usually referred to as “Corrective,” “Predictive Corrective” or “Implicit” method. Implicit method is based on the explicit method and the only difference is the addition of an extra procedure for evaluating the accuracy of the answer. In other words, if the answer does not lead to equilibrium in the system, some extra iteration (usually using Newton- Raphson method) will be applied to make the current configuration as close to the equilibrium state as possible (Figure 2-12). Because of the existence of a controlling mechanism for enforcing equilibrium in this system, larger step sizes can be chosen compared to the explicit method, however, Newton-Raphson iterations in each increment need calculation of Stiffness matrix, which is computationally expensive. Moreover, there is a possibility that the Newton- Raphson method does not converge. More details about the Finite Element Integrators used in the Abaqus package and specified modifications to make them more suitable for the woven unit cell deformation analysis is presented in the following chapter. 38 Figure 2-12: Implicit finite element procedure, prediction followed by correction (Felippa 2001). 2.9 Constitutive model for dry fabrics Finite strain analysis of woven fabrics most often utilizes rate constitutive (hypo-elastic) equations. The problem, however, is that these formulations use Green Naghdi (Green and Naghdi 1965; Dienes 1979) or Jaumann (1911) objective derivative, which cannot correctly describe the finite strains of fibrous materials (Badel et al. 2007). The alternative approach suggested by Badel et al. (2008a and c) is using a different objective derivative defined from the rotation of fibre frame. Moreover, they proved that using the fibre rotation objective derivative is the only approach which can be expected to result in accurate finite element models for fibrous materials. In the rate constitutive equations, a stress rate ߪ is related to the strain rate ܦ by a constitutive tensor ܥ that is a forth order tensor. To avoid rigid body rotations that can affect stress state, the derivative ߪ, called “objective derivative,” is the derivative for an observer who is fixed with the material. There are different ways to define the objective derivative. Most common objective derivatives are rotational objective derivatives that correspond to a rotation tensor, ܳ characterizing the rotation of the material. So, we have: ߪ ൌ ܥ: ܦ (2-50) with: 39 ߪ ൌ ܳ ቂ ௗ ௗ௧ ሺ்ܳ. ߪ ܳሻቃ ்ܳ (2-51) where the definition for ܳ depends on the selected objective derivative. Finite element procedures use the rate constitutive equations to update the stress in each increment after the nodal displacement and corresponding strain values are computed. The integration of equation (2-51) in each time increment can be done using the formula of Hughes and Winget (1980). Assuming the time increment to be Δݐ ൌ ݐାଵ െ ݐ: ሾߪାଵሿశభ ൌ ሾߪ ሿ ൣܥ ାଵ/ଶ൧ శభ/మሾΔߝሿశభ/మ (2-52) and: ሾΔߝሿశభ/మ ൌ ൣܦ ାଵ/ଶ൧ శభ/మ Δݐ (2-53) where matrix ሾܵሿ denotes the components of tensor ܵ in the rotated frame with unit vectors ሼ݁ ሽ at the n-th time increment. The objective derivative based on the fibre rotation can be found using ܳ ൌΦ where Φ is the rotation of the fibre. If the frame of fibre is presented by in which ଵ is the unit vector in the frame direction, then assuming the initial position of fibre is ൌ we have (Badel et al. 2008a and 2009): ଵ ൌ ۴.భ బ ห۴.భ బห (2-54) where ۴ is the deformation gradient tensor. The other basis vectors for this orthonormal frame can be obtained by the material transformation of ଶ: ଶ ൌ ۴.మିሺ۴.మ.భሻ భ |۴.మିሺ۴.మ.భሻ భ| (2-55) ଷ ൌ ଵ ൈ ଶ (2-56) Therefore, the rotation tensor Φ from the initial frame to the fibre frame is derived using: Φ ൌ ൫. ൯ ٔ (2-57) Eventually, the stress update formula appropriate for finite element modeling can be defined as: ሾߪାଵሿశభ ൌ ሾߪ ሿ ൣܥ ାଵ/ଶ൧ శభ/మሾΔߝሿశభ/మ (2-58) This subject has been discussed in more details in the subsequent chapter. 40 2.10 Summary The general procedures for creating geometry of a plain weave fabric, assigning material properties, defining the loading boundary condition, as well as the special type of boundary conditions known as periodic unit cell, were described in this chapter. The presented analysis framework helps analyze accurate representative models of woven fabrics under homogenous loadings. It is essential to apply non-linear finite element procedures in this type of analysis. Keeping in mind nonlinearity of the models, both predictive (explicit) and predictive-corrective (implicit) finite element solvers can be justifiable options depending on a given application. 41 Chapter 3: Creating FE model in Abaqus This chapter presents details of the finite model used in this thesis. The model is then used for the study of uncertainties at meso-level. It begins with the geometrical modeling of yarn as the basic element at meso-level. Next, an assembly of yarns are used for producing a unit cell model along with its boundary conditions. A special material constitutive model is described for implementation in Abaqus using user-defined subroutines. Finally, using an inverse method, material properties of yarns are identified to follow experimental measurements on a plain weave fabric. 3.1 Reference unit cell The selected fabric type in this study is a balanced plain weave, which is one of the most common types of woven fabrics with numerous applications in aerospace and transportation industries. Relations for defining constitutive equations of the fabric’s yarn material are found via an inverse method. This method is applied to construct a constitutive equation that forces the basic unit cell to follow the published experimental results for uni-axial, bi-axial (Buet-Gautier and Boisse 2001) and shear frame tests (Cao et al. 2008). Although, the unit cell dimensions for the above shear and axial tests have been different, both experimental works use glass fibres. Therefore, almost identical material properties are used for axial and shear loading cases. Indeed, slight differences in the material properties should not have a significant effect on the final goal of this research: identifying important uncertainty factors at meso-level under individual deformation modes. The geometrical definition used by Peng and Cao (2000), which was originally based on the study by McBride and Chen (1997), has been selected here (Figure 2-4a). The crucial parameters for a yarn are the yarn width (w), yarn spacing length (s), which is identical to the total length of the unit cell and yarn height (h). These parameters are geometrically shown in Figure 3-1. 42 Figure 3-1: Geometrical parameters for defining the basic yarn. 3.2 CAD modeling of yarns To produce the geometry of the yarn, a parametric CAD model has been made using the CATIA software. This model takes four independent input parameters to define the geometry of the unit cell (Figure 3-2). Yarn length, width, height and alignment angle are inserted to the formulas as an input. The cross sectional shape is defined by setting the first three parameters into formulas given by McBride and Chen (1997) (eqs. 2-1 to 2-4). These four formulas are used to locate 19 points on a plane. The cross section shape is defined by connecting these points to each other (Figure 3-3). In the same manner, the opposite cross section and yarn axis line has been generated. Consequently, in the final step the “Multi-section Solid” tool in CATIA is used to connect two sections using the equation of yarn axis as a guide (Figure 3-4). Yarn axis itself is also made from the connection of 11 points in a plane normal to the section of yarn. 43 Figure 3-2: Dialog box for defining the geometrical parameters of the yarn Figure 3-3: Defining the parametric cross section. 44 Figure 3-4: Creating solid yarn using two end cross sections and the connecting line. 3.3 Defining model in Abaqus FE software The yarn model created in the previous section is imported into Abaqus/CAE and assembled to create the unit cell. The assembling step in Abaqus defines the initial position of the part and does not imply any bonding on their deformations and translations. To define the internal reactions between the yarns, their interactions should be defined. In essence, defining interactions in Abaqus relates movements of a node set or reference point to other node sets or reference points. This can be via constraints to make nodes have the same displacement/rotation in a specific direction (coupling), assigning the displacement of one node set or reference point to be related to displacement of another node set or reference point through an arbitrary equation, or defining contact conditions between two surfaces or node sets. There is also possibility of defining kinematic connectors (hinge, slider, etc.) between nodes and reference points. Because of different kinematic conditions for applying the loading on the unit cell in axial and shear tests, separate sets of interactions may be defined for them. However, it should be noted that both of these interaction sets are designed in a way to satisfy the periodic boundary condition in applying the macroscopic field displacement on the unit cell. 45 Contact conditions are the same for both loading cases. Contact is defined between the yarn surfaces from the initial configuration. Tangential behavior at contact surfaces is defined using the penalty method with the friction coefficient. The tangential behavior at the contact surfaces depends on the value of normal pressure between two surfaces. Classical isotropic Coulomb friction implies no relative motion occurs if the equivalent frictional stress (߬) is less than ߬௧ ൌ ߤ , where ߤ is the friction coefficient and is the contact pressure. ߬ is defined to be ߬ ൌ ඥ߬ଵ ଶ ߬ଶ ଶ (3-1) where ߬ଵ,ଶ are frictional stresses at contact surface in two normal directions. With the penalty contact algorithm in Abaqus, the penalty stiffness term between two nodes with no relative displacement in contact is used to keep the relative displacements zero. Therefore, the relative motion in the absence of slip is equal to the friction force divided by the penalty stiffness, which is very close to zero value and ignorable (SIMULIA 2008). 3.3.1 Kinematic conditions for axial loading The periodicity boundary condition implies that the side planes should be kept plane and normal to the mid-plane. For axial loading there is an advantage, because during the analysis the planes stay parallel or normal to the global coordinate systems. Therefore, the “Equation” constraint in the Interaction module of Abaqus can be easily used for defining the periodic boundary condition. The four side surfaces are defined as sets, which are actually the boundaries and then the degree of freedom of nodes on each surface (set) is constrained by an equation to the degree of freedom of a reference point. This mechanical constraint makes all the nodes on the plane have the same displacement in the specified direction, which means they remain plane. To illustrate this procedure more clearly, it is shown graphically in Figures 3-5 and 3-6 for the boundary surface with the normal vector ො݊ ൌ ሼ1,0,0ሽ with respect to the global coordinate system. First, the boundary has been selected (highlighted part in Figure 3-5) and the nodes are assigned to a set named xP_Surf, which refers to the boundary surface with positive x component in its normal vector. Then, a reference point has been created separately and called RP-1. It should be noted that the reference point is defined as a point that does not belong to the main geometry of the unit cell. In other words, it is an external reference for adjusting movements on the yarn’s section surface. 46 Figure 3-5: Node sets created for defining periodic boundary condition on one of the unit cell surfaces. The arrows in Figure 3-5 shows the desired direction for constraining xP_Surf ,which is actually the positive x direction in the global coordinate system. Constraint is applied via defining the following equation for the above mentioned node sets. ݑଵ,௫_௦௨ െ ݑଵ,ோିଵ ൌ 0 (3-2) where ݑଵ,௫_௦௨ and ݑଵ,ோ_ଵ refer to the displacement in x direction for nodes in xP_Surf and RP-1, respectively. In Abaqus/CAE, this can be applied by the inputs shown in Figure 3-6 in the “Equation” creator dialog box of “Interaction module.” A similar procedure is also applied to the three other surrounding surfaces to ensure appropriate boundary condition. 47 Figure 3-6: Defining the constraint for surrounding surfaces in Abaqus/CAE After defining interactions in the unit cell to simulate the periodicity, the next step is to apply the loading boundary conditions, which are of displacement nature. Indeed, this is the actual boundary condition that applies to the fabric in a tension test. Figures 2-8 (a,b) show this boundary conditions schematically. However, for the sake of simplicity in applying boundary conditions and reading reaction forces, the prescribed displacements are on the reference point of each surface. This means that instead of assigning displacement for each node on that surface, displacement is assigned only to the reference point; however all the nodes on the surface will comply with the movements of that reference point. Overall, this means that in the post processing of the finite element simulation, there is no need to add up the reaction force from each node to find the total value of reaction forces. Hence, this is a simplification without losing any details of the problem. 3.3.2 Kinematic conditions for shear loading Based on the nature of the loading in the shear frame test, separate types of interactions and boundary conditions have been applied. This time, the peripheral lines do not stay normal to the global coordinates and, therefore, the equations defined for constraining the nodes on the yarn surfaces should be updated by rotation of the yarns. This can be done via defining user defined subroutines for constraining yarns. Another method, which is applied here, is setting up a frame composed of four rigid links and then forcing the nodes on the peripheral surfaces to follow a specific degree of freedom with respect to the coordinate system attached to the rigid bars. Figure 3-7 shows the configuration of the rigid bars and the coordinate systems attached to them, 48 before and after loading. Here the bars are pinned to each other at the corner points and they are only allowed rotational degrees of freedom, thus allowing them to rotate in the plane of deformation – the rest are constrained by boundary conditions to avoid rigid body motion. This type of boundary condition is shown schematically in Figure 2-8c and is applied to the rigid frame and then transferred to the unit cell using couplings. The coupling constraint is useful when a group of coupling nodes is constrained to the rigid body motion of a single node. Figure 3-7: Rigid frame and the coordinate system attached to each link before and after deformation. There are two other boundary conditions applied to the unit cell to make it follow the kinematic movement of yarns in the picture frame test. Coinciding nodes at the corner of the unit cell are pinned to each other (Figure 3-8a). Moreover, the nodes at the corner lines of the unit cell are constrained using an equation to remain perpendicular to the plane of deformation (Figure 3-8b). To keep the corner lines perpendicular, the displacement of all nodes are assigned to be the same in x and z directions. This is done again by defining the lines as node sets and assigning their degree of freedom to the degree of freedom of a reference point using an equation similar to the equation (3-2). 49 (a) (b) Figure 3-8: (a) Pinned nodes at the corner of the unit cell. (b) Node sets and reference points for keeping the corner lines perpendicular to the deformation plane. So far, the interactions between the bars and yarns of the unit cell have been defined. The next step is constructing interaction between the unit cell and the rigid frame. The “Coupling” constraint of Abaqus/CAE “interaction” module is used to define the periodic boundary condition on the unit cell. For each yarn, the reference node attached to the rigid bar beside that yarn is selected and then three boundary surfaces of the yarn, including two at each end and one side surface, are coupled in the appropriate direction to the movement of the reference node. Coupling here is quite similar to the equation constraint defined for the axial loading, however, the advantage is that because the reference node is attached to the rigid frame, updating of the coordinate system is done automatically in the solver based on the rotation of the rigid frame. Figure 3-9 shows the coupling used for one of the yarns. For the rigid bar shown with black color, the red circle indicates its reference node. The displacement of the nodes located on the surface highlighted in Figure 3-9 (a) is coupled in z-direction to the displacement of the reference node. This means that the highlighted surface acts like a rigid surface in the z-direction. Figure 3- 9 (b) shows the same subject for that yarn, however, this time end sections are selected and coupled to the x-direction. 50 Figure 3-9: Coupling of yarn (a) side surfce (b) end surfaces. Lastly, to make the boundary condition on the frame easier to apply and the reaction forces in the post processing easier to read, another interaction using “Equation” constraint is generated to make the corners of the rigid frame move with the reference point in the plane of deformation. Here, the boundary conditions can be applied to the reference point, which makes the reading of reaction force at post processing easier (Figure 3-10). Figure 3-10: Corners of the rigid frame constrained to the movement of a reference point. 3.4 Material modeling As previously mentioned, there are two ways for defining material properties for yarns in meso-level. There can be either the use homogenization results from micro scale to find the 51 equivalent material properties for the yarn or there can be an appropriate constitutive model for the yarn material from its physical properties. This model should be able to capture the fabric material response from actual test results. 3.4.1 Homogenized material properties Using micro to meso level homogenization, Peng and Cao (2002) applied the concept of homogenization in heterogeneous media (Bakhvalov and Panasenko 1989) to find the effective material properties of E-Glass/PP yarns from properties of pure E-Glass and Polypropylene. Defining a customized finite element procedure for the homogenization step they predicted the material properties of yarns as shown in Table 3-1. Table 3-1: Material properties of E-Glass/PP yarns predicted by Peng and Cao (2002) Longitudinal Stiffness (El) 51.92 (GPa) Transverse Stiffness (Et) 21.97 (GPa) Longitudinal Poisson’s Ratio (νlt) 0.2489 Transverse Poisson’s Ratio (νtt) 0.2143 Longitudinal Shear Stiffness (Glt) 8.856 (GPa) Transverse Shear Stiffness (Gtt) 6.250 (GPa) A unit cell model based on the boundary condition details given in this chapter and the material properties from Table 3-1 is analyzed under three basic loading modes. Three basic dimensional values of the yarns (Figure 3-1) are s=5.14 mm, w=3.72 mm and h=0.39 mm. This makes the model match the one used by Peng and Cao (2002) and provides a measure for validating our finite element model. Running the model using Abaqus/Standard 6.8-1 provided the following results that show a close correspondence to the results by Peng and Cao (2002). 52 Figure 3-11: Validation of the refernece finite element model 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.05 0.1 0.15 0.2 R ea ct io n F or ce ( N ) Displacement (mm)(a) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 0.05 0.1 0.15 0.2 R ea ct io n F or ce ( N ) Displacement (mm)(b) 0.0E+00 5.0E+04 1.0E+05 1.5E+05 2.0E+05 0 10 20 30 40 N or m al iz ed F or ce (N /m m 2) Shear angle (degree) Peng et al. (2004) Presented Work (c) 53 3.4.2 Postulating a suitable constitutive model in meso-level It has been proven that material models using homogenizations and averaging methods can capture and predict the fabric response reasonably well. However, as previously mentioned, there is a more advanced constitutive material model to characterize more specific behavior of fabrics (Gasser et al. 2000; Boisse et al. 2007 and Badel et al. 2008a). Mechanical behavior of yarns is very particular due to the fibrous nature of the yarns. This means that each yarn is made of a very large number of fibers that have a very small diameter so that they cannot bear any compressive force. Moreover, dry fabrics in yarns can easily slide on each other, which results in a substantial decrease in shear modulus and transverse young modulus. On the other hand, fabrics can have considerable stiffness during tension. The notable amount of difference between the stiffness in fibre direction and other stiffness values suggests that the appropriate constitutive model for yarns should be presented at the local frame coordinates of the fibres. In other words, the objective derivatives used by default in Abaqus (Green Naghdi or Jaumann) can lead to unsatisfactory results for the analysis of fibres (Badel et al. 2008a). As the fibre goes through deformation, the frame of rotation used by the software (e.g. Green Naghdi) rotates, however, because it is not obliged to always follow the actual fibre direction, it might have directions that do not coincide with the rotated fibre. This has been illustrated in Figure 3-12, where fi=1,2 shows the fibre frame and ei=1,2 shows any other possible frame, which is not obliged to follow the fibre direction. Keeping in mind that material properties are presented in the frame of fibre, as well as the huge difference of stiffness in the two local directions, the conventional rotation algorithm of the frame cannot be a superior choice. Figure 3-12: Comparison of Fibre Frame and other possible frame of rotation in 2D; (a) before deformation frames coincide, (b) the separation of frames after deformation. 54 The modification of the finite element procedure to handle the analysis in the correct frame can be done using two different methods (Badel et al 2008a): “Modified Green-Naghdi” and “Fibre Frame” formulations. The underlying mathematical procedures were briefly described in Chapter 2. In simple terms, this means that either the constitutive model should be presented and calculated in the fibre frame or material properties should be transformed to the current frame. In this project, Abaqus 6.8-1 is used, which employs the Green Naghdi stress rate in explicit calculations for solid (continuum) elements and Jaumann’s for Implicit (Abaqus/Standard). Nevertheless, differences between these two stress rates are significant only if finite rotation of a material point is accompanied by finite shear (SIMULIA 2008), which does not happen in elements for the studied loading types. Moreover, the rotation tensor for the Green-Naghdi stress rate is readily provided in every step as a parameter in user defined subroutines. Thus, the model is developed assuming the stress rate is updated by the Green-Naghdi objective derivative, and the rotation from the Green-Naghid frame is modified in order to find the actual rotation of material directions. The first method, which is called the “Modified Green-Naghdi,” begins with writing the stiffness matrix in the frame of fibre and then transforms it to the frame of software calculations. Thus, the only required modification is the rotating of the stiffness matrix. In mathematical terms: ሾሿ ൌ ሾሿ ሾሿ ሾ′ሿ (3-3) where ሾሿ, is the stiffness matrix in Green-Naghdi and the fibre frame respectively. ሾሿ and ሾ′ሿ are the base change matrices described in Appendix A. It must be emphasized that because of non tensorial nature of ሾሿ in the software (Abaqus uses a 6 by 6 matrix to store stiffness tensor), ሾ′ሿ is neither transpose of ሾሿ nor its inverse. Indeed ሾሿ is a singular matrix (|| ൌ 0). However, it can be proven that ሾሿ. ሾ′ሿ ൌ ሾሿ; where ሾሿ is the identity matrix. The second method is called “Fibre Frame,” which implies rotating stress and strain values from the software frame to the fibre frame, calculating the new stress values and then providing the software with new stress values, which are then transformed back into the software frame. In other words first we calculate: ሼሽ ൌ ሾ′ሿ ሼሽ (3-4a) ሼ∆ࢿሽ ൌ ሾ′ሿ ሼ∆ࢿሽ (3-4b) 55 Then the stress strain calculation for new stress values are done in the fibre frame: ሼାଵሽ ൌ ሼሽ ൣାଵ/ଶ൧ ൛∆ࢿ ାଵ/ଶൟ (3-5) However, the default frame for the software is different from the one that was done in the calculations, therefore, new stress components should be transformed back then sent back to the software: ሼାଵሽ ൌ ሾሿ ሼାଵሽ (3-6) 3.5 Implementing the material constitutive model in Abaqus Implementation of the material model described above in the Abaqus package can be done via defining a subroutine for material properties to express a particular constitutive behaviour of yarns. The analysis integration can also be done in two different manners, which are Explicit/Dynamic and Static/General; where the first method utilizes an explicit (predictive) solver and the second provides an implicit (predictive-corrective) solver. Abaqus accepts user defined subroutines written in FORTRAN programming language. Codes are compiled using Intel® FORTRAN compiler and linked to the program to be analyzed by Abaqus. Because of differences in procedure conducted by Implicit and Explicit solvers the subroutines accepted by Abaqus for these two solvers are different. 3.5.1 Explicit integrator In Abaqus, the explicit procedure is based on the implementation of a fully predictive integration rule together with the use of diagonal (“lumped”) element mass matrices. The equations of motion for the body are integrated using the explicit central-difference integration rule, which can be simply described as (SIMULIA 2008): ାଵ ൌ ∆ݐାଵ ሶ ାభ మ (3-7) ሶ ାభ మ ൌ ሶ ିభ మ ∆௧శభା ∆௧ ଶ ሷ (3-8) where ୨ and ∆ݐ୨ are the displacement vector and time incrimination in j-th increment, respectively. ሶ ୨ and ሷ ୨ are the first and second derivatives of the displacement vector which 56 obviously are velocity and accelerations. Moreover, half values for incriminations show the central difference method used in integration. Accelerations, which are the essential part to start the analysis are calculated via: ܝሷ ൌ ିሺ െ ሻ (3-9) where is the mass matrix and and are the applied load vector and internal force vector. For the sake of simplicity in computation of inverse mass matrix, Abaqus uses a lumped mass matrix. Internal forces are assembled from contributions from the individual elements. Thus, as it can be seen in these formulations in the explicit method no calculation of global stiffness matrix and its tangent matrix is required, which makes calculations much faster compared to the implicit integrator. However, to assure the stability, safe time incrimination should be selected. This incrimination is usually very small and makes the analysis go through many iterations. Nevertheless, the most prominent disadvantage of the explicit solver is that there are no mechanisms to guarantee the correctness of the answer and the final equilibrium state. 3.5.2 Explicit subroutine For the Explicit solver, a VUMAT subroutine should be used. This subroutine has a number of variables, which are assumed to be input values (read only parameters as called by the Abaqus Documentation) and some other variables, which are updated in the subroutine (write only parameters). For this analysis, the read values in the subroutine are strain increments, right stretch tensor and deformation gradient tensor. The parameters that are updated by the subroutine are the new stresses. Based on the needs of the analysis, there may be the possibility of using some state variables for recording internal parameters that are calculated in the code, but are not among the subroutine conventional parameters. Stress in the fibre direction and some failure criteria can be one of the possible reasons for using state variables. The general algorithm for this subroutine by using the Modified Green Naghdi or the Fibre Frame is shown in the following flowchart (Figure 3-13) and described below. 57 Figure 3-13: Flowchart for programming the user-defined subroutine in explicit solver The parameters provided for the explicit analysis are ሾ܃ሿబ , ሾ܃ାଵሿబ ; ሾ۴ሿబ , ሾ۴ାଵሿబ ; ሾોሿ and ሾ∆ઽሿశభ/మ where ሾ܃ሿ is the right stretch tensor, ሾ۴ሿ is the deformation gradient 58 tensor, ሾોሿ the stress tensor and ሾ∆ઽሿ represents the increment in the strain during current step of the analysis. The only output value for this subroutine is the stress value ሾોାଵሿశభ at the end of increment. It should be noted that the superscript for matrices shows their current increment number and subscript is the coordinate system that tensor is defined in. The input and output of this subroutine subscripts are ݁, which means they are in Green Naghdi (Abaqus default) frame at the j-th step. First, the rotation matrix used by the software to rotate the coordinate should be found; it is not readily given by the subroutine and should be calculated recalling that ۴ ൌ ܀ ܃, where ܀ is the rotation tensor: ሾ܀ାଵሿబ ൌ ሾ۴ ାଵሿబ. ሾ܃ ାଵሿబ ିଵ (3-10) Then the current work basis can be defined from the initial coordinate system by: ሼ ାଵሽ ൌ ሾ܀ାଵሿబ. ሼ ሽ (3-11) So far, the direction of the material calculations that Abaqus is using have been found. The next step is to find the fibre frame unit vector directions. Using (eqs. 2-54 to 2-56), the fibre frame unit vectors ሼାଵሽ can be calculated. Then the rotation matrix between ሼାଵሽ and ሼାଵሽ frames is calculated (Badel et al. 2008a): દ ൌ ٔ ൌ ൫൯ ٔ ൌ ൫ . ൯ ٔ (3-12) Subsequently, the transformation matrices ሾሿ and ሾ′ሿ are constructed and ሾሿ are generated from ሾሿ by applying the transformation shown in (3-3). In the case of nonlinear material properties, the strain in fibre directions should also be available to be accounted in the ሾሿ. Hence, there is one more transformation to be done, which is rotating strains from ሼାଵሽ frame to ሼାଵሽ frame. This is because strain dependant material properties, like longitudinal stress- strain relations and crushing between the yarns, are defined in the direction of fibre or its normal directions. This transformation is as follows. ۴ and ܃ are calculated at the end of each time increment so it is needed to work at the end of step configuration instead of mid step configuration. However, the problem is that the stress increments are given at mid-step configuration. Let ∆܀ ൌ ܀ାଵ. ൫܀ାଵ/ଶ൯ T be the rotation tensor between the end and mid-step configurations and ܂∆܀ the corresponding transformation such that for any two-dimensional tensor ܁ we have ܂∆܀ሺ܁ሻ ൌ ∆܀ . ܁ . ∆܀T. Then the following relation can be used (Badel et al 2008a): 59 ∆ݐ ൣ۲ାଵ/ଶ൧ శభ/మ ൌ ∆ݐ ൣ܂∆܀൫۲ାଵ/ଶ൯൧శభ/మ (3-13) while ∆ݐ ൣ۲ାଵ/ଶ൧ శభ/మ ൌ ሾ∆ઽሿశభ/మ hence, ሾ∆ઽሿశభ/మ ൌ ሾ ∆ܶ܀ሺ∆ઽሻሿశభ (3-14) This enables a direct use of the strain increment provided by the code. A similar logic can be applied to the constitutive matrix as a representation of a fourth-order tensor. For using Fibre Frame method the first parts of the procedures up to equation (3-12), which is includes calculating દ (rotation matrix) is the same. The only difference is that here, instead of rotating the constitutive matrix to the working frame, the stress and strain increment tensors are rotated to the frame of fibre for calculating the new stress values: ሾ܂∆܀ሺ∆ઽሻሿశభ ൌ ሾદሿశభ T ሾ܂∆܀ሺ∆ઽሻሿశభ ሾદሿశభ (3-15) ሾોାଵሿశభ ൌ ሾો ሿ ൣ ∆ܶ܀൫۱ାଵ/ଶ൯൧శభ . ሾ ∆ܶ܀ሺ∆ઽሻሿశభ (3-16) The procedure is complete by rotating the strain tensor back to the Abaqus working frame using equation (3-6). 3.5.3 Implicit integrator As discussed in the previous chapter, implicit analysis starts with a predictive increment followed by a series of iterations within each increment to attain equilibrium. Abaqus uses Newton’s method (modified versions are available as an option), where iterations are done for the calculation of global stiffness matrix and its tangent. Indeed, the matrix calculations regarding the global stiffness matrix and its tangents, as well as the fact that these calculations have to be done more than once during the iterations, makes the steps for implicit analysis run much slower than the steps for explicit analysis. However, the iterations, which are intended to force the system to attain its equilibrium, allows for larger steps to be chosen (∆ݐ୨ ) in the predictive part of the step. To summarize, it is evident that applying the explicit method makes Abaqus do more in faster increments, while using implicit integrator makes the computations for each increment slower, while using a fewer number of increments. Regarding the subroutine programming for 60 the constitutive model of these two analysis types, the implicit integrator needs to be provided with updated values of stress and stiffness matrix (tangent stiffness matrix is calculated internally by Abaqus). 3.5.4 Implicit subroutine The programming for implicit formulation is similar to explicit formulation, because the first level in each step done in explicit and implicit is the predictive level as described above. However, as previously mentioned, it needs the general stiffness matrix of the system in its corrective iterations. Therefore, the stiffness matrix of each element must also be provided for this subroutine and extra attention should be paid to define a stiffness matrix, which is consistent with the defined stress and strain relationship, otherwise the convergence would be violated. Regarding the fact that the stress-strain relations are readily defined in the frame of the fibre, the stiffness matrix should go though the transformation as well. Keeping in mind the fact that in Abaqus stiffness matrix is defined as a 6 by 6 matrix and shear strains for Implicit formulation are engineering shear strains ߛ and ߛ ൌ 2 ߝ; the transformation for stiffness can be derived by starting with: ሼሽ ൌ ሾሿ ሼࢽሽ (3-17) in which ሼࢽሽ ൌ ሾሿሼࢿሽ and ሼࢽሽ ൌ ሼߝଵଵ, ߝଶଶ, ߝଷଷ, ߛଵଶ, ߛଵଷ, ߛଶଷሽ ሾሿ ൌ ሾሿଷൈଷ ሾሿଷൈଷ ሾሿଷൈଷ 2 ሾሿଷൈଷ ൨ And from (3-4) we can write: ሾ′ሿ ሼሽ ൌ ሾሿ ሾሿ ሾ′ሿሼࢿሽ (3-18) So, one may try to find the stiffness matrix in frame of software by inverting ሾ′ሿ. However, this matrix is not a regular tensorial transformation and, moreover, it is not invertible. On the other hand, there is: ሼሽ ൌ ሾሿ ሼሽ ൌ ሾሿ ሾ′ሿ ሼሽ It is known that ሼሽ has a real answer, which is unique; Thus: ሾሿ ሾ′ሿ ൌ ሾሿ (3-19) where ሾሿ is the identity matrix. Therefore, by pre-multiplying two matrices there is: 61 ሼሽ ൌ ሾሿ ሾሿ ሾሿ ሾ′ሿ ሾሿିଵ ሼࢽሽ (3-20) It can be written in the following compacted form: ሼሽ ൌ ሾሿ ሼࢽሽ (3-21) ሾሿ ൌ ሾሿ ሾሿ ሾሿ ሾ′ሿ ሾሿିଵ (3-22) In the Implicit integrator there is great advantage for using the Modified Green-Naghdi method over the Fibre Frame method in subroutine definition. Indeed, the stiffness matrix in the default frame has to be provided for Abaqus regardless of using either the Modified Green- Naghdi or the Fibre Frame method. This means that the calculations must be done for the stiffness matrix rotation anyway. By using the Fibre Frame method, calculations to rotate the stress and strain tensor matrices must be done, as do the stress-strain calculations, and then they need to be transformed back to the default frame of Abaqus. In addition, the transformations on the stiffness matrix should be done for iterative process. On the other hand, the Modified Green- Naghdi method only needs the transforming stiffness matrix, which can be used both for stress- strain calculations and iterations. Hence, there is no need to do extra transformations on the stress and strain vectors to the Fibre Frame and inverse transformations to the work frame. At most, under the assumption of nonlinear material behavior, the strains must be transformed to transverse and longitudinal strains in the fibre direction. Badel et al. (2008a) and Badel et al. (2009) used the explicit integrator of Abaqus for shear frame test modeling of a unit cell. They implemented both the Fibre Frame and the Modified Green-Naghdi type of constitutive models and concluded that calculations using the Fibre Frame method are more stable, whereas the Modified Green-Naghdi formulation crashed during the analysis. However, this scenario is completely different using implicit integrators. The iteration mechanism in the implicit solvers assures that equilibrium is satisfied for the model. Moreover, doing the stress-strain calculations during an increment in one frame and iterations in another frame increase the risk of diverging, because of possible numerical round off errors. These errors might make the calculated stress values not match the global stiffness matrix. To summarize, it can be argued that using the Modified Green-Naghdi method in user defined material subroutine for constitutive material behavior of fabrics along with an implicit solver can provide a powerful, reliable and fast tool for finite element analysis of unit cells at meso-level. Apart from the above-mentioned facts, programming strategy for the implicit integrator using the Modified Green Naghdi and the Fibre Rotation method is very similar to what was 62 described for the explicit solver and can be illustrated in the flowchart as seen in Figure 3-14. There are, however, a couple of slight differences in the implicit subroutine: 1. Some of the tensors in subroutine variables are given in the rotated coordinate system, instead of the initial coordinate system. Hence, appropriate transformations should be applied to them. 2. The right stretch tensor is not provided readily and it should be calculated from the deformation gradient tensor. ሾܷሿଶ ൌ ሾܨሿ். ሾܨሿ, so ሾܷሿ can be calculated from the matrix square root of ሾܷሿଶ. It is worth noting that, ሾܷሿ cannot be calculated by simply finding the square root of each array in ሾܷሿଶ. A more detailed description on calculating the square root matrix is given in Appendix B. 63 Figure 3-14: Flow chart for programming the user-defined subroutine in implicit solver 64 3.6 Postulating material model constants Unfortunately, a comprehensive set of material properties consistent with the presented material constitutive model (rotation with the frame of fibre), along with its validation with experimental axial tension and shear tests, could not be found in the literature. Hence, an inverse method is applied to find the stress-strain relationship for constitutive models from results of tension and shear frame tests on similar materials. In other words, a postulation is made on the stiffness functions and model coefficients are determined by iterations. Acceptable results are those that make the response from the finite element simulations under all loading modes (uni- axial, bi-axial and shear frame test) close to the experiments. The experimental results were taken from Buet-Gautier and Boisse (2001) for tension tests and Cao et al. (2008) for shear tests. Fabrics are balanced plain weave with yarns made of glass fibres. Although materials are not taken from identical studies and the meso-level dimensions in these tests are different, an identical set of material properties for their yarns has been proposed and applied here. The main purpose of this research is identifying important uncertainty factors under each individual deformation mode. 3.6.1 Material model It is of high priority to have a model, which is both simple and general enough to be adapted to new materials. The proposed material model consists of two basic stiffness values, which are: 1. Longitudinal Stiffness (ܧଵଵ). 2. Transverse Stiffness (ܧଶଶ and ܧଷଷ). The remaining material properties in the stiffness matrix are zero or quasi-zero as discussed in Chapter 2. The stiffness for the yarn material is constructed as: ܥ ൌ ۏ ێ ێ ێ ێ ۍ ܧଵଵ 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 0 0 ܩ 0 0 0 ܩے ۑ ۑ ۑ ۑ ې (3-23) As it can be observed, besides the longitudinal and transverse stiffness values, shear modulus (ܩ) is the only non-zero value in this matrix. However, it should be very small compared to the 65 normal stiffness values. Its purpose is solely to assure non-singularity of stiffness matrix and stability in the numerical methods (Gasser et al. 2000). 3.6.2 Longitudinal yarn behavior It has been assumed that longitudinal stiffness is a function of longitudinal strain only. If the function is started as ܧଵଵ ൌ ܧଵଵሺߝଵଵሻ, the way to find ܧଵଵሺߝଵଵሻ can be via curve fitting using an inverse method that utilizes the results from finite element analysis and after a series of trial and errors finds the appropriate function which makes those results match experiments. There are a variety of functions that can be proposed for this purpose. Keeping in mind a minimum number of model factors, a piece-wise linear axial stress trend was implied by choosing: ܧଵଵሺߝଵଵሻ ൌ ൞ ܧଵଵ,ଵ ߝଵଵ ߝଵଵ,ଵ ܧଵଵ,ଶ ߝଵଵ,ଵ ߝଵଵ ൏ ߝଵଵ,ଶ… … … … ܧଵଵ, ߝଵଵ,ିଵ ൏ ߝଵଵ (3-24) 3.6.3 Transversal yarn behavior The transversal behavior of yarns has been described by the crushing law, wherein the transversal stiffness is related to some other factors, which usually depend on contact conditions between the yarns. This transverse stiffness can be described as a function of contact force (Kawabata et al. 1973a-c) or strain along the contact normal and longitudinal strain (Gasser et al. 2000) or other more complicated formulations (Badel et al. 2008b-c). The basic function that has been chosen here is the one suggested by Gasser et al. (2000) and is defined as: ܧଶଶ ൌ ܧ|ߝଶଶ |ߝଵଵ ܧଷଷ ൌ ܧ|ߝଷଷ |ߝଵଵ (3-25) where ܧ, ݊ and ݉ are model constants that relate to the material properties and can be identified by using an inverse method. However, it has been seen in previous works that a value of ݊ ؆ 1.0 and ݉ ؆ 2.0 can reflect crushing behavior quite well for glass fabrics (Badel et al. 2007; Badel et al. 2008a). Therefore, it can be presumed that these two parameter values before starting the inverse method are (݊ ൌ 1.0 and ݉ ൌ 2.0). The rest of inverse identification for the crushing phenomenon is about fitting ܧ. 66 3.7 Inverse identification Trial and error methods have been utilized here to find constants of stiffness functions in the previous section, as well as the minimum acceptable shear stiffness value to keep finite element runs stable. Finally, a mesh sensitivity analysis has been exploited to ensure the accuracy of the analysis. First, the axial tension is tested. The experimental data is taken from Buet-Gautier and Boisse (2001) for uni-axial and balanced bi-axial tension tests on balanced plain weave glass fabrics. The reported dimensions for the test samples are: Density (yarn/mm) = 0.22 ; Crimp (%) = 0.4 where ܥݎ݅݉ሺ%ሻ ൌ ݕܽݎ݊ ݈݁݊݃ݐ݄ െ ݈݁݊݃ݐ݄ ݂ ݂ܾܽݎ݅ܿ ݂ݎ݉ ݓ݄݁ݎ݁ ݐ݄݁ ݕܽݎ݊ ݅ݏ ݁ݔݐݎܽܿݐ݁݀ ݂ܾܽݎ݅ܿ ݈݁݊݃ݐ݄ ൈ 100 and from (2-5) the yarn length can be calculated as a function of yarn spacing (s) and height (h). Therefore, starting with the above constants we obtain s = 4.55 (mm) and h = 0.366 (mm). Using the image scale provided for the test samples (Figure 3-15), yarn width is estimated at 4.44 (mm). Figure 3-15: Photo from sample tested by Buet-Gautier and Boisse (2001) After finding the geometrical constants, the unit cell model is created in Abaqus/CAE, where the yarns were generated by a CATIA parametric design. Regarding the fact that stiffness along the yarn axis in the fabric is much higher than other directions, extra attention should be paid in defining the appropriate local material directions on each point of the yarn. Therefore, each yarn is split into different sections by a surface normal to its axis and material properties were defined along a local coordinate system on that section. Figure 3-16 shows separate sections in these c Zero problems problem, hourglass It ha the respo can also mentione the loadi waviness is no con yarns at importan stiffness propertie modulus The and local co oordinate sy Figure 3-16: and very cl in numeric the approa control) is s been obse nse is from be recogniz d in Chapte ng direction . At the sam siderable so contact surf t role in un function fi s from bi-ax is taken to b longitudinal ordinate sys stems lies a Local coord ose to zero al simulation ch develope used (Pian a rved that in axial stiffn ed from the r 2, in uni-a , they unde e time, yarn urce of exte aces are ve i-axial load rst by using ial tests and e 15 MPa. stiffness fu tems define long the lon inate system values of sti s and lead d for buildi nd Chen 19 some prelim ess, wherea mechanism xial loading rgo a type o s in the dire rnal force to ry low. Thu ing. This fa the result shear tests nction after d for assign gitudinal m for defining ffness param to spurious m ng finite el 83; Flanaga inary runs i s lateral com of deforma , because th f deformat ction of loa resist again s, effects f ct provides s from uni . It should b fitting was ing materia aterial direc material or eters, such odes (Gass ements with n and Belyts n uni-axial ponents do tion in uni- ere is no te ion in which ding become st this mech rom the late the opportu -axial loadi e noted tha l orientation tion. ientation alo as shear mo er at al. 200 reduced in chko 1981) loading, the not have m axial loadin nsion in the they tend straightene anism, inte ral stiffnes nity to fix ng and then t in all axia . The x-dire ng yarn dulus, can c 0). To solve tegration (c . major effec uch effect. g. As previo yarns norm to increase d. Because ractions bet s do not pla the longitu finding la l test runs, 67 ction ause this alled ts on This usly al to their there ween y an dinal teral shear 68 ܧଵଵሺߝଵଵሻ ൌ ቐ 100 MPa ߝଵଵ ൏ 1.0 ൈ 10ିଷ 5 GPa 1.0 ൈ 10ିଷ ߝଵଵ ൏ 1.6 ൈ 10ିଷ 50 GPa 1.6 ൈ 10ିଷ ߝଵଵ This yields the following response curves under uni-axial loading. Figure 3-17: Response force of fitted model for uni-axial test data As it can be observed in Figure 3-17, fitted longitudinal stiffness agreeably captures experimental results. Subsequently, a similar procedure was taken for finding ܧ in the crushing formula of the lateral stiffness. The ensuing curve shown in Figure 3-18 has been created by the following fitted function ܧ௧௧ ൌ 2.5 ൈ 10ହ|ߝ௧௧ |ߝଵଵ 3.0 MPa where ݐݐ ൌ ሼ22,33ሽ. Note that the 3.0 MPa constant added to the main function is a very small value to avoid zero stiffness at the beginning of numerical calculation; hence the possibility of divergence because of zero terms in stiffness matrix is decreased. 0 20 40 60 80 100 120 140 160 180 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 Re ac ti on F or ce (N /y ar n) Strain Buet and Boisse (2001) Fitted Model 69 Figure 3-18: Response force of fitted model for equi-biaxial test data To check the validity of the postulated material properties (model constants) and assure they work well for other cases, beside the fitted experiments, two more runs were conducted. These cases include proportional bi-axial loadings. Test specimens are loaded under the biaxial displacement load with strain ratios (k) of 0.5 and 2.0. Figure 3-19 shows the new results, which confirm the validity of applied model and material behavior. Figure 3-19: Response force of fitted model for unbalanced bi-axial test data Similarly, another inverse method was next performed to use the results of the shear frame test. Data was taken from a benchmark work on shear frame tests for glass fabrics (Cao et al. 2008). The latter work includes results of shear tests for similar samples tested by different laboratories around the world. Although the geometrical dimensions of fabrics for this series of 0 20 40 60 80 100 120 140 160 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Re ac ti on fo rc e (N /y ar n) Strain Buet and Boisse (2001) Fitted model 0 20 40 60 80 100 120 140 160 0 0.002 0.004 0.006 Re ac ti on fo rc e (N /y ar n) Strain k = 2 Buet and Boisse (2001) Fitted model 0 10 20 30 40 50 60 70 80 0 0.001 0.002 0.003 Strain k = 0.5 Buet and Boisse (2001) Fitted model 70 tests are different from those of axial tests above, they are both using similar materials. Thus, an attempt was made to use exactly the same material properties for this case. However, it was revealed that even though previously fitted longitudinal stiffness works quite well in shear, the crushing formula does not. This can be related to the nature of the interactions between yarns in axial loading and shear loading. In axial loading, interactions at contact surfaces are in the form of pressure from upper yarn to the yarn underneath. On the other hand, in shear loading there is side-to-side contact of yarns, which happens in larger shear strains (locking). Thus, the formula found for bi-axial loads is not capable of representing the locking correctly. There are two options for resolving this situation: either using a more complicated formula that can work for both shear and bi-axial loading or trying to find a new ܧ coefficient for the shear frame test. In this project, because the focus is on loading modes separately, the second method was chosen. By fitting the new ܧ coefficient the following formula has been proposed: ܧ௧௧ ൌ 1.0 ൈ 10ସ|ߝ௧௧ |ߝଵଵ 0.3 MPa which yields the fabric shear response shown in Figure 3-20. Using this ܧ coefficient the response curve falls between experiments reported in Cao et al. (2008). Figure 3-20 shows two of these experiments reported by Hong Kong University of Science and Technology (HKUST) and the University of Massachusetts Lowell (UML). Figure 3-20: Results of response curves from fitted model vs. two experiments from (Cao et al. 2001) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 10 20 30 40 50 N or m al iz ed r ea ct io n fo rc e (N /m m ) Shear angle (degree) Experiment (HKUST) Experiment (UML) Fitted model 71 3.8 Summary This chapter began by defining a geometrical parametric model of yarn in CATIA. Additionally, yarn geometrical models were imported into Abaqus and boundary conditions for defining loads and periodicity were defined based on the loading mode (axial or shear). Next, a special material behavior of yarns in dry fabrics was incorporated into Abaqus by FORTRAN subroutines, which were called during the analysis by the Abaqus solvers. Eventually, experimental results, along with an inverse method, were utilized to find the reference material properties, which are used in the following chapters for studying the effect of uncertainty factors. The identified material properties were based on two separate yarn crushing formulas depending on the given loading mode (axial or shear). 72 Chapter 4: Uncertainty factors and their modeling This chapter introduces the uncertainties and defects that often exist in woven fabrics at meso-level and explains their influence on the response of fabric at macro-level. To facilitate their incorporation into finite element models during factorial designs, the uncertainty factors are categorized into geometrical and material groups and a range of variation for factors of each group has been identified. 4.1 Introduction The previous chapters have discussed the modeling of the mechanical behavior of an ideal plain weave unit cell. As mentioned in Section 1-4, the main objective of the current study is to identify the sensitivity of this modeling to a variety of uncertainty factors that can exist in practice. Previous models, which have been used in finite element procedures, are based on an idealization of fabric geometry and its mechanical behavior. Having an estimation of the accuracy that may be reduced by making this simplifying idealization is of high importance for designers. There needs to be a comprehensive study on the models incorporating the uncertainties, as well as the identification of the most significant factors. The goal of this chapter is to introduce fabric uncertainty factors and explain how they can be incorporated into finite element models. Let us first begin with a discussion on the importance of these factors with reference to the experimental evidence. 4.2 Uncertainty of test results Experimental tests are the most reliable and common source to identify material properties. Test instruments like universal tension machines, shear frames and biaxial fixtures are among most common setups used for identifying material properties of woven fabrics. In addition, the developed strategy in meso-level modeling of fabrics to identify their mechanical behavior in meso-level can be an alternative method to identify effective material properties without conducting time consuming and expensive experimental tests. However, by developing a numerical model, deterministic results are obtained in output, whereas in actual experiments, there are usually different (probabilistic) results every time that the test is repeated. Although the 73 experimental results may be close to each other within acceptable standard deviations, their differences are still indicative of unavoidable noise factors in the response of fabrics. This may be related to the uncertainties that are intrinsic to the nature of fabric materials and test setups. By taking a closer look at a section of a typical woven fabric, it can be observed that it has several irregularities in yarn shapes, as well as voids and defects, which can be related to effective material properties (Figure 4-1). Figure 4-1: A tomographic image of a woven fabric and its section (Badel et al. 2009) These uncertainties in the response of fabrics present themselves both in axial and shear tests by means of non-repeatable data (Milani and Nemes 2004; Lussier and Chen 2002). Another interesting case study highlighting the effect of uncertainties is the work of Cao et al. (2008), in which three types of fabrics were provided from one producer and tested in different laboratories under shear frame and trellising tests. Combining all the normalized results from samples together showed surprisingly different results. Figure 4-2 reveals these results for shear frame and bias-extension (trellising) tests for plain weave. The legend on the right denotes the abbreviation of laboratory that where the samples were tested. The tests are separated into two categories: shear frame test and the bias-extension test. 74 Figure 4-2: Results of shear tests on similar woven fabric samples by different laboratories (Cao et al. 2008) The facts mentioned above clearly indicate the need for a model that can predict or at least highlight the effect of uncertainties on the fabric mechanical behavior. For this purpose, the first step should be to identify a set of uncertainty factors, then incorporate them into the finite element model and, finally, try to find their level of importance via factorial designs. Uncertainty factors studied here fall into two separate categories. They consist of geometrical and material uncertainty factors, where the former includes 4 factors and latter includes 3 factors, which lead to a 24 geometrical and a 23 material factorial design, respectively. Firstly, the reason for this is to decrease the number of computer experiments and, second, to compare similar types of factors together. For each category, a separate sensitivity analysis is conducted to identify the most significant factors. 4.3 Geometrical uncertainty factors Among the geometrical factors, the most important are the fabric unit cell dimensions: yarn spacing (s), yarn width (w) and yarn height (h). Furthermore, by taking a closer look at the yarn shape (Figure 4-3), another geometrical defect can be observed, which is the misalignment of the yarns angle. This can be caused by poor manufacturing procedures or pre-shearing, which can be 75 present in the fabric at unloaded configurations. As a result, these four parameters are considered in the geometrical 2-level factorial analysis. Figure 4-3: Misalignment in a typical woven fabric (Milani et al. 2007) The next step after identifying geometrical uncertainty factors is finding their variation intervals. This is considered an important step for reliable results to be achieved. Accordingly, there is no better means than using the reported actual values of geometrical parameters of fabric samples. In practical terms, the range of such variations depends highly upon the manufacturing process, type of fabrics, as well as the possibility of other uncontrollable noise factors. Here, attention is focused on the plain weave used in experiments by Cao et al. (2008), where Vetrotex Saint-Gobain Co. provided the samples made of Glass/PP. Table 1 from Cao et al. (2008) reports the measured dimensions of fabrics by different laboratories. A tolerance value of ±2.5% can be noticed in the reported range of yarn spacing (s). Similarly, ±3.55% exists for reported values of yarn width (w). Unfortunately, there is no sufficient data reported for the thickness of fabric. Therefore a ± 4.0 % deviation as a high estimation from the other two factor values has been applied, considering the fact that waviness in a manufactured fabric can be more random and uncontrollable compared to its spacing and width. Moreover, it has been recommended that in single replicate factorial designs for factors where the actual range of variation is not available, a higher estimation of the range should be applied (Montgomery 2000). For yarn misalignment it has been reported by Milani et al. (2007) to be up to of ±5.0 degrees. Table 1 shows the selected factors of uncertainty, their ranges of variation and the reason for selecting them. 76 Figure 4-4: The plain woven fabric used by Cao et al. (2008) Table 4-1: Geometrical uncertainty factors Parameter S (yarn spacing) w (yarn width) (h) Yarn thickness φ (misalignment) Change range ± 2.5 % ± 3.55 % ± 4.0 % ± 5° Reason data sources in Cao et al. (2008) data sources in Cao et al. (2008) upper estimation from previous data Milani et al. (2007) Different yarn spacing, width and height are modeled in the finite element unit cell by changing the dimensions in the CAD model followed by regenerating the finite element model and solving it for new results. However, it is important to note that two separate definitions of misalignments should be considered for shear and axial loading cases. This is based on the fact that in axial loading there is no difference between positive and negative misalignment due to symmetry. On the other hand, this is not the case for shear mode, where positive and negative misalignment can respectively delay or expedite the fiber locking (it is by definition when two adjacent yarns start a side-to-side contact). Figure 4-5 shows the definition of positive misalignment for the axial loading. Because of symmetry, the basic configuration (zero misalignment) in axial loading has been considered to be the lower bound of the misalignment case. Figu loading. occur soo 4.4 Mate By r constants importan shear mo shear mo actual pa shear stif re 4-6 show Positive or n ner (positiv Figur rial uncerta ecalling the in the st t material pa dulus are th dulus is not rameters tha fness is a v Figure 4- s the defin egative mis e) or later (n (a) e 4-6 : Negat inty factor constitutiv ress-strain r rameters ca ree factors a real facto t are measu alue that s 5: Positive m ed positive aligned case egative) tha ive (a) and p s e model use elationship n be recogn that exist in r here. In o red and dire hould be va isalignment and negat s are define n the no-mi ositive (b) m d for finite (e.g. long ized. Longi the materia ther words, ctly relate t nishing to for axial loa ive misalig d based on salignment isalignment element m itudinal and tudinal stiffn l constitutiv longitudina o the respon make the y ding nment angl whether it m case. (b) s in shear te odeling, esp lateral st ess, transve e model of l and transv se of the m arn material es for the akes the loc st ecially diff iffnesses), rse stiffnes yarns. How erse stiffnes aterial, whil model acc 77 shear king erent some s and ever, s are e the urate 78 (Gasser et al. 2000). If, however, it becomes absolutely zero, a divergence of the numerical procedure can occur, and thus it must hold a small non-zero value. Shear modulus in the yarn constitutive model is only a stabilizer value. Friction between the fabric yarns can be considered as another uncertainty factor, which can be related to the material. Unfortunately, there are no comprehensive studies that show valid and reliable ranges of variation of stiffness and friction coefficient values. Therefore, an estimation of these values is used based on other observations or reasoning as follows. For the stiffness values, the work by Milani et al. (2007) has been considered, where different samples were tested under axial and shear frame tests. According to the data, a range of ± 8% has been observed for stiffness of fibres under axial tension. Although this value was seen for the fabric and not for yarns, it was used as an approximation for the variation in longitudinal stiffness of, as well as the transverse stiffness. For the friction coefficient, the same problem prevails. Peng and Cao (2002) used a friction coefficient of 0.05 for glass/PP fabrics in the finite element modeling of a similar unit cell. However, the question is, “what is the appropriate range for this factor?” Here, for the sake of the numerical procedure, noticing that applying or not applying friction behavior can affect the numerical results, the zero friction factor (no friction) has been selected as the lower bound for this factor. To keep the symmetry of the levels in 2-level factorial design, a higher bound of 0.1 for this coefficient has been selected. Table 4-2 illustrates the factors and their levels used for the study of material uncertainties. Table 4-2: Material uncertainty factors Parameter E1 (Longitudinal Young Modulus) E0 (Transverse Young Modulus) (f) Friction Coeff. Change range ± 8 % ± 8 % 0 - 0.1 Reason Milani et al (2007) Same reason as E1 0 and double value by Peng and Cao (2002) 79 4.5 Summary Ideal unit cell models are incapable of representing uncertainties in the real fabric behavior. Results of experimental measurements show different uncertainty sources in woven fabrics that can lead to non-repeatability in test data and cause uncertainty in the performance of final products. In this chapter, the uncertainty sources and flaws at meso-level were introduced in two categories of geometrical and material uncertainties. A range of variation for each factor was extracted directly from the reported experimental measurements or indirectly from physical interpretation. Eventually, a table of uncertainty factors for each category was established for use in factorial designs in the following chapter. 80 Chapter 5: Sensitivity analysis on uncertainty factors After identifying the uncertainty factors and their variation range, the next step is to recognize the main factors and their interactions that play an important role in affecting the fabric response. Therefore, this chapter covers factorial designs with the help of computer experiments that are conducted on each series of geometrical and material uncertainty sets. It begins by defining a base unit cell, which uses the nominal values of uncertainty factors. Then parametric studies are conducted on this unit cell by varying parameters according to two-level factorial designs. Finally, results of factorial runs are fed into the Design Expert software to quantify the sensitivity of the unit cell response to the study factors and their interactions. 5.1 Base model and its validation Recalling Figure 2-4 (a), attention is focused towards a plain-weave fabric. Three nominal dimensions for this unit cell are s = 5.13 (mm), w = 4.44 (mm) and h = 0.5 (mm). These dimensions are selected to be close to the fabric tested in Cao et al. (2008). Thus, a comparison can be made between the variations of the results from the current design of experiment to those observed in the actual experiments. Material properties are as extracted in Chapter 3. The next step is defining each basic type of loading and a response for the unit cell to be used in the factorial design. For the uni-axial test, the strain values up to 5.848 ൈ 10ିଷ are selected, which needs a 0.3 mm displacement at the moving end of the boundary condition. This is close to the value of the strain measured in axial tests by Buet-Gautier and Boisse (2001). The same amount of displacement is selected for the basic unit cell under bi-axial loading. Figure 5-1 (a and b) shows the response of this nominal unit cell to the predefined axial loading, along with a case with finer mesh, and two cases with more material orientation sections and fewer sections along the yarn. As it can be observed in these figures, the original mesh size and number of sections are adequate for the basic model. 81 (a) (b) Figure 5-1: Sensitivity of basic model to mesh size and number of sections under axial loading The same unit cell model has been used for shear loading. The only differences are the type of boundary conditions that have been applied for simulating the shear frame test and the material coefficients, which are used for transverse stiffness. The applied loading parameter for the shear frame test is the stretch ratio (ߣ), which is defined as the ratio of diagonal length 0 5 10 15 20 25 30 0 0.001 0.002 0.003 0.004 0.005 0.006 Lo ad / U ni t c el l ( N ) Strain uni‐axial Basic 3 Sections Finer mesh 20 Sections 0 20 40 60 80 100 120 140 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 Lo ad / U ni t c el l ( N ) Strain bi‐axial Basic 20 Sec Finer Mesh 82 between deformed frame and un-deformed frame ߣ ൌ ܮଶ/ܮଵ (Figure 5-2). Stretch ratios up to 1.3 are considered for the factorial design. To check the mesh size and the number of sections for the basic shear model, the same unit cell has been run with a fewer number of sections and finer mesh. Figure 5-3 shows the results of this analysis, which proves the reliability of the basic model (mesh and section sensitivity). Figure 5-2: Diagonal length of the shear frame used for calculating stretch ratio (before and after deformation) Figure 5-3: Sensitivity of basic model to mesh size and number of sections under shear loading 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 N or m al iz ed fo rc e (N /m m ) Stretch ratio Basic 4 Section Finer Mesh 83 5.2 Factorial design: Fundamentals There are many circumstances where analysts deal with studying the effect of design factors. Generally, this task can be done effectively via factorial designs. In a factorial design, all possible combinations of factors at different levels are investigated. Then, results are analyzed to find the effect of each factor or interactions between them (Montgomery 2000). The effect of a factor is the change in response produced by a change in the level of the factor. However, sometimes the difference in response between the levels of one factor is not the same at all levels of the other factors. This occurs when there is an interaction between factors. The factorial design analysis provides knowledge of the average response changes that are expected by varying a factor or a combination of factors if there is an interaction. It can also highlight the importance of each factor with respect to others by calculating their percentage contribution on the output response. In order to keep the number of runs to a minimum and to have an optimum computational cost, it is preferred to define two levels for each factor. This is called a two-level factorial design or 2k factorial design, where k is the number of factors. If factors are affecting the response monotonically and their ranges of variation are not large enough to make a highly nonlinear response, it can be expected that this method will be highly effective. Two levels of high (positive) and low (negative) are defined for each factor; therefore, there are 2k runs in total. In physical experiments, the repeated tests in each combination of factorial design gives an internal estimation of error via analysis of variances (ANOVA), which can then be used to determine whether a factor has a considerable influence on the results (its contribution toward the response is higher than errors) or it is ineffective (Montgomery 2000). However, in computer experiments there is only one repeat for each run. Therefore, there should be other alternatives for assessing the effect of factors. There are different sensitivity analysis methods suggested for single replicate factorial designs. Daniel (1959) suggested using normal probability plots of the estimated effects. He argued, effects that are negligible tend to fall along the straight line on this plot. On the other hand, significant effects do not locate close to this line and are farther from that if they show more contribution towards the response. Similarly, half-normal plots can be used for single replicate designs. This is a plot of the absolute value of the effect estimates against their cumulative normal probabilities. Like the standard normal probability plots, in half-normal 84 probability plots the negligible factors are located along a straight line, whereas significant factors are away from this line relative to their contribution percentage. Although the most common method for the analysis of single replicate design is normal and half-normal probability plots, these methods can become subjective. This can lead to some difficulties in practice. In a comprehensive survey on other methods for single replicate design, Hamada and Balakrishnan (1998) found that the method proposed by Lenth (1989) could advantageously detect significant effects. This method is also easy to implement and, as a result, it appears in some software packages including “Design Expert” for analyzing data from single replicate factorials. The Lenth’s method is briefly described as follows. Suppose there are m contrasts (effects from factors and interactions) in a factorial design and are referred by c1, c2 , ... , cm ; hence for a 2k full factorial design, ݉ ൌ 2 െ 1. Lenth’s method starts with an estimation of the variance of the smallest contrast (in absolute value). Let ܲܵܧ ൌ 1.5 ൈ median൫ห ܿห: ห ܿห ൏ 2.5 ݏ൯ (5-1) ݏ ൌ 1.5 ൈ median൫ห ܿห൯ (5-2) Here PSE stands for “Pseudo Standard Error” and Lenth showed that it can be used as a reasonable estimation of contrast variance when there are not many active (significant) effects (Montgomery 2000). This parameter can be used to judge the significance of contrasts. To assess the significance of each contrast it is compared to the Margin of Error (ME), which is defined as: ܯܧ ൌ ݐఈ/ଶ,ௗ ܲܵܧ (5-3) where ݐఈ/ଶ,ௗ is the t-distribution with a significance level of ߙ and ݀ ൌ ݉/3 is the degree of freedom. Similarly, factors can be compared using another parameter suggested by Lenth called Simultaneous Margin of Error (SME): ܵܯܧ ൌ ݐఊ,ௗ ܲܵܧ (5-4) Here percentage point for the t-distribution is ߛ ൌ 1 െ ൫1 √1 െ ߙ ൯/2 (5-5) In the results to follow, Lenth’s parameters are reported from the “Design Expert” software by assuming ߙ ൌ 0.05. 85 5.3 Planning factorial design For the sake of simplicity, a letter is used to refer to each study factor (see Tables 5-1 and 5- 7). Moreover, for each run in the design of experiment, if that factor is in its higher level a positive sign is added after its letter, whereas if it is in its lower level a negative sign is used. Each run is identified by a code, which consists of symbols of factors in that experiment, followed by a sign that shows either that factor is in its higher level or lower level. For example, if A+B+C- is used as a code for an experiment, it means that A and B factors are in their higher levels and C is in its lower level. 5.3.1 Choosing response function It is highly important to ensure that an appropriate response and loading parameter is selected for a factorial design. Regarding the loading parameter, as it was previously mentioned, macroscopic strain on the unit cell sample in axial mode and stretch ratio in the shear mode are applied. It is evident that regardless of the fabric type and the uncertainties that exist in the test samples, strain and stretch ratio are two independent factors that are most often used in displacement control experiments. Normalized force can be selected as a good option for the response function. The detail of this normalization depends on the loading mode. For the axial loading, reaction force at the moving boundary per unit length of that boundary has been used as a normalized response parameter. For shear mode, however, the situation becomes more complicated. There are different formulas for the normalization of shear tests (Peng et al. 2004). Nevertheless, the best method for such normalization is using the amount of absorbed energy per unit area of fabric to find a normalized force (Peng et al. 2004). If it is assumed that there is the same configuration of deformation in two given samples and the external work in the first one is w1 and in the second one w2, then it is reasonable to say ௪భ భ ൌ ௪మ మ (5-6) where Ai is the deforming area of sample 1 or 2, respectively. If the displacement at the moving head of the frame is shown by d, reaction forces by P and length of fabrics with L, then it can be written as: 86 భ ௗభ మ ௗమ ൌ భ మ ൌ భ మ మ మ (5-7) However, the displacement in the moving head of frame can be related to the shear angle and side length of frame (l) via (Peng et al. 2004) ݀ ൌ ݈ ቂ2 cos ቀ π ସ െ ఏ ଶ ቁ െ √2ቃ (5-8) It can be easily checked that if ߣଵ ൌ ߣଶ then ߠଵ ൌ ߠଶ . After substituting (5-8) in (5-7) and simplifying there are భ మ ൌ భ మ/భ మ మ/మ (5-9) which assumes using the normalized force in the form of ܲ ൌ ܲ ݈/ܮ ଶ (5-10) For the misaligned unit cells, it can be shown that by defining an imaginary frame instead of the misaligned frame and using its length in equation (5-10) there will be an appropriate normalized force (Figure 5-4). Figure 5-4: Imaginary frame defined for shear frame test under misalignment, (a) positive misalignment and (b) negative misalignment Two indicators are selected for assessing the design of experiment responses. The first one is the final value of the normalized force at the end of deformation, which is shown by I∞. This indicator can give estimation of the reaction force that we have at the final step. However, it lacks any information about previous states and the history of deformation. Therefore, a second indicator I1 is defined as the area under response curve, which represents the normalized force 87 versus loading parameter (Figure 5-5). It is worth noting that I1 can be related to the energy absorbed by the material. Figure 5-5: Schematic of two indicators used for the assessment of case studies 5.4 Geometrical sensitivity analysis For geometrical factorial design, there are yarn spacing (s), yarn width (w), yarn height (h) and misalignment to be considered. Here, the letters A-D are used to refer to these factors as shown in Table 5-1. Table 5-1: Symbols for factors used in geometrical sensitivity analysis Factor Yarn spacing (s) Yarn width (w) Yarn height (h) Misalignment Symbol A B C D As explained in Section 4-3, because of differences that have been noticed on positive and negative misalignment definitions (mainly its symmetry in axial loading), two separate geometrical sensitivity designs are applied. In axial loading, misalignment is defined in the range of [0°, +5°], where 0° is the low and +5° is the high misalignment level, respectively. In shear mode, the misalignment range is [-5°, +5°]. Ranges for other factors are defined based on Table 4-1. Therefore, a table of runs for the full factorial design can be established as follows (Table 5- 2). 0 1 2 3 4 5 6 7 8 9 0 0.2 0.4 0.6 0.8 1 1.2 N or m al iz ed fo rc e Strain (axial mode) Stretch ratio (shear mode) I∞ I1 88 Table 5-2: The value of each factor in geometrical sensitivity analysis Run Number A B C D (Axial/Shear) Code of run 1 5 4.07 0.48 0/-5 A-B-C-D- 2 5.26 4.07 0.48 0/-5 A+B-C-D- 3 5 4.37 0.48 0-5 A-B+C-D- 4 5.26 4.37 0.48 0/-5 A+B+C-D- 5 5 4.07 0.52 0/-5 A-B-C+D- 6 5.26 4.07 0.52 0/-5 A+B-C+D- 7 5 4.37 0.52 0/-5 A-B+C+D- 8 5.26 4.37 0.52 0/-5 A+B+C+D- 9 5 4.07 0.48 5/5 A-B-C-D+ 10 5.26 4.07 0.48 5/5 A+B-C-D+ 11 5 4.37 0.48 5/5 A-B+C-D+ 12 5.26 4.37 0.48 5/5 A+B+C-D+ 13 5 4.07 0.52 5/5 A-B-C+D+ 14 5.26 4.07 0.52 5/5 A+B-C+D+ 15 5 4.37 0.52 5/5 A-B+C+D+ 16 5.26 4.37 0.52 5/5 A+B+C+D+ 5.4.1 Axial loading The results of all runs for uni-axial and bi-axial loadings are illustrated in Figure 5-6 in terms of normalized reaction force versus strain. Table 5-3 also shows the output of runs in terms of difference between indicators of each run and the basic unit cell (I∞,run- I∞,nominal and I1,run- I1,nominal). These results were inserted into Design Expert to analyze the effect of each factor and interaction between factors. Considering that four different outputs are reported for this design in Table 5-3, four separate sensitivity analyses are available. 89 Figure 5-6: Result of factorial design runs for axial mode (the corresponding values of indicators for each run are given in Table 5-3) 0 5 10 15 20 25 30 0 0.001 0.002 0.003 0.004 0.005 0.006 N or m al iz ed R ea ct io n Fo rc e (N /m m ) Strain bi‐axial mode A+B+C+D+ A‐B+C+D+ A+B‐C+D+ A‐B‐C+D+ A+B+C‐D+ A‐B+C‐D+ A+B‐C‐D+ A‐B‐C‐D+ A+B+C+D‐ A‐B+C+D‐ A+B‐C+D‐ A‐B‐C+D‐ A+B+C‐D‐ A‐B+C‐D‐ A+B‐C‐D‐ A‐B‐C‐D‐ Basic 0 1 2 3 4 5 6 7 8 9 10 0 0.001 0.002 0.003 0.004 0.005 0.006 N or m al iz ed R ea ct io n fo rc e (N /m m ) Strain uni‐axial mode A+B+C+D+ A‐B+C+D+ A+B‐C+D+ A‐B‐C+D+ A+B+C‐D+ A‐B+C‐D+ A+B‐C‐D+ A‐B‐C‐D+ A+B+C+D‐ A‐B+C+D‐ A+B‐C+D‐ A‐B‐C+D‐ A+B+C‐D‐ A‐B+C‐D‐ A+B‐C‐D‐ A‐B‐C‐D‐ Basic 90 Table 5-3: The values of response indicators in each run relative to those of the base unit cell (axial loading) Run Number Code of run I∞ uni-axial I1 uni-axial I∞ bi-axial I1 bi-axial 1 A-B-C-D- -0.1916 -0.0003 -0.57192 -0.00114 2 A+B-C-D- 0.51372 0.000656 0.490606 0.001171 3 A-B+C-D- 1.01068 0.001979 0.275882 0.001729 4 A+B+C-D- 1.69857 0.00284 2.05829 0.004923 5 A-B-C+D- -1.54142 -0.00241 -3.73812 -0.01003 6 A+B-C+D- -1.31748 -0.00194 -1.45825 -0.00471 7 A-B+C+D- -0.61916 -0.00059 -2.42092 -0.00461 8 A+B+C+D- -0.07373 -7.5E-05 0.264751 0.000391 9 A-B-C-D+ -0.99036 -0.00143 -5.15992 -0.00926 10 A+B-C-D+ -0.97457 -0.00154 -3.99685 -0.00626 11 A-B+C-D+ -0.26494 -7.5E-05 -4.20052 -0.00726 12 A+B+C-D+ -0.84478 -0.001 -4.0419 -0.00732 13 A-B-C+D+ -2.06018 -0.00316 -7.43372 -0.01564 14 A+B-C+D+ -2.16757 -0.0034 -5.28373 -0.00956 15 A-B+C+D+ -1.4076 -0.00194 -5.60492 -0.01072 16 A+B+C+D+ -1.69525 -0.00235 -4.45844 -0.00755 Figure 5-7 shows the half-normal probability plots for deviations of two indicators from the basic case in two modes of uni-axial and bi-axial loadings. Points on these plots are shown by two different colors in which orange represents positive and blue represents negative effects. Table 5-4 also shows the numerical illustration of these graphs in terms of percentage contribution and Lenth’s ME and SME with a 0.05 significance threshold. 91 Figure 5-7: Half normal probability plots for geometrical factorial design (axial loading) Table 5-4: Results of the 2-level material factorial design for shear tests; A (yarn spacing), B (yarn width), C (yarn height) and D (misalignment) uni-Axial bi-Axial I∞ I1 I∞ I1 Effect % Cont. Effect % Cont. Effect % Cont. Effect % Cont. A 0.150435 0.5 0.00014 0.2 1.553579 8.8 0.003501 10.8 B 0.816656 15.0 0.001539 21.0 1.128016 4.6 0.003125 8.6 C -1.35489 41.3 -0.00212 40.0 -1.87338 12.8 -0.00488 21.0 D -1.2356 34.4 -0.00188 31.4 -4.38504 70.0 -0.00766 51.9 AB -0.05898 0.1 -0.00013 0.2 -0.11028 0.0 -0.00067 0.4 AC -0.05685 0.1 -5.5E-05 0.0 0.511924 1.0 0.001391 1.7 AD -0.39021 3.4 -0.00056 2.8 -0.39904 0.6 -0.00046 0.2 BC 0.006072 0.0 -4.9E-05 0.0 0.295557 0.3 0.001234 1.3 BD -0.32163 2.3 -0.0005 2.2 -0.23591 0.2 -0.00116 1.2 CD 0.290899 1.9 0.000422 1.6 0.527972 1.0 0.001533 2.1 ABC 0.094285 0.2 9.5E-05 0.1 -0.03914 0.0 -0.00013 0.0 ABD -0.13499 0.4 -0.00012 0.1 -0.39171 0.6 -0.00082 0.6 ACD 0.099106 0.2 0.000152 0.2 -0.01823 0.0 0.000185 0.0 BCD 0.06135 0.1 0.000144 0.2 0.139378 0.1 0.000261 0.1 ABCD 0.009557 0.0 6.32E-05 0.0 0.039378 0.0 0.000168 0.0 Lenth’s ME 0.372846 0.000506 1.024626 0.003148 Lenth’s SME 0.756933 0.001028 2.08014 0.00639 Half Normal Plot 0.0000 0.0005 0.0011 0.0016 0.0021 0 10 20 30 50 70 80 90 95 99 A B C D AB AC AD BC BD CD ABCABD ACD BCD ABCD Half-Normal Plot 0.00 1.10 2.19 3.29 4.39 0 10 20 30 50 70 80 90 95 99 A B C D AB AC AD BC BD CD ABC ABD ACD BCD ABCD Half Normal Plot 0.000 0.002 0.004 0.006 0.008 0 10 20 30 50 70 80 90 95 99 A B C D AB AC AD BC BD CD ABC ABD ACDBCDABCD Half Normal Plot 0.00 0.34 0.68 1.02 1.35 0 10 20 30 50 70 80 90 95 99 A B C D ABAC AD BC BD CD ABC ABD ACD BCD ABCD I∞ I1 un i-a xi al m od e bi -a xi al m od e |Standardized Effect| |Standardized Effect| H al f n or m al % p ro ba bi lit y H al f n or m al % p ro ba bi lit y 92 From the illustrated plots and data provided in Table 5-4, the number of important factors and their level of their importance can be extracted. From Table 5-4 it can be concluded that for uni-axial loading B, C and D are important factors, because they exceed Lenth’s SME for both I∞ and I1. In addition, AD interaction may also be considered marginally important, because it exceeds Lenth’s ME, but it is still less than Lenth’s SME. These results have been graphically shown in Figure 5-7. Similarly, for the bi-axial loading mode, D is an absolutely significant factor as it exceeds Lenth’s ME and Lenth’s SME where A, B and C are factors which may be considered significant or insignificant based on the selected criteria, because they only exceed Lenth’s ME, not Lenth’s SME; especially for B, which is very close to Lenth’s ME for I∞ and barely fails to exceed it for I1. To give another visual representation of the main factors, Figure 5-8 shows the effect from each main factor by plotting the difference between I∞ for the basic case and the factorial design runs. The point that was previously noted can also be observed here. Among the interaction factors, only AD in uni-axial tests shows a noticeable influence on the results. This interaction effect is explicitly illustrated in Figure 5-9. For brevity the rest of interactions are not plotted here, however, they all show a parallel trend, indicating their very low importance compared to the other factors. I ∞ -I ∞ , n om in al I ∞ -I ∞ , n om in al uni‐axial bi‐axial Figure 5-8: Effects of main factors in their normalized range -1 0 1 -1.5 -1 -0.5 0 A: Yarn Spacing -1 0 1 -1.5 -1 -0.5 0 B: Yarn Width -1 0 1 -1.5 -1 -0.5 0 C: Yarn Height -1 0 1 -1.5 -1 -0.5 0 D: Misalignment -1 -0.5 0 0.5 1 -6 -4 -2 0 A: Yarn Spacing -1 -0.5 0 0.5 1 -6 -4 -2 0 B: Yarn Width -1 -0.5 0 0.5 1 -6 -4 -2 0 C: Yarn Height -1 -0.5 0 0.5 1 -6 -4 -2 0 D: Misalignment 93 Figure 5-9: Effect of AD interaction in normalized range of A (uni-axial mode) 5.4.2 Shear loading A similar set of factorial design to the one conducted in axial loading is conducted on the shear mode. As previously mentioned, the only difference with axial loading is the assigned range of variation of misalignment, which is [-5°, +5°]. Thus, for the case of D- we refer to -5° misalignment. Figure 5-10 shows the results of the runs. Moreover, the difference between I∞ and I1 for the basic case and factorial design runs is summarized in Table 5-4. A table of effects and half-normal probability plots are illustrated in Table 5-5 and Figure 5-11. -1 -0.5 0 0.5 1 -1.5 -1 -0.5 0 0.5 A: Yarn Spacing I ∞ - I ∞ , n om in al 94 Figure 5-10: Result of factorial design runs for shear mode (the corresponding values of indicators for each run are given in Table 5-5) Table 5-5: The difference between values of different indicators in each run using the basic unit cell (shear loading) Run Number Code of run I∞ I1 1 A-B-C-D- -0.10632 -0.0074 2 A+B-C-D- -0.11959 -0.00855 3 A-B+C-D- -0.08882 -0.00561 4 A+B+C-D- -0.10474 -0.00716 5 A-B-C+D- -0.10455 -0.0073 6 A+B-C+D- -0.11104 -0.008 7 A-B+C+D- -0.0736 -0.00452 8 A+B+C+D- -0.0932 -0.00635 9 A-B-C-D+ 0.68708 0.024762 10 A+B-C-D+ 0.573014 0.01848 11 A-B+C-D+ 1.32065 0.044763 12 A+B+C-D+ 0.87831 0.02975 13 A-B-C+D+ 0.884028 0.032952 14 A+B-C+D+ 0.64965 0.02325 15 A-B+C+D+ 1.2168 0.044561 16 A+B+C+D+ 1.04441 0.035066 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.05 0.1 0.15 0.2 0.25 0.3 N or m al iz ed R ea ct io n Fo rc e (N /m m ) Stretch Ratio Shear mode A+B+C+D+ A‐B+C+D+ A+B‐C+D+ A‐B‐C+D+ A+B+C‐D+ A‐B+C‐D+ A+B‐C‐D+ A‐B‐C‐D+ A+B+C+D‐ A‐B+C+D‐ A+B‐C+D‐ A‐B‐C+D‐ A+B+C‐D‐ A‐B+C‐D‐ A+B‐C‐D‐ A‐B‐C‐D‐ Basic1.0 1 1. 1.15 1.2 1.25 1.3 95 Figure 5-11: Half normal probability plots for geometrical factorial design (shear loading) Table 5-6: Results of the 2-level geometrical factorial design for shear tests; A (yarn spacing), B (yarn width), C (yarn height) and D (misalignment). I∞ I1 Effect % Cont. Effect % Cont. A -0.12731 1.4 -0.00572 2.0 B 0.218443 4.2 0.007788 3.7 C 0.046614 0.2 0.002577 0.4 D 1.006976 88.7 0.038557 90.0 AB -0.03526 0.1 -0.00126 0.1 AC 0.019094 0.0 0.000284 0.0 AD -0.11349 1.1 -0.00441 1.2 BC -0.02436 0.1 -0.00083 0.0 BD 0.198156 3.4 0.005886 2.1 CD 0.037345 0.1 0.001942 0.2 ABC 0.047477 0.2 0.001025 0.1 ABD -0.03132 0.1 -0.00087 0.0 ACD 0.018316 0.0 0.000241 0.0 BCD -0.02847 0.1 -0.00114 0.1 ABCD 0.050089 0.2 0.00121 0.1 Lenth’s ME 0.139969 0.004383 Lenth’s SME 0.284158 0.008898 As it can be noted from the plots in Figure 5-11 and Table 5-6, D is the most important factor when considering both criteria (Lenth’s ME and SME). It has such a large influence on the output results that it can even be considered as the only important factor for shear test results. Its percentage contribution of ~90% clearly illustrates this fact. After that, the main factor B and BD interactions can be considered of second-level importance, as they exceed Lenth’s ME, but not Lenth’s SME. Moreover, the main factor A and AD interactions show a slight contribution by Half Normal Plot 0.00 0.01 0.02 0.03 0.04 0 10 20 30 50 70 80 90 95 99 A B C D AB AC AD BC BD CD ABCABD ACD BCD ABCD Half Normal Plot 0.00 0.25 0.50 0.76 1.01 0 10 20 30 50 70 80 90 95 99 A B C D AB AC AD BC BD CD ABC ABD ACD BCD ABCD I1I∞ |Standardized Effect| |Standardized Effect| H al f n or m al % p ro ba bi lit y 96 barely passing the criteria for Lenth’s ME when considering the I1 indicator. To visually compare the effect of important factors (B and D) they are plotted in Figure 5-12 for the I∞ indicator. Figure 5-13 also shows the most important interaction (BD) in this loading mode. Figure 5-12: Effects of important main factors in their normalized range Figure 5-13: Effect of BD interaction in the normalized range of B -1 -0.5 0 0.5 1 -0.2 0 0.2 0.4 0.6 0.8 I ∞ - I ∞ , n om in al B: Yarn Width -1 -0.5 0 0.5 1 -0.2 0 0.2 0.4 0.6 0.8 D: Misalignment -1 -0.5 0 0.5 1 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 B: Yarn Width I ∞ - I ∞ , n om in al D- D+ 97 5.5 Material sensitivity analysis Sensitivity analysis on the material properties includes three factors of longitudinal stuffiness, transverse stiffness (crushing effect) and friction coefficient. The defined intervals for these factors were discussed in Section 4-4. Here again, the names of the factors are replaced by symbols to make the factorial design procedure more convenient in subsequent discussions. The symbol for each factor in the material sensitivity analysis is shown in Table 5-7. Table 5-8 also shows the code of each run in this factorial design and its meaning in terms of values of material properties. Table 5-7: Symbols for factors used in material sensitivity factorial analysis Factor Longitudinal stiffness fraction (E11) Transverse stiffness fraction (Ett) Friction coefficient (f) Symbol A B C Table 5-8: The value of each factor for material sensitivity analysis Run Number A B C Code of Run 1 0.92 E11 0.92 Ett 0.0 A-B-C- 2 1.08 E11 0.92 Ett 0.0 A+B-C- 3 0.92 E11 1.08 Ett 0.0 A-B+C- 4 1.08 E11 1.08 Ett 0.0 A+B+C- 5 0.92 E11 0.92 Ett 0.1 A-B-C+ 6 1.08 E11 0.92 Ett 0.1 A+B-C+ 7 0.92 E11 1.08 Ett 0.1 A-B+C+ 8 1.08 E11 1.08 Ett 0.1 A+B+C+ As it can be observed from Table 5-8, a two-level full factorial design with three factors, which has eight test samples (unit cells with different material properties) is dealt with. In addition, there are three tests for each case, which includes uni-axial, bi-axial and shear. Moreover, two indicators for each test are considered. Therefore, for each run of the factorial design, six output results should be calculated. These results are summarized in Figure 5-14 and Table 5-9. 98 Figure 5-14: Result of factorial design runs (the corresponding values of indicators for each run are given in Table 5-9) 0 1 2 3 4 5 6 7 8 9 0 0.001 0.002 0.003 0.004 0.005 0.006 N or m al iz ed R ac ti on F or ce (N /m m ) Strain uni‐axial mode Basic A+B+C+ A‐B+C+ A+B‐C+ A‐B‐C+ A+B+C‐ A‐B+C‐ A+B‐C‐ A‐B‐C‐ 0 5 10 15 20 25 30 35 0 0.001 0.002 0.003 0.004 0.005 0.006 N or m al iz ed R ea ct io n Fo rc e (N /m m ) Strain bi‐axial mode Basic A+B+C+ A‐B+C+ A+B‐C+ A‐B‐C+ A+B+C‐ A‐B+C‐ A+B‐C‐ A‐B‐C‐ 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.05 0.1 0.15 0.2 0.25 0.3 N or m al iz ed R ea ct io n Fo rc e (n /m m ) Stretch ratio shear mode Basic A+B+C+ A‐B+C+ A+B‐C+ A‐B‐C+ A+B+C‐ A‐B+C‐ A+B‐C‐ A‐B‐C‐ 1.0 1.05 1.1 1.15 1.2 1.25 1.3 99 Table 5-9: The values of response indicators in each run relative to those of the base unit cell (shear loading) Run Number Code of run I∞ uniAxial I1 uniAxial I∞ biAxial I1 biAxial I∞ shear I1 shear 1 A-B-C- -0.43289 -0.0007 -6.7386 -0.0146 -0.01848 -0.00138 2 A+B-C- -0.10957 -0.00024 -4.48031 -0.01116 -0.01856 -0.0014 3 A-B+C- -0.81228 -0.0013 -6.08207 -0.01249 -0.00582 -0.00054 4 A+B+C- -0.46049 -0.00074 -3.64834 -0.00852 -0.00665 -0.00056 5 A-B-C+ -0.1276 -0.00014 0.540936 -0.00047 -0.00157 0.000253 6 A+B-C+ 0.352651 0.000647 3.44522 0.003594 -0.00311 0.000211 7 A-B+C+ -0.12103 0.000182 1.35965 0.001694 0.016895 0.001408 8 A+B+C+ 0.828635 0.001643 4.81618 0.00747 0.015386 0.001364 Figure 5-15 shows the half-normal probability plots for the material factorial design. Table 5-10 demonstrates the final output of post analysis for finding the effect of factors and their contribution percentages. 100 Figure 5-15: Half-Normal probability plots of material factorial design 0.00 1.94 3.89 5.83 7.78 0 10 20 30 50 70 80 90 95 99 A B C AB AC BC ABC Half Normal Plot 0.0000 0.0004 0.0009 0.0013 0.0018 0 10 20 30 50 70 80 90 95 99 A B C AB AC BC ABC Half Normal Plot 0.00 0.17 0.34 0.52 0.69 0 10 20 30 50 70 80 90 95 99 A B C AB AC BC ABC Half Normal Plot 0.0000 0.0003 0.0007 0.0010 0.0013 0 10 20 30 50 70 80 90 95 99 A B C AB AC BC ABC Half Normal Plot 0.000 0.004 0.007 0.011 0.015 0 10 20 30 50 70 80 90 95 99 A B C AB AC BC ABC Half Normal Plot 0.000 0.005 0.010 0.014 0.019 0 10 20 30 50 70 80 90 95 99 A B C AB AC BC ABC I∞ I1 un i-a xi al m od e bi -a xi al m od e sh ea r m od e |Standardized Effect| |Standardized Effect| H al f n or m al % p ro ba bi lit y H al f n or m al % p ro ba bi lit y H al f n or m al % p ro ba bi lit y 101 Table 5-10: Results of the 2-level factorial design for geometrical shear tests; A (longitudinal stiffness), B (transverse stiffness) and C (friction coefficient) uni-axial Bi-axial Shear I∞ I1 I∞ I1 I∞ I1 Effect % Cont. Effect % Cont. Effect % Cont. Effect % Cont. Effect % Cont. Effect % Cont. A 0.5263 30.5 .000818 22.7 2.763 11.0 0.00431 7.6 -.00099 0.2 -.00003 0.0 B -0.062 0.4 0.0001 0.1 0.920 1.2 0.00270 3.0 .01538 38.2 0.00100 23.8 C 0.6870 52.0 0.0013 59.7 7.778 87.4 0.01476 89.1 .01928 60.0 0.00178 75.6 AB 0.1245 1.7 0.0002 1.3 0.182 0.0 0.00056 0.1 -.00018 0.0 0.00000 0.0 AC 0.1887 3.9 0.0003 3.1 0.417 0.3 0.00061 0.2 -.00054 0.0 -.00001 0.0 BC 0.3032 10.1 0.0006 12.4 0.175 0.0 0.00032 0.0 .00310 1.6 0.00016 0.6 ABC 0.1102 1.3 0.0001 0.7 0.094 0.0 0.00030 0.0 .00019 0.0 0.00000 0.0 ME 1.0655 0.0014 1.027 0.00249 .00303 0.00004 SME 2.5499 0.0034 2.458 0.00596 0.00725 0.00009 Analyzing the results of Table 5-10 and Figure 5-15 clearly indicates the contribution of each of the material uncertainty factors to the response of the fabric. For uni-axial loading, there is a high percentage of contribution from A, C and BC. However, comparing their effects with Lenth’s ME and SME shows that they do not exceed any of these criteria. This can also be verified in the half-normal probability plots of the uni-axial mode by noting that all the points are almost in a straight line. Thus, it can be concluded that in the uni-axial loading, although there is a variation of response, none of the material factors and interactions are statistically important enough to substantially change the output. Studying Figure 5-14 can provide an example of how this is possible. All the graphs are concentrated around the basic line and they show less variation compared to bi-axial mode. Nonetheless, it should be noted here that C (friction), A (longitudinal stiffness) and BC interaction (transverse stiffness and friction) are factors that could play a more important role if a higher value of significance threshold α was applied. For bi-axial loading, the percentage contributions of A and C show that they can be highlighted as the main important factors. Half-normal probability plots also confirm this fact. However, from SME and ME parameters it is observed that C is a very important factor, while A can be rejected using a more restrict criterion. Another important fact to be noted here is that effects from interactions are at such a low level that they are not comparable to the main effects. In shear loading, high percentage contributions from B and C and a very slight contribution from their interaction can be noticed from Table 5-10. Moreover, interpretations from half- normal plots show that these two factors are dominant. Furthermore, their effects exceed Lenth’s 102 ME and SME by a considerable margin, which insists on their high level of importance and influence on the response of fabrics. Figure 5-16 graphically illustrates the effects of the main factors in the performed material factorial design. As it can be observed under all loading modes, the effect of friction is highly considerable. In addition, the effect of longitudinal stiffness in the axial loading modes and the transverse stiffness in the shear mode are other interesting points to be observed. Figure 5-16: Effect of main factors in material factorial design -1 0 1 -1.5 -1 -0.5 0 0.5 un i-a xi al m od e I ∞ - I ∞ , n om in al -1 0 1 -1.5 -1 -0.5 0 0.5 -1 0 1 -1.5 -1 -0.5 0 0.5 -1 0 1 -6 -4 -2 0 2 4 bi -a xi al m od e I ∞ - I ∞ , n om in al -1 0 1 -6 -4 -2 0 2 4 -1 0 1 -6 -4 -2 0 2 4 -1 0 1 -0.02 0 0.02 A: E11 sh ea r m od e I ∞ - I ∞ , n om in al -1 0 1 -0.02 0 0.02 B: E22, E33 -1 0 1 -0.02 0 0.02 C: friction 103 5.6 Non-repeatability of test results Figure 4-2 illustrates the non-repeatability of experimental results in the actual shear tests. It can provide a good example to compare the results from the factorial design with the test data that already reflect uncertainty effects. Figure 5-17 shows response curves from all the geometrical and material factorial designs (grey lines), the nominal (basic) unit cell (red line), along with two selected experimental results from Cao et al., which indicates the upper and lower limits of measured forces (green lines). It should be noted that in real samples, all types of uncertainty factors are randomly distributed in the fabric, whereas in the modeled unit cell they are distributed uniformly. For instance, the assumption of having +5° misalignment everywhere in a fabric is improbable. Although there might be some regions with this amount of imperfection, their overall effect becomes less by contribution of other misaligned regions, which may even act in opposite ways. The fibre misalignment had been previously estimated as an important factor using macro-level models (Milani et al. 2010). Indeed, in the present meso- level study, according to Figure 5-17, all the runs which belong to +5° shear misalignment are higher than the nominal case (with zero misalignment), whereas those with -5° misalignment are lower than the nominal case. It can be directly related to the adjacent yarns locking in shear deformation. An interesting observation is the configuration of yarns at the end of the loading, which is completely different for positive or negative misalignment. These configurations are illustrated in Figure 5-17. As it can be observed for positive misalignment, adjacent yarns are locked into each other, while for the negative case, there is still a hollow space (trellis) between yarns and they can rotate on each other without a side contact. Thus, higher reaction forces are expected in the former case due to severe contact compared to the latter case. By focusing at lower stretch ratios, there is a better agreement between the factorial runs, because the locking has not yet occurred. However, after the occurrence of locking, the deformation mechanism changes and higher reaction forces come to action. 104 Figure 5-17: Non-repeatability of test data and comparing it to the factorial design results. Red line: nominal unit cell; Grey lines: factorial designs; Green lines: experiment 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 0.1 0.2 0.3 0 0.05 0.1 0.15 0.2 0.25 0 0.05 0.1 0.15 0.2 0.25 Stretch ratio N or m al iz ed R ea ct io n Fo rc e (N /m m ) 1.0 1.1 1.2 1.3 1.0 1 1.1 1.15 1. 1 105 5.7 Some more discussions on the results 5.7.1 Important effect of misalignment in axial loading In another work by Komeili and Milani (2009), an analysis on the effect of a set of uncertainty factors in the response of a unit cell under the same loading modes was conducted. The results in that work showed a negligible influence of misalignment for axial loading compared to other factors. However, as mentioned in Section 5-6, here the misalignment turned out to be a very important factor in axial loading. The reason for this contradiction can be closely related to the difference between the material models that have been used. Komeili and Milani (2009) use the homogenized material properties by Peng and Cao (2002), in which the magnitude of shear and transverse stiffness is in the same order of magnitude as the longitudinal stiffness (Table 3-1). On the other hand, the current model uses the Gasser et al. (2000) approach, whereby the shear stiffness vanishes to zero and transverse stiffness is nonlinearly related to axial and transverse strains where it has a very low magnitude at the beginning of deformation. The high value of shear stiffness allows for the yarns to act more like a beam and gives them the ability to bear shearing during axial loading of the unit cell. Even though they are in an inclined angle with the direction of loading, they can still show a high resistance against undulation. On the other hand, for the current model, because of low shear stiffness, yarns tend to first be aligned along the loading direction with no high resistance and then other material factors come to effect. To prove this point, a new set of runs was conducted. Namely, two bi- axial runs with A-B-C+D- and A-B-C-D- levels were selected to study the effect of misalignment under different shear stiffness values. The new set of runs was implemented under the circumstance of very low (G = 15 MPa) and very high (G = 7500 MPa) shear stiffness; the latter is almost in the order of magnitude of longitudinal stiffness and very close to the value used by Komeili and Milani (2009). Results show a higher domination of misalignment in the presence of lower shear stiffness (Figure 5-18). Interestingly, the misalignment effect almost vanishes for the high value of shear stiffness and this proves that the choice of the material model can significantly affect the results of parametric studies at meso-level. 106 Figure 5-18: Effect of shear stiffness on the importance of misalignment in axial loading 5.7.2 Effect of fricion in axial loading Another surprising result from the performed sensitivity analysis on the material factors is the point that friction turned out to be a crucial factor under all of the loading modes. For shear loading there was no difficulty in the physical interpretation of this, because in the shear frame test, yarns tended to rotate on each other rather than elongating or crushing. This is the dominant 0 5 10 15 20 25 30 35 40 45 0 0.002 0.004 0.006 N or m al iz ed re ac tio n Fo rc e (N /m m ) G = 7500 MPa A-B-C+D- A-B-C+D+ 0 5 10 15 20 25 0 0.002 0.004 0.006 N or m al iz ed R ea ct io n Fo rc e (N /m m ) Strain G = 15.0 MPa A-B-C+D- A-B-C+D+ 107 source of resisting force up to the locking point. However, for the axial mode it is expected that the longitudinal stiffness is the leading factor. This can be related back to the applied loading magnitude in the axial tests. Indeed, the magnitude of strain that was used here is close to what is reported in tensile tests by Buet-Gautier and Boisse (2001). At this level of strain (ߝ ൎ 0.006) yarns may not have come to full tension yet, hence the undulation and relative movements of the yarns are two dominant mechanisms for reaction force (Figure 5-19). In this case, friction can play an important role through the reaction forces at the contact surfaces while yarns moving on each other. Figure 5-19: Superposition of deformed and un-deformed unit cell in bi-axial loading. Solid shape is the deformed and shadow shape is the un-deformed unit cell Nevertheless, in higher strain values there should be a higher contribution from the longitudinal stiffness. In this case, the effect of friction almost vanishes and the longitudinal stiffness becomes the dominant factor. Figure 5-20 shows the result of a new factorial design that was conducted with two factors: longitudinal stiffness (A) and friction (B). This time, strain magnitude in the unit cell is assigned to be 0.0155, which is almost three times the value of the original. Plots in this figure prove that although with lower strain values as the applied load, friction is a considerably significant factor, in larger strain values it loses its significance and the longitudinal stiffness become the dominant factor. Nonetheless, because woven fabrics’ main application is in fabrication processes where large shear deformations and less axial deformations are present (like stamping), the range of strain used in the original factorial design is reasonable for practical use. 108 Figure 5-20: Effect of friction and longitudinal stiffness in larger strain values; A (longitudinal stiffness) and B (friction coefficient) 5.8 Summary First, a nominal unit cell model with material properties extracted from experimental results in Chapter 3 was constructed. Then, the geometrical and material parameters of the model were varied to run a series of two-level factorial designs on geometrical and material uncertainties. For each run, two indicators were defined that show the final values of normalized force at the highest value of deformation and the area under the force-deformation curve. The difference between these two indicators for each run with those of nominal run was used in the factorial design analyses. Finally, the main effects and parameter interactions were analyzed using percentage contributions, half-normal probability plots and Lenth’s criteria. A more detailed analysis on the effects of misalignment and friction in axial loading was presented and it was concludes that the choice of yarn constitutive model as well as the range of deformation can closely influence these effects. 0 20 40 60 80 100 120 0 0.005 0.01 0.015 N or m al iz ed R ea ct io n Fo rc e (N /m m ) Strain Basic A+B+ A+B- A-B+ A-B- 109 Chapter 6: Conclusions and remarks This chapter begins with a summary of the methodology used in this study. Afterwards, the main conclusions observed in preceding chapters are presented. In the final section of the chapter, suggestions on potential future work in the area of fabric meso-modeling and design are presented. 6.1 Summary of methodology A special finite element model using the explicit and implicit integrator of the Abaqus package was developed for studying the mechanical behaviour of fabric unit cells under three basic loading modes: uni-axial, bi-axial and pure shear. The material constitutive model was specifically developed for dry fabrics at meso-level. It uses a modified type of objective derivative, which is based on the direction of fibres for updating stress tensor, whereas Abaqus itself uses the Green-Naghdi or Jaumann objective derivative by default. The user defined material constitutive model has been incorporated into Abaqus by a FORTRAN subroutine to represent correct yarn behaviour. Subsequently, by applying an inverse method, the material properties are found to make the finite element model capture the published experimental results. Two separate series of uncertainty sensitivity analyses have been conducted on the geometrical and material factors in meso-level. Four geometrical factors, the yarn spacing, width, height and fibre misalignment have been considered in a factorial design and the most significant factors on the response of a plain weave unit cell were identified. Similarly, a factorial design on the yarn longitudinal and transverse stiffnesses, as well as the friction coefficient between yarns has been conducted. The effect and percentage contribution of factors and the interactions have been identified using the Design Expert software. Eventually, these effects and interactions were analysed using half-normal probability plots and Lenth’s ME and SME criteria. Because of the symmetry of misalignment under axial loading modes, there is no difference between +5° and -5° misalignment, as they both represent the same unit cell behaviour. Therefore, the lower bound of misalignment in axial loading is considered to be 0° during factorial design. In other words, the case of no misalignment is considered to be the value of the lower bound for this factor in axial modes. 110 6.2 Summary of inverse method and sensitivity analysis results An inverse method was used for finding the material properties of yarns on meso-level from experimental tests in macro-level samples. The fact that fitted material properties agreeably matched experimental curves proves that experimental data on fabrics can be advantageously used along with an inverse method to find the unit cell material properties in meso-level, which is hard to measure directly. Moreover, the observations from the results of factorial designs provided valuable insight on the effects of uncertainty factors on the response of the fabric, and also on the obstacles of numerical modeling of unit cells. Two-level factorial design has been conducted using finite element modeling of a plain weave unit cell and incorporating uncertainty factors on yarns geometry and material properties. For axial modes, all finite element samples were tested up to an identical value of strain (5.5848 ൈ 10ିଷ). Similarly, for shear, samples were tested under an equal stretch value (1.3). A normalized force value was used for assessing the response of the material under different uncertainty factor conditions. In axial modes, the load was normalized by the length of the sample and, for the shear mode, the absorbed energy per unit area of the fabric was used. Following from this, the normalized curves versus strain (for axial modes) or stretch (for shear mode) were plotted and the final values of normalized force (I∞) and the area under each curve (I1) were considered as outputs for the factorial design runs. Finally, the difference between I∞ and I1 of each run with the nominal unit cell was fed into the Design Expert package for analysis. The summary of the analysis results is as follows. It is worth noting that all the following conclusions are limited to a balanced plain weave glass fabric under the three basic loading modes (uni-axial tension, equi-biaxial tension and picture frame shear) with specified range of deformations. For all factorial design cases, the yarn material model has been based on the framework discussed in Section 3.4.2. 6.2.1 Geometrical sensitivity Geometrical sensitivity analysis results show that misalignment is a very important factor regardless of the loading mode. This can be physically described by considering the kinematics of deformation in fabrics. In axial loading, if yarns are located along the direction of loading, from the initiation of deformation, the main resisting force is from the stretching of the yarns. On 111 the other hand, in the presence of misalignment, yarns have to move and become parallel to the direction of loading before any stretching happens. In shear mode, kinematic of deformation can be separated into two stages. The first step is from the beginning of deformation to the point that adjacent yarns come to a side to side contact (also known as “locking”). The second stage is after the initiation of the above contact when pressure between the side faces of yarns develops. During the first step, the movement of yarns on each other and their friction is the main source of reaction force, whereas in the latter step, the pressure between the locked yarns plays the key role. It should be mentioned that much higher reaction forces are expected in the second step because of the locking. Considering the configuration of yarns in a unit cell and the kinematics of picture frame test, misalignment is the most important factor in fabric shear loading because it can expedite or postpone the locking. Yarn height and width are also found to be important factors in axial modes. Physically, increasing the yarns height or width means increasing its cross sectional area, which leads to a higher reaction force, due to the stretching of a bigger section of the yarn. Although the numerical value of the area change by varying these factors is very small, a noticeable difference between their contribution percentages during the parametric study on the response of the unit cell was observed and it suggested another influential mechanism. This mechanism relates to the undulation of the yarns. A larger height in a yarn means the curvier the axial line of the yarn. Accordingly, yarns with higher heights have to undergo higher strain before they start stretching in the yarn and this can closely affect the reaction force of the unit cell. Yarn spacing presented itself as a factor with marginal importance. Its influence on the fabric response can be related to the undulation and rotation of misaligned yarns along the loading direction (this is why spacing-misalignment interaction was seen to be an important factor in uni-axial mode). Nonetheless, its influence is not dramatically high. In the shear mode, the yarn width is recognized to be the second important factor after misalignment. This can again be related to its effect on the locking phenomenon. Wider yarns mean that the side contact between them occur sooner. Thus, it is expected to be significant in the variation of reaction force. Moreover, the interaction from misalignment and yarn width in shear mode is observed to be important. This suggests that these two factors could be the only significant geometrical uncertainty factors (statistically) under pure shear loading of a plain weave specimen. 112 6.2.2 Material sensitivity From the percentage contribution of factors obtained in factorial design, it is clear that friction is the dominant material factor in all modes. Moreover, longitudinal stiffness in axial mode, as well as the transverse stiffness in shear mode, are the second and third most important factors. However, in the uni-axial mode, the effects exceed neither Lenth’s ME nor SME criterion. It is also evident that half-normal probability plots show that factors and interactions are almost on a straight line. Thus, it may be concluded that although a noticeable percentage contribution from some factors in uni-axial loading exists, no material factor or interaction is statistically influential in the fabric response and they may be considered as random noise factors in the system. There is no doubt on the significance of longitudinal stiffness in both axial modes and transverse stiffness in bi-axial loading, because of the resistance from yarns against undulation through the pressure in their side contact surfaces. However, friction does not seem to be very influential against the stretching force of yarns. Indeed, at the applied level of strain in the conducted design of the experiment, yarns are still in their undulation step and contact interaction is the dominant source of reaction force. Nonetheless, increasing the strain magnitude to higher values showed a vanishing effect of friction and an increasing effect of the longitudinal stiffness. In the shear mode, friction clearly plays an important role, because of its effect before and after the formation of side contacts (locking) in the yarns. In the first stage of deformation, when yarns are rotating on each other, friction plays an important role in creating reaction force. In second stage, when yarns are in side contact and severe pressure exists between them, friction becomes even more important on contact surfaces. Additionally, in the second stage, it can be observed that due to the contact condition, the transverse stiffness becomes highly important. 6.3 Recommended future work This work presented an analysis on the finite element of fabric unit cell under uncertainty factors and a design of experiment analysis was conducted to find the influential factors. However, clearly there is still a very broad range of studies and research in this field that could not be covered here. The following are some suggested future research topics. 113 • Incorporate the current results into a macro-level homogenized material constitutive model. Factorial designs in the work revealed the significant factors in meso-level and their contributions toward the response. This can effectively be used in homogenization techniques, which use fabric models from meso-level to postulate constitutive models at macro-level. Defining response surfaces and deriving material models in macro-level as a function of uncertainty factors can result in design equations, which can advantageously capture uncertainties in large scale macro simulations. Subsequently, these set of homogenized macro-level material properties may advantageously be used in related simulations of composite parts. They can also be incorporated into commercial finite element packages for simulating woven fabric composite forming processes. • In this work the focus was on three basic modes of loading in a predefined deformation range. Similarly, however, a more general study can be conducted to consider the effect of combined loading. It could also be worthwhile to study the effect of loading magnitude on the sensitivity analysis results, because, as it was observed here, some factors can either loosen or strengthen their significance by increasing the loading magnitude (here the applied displacement range). • Although the basic results for the nominal unit cell were compared to experimental data and they coincided with the published results, no comprehensive experimental test data is available in the literature specifically on uncertainty analysis. There is a major need for more experimental measurements in this field both for non-repeatability of raw data and validation of sensitivity analysis results. This is especially the case if a comprehensive macro-level constitutive model is to be generated for fabrics and these results can be practically used for manufacturing process simulations or quality control applications. Thus, an experimental survey on the range of uncertainty factors is of high priority. • Most of the observations and conclusions in this work were limited to the 2D type of textile composites. Even the conclusions that may pertain to a similar fabric (e.g. two- dimensional twill weave) should be verified by a similar method. Thus, for new fabric types and material properties (especially if a very different yarn material has been used), similar factorial designs can be conducted to find significant uncertainty factors and their levels of importance. It would be of interest to investigate how the same geometrical and material uncertainty factors play roles in different fabric architectures. 114 Bibliography Aitharaju, V, & Averill, R (1999) Three-dimensional properties of woven-fabric composites. Composites science and technology 59. Badel, P, Gauthier, S, Vidal-Sallé, E, & Boisse, P (2009) Rate constitutive equations for computational analyses of textile composite reinforcement mechanical behaviour during forming. Composites Part A: Applied Science and Manufacturing 40: 997-1007. Badel, P, Maire, E, Vidal-Sallé, E, & Boisse, P (2008) Computational determination of the mechanical behavior of textile composite reinforcement. Validation with x-ray tomography. International Journal of Material Forming 1: 823-826. Badel, P, Vidal-Sallé, E, & Boisse, P (2007) Computational analysis of the mechanical behavior of textile composite reinforcements. ESAFORM Conference on Material Forming, (Cueto, E & Chinesta, F), pp. 1023-1026. American Institute of Physics. Badel, P, Vidalsalle, E, & Boisse, P (2008) Large deformation analysis of fibrous materials using rate constitutive equations. Computers & Structures 86: 1164-1175. Badel, P, Vidalsalle, E, Maire, E, & Boisse, P (2008) Simulation and tomography analysis of textile composite reinforcement deformation at the mesoscopic scale. Composites Science and Technology 68: 2433-2440. Bakhvalov, N, & Panasenko, G (1989) Homogenisation: averaging processes in periodic media: mathematical problems in the mechanics of composite materials. Springer. Ballhause, D, Konig, M, Kroplin, B, & Source, C (2006) Modelling of woven fabrics with the Discrete Element Method. European Conference on Computational Mechanics. Lisbon, Portugal. Ballhause, D, König, M, & Kröplin, B (2005) A microstructure model for fabric- reinforced membranes based on discrete element modelling. International conference on Textile Composites and Inflatable Structures, (Onate, E & Kroplin, B), pp. 1-10. CIMNE, Barcelona. Bigaud, D, & Hamelin, P (1997) Mechanical properties prediction of textile-reinforced composite materials using a multiscale energetic approach. Composite structures 38: 361–371. Boisse, P, Akkerman, R, Cao, J, Chen, J, Lomov, S, & Long, A (2007) Composites Forming. Advances in Material Forming - Esaform 10 years on. Springer, Paris. Boisse, P, Badel, P, & Vidal-salle, E (2007) Computational determination of in-plane shear mechanical behaviour of textile composite reinforcements. Computational Materials Science 40: 439-448. 115 Boisse, P, Borr, M, Buet, K, & Cherouat, A (1997) Finite element simulations of textile composite forming including the biaxial fabric behaviour. Composites. Part B, Engineering 28: 453–464. Boisse, P, Buet, K, Gasser, A, & Launay, J (2001) Meso/macro-mechanical behaviour of textile reinforcements for thin composites. Composites Science and Technology 61: 395–401. Boisse, P, Gasser, A, Hagege, B, & Billoet, J (2005) Analysis of the mechanical behavior of woven fibrous material using virtual tests at the unit cell level. Journal of Materials Science 40: 5955-5962. Boisse, P, Zouari, B, & Gasser, A (2005) A mesoscopic approach for the simulation of woven fibre composite forming. Composites Science and Technology 65: 429-436. Boubaker, B, Haussy, B, & Ganghoffer, J (2007) Consideration of the yarn–yarn interactions in meso/macro discrete model of fabric. Part I: Single yarn behaviour. Mechanics Research Communications 34: 359-370. Boubaker, B, Haussy, B, & Ganghoffer, J (2007) Consideration of the yarn–yarn interactions in meso/macro discrete model of fabric. Part II: Woven fabric under uniaxial and biaxial extension. Mechanics Research Communications 34: 371-378. Bu, H, Wang, J, & Huang, X (2009) Fabric defect detection based on multiple fractal features and support vector data description. Engineering Applications of Artificial Intelligence 22: 224-235. Buet-Gautier, K, & Boisse, P (2001) Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements. Experimental Mechanics 41: 260–269. Cao, J et al. (2008) Characterization of mechanical behavior of woven fabrics: Experimental methods and benchmark results. Composites Part A: Applied Science and Manufacturing 39: 1037-1053. Carvelli, V, & Poggi, C (2001) A homogenization procedure for the numerical analysis of woven fabric composites. Composites Part A 32: 1425–1432. Chen, D, & Cheng, S (1993) Analysis of composite materials: A micromechanical approach. Journal of Reinforced Plastics and Composites 12: 1323-1338. Chen, J, Lussier, D, Cao, J, & Peng, X (2001) Materials characterization methods and material models for stamping of plain woven composites. International Journal of Forming Processes 4: 269–284. Chung, PW, & Tamma, KK (1999) Woven fabric composites—developments in engineering bounds, homogenization and applications. International Journal for Numerical Methods in Engineering 45: 1757-1790. 116 Cox, B, Carter, W, & Fleck, N (1994) Binary model of textile composites - I. Formulation. Acta metallurgica et materialia 42: 3463-3479. Creighton, C, Sutcliffe, M, & Clyne, T (2001) A multiple field image analysis procedure for characterisation of fibre alignment in composites. Composites Part A 32: 221–229. Daniel, C (1959) Use of Half-Normal Plots in Interpreting Factorial Two Level Experiments. Technometrics 1: 311-342. Dienes, JK (1979) On the analysis of rotation and stress rate in deforming bodies. Acta Mechanica 32: 217-232. Durville, D (2005) Approach of the constitutive material behaviour of textile composites through simulation. International conference on Textile Composites and Inflatable Structures, pp. 1-10. France. Elleithy, R (2000) The hierarchical structure and flexure behavior of woven carbon fiber epoxy composite. Polymer Composites 27. Felippa, CA (2001) Intorduction to finite element methods. Depertment of Aerospace Engineering Science, University of Colorado, Boulder. Flanagan, D, & Belytschko, T (1981) A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. International Journal for Numerical Methods in Engineering. Gasser, A, Boisse, P, & Hanklar, S (2000) Mechanical behaviour of dry fabric reinforcements. 3D simulations versus biaxial tests. Computational Materials Science 17: 7–20. Green, A, & Naghdi, P (1965) A general theory of an elastic-plastic continuum. Archive for rational mechanics and analysis 18. Guagliano, M, & Riva, E (2001) Mechanical behaviour prediction in plain weave composites. Journal of strain analysis for engineering design 36: 153-162. Halpin, JC, & Tsai, SW (1969) Effective environmental factors on composite materials. Hamada, M, & Balakrishnan, N (1998) Analyzing unreplicated factorial experiments: A review with some new proposals. Statistica Sinica 8: 1–28. Hamila, N, & Boisse, P (2007) A Meso–Macro Three Node Finite Element for Draping of Textile Composite Preforms. Applied Composite Materials 14: 235-250. Harrison, P, & Johnson, E (1996) A mixed variational formulation for interlaminar stresses in thickness-tapered composite laminates. International Journal of Solids and Structures 33: 2377-2399. 117 Hivet, G, & Boisse, P (2005) Consistent 3D geometrical model of fabric elementary cell. Application to a meshing preprocessor for 3D finite element analysis. Finite Elements in Analysis and Design 42: 25-49. Hivet, G, & Boisse, P (2008) Consistent mesoscopic mechanical behaviour model for woven composite reinforcements in biaxial tension. Composites Part B: Engineering 39: 345- 361. Hivet, G, Launay, J, Gasser, A, Daniel, JL, & Boisse, P (2002) Mechanical Behavior of Woven Composite Reinforcements While Forming. Journal of Thermoplastic Composite Materials 15: 545-555. Hughes, T, & Winget, J (1980) Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. International journal for numerical methods in. Ishikawa, T, & Chou, T (1983) On-dimensional micromechanical analysis of woven fabric composites. AIAA Journal 21: 1714-1721. Ishikawa, T, & Chou, T (1982) Stiffness and strength behaviour of woven fabric composites. Journal of Materials Science 17: 3211–3220. Ishikawa, T, Matsushima, M, Hayashi, Y, & Chou, T (1985) Experimental Confirmation of the Theory of Elastic Moduli of Fabric Composites. Journal of Composite Materials 19: 443- 458. Ivanov, I, & Tabiei, A (2002) Flexible Woven Fabric Micromechanical Material Model with Fiber Reorientation. Mechanics of Advanced Materials and Structures 37 - 51. Jaumann, G (1911) Geschlossenes System physikalischer und chemischer Differentialgesetze. Sitzgsber. Akad. Wiss. Wien (IIa). Jiang, Y, Tabiei, A, & Stimitses, J (2000) Mircomechanics-based plain weave fabric composites constitutuve equations for local and global analysis. Journal of Composite Science and Technology 60: 553-560. Kaw, AK (1997) Mechanics of composite materials. New York: CRC Press. Kawabata, S, Niwa, M, & Kawai, H (1973) Finite-deformation theory of plain-weave fabrics - 1. The biaxial-deformation theory. Journal of the Textile Institute 64: 21-46. Kawabata, S, Niwa, M, & Kawai, H (1973) Finite-deformation theory of plain-weave fabrics - 2. The uniaxial-deformation theory. Journal of the Textile Institute 64: 47-61. Kawabata, S, Niwa, M, & Kawai, H (1973) Finite-deformation theory of plain-weave fabrics - 3. The shear-deformation theory. Journal of the Textile Institute 64: 62-85. 118 Kawabata, S (1989) Nonlinear mechanics of woven and knitted materials. Elsevier Science Publishers, Textile Structural 46: 285-293. King, M, Jearanaisilawong, P, & Socrate, S (2005) A continuum constitutive model for the mechanical behavior of woven fabrics. International Journal of Solids and Structures 42: 3867-3896. Komeili, M, & Milani, A (2009) On the Effect of Uncertainty Factors on Mechanical Behavior of Woven Fabric Composites at Meso-Level. ASME International Mechanical Engineering Congress & Exposition. Lake Buena Vista. Kuhn, JL, Haan, SI, & Charalambides, PG (1999) A Semi-Analytical Method for the Calculation of the Elastic Micro-Fields in Plain Weave Fabric Composites Subjected to In-Plane Loading. Journal of Composite Materials 33: 221-266. Lancaster, P, & Tismenetsky, M (1985) The theory of matrices: with applications. Lebrun, G, Bureau, M, & Denault, J (2003) Evaluation of bias-extension and picture- frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics. Composite Structures 61: 341-352. Lenth, R (1989) Quick and easy analysis of unreplicated factorials. Technometrics 469– 473. Lim, T, & Ramakrishna, S (2002) Modelling of composite sheet forming: a review. Composites Part A 33: 515–537. Lin, H, Clifford, MJ, Long, AC, & Sherburn, M (2009) Finite element modelling of fabric shear. Modelling and Simulation in Materials Science and Engineering 17: 015008. Lomov, S, & Verpoest, I (2006) Model of shear of woven fabric and parametric description of shear resistance of glass woven reinforcements. Composites Science and Technology 66: 919-933. Lomov, S, Huysmans, G, Luo, Y, Parnas, R, Prodromou, A, Verpoest, I, & Phelan, F (2001) Textile composites: modelling strategies. Composites part A 32: 1379–1394. Lomov, S, Ivanov, D, Verpoest, I, Zako, M, Kurashiki, T, Nakai, H, & Hirosawa, S (2007) Meso-FE modelling of textile composites: Road map, data flow and algorithms. Composites Science and Technology 67: 1870-1891. Lua, J (2007) Thermal–mechanical cell model for unbalanced plain weave woven fabric composites. Composites Part A: Applied Science and Manufacturing 38: 1019-1037. Lussier, D, & Chen, J (2002) Material characterization of woven fabrics for thermoforming of composites. Journal of Thermoplastic Composite Materials 15: 497. 119 Mase, GT, & Mase, GE (1999) Continuum mechanics for engineers. CRC Press LLC. McGlockton, M, Cox, B, & McMeeking, R (2003) A Binary Model of textile composites: III high failure strain and work of fracture in 3D weaves. Journal of the Mechanics and Physics of Solids 51: 1573-1600. Mcbride, TM, & Chen, J (1997) Unit-cell geometry in plain-weave during shear deformations fabrics. Composites Science and Technology 57: 345-351. Milani, A, & Nemes, J (2004) An intelligent inverse method for characterization of textile reinforced thermoplastic composites using a hyperelastic constitutive model. Composites Science and Technology 64: 1565-1576. Milani, A, Nemes, J, Abeyaratne, R, & Holzapfel, G (2007) A method for the approximation of non-uniform fiber misalignment in textile composites using picture frame test. Composites Part A: Applied Science and Manufacturing 38: 1493-1501. Milani, A, Nemes, J, Lebrun, G, & Bureau, M (2009) A comparative analysis of a modified picture frame test for characterization of woven fabrics. Polymer Composites 31: 561– 568. Montgomery, DC (2000) Design and Analysis of Experiments. Wiley. Naik, N, & Shembekar, P (1992) Elastic Behavior of Woven Fabric Composites: I-- Lamina Analysis. Journal of Composite Materials 26: 2196-2225. Pagano, N (1978) Stress fields in composite laminates. International Journal of Solids and Structures 14: 385-400. Page, J, & Wang, J (2002) Prediction of shear force using 3D non-linear FEM analyses for a plain weave carbon fabric in a bias extension state. Finite Elements in Analysis and Design 38: 755-764. Peng, X, & Cao, J (2002) A dual homogenization and finite element approach for material characterization of textile composites. Composites Part B 33: 45–56. Peng, X, & Cao, J (2000) Numerical determination of mechanical elastic constants of textile composites. PROCEEDINGS-AMERICAN SOCIETY FOR COMPOSITES, p. 677–688. Peng, X, Cao, J, Chen, J, Xue, P, Lussier, D, & Liu, L (2004) Experimental and numerical analysis on normalization of picture frame tests for composite materials. Composites Science and Technology 64: 11-21. Pian, T, & Chen, D (1983) On the suppression of zero energy deformation modes. International Journal for Numerical Methods in Engineering. 120 Potluri, P, & Manan, A (2007) Mechanics of non-orthogonally interlaced textile composites. Composites Part A: Applied Science and Manufacturing 38: 1216-1226. Potluri, P, Perezciurezu, D, & Ramagulam, R (2006) Measurement of meso-scale shear deformations for modelling textile composites. Composites Part A: Applied Science and Manufacturing 37: 303-314. Reissner, E (1950) On a variational theorem in elasticity. J. Math. Phys 29: 90-95. Roy, AK, & Sihn, S (2001) Development of a three-dimensional mixed variational model for woven composites. I. Mathematical formulation. Science 38: 5935-5947. SIMULIA (2008) Abaqus 6.8-1 Documentation. Sejnoha, M, & Zeman, J (2008) Micromechanical modeling of imperfect textile composites. International Journal of Engineering Science 46: 513-526. Seng, LK, Earn, TT, & Shim, VP (1997) Hygrothermal analysis of woven-fabric composite plates. Composites. Part B, Engineering 28: 573–581. Sherburn, M (2007) Geometric and Mechanical Modelling of Textiles. Microscopy. Sihn, S, & Roy, AK (2001) Development of a three-dimensional mixed variational model for woven composites. II. Numerical solution and validation. International Journal of Solids and Structures 38: 5949-5962. Sihn, S, Iarve, EV, & Roy, AK (2004) Three-dimensional stress analysis of textile composites: Part I. Numerical analysis. International Journal of Solids and Structures 41: 1377- 1393. Sihn, S, Iarve, EV, & Roy, AK (2004) Three-dimensional stress analysis of textile composites: Part II. Asymptotic analysis. International Journal of Solids and Structures 41: 1395-1410. Skordos, A, & Sutcliffe, M (2008) Stochastic simulation of woven composites forming. Composites Science and Technology 68: 283-296. Tabiei, A, & Ivanov, I (2003) Fiber Reorientation in Laminated and Woven Composites for Finite Element Simulations. Journal of Thermoplastic Composite Materials 16: 457-474. Tabiei, A, & Jiang, Y (1999) Woven fabric composite material model with material nonlinearity for nonlinear finite element simulation. International Journal of Solids and Structures. Tabiei, A, & Yi, W (2002) Comparative study of predictive methods for woven fabric composite elastic properties. Composite Structures 58: 149–164. 121 Tabiei, A, Jiang, Y, & Yi, W (2000) A novel micromechanics-based plain weave fabric composite constitutive model with material nonlinear behavior. AIAA Journal 38. Takano, N, Uetsuji, Y, Kashiwagi, Y, & Zako, M (1999) Hierarchical modelling of textile composite materials and structures by the homogenization method. Modelling and Simulation in Materials Science and Engineering 7: 207-231. Tan, P, Tong, L, & Steven, G (1997) Modelling for predicting the mechanical properties of textile composites—a review. Composites Part A 28: 903–922. Tanov, R, & Tabiei, A (2001) Computationally efficient micro-mechanical woven fabric constitutive. J Models Appl Mech 68: 553-560. Vandeurzen, P, Ivens, J, & Verpoest, I (1996) A three-dimensional micromechanical analysis of woven-fabric composites: II. elastic analysis. Composites Science and Technology 56: 1317-1327. Vandeurzen, P, Ivens, J, & Verpoest, I (1997) Microstress analysis in woven fabric composites using variational principles. Plastics, Rubber and Composites Processing and Applications 26: 438-446. Vandeurzen, P, Vandeurzen, P, Ivens, J, & Verpoest, I (1996) A three-dimensional micromechanical analysis of woven-fabric composites: I. Geometric analysis. Composites Science and Technology 56: 1303-1315. Whitcomb, J, Srirengan, K, & Chapman, CS (1994) Evaluation of homogenization for global/local stress analysis of textile composites. Collection of Technical Papers - AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 3: 1649-1663. Whitcomb, J (1989) Three-dimensional stress analysis of plain weave composites. Composite Materials: Fatigue and Fracture (Third Volume). Woo, K, & Whitcomb, JD (1993) Macro finite element using subdomain integration. Communications in Numerical Methods in Engineering, 9: 937-949. Xu, J, Cox, B, McGlockton, M, & Carter, W (1995) Binary model of textile composites - II. The elastic regime. Acta metallurgica et materialia 43: 3511-3524. Xue, P, Cao, J, & Chen, J (2005) Integrated micro/macro-mechanical model of woven fabric composites under large deformation. Composite Structures 70: 69-80. Xue, P, Peng, X, & Cao, J (2003) A non-orthogonal constitutive model for characterizing woven composites. Composites Part A: Applied Science and Manufacturing 34: 183-193. 122 Zeman, J, & Sejnoha, M (2004) Homogenization of balanced plain weave composites with imperfect microstructure: Part I-Theoretical formulation. International Journal of Solids and Structures 41: 6549-6571. 123 Appendix A: Transformation of second order tensors in vector form Transforming stress and strain vectors from one coordinate system to another via a rotation is theoretically derived from formulation of transformation of stress or strain tensors. Assume ۯ to be a second order tensor in three-dimensional space and દ to be the transformation tensor between base vectors of a current coordinate and another coordinate of interest. Then the stress or strain tensor in the rotated coordinate can be found by (Mase and Mase 1999): ۯԢ ൌ દ ۯ દ܂ ۯ ൌ દ܂ ۯԢ દ (A-1) in which, for this case, because the focus is on the transformation from the fibre frame to the software frame, it can be verified that: દ ൌ ൫. ൯ ٔ ൌ ൫൯ ٔ (A-2) Considering the symmetry in ۯ and ۯԢ, they are stored in the software as 6 by 1 vectors to reduce the computational cost. In that case, the conventional formula for tensor transformation cannot be used. However, a new transformation matrix based on the original relationship in tensorial transformation has been found as follows (Badel et al 2008a). ۏ ێ ێ ێ ێ ۍ ܣଵଵ ܣଶଶ ܣଷଷ ܣଵଶ ܣଵଷ ܣଶଷے ۑ ۑ ۑ ۑ ې ൌ ሾܶሿ. ۏ ێ ێ ێ ێ ێ ۍ ܣԢଵଵ ܣԢଶଶ ܣԢଷଷ ܣԢଵଶ ܣԢଵଷ ܣԢଶଷے ۑ ۑ ۑ ۑ ۑ ې (A-3) in which ሾܶሿ ൌ ሾܶሿ ሾܶሿ ሾܶሿ ሾܶሿ ൨ (A-4) ሾܶሿ ൌ ሺ ଵ݂ሻଵଶ ሺ ଶ݂ሻଵଶ ሺ ଷ݂ሻଵଶ ሺ ଵ݂ሻଶଶ ሺ ଶ݂ሻଶଶ ሺ ଷ݂ሻଶଶ ሺ ଵ݂ሻଷଶ ሺ ଶ݂ሻଷଶ ሺ ଷ݂ሻଷଶ ሾܶሿ ൌ 2 ሺ ଵ݂ሻଵ ሺ ଶ݂ሻଵ 2 ሺ ଵ݂ሻଵ ሺ ଷ݂ሻଵ 2 ሺ ଶ݂ሻଵ ሺ ଷ݂ሻଵ 2 ሺ ଵ݂ሻଶ ሺ ଶ݂ሻଶ 2 ሺ ଵ݂ሻଶ ሺ ଷ݂ሻଶ 2 ሺ ଶ݂ሻଶ ሺ ଷ݂ሻଶ 2 ሺ ଵ݂ሻଷ ሺ ଶ݂ሻଷ 2 ሺ ଵ݂ሻଷ ሺ ଷ݂ሻଷ 2 ሺ ଶ݂ሻଷ ሺ ଷ݂ሻଷ ሾܶሿ ൌ ሺ ଵ݂ሻଵ ሺ ଵ݂ሻଶ ሺ ଶ݂ሻଵ ሺ ଶ݂ሻଶ ሺ ଷ݂ሻଵ ሺ ଷ݂ሻଶ ሺ ଵ݂ሻଵ ሺ ଵ݂ሻଷ ሺ ଶ݂ሻଵ ሺ ଶ݂ሻଷ ሺ ଷ݂ሻଵ ሺ ଷ݂ሻଷ ሺ ଵ݂ሻଶ ሺ ଵ݂ሻଷ ሺ ଶ݂ሻଶ ሺ ଶ݂ሻଷ ሺ ଷ݂ሻଶ ሺ ଷ݂ሻଷ 124 ሾܶሿ ൌ ሺ ଵ݂ሻଶ ሺ ଶ݂ሻଵ ሺ ଵ݂ሻଵ ሺ ଶ݂ሻଶ ሺ ଵ݂ሻଶ ሺ ଷ݂ሻଵ ሺ ଵ݂ሻଵ ሺ ଷ݂ሻଶ ሺ ଶ݂ሻଶ ሺ ଷ݂ሻଵ ሺ ଶ݂ሻଵ ሺ ଷ݂ሻଶ ሺ ଵ݂ሻଷ ሺ ଶ݂ሻଵ ሺ ଵ݂ሻଵ ሺ ଶ݂ሻଷ ሺ ଵ݂ሻଷ ሺ ଷ݂ሻଵ ሺ ଵ݂ሻଵ ሺ ଷ݂ሻଷ ሺ ଶ݂ሻଷ ሺ ଷ݂ሻଵ ሺ ଶ݂ሻଵ ሺ ଷ݂ሻଷ ሺ ଵ݂ሻଷ ሺ ଶ݂ሻଶ ሺ ଵ݂ሻଶ ሺ ଶ݂ሻଷ ሺ ଵ݂ሻଷ ሺ ଷ݂ሻଶ ሺ ଵ݂ሻଶ ሺ ଷ݂ሻଷ ሺ ଶ݂ሻଷ ሺ ଷ݂ሻଶ ሺ ଶ݂ሻଶ ሺ ଷ݂ሻଷ and similarly, ۏ ێ ێ ێ ێ ێ ۍ ܣԢଵଵ ܣԢଶଶ ܣԢଷଷ ܣԢଵଶ ܣԢଵଷ ܣԢଶଷے ۑ ۑ ۑ ۑ ۑ ې ൌ ሾܶԢሿ. ۏ ێ ێ ێ ێ ۍ ܣଵଵ ܣଶଶ ܣଷଷ ܣଵଶ ܣଵଷ ܣଶଷے ۑ ۑ ۑ ۑ ې (A-5) in which ሾܶԢሿ ൌ ሾܶԢሿ ሾܶԢሿ ሾܶԢሿ ሾܶԢሿ ൨ (A-6) ሾܶԢሿ ൌ ሾܶሿ் ሾܶԢሿ ൌ 2 ሺ ଵ݂ሻଵ ሺ ଵ݂ሻଶ 2 ሺ ଵ݂ሻଵ ሺ ଵ݂ሻଷ 2 ሺ ଵ݂ሻଶ ሺ ଵ݂ሻଷ 2 ሺ ଶ݂ሻଵ ሺ ଶ݂ሻଶ 2 ሺ ଶ݂ሻଵ ሺ ଶ݂ሻଷ 2 ሺ ଶ݂ሻଶ ሺ ଶ݂ሻଷ 2 ሺ ଷ݂ሻଵ ሺ ଷ݂ሻଶ 2 ሺ ଷ݂ሻଵ ሺ ଷ݂ሻଷ 2 ሺ ଷ݂ሻଶ ሺ ଷ݂ሻଷ ሾܶԢሿ ൌ ሺ ଵ݂ሻଵ ሺ ଶ݂ሻଵ ሺ ଵ݂ሻଶ ሺ ଶ݂ሻଶ ሺ ଵ݂ሻଷ ሺ ଶ݂ሻଷ ሺ ଵ݂ሻଵ ሺ ଷ݂ሻଵ ሺ ଵ݂ሻଶ ሺ ଷ݂ሻଶ ሺ ଵ݂ሻଷ ሺ ଷ݂ሻଷ ሺ ଶ݂ሻଵ ሺ ଷ݂ሻଵ ሺ ଶ݂ሻଶ ሺ ଷ݂ሻଶ ሺ ଶ݂ሻଷ ሺ ଷ݂ሻଷ ሾܶԢሿ ൌ ሾܶሿ் 125 Appendix B: Matrix form of square root Matrix ሾ܁ሿ is said to be the square root of ሾ۰ሿ if ሾ܁ሿ. ሾ܁ሿ ൌ ሾ۰ሿ. For a diagonal matrix the square root can be found by taking square roots of its components. However, this does not pertain for a matrix of general form. If ሾ۰ሿ is a positive definite matrix there is a diagonal matrix ሾ۲ሿ that can satisfy the following equation (Lancaster and Tismenetsky 1985) ሾ۰ሿ ൌ ሾ܁ሿ ൌ ሾ܃ሿ. ሾ۲ሿ. ሾ܃ሿT (B-1) where ሾ܃ሿ is the matrix having eigenvectors of ሾ۰ሿ as columns. Thus, for finding the square root of ሾ۰ሿ it can be diagonalized using its eigenvectors and then transforming it back after finding the square root of the diagonlized matrix. In other words: ሾ۲ሿ ൌ ሾ܃ሿT. ሾ۰ሿ. ሾ܃ሿ (B-2) ሾ۲ሿ ൌ diagሾඥሾ۲ሿଵଵ ଶ , ඥሾ۲ሿଶଶ ଶ , … , ඥሾ۲ሿ୬୬ଶ ሿ (B-3) ሾ܁ሿ ൌ ሾ܃ሿ. ሾ۲ሿ. ሾ܃ሿT (B-4) This method can be simply applied by programming in a subroutine as for the current project it was a FORTRAN code in the Abaqus UMAT subroutine (available via contacting the authors).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Mechanical behavior of woven fabric composites under...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Mechanical behavior of woven fabric composites under meso-level uncertainties : modeling and sensitivity… Komeili, Mojtaba 2010
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Mechanical behavior of woven fabric composites under meso-level uncertainties : modeling and sensitivity analysis |
Creator |
Komeili, Mojtaba |
Publisher | University of British Columbia |
Date Issued | 2010 |
Description | Unit cell modeling of the woven fabric composites is a strong tool for studying fabric behavior at meso-level. Among the suggested methods the three-dimensional finite element analysis is found to be promising and of greater interest to researchers in the field. However, due to the multi-scale nature and particular behavior of fibrous yarns, numerical procedures applicable to woven fabrics differ from conventional finite element routines. Moreover, from a practical point of view, most of the models in the literature focus on the ideal unit cells and excludes the meso-level uncertainty factors. On the other hand, experimental measurements indicate non-repeatability of test data which is a result of embedded inherent uncertainties. In this study, a finite element model capable of defining the representative volume element of a plain weave under homogenous loading has been created. Because of particular behavior of dry fabrics, a special constitutive material behaviour is defined via user-defined subroutines, which links to the solver of a finite element package (Abaqus). Material properties of yarns for dry glass fabrics are extracted by fitting their numerical responses to results of experiments. Subsequently, geometrical and material meso-level uncertainty factors are studied separately using two-level factorial designs. The output of runs are inserted into a commercial design package (Design Expert) for producing results in terms of probability plots, effects and percentage contribution of each factor and their interactions. Result shows that depending on the loading type, there are factors that show significant contribution toward the final response, whereas others are quite negligible. More elaboration on the design results provides informative conclusions regarding modeling of unit cells, as well as their behavior during different loading steps. These conclusions can be used in more comprehensive unit-cell homogenization formulations, as well as for defining appropriate tolerances for meso-level defects in woven fabric composites. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0071372 |
URI | http://hdl.handle.net/2429/29237 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Engineering, School of (Okanagan) |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCO |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_Spring_Komeili_Mojtaba.pdf [ 3.6MB ]
- Metadata
- JSON: 24-1.0071372.json
- JSON-LD: 24-1.0071372-ld.json
- RDF/XML (Pretty): 24-1.0071372-rdf.xml
- RDF/JSON: 24-1.0071372-rdf.json
- Turtle: 24-1.0071372-turtle.txt
- N-Triples: 24-1.0071372-rdf-ntriples.txt
- Original Record: 24-1.0071372-source.json
- Full Text
- 24-1.0071372-fulltext.txt
- Citation
- 24-1.0071372.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0071372/manifest