UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A multiscale analysis of forming induced wrinkles in woven composite preforms Hosseini, Abbas 2018

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2018_may_hosseini_abbas.pdf [ 3.42MB ]
Metadata
JSON: 24-1.0363919.json
JSON-LD: 24-1.0363919-ld.json
RDF/XML (Pretty): 24-1.0363919-rdf.xml
RDF/JSON: 24-1.0363919-rdf.json
Turtle: 24-1.0363919-turtle.txt
N-Triples: 24-1.0363919-rdf-ntriples.txt
Original Record: 24-1.0363919-source.json
Full Text
24-1.0363919-fulltext.txt
Citation
24-1.0363919.ris

Full Text

I  A MULTISCALE ANALYSIS OF FORMING INDUCED WRINKLES IN WOVEN COMPOSITE PREFORMS  by  ABBAS HOSSEINI   A DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY  in   THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES  (Mechanical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  February 2018  © Abbas Hosseini, 2018  ii  Abstract A multi-scale analysis of forming induced wrinkles in woven composite preforms is conducted. The main objective of the present thesis is to propose an analytical model for the prediction of a class of forming induced wrinkles. In this research, wrinkling phenomenon in woven fabrics is divided into two distinct classes based on the nature and mechanism of formation, namely shear wrinkling (i.e., developed as a result of large shear deformation) and non-shear wrinkling (i.e., developed due to longitudinal contraction of the length of the fabric strip). Shear wrinkling as one of the most critical defects in the draping of woven fabrics is investigated in this study. To develop an analytical model for the prediction of shear wrinkling, an in-depth investigation is conducted to scrutinize the mesoscale mechanisms of woven fabrics during shear deformation. It is revealed that shear of woven fabrics occurs through two distinct deformation modes (i.e., trellising mode and pure multi-scale shear mode) depending on the boundary conditions applied to the fabric. In accordance with the deformation mode of woven fabrics, analytical equations are derived to predict the onset of shear wrinkling of the fabrics. Experimental investigations are conducted to assess the accuracy of the analytical model. Afterward, the analytical model for the prediction of shear wrinkling is implemented into ABAQUS Finite Element (FE) package as wrinkling indicator/predictor to improve the accuracy of the FE results. Furthermore, a parametric study is carried out to determine the effect of fabric and process parameters on the resistance of woven fabrics against wrinkling.   iii  Lay Summary Textile Structural Composites (TSC) have extensively gained attention in a wide variety of applications as enhanced multi-purpose material systems. Especially, woven fabric reinforced composites, as a class of TSCs, are commonly used in industrial sectors ranging from aerospace to marine and energy, owing to their superior formability and high out-of-plane stiffness and strength. Despite the satisfactory formability of woven fabrics compared to their unidirectional counterparts, there still exist a number of modes of a defect such as wrinkling that restrict the effortless forming of this class of textiles. There are two different types of forming induced wrinkles in woven composite preforms, namely non-shear wrinkles and shear wrinkles. Shear wrinkling, which is caused due to large fabric shear deformation, is receiving much attention in the current literature as one of the most concerning modes of defect in composites manufacturing. Accordingly, the present study is conducted to specifically investigate shear wrinkling behavior of woven fabrics.      iv  Preface This thesis is comprised of a number of collaborative research projects. Chapters 1 and 2 of the thesis are a result of a close collaboration between our team at UBC and Prof. Michael Sutcliffe at University of Cambridge. This collaboration took place in the Cambridge Center for Micromechanics, where I worked on the project under the supervision of Prof. Sutcliffe. These two chapters have led to the following paper: - A. Hosseini, M.P.F. Sutcliffe, F. Sassani, F.K. Ko “Forming Induced Wrinkles in Composite Preforms: A Review”. Under Submission, 2018. To prepare the manuscript of the above-mentioned paper, the professors were more focused on the framework and the topics that needed to be discussed, while I conducted the literature review and wrote the draft of the paper. Meanwhile, our team at UBC Vancouver was in a very close collaboration with Prof. Abbas Milani’s research group at UBC Okanagan. This collaboration started from CANCOM Conference held in Edmonton in August 2015, where I met Dr. Masoud Haghi Kashani. At the time, he was a PhD student at UBC Okanagan, working on shear deformation of woven fabrics. In his research, his primary focus was on the shear behavior of woven fabrics prior to the locking point, while that of mine was predominantly concentrated on the shear behavior of woven fabrics in the post-locking stage – from locking to wrinkling. Our discussion in Edmonton revealed that we have a cornerstone in our research that could complement one another’s research. Thereafter, through the Composites Research Network at the Vancouver and Okanagan nodes, we started working collaboratively under joint supervisions of Prof. Sassani and Prof. Ko (my co-supervisors) and Prof. Milani (Dr. Kashani’s supervisor). v  Eventually, the collaboration led to developing a new shear model for plain-woven fabrics. The shear model was employed in our study to predict the locking angle in the fabric (i.e., the most important trigger for shear wrinkling). One of the challenges we had at the time was that the shear model could not precisely predict locking angle of woven fabrics based on the assumptions available in the literature. I proposed to update the assumption by defining locking as the instant when the distance between adjacent parallel yarns is bounded by the thickness of interlacing yarn (i.e., in the past studies, locking was defined as the instant at when the distance between adjacent parallel yarns becomes zero). This proposed assumption notably improved the accuracy of the shear model for prediction of locking. The results of this part of the collaboration formed Chapter 3 of my thesis. Moreover, the results were published in the following journal: - M. Haghi Kashani, A. Hosseini, F. Sassani, F.K. Ko, A.S. Milani.  “The role of intra-yarn shear in integrated multi-scale deformation analyses of woven fabrics: A critical review”. Critical Reviews in Solid State and Materials Sciences, 2017 (5-Year Impact Factor: 5.453). http://dx.doi.org/10.1080/10408436.2017.1342597. After this stage, our collaboration with Prof. Milani and Dr. Kashani continued. We worked together on an analytical model to predict the onset of shear wrinkling by relying on the shear model developed in the former stage. I proposed an equivalent structure composed of bars and different types of springs which closely mimics the behavior of woven fabrics in the post-locking stage. Building on this structure, I derived an analytical equation to predict the onset of shear wrinkling of woven fabrics under trellising deformation mode (i.e., one of the possible shear deformation modes). Dr. Kashani actively contributed to the discussions that eventually, led to improving the presentation of the final equation. Moreover, he significantly vi  took part in the experimental validation of the proposed model. The results of this part of the project were the basis for Chapter 4 of my thesis. Moreover, the results were published in: - A. Hosseini, M. Haghi Kashani, F. Sassani, A.S. Milani, F.K. Ko. “A Mesoscopic Analytical Model to Predict the Onset of Wrinkling in Plain Woven Preforms under Bias Extension Shear Deformation”. Materials, 2017 (5-Year Impact Factor: 3.236). http://dx.doi.org/10.3390/ma10101184. - A. Hosseini, M. Haghi Kashani, F. Sassani, A.S. Milani, F.K. Ko. “A mechanism-based analytical model for predicting the onset of wrinkling in plain-woven composite preforms – A bridge from micro-scale to macro-scale”. The 21st annual Pacific Centre for Advanced Materials and Microstructures (PCAMM 2016) Conference, December 2016, Vancouver, Canada. To develop the analytical model for the prediction of shear wrinkling, we were in need of a correct understanding of the coupling effects between the mechanical properties of woven fabrics. Accordingly, an experimental study was conducted in which Dr. Kashani was the primary contributor, while I actively contributed to the discussions. The results of this part of the project demonstrated that there exists a coupling between tension-shear properties of woven fabrics which markedly affect shear wrinkling behavior of the fabrics. Noticing the inherent couplings in woven fabrics led us to the fact that the analytical model for predicting shear wrinkling should be a tension-shear-coupling-integrated equation. Integrating this effect into the analytical model remarkably improved the effectiveness of the model. The following paper is the outcome of this phase of the project: - M. Haghi Kashani, A. Hosseini, F. Sassani, F.K. Ko, A.S. Milani. “Understanding Different Types of Coupling in Mechanical Behaviour of Woven Fabric Reinforcements: A Critical Review and Analysis”. Composite Structures, 2017 (5-Year Impact Factor: 4.193). https://doi.org/10.1016/j.compstruct.2017.06.069. vii  Subsequently, I generalized the analytical model introduced in Chapter 4 in a manner that could be used for the prediction of wrinkling of woven fabrics under another possible shear deformation mode (i.e., pure multi-scale shear deformation mode). There was a source of controversy in the literature for quite a long time that why the onset of wrinkling in an identical woven fabric but under bias extension test and picture frame test is asynchronous. Our analytical study unveiled the rationale behind this experimental observation. This part of the project formed Chapter 5 of my thesis and led to the following publications: - A. Hosseini, M. Haghi Kashani, F. Sassani, A.S. Milani, F.K. Ko. “Identifying Distinct Shear Wrinkling Behavior of Woven Composite Preforms under Bias Extension and Picture Frame Tests”. Composite Structures, 2018 (5-Year Impact Factor: 4.193). https://doi.org/10.1016/j.compstruct.2017.11.033. - M. Haghi Kashani, A. Hosseini, F. Sassani, F.K. Ko, A.S. Milani. “Bias extension test or picture frame test? Toward a best-practice for characterizing the true shear behavior of woven fabrics for forming simulations”.  The 21st Annual Pacific Centre for Advanced Materials and Microstructures (PCAMM 2016) Conference, December 2016, Vancouver, Canada. Next phase of the project was to implement the shear model of woven fabrics and the analytical model for prediction of shear wrinkling into ABAQUS FE package. To do so, we collaborated with Mr. Masoud Hejazi who was an M.ASc. student in our laboratory. He helped us write a User Material Subroutine (UMAT) in ABAQUS. In this phase, Dr. Kashani provided the constitutive equations of the woven fabric, I presented the analytical wrinkling indicator and the algorithm for detection of shear wrinkling and Mr. Hejazi implemented them into ABAQUS. The written UMAT code was submitted to the 2017 American Society of Composites conference at Purdue University to participate in the material model code viii  competition. We were selected as one of the scholarship winners. The outcome of this part of the project is as follows: - M. Haghi Kashani, A. Hosseini, F. Sassani, F.K. Ko, A.S. Milani. “A coupled non-orthogonal hypoelastic model for woven fabrics”. Under Submission, 2018. - M. Haghi Kashani, M. Hejazi, A. Hosseini, F. Sassani, F.K. Ko, A.S. Milani. “Toward enhanced forming simulation of woven fabrics using a coupled non-orthogonal hypoelastic constitutive model, integrated with a new wrinkling onset criterion”. The 32nd American Society of Composites (ASC) Technical Conference, October 2017, Purdue University, USA. Finally, a parametric study was conducted in which the effect of fabric and process parameters on wrinkling phenomenon was investigated. I was the primary contributor in this part of the project and our collaboration with Prof. Milani’s laboratory is highly recognized. The results of this part of the project are under preparation for publication in a journal: - A. Hosseini, M. Haghi Kashani, F. Sassani, A.S. Milani, F.K. Ko. “A parametric study of forming induced wrinkles in woven composite preforms”. Under Preparation, 2018.   Chapter 6 of my thesis is composed of “the implementation of wrinkling model into ABAQUS” and “parametric study of the wrinkling phenomenon”.   ix  Table of Contents Abstract ..................................................................................................................................... ii Lay Summary ........................................................................................................................... iii Preface...................................................................................................................................... iv Table of Contents ..................................................................................................................... ix List of Tables .......................................................................................................................... xii List of Figures ........................................................................................................................ xiii List of Symbols ..................................................................................................................... xxii List of Abbreviations ........................................................................................................... xxvi Acknowledgements ............................................................................................................. xxvii Chapter One: Introduction ........................................................................................................ 1 1.1 Nature of Forming Induced Wrinkles ........................................................................ 3 1.1.1 A General Overview of Wrinkling ..................................................................... 3 1.2 Forming Induced Wrinkles in Woven Composite Preforms ...................................... 5 1.3 Summary .................................................................................................................. 10 Chapter Two: Background Review and Research Objectives ................................................ 11 2.1 Overview .................................................................................................................. 11 2.2 Methods for the Analysis of Wrinkling.................................................................... 11 2.2.1 Analytical Methods ........................................................................................... 11 2.2.2 Finite Element Analysis .................................................................................... 18 2.3 Research Objectives and Thesis Framework ........................................................... 26 Chapter Three: Meso-Scale Mechanism of Shear Deformation of Woven Fabrics ............... 29 3.1 Overview .................................................................................................................. 29 3.2 Shear Behavior of Woven Fabrics ........................................................................... 30 3.2.1 Shear Deformation of Woven Fabrics under Trellising Mode ......................... 30 x  3.2.2 Characterization of Shear Behavior of Woven Fabrics .................................... 34 3.2.3 Proposing a New Shear Deformation Mode for Woven Fabrics ...................... 44 3.2.4 Modifying Prodromou’s Equation for Prediction of Locking Angle of Woven Fabrics under Trellising Mode ........................................................................................ 48 3.3 Summary .................................................................................................................. 50 Chapter Four: A Mesoscopic Analytical Model to Predict the Onset of Wrinkling in Plain- ation......................................................................................................................................... 52 4.1 Overview .................................................................................................................. 52 4.2 Mechanism of Wrinkle Formation under Trellising Shear Mode ............................ 53 4.3 Prediction of Fabric Wrinkling Based on an Instability Analysis ............................ 56 4.3.1 Analogy-Based Equivalent Structure ................................................................ 56 4.3.2 Characterization of Stiffness Elements of the PWCPE Structure ..................... 59 4.3.3 Determination of Critical Compressive Force (Pcr) and Shear Angle at the Onset of Wrinkling .................................................................................................................... 62 4.4 Experimental Evaluation .......................................................................................... 66 4.4.1 Geometric Characterization .............................................................................. 66 4.4.2 Meso-Mechanical Characterization of the Fabric ............................................. 67 4.4.3 Bias Extension Test for Validating the Analytical Wrinkling Model............... 72 4.5 Summary .................................................................................................................. 74 Chapter Five: Prediction of Wrinkling Initiation in PWCPs under Pure Multi-Scale Shear Deformation Mode  ................................................................................................................. 76 5.1 Overview .................................................................................................................. 76 5.2 Mechanism of Wrinkle Formation under Pure Multi-Scale Shear Deformation Mode   .................................................................................................................................. 77 5.3 Modified Analytical Model for the Prediction of Wrinkling under Pure Multi-Scale Shear Deformation Mode .................................................................................................... 79 Woven Preforms under Trellising Shear Deformation Onset of Wrinkling Deformation Mode Shear Deformation Mode Fabrics under Trellising Mode xi  5.4 Why Are the Onsets of Wrinkling in Bias Extension Test and Picture Frame Test Asynchronous? .................................................................................................................... 83 5.4.1 Common Understanding in the Literature ........................................................ 83 5.4.2 New Understanding .......................................................................................... 84 5.5 Experimental Validation .......................................................................................... 87 5.6 Summary .................................................................................................................. 92 Chapter Six: Potential Applications of the Analytical Model of Wrinkling Prediction  ........ 94 6.1 Overview .................................................................................................................. 94 6.2 Implementation of the Analytical Model in ABAQUS FE Package ........................ 94 6.3 Parametric Study of Wrinkling Initiation ............................................................... 100 6.3.1 Fabric Parameters............................................................................................ 102 6.3.2 Process Parameters.......................................................................................... 108 6.4 Summary ................................................................................................................ 112 Chapter Seven: Conclusions and Future Works  .................................................................. 113 References ............................................................................................................................. 119 Appendix: A Manual for ABAQUS User Subroutine UMAT_Woven Fabric ..................... 126           Asynchronous? xii  List of Tables Table 2-1) A summary of the analytical studies available in the literatue for wrinkling prediction. ............................................................................................................................... 14 Table 2-2)   A summary of efforts to model wrinkling in finite element software. ................ 27 Table 4-1)     Summary of the measured geometric and mechanical characteristics of the PWCP.................................................................................................................................................. 71                  prediction PWCP. xiii  List of Figures Classification of manufacturing processes (the shaded boxes highlight the path of interest for this study). .................................................................................. 1 Sheet stamp forming process in different materials. ......................................... 2 Some examples of sheet forming induced defects in preforms (a) wrinkling in an engineering fabric [6], (b) in-plane waviness in a woven fabric [7], (c) yarn sliding in an NCF [8], (d) slippage of a set of crossovers in a woven fabric [9]............................................................................................................................ 3 A sheet metal under quasi-static compressive load (a) slightly below buckling load, (b) slightly above buckling load............................................................... 4 Schematic of a sheet metal during hemisphere stamping process (a) front view, (b) top view. ...................................................................................................... 5 Mesoscopic deformation pattern of a woven fabric during the forming process............................................................................................................................ 7 (a) Formed specimen, (b) details of wrinkling [25]. ......................................... 8 A plain-weave fabric element (a) before shear, (b) during shear and prior to locking, (c) at the lock-up, (d) in post-locking [23]. ......................................... 8 Mechanism of wrinkle formation in a plain-weave element due to large in-plane shear deformation (T and Ms denote tension along yarn and in-plane shear coupling, respectively). ..................................................................................... 9 Wrinkling due to large shear deformation in (a) plain-weave fabric, (b) twill-weave fabric, (c) 5H satin fabric [32]. .............................................................. 9 Figure 1-1) Figure 1-2) Figure 1-3) Figure 1-4) Figure 1-5) Figure 1-6) Figure 1-7) Figure 1-8) Figure 1-9) Figure 1-10) xiv  Flowchart of detecting the onset of wrinkling in FEM simulation using the proposed wrinkling indicator/predictor [15]. .................................................. 12 Load-deformation behavior for an in-plane loaded plate. .............................. 13 (a) A plain woven fabric element, (b) geometry of the fabric element before and after shear (reproduced from [39] for quality). ............................................... 16 Effective factors of an FE model in the simulation of wrinkling. .................. 19 Discrete model of a plain-woven fabric [45]. ................................................. 20 A four-node element composed of woven unit cells [51]. .............................. 21 Draping of a woven fabric with considering (a) only tension term, (b) tension and trellising shear terms, (c) tension, trellising shear and bending terms [10].......................................................................................................................... 22 Methods for modeling wrinkles in sheet metal forming processes................. 23 A simply supported, in-plane loaded, thin elastic plate in (a) ABAQUS Standard, (b) ABAQUS Explicit..................................................................... 26 A typical force-shear angle response of woven fabrics at macro-level, showing four phases: 1- shear with static friction, 2- shear with dynamic friction, 3- locking, and 4- wrinkling. ............................................................................... 30 (a) The motion field of a representative fabric element under trellising mode, (b) velocity vectors during deformation. The centers of rotation are shown as CRi. l is the initial distance between the central axes of the adjacent parallel yarns, and S represents the length of yarn in the fabric element. ................... 31 In-detail meso-level deformation of a cross-over (a) before deformation; (b) after trellising shear deformation. ................................................................... 32 Figure 2-1) Figure 2-2) Figure 2-3) Figure 2-4) Figure 2-5) Figure 2-6) Figure 2-7) Figure 2-8) Figure 2-9) Figure 3-1) Figure 3-2) Figure 3-3) xv  A plain-woven fabric element (a) before shear, (b) during shear and prior to locking, (c) at the lock-up. .............................................................................. 33 Two different test set-ups commonly used for the shear characterization of woven fabrics: (a) bias extension test, (b) picture frame test [29]; in both tests the yarns are initially oriented at ±45̊ with respected to the vertical loading direction. ......................................................................................................... 34 A sample (a) at the beginning of BET, (b) during BET.................................. 35 A typical picture frame test set up. ................................................................. 36 In-detail meso-level deformation of a cross-over (a) before deformation; (b) after shear deformation in picture frame test. γ, w0 and w and are the shear angle of the fabric, the initial width and the instantaneous width of the yarns, respectively. .................................................................................................... 37 (a) Shear force – shear angle diagram, (b) fabric element under trellising shear mode (only yarn rotation causes the shear of the fabric element), (c) fabric element under pure shear mode (taking the intra-yarn shear into account) [45].......................................................................................................................... 38 Comparison between intra-yarn shear (i.e., in this figure it is named by local shear angle) and the applied frame shear angle [66]. RR1, RR2, and RR3 are three different fabrics, the stiffest of which is RR3. A and B represent the shear angle within yarns spacing and the yarn itself. ............................................... 39 Yarn width change during shear deformation of the fabric in PFT [67]. ....... 40 Schematic of a fabric element of a typical plain-woven fabric modeled in TexGen (a) top view, (b) 3D presentation of a single yarn in the fabric element Figure 3-4) Figure 3-5) Figure 3-6) Figure 3-7) Figure 3-8) Figure 3-9) Figure 3-10) Figure 3-11) Figure 3-12) xvi  subject to in-plane shear force, (c) a simplified structural beam representation. Notice that the in-plane shear force is along with the width of yarns, which is generally not substantially less than the length of a yarn between two crossovers. ....................................................................................................... 41 Coordinate systems to demonstrate the undeformed and deformed shape of a fabric sample subject to an ideal picture frame test. ....................................... 44 (a) A representative fabric element under picture frame test, (b) deformation of a yarn segment during the picture frame test, showing intra-yarn shear and hence a decrease in the instantaneous yarn width, (c) deformation of the representative fabric element. ......................................................................... 45 Trajectory of the spacing between the adjacent parallel yarns over pure multi-scale shear deformation mode, which is drawn using equation (3-6). ............ 47 The effect of interlacing yarn on locking angle by considering 3D schematic of the fabric. ........................................................................................................ 48 A fabric element at the beginning of trellising shear deformation mode. ...... 49 Predictions of the onset of locking in a fabric element under trellising mode using equation (3-9) and equation (2-2) [39]. The width and thickness of the yarns and the spacing between the adjacent parallel yarns are assumed 1.4 mm, 0.355 mm and 0.65 mm, respectively. ............................................................. 49 A comparison of the predictions of locking angle by relying on trellising shear mode and pure multi-scale shear mode. The width and thickness of the yarns and the spacing between the adjacent parallel yarns are assumed 1.4 mm, 0.355 mm and 0.65 mm, respectively. ....................................................................... 50 Figure 3-13) Figure 3-14) Figure 3-15) Figure 3-16) Figure 3-17) Figure 3-18) Figure 3-19) xvii  Wrinkles formed due to (a) in-plane compression, (b) out-of-plane bending, and (c) in-plane shear............................................................................................. 53 A PWCP in the course of large in-plane shear deformation. .......................... 54 Compressive forces induced between adjacent yarns during trellising shear deformation of PWCP; (a) simulation, (b) experimental investigation of the yarn strain field in large shear deformation [3]. ............................................. 55 Mechanism of wrinkle formation in a PWCP element due to large trellising shear deformation (T and Ms denote the tension along yarn and in-plane shear coupling, respectively). ................................................................................... 56 A four-link equivalent structure mimicking the behavior of a representative PWCP element in the post-locking stage. ....................................................... 57 Equivalent structure (a) before equilibrium deviation, (b) after equilibrium deviation. ......................................................................................................... 58 Equivalent torsion spring Kt, mimicking the bending effect of interlacing yarn Y-3 on the out of plane deformation of Y-2, and vice versa. ......................... 59 Modeling of yarn as a bending beam (made of a linear, transversely isotropic, homogeneous material with an effective flexural rigidity in the longitudinal direction Qb), and an equivalent torsion spring coefficient, Kt. ...................... 60 The PWCPE structure in an imaginary deviated mode; note that the P forces remain in the original directions. .................................................................... 64 The geometric characteristics of the tested PWCP. Mean value of yarn width along with its standard deviation are, in turn, 1.4 mm and 0.04 mm. Mean value of spacing between yarns along with its standard deviation are, in turn, 0.65 mm Figure 4-1) Figure 4-2) Figure 4-3) Figure 4-4) Figure 4-5) Figure 4-6) Figure 4-7) Figure 4-8) Figure 4-9) Figure 4-10) xviii  and 0.05 mm. Mean value of yarn thickness along with its standard deviation are, in turn, 0.355 mm and 0.04 mm. ............................................................... 67 The experimental set-up used for measuring effective bending rigidity of the yarn. ................................................................................................................ 68 Test set-up used to measure the lateral stiffness of the yarn. .......................... 69 The force-displacement curve for a single yarn under compression in the width direction (Force-displacement curve without any yarn in the fixture was first drawn to determine the resistance of the fixture; then this resistance was subtracted from force-displacement curve of the yarn under lateral compression and the net value is shown here). .................................................................... 70 The PWCP (a) before shear deformation, (b) at the onset of locking (𝛾 =28.7±0.8 ̊), (c) at the onset of wrinkling (𝛾 =38±0.5 ̊). ........................ 73 Experimental measurements of locking and wrinkling angles versus the predictions of the proposed analytical model. ................................................ 73 Reduction rate of the instantaneous width of the fabric element (L) and the total instantaneous width of two parallel adjacent yarns (2w) under pure multi-scale shear mode  ..................................................................................................... 78 The contact forces due to large pure multi-scale shear deformation exerted on (a) the weft yarns, (b) the warp yarns (the black line between the adjacent parallel yarns represents the interlacing yarn from top view)......................... 78 (a) A fabric element deformed and locked under pure multi-scale shear deformation, (b) an equivalent representative structure mimicking the behavior of the fabric element in the post-locking stage. .............................................. 80 Figure 4-11) Figure 4-12) Figure 4-13) Figure 4-14) Figure 4-15) Figure 5-1) Figure 5-2) Figure 5-3) xix  A comparison of the predictions of wrinkling angle (i.e., the shear angle at which wrinkling initiates) by relying on trellising shear mode and pure multi-scale shear mode. The width and thickness of the yarns and the spacing between the adjacent parallel yarns are assumed 1.4 mm, 0.355 mm and 0.65 mm, respectively. .................................................................................................... 82 Analytical representation of the yarn spacing of the fabric element under trellising and pure multi-scale shear deformation modes. .............................. 85 Analytical representation of the yarn width of the fabric element under trellising and pure multi-scale shear deformation mode. ............................................... 86 The picture frame test set-up used in the experiments. ................................... 88 (a) Using needles in the PFT boundary conditions, (b) a folded sample in the grips to increase the interaction between the fabric ends and the needles. ..... 90 (a) Wrinkling initiation in the PFT arm regions due to the presence of transverse yarns, (b) removing most of the transverse yarns, (c) the shear deformation at 64o with no wrinkles. ...................................................................................... 90 Experimental measurements of locking and wrinkling angles versus the predictions proposed by the analytical model. ................................................ 92 Validation of the UMAT code implemented in ABAQUS, (a) accurate prediction of longitudinal stress under sequential tensile test with 2% transverse pre-tension, (b) precise anticipation of transverse stress under sequential tensile test with 2% transverse pre-tension. These graphs highlight that how significant is the effect of coupling on the fabric behavior.  ............................................ 96 Figure 5-4) Figure 5-5) Figure 5-6) Figure 5-7) Figure 5-8) Figure 5-9) Figure 5-10) Figure 6-1) xx  The flowchart based on which wrinkling initiation in the finite element code is detected. .......................................................................................................... 98 (a) FE simulation results of the woven fabric equipped with the wrinkling detection flowchart (i.e., the red spot shows the region in where wrinkling is about to occur), (b) a real sample under bias extension test. ........................ 100 FE simulations of woven fabric forming process, (a) [62], (b) [27]. ............ 100 Parameters affecting wrinkling of the fabrics. .............................................. 101 Qualitative analysis of the effect of yarn bending rigidity on the critical shear angle of the plain-woven fabric under (a) trellising mode, (b) pure multi-scale shear mode. ................................................................................................... 103 Revisited figure from Chapter 3. .................................................................. 104 Qualitative analysis of the effect of yarn lateral stiffness on critical shear angle of the plain-woven fabric under (a) trellising mode, (b) pure multi-scale shear mode. ............................................................................................................. 105 Qualitative analysis of the effect of spacing between yarns on critical shear angle of the plain-woven fabric under (a) trellising mode, (b) pure multi-scale shear mode. ................................................................................................... 106 Draped woven fabric over a hemispherical punch, (a) a loose plain-woven fabric, (b) a tight plain-woven fabric [32]. ................................................... 107 Qualitative analysis of the effect of thickness of interlacing yarn on locking angle of the plain-woven fabric under trellising mode. ................................ 108 The schematic of a draping process. ............................................................. 109 Figure 6-2) Figure 6-3) Figure 6-4) Figure 6-5) Figure 6-6) Figure 6-7) Figure 6-8) Figure 6-9) Figure 6-10) Figure 6-11) Figure 6-12) xxi  Qualitative analysis of the effect of tension in yarns on the critical shear angle of the plain-woven fabric under, (a) trellising mode, (b) pure multi-scale shear mode. ............................................................................................................. 110 The effect of BHF on wrinkling of fabrics, (a) BHF = 0, (b) BHF = 5000 N [81]........................................................................................................................ 111 Effect of punch displacement and magnitude of BHF on the amplitude of wrinkles [6]. .................................................................................................. 111                             Figure 6-13) Figure 6-14) Figure 6-15) xxii  List of Symbols  English Symbols  A   a lumped parameter A' cross-sectional area of a yarn A  the length of the plate b the width of the plate C the right Cauchy–Green tensor 𝑐𝑔 a correction factor representing the impact of gaps between filaments within a dry yarn dWshear the increment of the work done on the fabric element dV the volume of fabric element D the bending stiffness of the plate Dʹ the diameter of fiber cross-section E Young’s modulus of the yarn E strain tensor Fshear the shear load per unit length G shear modulus of a yarn I bending moment of inertia of a yarn J the polar moment of inertia of a yarn l initial distance between the central axes of the adjacent parallel yarns xxiii  L the instantaneous length of the fabric element L0 the initial length of the fabric element L (, Δ) total potential energy of the PWCPE structure 𝑀𝑦 the resultant bending moment n the number of fibers in a yarn P the compressive force between the adjacent parallel yarns ?̅?  the lateral compressive force between yarns in Zhu’s model Pcr the critical compressive force between the adjacent parallel yarn under trellising mode (Pcr)pure multi-scale shear mode or (Pcr)PMSM the critical compressive force between the adjacent parallel yarn under pure multiscale shear mode 𝑄𝑏 the effective flexural rigidity of yarn in its longitudinal direction S the length of yarn in the representative PWCP element SP spacing between yarns T tension along the yarn t yarn thickness 𝑇𝑥 the resultant torsional torque 𝑉𝑧 shear force at each point of the beam w0 the initial yarn width w the instantaneous width of yarn w0 initial width of yarn Wint(η) the internal virtual work xxiv  Wext(η) the external virtual work Wacc(η) the virtual work by acceleration 𝑊𝑖𝑛𝑡𝑇 (η) energy dissipated for tensioning the yarns 𝑊𝑖𝑛𝑡𝑆 (η) energy dissipated for trellising shear between the yarns 𝑊𝑖𝑛𝑡𝑏 (η) energy dissipated for bending the yarns (Wshear)locking-to-wrinkling the area under force-displacement of the BET X the undeformed vector of a material point within the fabric subjected to PFT 𝑥 the deformed vector of a material point within the fabric subjected to PFT   Greek Symbols  𝛼  cross-section shape coefficient of the equivalent beam γ  the inter-yarn shear angle 𝛾Locking  locking angle 𝛾Wrinkling  wrinkling  angle Δ compaction of the yarn due to lateral compressive force or the compaction of the slide springs Δcr  the critical yarn compaction in width direction just before wrinkling under trellising mode (∆𝑐𝑟)PMSM  the critical compaction of yarns in the width direction just before the wrinkling under pure multi-scale shear mode xxv  𝜀1  the longitudinal strains along the warp 𝜀2  the longitudinal strains along the weft yarn 𝜂 a virtual displacement field   the angle of rods with respect to the plane of initial mode                   xxvi  List of Abbreviations BET bias extension test BHF blank holder force CR centers of rotation FE finite element FEA finite element analysis FEM finite element modeling FRP fiber-reinforced plastic NCF non-crimp fabric PJT pin-joint theory PFT picture frame test PWCP plain-woven composite preform PWCPE PWCP Equivalent structure UMAT user material subroutine           xxvii  Acknowledgements Firstly, I would like to express my special appreciation and thanks to my co-advisors Prof. Farrokh Sassani and Prof. Frank K. Ko for the support of my Ph.D. study and related research. You have been truly a tremendous mentor for me with immense knowledge and endless patience. I would like to thank you for encouraging my research and for allowing me to grow as a research scientist. In addition to my co-advisors, I would like to express my sincere gratitude to the rest of my thesis committee: Prof. Steve Feng and Prof. Srikantha Phani for their insightful comments. My special thanks also go to Prof. Abbas Milani who provided me with an opportunity to have collaboration with his team members and who gave me access to the laboratory and research facilities. Without his precious support, it would not be possible to conduct the experimental part of the research. Furthermore, the support of colleagues and stimulating discussions at the Composites Research Network (CRN) are greatly valued. I am extremely thankful to Prof. Michael Sutcliffe who hosted me at the University of Cambridge and allowed me to do research into the analysis of wrinkling in non-crimp fabrics. This precious opportunity gave me an insight into forming behavior of NCFs and paved the way for better analysis of this type of composite reinforcement. Moreover, I thank my fellow labmates – Dr. Behnam Razavi, Dr. Atefeh Einafshar, Mr. Morteza Taiebat, Mr. Sina Nezafatkhah, Mr. Mohammad Saryazdi, Ms. Sara Azizkhani and Ms. Hamide Torabi - for all the fun we have had in the last four years. In particular, I am grateful to Dr. Masoud Haghi Kashani and Mr. Masoud Hejazi for the stimulating discussions and all the beautiful time we have had together. xxviii  Last but not the least, I would like to thank immensely: - my parents who devoted their lives to me. - my siblings who supported me spiritually. - my wife who literally motivated me whenever I faced a challenge. She has been always inspiring in my life and without her patience and support; taking even one step was not possible. - my adorable son who made my life colorful and gave me the impetus to be stronger in encountering hardships and challenges.              xxix     To My Beloved Parents– Ashraf and Mohammad   To My Lovely Wife– Atena   To My Adorable Son – Adrian     1  1 Chapter One Introduction The word manufacturing first appeared in English in 1567 and is derived from the Latin manu factus, meaning “made by hand” [1]. This term, however, has evolved over time and now it means the transformation of raw materials into useful products through the use of least expensive methods and not necessarily made by hand [2]. Manufacturing processes have been classified in a wide variety of ways, one of the most comprehensive of which is presented by El Wakil [2]. Figure 1-1 shows a summary of the classification of the manufacturing processes. Forming, as the manufacturing process of interest in this study, is defined as changing the shape of the raw stock without adding material to it or taking material away from it. Forming processes are performed through two distinct methods: bulk forming and sheet forming. In sheet forming processes, stamping or deep drawing is one of the most practical techniques to shape the raw material [3][4][5].  Figure 1-1) Classification of manufacturing processes (the shaded boxes highlight the path of interest for this study). Manufacturing ProcessesFormingBulk Forming Sheet FormingIncremental Forming SpinningStamping/Deep DrawingMisc.Joining Molding Material Removing2  The raw material for sheet stamping process could be a sheet metal, a fiber-reinforced plastic (FRP) sheet, a fiber-metal laminate, etc. In the case of sheet stamping of FRPs, there exist two commonly used techniques, namely preforming along with resin impregnation and stamping of pre-impregnated sheets (i.e., where the reinforcement is usually pre-impregnated with a thermoset resin; however, pre-consolidated thermoplastic sheets are sometimes classified in this group). Figure 1-2 presents some of the material system used for sheet stamping process and highlights the path of interest of this study with shading. As the figure shows, this thesis is focused on preforming of a class of 2D fabrics called woven fabrics. This process is also known as draping or sheet stamp forming (i.e., sometimes its short form – sheet forming – is used) of engineering fabrics. Despite the fact that there are a large number of various engineering fabrics used in the industry, we selected woven fabrics because of their substantially extensive application in different industrial sectors, such as aerospace, automotive, marine and energy (e.g., wind turbine).  Figure 1-2) Sheet stamp forming process in different materials. Sheet Stamp FormingSheet Metals FRP CompositesPreforms1DUnidirectional2DWoven Fabrics Non-Crimp Fabrics3DInterlock WeavePre-pregs  Fiber – Metal Laminate Composites 3  The previous studies have revealed that several forms of defects would appear during the course of preforming process, which could affect both the appearance and the performance of the final product [6][7][8]. Some examples of these defects are shown in Figure 1-3. Among the forming induced defects, wrinkling has been identified as one of the most critical defects [9][10]. Accordingly, the current study focuses explicitly on the wrinkling behavior of woven fabrics.  Figure 1-3) Some examples of sheet forming induced defects in preforms (a) wrinkling in an engineering fabric [6], (b) in-plane waviness in a woven fabric [7], (c) yarn sliding in an NCF [8], (d) slippage of a set of crossovers in a woven fabric [9].  1.1 Nature of Forming Induced Wrinkles 1.1.1 A General Overview of Wrinkling Compared to sheet metal forming, the field of sheet forming of composite preforms and the analysis of wrinkling is in its infancy dating back only to the late-1980s [11][12]. While in sheet metal forming the wrinkling phenomenon has been extensively studied for a relatively long time, with much research conducted in the 1950s [13]. In sheet metal forming processes, wrinkling is treated as a local buckling problem, which appears in some regions of the formed sheet metal [14][15]. Buckling in classical mechanics is defined as the deformation state at 4  which a deviation from in-plane deformation mode to out-of-plane deformation mode arises due to an instability [16][17]. The instability associated with buckling is developed owing to excessive compressive force within the material. Figure 1-4 illustrates the schematic of one mode of plate buckling under in-plane compression. As can be seen in the figure, the sheet metal suddenly changes its deformation mode from in-plane to out-of-plane when quasi-static compressive force P reaches a critical value which is referred to as buckling load. The magnitude of this load depends on a number of parameters such as material properties (i.e., material moduli and Poisson’s ratios), plate geometry (i.e., width, length and thickness of the plate), boundary conditions, etc. [16].   Figure 1-4) A sheet metal under quasi-static compressive load (a) slightly below buckling load, (b) slightly above buckling load.  For sheet metals, compression within the material could be created either by applying an external compressive load or by imposing displacements which generate compressive strains. In any case, the net result is local buckling or wrinkling. The compression induced during sheet metal forming processes predominantly belongs to the latter case of imposed displacements. Figure 1-5 illustrates this point, showing how material initially located at strip I is deformed to location strip II by the drawing process. The reduction in the strip circumference causes in-plane compression within the material [18]. If the compression at a given position within the 5  material exceeds the buckling load, which depends on several factors [16][19], localized buckling or wrinkling in that region of sheet metal will occur [14][20][21][22].   Figure 1-5) Schematic of a sheet metal during hemisphere stamping process (a) front view, (b) top view.  1.2 Forming Induced Wrinkles in Woven Composite Preforms In spite of sheet metals that constriction in circumference during forming process develops large compressive stress within the material, the story is quite different for woven composite preforms in which shear stiffness is relatively low in the early stage of shear deformation. The inherently low shear stiffness of the fabrics, which roots in their meso-structure, makes it possible to reduce the length of the strip in some regions of the fabric without causing high in-plane compression. This issue significantly impacts wrinkling behavior of the fabrics. Accordingly, it is of great importance to delve into the mesoscale deformation mechanisms of the fabrics to better understand the nature of wrinkling in the woven reinforcements. As stated earlier, the behavior of woven fabrics in sheet forming processes is not identical with that of sheet metals due to their initially low shear stiffness and quasi-inextensibility of the yarns. These two important features of woven fabrics lead to the issue that intra-ply shear 6  becomes the main deformation mode of the fabrics for conforming a shape [23]. In addition to intra-ply shear, however, there exist some other mechanisms which contribute to the deformation of the fabric. The occurrence of each mechanism at any local point of the deforming fabric significantly depends on the geometry of the punch, boundary conditions, etc. This explanation would be more comprehensible by studying an example in which a woven fabric is draped into a hemispherical shape. Hemispherical stamping is selected not only on account of its simplicity but also its capability to create almost entire deformation mechanisms occurring in the forming of complex geometries [24]. If we consider two virtual strips on a woven fabric as shown in Figure 1-6, by the punch penetration, the strip I is forced into strip II. In woven fabrics, the length constriction is a function of mesoscale deformation mechanisms of the fabric. In Region A, for example, the intra-ply shear brings about constriction in length of the strip (i.e., intra-ply shear occurs in this region of the fabric since the force exerted to the fabric from the punch is in the bias direction of the warp and weft yarns), while the reduction in length of the strip in Region B can be accommodated only by compression (i.e., since the punch force is parallel to one group of yarns and perpendicular to the other group of yarns, shear does not occur in this region). The length of the strip in Region C could be reduced by a combination of these deformation mechanisms. If the compression induced in Region B or C becomes larger than a critical value, wrinkling occurs. This type of wrinkles, from the formation mechanism point of view, is analogous to those appearing in the sheet metal forming processes. Only a few papers could be found in the past research works that reported this class of wrinkling in the forming of woven fabrics [25][10][26]. Among those, Skordos et al. elaborately looked into this phenomenon using experimental investigations and finite element 7  simulations [25]. Figure 1-7 shows an example of forming induced wrinkles in a part of a woven preform with mesoscale architecture analogous to Region B in Figure 1-6.   Figure 1-6) Mesoscopic deformation pattern of a woven fabric during the forming process.  Besides this type of wrinkling, there exists another mechanism for the formation of wrinkles in woven fabrics. If we return to the example presented in Figure 1-6, fabric elements in Region A take the shape of the punch – reduce the length of the strip segment – through intra-ply shear deformation. In this region, the shear angle of fabric elements increases by the punch movement. In the initial stage of shear deformation, the shear resistance of the fabric is low due to low frictional forces between overlapping yarns [27]. Afterward, the shear resistance of the fabric increases more rapidly when inter-yarn gaps vanish – locking or 8  jamming initiates – and adjacent parallel yarns begin to press each other transversely (Figure 1-8 (a–c)).  Figure 1-7) (a) Formed specimen, (b) details of wrinkling [25].   Figure 1-8) A plain-weave fabric element (a) before shear, (b) during shear and prior to locking, (c) at the lock-up, (d) in post-locking [23].  Further intra-ply shear of the fabric in the post-locking stage is dominantly by yarns lateral compaction (Figure 1-8 (d)) [23][10][28][29]. By growing the shear deformation of woven fabrics, the pressure between adjacent parallel yarns – in-plane compression – rises, triggering the formation of another type of wrinkles, which is referred to as shear wrinkles (Figure 1-9) [30][31]. 9   Figure 1-9) Mechanism of wrinkle formation in a plain-weave element due to large in-plane shear deformation (T and Ms denote tension along yarn and in-plane shear coupling, respectively).  This type of wrinkling is analogous to the wrinkles formed in bias extension test and picture frame test (a detailed review of the mesoscale deformation mechanisms of these two tests is presented in Chapter 3). Figure 1-10 illustrates some examples of the wrinkles formed due to large intra-ply shear deformation. Hereafter, in order to distinguish different types of wrinkles, terms “non-shear wrinkling” and “shear wrinkling” are used for the former and latter types of the wrinkling.  Figure 1-10) Wrinkling due to large shear deformation in (a) plain-weave fabric, (b) twill-weave fabric, (c) 5H satin fabric [32]. 10   1.3 Summary To recapitulate, forming induced wrinkles in woven fabrics could be classified into two distinct categories, namely non-shear and shear wrinkles. Although the nature of both types of wrinkles is identical (i.e., both wrinkles are developed due to excessive in-plane compression), the mechanisms contributing to initiate this phenomenon are quite different. That is to say, the in-plane compression induced in the former wrinkling case is due to the longitudinal contraction (of the length) of the fabric segment, whereas in the latter case, this compression is initiated because of the large intra-ply shear deformation of the fabric (i.e., compression is as a result of locking).            11  2 Chapter Two Background Review and Research Objectives 2.1 Overview This chapter presents a detailed discussion of the methods for the analysis of wrinkling. First, analytical methods for the prediction of wrinkling are reviewed. Subsequently, finite element simulations are introduced as a useful tool for the study of wrinkling phenomenon. The integrated review of analytical methods along with FE simulations answers an important question that how having an accurate analytical criterion for wrinkling prediction could improve the accuracy of the finite element simulations. Moreover, the gaps in the literature and the needs of the industry are highlighted, based on which the research objective is defined. 2.2 Methods for the Analysis of Wrinkling 2.2.1 Analytical Methods 2.2.1.1 The Importance of Analytical Methods Analytical methods have always been a powerful tool to either establish criteria for wrinkling prediction or identify the contributing parameters in the wrinkling phenomenon. These methods have been extensively employed in the field of sheet metal forming. Considering the fact that the literature of analytical methods in sheet metal forming is far more established than that of the preforming field, it would be advantageous to briefly review some of the studies available in the field of sheet metal forming. This would lead to acquiring a better insight on how to effectively use analytical methods in studying the wrinkling of composite 12  preforms. In this section, therefore, a couple of potential applications of analytical methods in sheet metal forming are discussed. - Application one (establishing the criteria for wrinkling prediction): The majority of analytical criteria for the prediction of wrinkling of sheet metals have been proposed by relying on the fact that wrinkling is a local buckling problem [14][15][33][34]. As an important application, these criteria are usually linked to finite element packages to streamline the capability of the software for the prediction of wrinkling initiation [14][15][35]. Wang and Cao [15], as an example, developed an analytical criterion to predict wrinkling in sheet metals and then, incorporated it in ABAQUS FE package as a wrinkling detector. Using the flowchart shown in Figure 2-1, they could identify the regions of the sheet metal in where wrinkling would occur during sheet forming process.  Figure 2-1) Flowchart of detecting the onset of wrinkling in FE simulation using the proposed wrinkling indicator/predictor [15]. 13  - Application two (identifying the contributing parameters in the wrinkling phenomenon): Besides the application of analytical methods in predicting the onset of wrinkling, these methods could also be used for identifying the contributing parameters in wrinkling. This statement could be more tangible by studying a simple example of wrinkling in a sheet material. Figure 2-2 illustrates the load-deformation behavior of a simply supported, perfectly flat, isotropic, elastic plate uniformly loaded in compression.                   Figure 2-2) Load-deformation behavior for an in-plane loaded plate.  As illustrated in the load-deformation response of Figure 2-2, initially the load increases linearly with the applied displacement, as the plate simply shortens in the load direction while remaining flat. Then, at a critical load Ncr, there is a bifurcation in the response as the deformation switches from the flat shape to the buckled shape. The critical load per unit width at which the plate buckles – the buckling load – is given by:  𝑁𝑐𝑟 = 𝐷 [𝜋2 + (𝜋𝑎𝑏)2]  (2-1)   14  where a, b and D are the length, the width and the bending stiffness of the plate, respectively [16]. The bending stiffness of the plate depends on the material properties and the thickness of the plate. Hence, the model identifies the relevant material and geometric parameters controlling buckling. Evidently, the analysis can be extended to more sophisticated examples to understand more about the factors affecting buckling. In summary, this section illustrates another potential application of analytical models for identifying the parameters affecting wrinkling phenomenon. 2.2.1.2 Non-Shear Wrinkling Unlike the sheet metals in which the efforts toward developing analytical wrinkling criteria are noticeable, the same trend is not observed in the composite preforms. In the sense that only a handful of papers, as presented in Table 2-1, could be found that analytically dealt with wrinkling of the fabrics. Table 2-1) A summary of the analytical studies available in the literature for wrinkling prediction. Wrinkling Type  Woven Fabric  Non-Shear Wrinkles  Clapp et al. [36] Anandjiwala et al. [37] Clapp et al. [38]  Shear Wrinkles  Prodromou et al. [39] Zhu et al. [40] Rozant et al. [41] Liu et al. [42]    As mentioned in Chapter 1, non-shear wrinkling has not been widely considered in the fabric forming literature, and most of the work in this area relates to the apparel industry 15  [36][37][38]. The mechanism controlling non-shear wrinkling of fabrics is analogous to wrinkling in sheet metal; thus, the analytical methods developed for the analysis of wrinkling of sheet metals can be modified and used for investigating non-shear wrinkling problems [36][37]. 2.2.1.3 Shear Wrinkling In contrast to non-shear wrinkling, analytical modelling of shear wrinkling during preforming of fabrics has attracted significant attention, see for example [39][40][41][42]. In shear wrinkling, shear angles at which locking and wrinkling take place are of great importance, and the analytical studies available in the literature have focused on the methods that could predict these angles. Among the research toward establishing analytical criteria for shear wrinkling prediction in woven fabrics, Prodromou’s work could be considered as one of the very first efforts [39]. As a matter of fact, they were among the pioneers who recognized the importance of a model to determine locking and the onset of wrinkling in the woven fabrics. In their model, the locking was defined as the shear angle at which the distance between the warp and weft yarns becomes zero. To develop the approach, they modified the pin-joint theory (PJT), which was solely based on the kinematics of the fabric structure, and included fabric construction parameters, such as yarn spacing and yarn size, in the analytical calculations. Figure 2-3 (a) illustrates a plain-woven fabric element. According to the figure, locking was defined as a moment at when the orthogonal distance between the center of the fabric element and the mid-line of the yarn (i.e., DE in Figure 2-3 (b)) becomes equal to the half-width of the yarn. That is to say, the distance between the adjacent parallel yarns becomes zero. 16  An equation in accordance with this geometric assumption was obtained based on which locking of the fabric could be determined. Equation (2-2) yields locking angle of a plain-woven fabric [39]: 𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔 = 90° − 𝑠𝑖𝑛−1(𝑤0𝑙)  (2-2)  where 𝛾Locking, l, and w0 are locking angle, initial distance between the central axes of the adjacent parallel yarns, and the width of the yarn, respectively. Another assumption they made is that the onsets of locking and shear wrinkling are coincidental, while the fabric structure undergoes further shear deformation even after locking due to lateral compaction of yarns. This provides more room for shear deformation and the onset of shear wrinkling is markedly different than that of the locking.         Figure 2-3) (a) A plain-woven fabric element, (b) geometry of the fabric element before and after shear (reproduced from [39] for quality).  In addition to the abovementioned research, an energy-based model to predict the shear wrinkle initiation in the plain-woven fabrics has been proposed. Zhu et al. [40] developed an analytical criterion relying on the assumption that wrinkles initiate whenever the energy required for in-plane shear deformation mode becomes greater than that of the out-of-plane Load Direction Load Direction Weft Yarn Warp Yarn 17  deformation mode. The in-plane shear deformation energy in their model consists of three terms: - The energy dissipated by the rotational friction of the yarns. - The energy dissipated by sliding friction between the adjacent parallel yarns. - The energy dissipated by lateral compaction of the yarns. Likewise, they determined the energy required for out-of-plane deformation (i.e., wrinkling) of the fabric by using the following terms: - The energy needed for bending the yarns. - The energy stored due to tension in the yarns. By equating the energy terms for in-plane and out-of-plane deformation modes, they could propose the following analytical equation [43]:  𝑊𝑣 = 𝐸. 𝑐𝑜𝑠𝛾. [4𝐽𝐷ʹ+𝑛𝜋𝐷ʹ2𝑤0(𝜀1 + 𝜀2)16] − 2𝑤0𝑠𝑖𝑛𝛾. ?̅?(𝛾) (2-3)  where E is Young’s modulus of the yarn, γ is the inter-yarn shear angle, Dʹ is the diameter of fiber cross-section, J is the polar moment of inertia of a yarn, n is the number of fibers in a yarn, w0 is the initial yarn width, 𝜀1 and 𝜀2 are the longitudinal strains along the warp and weft yarn directions, respectively, and ?̅? is the lateral compressive force between yarns. The condition under which wrinkling occurs is as follows: {𝑊𝑣 > 0𝑊𝑣 = 0𝑊𝑣 < 0 (2-4)  Wrinkling does not occurs Critical State Wrinkling occurs 18  The model proposed by Zhu could be challenged given the fact that the comparison was made between two terms of energy that are both internal, and their sum always equals the external work done by the shear forces. As a matter of fact, they tried to make use of stable equilibrium approach to determine the onset of deviation in deformation modes [17]. However, they did not take this key issue into account that the equilibrium destabilizing energy, which is caused by in-plane deformation mode, should be compared to the equilibrium restoring energy [16][17]. Based on the latter proven explained definition, another term which is the work done by external in-plane forces is overlooked that would impact the prediction of the onset of wrinkling. Some other papers are also available in the literature [41][42] that could be useful for the analytical investigation of shear wrinkling phenomenon. It is worth mentioning that the analytical criteria reviewed in this section could have the same applications as those of the sheet metals, namely being implemented in FE packages and being used for the parametric study of wrinkling phenomenon. In the sections ahead, these issues will be more discussed. 2.2.2 Finite Element Analysis Another useful tool to study the forming induced wrinkles of the fabrics is the Finite Element Analysis (FEA). In FEA, the forming process is simulated and the behavior of the fabric during the process is analyzed. The success of the simulation to precisely predict the forming induced wrinkles depends on several factors. Figure 2-4 presents some of the most influential factors having an impact on the accuracy of the FEA. 19   Figure 2-4) Effective factors of an FE model in the simulation of wrinkling.  In accordance with the literature, the major efforts in FE analysis of preforming process – in particular, the forming induced wrinkles – have been devoted to the material modeling strategy; and the other factors have been less comprehensively addressed. While the studies reported in the field of sheet metal forming unveil the fact that FE solver coupled with the methods to model wrinkling are of great significance [21][44]. Hence, this section of the thesis aims to review some of the studies focused on the factors introduced in Figure 2-4. The section is organized in a manner that “Section 2.2.2.1” deals with the material modeling strategy and briefly describes discrete, semi-discrete and continuous approaches; and “Section 2.2.2.2” tackles FE solvers and the methods for modeling the wrinkles. 2.2.2.1 Material Modeling Strategy As stated earlier, the research on material modeling of the fabrics is remarkably broad and above of all, there exists a review paper in which the material modeling of the fabrics is comprehensively discussed [24]. Accordingly, revisiting the topic in the thesis would seem FEA FactorsMaterial Modeling StrategyDiscreteSemi-discreteContinuousFE SolverExplicit ImplicitMethods of Modeling Wrinkles20  redundant. However, for its relevance and completeness, we review a few papers on the subject of fabric modeling strategy to set the stage for discussing the other contributing factors. In the past studies, three major strategies were introduced to model woven fabrics, namely discrete, semi-discrete and continuous. - Discrete Approach: In this method, the components of the fabric including yarns and interlacements are all modeled, separately. Figure 2-5 shows an example of the discrete model of a woven fabric.  Figure 2-5) Discrete model of a plain-woven fabric [45].  In view of the extent of the literature, it is seen that almost the majority of discrete models show a reasonable capability of modeling wrinkling as a result of taking the sub-scale mechanisms of the fabrics into account [46][47][48]; however, the overall computational cost and the level of complexity associated with these models are relatively high [49][50]. - Semi-Discrete Approach: The basic idea of this method is to develop a customized macro-scale element, which represents the behavior of a number of unit cells of the fabric. Figure 2-6 illustrates an example of a four-node element with classic bilinear interpolation functions. 21   Figure 2-6) A four-node element composed of woven unit cells [51].   If a virtual displacement field 𝜂 is applied to the presented element, the virtual work theorem relates the internal (Wint(η)), external (Wext(η)), and acceleration (Wacc(η)) virtual works as follows [52]:  𝑊𝑒𝑥𝑡(𝜂) − 𝑊𝑖𝑛𝑡(𝜂) = 𝑊𝑎𝑐𝑐(𝜂) (2-5)  The first and the last terms of the equation are similar to those of the conventional FE methods [52]; however, the internal virtual work in this element is obtained by the following expression:  𝑊𝑖𝑛𝑡(𝜂) = 𝑊𝑖𝑛𝑡𝑇 (η)+ 𝑊𝑖𝑛𝑡𝑆 (η)+ 𝑊𝑖𝑛𝑡𝑏 (η) (2-6)  where 𝑊𝑖𝑛𝑡𝑇 (η), 𝑊𝑖𝑛𝑡𝑆 (η) and 𝑊𝑖𝑛𝑡𝑏 (η) are, in turn, the energies dissipated for tensioning the yarns, trellising shear between the yarns, and bending the yarns [52]. These terms of internal energy are crucially important since the results obtained for wrinkling prediction are sensitive to them [10][53][54]. The point is graphically illustrated in Figure 2-7 in which the behavior of a woven fabric is modeled by including the various internal energy terms. 22   Figure 2-7) Draping of a woven fabric with considering (a) only tension term, (b) tension and trellising shear terms, (c) tension, trellising shear and bending terms [10].  Despite the advantages of the semi-discrete approach such as comparably low computational cost, the coupling effects as an influential factor on wrinkling have not been correctly implemented in this approach. From a mechanistic point of view, coupling means that applying a load in a given direction changes the meso-structure, and subsequently, the mechanical properties of the fabric. Hence, the terms of internal energy should be functions of each other, while in the available studies, they are assumed independent [55].  - Continuous Approach: This type of fabric modeling is the most convenient option for implementation into a finite element package and minimizing the computational cost [49]. In this approach, the fabric  is modeled as a continuous material, while exhibits anisotropic properties while undergoing large shear and bending deformations [24]. To put the approach efficiently into practice, constitutive equations should be developed which closely mimic the mesoscale behaviour of the fabric. To this end, numerous studies have been conducted aiming to propose accurate and efficient approaches for establishing constitutive equations. Yet again, a detailed review of these approaches could be found in [24]. In almost all conducted studies, it has been demonstrated that the accuracy of the FE 23  models in making predictions of wrinkling strongly depends on the constitutive equations [56][57][58][59]. As such, deriving accurate constitutive equations is the most fundamental step to properly analyze wrinkling behavior of the fabrics. 2.2.2.2 Methods of Modeling Wrinkles and FE Solver Technique Different techniques can be used within FEA to address buckling and instabilities. While the literature in this area in fabric forming is relatively modest, the world of sheet metal forming provides a richer stream of work [21][44], which can be applied to fabric forming. Based on a summary of research works conducted by Henriques et al. [21], Liu et al. [44], and Wang & Cao [15], the methods of modelling and predicting wrinkling in FEA of sheet metals can be classified into four main classes, namely non-bifurcation, bifurcation, eigenvalue, and analytical wrinkling indicator method (Figure 2-8).  Figure 2-8) Methods for modelling wrinkles in sheet metal forming processes. Finite Element AnalysisNon-Bifurcation MethodInitial ImperfectionNumerical Round-off ErrorsBifurcation Method Determinant of Stiffness MatrixEigenvalue MethodWrinkling Indicator MethodPlastic Bifurcation ApproachDynamic Relaxation ApproachEnergy Approach24  The non-bifurcation method is based on the fact that a perfect material under perfect in-plane loading would not buckle at any value [14]. Therefore, some imperfections, for example geometric or material imperfections or load eccentricity, must be introduced to trigger wrinkle formation in an FE analysis. Building on this notion, FE models have been developed with initially embedded imperfections [14][60]. Wrinkling prediction results using this approach are reasonable because of the existence of inherent imperfections in the real material. Having said that, the predictions can be sensitive to the pre-defined imperfections [21]. The bifurcation analysis does not need to have pre-defined imperfections. Instead, the determinant of the stiffness matrix of the structure at each iteration of the increment in the implicit solver is checked and when its sign changes abruptly from positive definite to negative, wrinkling is detected [44]. The eigenvalue buckling analysis identifies the singular points of a stiffness matrix of a structure under a linear perturbation. This method can only be applied to elastic problems. Moreover, effects involving time or strain rate are ignored in the procedure [44]. The analytical wrinkling indicator approach is one of the most practical methods to predict wrinkling in FEA. Wrinkling criteria, which can be developed using a plastic bifurcation approach, a dynamic relaxation method or an energy approach, are integrated into the FE software to locate the zone and the onset of wrinkling. There are two commonly used types of finite element solver, implicit and explicit, both of which can be employed for the analysis of mechanical problems depending on the nature and the conditions of the problem. The implicit method is usually used for solving linear or non-linear static problems or low-speed quasi-static processes. The explicit method is a true 25  dynamic procedure, originally developed to model high-speed impact events in which inertia plays a dominant role in the solution [61]. However, this solver method has also proved valuable in solving quasi-static problems as well as the dynamic problems that it was originally developed for. In forming simulations, therefore, both solvers could potentially be used because of the quasi-static nature of the process. A key question is whether both types of solver give the same wrinkling predictions. If we model a simply supported, thin elastic plate in ABAQUS, using this program's implicit solver (ABAQUS Standard), then the plate stays flat under the action of an in-plane compressive displacement, see Figure 2-9 (a). However, the same analysis using the corresponding explicit solver (ABAQUS Explicit), gives a buckling response for the same loading and boundary conditions, Figure 2-9 (b). Modelling of wrinkles using an explicit solver is considered as a type of non-bifurcation method. That is to say, the accumulation of numerical errors in the explicit solvers act as the imperfections in the model, and wrinkling is simply modelled without introducing pre-defined imperfections or using any other applicable methods [15][21][44]. Although the explicit method can automatically generate deformed shapes with wrinkles, the onset and growth of the buckling obtained from the explicit code is sensitive to the input parameters of the FE model, such as element type, mesh density, simulation speed, etc. [44]. The existence of these types of sensitivity means that FE models with explicit integration method can be unreliable in terms of prediction of wrinkling. Thus, a sensitivity analysis of the FE parameters is highly recommended when using explicit solvers to model buckling [15].  26   Figure 2-9) A simply supported, in-plane loaded, thin elastic plate in (a) ABAQUS Standard, (b) ABAQUS Explicit.  In contrast to the solution using an explicit solver, the implicit technique does not show any wrinkles in the FE model illustrated in Figure 2-9 (a), on account of its accuracy in solving the problem. Accordingly, the FE model with implicit solver should be equipped either with pre-defined imperfections, a stiffness matrix bifurcation approach or analytical wrinkling indicators to predict/model wrinkling. Linking implicit solvers to the abovementioned approaches means that we can take advantage of the reliability of this solver in complicated problems on one hand, whilst overcoming the numerical difficulty of the instability analysis – buckling – on the other hand [15][44]. Specifically, Wang and Cao [15] demonstrated that the predictions of wrinkling using an implicit solver linked with an analytical wrinkling indicator gives “excellent” agreement with the experimental results. 2.3 Research Objectives and Thesis Framework It seems that almost all the FE modelling of the fabric forming process, to the best knowledge of the author, have adopted an explicit solver approach (Table 2.2). Hence, the 27  accumulation of numerical errors of the solver could be a key factor in triggering wrinkle formation. This observation suggests that it would be worth applying the analytical wrinkling techniques described briefly in Section 2.2.1.3 to fabric forming modelling. The analytical wrinkling criteria introduced in Section 2.2.1.3 could be used as wrinkling indicators in an FE software with an implicit solver to predict wrinkles in a more precise manner. On the other hand, reviewing the literature revealed that the previous efforts in this regard are insufficient and there still exists a need for precise prediction of wrinkling (Table 2-1). Table 2-2) A summary of previous efforts to model wrinkling in finite element software. Methods of Modeling Wrinkles in FE Software Fabric Sheet metal Non-Bifurcation Method Initial Imperfection - Cao et al. [14] Numerical Round-Off Errors  Allaoui et al. [62] Khan et al. [63] Boisse et al. [54] Skordos et al. [25]   Wang et al. [15] Bifurcation Method - Wong et al. [64] Eigenvalue Method - Wong et al. [64] Wrinkling Indicator Method -  Wang et al. [15]  An accurate analytical model for the prediction of wrinkling not only can be used as a wrinkling initiation criterion (i.e., wrinkling indicator in FE software) but also for the parametric study purposes. As such, the main objective of the present thesis is as follows: Main Objective: Propose an accurate analytical model to predict the onset of shear wrinkling in woven composite preforms – the case of interest in preforming field. To achieve the main objective, the thesis is structured as follows: 28  Chapter 3 studies the mesoscale mechanisms of shear deformation of woven fabrics and reveals that there exist two distinct shear deformation modes in woven fabrics – trellising mode and pure multi-scale shear mode. Chapter 4 investigates the nature of wrinkle formation and proposes an analytical model to predict the onset of shear wrinkling of woven fabrics under trellising mode. Chapter 5 generalizes the analytical model introduced in Chapter 4 in a manner to be used for woven fabrics under pure multi-scale shear deformation mode. Chapter 6 deals with some of the potential applications of the proposed analytical model. The proposed analytical model is used as a wrinkling indicator for the prediction of shear wrinkling in an FE software and also, identifies the contributing parameters in shear wrinkling to improve the capability of the fabrics against this phenomenon.       29  3 Chapter Three Mesoscale Mechanisms of Shear Deformation of Woven Fabrics 3.1 Overview This chapter attempts to bring an enhanced insight into the analysis of in-plane (i.e., intra-ply) shear behavior of woven fabrics. The analysis of this type of deformation is important in our study because shear wrinkling, which is the case of interest in preforming field, is initiated as a result of a large in-plane shear deformation. Hence, a comprehensive study of the shear deformation of woven fabrics is a fundamental step toward proposing an analytical criterion for the prediction of shear wrinkling. It has been known in the literature for quite a long time that shear deformation of woven fabrics occurs only as a result of trellising mode, and Pin-Joint Theory (PJT) has been a very well-established theory to describe the kinematics of this deformation mode of the fabrics. Through this chapter, however, this common understanding is challenged by undertaking a multi-scale analysis along with a critical review and integration of the previous experimental, analytical and numerical studies. It is demonstrated that the trellising mode is not the only deformation mode which contributes into shear deformation of woven fabrics and there exists another potential mode which could cause intra-ply shear deformation of the fabrics. In the present study, the latter deformation mode is named pure multi-scale shear deformation mode. The occurrence of trellising and pure multi-scale shear deformation modes highly depends on 30  the boundary conditions applied to the fabric. In this chapter, a new theory is proposed to better represent the shear deformation of woven fabrics under pure multi-scale shear deformation mode. Finally, new analytical equations are derived to predict locking angle (i.e., the most important trigger of shear wrinkling) of woven fabrics under trellising and pure multi-scale shear deformation modes. The newly proposed equations notably improve the theoretical prediction of locking angle compared to the equations available in the literature. 3.2 Shear Behavior of Woven Fabrics 3.2.1 Shear Deformation of Woven Fabrics under Trellising Mode A typical shear response of woven fabrics is illustrated in Figure 3-1. According to this figure, the fabric can experience four different regimes during shear deformation, which are referred to as shear with static friction, dynamic friction, locking, and wrinkling.  Figure 3-1) A typical force-shear angle response of woven fabrics at macro-level, showing four phases: 1- shear with static friction, 2- shear with dynamic friction, 3- locking, and 4- wrinkling.  To describe the shear behavior of woven fabrics under trellising mode, the pin-joint theory has been extensively used in the literature [23]. Relying on the kinematics of the structures 31  governed by the PJT, yarns are analogous to rigid trusses pinned to each other at crossovers, i.e., experiencing fully rigid body rotation with respect to each other and creating a trellis structure. This implies that the width of yarns remains constant over the trellising shear deformation prior to locking. Figure 3-2 (a) shows a representative element of a plain-woven fabric. To apply a trellising shear deformation to this fabric element – changing the angle between the warp and the weft yarns – imaginary pins A, B, C, D must undergo the displacements illustrated in Figure 3-2 (a).  Figure 3-2) (a) The motion field of a representative fabric element under trellising mode, (b) velocity vectors during deformation. The centers of rotation are shown as CRi. l is the initial distance between the central axes of the adjacent parallel yarns, and S represents the length of yarn in the fabric element.  In fact, the displacement of the pins leads to rotation of yarns with respect to each other about the interlacing pins. For instance, Yarn 4 experiences relative rigid rotation to Yarn 1 about pin D. However, the actual center of rotation of each yarn as indicated in Figure 3-2 is CRi. Using the concept of dynamics of mechanisms, the center of rotation for each yarn can be identified through a geometrical approach by drawing perpendicular lines to the velocity 32  vectors of two pins of each yarn. Figure 3-2 (b) is presented to evidently explain why it is presumed that the yarns are pinned at the crossovers. This figure depicts the different velocities of the points at the same location but within two distinct yarns. To illustrate, points E1 and E4 which in turn belong to the Yarn 1 and Yarn 4 have the same velocity magnitude, but in different directions. On the other hand, the direction of the velocity vector of points F1 and F4 is identical, but the magnitude is different because of different distances from the corresponding centers of rotation. The exclusive location at which two overlapping points from two distinct yarns have the same velocity, is the pin’s location, implying that points D1 and D4 remain on each other over trellising shear deformation. Figure 3-3 (a) and (b) illustrate the resultant relative displacement between yarns 1 and 4 at crossover D before and after trellising shear deformation, respectively.  Figure 3-3) In-detail meso-level deformation of a cross-over (a) before deformation; (b) after trellising shear deformation.  Given the detailed explanation of the trellising mode, as the shear deformation in the fabric element proceeds, the yarns continue to rotate about their corresponding instantaneous centers of rotation (CR). This deformation mechanism is illustrated in Figure 3-4. 33    Figure 3-4) A plain-woven fabric element (a) before shear, (b) during shear and prior to locking, (c) at the lock-up.  By taking this figure into account, the understanding of different regimes in Figure 3-1 becomes easier. In the sense that at the beginning of the shear, there is static friction between overlapping yarns which resists against yarn rotation and then, by the growth of shear deformation, the static friction turns into dynamic friction that causes the compliant response of the fabric. As shear deformation continues (i.e., the yarns continue to rotate), the distance between the adjacent parallel yarns decreases and at some point, fiber jamming or locking occurs. After the occurrence of this instant, the growth of shear deformation establishes side-to-side (i.e., lateral) compressive force between the adjacent parallel yarns. This compressive force ramps up by the increase in shear deformation and eventually, causes shear wrinkling of the fabric. Accordingly, determining the shear angle at which locking initiates (i.e., hereafter is called locking angle) is of great importance. Based on a theory proposed by Prodromou and Chen [39], locking can be defined as an instant at which the distance between the adjacent parallel yarns becomes zero. Relying on this assumption, they derived equation (2-2) of the present thesis for the prediction of locking angle of the fabric. However, the accuracy of this theory will be challenged later in this chapter. 34  3.2.2 Characterization of Shear Behavior of Woven Fabrics In the literature, there exist two different experimental methods for the characterization of the shear behavior of woven fabrics, namely Bias Extension Test (BET) and Picture Frame Test (PFT). Figure 3-5 illustrates the general set-up for the BET and PFT. In addition to offering a simpler set-up, the BET is known to be more analogous to real stamping process of fabrics with respect to the possibility of inducing both wrinkling and translational slippage defects, depending on the chosen dimensions of samples. On the other hand, this test suffers from a heterogeneous deformation field within the sample, whereas the PFT is capable of offering more uniform deformation fields without translational slippage of yarns at crossovers; hence, simplifying the test characterization. As a disadvantage, the sample installation and fixture boundary effect in the PFT can cause very large deviations in the characterization results. In the following section, these two shear characterization procedures are explained in more details.  Figure 3-5) Two different test set-ups commonly used for the shear characterization of woven fabrics, (a) bias extension test, (b) picture frame test [29]; in both tests, the yarns are initially oriented at ±45̊ with respected to the vertical loading direction.                                            (a)                    (b)  𝐒𝐡𝐞𝐚𝐫 𝐚𝐧𝐠𝐥𝐞 = 𝜽𝑷𝑭) 35  3.2.2.1 Bias Extension Test In bias extension test, a rectangular sample is cut and placed into a tensile machine in a manner that the direction of applied tension is in the bias direction of the warp and weft yarns. The rectangular sample in BET consists of three distinct regions (I, II and III) each of which experiences different shear patterns (Figure 3-6) [65]. If the initial length of the specimen is equal to or larger than twice the width, the yarns in the central zone III are free at their both ends. By making the assumptions that the yarns are inextensible in their longitudinal direction and that there exists no slippage between the warp and the weft yarns, the zone III undergoes shear, and its corresponding deformation mode is the ‘trellising’ of yarns – yarn rotation.  Figure 3-6) A sample (a) at the beginning of BET, (b) during BET.  3.2.2.2 Picture Frame Test In picture frame test, a sample is placed into a hinged frame with four rigid arms of equal length as shown in Figure 3-7. If a tensile force is applied across diagonally opposing corners 36  of the picture frame, the frame changes its shape from a square into a lozenge [23]. Hence, the fabric sample placed inside the picture frame – oriented same as the element depicted in Figure 3-7 – undergoes pure in-plane shear deformation. Although the macroscale behavior of the fabrics under BET and PFT appears identical (i.e., in both cases, the fabric undergoes shear deformation), the mesoscopic investigation of the tests reveals another fact. In the majority of past studies, trellising mode, in which the yarns are assumed free at their both ends, and their width remains constant prior to locking, was identified as the only deformation mode occurring in the fabric under both BET and PFT. However, by taking a close look at the fabric elements at the boundaries in PFT and considering the repetitive pattern of the elements, it becomes evident that the mesoscale deformation mechanism of the fabric under picture frame test not only is not identical with that of the bias extension test but also is quite distinct.  Figure 3-7) A typical picture frame test set up. 37  Figure 3-8 illustrates a cross-over of woven fabric under picture frame test. As can be seen in the figure, shear deformation applied to the fabric causes shear deformation within the yarns (i.e., this is known as intra-yarn shear). Hence, the instantaneous width of yarn – w – does not remain unchanged prior to locking and follows a theoretical cosine formula (i.e., w = w0 cosγ, where γ and w0 are the shear angle of the fabric and initial yarn width, respectively). By this observation, the assumption of using trellising shear mode for PFT is no longer valid and another deformation mode is needed for interpreting the behavior of woven fabrics in this test. In the following sections, we will first prove the existence of intra-yarn shear in picture frame test and then, propose a new deformation mode to better explain the behavior of woven fabrics under PFT.  Figure 3-8) In-detail meso-level deformation of a cross-over, (a) before deformation; (b) after shear deformation in picture frame test. γ, w0 and w and are the shear angle of the fabric, the initial width and the instantaneous width of the yarns, respectively.  3.2.2.2.1 Existence of Intra-Yarn Shear in Picture Frame Test It can be numerically, experimentally and analytically demonstrated that in-plane shear applied to the fabric in picture frame test causes intra-yarn shear which impacts the response of fabric, notably. In the PFT, the intra-yarn shear is highly possible to arise, compared to the 38  BET, because this phenomenon originates from the boundary conditions applied at the yarn ends. Numerical Evidence - The effect of boundary conditions to create intra-yarn shear was first introduced numerically by Lin et al. [45]. In their numerical simulations where yarns rotation was merely taken into account, the simulations could not properly predict the response of the fabric in picture frame test. Hence, they modified their numerical simulations by taking the intra-yarn shear into account and obtained some results with higher accuracy. Figure 3-9 (adapted from [45]) shows the response of a woven fabric in picture frame test versus the results obtained from finite element simulations by considering trellising mode in which there exists no intra-yarn shear and pure shear mode which takes the intra-yarn shear into account. As figure shows, the model, which considers intra-yarn shear in the analysis, yields the results closer to the experimental measurements, proving the existence of intra-yarn shear in the picture frame test.  Figure 3-9) (a) Shear force – shear angle diagram, (b) fabric element under trellising shear mode (only yarn rotation causes the shear of the fabric element), (c) fabric element under pure shear mode (taking the intra-yarn shear into account) [45]. 39  Experimental Evidence - The performed micro-level optical measurements of strain field within yarns of a plain-woven fabric under picture frame test revealed the presence of shear deformation within yarns [66]. Figure 3-10 (adapted from [66]) shows strain field within a yarn during a picture frame test, proving the existence of intra-yarn shear in PFT.   Figure 3-10) Comparison between intra-yarn shear (i.e., in this figure, it is named by local shear angle) and the applied frame shear angle [66]. RR1, RR2, and RR3 are three different fabrics, the stiffest of which is RR3. A and B represent the shear angle within yarns spacing and the yarn itself.  Another experimental investigation implicitly highlighted the occurrence of intra-yarn shear mechanism, by tracking the width of yarns over shear deformation [67]. If the trellising shear mode, which is based on PJT, would happen in the PFT, the width of yarns should have remained constant up to the locking point. However, results in Zhu et al. [67] clearly reveal that the width of yarns reduces during shear deformation even before reaching the locking angle (Figure 3-11), offering another evidence for the presence of intra-yarn shear in PFT.  40   Figure 3-11) Yarn width change during shear deformation of the fabric in PFT [67].  Theoretical Evidence – Regarding theoretical evidence for the existence of shear deformation within yarns, although brief explanations – using fabric deformation continuity at fixture boundaries – has been provided in [68], a detailed theoretical explanation is presented here to better explain the existence of intra-yarn shear in PFT. If we assume a yarn under PFT as a curved beam which undergoes transverse shear as indicated in Figure 3-12, and ignore intra-yarn shear, then only bending and torsion energies contribute into the shear force-shear angle curve illustrated in Figure 3-1. Considering a beam subjected to a force along its thickness direction, it is common to overlook the contribution of transverse shear energy in the total strain energy of the beam, providing its length to thickness ratio is greater than 16 (ASTM D7264). However, the ignorance of the transverse shear energy in the yarns of woven fabrics subjected to in-plane shear loading is highly questionable, in that the ratio of the length to the thickness of an equivalent beam, representing a typical yarn in the fabric element, is only slightly greater than 2. To clarify, Figure 3-12 (a) demonstrates a typical woven fabric which is modeled in TexGen. According to this figure, the length of the yarn is the length of the modelled cantilever beam, while the thickness of the cantilever beam is the 41  width of yarn, because the in-plane shear force is along the width of yarn (Figure 3-12 (b)). Consequently, for related meso-level deformation analysis purposes, the ratio of length to the thickness of the proposed equivalent beam is the ratio of the length to the width of yarn in a representative fabric element, whose value is nearly 2 given the fact that length of a yarn in a fabric element of a balanced woven fabric is around twice the width of yarn plus the spacing between yarns as shown in Figure 3-12 (a).    (a) (b) (c) Figure 3-12) Schematic of a fabric element of a typical plain-woven fabric modeled in TexGen (a) top view, (b) 3D presentation of a single yarn in the fabric element subjected to in-plane shear force, (c) a simplified structural beam representation. Notice that the in-plane shear force is along the width of yarns, which is generally not substantially less than the length of a yarn between two crossovers.  Accordingly, instead of Euler-Bernoulli beam theory, the High-Order beam theory or the First-Order Beam Theory (Timoshenko Beam Theory) should be employed. Under such theories, the stored shear energy in an arc beam shown in Figure 3-12 (c) should be added to the bending and torsional deformation energies (i.e., it means that the shear energy in the beam is not negligible compared to the bending energy; therefore, it needs to be taken into account in the analyses) as follows: 42  𝑈 = ∫𝑀𝑦22𝐸𝐼𝑦𝑦𝑑𝑥 + ∫𝑇𝑥2𝐺𝐽𝑑𝑥 + ∫𝑉𝑧(𝑥)22𝐺𝐴′/𝛼𝑑𝑥 (3-1) where E, I, G, J, A', and 𝛼 represent the Young’s modulus, bending moment of inertia, shear modulus, polar moment of inertia, cross-sectional area, and cross-section shape coefficient of the equivalent beam, respectively. Also, the resultant bending moment, torsional torque, and shear force at each point of the beam are regarded as 𝑀𝑦, 𝑇𝑥, and 𝑉𝑧, respectively. As a consequence, it is now evident that during the course of picture frame test, the yarns also undergo shear deformation that proves the presence of intra-yarn shear. Another theoretical argument to substantiate the existence of intra-yarn shear in the PFT can be offered in regards to kinemics of such tests from a continuum mechanics point of view. Milani et al. using a continuum mechanics based model proposed an explicit predictive equation for the shear response of woven fabrics under fiber misalignment [69]. If X and 𝑥 vectors denote the undeformed and deformed vector of a material point within the fabric subjected to PFT (Figure 3-13), the relationship between components of these vectors under a fully uniform picture frame test without local bending close to grips – ideal PFT – can be attained as: {𝑥1 = 𝑋1 cos𝛾2+ 𝑋2 sin𝛾2𝑥2 = 𝑋1 sin𝛾2+ 𝑋2 cos𝛾2 (3-2) Thereafter, the deformation gradient tensor can be calculated as: 43  𝐹 = [𝜕𝑥𝑖𝜕𝑋𝑗] =  [     cos𝛾2sin𝛾20sin𝛾2cos𝛾200 0𝜕𝑥3𝜕𝑋3]      (3-3) Given the incompressibility assumption, det[F]=1, and it leads to determination of the unknown component 𝜕𝑥3𝜕𝑋3 as 1cos𝛾. Existence of gaps between filaments within yarns, however, challenges the incompressible behavior assumption for woven fabrics. Comparison between the above presented expression for change in thickness (1cos𝛾) with the experimental results observed in [67] suggests that yarns are slightly thickened during shear deformation, but not as significantly as the proposed formula by the continuum mechanics based model in [69]. Accordingly, 𝑐𝑔, a correction factor representing the impact of gaps between filaments within a dry yarn on the incompressibility assumption, is introduced here. Subsequently, the right Cauchy–Green tensor and the strain tensor, known as C and E, respectively, can be calculated as: 𝑪 = [𝑭]𝑇[𝑭] = [1 sin 𝛾 0sin 𝛾 1 00 0 (𝑐𝑔 sec 𝛾)2] (3-4) 𝑬 =12[𝑪 − 𝑰]=12[0 sin 𝛾 0sin 𝛾 0 00 0 (𝑐𝑔 sec 𝛾)2 − 1] (3-5) Equation (3-5) reveals that the fabric sample subjected to uniform picture frame loading must experience pure shear deformation at every single point. To elucidate more, thanks to 44  arbitrarily chosen point P, it can be located within yarns spacing regions as well as within yarns. Hence, shear deformation field is created not only between yarns, but also inside the yarns.   (a) Before deformation (b) After deformation  Figure 3-13) Coordinate systems to demonstrate the undeformed and deformed shape of a fabric sample subject to an ideal picture frame test.  In view of the numerical and experimental studies reviewed above and the conducted theoretical investigation, it is apparent that trellising shear mode in which the yarn width is assumed unchanged prior to locking in shear deformation (i.e., there exists no intra-yarn shear) does not properly fit the case of picture frame test. As such, a new deformation mode which closely mimics the mesoscopic behavior of woven fabrics in PFT seems necessary. 3.2.3 Proposing a New Shear Deformation Mode for Woven Fabrics As mentioned earlier, the trellising mode in which yarns experience rigid body rotation to create shear in the fabric is not valid for the picture test because of the boundary conditions applied to ends of yarns; instead pure multi-scale shear deformation mode would be more realistic. This deformation mode is analogous to pure shear deformation of a continuum 45  element in which every single point of the element undergoes shear deformation. This simply means that the yarns in the fabric element go through shear deformation (i.e., the fabric experiences intra-yarn shear). Although the existence of intra-yarn shear in PFT has been observed and reported in the literature, its impact on the shear behavior of the fabric has not been investigated yet. Accordingly, this section of the thesis aims to demonstrate that how intra-yarn shear can significantly affect fabric behavior and why this mesoscale mechanism must be considered in the analyses of woven fabrics. Figure 3-14 (a) shows a representative fabric element in picture frame test. The arrows in the figure depict the shear force applied to the fabric element. By applying a shear force to the ends of yarns, the fabric element experiences pure shear deformation. As can be seen in Figure 3-14 (b) and (c), with the increase in the shear deformation of the fabric element, the instantaneous yarn width (w) along with the instantaneous length of the considered element (𝐿) decrease.  Figure 3-14) (a) A representative fabric element under pure multi-scale shear deformation mode, (b) deformation of a yarn segment during the picture frame test, showing intra-yarn shear and hence a decrease in the instantaneous yarn width, (c) deformation of the representative fabric element. 46  The changes in the yarn width and the length of fabric element lead to the change in spacing between the adjacent parallel yarns. The spacing between the adjacent parallel yarns (SP) in the fabric element is determined by the following mathematical equation:  𝑆𝑃 = 𝐿 − 2𝑤 = 𝐿0. cos(𝛾) − 2(𝑤0 cos(𝛾)) (3-6)  in which L0, w0 and 𝛾 are the initial length of the fabric element, the initial yarn width and the shear angle of the fabric element, respectively. Figure 3-15 illustrates the changes of the spacing between adjacent parallel yarns versus shear deformation of the fabric element. It is seen that the spacing between the adjacent parallel yarns tapers off, as shear deformation of the fabric element increases. Since the yarn width continuously changes over shear deformation, it is expected that the locking angle of the fabric under pure multi-scale shear deformation mode is notably different than that of the trellising mode.  Referring to the definition of locking provided by Prodromou et al. [67] in which the spacing between the adjacent parallel yarns should become zero (SP = 0), we expect to obtain the locking angle of the fabric under pure multi-scale shear deformation mode, analytically. However, the equation below does not provide any solution unless 2w0 and L0 are equal:  𝑆𝑃 =  𝐿0. cos(𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔) − 2(𝑤0 cos(𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔)) = 0 (3-7)  This geometrically means that there should not be any spacing between the yarns in a given fabric, which is not practical. The key point resolving the dispute presented above is hidden in the assumption given in Prodromou’s theory in which “locking occurs when the distance between the adjacent yarns becomes zero”. Their erroneous assumption stems from looking at the fabric element in a 2D schematic. Instead, if we consider a 3D fabric element similar to Figure 3-16, locking occurs when each yarn of a given family is sandwiched between 47  two adjacent yarns of the other family. This in turn implies that the critical distance between two parallel yarns to initiate locking is the ‘effective’ thickness of the other family of yarn (i.e., interlacing yarn), rather than zero. Mathematically, this point adds one more term to equation (3-7), leading to equation (3-8):  𝑆𝑃 = 𝐿0. cos(𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔) − 2(𝑤0 cos(𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔)) = 𝑡 (3-8)  where t is the yarn thickness. Now, based on our newly proposed definition of locking instant, the locking angle of a given woven fabric under pure multi-scale shear deformation mode can be analytically determined as shown in equation (3-8). Moreover, the analytical equation for the prediction of locking angle of woven fabrics under trellising mode should be updated by incorporating the effect of interlacing yarn as indicated in next section.  Figure 3-15) Trajectory of the spacing between the adjacent parallel yarns over pure multi-scale shear deformation mode, which is drawn using equation (3-6).  SP0 48    Figure 3-16) The effect of interlacing yarn on locking angle by considering 3D schematic of the fabric.  3.2.4 Modifying Prodromou’s Equation for Prediction of Locking Angle of Woven Fabrics under Trellising Mode Although the interpretation of locking of woven fabrics under trellising mode proposed by Prodromou et al. [67] was a good starting point, they have ignored the fact that locking occurs when the distance between the adjacent parallel yarns is bounded by the thickness of the interlacing yarn. To elaborate on, if we consider the fabric element illustrated in Figure 3-17, locking of the fabric occurs when OD, which is initially equal to 𝑦𝑎𝑟𝑛 𝑤𝑖𝑑𝑡ℎ + 𝑔𝑎𝑝2, becomes equal to 𝑦𝑎𝑟𝑛 𝑤𝑖𝑑𝑡ℎ + 𝑡ℎ𝑖𝑐𝑘𝑛𝑒𝑠𝑠2. Built on this newly proposed geometric interpretation of locking, we derive another equation to predict the locking angle of plain-woven fabrics under trellising shear deformation mode: 𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔 = 𝑐𝑜𝑠−1 (𝑤0 + 𝑡𝑙) (3-9) where γlocking, w0, t, and l are the locking angle, initial yarn width, yarn thickness and initial distance between the central axes of the adjacent parallel yarns, respectively (Figure 3-17). Spacing Interlacing Yarn 49   Figure 3-17) A fabric element at the beginning of trellising shear deformation mode.  Figure 3-18 shows that how significant is the incorporation of the thickness of interlacing yarn on the prediction of locking angle. An experimental comparison is also drawn in Chapter 4 of the present thesis, demonstrating that the predictions by our newly proposed equation are closer to the experimental measurements compared to those predicted by the model introduced in [39].   Figure 3-18) Predictions of the onset of locking in a fabric element under trellising mode using equation (3-9) and equation (2-2) [39]. The width and thickness of the yarns and the spacing between the adjacent parallel yarns are assumed 1.4 mm, 0.355 mm and 0.65 mm, respectively.  Furthermore, a comparison between the analytical predictions of locking angle of an identical plain-woven fabric under trellising shear mode and the pure multi-scale shear mode is drawn in Figure 3-19. This compassion reveals that how different are the predictions of Equation (2-2) Equation (3-9) 50  locking angle under different shear deformation modes. Hence, it is of great importance to distinguish different shear modes of woven fabrics for the further analyses. Note that the experimental tests are also conducted in Chapter 4 and Chapter 5 of the thesis to validate the theoretical predictions.      Figure 3-19) A comparison of the predictions of locking angle by relying on trellising shear mode and pure multi-scale shear mode. The width and thickness of the yarns and the spacing between the adjacent parallel yarns are assumed 1.4 mm, 0.355 mm and 0.65 mm, respectively.  3.3 Summary The mesoscopic analysis of BET and PFT revealed that there exist two different shear deformation modes in woven fabrics, namely trellising shear mode and pure multi-scale shear mode. Despite the fact that in a few studies the existence of intra-yarn shear in PFT had been observed, there was not a clear understanding that how this mechanism would affect the shear behavior of woven fabrics. Built on the existence of intra-yarn shear in PFT, this chapter clearly demonstrated that trellising mode is not appropriate to describe the shear behavior of woven fabrics under picture frame test. Instead, another mode, which is called pure multi-scale shear deformation, was used to better represent the behavior of the deforming fabrics under PFT. Moreover, through this chapter new analytical equations were provided to predict the onset of Trellising Mode Pure Multi-scale Shear Mode    51  locking of plain-woven fabrics under trellising and pure multi-scale shear deformation modes. The predictions revealed that the difference between the locking angles of an identical fabric under these two modes is substantially large. This conclusion is important since locking is the most influential trigger in shear wrinkling phenomenon. That is, the larger locking angle, the later shear wrinkling of the fabric. Accordingly, it is very crucial to understand the nature of shear deformation of woven fabrics to be able to precisely predict shear wrinkling. What’s more, in the real fabric forming process, shear deformation mode of woven fabrics could be either trellising or pure multi-scale shear depending on the boundary conditions applied to the fabric. To elucidate, if the fabric is free at its ends – there is no or little blank holder force, then the trellising is the dominant shear deformation mode. However, if the fabric is constrained around its periphery by a blank holding system, then the pure multi-scale shear mode is the dominant mode of deformation. As such, shear wrinkling prediction based on both deformation modes would be of great importance in the academia and the industry. In the following sections, therefore, analytical equations are derived for the prediction of shear wrinkling of woven fabrics based on trellising mode and pure multi-scale shear mode.    52  4 Chapter Four A Mesoscopic Analytical Model to Predict the Onset of Wrinkling in Plain-Woven Preforms under Trellising Shear Deformation Mode 4.1 Overview This chapter presents a mesoscopic analytical model for the prediction of wrinkling of Plain-Woven Composite Preforms (PWCPs) under trellising shear mode, based on a new instability analysis point of view. The analysis aims to pave the way toward better understanding of the nature of wrinkle formation in woven fabrics due to trellising shear mode and identify the effect of contributing fabric and process parameters on the onset of wrinkling. To this end, the mechanism of wrinkle formation in PWCPs in mesoscale is first scrutinized, and an equivalent structure composed of bars and different types of springs is proposed, mimicking the behavior of a representative PWCP element at the post-locking stage. The parameters of this equivalent structure are derived based on geometric and mechanical characteristics of the PWCP. Among different instability analysis techniques, the principle of minimum total potential energy is employed, and experimental validation is carried out to reveal the effectiveness of the derived prediction equation. 53  4.2 Mechanism of Wrinkle Formation under Trellising Shear Mode Wrinkles in the PWCPs can be formed through three distinct loading conditions namely in-plane compression, out-of-plane bending and large in-plane shear deformation (Figure 4-1). As stated earlier, in-plane shear deformation of woven fabrics can occur through two distinct deformation modes, namely trellising and pure multi-scale shear deformation. Trellising was described based on the pin-joint theory in which in the initial stage of shear deformation of PWCPs, yarns rotate with little resistance, which is mainly due to frictional forces between overlapping families of yarns. Afterward, the shear resistance of PWCPs increases more rapidly when inter-yarn gaps vanish – locking initiates – and adjacent parallel yarns begin to press each other transversely (Figure 4-2). Further shear deformation of PWCPs at this stage is dominantly by yarns lateral compaction. By growing the shear deformation of PWCPs, the pressure between adjacent parallel yarns rises, triggering the formation of wrinkles.         Figure 4-1) Wrinkles formed due to (a) in-plane compression, (b) out-of-plane bending, and (c) in-plane shear.   (a) (b) (c) 54        Figure 4-2) A PWCP in the course of large in-plane shear deformation.  In light of the above discussion, the PWCPs in advance of wrinkling experience two main stages, “pre-locking” and “post-locking”. As demonstrated in Chapter 3, the locking of the PWCPs mainly depends on kinematics and the geometry of the yarns in the woven fabric structure: it is defined as an instant at which the distance between the adjacent parallel yarns is bounded by the thickness of the interlacing yarn. Contrary to the pre-locking stage, the behavior of PWCPs in the post-locking stage is dominantly governed by the material properties of the yarn (e.g., its transverse compaction resistance and flexural rigidity). Hence, the onset of wrinkling in PWCPs as a whole depends not only on kinematics and geometry of the fabric elements but also the material properties of the yarns. Although much earlier research on the behavior of PWCPs in the pre-locking stage has been conducted [39], the mechanism of deformation in the post-locking stage has not yet fully understood. In other related references [9][42], it was reported that the smallest fabric sub-structure, in which wrinkles can potentially initiate, is the representative element illustrated in Figure 3-17. Accordingly, this element is analyzed as the source of wrinkling initiation in the PWCP Undeformed Prior to locking Locked position 55  structure to identify the onset of shear wrinkling [70]. Figure 4-3 illustrates the compressive forces establishing between adjacent parallel yarns during the course of lock-up, and thereof Figure 4-4 depicts the mechanism of wrinkle formation in the representative element of PWCP.        Figure 4-3) Compressive forces induced between adjacent yarns during trellising shear deformation of PWCP; (a) simulation, (b) experimental investigation of the yarn strain field in large shear deformation [3].  In mechanical structures under compressive forces, there are discrete values of the load at which secondary equilibrium configurations may appear in the neighborhood of the initial equilibrium position. To put it differently, a system under critical compressive loads admits additional but adjacent equilibrium states with different deformation patterns called structural modes [17]. The nature of shear wrinkling in PWCPs can be analogous to a deviation in deformation mode of the material structure. This analogy opens up a window toward theoretical analysis of shear wrinkling phenomenon in PWCPs from an instability point of view.  Top View Bottom View (a) (b) Y-3 Y-4 Y-1 Y-2 Y-1 Y-2 Y-3 Y-4 56          Figure 4-4) Mechanism of wrinkle formation in a PWCP element due to large trellising shear deformation (T and Ms denote the tension along yarn and in-plane shear coupling, respectively).  In the next section, an analogy-based equivalent structure for the representative element of PWCP, which mimics the behavior of the fabric at the post-locking stage, is introduced. Subsequently, the analytical criterion at which deviation of deformation mode from the in-plane to out-of-plane mode becomes possible is derived and validated experimentally. 4.3 Prediction of Fabric Wrinkling Based on an Instability Analysis 4.3.1 Analogy-Based Equivalent Structure As alluded in Figure 4-3 and Figure 4-4, the adjacent parallel yarns of the fabric element undergo compressive forces during shear deformation, and subsequently, at a certain instant termed the onset of wrinkling, the fabric structure alters its deformation pattern from the in-plane mode to the out-of-plane mode. Given this definition, a simple four-link equivalent Before Wrinkling After Wrinkling T T T T T T T T MS MS 57  structure mimicking the behavior of a representative PWCP element in the post-locking stage is proposed as shown in Figure 4-5.  Figure 4-5) A four-link equivalent structure mimicking the behavior of a representative PWCP element in the post-locking stage.  It is assumed that the instant of deviation in equilibrium mode of this equivalent structure is coincidental with the onset of wrinkling in the fabric (Figure 4-6). This structure, hereafter called PWCP Equivalent (PWCPE), comprises of four rigid bars of length w0 (initial width of yarns). These bars are clamped elastically at one end against rotation – about the dash-lined axes illustrated in Figure 4-5 – using four torsion springs of stiffness Kt, and linked to ball-joints at the other end (Point O). The torsion springs imitate the role of interlacing yarns in the fabric (Figure 4-7). To elucidate, yarn Y-2 at the onset of wrinkling tends to rotate about its longitudinal outer edge (i.e., about the JJʹ axis in Figure 4-7) and form a wrinkle. To accommodate such deformation, the transverse bending of the interlacing yarn Y-3 is necessitated as shown in Figure 4-8. That is, the interlacing yarn Y-3 resist against wrinkling. Hence, this out-of-plane resistance of interlacing yarn is modeled by an equivalent torsion 58  spring (Kt). The other yarns in the representative PWCP element simultaneously experience the same sequence of events.         Figure 4-6) Equivalent structure (a) before equilibrium deviation, (b) after equilibrium deviation.  Four sliding springs with the stiffness coefficient K are placed around the bars to capture the lateral compaction stiffness of yarns – in the width direction – and their stiffness coefficient needs to be characterized via experimentation (to be discussed in more details in Section 4.4.2.2). The bars are treated as the cores of the sliding springs. The quasi-static compressive forces P act on the system as illustrated in Figure 4-6. These forces are to model the compaction forces between adjacent parallel yarns in the real fabric. If the compressive forces P in the PWCPE structure reach a critical value termed as Pcr, the equivalent structure will become unstable. At the onset of instability, the PWCPE structure will change its deformation mode from in-plane to out-of-plane (Figure 4-6). Accordingly, obtaining Pcr is a significant step to O A C B D (a) (b) P P P P MS T T T T T T T T Pcr Pcr Pcr Pcr O' A' C' B' D' MS 59  predict the onset of wrinkling in the fabric structure. Before that, however, we first need to characterize the stiffness properties of the structure as follows.   Figure 4-7) Equivalent torsion spring Kt, mimicking the bending effect of interlacing yarn Y-3 on the out of plane deformation of Y-2, and vice versa.  4.3.2 Characterization of Stiffness Elements of the PWCPE Structure As stated earlier, the stiffness coefficient of the sliding springs, which model the lateral compaction of the yarns, needs to be characterized via experimentation. However, the stiffness of torsion springs is characterized analytically for a given fabric. To derive the equivalent torsion spring coefficients, it is presumed that yarn Y-3 is a straight beam (Figure 4-8) made of a linear, transversely isotropic, homogeneous material with an effective flexural rigidity in the longitudinal direction Qb [71]. 60   Figure 4-8) Modeling of yarn as a bending beam (made of a linear, transversely isotropic, homogeneous material with an effective flexural rigidity in the longitudinal direction Qb), and an equivalent torsion spring coefficient, Kt.  The Euler-Bernoulli theory is used to analyze such hypothesized beam [71][67]. Moreover, as the proposed model takes the effect of tension on the yarns into account, instead of a first-order analysis which is based on the superposition principle and considers the deflection of the beam independent of the applied tension, the second-order analysis [72] is undertaken which accounts for the effect of tension on the overall flexural rigidity of yarns, F YS/2. Considering the assumptions made, the deflection of the yarn as a beam at its mid-span and under arbitrary vertical force F can be written as [18]: 𝑌𝑆/2 =𝐹 2𝑄𝑏  (𝑇𝑄𝑏)32( 𝑆√𝑇𝑄𝑏2− tanh( 𝑆√𝑇𝑄𝑏2) )  (4-1)  where T, 𝑄𝑏, and S are tension along the yarn, the effective flexural rigidity of yarn in its longitudinal direction, and the length of yarn in the representative PWCP element (Figure 3-17), respectively. Note that the second order analysis determines the extent of increase in bending stiffness for an integrated beam. However, yarns are comprised of thousands of filaments. In simple words, the effective bending rigidity of a yarn, Qb, where filaments are only in contact with each other and can slide over one another is much less than when the 61  filaments are perfectly jointed to each other and make an integrated beam. This conflict may be resolved by the fact that measuring the flexural rigidity using the setup explained in Section 4.4.2.1 indirectly takes the discrete nature of the yarns into account. However, one point has still been overlooked. To elaborate, tension in the yarns increases the overall flexural rigidity of yarns, F YS/2, through two mechanisms:  Taking the resultant bending moment from tension into account (second-order analysis) reduces the quantity YS/2 and in turn, increases the overall flexural rigidity. This issue is incorporated in equation (4-1).  Increase in the internal bending rigidity of the yarns by making more interaction/friction between filaments, meaning that they are virtually jointed and make an integrated beam with the total cross-section of sub-beams. Mathematically, Qb in equation (4-1) should be presumed a function of tension. This mechanism has not been a part of the present study. In addition to the assumptions stated above, it should be highlighted that in the current study, in spite of nonlinear tensile behavior of yarns, the linearity assumption does not affect the results since the level of tension in our experimental investigation is kept at zero.  Now, by equating the deflection of the mid-span of the beam (𝑌𝑆/2), which is caused by the arbitrary force F, to the nodal displacement of the end of the link caused by the same force (𝑶𝑶′ in Figure 4-8), the equivalent stiffness of the torsion spring is determined as: 62  𝐾𝑡 =2𝑄𝑏 (𝑇𝑄𝑏)32𝑤2( 𝑆√𝑇𝑄𝑏2) − tanh( 𝑆√𝑇𝑄𝑏2)  =  2𝑄𝑏 (𝑇𝑄𝑏)32(𝑤0 − ∆)2( 𝑆√𝑇𝑄𝑏2) − tanh( 𝑆√𝑇𝑄𝑏2)  (4-2)  where w, w0 and Δ are, in turn, the instantaneous width of yarn, initial width of yarn and compaction of the yarn due to lateral compressive force. For subsequent calculations, it is convenient to define the following lumped parameter notation: 𝐴 =  2𝑄𝑏 (𝑇𝑄𝑏)32( 𝑆√𝑇𝑄𝑏2) − tanh( 𝑆√𝑇𝑄𝑏2)  (4-3)   4.3.3 Determination of Critical Compressive Force (Pcr) and Shear Angle at the Onset of Wrinkling As discussed in Section 4.3.1, the onset of wrinkling in the fabric element is coincident with the instant at which a deviation of deformation mode in the PWCPE structure occurs. As such, the methods developed in the field of structural analysis can be borrowed to estimate the onset of deviation in the equilibrium mode (wrinkling) of PWCPE. One of the de facto analytical approaches is the principle of minimum total potential energy. This approach is based on the fact that the mechanical system is conservative (i.e., there is no friction in the system). However, we are aware that there exist two sources of energy loss due to friction. The first is the energy dissipated through the friction at crossovers. Because this is substantially less than the lateral compaction energy [9], it does not violate the assumption of the approach. The second source of energy dissipation is the filament sliding within the yarns which is under 63  lateral compaction. This mechanism causes an elastoplastic behavior of the yarns in the lateral (compaction) direction. Since the compaction of filaments is the main source of resistance of yarns against lateral compression, the sliding friction between filaments can be also ignored. Moreover, K is assumed linear whose reason is explained in the next section. By taking Figure 4-9 into consideration, the total potential energy of the system in an imaginary deviated mode is:   𝐿(𝜃, ∆) = 4(𝐾𝑡𝜃22+𝐾∆22− 𝑃(𝑤0 − (𝑤0 − ∆)𝑐𝑜𝑠𝜃)) (4-4)   where   is the angle of rods with respect to the plane of initial mode, Δ is the compaction of the slide springs, and P is the compressive forces exerted on the sliding springs. In equation (4-4), term 1 and term 2 are in turn the strain energy of the system – negative of the work done by the internal forces – and potential energy of the external forces – negative of work done by external forces [16]. Since the bars in the equivalent structure are assumed rigid (merely as the cores of the sliding springs), there exists no term of strain energy for these members. Owing to the assumption of conservation of energy in the system, the equilibrium of the PWCPE structure in the initial mode is: {𝜕𝐿𝜕𝜃= 0          ⇒          𝐾𝑡𝜃 − 𝑃(𝑤0 − ∆)𝑠𝑖𝑛𝜃 = 0          ⇒           𝜃 = 0 𝜕𝐿𝜕∆= 0          ⇒           𝐾∆ − 𝑃𝑐𝑜𝑠𝜃 = 0                           ⇒           ∆=𝑃𝐾 (4-5)  Term 2 Term 1 64   Figure 4-9) The PWCPE structure in an imaginary deviated mode; note that the P forces remain in the original directions.  On account of the fact that the PWCPE structure is a system with multiple degrees of freedom, the type of extrema in total potential energy function, L (, Δ), can be obtained with the aid of the Hessian matrix. If the determinant of the matrix at the equilibrium points is positive definite and 𝜕2𝐿𝜕𝜃2> 0, then the total potential energy of the system has a relative minimum, inferring that the equilibrium mode of the structure does not intend to deviate to an adjacent equilibrium state – the Lagrange-Dirichlet theorem [73]. However, if the matrix is either negative definite or is not a sign definite, the equilibrium mode of the structure opts to deviate to an adjacent equilibrium mode – first Lyapunov theorem [73]. Hence, it is inferred that the instant at which the determinant of the Hessian matrix of the structure changes its sign from positive definite to negative definite, may be considered as the onset of equilibrium mode deviation. The determinant of Hessian matrix for the PWCPE structure is:  65  𝐻 = |[4𝐾𝑡 − 4𝑃(𝑤0 − ∆)𝑐𝑜𝑠𝜃 4𝑃𝑠𝑖𝑛𝜃4𝑃𝑠𝑖𝑛𝜃 4𝐾]|𝜃=0∆=𝑃𝐾 (4-6)  The instant at which deviation in deformation mode occurs is defined by H=0. It should be noted that in this structure if the determinant of the Hessian matrix is positive, 𝜕2𝐿𝜕𝜃2 is always positive. By considering equations (4-2) and (4-3), it is evident that Kt depends on Δ and A. Hence, we have:  𝐻 = |[4𝐴(𝑤0 − ∆)2 − 4𝑃(𝑤0 − ∆)𝑐𝑜𝑠𝜃 4𝑃𝑠𝑖𝑛𝜃4𝑃𝑠𝑖𝑛𝜃 4𝐾]| 𝜃=0∆=𝑃𝐾= 0 (4-7)  Substituting Δ by P/K and  by zero, a quadratic equation with respect to the compressive loads is found, whose roots could identify the critical compressive load Pcr:  (1 +𝐴𝐾)𝑃2 − (2𝐴𝑤0 + 𝐾𝑤0)𝑃 + 𝐾𝐴𝑤02 = 0 (4-8)  And hence,  𝑃𝑐𝑟 =2𝐴𝑤0 + 𝐾𝑤0 ± 𝐾𝑤02(1 +𝐴𝐾) (4-9)  Only the minus sign of equation (4-9) should be considered as the critical load since the sign of equation (4-8) changes from positive to negative at:  𝑃𝑐𝑟 =𝐴𝐾𝑤0𝐴 + 𝐾 (4-10)  Considering equation (4-5), the critical yarn compaction, Δcr, can be determined as: 66  𝛥𝑐𝑟 =𝑃𝑐𝑟𝐾=𝐴𝑤0𝐴 + 𝐾=2𝑄𝑏 (𝑇𝑄𝑏)32𝑤02𝑄𝑏 (𝑇𝑄𝑏)32+ 𝐾(  ( 𝑆√𝑇𝑄𝑏2) − tanh( 𝑆√𝑇𝑄𝑏2) )   (4-11)  Finally, referring to Figure 3-17 and the discussion presented in Section 3.2.1, the shear angle at which wrinkling occurs – critical shear angle – is:  𝛾𝑊𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔  = cos−1 (𝑤0 + 𝑡 − ∆𝑐𝑟𝑙) (4-12)   This implies that as the shear deformation grows, at some point the distance between adjacent parallel yarns reaches the thickness of the interlacing yarn – locking begins; subsequently, the adjacent parallel yarns start compacting each other up to Δcr when eventually wrinkle sets in (equation (4-12)). It should be mentioned that in this analytical model, the compaction of the yarns in thickness direction, which is substantially less than that in the width direction, is ignored. 4.4 Experimental Evaluation 4.4.1 Geometric Characterization A commercially available PWCP made of carbon fiber was acquired from APC Composite Inc. The geometric specifications of the fabric were obtained using a Nikon optical microscope. For geometric characterization of yarns, six measurements at different locations of the fabric were taken, and the measured values were averaged as indicated in Figure 4-10. 67  To measure the yarn width and the spacing between two adjacent parallel yarns, the PWCP sample was laid under the microscope and the measurements were taken. However, measuring the thickness of the yarn was not as straightforward. Six yarns were dipped into PDMS which is a transparent and soft thermoset polymer. After curing the polymer, a section-cut was passed perpendicular to the longitudinal axis of the yarn. Subsequently, the thickness of the yarn was measured under the microscope (Figure 4-10).  Figure 4-10) The geometric characteristics of the tested PWCP. Mean value of yarn width along with its standard deviation are, in turn, 1.4 mm and 0.04 mm. Mean value of spacing between yarns along with its standard deviation are, in turn, 0.65 mm and 0.05 mm. Mean value of yarn thickness along with its standard deviation are, in turn, 0.355 mm and 0.04 mm.  4.4.2 Meso-Mechanical Characterization of the Fabric Mechanical properties of the yarns comprising flexural rigidity along with lateral compaction stiffness of yarns – in the width direction – needed to be experimentally measured. There are some studies targeting characterization of the effective flexural rigidity of the yarns in longitudinal direction [71][74], whereas measuring the lateral compaction stiffness of yarns has not received much attention in the literature. 68  4.4.2.1 Effective Flexural Rigidity of Yarn in Longitudinal Direction (Qb) There exists a common test method for determining the effective bending rigidity of a yarn. The experiment is based on Peirce’s cantilever testing procedure, in which a yarn bends under its own weight [71][74]. Figure 4-11 illustrates the experimental set-up used in this study. Five measurements were taken and the averaged value along with the standard deviation for the effective bending rigidity of the yarn was found to be 2.75 N.mm2 and 0.21 N.mm2, respectively.    Figure 4-11) The experimental set-up used for measuring effective bending rigidity of the yarn.  4.4.2.2 Effective Lateral Stiffness of Yarn (K) 4.4.2.2.1 Direct Method The first approach for determining the lateral stiffness of a yarn would be to directly measure the lateral force-lateral displacement of a single yarn – or a stack of yarns – using a compression test as described in [9]. To do so, an experimental set-up was prepared. An aluminum fixture was designed and fabricated to transform the tensile force provided by a KES-G1 tensile machine (Kato Tech, Kyoto, Japan) to a compressive force. Moreover, a 69  punch, made of brass, was utilized to apply a compressive force to the yarns (Figure 4-12). Two loose tapes were attached to the yarn stack to mimic the effect of yarn interlacement. A compressive force was applied to the stack of yarns and subsequently, the force-displacement curve for a single yarn was extracted based on the procedure given in section 3.3.2 of [9]. Figure 4-13 shows the force-displacement curve obtained for a single yarn. As can be seen, the graph shows two distinct trends, first of which is related to the elastoplastic behavior of the yarn in which the filaments slide over each other. The amount of energy dissipated by this mechanism is negligible compared to the second part of the response in which the filaments are laterally under compression. Despite the fact that the second part of the graph has a non-linear behavior with fluctuations, the R-squared (coefficient of determination) of the fitted green line was found to be 0.93. Hence, for the tested fabric, by considering an equivalent linear behavior for the compression of yarn, the slope of the fitted line (10.6 N/mm) was used to estimate K.  Figure 4-12) Test set-up used to measure the lateral stiffness of the yarn. 70   Figure 4-13) The force-displacement curve for a single yarn under compression in the width direction (Force-displacement curve without any yarn in the fixture was first drawn to determine the resistance of the fixture; then this resistance was subtracted from force-displacement curve of the yarn under lateral compression and the net value is shown here).  4.4.2.2.2 Inverse Method Another practical approach for determining K is to use an inverse method based on a bias extension test. Launay et al. [29] introduced an equation that relates the tensile load to the shear load per unit length (Fshear). By measuring the BET shear load per unit length, the increment of the work done on the fabric element (dWshear) through shear stress can be obtained as: 𝑑𝑊𝑠ℎ𝑒𝑎𝑟 = (𝑑𝑉). 𝜏𝑠ℎ𝑒𝑎𝑟 . (𝑑𝛾) (4-13)  where dV is the volume of fabric element. Therefore, the work done on this element from the onset of locking to wrinkling is given by: (𝑊𝑠ℎ𝑒𝑎𝑟)𝐿𝑜𝑐𝑘𝑖𝑛𝑔−𝑡𝑜−𝑤𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔 = 𝑆2𝑡 ∫ 𝜏𝑠ℎ𝑒𝑎𝑟𝛾𝑤𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔𝛾𝑙𝑜𝑐𝑘𝑖𝑛𝑔𝑑𝛾                                                   = 𝑆2 ∫ 𝐹𝑠ℎ𝑒𝑎𝑟 . 𝑑𝛾𝛾𝑊𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔 (4-14) 71  If we assume that the shear work from locking to wrinkling mainly causes the transverse compaction of the yarns in the fabric element, for the PWCPE structure we can write: (𝑊𝑠ℎ𝑒𝑎𝑟)𝐿𝑜𝑐𝑘𝑖𝑛𝑔−𝑡𝑜−𝑤𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔 = 4 × (12𝐾∆2) (4-15)  The value of Δ at wrinkling onset is ∆𝑐𝑟 and can be determined from equation (4-12). Then, 𝐾 = 12 (𝑤0 + 𝑡 − 𝑙 𝑐𝑜𝑠(𝛾𝑤𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔))2 × ((𝑊𝑠ℎ𝑒𝑎𝑟)𝑙𝑜𝑐𝑘𝑖𝑛𝑔−𝑡𝑜−𝑤𝑟𝑖𝑛𝑘𝑙𝑖𝑛𝑔) (4-16)  where (Wshear)locking-to-wrinkling is obtained from the area under force-displacement of the BET. The value obtained for K using this method is 11.53 N/mm. Although the value of K extracted using the inverse method is comparable to that of the direct measurement approach, the characterization challenge in the inverse method is deemed to be substantially less, though more subjective. To minimize modeling error, an average value between the direct and inverse approaches was eventually used for the final lateral compaction stiffness of yarn (Table 4-1). In addition to these two approaches, a micro-level FE simulation may be considered in future studies as a viable alternative. Table 4-1) Summary of the measured geometric and mechanical characteristics of the PWCP. Yarn Material  Geometric Characteristics  Mechanical Characteristics w0 (mm) S (mm) l (mm) t (mm) Flexural Rigidity - 𝑄𝑏 (N.mm2) Yarn Lateral Stiffness - K (N/mm) Carbon fiber 1.4±0.04 4.1 2.05 0.355±0.04 2.75±0.21 11.06 72  4.4.3 Bias Extension Test for Validating the Analytical Wrinkling Model  By substituting the geometric and mechanical properties of the fabric – presented in Table 4-1 – into equations (3-9) and (4-12), the locking and critical (wrinkling) shear angles of the woven fabric was determined to be 31.1o and 40.9o, respectively. To evaluate the validity of the proposed equations, a BET was carried out (Figure 4-14). Three repeats were conducted, and in each repeat, the shear angle of the fabric was recorded using a digital camera mounted in front of the test set-up. Afterward, the pictures were exported to Autodesk Inventor CAD software and the shear angle of the central fabric elements was determined. Figure 4-15 shows the values for the measured shear angles corresponding to locking and wrinkling versus those calculated from the analytical approach. The experimental results revealed an acceptable agreement with the analytical predictions (i.e., the difference as a percentage of the average (i.e., |𝑉𝑎𝑙𝑢𝑒 𝐼 − 𝑉𝑎𝑙𝑢𝑒 𝐼𝐼|(𝑉𝑎𝑙𝑢𝑒 𝐼+𝑉𝑎𝑙𝑢𝑒 𝐼𝐼2)) between the experimental result of wrinkling and analytical prediction is about 7.35%). This difference could stem from a number of potential sources, such as the assumptions made in deriving equivalent stiffness of torsion springs (Kt), linearity assumption for yarn behavior under compression in width direction (K), errors in mechanical and geometric characterization of the fabric, the lack of a complete uniformity in the fabric properties, etc. It should be highlighted that the mentioned potential errors in the analytical model would sometimes mask the negative effects of one another and make the comparison between analytical predictions and experimental results very difficult. The issue is inherent to all analytical models.   73          Figure 4-14) The PWCP (a) before shear deformation, (b) at the onset of locking (𝛾 =28.7±0.8 ̊), (c) at the onset of wrinkling (𝛾 =38±0.5 ̊).   Figure 4-15) Experimental measurements of locking and wrinkling angles versus the predictions of the proposed analytical model.  (a) (b) (c) 74  4.5 Summary A mesoscopic study on wrinkling of the PWCPs under trellising shear mode was presented based on an instability analysis point of view. By considering the mechanism of wrinkle formation in mesoscale, an equivalent structure mimicking the behavior of a representative PWCP element in the post-locking stage was proposed. The stability of the equivalent structure was assessed using the principle of minimum total potential energy. The onset of instability in the equivalent structure, which is coincident with the onset of wrinkling in the PWCP, was analytically determined. The proposed equations (4-11) and (4-12) imply that a number of parameters influence the extent of the critical shear angle of the PWCPs: the lateral compaction stiffness of yarns (K), tension along yarns (T), effective flexural rigidity of yarns (𝑄𝑏), initial yarn width (w0), initial yarn thickness (t), the yarn length in the representative PWCP element (S), and the initial distance between the central axes of the adjacent parallel yarns (l). The analytical model introduced in this paper may be helpful for the following applications. Fabric Design for Superior Wrinkling Behavior: Occurrence of wrinkling in the course of draping of the PWCPs is currently a major concern and can be of the types illustrated in Figure 4-1. Among different types, wrinkling due to large in-plane shear deformation is deemed critical in draping around parts with sharp corners. The presented model may be considered as a first step toward better analytically comprehending the mechanics of wrinkle formation in the PWCPs due to large shear deformation. One of the main merits of the approach is incorporating both geometric and mechanical characteristics of the fabrics. Since the model comprises of mesoscale parameters of the PWCPs, it opens a window toward analyzing the impact of sub-scale 75  contributing parameters on the shear wrinkling behavior of the woven composite preforms. More specifically, upon a parametric study of the mesoscale specifications of different fabrics, optimum PWCPs can be analytically designed/selected to offer a superior resistance to shear wrinkling. Minimum Tensioning of Yarns during Forming Process: Determining the minimum tension along yarns to avoid wrinkling is one of the most practical challenges [55], and the proposed model has the potential to determine the minimum tension. Equation (4-12) explicitly reveals how the critical shear angle can be affected by adjusting tension along yarns. Referring to equation (4-3), it is evident that A is increased if tension along yarns T is ramped up. By taking equation (4-11) into account, Δcr becomes larger by the growth of T, inferring that wrinkling can be postponed or even entirely avoided. Substituting the values of the parameters of a given PWCP and the rough approximation of the expected maximum shear angle into equation (4-12), e.g. using available kinematics based draping software, the minimum tension required along yarns to avoid shear wrinkling is obtainable. This cost-effective solution would be useful for the industrial optimization of forming processes of large PWCP parts. In the current study, a linear behavior of yarns under tension was presumed, while for more precise predictions, a model based on a non-linear behavior of the yarns should be developed.    76  5 Chapter Five Prediction of Wrinkling Initiation in PWCPs under Pure Multi-Scale Shear Deformation Mode 5.1 Overview As discussed earlier in the previous chapters, shear deformation of woven fabrics occurs through two distinct deformation modes, namely trellising and pure multi-scale shear mode. A mesoscopic analytical model to predict the onset of wrinkling of PWCPs under trellising shear mode was presented in Chapter 4. However, this model is not capable of predicting wrinkling initiation in PWCPs under pure multi-scale shear deformation mode owing to the distinct mesoscale deformation mechanisms. Hence, this chapter aims to introduce a modified analytical criterion for the prediction of the onset of wrinkling of PWCPs under pure multi-scale shear deformation mode. To this end, the mechanism of wrinkle formation of PWCPs under pure multi-scale shear deformation mode is first discussed. Afterward, a modified analytical model to predict wrinkling initiation is proposed. Furthermore, a detailed discussion is presented to explain the reason of asynchronous onsets of wrinkling of PWCPs under bias extension and picture frame tests. The discussion challenges the currently agreed justification in the literature. 77  5.2 Mechanism of Wrinkle Formation under Pure Multi-Scale Shear Deformation Mode There exists an inherent difference between trellising mode and pure multi-scale shear deformation mode which stems from different boundary conditions applied to the ends of yarns. If we revisit Figure 3-14, it is evident that as shear deformation of the fabric element increases, the instantaneous yarn width (w) along with the instantaneous length of the element (L) decreases. This causes reduction of the spacing between the adjacent parallel yarns (Figure 3-15). The reduction in gaps between the yarns continues up to the locking point. After locking initiates, the yarn compaction starts and increases by the progress of the shear deformation of the fabric element. This phenomenon takes place because the reduction rate of the instantaneous length of the element (L) over shear deformation is larger than that of the width of the two yarns engaged (2w). This fact is demonstrated in Figure 5-1, in which the derivatives of the instantaneous length of the element, 𝑑𝑑𝛾𝐿0. cos(𝛾), and the instantaneous width of the two yarns, 𝑑𝑑𝛾2(𝑤0 cos(𝛾)), with respect to shear angle (𝛾) are plotted. Similar to trellising, as the shear deformation of the fabric element under pure multi-scale shear mode continues, the contact forces between the adjacent parallel yarns grow to a certain limit termed the critical compressive force that causes wrinkling of the fabric. Figure 5-2 illustrates the post-locking stage of the adjacent parallel yarns of the fabric element under pure multi-scale shear deformation mode. These compressive forces press the yarns in the width direction and provide more room for shear deformation of the element. Hence, by having the value of the reduction in the yarn width just before wrinkling, the corresponding critical shear angle (i.e., the onset 78  of shear wrinkling) under pure multi-scale shear mode can be determined as shown in the next section.     Figure 5-1) Reduction rate of the instantaneous width of the fabric element (L) and the total instantaneous width of two parallel adjacent yarns (2w) under pure multi-scale shear mode.   Figure 5-2) The contact forces due to large pure multi-scale shear deformation exerted on (a) the weft yarns, (b) the warp yarns (the black line between the adjacent parallel yarns represents the interlacing yarn from top view). (b) (a) 10 30 50 70 - 3 - 2 - 1 Shear angle of the fabric element (degree) Instantaneous width of fabric element (L) Instantaneous width of two yarns (2w) The reduction rate (mm/rad) 0 79  5.3 Modified Analytical Model for the Prediction of Wrinkling under Pure Multi-Scale Shear Deformation Mode We know that at the “onset of wrinkling”, the fabric deformation pattern changes from in-plane mode to out-of-plane mode. The four-link equivalent structure introduced in Chapter 4, which mimics the behavior of the representative PWCPs in the post-locking stage, is employed in this analysis with some modifications. In trellising deformation mode, the intra-yarn shear was ignored, inferring that the width of the yarns remains unchanged prior to locking. However, due to the boundary condition of the ends of the yarns in pure multi-scale shear deformation mode, the yarn width decreases as a result of the intra-yarn shear. In this section, the four-link equivalent structure introduced in Chapter 4 is modified to take the yarn width change into account. The modified structure is comprised of four bars of length w (the instantaneous width of the yarns by considering the intra-yarn shear). The bars are clamped elastically at one end against rotation – about the dash-lined axes shown in Figure 5-3 – using four torsion springs of stiffness Kt, and linked to ball-joints at the other end (Point O). These torsion springs represent the equivalent resistance of the fabric element against out-of-plane deformation. Four sliding springs with stiffness K are placed around the bars to capture the lateral compaction resistance of the yarns. As the shear deformation of the element proceeds, the neutral length of the sliding springs is updated to match the length of the bars. The bars are treated as the core of the sliding springs. The quasi-static equivalent compressive force P (i.e., the distributed contact force illustrated in Figure 5-2 times by L0 shown in Figure 3-14) acts on the system as illustrated in Figure 5-3. This force represents the compaction force between adjacent parallel yarns in the fabric 80  element. When the compressive force P reaches a critical value (Pcr)pure multi-scale shear mode or (Pcr)PMSM, the equivalent structure becomes unstable. At the onset of instability, the deformation mode of the equivalent structure changes from in-plane to out-of-plane.    Figure 5-3) (a) A fabric element deformed and locked under pure multi-scale shear deformation, (b) an equivalent representative structure mimicking the behavior of the fabric element in the post-locking stage.  All the calculations in Chapter 4 for the prediction of critical compressive force under trellising mode need to be repeated with the difference that the yarn width in the calculation is no longer w0; instead 𝑤 = 𝑤0 cos 𝛾 which is the instantaneous yarn width should be used. By updating all the equations, equation (5-1) is derived which shows the critical compressive force by taking the intra-yarn shear into consideration:  (𝑃𝑐𝑟)𝑃𝑀𝑆𝑀 =𝐴𝐾𝑤0 cos((𝛾𝑐𝑟)𝑃𝑀𝑆𝑀)𝐴 + 𝐾 (5-1)  81  where K is the stiffness coefficient of the sliding springs, (𝛾𝑐𝑟)𝑃𝑀𝑆𝑀 is the shear angle of the fabric element under pure multi-scale shear deformation mode at the onset of wrinkling and A is a lumped parameter notation as follows: 𝐴 =  2𝑄𝑏 (𝑇𝑄𝑏)32( 𝑆√𝑇𝑄𝑏2) − tanh( 𝑆√𝑇𝑄𝑏2)  (5-2)  where T, 𝑄𝑏, and S are the tension along the yarn, the flexural rigidity of yarn, and the length of the yarn in the fabric element, respectively. By calculating the critical compressive force between the adjacent parallel yarns, the critical compaction of yarns in the width direction just before the wrinkling (∆𝑐𝑟)PMSM can be determined as:  (𝛥𝑐𝑟)𝑃𝑀𝑆𝑀 =(𝑃𝑐𝑟)𝑃𝑀𝑆𝑀𝐾=𝐴𝑤0cos (𝛾𝑐𝑟)𝑃𝑀𝑆𝑀𝐴 + 𝐾 (5-3)  To determine locking angle of the fabric element in pure multi-scale shear deformation mode, the following geometric equation was derived in Chapter 3:  𝑆𝑃 = 𝐿 − 2𝑤 = 𝐿0. cos(𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔) − 2(𝑤0 cos(𝛾𝐿𝑜𝑐𝑘𝑖𝑛𝑔)) = 𝑡 (5-4)  We know that after locking, the yarns will undergo lateral compaction until wrinkling sets in. Hence, the mathematical interpretation of the conditions at which wrinkling begins under pure multi-scale shear deformation mode is:  𝐿 − 2(𝑤 − (∆𝑐𝑟)𝑃𝑀𝑆𝑀) =                                = 𝐿0. cos((𝛾𝑐𝑟)𝑃𝑀𝑆𝑀) − 2(𝑤0 cos((𝛾𝑐𝑟)𝑃𝑀𝑆𝑀) − (∆𝑐𝑟)𝑃𝑀𝑆𝑀) = 𝑡 (5-5) 82  Finally, the critical shear angle can be determined by substituting (∆𝑐𝑟)𝑃𝑀𝑆𝑀 obtained from equation (5-3) into equation (5-5) and solving for (𝛾𝑐𝑟)𝑃𝑀𝑆𝑀: (𝛾𝑐𝑟)𝑃𝑀𝑆𝑀 = 𝑐𝑜𝑠−1 (𝑡𝐿0 − 2𝑤0 +2𝐴𝑤0𝐴 + 𝐾) (5-6)  It should be noted that in this analytical model, the compaction of yarns in the thickness direction is ignored since its value for a typical fabric is much smaller than that of the width direction. In order to see how significantly different are the onsets of wrinkling of an identical woven fabric under trellising mode and pure multi-scale shear deformation mode, the analytical predictions are compared in Figure 5-4. As can be seen in the figure, the difference as a percentage of the average (i.e., |𝑉𝑎𝑙𝑢𝑒 𝐼 − 𝑉𝑎𝑙𝑢𝑒 𝐼𝐼|(𝑉𝑎𝑙𝑢𝑒 𝐼+𝑉𝑎𝑙𝑢𝑒 𝐼𝐼2)) between the onsets of wrinkling is about 52.2%, which is substantially large.   Figure 5-4) A comparison of the predictions of wrinkling angle (i.e., the shear angle at which wrinkling initiates) by relying on trellising shear mode and pure multi-scale shear mode. The width and thickness of the yarns and the spacing between the adjacent parallel yarns are assumed 1.4 mm, 0.355 mm and 0.65 mm, respectively. Trellising Mode Pure Multi-Scale Shear Mode    83  This difference has been experimentally observed and reported in the literature by comparing the onsets of wrinkling of an identical fabric but under two different tests, namely bias extension test (i.e., whose deformation mode is trellising) and picture frame test (i.e., whose deformation mode is pure multi-scale shear). A large portion of the literature has been devoted to this question that why the onsets of wrinkling in these two tests are asynchronous. In the following section, the rationale behind this observation is discussed. 5.4 Why Are the Onsets of Wrinkling in Bias Extension Test and Picture Frame Test Asynchronous? 5.4.1 Common Understanding in the Literature In the past studies, the only deformation mode to describe the shear behavior of woven fabrics was trellising. Based on this assumption, the asynchronous onsets of wrinkling of an identical fabric under BET and PFT were really puzzling. To justify this observation, a number of research was undertaken, in the majority of which tension in yarns during picture frame test was introduced as the source of this phenomenon. As an example, Lanuany et al. [29] justified the distinct response of a woven fabric in the PFT and BET tests – asynchronous onsets of wrinkling – based on the existence of spurious tension in the yarns during picture frame testing. Namely, they explained that tension along the yarns, which would be a result of misalignment of the yarns, can noticeably delay the onset of wrinkling in the PFT. Hivet and Duong further tackled this issue using the same assumption and correlated different responses of the fabric to the tension induced in the yarns. However, the source of tension was not deemed to be misalignment, rather it was assumed to be due to the increase in the crimping of yarns [75]. It has been shown both numerically and analytically (using 84  continuum mechanics) that the crimping of yarns rises over the shear deformation [45]. Experimentally, this issue was observed in the BET in which the yarns are free at their ends [7]. However, the same observation was not made in the PFT owing to the boundary conditions imposed on the ends of yarns. In simple words, on one hand, the effective length of the yarns – the length of the projection of the yarn onto a plane that contains the fabric – over shear deformation in the PFT must remain constant due to the frame kinematics [5]; on the other hand, the growth in crimping tends to reduce the effective length. Accordingly, some tension has to be produced in the yarns under the PFT to maintain the effective length unchanged. By this point, they concluded that in the PFT some tension is generated in the yarns and it delays the wrinkling of the woven preforms. 5.4.2 New Understanding Although the interpretations given in the literature could be correct, an important fact that the deformation modes of the fabrics under these two tests are different was not taken into account. Based on the hypothesis proposed in Chapter 3 (i.e., distinct deformation modes), the rate of vanishing the gaps between the adjacent parallel yarns due to shear deformation is different in woven fabrics under trellising mode and pure multi-scale shear deformation mode. Figure 5-5 shows the trajectory of spacing between the adjacent parallel yarns in the trellising and pure multi-scale shear modes. It is seen that locking in the fabric under trellising mode occurs far ahead of the one under pure multi-scale shear mode. Since intra-yarn shear delays locking in the latter case, wrinkling is also automatically delayed.  85     Figure 5-5) Analytical representation of the yarn spacing of the fabric element under trellising and pure multi-scale shear deformation modes.  Another reason would be hidden in the trajectory of the yarn width under trellising and pure multi-scale shear modes. Let’s define a local coordinate system in a yarn whose axes-1, 2 and 3 are in the longitudinal direction, width direction and thickness direction of the yarn, respectively. As Figure 5-6 shows, the trajectory of the yarn width (i.e., the deformation of the yarn in direction 2 of the local coordinate system) in trellising mode and pure multi-scale shear mode is not identical due to intra-yarn shear induced in the latter case. In trellising, since there would be no intra-yarn shear, the yarn width remains unchanged prior to locking (the stage between points (1) and (2) in Figure 5-6). After locking, however, the yarns become compressed in direction 2 due to compressive contact forces and this continues until wrinkles appear in the fabric (the stage between points (2) and (3) in Figure 5-6). On the contrary, the yarn width in pure multi-scale shear mode begins to change upon the shear deformation, and it continues up to locking (the stage between points (1) and (4)). In the pre-locking stage, the yarn width changes according to 𝑤 = 𝑤0𝑐𝑜𝑠𝛾 and it would continue until wrinkles appear (the Trellising mode Pure multi-scale mode Locking angle under trellising mode Locking angle under pure multi-scale mode Weft yarn thickness 0.0 0.2 0.4 0.6 Spacing between the warp yarns (mm) 10 20 30 40 50 60 Shear angle of the fabric element (degree) 70 0.8 86  stage between points (4) and (5) shown by the red dashed-line). However, as locking begins, the compaction of the yarns due to compressive contact forces is initiated, and hence the reduction in the yarn width – in direction 2 – becomes faster (the stage between points (4) and (5) shown by green solid-line in Figure 5-6). Eventually, the fabric structure wrinkles due to the critical compressive contact forces. Figure 5-6 indicates that the span of locking-to-wrinkling in the pure multi-scale shear mode is larger than that of the trellising, theoretically owing to the smaller rate in increase of the compressive forces between the adjacent parallel yarns.   Figure 5-6) Analytical representation of the yarn width of the fabric element under trellising and pure multi-scale shear deformation mode.  In trellising, the shear deformation of the fabric element in the post-locking stage is not possible except by compaction of the yarns in direction 2. This causes a sharp increase in the compressive contact forces between the adjacent parallel yarns that eventually lead to instability of the element. By contrast, the rate of increase in the compressive forces between 0.2 0.6 1.0 1.4 The yarn width (mm) 10 20 30 40 50 60 Shear angle of the fabric element (degree) 70      1 2 3 4 5 Trellising mode Pure multi-scale mode 87  the yarns in the pure multi-scale shear mode is substantially less since the shear deformation of the fabric element in post-locking stage becomes possible not only by the compaction of the yarns in the width direction – direction 2 – but also by the intra-yarn shear that reduces the width of the yarns. As a summary, there exist two reasons for the delayed wrinkling of woven fabrics under PFT (i.e., whose deformation mode is pure multi-scale shear) to BET (i.e., whose deformation mode is trellising): I- Locking, which is the most important trigger in shear wrinkling, occurs far later in the fabric under PFT compared to BET. II- Compressive forces induced between the adjacent parallel yarns after locking is less in the fabric under picture frame test compared to bias extension test. To assess the effectiveness of the proposed hypothesis and the analytical model developed for the prediction of the onset of wrinkling of PWCPs under pure multi-scale shear deformation mode, an experimental investigation is carried out in which a woven fabric is tested by the picture frame test. The results obtained from PFT are compared with those obtained from the bias extension test. 5.5 Experimental Validation In spite of the simplicity in conducting the BET, the picture frame testing whose deformation mode is pure multi-scale shear requires more sophisticated equipment, in particular, if the yarn tensioning needs to be controlled/monitored. A customized biaxial tensile-picture frame fixture was designed and fabricated and was vertically installed on the 88  Instron machine. The fixture was comprised of four steel bars with equal length (23.5 cm) joint to each other with pins, as demonstrated in Figure 5-7. The tensile force was applied along the diagonal of the square frame, generating shear deformation. In addition to the Instron motor to provide the shear loading, four servomotors were perpendicularly attached to each of the frame arms to provide tensile/compressive loading on the warp and weft yarns of the fabric. The load cell data were transferred to data acquisition system, and the LABVIEW was used to control and maintain the tension in the yarns at zero level (i.e., to ensure a pure shear condition).  Figure 5-7) The picture frame test set-up used in the experiments.  89  Using appropriate boundary conditions to generate a uniform deformation field not only in the region of interest but also close to the grips, is regarded as one of the main challenges of the picture frame testing. Non-uniformity in the deformation field, which is a result of e.g., misalignment during installation, can cause fabric samples to experience tension [29], leading to a discrepancy in the results. To cope with this challenge, needles were used to facilitate the yarns settlement at the grips, resulting in more uniform deformation fields [68][76][77]. Haghi Kashani et al. [68] introduced two jaws with needles and slight pressure as a modified boundary condition, which is a combination of a perfectly clamped boundary condition and a free boundary condition. The implementation of the proposed boundary condition indicated that tension in the yarns during the shear deformation is notably less than the simply-clamped boundary condition [13]. Accordingly, in the current study, this boundary condition was adopted (Figure 5-8 (a)). The width of each strip of the sample was set to be 12 cm. Although the effective distance between the grips is 23.5 cm, the length of the samples was chosen 30 cm to provide some room for folding the ends of the sample within the grips and creating good hold between the fabrics and the grips (Figure 5-8 (b)). It has been recommended [78] that the transverse yarns in the arm regions of the PFT sample (see also Figure 3-7) should be removed; otherwise, wrinkling would initiate in the arm regions rather than in the region of interest. The same observation was made in this study as shown in Figure 5-9 (a). This observation can be justified by considering the mesoscale deformation mechanism participating in different regions of a PFT sample. To explain, the mesoscale deformation field in the region of interest follows the pure shear deformation, whereas, in the arm regions with one family of yarns unclamped, the deformation pattern would be something between trellising and pure shear. As a result, locking and wrinkling appear in 90  the arm regions earlier than in the region of interest. To alleviate the chance of wrinkling in the arm regions, most of the transverse yarns were removed from the sample as illustrated in Figure 5-9 (b).  Figure 5-9 (c) shows the latter sample set-up at 64o shear, while no wrinkling initiated neither in the arm regions nor in the region of interest.  Figure 5-8) (a) Using needles in the PFT boundary conditions, (b) a folded sample in the grips to increase the interaction between the fabric ends and the needles.   Figure 5-9) (a) Wrinkling initiation in the PFT arm regions due to the presence of transverse yarns, (b) removing most of the transverse yarns, (c) the shear deformation at 64o with no wrinkles.  91  Despite using the needles at the fabric boundary conditions, a slight increase in yarns tension was recorded by the load cells during the test (up to ~40 N or 0.8 N per yarn). As addressed in Section 5.4.1, this increase in tension can be due to the increase in the crimping of yarns during the shear deformation [75] which reduces the effective length of the yarns. Accordingly, to remove the induced tension, the LABVIEW code was programmed in a way that once the load cells show nonzero values, it turns on the biaxial motors in the compressive direction for a few seconds, keeping the tension within the sample at zero. Doing so eliminated the possibility of the tension along yarns. Figure 5-10 (b) shows the experimentally measured locking angle in the picture frame test versus analytically predicted values. Due to the limitation of the device, the shear angle of the frame could not exceed 64o. Wrinkling of the sample was not observed up to this angle, and this is consistent with the analytically predicted value of 69o through equation (5-6). The analytical prediction for the locking angle in the PFT (i.e., presented in Chapter 3) was also reasonably close to that of experimentally measured angle in Figure 5-10. Hence, despite the fact that tension in yarns was kept at zero, the onsets of wrinkling of woven fabrics under bias extension (i.e., whose deformation mode is trellising) and picture frame test (i.e., whose deformation mode is pure multi-scale shear) were asynchronous, proving the hypothesis proposed in this study. 92   Figure 5-10) Experimental measurements of locking and wrinkling angles versus the predictions proposed by the analytical model.  5.6 Summary The mesoscopic analytical model for the prediction of the onset of wrinkling of PWCPs under trellising shear deformation mode was modified in this chapter to be used for the analysis of wrinkling due to pure multi-scale shear deformation. To do so, the effect of intra-yarn shear, which impacts the yarn width, was incorporated to the model. The analytical results revealed that the onsets of wrinkling of an identical woven fabric under trellising and pure multi-scale shear deformation modes could be substantially different. This conclusion opened up a window toward a better explanation for the question that why the onsets of wrinkling in BET and PFT are asynchronous. For a long time, the asynchronous onset of wrinkling in the BET and PFT attracted much attention in the literature. The different response of the fabric in these two tests, particularly the onset wrinkling, has been dominantly explained by the assumption of the 93  existence of some tension in the yarns during the PFT. In this study, however, it was demonstrated that distinct deformation modes during BET and PFT are the main source of discrepancy in the response of woven fabrics to wrinkling phenomenon. It was experimentally shown that the onset of wrinkling in the BET and PFT could be as different as at least 46%, even when the tension in yarns is kept at zero through the test.                94  6 Chapter Six Potential Applications of the Analytical Model of Wrinkling Prediction 6.1 Overview Analytical wrinkling criteria have always been a powerful tool to either be used in finite element packages as a wrinkling indicator (predictor) or identify the contributing parameters in the wrinkling phenomenon. In the field of sheet metal forming, the efforts toward employing the analytical wrinkling criteria for the abovementioned purposes are numerous; however, the same trend cannot be found in the literature of fabric forming. Hence, it would be advantageous to examine the capability of the analytical model developed in this thesis for predicting wrinkling initiation in FE packages and determining the effect of contributing parameters on this phenomenon. 6.2 Implementation of the Analytical Model in ABAQUS FE Package Per the discussion presented in Chapter 2, continuum modeling of woven fabrics is currently one of the most convenient and common strategies for the modeling of the fibrous materials thanks to the ease of implementation into the finite element packages and minimization of the computational cost. Almost all the attention of the academia and the industry in the finite element modeling of fabric forming processes was devoted only to precise constitutive equations which closely model the behavior of a deforming fabric. However, they did not consider an important fact that how wrinkles initiate in an FE model. Ignoring this 95  question resulted in discrepancies between the finite element predictions of wrinkling and the experimental results. To assess how vital is this question in the fabric forming process, a simple problem – bias extension test – is modeled in ABAQUS FE package. To do so, having an accurate constitutive equation for defining the behavior of plain-woven fabrics is of great importance. In another study conducted in parallel to the current thesis (i.e., the author was part of that study and the details of the fabric modeling can be found in [79]), the behavior of a plain-woven fabric was modeled using a hypoelastic constitutive model in which the non-orthogonality and coupling effects of the fabric (i.e., coupling means that the effective material properties of a given direction of the fabric change due to the applied deformation in other directions) were taken into account. The hypoelastic approach was put into practice since the experimental investigation carried out on the fabric revealed that the material behavior is notably dependent on load history.  After defining the constitutive equations of the woven fabric, the material model is implemented in ABAQUS FE package via User Material (UMAT) subroutine (please refer to the appendix for more information regarding the UMAT subroutine). Four-node shell elements (S4) are used in this simulation. To validate the accurate implementation of the defined material model in UMAT, a number of loading scenarios are simulated in the software [79]. Then, the FE results are compared with the results of the experiments and those of the constitutive model directly drawn by MATLAB. Figure 6-1 shows two examples of the examined loading scenarios to validate the accurate implementation of the constitutive equations. As the figure indicates, the FE results are reasonably close to those obtained from the experiments and the constitutive model directly drawn by MATLAB, proving the accurate implementation of the constitutive equations into ABAQUS. 96        Figure 6-1) Validation of the UMAT code implemented in ABAQUS, (a) accurate prediction of longitudinal stress under sequential tensile test with 2% transverse pre-tension, (b) precise anticipation of transverse stress under sequential tensile test with 2% transverse pre-tension. These graphs highlight that how significant is the effect of coupling on the fabric behavior.  After ensuring that the material model is correctly implemented in ABAQUS, the bias extension test is simulated. In spite of large shear angles of the fabric in the middle of the sample (i.e., 70̊), the finite element model does not show any wrinkle in region III of the sample (for definition of region III, refer to Figure 3-6). The most important reason of observing such phenomenon is that since there is no contact in this example and also implicit solver is used, 0 0.01 0.02 0.03 Longitudinal strain 0.04 5 10 15 20 Longitudinal stress (MPa) 25 (a) Experiment Uncoupled model Coupled model ABAQUS 0 0.01 0.02 0.03 0.04 2.5 3.5 4.5 5.5 Longitudinal strain Transverse stress (MPa) (b) Experiment Uncoupled model Coupled model ABAQUS 97  the finite element software could accurately solve the problem. This infers that the numerical errors are too small to act as imperfections in the model. On the contrary, the observations made during the bias extension test as reported in Chapter 4 unveil another fact. In the sense that the fabric wrinkles when the shear angle exceeds 38̊. This evidently implies that even having an accurate material model does not guarantee accurate prediction of wrinkling in the finite element software. However, employing a constitutive model, which takes the coupling and non-orthogonality into account, is considerably essential to accurately predict the stress field, and in turn the level of tension which yarns experience over the forming process. The calculated tension level will be further used in the wrinkling prediction, as tension is one of the significant parameters in the derived equations for wrinkling in Chapters 4 and 5. To deal with the issue discussed above (i.e., the lack of ability of the FE model to capture wrinkles in bias extension simulation), the finite element code was equipped/integrated with the analytical model developed in Chapter 4. To this end, the flowchart presented in Figure 6-2 is used. The whole idea of the flowchart is to detect the elements in the model which satisfy the condition of wrinkling criterion. To benefit from the analytical model, the mesoscopic geometric and mechanical characteristics of the fabric coupled with the stresses in yarn direction are needed. According to the flowchart, the fabric characteristics are first defined and then, simulation is started. The incremental deformation is applied to the fabric and in each increment, critical shear angle of all elements – the shear angle at which wrinkling initiates – are calculated in terms of material characteristics and the stress field (i.e., tension level).   98                           Figure 6-2) The flowchart based on which wrinkling initiation in the finite element code is detected. Apply incremental deformation Start FE analysis  Determine stress-strain fields using the developed coupled non-orthogonal hypoelastic constitutive model for each element.  Use the analytical model developed in Chapter 4 to obtain the critical shear angle of each element based on the mesoscopic properties of the fabric and stress field of the element.  Is the shear strain of the element larger than the calculated critical value?  Apply an out of plane displacement to the center node of the element. Wrinkling is detected in the model.  Read the mesoscale geometric and mechanical properties of the fabric. No 99  Afterward, the shear strain of each element determined by the FE model is compared with its corresponding critical shear angle. If the shear strain of the element exceeds the calculated critical shear angle, then wrinkling is initiated in that element. After locating the element(s) at which wrinkling is about to occur, a geometric imperfection – out of plane nodal displacement – is applied to the element and the simulation continues. It should be highlighted that the proposed flowchart is a fundamental step toward improving the accuracy of the finite element codes and needs to be more elaborated in the future studies. The approach is inspired by the current methods available in the field of sheet metal forming. Figure 6-3 (a) and (b) illustrate the results of finite element modeling of a woven fabric under bias extension test which is equipped with the analytical wrinkling criterion along with a real woven fabric sample, respectively. As illustrated in the figure, the new FE model is capable of predicting wrinkling in region III of the sample which is consistent with the experimental observations. In spite of the discussion presented above, there are a number of studies in the literature that present finite element models capable of simulating shear wrinkles in the fabric forming process without any wrinkling indicator [27][54][62][80]. Figure 6-4 shows two of those models. This observation contradicts the discussion we just presented (i.e., our FE model could not model wrinkling in bias extension test). The answer would be hidden in the fact that in our model, which simulates bias extension test, the implicit solver causes negligible numerical errors (i.e., accurate predictions). While in the models shown in Figure 6-4, there are some contacts (e.g., fabric-die, fabric-punch and fabric-holder) that make the explicit solver yield results with relatively large round-up numerical errors. The numerical errors act as geometric imperfections in the model and cause wrinkling initiation in the FE software. Hence, it could 100  be concluded that the wrinkles in those models would not be accurate and be as a result of some unwanted numerical errors, whereas wrinkles simulated with our proposed method (i.e., wrinkling indicator method) are as a result of the true wrinkling initiation conditions.   Figure 6-3) (a) FE simulation results of the woven fabric equipped with the wrinkling detection flowchart (the red zone shows the region where wrinkling is about to occur), (b) a real sample under bias extension test.                   Figure 6-4) FE simulations of woven fabric forming process, (a) [62], (b) [27].  6.3 Parametric Study of Wrinkling Initiation The parameters affecting wrinkling behavior of composite preforms could be divided into two main groups, namely fabric parameters and process parameters. The former group comprises those parameters which are related to the fabric such as cross-section, bending (a) (b) 101  rigidity and lateral stiffness of yarn, fiber architecture, stitch pattern (i.e., for non-crimp fabrics), etc. The latter group, however, pertains to the parameters at the processing level, such as blank holding system [81][82], tooling surface properties [81][83], blank characteristics [84] and the geometry of the punch [43]. Figure 6-5 presents a summary of the parameters, which have an impact on the forming induced wrinkles of the fabrics.  Figure 6-5) Parameters affecting wrinkling of the fabrics.  On the contrary to the process parameters which attracted a great deal of attention from the academia and the industry, the fabric parameters were practically untouched and there could not be found a major research in the literature, which targets explicitly the effect of fabric parameters on wrinkling of woven fabrics. In this thesis, therefore, we approach the problem from a new aspect and conduct a parametric study at the fabric level to demonstrate how significant are the effect of fabric parameters on wrinkling behavior of woven fabrics. In addition, the effect of one of the most influential process parameters (i.e., blank holder force which causes tension in yarns) on wrinkling is studied. Accordingly, the section is organized Contributing ParametersFabric ParametersYarn Cross-SectionYarn Bending RigidityYarn Lateral StiffnessFiber ArchitectureStitch PatternProcess ParametersBlank HolderTooling Surface PropertiesBlank Holder TypeBlank CharactirsticsBlank OrientationBlank SizeGeometry of the Punch102  in a manner that “Section 6.3.1” deals with fabric parameters; and “Section 6.3.2” investigates the effect of blank holder force on wrinkling phenomenon. 6.3.1 Fabric Parameters Through this section, it is demonstrated that the fabric parameters can notably impact both types of wrinkling, namely shear and non-shear wrinkling. As an example, if we revisit equation (2-1), we note that the bending stiffness of the fabrics plays a pivotal role in non-shear wrinkling. To explain, a fabric with higher bending stiffness requires a larger in-plane compression for buckling. This infers that the capability of the fabric against non-shear wrinkling would be enhanced as the bending stiffness of the fabric increases. Greater bending stiffness of the fabric is achieved either by stiffer yarns or (and) optimized fiber architecture. Hence, it is evident that the fabric parameters have an influence on the performance of the fabrics in non-shear wrinkling. Since non-shear wrinkling is not in the scope of the present research, it is not elaborated further in the thesis. A similar discussion is valid for shear wrinkling. That is, the fabric parameters could be effective in the capability of the fabric in resisting against shear wrinkling. As stated earlier, locking and wrinkling angles are the most important factors in shear wrinkling. If the fabric parameters are designed such that the values of these angles become larger, then the fabric becomes more resistant against shear wrinkling. The analytical equations introduced in this thesis for the prediction of the onset of wrinkling are useful sources for identifying the effect of fabric parameters on shear wrinkling. For instance, equations (4-12) & (5-6) of the present thesis introduce a number of parameters affecting the onset of shear wrinkling. These parameters are lateral stiffness, bending rigidity, the geometry of the cross-section of yarns, 103  and the spacing between the adjacent parallel yarns. In the following sections, the effect of these parameters on the onset of shear wrinkling is qualitatively studied and it will be demonstrated that how shear wrinkling can be affected by changing these parameters.  6.3.1.1 Mesoscopic Mechanical Parameters Referring to equations (4-12) and (5-6), it is evident that bending rigidity of yarns (Qb) and yarn lateral stiffness (K) are two most influential mesoscopic mechanical parameters on the onset of wrinkling of plain-woven fabrics. Figure 6-6 illustrates the qualitative effect of yarn bending rigidity on the critical shear angle of woven fabrics.      Figure 6-6) Qualitative analysis of the effect of yarn bending rigidity on the critical shear angle of the plain-woven fabric under (a) trellising mode, (b) pure multi-scale shear mode.   As the figure shows, by the increase in bending rigidity of yarns, the onset of shear wrinkling is delayed (i.e., larger critical shear angle). In the analytical model presented in Chapter 4, it has been discussed that the interlacing yarn acts as a resistant factor to out of plane deviation of the two adjacent parallel yarns which are under lateral compression. To elaborate the former statement, we revisit Figure 4-7 here (Figure 6-7). As can be seen, at the onset of wrinkling, yarn Y-2 tends to rotate about the edge JJ’ to form a wrinkle, whereas yarn Y-3 2.0 2.4 2.8 3.2 3.6 39 40 41 42 43 Critical shear angle (degree) Bending rigidity of the yarn – Qb (N.mm2) (a) 2.0 2.4 2.8 3.2 3.6 68.5 69.5 70.5 71.5 72.5 Critical shear angle (degree) Bending rigidity of the yarn – Qb (N.mm2) (b) 104  resists against this rotation. Hence, the more bending rigidity of the interlacing yarn (yarn Y-3), the later shear wrinkling initiation in the fabric element.  Figure 6-7) Revisited figure from Chapter 3.  In contrast to bending rigidity of the yarns, lateral stiffness of the yarns has an opposite effect on wrinkling initiation of woven fabrics. The qualitative effect of lateral stiffness of the yarns on the critical shear angle of woven fabrics is illustrated in Figure 6-8. In accordance with the figure, as the lateral stiffness of the yarns increases, the possibility of wrinkling initiation in the fabric ramps up. This would be justified by considering the fact that stiffer yarns – in the width direction – in the fabric cause larger in-plane compression due to shear deformation. To describe, as shear deformation of the fabric increases, locking initiates. After locking point, the more shear deformation leads to the compaction of the yarns in the width direction. Yarns possessing greater lateral stiffness (K) induce more in-plane compression within the fabric, raising the possibility of instability in the fabric structure. From a macroscale point of view, it means that the shear modulus of the fabric after locking becomes larger by having stiffer yarns in width direction, which results in establishing larger in-plane compression. 105                  Figure 6-8) Qualitative analysis of the effect of yarn lateral stiffness on critical shear angle of the plain-woven fabric under (a) trellising mode, (b) pure multi-scale shear mode.  Although the discussion presented above provides us with a qualitative understanding of the effect of mechanical properties of yarns on the onset of wrinkling of woven fabrics, it should be highlighted that all the discussions presented above are based on the assumption that changing one of the mechanical properties of the yarns does not affect the other properties (i.e., the coupling effects between the yarn properties are not taken into account). In terms of practical implementation of the analytical results obtained, this section of the thesis would open up a discussion in the field of textile engineering on how to make a yarn with minimum lateral stiffness while its flexural rigidity is maximum. 6.3.1.2 Mesoscopic Geometric Parameters Besides the mesoscopic mechanical characteristics of woven fabrics, geometric parameters are also of great importance in wrinkling behavior of the fabrics. The geometric parameters are classified into two levels, namely the geometry of the cross-section of the yarns (i.e., width and thickness of yarn) and spacing between the adjacent parallel yarns. Discussion on the effect of width and thickness of yarns on shear wrinkling is very complicated. However, 69 70 71 72 73 Critical shear angle (degree) (b) 8.0 10 12 14 Lateral stiffness of yarn – K (N/mm) 8.0 10 12 14 Critical shear angle (degree) Lateral stiffness of yarn – K (N/mm) (a) 39.5 40.5 41.5 42.5 43.5 106  a parametric study of spacing between the adjacent parallel yarns seems straightforward. Figure 6-9 indicates the qualitative impact of spacing between yarns on wrinkling initiation. As can be seen in the figure, larger gap between yarns – less interlacing density – causes delayed shear wrinkling in woven fabrics. The reason of this observation roots in the fact that the larger spacing between the adjacent parallel yarns results in the later locking initiation. Since locking, which is the most important trigger in shear wrinkling, occurs later, the shear wrinkling is also delayed.         Figure 6-9) Qualitative analysis of the effect of spacing between yarns on critical shear angle of the plain-woven fabric under (a) trellising mode, (b) pure multi-scale shear mode.  An experimental study has been carried out in the literature in which two different fabrics – a loose and a tight plain-woven fabric – are draped over a hemispherical punch [32]. The experimental results reveal that the loose plain-woven fabric is draped over the punch without the formation of wrinkling, while in the tight fabric, there appear several wrinkles (Figure 6-10). This experimental observation is consistent with the parametric study conducted in this chapter. 0.9 66 68 70 72 74 Critical shear angle (degree) (b) Spacing between yarns – SP (mm) 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 Critical shear angle (degree) Spacing between yarns – SP (mm) 34 38 42 46 0.9 (a) 107   Figure 6-10) Draped woven fabric over a hemispherical punch, (a) a loose plain-woven fabric, (b) a tight plain-woven fabric [32].  As stated earlier, the parametric analysis of the cross-sectional geometry of the yarn is very challenging. As an example, referring to equation (3-9), it is evident that the increase in thickness of yarn causes earlier locking of the fabric (Figure 6-11). Since locking is the most influential trigger of shear wrinkling, it is expected that shear wrinkling sets in earlier in the fabric. However, there exists another fact that should not be overlooked. That is, as the thickness of yarn becomes larger, its flexural rigidity increases which causes more resistance against wrinkling initiation. Accordingly, there exists an interaction between the positive and negative effects of yarn thickness on shear wrinkling which makes the parametric study case sensitive. The same story can be mentioned for the effect of yarn width on shear wrinkling. On one hand, the increase in yarn width makes the fabric structure susceptible to earlier instability (i.e., by taking the equations (4-12) and (5-6) into account); on the other hand, it increases flexural rigidity of the yarn. Yet again, both positive and negative interactional effects can be observed in this parameter. 108  In light of the discussion presented here, it is now obvious that the fabric parameters can notably influence the behavior of woven fabrics in wrinkling phenomenon. Hence, it would be really advantageous to efficiently design the fabric parameters in addition to the process parameters. This study unveils the importance of the fabric parameters on wrinkling and calls for more attention from the field of textile engineering to focus more on the mesoscopic parameters of the fabric.     Figure 6-11) Qualitative analysis of the effect of thickness of interlacing yarn on locking angle of the plain-woven fabric under trellising mode.  6.3.2 Process Parameters As shown in Figure 6-5, there are several process parameters which affect wrinkling behavior of woven fabrics, including blank holder force, tooling surface properties, the geometry of the punch, etc. Among these parameters, blank holder force has been identified as the most influential parameter. Therefore, in the following section, the effect of blank holder force on wrinkling is studied. 0.15 0.25 0.35 0.45 Locking angle (degree) Thickness of interlacing yarns – t (mm) 26 30 34 38 42 109  6.3.2.1 Blank Holder Force In draping of composite preforms, Blank Holder Force (BHF) plays a significant role in controlling the formation of both types of wrinkles. During the forming process, as punch travels, the fabric is drawn into the die cavity. Blank holder force, as one of the most important process parameters, causes frictional forces on the interfaces of fabric-die and fabric-blank holder as shown in Figure 6-12. These tangential forces result in tension in the yarns [81] whose magnitude depends on tool surface properties [85], the magnitude of BHF [6], yarns orientation, etc. It has also been experimentally proven that if the tension in the yarns increases, the formation of wrinkles (i.e., shear and non-shear wrinkles) is delayed [40][43]. Hence, the role of BHF in delaying or even preventing wrinkling – by inducing tension in yarns – is evident. Studies on this end are numerous [84][6][81][82][86][87], most of which investigated the effect of this parameter experimentally and some of which dealt with the problem using finite element simulations. However, the analytical analysis of the effect of this parameter is relatively limited.    Figure 6-12) The schematic of a draping process.  In this section of the thesis, we aim to investigate the effect of this parameter on wrinkling using the analytical model developed in Chapters 4 and 5. If we revisit the equations (4-12) and (5-6), which predict in turn the onset of shear wrinkling in plain-woven fabrics under trellising and pure multi-scale shear deformation mode, it is seen that as tension in yarns BHF BHF Punch Die Die Fabric blank Frictional Force Frictional Force 110  increases, wrinkling initiation is notably delayed (Figure 6-13). This observation is a result of the second order analysis conducted in deriving equation (4-2) for determination of Kt (for more discussion, please refer to Section 4.3.2 of the present thesis). Also, we know that tension in yarns is induced as a result of blank holder force. Consequently, it means that larger BHF causes more tension in the yarns and eventually, delayed wrinkling.       Figure 6-13) Qualitative analysis of the effect of tension in yarns on the critical shear angle of the plain-woven fabric under, (a) trellising mode, (b) pure multi-scale shear mode.  This analytical prediction is consistent with the FE results reported in Yu et al. [81]. They demonstrated that wrinkling initiation is markedly delayed as blank holder force increases. Figure 6-14 (a) and (b) illustrate how wrinkles are reduced by increasing the magnitude of the BHF. A more explicit correlation between the magnitude of BHF and wrinkling could be observed in a research conducted by Arnold et al. [6] in which the effect of BHF on the amplitude of wrinkles is studied. Figure 6-15 illustrates that as the magnitude of BHF increases, wrinkle amplitude decreases. This experimental observation is also in agreement with the predictions of our proposed analytical model. 1 3 5 7 Critical shear angle (degree) Tensions in yarns – T (N) 72 74 76 78 80 9 82 (b) Critical shear angle (degree) 65 1 3 5 7 Tension in yarns – T (N) 45 50 55 60 9 (a) 111    Figure 6-14) The effect of BHF on wrinkling of fabrics, (a) BHF = 0, (b) BHF = 5000 N [81].    Figure 6-15) Effect of punch displacement and magnitude of BHF on the amplitude of wrinkles [6].  One of the industrial applications of the proposed analytical model is to recursively determine the minimum blank holder force required to avoid wrinkling for a given geometry. To explain, if one substitutes the values of the parameters of a given fabric and the rough approximation of the expected maximum shear angle – e.g., using available kinematics based 112  draping software – into equation (5-6), the minimum tension required along yarns to avoid shear wrinkling is obtainable. Then, by having the frictional properties of the fabric-tool, the minimum BHF to avoid shear wrinkling could be determined. Determining minimum BHF to avoid wrinkling has always been of interest in the industry since large blank holder force causes fiber breakage in woven fabrics which is considered as a critical defect in the composites industry.  6.4 Summary The analytical model for the prediction of shear wrinkling of woven fabrics was employed in this chapter as a wrinkling predictor in ABAQUS FE package and a useful tool for the parametric analysis of wrinkling phenomenon. To make use of the analytical model in FE software, an algorithm was proposed based on which shear wrinkling of woven fabrics could be accurately predicted. Moreover, it was demonstrated that fabric parameters including cross-section, bending rigidity and lateral stiffness of yarn, along with interlacing density have a notable impact on shear wrinkling. Accordingly, designing these parameters in a manner that cause delayed wrinkling was strongly recommended. Finally, the application of the analytical model to determine a rough approximation of minimum BHF to avoid wrinkling was explained.      113  7 Chapter Seven Conclusions and Future Works Forming induced wrinkles in woven composite preforms have always been a controversial issue in both the academia and the composite processing industry. Accordingly, a large portion of the literature of composite manufacturing has been devoted to the study of this phenomenon – forming induced wrinkle. In spite of the considerable amount of research into forming induced wrinkles in woven fabrics, there was not a general agreement on the mechanism of wrinkle formation. In almost all studies, all wrinkles were treated similarly; however, through this research, we classified wrinkles into two groups: non-shear wrinkle and shear wrinkle. Although the nature of both types of wrinkles is identical (i.e., both wrinkles are developed due to excessive in-plane compression), the mechanisms contributing to the initiation of this phenomenon are entirely different. To put it another way, the in-plane compression induced in the former wrinkling case is due to the longitudinal contraction of the fabric segment, whereas in the latter case, this compression is initiated because of the large intra-ply shear deformation of the fabric (i.e., compression is as a result of locking). In the field of fabric forming, shear wrinkling has been of much interest since the possibility of formation of this type of wrinkling is relatively high in intricate shapes with sharp corners. As such, the present thesis was aimed to specifically study shear wrinkling in woven fabrics. In the literature of shear wrinkling, the lack of an accurate analytical method for the prediction of the onset of this phenomenon was evident. The analytical models are of significant importance since they can be used as either a shear wrinkling indicator in the finite 114  element software for improving the accuracy of the wrinkling predictions and/or a valuable source for parametric study of this phenomenon. Given the discussion presented above, we undertook a fundamental study to develop an analytical model for the prediction of shear wrinkling. To this end, a mesoscopic investigation of the behavior of woven fabric in shear deformation was conducted and it was revealed that shear deformation of woven fabrics could occur through two distinct deformation modes: trellising shear mode (i.e., the mode occurring in bias extension test) and pure multi-scale shear deformation mode (i.e., the mode occurring in picture frame test). For quite a long time, trellising was the only identified deformation mode to describe the behavior of woven fabrics during fabric shearing. However, our detailed analysis of the fabric deformation field in the course of picture frame test revealed that there exists another deformation mode – pure multi-scale shear deformation mode – in the fabric, in which yarn width does not remain unchanged prior to locking as a result of intra-yarn shear. That is, the trajectories of yarn width in trellising and pure multi-scale shear deformation modes are quite distinct. This inherent difference in turn results in asynchronous onsets of locking and eventually wrinkling. Acquiring this understanding sheds light on an area of the literature, which was a source of controversy for years. To explain, the source of discrepancies between the wrinkling behavior of woven fabrics under bias extension test and picture frame test was to some extent ambiguous in the literature. This research revealed that this discrepancy stems from the distinct deformation modes of woven fabrics in the course of shear deformation. Building on this understanding, we studied the mesoscopic nature of shear wrinkling and then, proposed an analytical model to predict the onset of this phenomenon with the aid of the instability analysis of the fabric structure. This approach lies among the first efforts to predict 115  the onset of wrinkling in woven fabrics using the concept of instability. The analytical model developed in this study includes both fabric parameters, such as the lateral compaction stiffness of yarns (K), effective flexural rigidity of yarns (𝑄𝑏), initial yarn width (w0), initial yarn thickness (t), the yarn length in the representative PWCP element (S), and the initial distance between the central axes of the adjacent parallel yarns (l); and process parameters, such as blank holder force. The analytical model implies that geometric characteristics of the fabric along with its mechanical properties can notably impact the onset of wrinkling. For example, an increase in spacing between the adjacent parallel yarns can cause delayed locking and shear wrinkling. Likewise, rising bending rigidity of the yarn can improve the resistance of the fabric against wrinkling. In opposite, having a stiffer yarn in the width direction – direction 2 in local coordinate system define in Section 5.4.2 – diminishes the performance of the fabric in wrinkling resistance. All things considered, it would be advantageous to sometimes design the fabric parameters in a manner to have an improved performance in wrinkling phenomenon. In addition, the analytical model developed in this study can be used to study the effect of yarn tensioning on the onset of wrinkling. Determining the minimum tension along yarns to avoid wrinkling is one of the most practical challenges, and the proposed model has the potential to determine the minimum tension. Equation (4-12) explicitly reveals how the critical shear angle can be affected by adjusting tension along yarns. Referring to equation (4-3), it is evident that A is increased if tension along yarns T is ramped up. By taking equation (4-11) into account, Δcr becomes larger by the growth of T, inferring that wrinkling can be postponed or even entirely avoided. Substituting the values of the parameters of a given PWCP and the rough approximation of the expected maximum shear angle into equation (4-12), e.g. using 116  available kinematics based draping software, the minimum tension required along yarns to avoid shear wrinkling is obtainable. This cost-effective solution would be useful for the industrial optimization of forming processes of large PWCP parts. This approach is among a very few studies that tried to analytically investigate the effect of yarn tensioning on wrinkling phenomenon.  Finally, the analytical model developed in the thesis can be used in finite element packages (e.g., ABAQUS) as a wrinkling indicator. In this study as the first effort in the literature of fabric forming simulation, we integrated the analytical wrinkling indicator in ABAQUS FE package and demonstrated that the accuracy of the software in predicting shear wrinkling could be notably improved. Since the current research is one of a few studies that analytically approached shear wrinkling problem, there still exist a number of potential research initiatives that could be studied in the future: - The analytical model developed in this study was only for plain-woven fabrics whereas there are other weave patterns which are of interest in the industry such as 8H satin, twill, etc. Hence, it would be interesting to update the analytical model in a manner that can be applied to other fiber architectures.  - To simplify the calculations in this study, the stress-strain behavior of the yarn was assumed linear (i.e., it has been discussed that why it did not affect our results notably). However, it is believed that the quality and accuracy of the predictions will be improved if the non-linear behavior of yarn is considered in the analyses.  117  - The effectiveness of the analytical model to predict wrinkling in FE simulations was only assessed in bias extension test. However, it is recommended that the same approach (flowchart presented in Figure 6-2) to be used for the analysis of a more realistic problem such as fabric forming process.  - It should be highlighted that all the discussions presented in Section 6.3.1 (i.e., parametric study) were based on the assumption that changing one of the properties of the yarns does not cause the change in the other properties (i.e., coupling effects between yarn properties were ignored). While in the real conditions there exist coupling effects between yarn properties. Accordingly, a comprehensive study in textile engineering field would be helpful for designing yarn properties (i.e., based on the results of the parametric study). This research would lead to a bridge between the micro/meso level properties of the fabric to its macro level properties such as wrinkling resistance behavior.  - For quite a long time, BHF was one of the most effective parameters to control wrinkling initiation. On one hand, BHF decreases the possibility of wrinkling formation and on the other hand, it raises the risk of yarn breakage. Hence, there is a trade-off between its positive and negative effects. After conducting the present study, it is evident that intra-yarn shear can be another potential parameter to delay shear wrinkling. As such, a new blank holder system can be designed to only create intra-yarn shear without causing large tension in yarns.  - In Chapter 3, it was demonstrated that the mesoscale deformation mechanisms of woven fabrics under bias extension test and picture frame test are distinct. In the fabric forming process, therefore, if blank holder system is used – boundary conditions applied to the ends of yarns are analogous to PFT, the results of picture frame test should be used in the 118  analyses. However, if blank holder system is not employed in the forming process, the results of bias extension test should be utilized.                                  119  References [1] S. Kalapakjian, S. Schmid, Manufacturing Engineering and Technology, Sixth Edit, Pearson Education, 2010. [2] S. El Wakil, Processes and Design for Manufacturing, Second Edi, PWS Publishing Company, 1998. [3] A. Hosseini, M. Kadkhodayan, A hybrid NN-FE approach to adjust blank holder gap over punch stroke in deep drawing process, Int. J. Adv. Manuf. Technol. 71 (2014) 337–355. doi:10.1007/s00170-013-5479-7. [4] M. Kadkhodayan, A. Hosseini, The Effect of Blank Holder Gap Profile upon Section Thickness Distribution in Deep Drawing Process, in: ICME, Tehran, 2011. [5] M. Kadkhodayan, A. Hosseini, The Optimization of Blank Holder Gap Profile in Deep Drawing Process, in: ICME, Tehran, 2011. [6] S.E. Arnold, M.P.F. Sutcliffe, W.L.A. Oram, Experimental measurement of wrinkle formation during draping of non-crimp fabric, Compos. Part A Appl. Sci. Manuf. 82 (2016) 159–169. doi:10.1016/j.compositesa.2015.12.011. [7] J.S. Lightfoot, M.R. Wisnom, K. Potter, Defects in woven preforms: Formation mechanisms and the effects of laminate design and layup protocol, Compos. Part A Appl. Sci. Manuf. 51 (2013) 99–107. doi:10.1016/j.compositesa.2013.04.004. [8] S. Bel, P. Boisse, F. Dumont, Analyses of the deformation mechanisms of non-crimp fabric composite reinforcements during preforming, Appl. Compos. Mater. 19 (2012) 513–528. doi:10.1007/s10443-011-9207-x. [9] B. Zhu, Sheet Forming of Woven Textile Composite Preforms: Formability and Wrinkling, 2007. [10] P. Boisse, N. Hamila, E. Vidal-Sallé, F. Dumont, Simulation of wrinkling during textile composite reinforcement forming. Influence of tensile, in-plane shear and bending stiffnesses, Compos. Sci. Technol. 71 (2011) 683–692. doi:10.1016/j.compscitech.2011.01.011. [11] J.B. Cattanach, G. Guff, F.N. Cogswell, The Processing of Thermoplastics Containing High Loadings of Long and Continuous Reinforcing Fibers, J. Polym. Eng. 6 (1986) 345–362. [12] J.A. Barnes, F.N. Cogswell, Transverse flow processes in continuous fibre-reinforced thermoplastic composites, Composites. 20 (1989) 38–42. doi:10.1016/0010-4361(89)90680-0. [13] R. Hill, A general theory of uniqueness and stability in elastic-plastic solids, J. Mech. Phys. Solids. 6 (1958) 236–249. doi:10.1016/0022-5096(58)90029-2. 120  [14] J. Cao, M.C. Boyce, Wrinkling behavior of rectangular plates under lateral constraint, Int. J. Solids Struct. 34 (1997) 153–176. doi:10.1016/S0020-7683(96)00008-X. [15] X. Wang, J. Cao, On the prediction of side-wall wrinkling in sheet metal forming processes, Int. J. Mech. Sci. 42 (2000) 2369–2394. doi:10.1016/S0020-7403(99)00078-8. [16] R.M. Jones, Buckling of Bars, Plates and Shells, Bull Ridge Publishing, Ann Arbor, 2006. [17] M.L. Gambhir, Stability Analysis and Design of Structures, Springer-Verlag Berlin Heidelberg, Berlin, 2004. doi:10.1061/(ASCE)ST.1943-541X.0001434. [18] Deep Drawing, (n.d.). http://thelibraryofmanufacturing.com/deep_drawing.html (accessed January 1, 2017). [19] B. Budiansky, J.W. Hutchinson, Buckling: Progress and Challenge, in: J.F. Besseling, A.M.A. van der Heijden (Eds.), Trends Solid Mech., Delft University Press, 1979. [20] Z. Marciniak, J.L. Duncan, S.J. Hu, Sheet deformation processes, 2002. doi:10.1007/BF02330250. [21] M.P. Henriques, T.J. Grilo, R.J.A. De Sousa, R.A.F. Valente, Numerical simulation and prediction of wrinkling defects in sheet metal forming, in: Stat. Comput. Tech. Manuf., 2013. doi:10.1007/978-3-642-25859-6_6. [22] W.F. Hosford, CaddellRobert, Metal forming: Mechanics and metallurgy, 2007. doi:10.1016/0378-3804(85)90124-X. [23] P. Boisse, N. Hamila, E. Guzman-Maldonado, A. Madeo, G. Hivet, F. Dell’Isola, The bias-extension test for the analysis of in-plane shear properties of textile composite reinforcements and prepregs: a review, Int. J. Mater. Form. (2016) 1–20. doi:10.1007/s12289-016-1294-7. [24] T. Gereke, O. Döbrich, M. Hübner, C. Cherif, Experimental and computational composite textile reinforcement forming: A review, Compos. Part A Appl. Sci. Manuf. 46 (2013) 1–10. doi:10.1016/j.compositesa.2012.10.004. [25] A.A. Skordos, C. Monroy Aceves, M.P.F. Sutcliffe, A simplified rate dependent model of forming and wrinkling of pre-impregnated woven composites, Compos. Part A Appl. Sci. Manuf. 38 (2007) 1318–1330. doi:10.1016/j.compositesa.2006.11.005. [26] G. Lebrun, M.N. Bureau, J. Denault, Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics, Compos. Struct. 61 (2003) 341–352. doi:10.1016/S0263-8223(03)00057-6. [27] M.R. Garnich, N.A. Klymyshyn, Multiscale analysis of stamp forming of a woven composite, J. Thermoplast. Compos. Mater. 26 (2013) 640–662. doi:10.1177/0892705711428654. 121  [28] S. V. Lomov, P. Boisse, E. Deluycker, F. Morestin, K. Vanclooster, D. Vandepitte, I. Verpoest, A. Willems, Full-field strain measurements in textile deformability studies, Compos. Part A Appl. Sci. Manuf. 39 (2008) 1232–1244. doi:10.1016/j.compositesa.2007.09.014. [29] J. Launay, G. Hivet, A. V. Duong, P. Boisse, Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements, Compos. Sci. Technol. 68 (2008) 506–515. doi:10.1016/j.compscitech.2007.06.021. [30] C. Mack, H.M. Taylor, The Fitting of Woven Cloth to Surfaces, J. Text. Inst. Trans. 47 (1956) T477–T488. doi:10.1080/19447027.1956.10750433. [31] J. Cao, R. Akkerman, P. Boisse, J. Chen, H.S. Cheng, E.F. de Graaf, J.L. Gorczyca, P. Harrison, G. Hivet, J. Launay, W. Lee, L. Liu, S. V. Lomov, A. Long, E. de Luycker, F. Morestin, J. Padvoiskis, X.Q. Peng, J. Sherwood, T. Stoilova, X.M. Tao, I. Verpoest, A. Willems, J. Wiggers, T.X. Yu, B. Zhu, Characterization of mechanical behavior of woven fabrics: Experimental methods and benchmark results, Compos. Part A Appl. Sci. Manuf. 39 (2008) 1037–1053. doi:10.1016/j.compositesa.2008.02.016. [32] U. Mohammed, C. Lekakou, M.G. Bader, Experimental studies and analysis of the draping of woven fabrics, Compos. Part A Appl. Sci. Manuf. 31 (2000) 1409–1420. doi:10.1016/S1359-835X(00)00080-4. [33] B.W. Senior, Flange wrinkling in deep-drawing operations, J. Mech. Phys. Solids. 4 (1956) 235–246. doi:10.1016/0022-5096(56)90032-1. [34] M.A. Shafaat, M. Abbasi, M. Ketabchi, Investigation into wall wrinkling in deep drawing process of conical cups, J. Mater. Process. Technol. 211 (2011) 1783–1795. doi:10.1016/j.jmatprotec.2011.05.026. [35] J. Cao, S.H. Cheng, H.P. Wang, C.T. Wang, Buckling of sheet metals in contact with tool surfaces, CIRP Ann. - Manuf. Technol. 56 (2007) 253–256. doi:10.1016/j.cirp.2007.05.059. [36] T.G. Clapp, H. Peng, Buckling of Woven Fabrics, Text. Res. J. 60 (1990) 641–645. doi:10.1177/004051759006001103. [37] R.D. Anandjiwala, Nonlinear Buckling of Woven Fabrics Part I: Elastic and Nonelastic Cases, Text. Res. J. 76 (2006) 160–168. doi:10.1177/0040517506057957. [38] T.G. Clapp, H. Peng, Part III : Experimental Validation of Theoretical Models, (n.d.) 641–645. [39] A.G. Prodromou, J. Chen, On the relationship between shear angle and wrinkling of textile composite preforms, Compos. Part A Appl. Sci. Manuf. 28A (1997) 491–503. [40] B. Zhu, T.X. Yu, J. Teng, X.M. Tao, Theoretical Modeling of Large Shear Deformation and Wrinkling of Plain Woven Composite, J. Compos. Mater. 43 (2008) 125–138. doi:10.1177/0021998308098237. 122  [41] O. Rozant, P.E. Bourban, J.A.E. Månson, Drapability of dry textile fabrics for stampable thermoplastic preforms, Compos. Part A Appl. Sci. Manuf. 31 (2000) 1167–1177. doi:10.1016/S1359-835X(00)00100-7. [42] L. Liu, J. Chen, X. Li, J. Sherwood, Two-dimensional macro-mechanics shear models of woven fabrics, Compos. Part A Appl. Sci. Manuf. 36 (2005) 105–114. doi:10.1016/j.compositesa.2004.07.004. [43] B. Zhu, T.X. Yu, H. Zhang, X.M. Tao, Experimental investigation of formability of commingled woven composite preform in stamping operation, Compos. Part B. 42 (2011) 289–295. [44] N. Liu, H. Yang, H. Li, S. Yan, Plastic wrinkling prediction in thin-walled part forming process: A review, Chinese J. Aeronaut. 29 (2016) 1–14. doi:10.1016/j.cja.2015.09.004. [45] H. Lin, M. Clifford, A.C. Long, M. Sherburn, Finite element modelling of fabric shear, Model. Simul. Mater. Sci. Eng. 17 (2008) 15008. doi:10.1088/0965-0393/17/1/015008. [46] J. Sirtautas, A.K. Pickett, P. L??picier, A mesoscopic model for coupled drape-infusion simulation of biaxial Non-Crimp Fabric, Compos. Part B Eng. 47 (2013) 48–57. doi:10.1016/j.compositesb.2012.09.088. [47] D. Ballhause, M. König, B. Kröplin, Modelling Fabric-Reinforced Membranes with the Discrete Element Method, in: E. Oñate, B. Kröplin (Eds.), Text. Compos. Inflatable Struct. II, Springer Netherlands, Dordrecht, 2008: pp. 51–67. doi:10.1007/978-1-4020-6856-0_4. [48] S. Gatouillat, E. Vidal-Sallé, P. Boisse, Advantages of the meso/macro approach for the simulation of fibre composite reinforcements, Int. J. Mater. Form. 3 (2010) 643–646. doi:10.1007/s12289-010-0852-7. [49] M. Komeili, Multi-Scale Characterization and Modeling of Shear-Tension Interaction in Woven Fabrics for Composite Forming and Structural Applications, 2014. [50] M. Komeili, A.S. Milani, On effect of shear-tension coupling in forming simulation of woven fabric reinforcements, Compos. Part B Eng. 99 (2016) 17–29. doi:10.1016/j.compositesb.2016.05.004. [51] P. Boisse, Meso-macro approach for composites forming simulation, J. Mater. Sci. 41 (2006) 6591–6598. doi:10.1007/s10853-006-0198-1. [52] N. Hamila, P. Boisse, F. Sabourin, M. Brunet, A semi-discrete shell finite element for textile composite reinforcement forming simulation, Int. J. Numer. Methods Eng. (2009) 1443–1466. doi:10.1002/nme. [53] N. Hamila, P. Boisse, Simulations of textile composite reinforcement draping using a new semi-discrete three node finite element, Compos. Part B Eng. 39 (2008) 999–1010. doi:10.1016/j.compositesb.2007.11.008. [54] P. Boisse, N. Hamila, A. Madeo, Modelling the development of defects during 123  composite reinforcements and prepreg forming., Philos. Trans. A. Math. Phys. Eng. Sci. 374 (2016) 20150269. doi:http://dx.doi.org/10.1098/rsta.2015.0269. [55] M. Haghi Kashani, A. Hosseini, F. Sassani, F.K. Ko, A. Milani, Understanding Different Types of Coupling in Mechanical Behavior of Woven Fabric Reinforcement: A Critical Review and Analysis, Compos. Struct. 179 (2017) 558–567. [56] P. Boisse, B. Zouari, A. Gasser, A mesoscopic approach for the simulation of woven fibre composite forming, Compos. Sci. Technol. 65 (2005) 429–436. doi:10.1016/j.compscitech.2004.09.024. [57] P. Badel, S. Gauthier, E. Vidal-Sallé, P. Boisse, Rate constitutive equations for computational analyses of textile composite reinforcement mechanical behaviour during forming, Compos. Part A Appl. Sci. Manuf. 40 (2009) 997–1007. doi:10.1016/j.compositesa.2008.04.015. [58] W.. Yu, F. Pourboghrat, K. Chung, M. Zampaloni, T.. Kang, Non-orthogonal constitutive equation for woven fabric reinforced thermoplastic composites, Compos. Part A Appl. Sci. Manuf. 33 (2002) 1095–1105. doi:10.1016/S1359-835X(02)00053-2. [59] Y. Aimène, E. Vidal-Sallé, B. Hagège, F. Sidoroff, P. Boisse, A Hyperelastic Approach for Composite Reinforcement Large Deformation Analysis, J. Compos. Mater. 44 (2010) 5–26. doi:10.1177/0021998309345348. [60] J. Cao, Prediction of plastic wrinkling using the energy method, J. Appl. Mech. - Trans. ASME. 66 (1999) 646–652. [61] Dassault Systèmes Simulia, Getting Started With ABAQUS, 2010. https://www.sharcnet.ca/Software/Abaqus610/Documentation/docs/v6.10/pdf_books/GET_STARTED_K.pdf. [62] S. Allaoui, P. Boisse, S. Chatel, N. Hamila, G. Hivet, D. Soulat, E. Vidal-Salle, Experimental and numerical analyses of textile reinforcement forming of a tetrahedral shape, Compos. Part A Appl. Sci. Manuf. 42 (2011) 612–622. doi:10.1016/j.compositesa.2011.02.001. [63] M.A. Khan, T. Mabrouki, E. Vidal-Sallé, P. Boisse, Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark, J. Mater. Process. Technol. 210 (2010) 378–388. doi:10.1016/j.jmatprotec.2009.09.027. [64] W. W, P. S, Wrinkled membranes III: numerical simulations, J Mech Mater Struct. 1 (2006) 63–95. [65] P. Harrison, M.J. Clifford, A.C. Long, Shear characterisation of viscous woven textile composites: A comparison between picture frame and bias extension experiments, Compos. Sci. Technol. 64 (2004) 1453–1465. doi:10.1016/j.compscitech.2003.10.015. [66] S. V. Lomov, Picture Frame Test of Woven Composite Reinforcements with a Full-Field Strain Registration, Text. Res. J. 76 (2006) 243–252. 124  doi:10.1177/0040517506061032. [67] B. Zhu, T.X. Yu, X.M. Tao, An experimental study of in-plane large shear deformation of woven fabric composite, Compos. Sci. Technol. 67 (2007) 252–261. doi:10.1016/j.compscitech.2006.08.011. [68] M. Haghi Kashani, A Coupled Non-Orthogonal Hypoelastic Constitutive Model for Simulation of Woven Fabrics, University of British Columbia, 2017. [69] A.S. Milani, J.A. Nemes, R.C. Abeyaratne, G.A. Holzapfel, A method for the approximation of non-uniform fiber misalignment in textile composites using picture frame test, Compos. Part A Appl. Sci. Manuf. 38 (2007) 1493–1501. doi:10.1016/j.compositesa.2007.01.008. [70] N. Kaynia, Instability-Induced Transformation of Interfacial Layers in Composites and its Multifunctional Applications, Massachusetts Institute of Technology, 2016. [71] B. Cornelissen, R. Akkerman, Analysis of yarn bending behaviour, in: ICCM-17 17th Int. Conf. Compos. Mater., 2009. http://doc.utwente.nl/74055/. [72] A.R. Zamani, S.O. Oyadiji, Analytical modelling of Kirschner wires in Ilizarov circular external fixator as pretensioned slender beams., J. R. Soc. Interface. 6 (2009) 243–56. doi:10.1098/rsif.2008.0251. [73] A.V. Perelmuter, V. Slivker, Handbook of Mechanical Stability in Engineering, World Science Publishing, Singapore, 2013. [74] E. Syerko, S. Comas-Cardona, C. Binetruy, Models of mechanical properties/behavior of dry fibrous materials at various scales in bending and tension: A review, Compos. Part A Appl. Sci. Manuf. 43 (2012) 1365–1388. doi:10.1016/j.compositesa.2012.03.012. [75] G. Hivet, A. V. Duong, A contribution to the analysis of the intrinsic shear behavior of fabrics, J. Compos. Mater. 45 (2010) 695–716. doi:10.1177/0021998310382315. [76] A. Hosseini, M. Haghi Kashani, F. Sassani, A. Milani, F.K. Ko, Identifying the Distinct Shear Wrinkling Behavior of Woven Composite Preforms under Bias Extension and Picture Frame Tests, Compos. Struct. (2017). doi:https://doi.org/10.1016/j.compstruct.2017.11.033. [77] F. Nosrat-Nezami, T. Gereke, C. Eberdt, C. Cherif, Characterisation of the shear-tension coupling of carbon-fibre fabric under controlled membrane tensions for precise simulative predictions of industrial preforming processes, Compos. Part A Appl. Sci. Manuf. 67 (2014) 131–139. doi:10.1016/j.compositesa.2014.08.030. [78] L. Li, Y. Zhao, H. gia nam Vuong, Y. Chen, J. Yang, Y. Duan, In-plane shear investigation of biaxial carbon non-crimp fabrics with experimental tests and finite element modeling, Mater. Des. 63 (2014) 757–765. doi:10.1016/j.matdes.2014.07.007. [79] M. H. Kashani, A. Hosseini, F. Sassani, F.K. Ko, A. Milani, A coupled non-orthogonal 125  hypoelastic model for woven fabrics, Submitted (2018). [80] B. El Abed, S. Msahli, H. Bel Hadj Salah, R. Zaouali, F. Sakli, Numerical simulation of woven fabric wrinkling, J. Text. Inst. 102 (2011) 77–86. doi:10.1080/00405000903509686. [81] W.R. Yu, P. Harrison, A. Long, Finite element forming simulation for non-crimp fabrics using a non-orthogonal constitutive equation, Compos. Part A Appl. Sci. Manuf. 36 (2005) 1079–1093. doi:10.1016/j.compositesa.2005.01.007. [82] A. Mallach, F. Ha rtel, F. Heieck, J.-P. Fuhr, P. Middendorf, M. Gude, Experimental comparison of a macroscopic draping simulation for dry non-crimp fabric preforming on a complex geometry by means of optical measurement, J. Compos. Mater. (2016) 1–13. doi:10.1177/0021998316670477. [83] M.S. Amabile, V. Eckers, T. Gries, Draping of non-crimp fabrics for fibre reinforced composites, Int. J. Mater. Form. 3 (2010) 647–650. doi:10.1007/s12289-010-0853-6. [84] H. Lin, J. Wang, A.C. Long, M.J. Clifford, P. Harrison, Predictive modelling for optimization of textile composite forming, Compos. Sci. Technol. 67 (2007) 3242–3252. doi:10.1016/j.compscitech.2007.03.040. [85] D.M. Mulvihill, M.P.F. Sutcliffe, Composites : Part A Effect of tool surface topography on friction with carbon fibre tows for composite fabric forming, Compos. Part A. 93 (2017) 199–206. doi:10.1016/j.compositesa.2016.10.017. [86] W.-R. Yu, P. Harrison, A.C. Long, Ideal Forming of Non-Crimp Fabric Preforms through Optimization of Blank Shape and Blank Holding Force, (2004). [87] J.S. Lee, S.J. Hong, W.R. Yu, T.J. Kang, The effect of blank holder force on the stamp forming behavior of non-crimp fabric with a chain stitch, Compos. Sci. Technol. 67 (2007) 357–366. doi:10.1016/j.compscitech.2006.09.009. [88] M. Kashani, M. Hejazi, A. Hosseini, F. Sassani, F. Ko, A. Milani, Toward enhanced forming simulation of woven fabrics using a coupled non-orthogonal hypoelastic constitutive model, integrated with a new wrinkling onset criterion, in: 32nd Am. Soc. Compos. Tech. Conf., West Lafayette - Purdue University, 2017. [89] A. Hosseini, M.H. Kashani, F. Sassani, A. Milani, F.K. Ko, A Mesoscopic Analytical Model to Predict the Onset of Wrinkling in Plain Woven Preforms under Bias Extension Shear Deformation, Materials (Basel). 10 (2017). doi: 10.3390/ma10101184.     126  8 Appendix: A Manual for ABAQUS User Subroutine UMAT_Woven Fabric 1. Introduction This manual contains review and evaluation of the ABAQUS User Subroutine UMAT_Woven Fabric. UMAT_Woven Fabric is a Fortran based code for the purpose of the numerical implementation of a new material model for woven fabrics in ABAQUS. This version of UMAT_Woven Fabric is based on the material model presented in [88]. Section 2 of this manual provides context for the subroutine function along with input and output parameters. Furthermore, an example is provided in Section 3 of this user manual to illustrate the effectiveness and validity of the newly developed UMAT compared to currently available material model for woven fabrics in ABAQUS. Eventually, the potential capabilities of the subroutine are discussed in Section 3. 2. Description of the Subroutine 2.1. Overview One of the advantages of woven fabrics over unidirectional prepregs is their superior formability thanks to the large shear deformation capability. There exists, however, a limit on the shear deformation of woven fabrics, namely the wrinkling. Applying tension to delay wrinkling during forming processes, a consequence of the inherent coupling in woven fabrics, is widely known to the industry. Yet, inherent coupling – change in the effective material properties of a given direction of the fabric due to the applied deformation in other directions – has not been fully understood and implemented in the forming simulations of fabric 127  reinforcements. Coupling should be incorporated in numerical optimization routines to accurately predict the deformation of the material under complex forming set-ups, and more importantly to predict a realistic yarn tension level that can suppress wrinkles. Toward this goal, a new coupled non-orthogonal model which predicts not only the stress-strain path, but also the critical point (shear wrinkling) of the woven fabrics should be proposed and implemented in a numerical simulation. The theory of the material model is thoroughly discussed in [55][88]. In summary, a coupled non-orthogonal hypoelastic constitutive model along with a criterion for wrinkling in terms of coupling are offered. To show its application, the model is implemented in ABAQUS via a UMAT code to predict the stress and strain fields as well as the onset of wrinkling under large shear deformation. The popular user subroutine to study the solid mechanics of woven fabrics is VFabric subroutine, which is deeply explained in ABAQUS (DS Simulia) user manual. In spite of taking non-orthogonality into account, the inherent coupling is ignored in the VFabric subroutine. The code provided here is able to take the inherent coupling into consideration, resulting in more accurate prediction of the response of woven fabrics under forming processes. Moreover, the VFabric subroutine is efficient for dynamic explicit analyses and could not be used for implicit solver. Although the explicit analysis is more practical for studying forming process, its precision is lower than implicit analysis. Hence, the presented subroutine has been written for the implicit solver in ABAQUS. When ABAQUS calls the UMAT_Woven Fabrics, the subroutine is provided with the inputs such as current stress as well as strain components and solution dependent state variables (SDV). The Subroutine undertakes the stress analysis and predicts the response of woven fabrics in the current time increment by providing outputs 128  at the end of the increment. These outputs (stress, strain and SDV) will be used as next increment inputs. 2.2. Definitions of Subroutine Functions, Inputs and Outputs The subroutine includes four basic modules. In the first module, the stress and strain arrays are transformed from the orthogonal global coordinate system to the local non-orthogonal coordinate system. The material stiffness matrix will be determined and assembled in the second module based on local strain array. The third module transfers the local stiffness matrix to the global stiffness matrix. Finally, the last module computes the global stress array and SDVs for next increment. 2.2.1. Subroutine Parameters The subroutine was written in accordance with ABAQUS User Material Manual Guide. The following lines were imported from ABAQUS to provide proper context for UMAT coding:     SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,      1 RPL,DDSDDT,DRPLDE,DRPLDT,        2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,    3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,      4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC) C       INCLUDE 'ABA_PARAM.INC' C       CHARACTER*80 CMNAME       DIMENSION STRESS(NTENS),STATEV(NSTATV),      1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),      2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),      3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),     4 JSTEP(4)  129  We also need to define the dimension of arrays used in the subroutine code (the STRANL and STRESSL are local strain and stress arrays, the rest of the arrays have been defined in the code, Section 3.3): DIMENSION DOBAR(4), DCBAR(4), DOSBAR(1), DCSBAR(1),      1STRANL(NTENS),STRESSL(NTENS),TRANSF(NTENS,NTENS),      2 TRANSFT(NTENS,NTENS),      3DBAR(NTENS,NTENS),DSTRANL(NTENS),DENTAL(NTENS,NTENS),     4 DDSDDET(NTENS,NTENS), DELTA(1) The coeficint and constant numbers that are used to define the stiffness function are coded as follows: PARAMETER(ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0, FIVE=5.D0, 1 AII=0.4518D12, BII=3.2748D10, CII=7.3238D8, DII=6.6648.D0, ………………. The basic arrays in this subroutine are: TRANSF - the transformation matrix form local coordinate to global  TRANFT - the transformation matrix from global coordinates to local  DDSDDET - the local material stiffness matrix STATEV(1) - Є1 STATEV(2) - Є2 STATEV(3) - Є12 STATEV(4) - 𝜎1 STATEV(5) - 𝜎2 STATEV(6) - 𝜎12 STATEV (NSTATV) - the wrinkling indicator (1=wrinkled, 0=Normal)   130  2.2.2. Transformation from global orthogonal coordinate system to local non-orthogonal coordinate system The details of the transformation between the orthogonal and non-orthogonal coordinate systems are presented in [55]. The transformation matrix can be written as:  This matrix implemented in the code as: TRANSF(1,1)=COS(STRAN(3)/TWO)**TWO TRANSF(1,2)=SIN(STRAN(3)/TWO)**TWO TRANSF(1,3)=SIN(STRAN(3)) TRANSF(2,1)=SIN(STRAN(3)/TWO)**TWO TRANSF(2,2)=COS(STRAN(3)/TWO)**TWO TRANSF(2,3)=SIN(STRAN(3)) TRANSF(3,1)=SIN(STRAN(3)/TWO)*COS(STRAN(3)/TWO) TRANSF(3,2)=SIN(STRAN(3)/TWO)*COS(STRAN(3)/TWO) TRANSF(3,3)=ONE    2.2.3. Stiffness functions in the non-orthogonal coordinate system The coupled stiffness functions in the local non-orthogonal coordinate system are determined based on the local strains. For a new woven fabric, a user needs only to modify this section based on the characterization results. As a matter of fact, this section of the code is the 131  main body of the code which presents the coupled constitutive model in the non-orthogonal coordinate system. DOBAR(1)=-AII*STATEV(1)**FIVE+BII*STATEV(1)**FOUR 1 -CII*STATEV(1)**THREE+DII*STATEV(1)**TWO-EII*STATEV(1)+FII       IF(STATEV(1).GT.EPSILONII) THEN         DOBAR(1)=-AII*EPSILONII**FIVE+BII*EPSILONII**FOUR 1 -CII*EPSILONII**THREE+DII*EPSILONII**TWO-EII*EPSILONII+FII       END IF       IF(STATEV(1).LT.ZERO) THEN         DOBAR(1)=FII       END IF …………. . . . . . . . . . C            IF(STATEV(1).LT.ZERO.OR.STATEV(1).EQ.ZERO) THEN         IF(STATEV(2).LT.ZERO.OR.STATEV(2).EQ.ZERO) THEN           DOSBAR(1)=ZPTZS           IF(STATEV(3).GT.ZPZS)Then             DOSBAR(1)=DOSBAR(1)-ZPOZE           END IF           DCSBAR(1)=ONE         END IF 132        END IF C 2.2.4. Generating the local stiffness matrix The local stifness matrix is assmebled as follows : DBAR(1,1)=DOBAR(1)*DCBAR(1)       DBAR(1,2)=DOBAR(3)*DCBAR(3)        DBAR(1,3)=ZERO        DBAR(2,1)=DOBAR(4)*DCBAR(4)        DBAR(2,2)=DOBAR(2)*DCBAR(2)        DBAR(2,3)=ZERO        DBAR(3,1)=ZERO        DBAR(3,2)=ZERO        DBAR(3,3)=DCSBAR(1)*DOSBAR(1)  2.2.5. Stress analysis in Global Coordinate System After trasnfering stiffness matrix form local to global coordinate system (by using TRANSF array), the stress analysis and updating the state varibales are as follows:   DO K1=1, NTENS         DO K2=1, NTENS           DDSDDE(K1,K2)=DDSDDET(K1,K2)         END DO       END DO          DO K1=1, NTENS         DO K2=1, NTENS           STRESS(K1)=STRESS(K1)+DDSDDE(K1,K2)*DSTRAN(K2)         END DO       END DO Updating SDVs: DO K1=1, NTENS         STATEV(K1)=STRANL(K1) 133          STATEV(NTENS+K1)=STRESSL(K1)       END DO    2.2.6. Wrinkling prediction The criterion for prdiction of wrinkling occurrence is developed in [76][89]. The last step  of the material subroutine is to predict wrinkling initiation which is coded as: DELTA(1)=ABS(STRESSL(1))**(ONE+HALF)/(ABS(STRANL(1))**(ONE+HALF)+      1 STRANL(1)**HALF-TANH(ABS(STRANL(1))**HALF))       IF(STATEV(3).LT.-ACOS(ANGS).OR.STATEV(3).GT.DELTA(1)) THEN         STATEV(3*NTENS)=ONE       END IF       STATEV(NSTATV)=DELTA(1)  3. Use of the Subroutine in ABAQUS To use the UMAT_Woven Fabrics, it is required to select the user material in the defining properties section. Afterward, when it comes to define the job, we should provide the path for the Fortran file of UMAT_Woven Fabrics in the special section of the Job bar. 3.1. An example to demonstrate the advantage of the new subroutine To prove the effectiveness and validity of the written subroutine, an example is provided below. Firstly, 2% tension in the transverse direction is applied and then tension in the main direction is applied up to 4%. Figure A-1 demonstrates the higher accuracy of the coupled model (the new model) compared to the uncoupled model (the currently available woven fabrics subroutine in ABAQUS - VFabric). 134   Figure A-1. Comparison between the newly dveloped material model for woven fabrics and the currently available woven fabrics model – Vfabric - in ABAQUS.   3.2. Potential Application The results showed the advantages of the presented UMAT for woven fabrics over the VFabric subroutine which is available in ABAQUS. In fact, the new subroutine predicts the behavior of woven fabrics under forming processes with higher precision in comparison with VFabric subroutine. Moreover, thanks to taking the inherent coupling into account, the manufacturing process parameters such as the applied pressure on the blank holder can be optimized using this new UMAT to prevent wrinkling. Furthermore, since the presented code was written for implicit solutions, coupled thermal-mechanical analyses to capture the effect of temperature on the mechanical properties is also possible. The last, but not least, it can be used for various types of woven fabrics by just modifying the stiffness functions section based on the characterization results of the given woven fabrics. 3.3. UMAT_WOVEN FABRICS C **************************************************************** C **************************************************************** C                        UMAT FOR Woven Fabircs C                        UMAT_WOVEN FABRICS 135  C             User Material Subroutine for Coupled Non-Orthogonal  C             Constitutive Modeling of Woven Fabrics  C    C    January, 2017 C    Version  1.1  C **************************************************************** C    Applications: C                2-D Planer Stress/Strain C                Conventional Shell Element S4 C **************************************************************** C    Authors: C                Masoud Hejazi C        Masoud Haghi Kashani C                Abbas Hosseini C                Farrokh Sassani C                Frank ko C                Abbas Milani C C    Department of Mechanical Engineering, C    University Of British Columbia, Canada  C            C  C ---------------------------------------------------------------- C       SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,      1 RPL,DDSDDT,DRPLDE,DRPLDT,      2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,      3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,      4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC) C       INCLUDE 'ABA_PARAM.INC' C       CHARACTER*80 CMNAME       DIMENSION STRESS(NTENS),STATEV(NSTATV),      1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),      2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),      3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),      4 JSTEP(4) C C LOCAL ARRAYS C ---------------------------------------------------------------- C DOBAR -  PURE TENSILE STIFNESS MODULES  C DCBAR -  COUPLING TENSILE STIFNESS INDUCED MODULUS C DOSBAR-  PURE SHEAR STIFNESS MODULES  C DCSBAR-  COUPLING SHEAR STIFNESS INDUCED MODULUS C TRANSF-  TRANSFORMATION MATRIX 136  C DBAR-    JACOBIAN NON-ORTHOGONAL MATRIX C STRANL-  NON-ORTHOGONAL STRAIN ARRAY C STREESL- NON-ORTHOGONAL STRESS ARRAY C       DIMENSION DOBAR(4), DCBAR(4), DOSBAR(1), DCSBAR(1),      1 STRANL(NTENS), STRESSL(NTENS), TRANSF(NTENS,NTENS),      2 TRANSFT(NTENS,NTENS),      3 DBAR(NTENS,NTENS), DSTRANL(NTENS), DENTAL(NTENS,NTENS),      4 DDSDDET(NTENS,NTENS), DELTA(1) C C       PARAMETER(ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0, FIVE=5.D0,      1 AII=0.4518D12, BII=3.2748D10, CII=7.3238D8, DII=6.6648D6,      2 EII=1.3571D4, FII=0.3344D2, GII=27.47D0,      3 EPSILONII=0.03D0, OPTH=1.3D0, ZPSF=0.83D0, ZPTE=0.28D0,      4 ZPTF=0.025D0, FIFFOUR=54.D0, TPS=2.6D0, THUND=1200.D0,      5 ZPZT=0.02D0, THOO=211.D0, ZPTZS=0.1587D0, ZPOZE=0.1071D0,      6 ZPZS=0.07D0, SVHU=1700.D0, TAH=2.5D0, OFFN=14597.D0,      7 OPTT=1.22D0, TWTH=12.D3, HALF=0.5D0, OFSZ=1570.D0,      8 OFSE=15685.D0, OOTF=11353.D0, ZPZOS=0.017D0, ZERO=0.D0,       9 TOL=1.D-10, ANGS=0.921D0) C C TRANSFORMATION MATRIX, ORTHOGONAL TO NONORTHOGONAL C-----------------------------------------------------------------  C       TRANSF(1,1)=COS(STRAN(3)/TWO)**TWO       TRANSF(1,2)=SIN(STRAN(3)/TWO)**TWO       TRANSF(1,3)=SIN(STRAN(3))       TRANSF(2,1)=SIN(STRAN(3)/TWO)**TWO       TRANSF(2,2)=COS(STRAN(3)/TWO)**TWO       TRANSF(2,3)=SIN(STRAN(3))       TRANSF(3,1)=SIN(STRAN(3)/TWO)*COS(STRAN(3)/TWO)       TRANSF(3,2)=SIN(STRAN(3)/TWO)*COS(STRAN(3)/TWO)       TRANSF(3,3)=ONE                        DO K1=1, NTENS         DO K2=1, NTENS         TRANSFT(K1,K2)=TRANSF(K2,K1)         END DO       END DO C C TRANSFORMING TO NONORTHOGONAL STRANL AND DSTRANL C --------------------------------------------------------------- C 137        DO K1=1, NTENS         STRANL(K1)=ZERO         DSTRANL(K1)=ZERO       END DO              DO K1=1, NTENS         DO K2=1, NTENS           STRANL(K1)=STRANL(K1)+TRANSFT(K1,K2)*(STRAN(K2)+DSTRAN(K2))         END DO       END DO              DO K1=1, NTENS         DO K2=1, NTENS           DSTRANL(K1)=DSTRANL(K1)+TRANSFT(K1,K2)*DSTRAN(K2)         END DO       END DO C C C---------------------------------------------------------------- C --------------------------------------------------------------- C --------------------------------------------------------------- C FUNCTIONS DOBAR DCBAR DOSBAR DC C      -------------------------------------------------- C DOBAR.. C -----------------------------------------------------------------        DOBAR(1)=-AII*STATEV(1)**FIVE+BII*STATEV(1)**FOUR      1 -CII*STATEV(1)**THREE+DII*STATEV(1)**TWO-EII*STATEV(1)+FII       IF(STATEV(1).GT.EPSILONII) THEN         DOBAR(1)=-AII*EPSILONII**FIVE+BII*EPSILONII**FOUR      1 -CII*EPSILONII**THREE+DII*EPSILONII**TWO-EII*EPSILONII+FII       END IF       IF(STATEV(1).LT.ZERO) THEN         DOBAR(1)=FII       END IF       C            DOBAR(2)=-AII*STATEV(2)**FIVE+BII*STATEV(2)**FOUR      1 -CII*STATEV(2)**THREE+DII*STATEV(2)**TWO-EII*STATEV(2)+FII       IF(STATEV(2).GT.EPSILONII) THEN         DOBAR(2)=-AII*EPSILONII**FIVE+BII*EPSILONII**FOUR      1 -CII*EPSILONII**THREE+DII*EPSILONII**TWO-EII*EPSILONII+FII       END IF       IF(STATEV(2).LT.ZERO) THEN         DOBAR(2)=FII       END IF 138   C            DOBAR(3)=(TPS+THUND*STATEV(2)-THUND*(STATEV(2)-ZPZT))/OPTH       IF(STATEV(2).LT.ZPZT) THEN         DOBAR(3)=(TPS+THUND*STATEV(2))/OPTH       END IF       IF(STATEV(2).LT.ZERO) THEN         DOBAR(3)=TPS       END IF C             DOBAR(4)=(TPS+THUND*STATEV(1)-THUND*(STATEV(1)-ZPZT))/OPTH       IF(STATEV(1).LT.ZPZT) THEN         DOBAR(4)=(TPS+THUND*STATEV(1))/OPTH       END IF       IF(STATEV(1).LT.ZERO) THEN         DOBAR(4)=TPS       END IF C ----------------------------------------------------------------- C C DCBAR C -----------------------------------------------------------------       DCBAR(1)=ONE-GII*(STATEV(2)+TOL)**OPTH       IF(STATEV(1).GT.ZPTF) THEN         DCBAR(1)=DCBAR(1)-ZPSF*(STATEV(2)+TOL)**ZPTE       END IF       IF(STATEV(2).GT.ZPTF) THEN         IF(STATEV(1).LT.EPSILONII.OR.STATEV(1).EQ.EPSILONII) THEN           DCBAR(1)=DCBAR(1)+(TWO-FIFFOUR*STATEV(1))         END IF         IF(STATEV(1).GT.EPSILONII) THEN         DCBAR(1)=DCBAR(1)+(TWO-FIFFOUR*EPSILONII)         END IF       END IF                     DCBAR(2)=ONE-GII*(STATEV(1)+TOL)**OPTH       IF(STATEV(2).GT.ZPTF) THEN         DCBAR(2)=DCBAR(2)-ZPSF*(STATEV(1)+TOL)**ZPTE       END IF       IF(STATEV(1).GT.ZPTF) THEN         IF(STATEV(2).LT.EPSILONII.OR.STATEV(2).EQ.EPSILONII) THEN           DCBAR(2)=DCBAR(2)+(TWO-FIFFOUR*STATEV(2))         END IF         IF(STATEV(2).GT.EPSILONII) THEN         DCBAR(2)=DCBAR(2)+(TWO-FIFFOUR*EPSILONII)         END IF 139        END IF            DCBAR(3)=ONE+THOO*STATEV(1) C            DCBAR(4)=ONE+THOO*STATEV(2)              IF(STATEV(2).LT.ZERO.OR.STATEV(1).LT.ZERO) THEN        DO K1=1, 4          DCBAR(K1)=ONE        END DO       END IF C ----------------------------------------------------------------- C ----------------------------------------------------------------- C C DOSBAR AND DCSBAR C -----------------------------------------------------------------       DOSBAR(1)=ZPTZS C             DCSBAR(1)=ONE+SVHU*(STATEV(1)+STATEV(2))       If(STATEV(1).GT.ZPTF) THEN         DCSBAR(1)=DCSBAR(1)+OFFN*(STATEV(1)-ZPTF)**OPTT       END IF C             IF(STATEV(2).GT.ZPTF) THEN         DCSBAR(1)=DCSBAR(1)+OFFN*(STATEV(2)-ZPTF)**OPTT C                 IF(STATEV(1).GT.ZPTF) THEN           DCSBAR(1)=DCSBAR(1)      1    +TWTH*((STATEV(1)-ZPTF)*(STATEV(2)-ZPTF))**HALF         END IF C                 END IF C       C             IF(STATEV(3).GT.ZPZOS) THEN         DCSBAR(1)=DCSBAR(1)-OFSZ*(STATEV(1)+STATEV(2)) C                 If(STATEV(1).GT.ZPTF) THEN           DCSBAR(1)=DCSBAR(1)-OFSE*(STATEV(1)-ZPTF)**OPTH         END IF C               IF(STATEV(2).GT.ZPTF) THEN           DCSBAR(1)=DCSBAR(1)-OFSE*(STATEV(2)-ZPTF)**OPTH C                   IF(STATEV(1).GT.ZPTF) THEN             DCSBAR(1)=DCSBAR(1) 140       1      -OOTF*((STATEV(1)-ZPTF)*(STATEV(2)-ZPTF))**HALF C                     END IF         END IF       END IF C            IF(STATEV(1).LT.ZERO.OR.STATEV(1).EQ.ZERO) THEN         IF(STATEV(2).LT.ZERO.OR.STATEV(2).EQ.ZERO) THEN           DOSBAR(1)=ZPTZS           IF(STATEV(3).GT.ZPZS)Then             DOSBAR(1)=DOSBAR(1)-ZPOZE           END IF           DCSBAR(1)=ONE         END IF       END IF C          C -------------------------------------------------------------------------- C -------------------------------------------------------------------------- C C C ASSEMBELING NONORTHOGONAL JACOBIAN MATRIX  C --------------------------------------------------------------------------       DBAR(1,1)=DOBAR(1)*DCBAR(1)       DBAR(1,2)=DOBAR(3)*DCBAR(3)       DBAR(1,3)=ZERO       DBAR(2,1)=DOBAR(4)*DCBAR(4)       DBAR(2,2)=DOBAR(2)*DCBAR(2)       DBAR(2,3)=ZERO       DBAR(3,1)=ZERO       DBAR(3,2)=ZERO       DBAR(3,3)=DCSBAR(1)*DOSBAR(1)  C--------------------------------------------------------------------------- C -------------------------------------------------------------------------- C --------------------------------------------------------------------------  C CALCULATING ORTHOGONAL JACOBIAN MATRIX AND STRESS  C ------------------------------------------------------------------------       C       DO K1=1, NTENS C        STRESS(K1)=ZERO         DO K2=1, NTENS          DDSDDET(K1,K2)=ZERO          DENTAL(K1,K2)=ZERO         END DO 141        END DO       DO K1=1, NTENS         DO K2=1, NTENS           DO K3=1, NTENS            DENTAL(K1,K2)=DENTAL(K1,K2)+TRANSF(K1,K3)*DBAR(K3,K2)           END DO         END DO       END DO                 DO K1=1, NTENS         DO K2=1, NTENS           DO K3=1, NTENS             DDSDDET(K1,K2)=DDSDDET(K1,K2)+DENTAL(K1,K3)*TRANSFT(K3,K2)           END DO         END DO       END DO               DO K1=1, NTENS         DO K2=1, NTENS           DDSDDE(K1,K2)=DDSDDET(K1,K2)         END DO       END DO              DO K1=1, NTENS         DO K2=1, NTENS           STRESS(K1)=STRESS(K1)+DDSDDE(K1,K2)*DSTRAN(K2)         END DO       END DO C CALCULATING NONORTHOGONAL STRESSES C---------------------------------------------------------------------------       DO K1=1, NTENS         STRESSL(K1)=ZERO       END DO       DO K1=1, NTENS         DO K2=1, NTENS           STRESSL(K1)=STRESSL(K1)+TRANSFT(K1,K2)*STRESS(K2)         END DO       END DO C-------------------------------------------------------------------------- C-------------------------------------------------------------------------- C        DO K1=1, NTENS         STATEV(K1)=STRANL(K1)         STATEV(NTENS+K1)=STRESSL(K1) 142        END DO   C WRINKLE PREDECTION C--------------------------------------------------------------------------       DELTA(1)=ABS(STRESSL(1))**(ONE+HALF)/(ABS(STRANL(1))**(ONE+HALF)+      1 STRANL(1)**HALF-TANH(ABS(STRANL(1))**HALF))       IF(STATEV(3).LT.-ACOS(ANGS).OR.STATEV(3).GT.DELTA(1)) THEN         STATEV(3*NTENS)=ONE       END IF       STATEV(NSTATV)=DELTA(1)    C ------------------------------------------------------------------------- C--------------------------------------------------------------------------       C             RETURN       END  3.4. The Input File of the Example  *Heading ** Job name: picturerame Model name: Model-1 ** Generated by: Abaqus/CAE 6.14-2 *Preprint, echo=NO, model=NO, history=NO, contact=NO ** ** PARTS ** *Part, name=Part-1 *Node       1,           0.,           0.       2,           1.,           0.       3,           2.,           0.       4,           3.,           0.       5,           4.,           0.       6,           5.,           0.       7,           0.,           1.       8,           1.,           1.       9,           2.,           1.      10,           3.,           1.      11,           4.,           1.      12,           5.,           1.      13,           0.,           2.      14,           1.,           2.      15,           2.,           2.      16,           3.,           2.      17,           4.,           2.      18,           5.,           2. 143       19,           0.,           3.      20,           1.,           3.      21,           2.,           3.      22,           3.,           3.      23,           4.,           3.      24,           5.,           3.      25,           0.,           4.      26,           1.,           4.      27,           2.,           4.      28,           3.,           4.      29,           4.,           4.      30,           5.,           4.      31,           0.,           5.      32,           1.,           5.      33,           2.,           5.      34,           3.,           5.      35,           4.,           5.      36,           5.,           5. *Element, type=CPS4  1,  1,  2,  8,  7  2,  2,  3,  9,  8  3,  3,  4, 10,  9  4,  4,  5, 11, 10  5,  5,  6, 12, 11  6,  7,  8, 14, 13  7,  8,  9, 15, 14  8,  9, 10, 16, 15  9, 10, 11, 17, 16 10, 11, 12, 18, 17 11, 13, 14, 20, 19 12, 14, 15, 21, 20 13, 15, 16, 22, 21 14, 16, 17, 23, 22 15, 17, 18, 24, 23 16, 19, 20, 26, 25 17, 20, 21, 27, 26 18, 21, 22, 28, 27 19, 22, 23, 29, 28 20, 23, 24, 30, 29 21, 25, 26, 32, 31 22, 26, 27, 33, 32 23, 27, 28, 34, 33 24, 28, 29, 35, 34 25, 29, 30, 36, 35 *Nset, nset=Set-1, generate   1,  36,   1 144  *Elset, elset=Set-1, generate   1,  25,   1 ** Section: Section-1 *Solid Section, elset=Set-1, material=Material-1 1.3, *End Part **   ** ** ASSEMBLY ** *Assembly, name=Assembly **   *Instance, name=Part-1-1, part=Part-1 *End Instance **   *Nset, nset=Set-1, instance=Part-1-1, generate   1,  31,   6 *Elset, elset=Set-1, instance=Part-1-1, generate   1,  21,   5 *Nset, nset=Set-2, instance=Part-1-1, generate  1,  6,  1 *Elset, elset=Set-2, instance=Part-1-1, generate  1,  5,  1 *Nset, nset=Set-3, instance=Part-1-1, generate   6,  36,   6 *Elset, elset=Set-3, instance=Part-1-1, generate   5,  25,   5 *Nset, nset=Set-4, instance=Part-1-1, generate  31,  36,   1 *Elset, elset=Set-4, instance=Part-1-1, generate  21,  25,   1 *Nset, nset=Set-5, instance=Part-1-1  1, *End Assembly **  ** MATERIALS **  *Material, name=Material-1 *Depvar      30, *User Material, constants=1 1., ** ---------------------------------------------------------------- **  ** STEP: Step-1 **  145  *Step, name=Step-1, nlgeom=NO *Static 0.0625, 1., 1e-05, 0.0625 **  ** BOUNDARY CONDITIONS **  ** Name: BC-3 Type: Displacement/Rotation *Boundary Set-3, 2, 2, 1. ** Name: BC-4 Type: Displacement/Rotation *Boundary Set-4, 1, 1, 1. ** Name: BC-5 Type: Displacement/Rotation *Boundary Set-5, 1, 1 Set-5, 2, 2 **  ** OUTPUT REQUESTS **  *Restart, write, frequency=0 **  ** FIELD OUTPUT: F-Output-1 **  *Output, field, variable=PRESELECT **  ** HISTORY OUTPUT: H-Output-3 **  *Output, history *Element Output, elset=Part-1-1.Set-1 E11, E12, E22 **  ** HISTORY OUTPUT: H-Output-2 **  *Element Output, elset=Part-1-1.Set-1 S11, S12, S22, SDV **  ** HISTORY OUTPUT: H-Output-1 **  *Output, history, variable=PRESELECT *End Step 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0363919/manifest

Comment

Related Items