UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A three-phase integrated flow-stress framework for process modelling of composite materials Amini Niaki, Sina 2017

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

Item Metadata

Download

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

Full Text

A THREE-PHASE INTEGRATED FLOW-STRESS FRAMEWORK FOR PROCESS MODELLING OF COMPOSITE MATERIALS by  Sina Amini Niaki  B.Sc., Mechanical Engineering, Sharif University of Technology, 2009 M.A.Sc., Mechanical Engineering, Sharif University of Technology, 2011  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES  (Civil Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  September 2017  © Sina Amini Niaki, 2017 ii  Abstract  An accurate and efficient predictive tool for process modelling of fibre-reinforced polymeric composite materials is highly desirable because of high production cost and substantial risk involved in their manufacturing. Process modelling of composite materials is complex due to multitude of physical and chemical processes such as flow of resin and gaseous volatiles through the fibre-bed, thermochemical changes, and stress development. These phenomena are usually simulated sequentially in a so-called “integrated sub-model” framework.  In the sequential method, (i) mapping of the results from one state to another is performed in a tedious and inefficient process, (ii) the concurrent occurrence of resin flow and stress development is ignored, (iii) the history of pressure during the flow regime of processing is relinquished, and (iv) the local spatial and temporal variations in resin gelation is not captured. Incorporating the main steps of process modeling into a unified module that captures the various phenomena in a continuous manner would help overcome the foregoing drawbacks of the sequential approach.  In this thesis, a new methodology is presented to integrate the simulation of resin and gas flow and stress development into a unified computational modelling framework. The governing equations are developed for the general case of a composite system that initially is regarded as a three-phase liquid-gas-solid system and, as a consequence of curing, the resin undergoes a transition from a fluid-like state into an elastic solid forming a solid cured composite material. The employed constitutive equations provide a continuous representation of the evolving iii  material behaviour while maintaining consistency with the formulations that are typically used to represent the material at each of the two processing extremes. The model is implemented in a 2D plane strain displacement-velocity-pressure (u-v-P) finite element code developed in MATLAB®. Various numerical examples are presented to demonstrate the capability of the Integrated Flow-Stress (IFS) model to predict the flow-compaction and stress development throughout the curing process of thermoset composite materials. The interactive effects of resin flow, gas flow, and stress development are investigated and comparisons are made with the predicted results obtained from the application of the stress model alone.     iv  Lay Summary  A composite material is made of several constituents with significantly different properties, when combined result in a material with all desired properties for applications in different key economic sectors such as renewable energy, aerospace, and transportation. Modelling tools are used to simulate their manufacturing virtually before the real risky and costly process. Because of their complex behavior during manufacturing, the available simulation tools are unable to accurately predict the outcome of processing when applied to the entire process. In this thesis, an accurate mechanistic model is developed such that a single numerical tool can be used to model the whole manufacturing process. All physical states of constituent materials as well as the evolution of their properties are considered. Consequently, the need to have separate models for different stages of processing is avoided and a fast, convenient, and accurate simulation tool is provided for engineers and researchers in composite manufacturing industry.   v  Preface  This thesis presents the research performed by Sina Amini Niaki. The research was supervised by Dr. Reza Vaziri and co-supervised by Dr. Anoush Poursartip at The University of British Columbia.  The two-phase isotropic model in Chapter 3 (section 3.1), and the related examples in Chapter 5 (section 5.1 and 5.2) has been published. S. A. Niaki, A. Forghani, R. Vaziri, and A. Poursartip (2017). A two-phase integrated flow-stress process model for composites with application to highly compressible phases. Mechanics of Materials, 109, 51-66. Most parts of this paper were written by Sina Amini Niaki.  The three-phase isotropic model in Chapter 3 (section 3.2), and the related examples in Chapter 5 (section 5.3) have been submitted for publication as a journal paper. S. A. Niaki, A. Forghani, R. Vaziri, and A. Poursartip (2017). A Three-Phase Integrated Flow-Stress Model for Processing of Composites. Most parts of this paper were written by Sina Amini Niaki.  The two-phase orthotropic model in Chapter 4 and the related examples in Chapter 5 (section 5.2) is being submitted for publication as a journal paper. S. A. Niaki, A. Forghani, R. Vaziri, and A. Poursartip (2017). A Two-Phase Orthotropic Integrated Flow-Stress Model for Process Simulation of Composite Materials.  Most parts of this paper were written by Sina Amini Niaki.  vi  Portions of Chapters 3 and 5 have been presented at the 10th Canadian-International Conference on Composites (CANCOM2017) by Sina Amini Niaki. S. A. Niaki, A. Forghani, R. Vaziri, and A. Poursartip “A Two-Phase Integrated Flow-Stress Process Model for Polymeric Composite Materials” 10th Canadian-International Conference on Composites (CANCOM2017), Ottawa, Canada, July 2017.  Portions of Chapters 4 and 5 have been submitted to the 21st International Conference on Composite Materials (ICCM21) and will be presented by Sina Amini Niaki. S. A. Niaki, A. Forghani, R. Vaziri, and A. Poursartip “A Two-Phase Orthotropic Model for Integrating Flow and Stress Development in Processing of Composite Materials” 21st International Conference on Composite Materials (ICCM21), Xi’an, China, August 2017.  vii  Table of Contents  Abstract .......................................................................................................................................... ii Lay Summary ............................................................................................................................... iv Preface .............................................................................................................................................v Table of Contents ........................................................................................................................ vii List of Tables ............................................................................................................................... xii List of Figures ............................................................................................................................. xiii Acknowledgements ......................................................................................................................xx Dedication ................................................................................................................................... xxi Chapter 1: Introduction ................................................................................................................1 1.1 Motivations ..................................................................................................................... 2 1.2 Objectives and Outline of Thesis .................................................................................... 3 1.3 Figures............................................................................................................................. 6 Chapter 2: Literature Review .......................................................................................................7 2.1 Flow Models ................................................................................................................... 8 2.1.1 Composite Materials Compaction............................................................................... 9 2.1.2 Liquid Composite Molding (LCM) .......................................................................... 12 2.1.3 Three-Phase Flow Models ........................................................................................ 15 2.2 Stress Models ................................................................................................................ 19 2.3 Integrated Approaches .................................................................................................. 21 2.4 Porosity ......................................................................................................................... 22 2.5 Summary ....................................................................................................................... 24 viii  2.6 Figures........................................................................................................................... 26 Chapter 3: Isotropic Integrated Flow-Stress Model (IFS-Isotropic) ......................................29 3.1 Formulation of Two-Phase Model (2IFS-Isotropic) ..................................................... 30 3.1.1 Conservation Equations ............................................................................................ 30 3.1.1.1 Mass Conservation ............................................................................................ 31 3.1.1.2 Equilibrium of the Matrix ................................................................................. 32 3.1.1.3 Equilibrium of the System ................................................................................ 33 3.1.2 Constitutive Equations .............................................................................................. 33 3.1.2.1 Bulk Response .................................................................................................. 34  Unsolidified Matrix Phase ........................................................................... 34  Solidifying Matrix Phase ............................................................................. 36 3.1.2.2 Effective Stress (Non-Bulk) Response ............................................................. 39  Unsolidified Matrix Phase ........................................................................... 40  Solidifying Matrix Phase ............................................................................. 41 3.1.3 Summary of the Governing Equations (2IFS-Isotropic) ........................................... 43 3.1.3.1 Special Case: Unsolidified Material (𝝀 = 𝟎) .................................................... 44 3.1.3.2 Special Case: Fully Solidified Material (𝝀 = 𝟏) ............................................... 46 3.1.4 Finite Element Implementation (2IFS-Isotropic)...................................................... 47 3.1.4.1 Volume Fraction Calculation ............................................................................ 49 3.2 Formulation of the Three-Phase Model (3IFS-Isotropic) ............................................. 53 3.2.1 Conservation Equations ............................................................................................ 53 3.2.1.1 Mass Conservation ............................................................................................ 53 3.2.1.2 Equilibrium of the Fluid Phases ........................................................................ 54 ix  3.2.1.3 Pressure Equilibrium ......................................................................................... 55 3.2.1.4 Equilibrium of the Three-Phase System ........................................................... 55 3.2.2 Constitutive Equations .............................................................................................. 56 3.2.2.1 Rate of Compressibility .................................................................................... 56 3.2.2.2 Effective Stress ................................................................................................. 58 3.2.3 Summary of the Governing Equations (3IFS-Isotropic) ........................................... 59 3.2.4 Finite Element Implementation (3IFS-Isotropic)...................................................... 60 3.2.4.1 Volume Fraction Calculation ............................................................................ 62 3.3 Summary ....................................................................................................................... 62 3.4 Tables ............................................................................................................................ 64 3.5 Figures........................................................................................................................... 66 Chapter 4: Orthotropic Integrated Flow-Stress Model (IFS-Orthotropic)............................74 4.1 Biot’s Coefficient for Orthotropic System .................................................................... 76 4.1.1 Special Case of Isotropic Materials .......................................................................... 77 4.2 Rate of Change of Compressibility, 𝜿, for Orthotopic System..................................... 78 4.2.1 Special Case of 𝜿 for Isotropic Materials ................................................................. 81 4.2.2 Special Case of Solid Density Equation for Isotropic Materials .............................. 82 4.3 First Invariant of 𝛔′ ....................................................................................................... 83 4.3.1 Special Case for Isotropic Materials ......................................................................... 83 4.4 Evolution of Material Properties ................................................................................... 84 4.5 Summary of Governing Equations................................................................................ 85 4.5.1 Two-Phase Model ..................................................................................................... 85 4.5.2 Three-Phase Model ................................................................................................... 86 x  4.6 Summary ....................................................................................................................... 87 4.7 Tables ............................................................................................................................ 88 Chapter 5: Numerical Applications ...........................................................................................89 5.1 Assumptions .................................................................................................................. 89 5.2 Fully Solidified Material Example ................................................................................ 92 5.3 Two-Phase Examples .................................................................................................... 93 5.3.1 Unsolidified Isotropic Material - Mechanical Loading ............................................ 93 5.3.2 Curing Isotropic Composite Material ....................................................................... 96 5.3.3 Composite with Highly Compressible Fluid Phase under Thermal Loading ......... 100 5.3.4 Curing Orthotropic Composite Material ................................................................. 103 5.3.5 Unidirectional [0°] Angle Laminate Undergoing Compaction and Cure ............... 105 5.4 Three-Phase Examples ................................................................................................ 108 5.4.1 Debulking of a Column of Partially Saturated Composite ..................................... 108 5.4.2 Curing Isotropic Composite Material – Permeable B.C. ........................................ 110 5.4.3 Curing Isotropic Composite Material – Impermeable B.C. .................................... 112 5.4.4 Unsaturated Unidirectional Angle Laminate Undergoing Compaction and Cure .. 113 5.5 Tables .......................................................................................................................... 115 5.6 Figures......................................................................................................................... 118 Chapter 6: Summary, Conclusions, and Suggestions for Future Works .............................147 6.1 Summary and Conclusions ......................................................................................... 147 6.2 Suggestions for Future Works .................................................................................... 149 6.2.1 Incorporation of Fluids Phase Exchange ................................................................ 150 6.2.2 Full Form of Governing Equations ......................................................................... 150 xi  6.2.3 Adding More Phases to the Framework .................................................................. 150 6.2.4 Use of Coupled Flow Models ................................................................................. 151 6.2.5 Full Viscoelastic Material Models .......................................................................... 151 6.2.6 𝒖-𝑷 Formualtion ..................................................................................................... 151 6.2.7 Large Strain Formulation ........................................................................................ 152 6.2.8 Other Type of Elements .......................................................................................... 152 6.2.9 Solidification Factor Investigation .......................................................................... 153 Bibliography ...............................................................................................................................154 Appendices ..................................................................................................................................161 Appendix A Iterative Solution for System of Governing Equations ...................................... 161 A.1 Two-Phase System .................................................................................................. 161 A.2 Three-Phase System ................................................................................................ 165  xii  List of Tables  Table 3.1 Definition of matrices and vectors used in the two-phase isotropic FE model ............ 64 Table 3.2 Definition of matrices and vectors used in the three-phase isotropic FE model .......... 65 Table 4.1 Definition of specific matrices and vectors used in the two-phase orthotropic FE model that are different from those defined for the two-phase isotropic model (listed in Table 3.1) ..... 88 Table 4.2 Definition of specific matrices and vectors used in the three-phase orthotropic FE model that are different from those defined for the two-phase isotropic model (listed in Table 3.2) ...................................................................................................................................... 88 Table 5.1 Mechanical properties of the fibre, resin and atmospheric air .................................... 115 Table 5.2 Cure kinetics and viscosity models for 3501-6 resin and atmospheric air ................. 115 Table 5.3 Fibre-bed properties for AS4/3501-6 CFRP ............................................................... 116 Table 5.4 Fibre-bed Pseudo-Viscoelastic (PVE) formulation used for the 3501-6 resin (Zobeiry et al. [51]) .................................................................................................................................... 117 Table 5.5 Relaxation times and weight factor for 3501-6 resin for 𝝌𝟎 = 𝟎. 𝟗𝟖 (Kim and White [76])............................................................................................................................................. 117  xiii  List of Figures  Figure 1.1 Schematic of a controlled loop in an autoclave processing [1] ..................................... 6 Figure 1.2 Schematic illustration of resin polymerization during cure: a) pre-cure stage, b) start of curing as molecular size increases, c) continuous network at gelation, d) fully cured solid resin {{158 Berglunds, L. A. 1991; }} .................................................................................................... 6 Figure 2.1 Schematic showing the sub-model structure used in COMPRO [8] ........................... 26 Figure 2.2 Spring-piston analogy for effective stress [68] ........................................................... 27 Figure 2.3 CV-FEM Discretization of the domain ....................................................................... 27 Figure 2.4 Schematic showing an integrated model structures..................................................... 28 Figure 3.1 A representation of the two-phase composite system and its constituents .................. 66 Figure 3.2 A representation of axial and transverse flow of a single fluid through fibrous porous medium ......................................................................................................................................... 66 Figure 3.3 Evolution of the properties of the solid/solid-skeleton and fluid phases with increasing solidity of the fluid phase .............................................................................................................. 67 Figure 3.4 Conceptual example showing the effect of liquid resin pressure on the stress in the solidified resin after curing ........................................................................................................... 67 Figure 3.5 Schematic representation of 4-1 element .................................................................... 68 Figure 3.6 Algorithm for the finite element implementation of the two-phase IFS model .......... 69 Figure 3.7 Comparison of volume fraction equations based on small or large strain assumptions: a) expansion of fluid phase b) expansion of solid phase c) flow of the fluid into the system ...... 70 Figure 3.8 A representation of three-phase composite system and its constituents ...................... 71 xiv  Figure 3.9 A representation of decoupled multiphase fluids axial and transverse flow through fibrous porous medium ................................................................................................................. 71 Figure 3.10 Representation of three-phase composite system and the fluids flow mechanism as fluid 𝐹 solidifies ............................................................................................................................ 72 Figure 3.11 Evolution of the properties of the solid/solid-skeleton (𝑆/𝑆𝐾) and multiphase matrix phase (combined of fluid phases 𝐹 and 𝐹) with increasing solidity of the fluid phase ................ 72 Figure 3.12 Algorithm for the finite element implementation of the three-phase IFS model ...... 73 Figure 5.1 Assumed variation of the solidification factor with the degree of cure .................... 118 Figure 5.2 Assumed variation of solid and solid-skeleton properties as a function of the solidification factor (shown for bulk and shear moduli as example) .......................................... 118 Figure 5.3 Schematic of the miromechanics used in this study (a) two-phase (b) three-phase .. 119 Figure 5.4 A solid material under mechanical and thermal loading ........................................... 119 Figure 5.5 Comparison of the strains obtained using the current IFS model and Abaqus 4-noded plane strain element for fully solidified material ........................................................................ 120 Figure 5.6 Comparison of the stresses obtained using the current IFS model and Abaqus 4-noded plane strain element for fully solidified material ........................................................................ 120 Figure 5.7 A solid material under mechanical and thermal loading to investigate the effect of time step on the the IFS model predictions ................................................................................. 121 Figure 5.8 Load-displacement for top surface of the example shown in Figure 5.7 .................. 121 Figure 5.9 An unsolidified composite system subjected to uniaxial strain loading (a) permeable boundary condition at the top surface only, and (b) impermeable boundary condition at all surfaces ....................................................................................................................................... 122 xv  Figure 5.10 Predicted time-histories of the resin and air pressures for both examples shown in Figure 5.9(a) (permeable top surface) and Figure 5.9(b) (impermeable top surface) ................. 122 Figure 5.11 Predicted time-histories of the resin and air velocities for both examples shown in Figure 5.9(a) (permeable top surface) ......................................................................................... 123 Figure 5.12 Predicted time-histories of the sample vertical strain containing resin and air as matrix for both examples shown in Figure 5.9(a) (permeable top surface) and Figure 5.9(b) (impermeable top surface) .......................................................................................................... 123 Figure 5.13 Predicted time-histories of the resin and air volume fractions for both examples shown in Figure 5.9(a) (permeable top surface) and Figure 5.9(b) (impermeable top surface) . 124 Figure 5.14 Geometry and different boundary conditions of an isotropic (e.g. a unidirectional composite with fibre directions normal to the plane being analyzed) sample undergoing curing..................................................................................................................................................... 124 Figure 5.15 Possible location of boundary conditions for the examples shown in Figure 5.14 within a representative practical geometry ................................................................................. 125 Figure 5.16 (a) Time-history of applied temperature and predicted viscosity, (b) predicted degree of cure, and solidification factor of material for the case studies shown in Figure 5.14 ............ 125 Figure 5.17 Predicted time-histories of the resin pressure for all examples shown in Figure 5.14 .................................................................................................................................. 126 Figure 5.18 Predicted time-histories of the resin velocity in the vertical direction at the top surface for examples shown in Figure 5.14 with permeable top surface (two top cases) .......... 126 Figure 5.19 Predicted time-histories of the resin volume fraction for all examples shown in Figure 5.14 .................................................................................................................................. 127 xvi  Figure 5.20 Predicted time-histories of stresses for the examples shown in Figure 5.14 with resin matrix and constrained top surface (two left cases) ........................................................... 127 Figure 5.21 Predicted time-histories of vertical strain for the examples shown in Figure 5.14 with resin matrix and free top surface (two right cases) ............................................................. 128 Figure 5.22 Time-history of applied temperature and predicted viscosity of air ........................ 128 Figure 5.23 Predicted time-histories of the air pressure for all examples shown in Figure 5.14..................................................................................................................................................... 129 Figure 5.24 Predicted time-histories of the air velocity in the vertical direction at the top surface for examples shown in Figure 5.14 with permeable top surface (two top cases) .......... 129 Figure 5.25 Predicted time-histories of the air volume fraction for all examples shown in Figure 5.14 .................................................................................................................................. 130 Figure 5.26 Predicted time-histories of stresses for examples shown in Figure 5.14 with air matrix and constrained top surface (two left cases) .................................................................... 130 Figure 5.27 Predicted time-histories of vertical strain for free samples at top surface for examples shown in Figure 5.14 with air matrix and free top surface (two right cases) ............. 131 Figure 5.28 Geometry and boundary condition of orthotropic sample undergoing curing ........ 131 Figure 5.29 Final deformed shape of the orthotropic sample shown in Figure 5.28 for three different fibre angles (deformations are scaled by a factor of 10 for illustration purpose) ........ 132 Figure 5.30 Predicted time-histories of the resin velocity through the (a) right vertical surface (b) top horizontal surface for the example shown in Figure 5.28 ........................................... 132 Figure 5.31 Predicted time-histories of the (a) normal strain for 𝛽 = 0° and 𝛽 = 90° (b) shear strain for 𝛽 = 45° for the example shown in Figure 5.28 ....................................................... 133 xvii  Figure 5.32 Predicted time-histories of the stress for 𝛽 = 45° (a) normal stress (b) shear stress for the example shown in Figure 5.28 ............................................................................ 133 Figure 5.33 (a) Geometry and boundary conditions of the convex angle laminate; and (b) location of key points and the centre-line of the angle laminate ................................................ 134 Figure 5.34 (a) Time history of autoclave and part temperature and predicted resin viscosity (b) predicted degree of cure and solidification factor for the case study shown in Figure 5.33 ...... 134 Figure 5.35 Predicted time history of resin bulk and shear moduli based on pseudo-viscoelastic model for the case study shown in Figure 5.33 .......................................................................... 135 Figure 5.36 FE meshes used in the analysis of angle laminates ................................................. 135 Figure 5.37 Predicted deformed shape of the laminate for the case study shown in Figure 5.33 (deformations are scaled up by a factor of 2 for illustration purposes) ...................................... 136 Figure 5.38 Photograph of the cross-section of [0°] laminate on a convex tool with permeable boundary (bleeding system) showing the corner thickening (with permission from Hubert [1])..................................................................................................................................................... 136 Figure 5.39 Predicted deformed shape of a [0°] laminate on a convex tool with permeable boundary (bleeding system) at end of compaction using flow model (with permission from Hubert [1]) .................................................................................................................................. 137 Figure 5.40 (a) Normal displacement of critical points A, B, and C through the curing process (b) tangential displacement (in 𝜉 direction in Figure 5.37) for different through thickness position (in 𝜂 direction in Figure 5.37) of the laminate edge at different times for the case study shown in Figure 5.33 .................................................................................................................................. 137 Figure 5.41 Predicted fibre volume fraction, (1 − 𝜑), distribution at different instances of time for case study shown in Figure 5.33 ......................................................................................... 138 xviii  Figure 5.42 Predicted distribution of the resultant (a) axial force (b) bending moment, along the mid-plane of the laminate at three different times for the case study shown in Figure 5.33 .... 138 Figure 5.43 Geometry and boundary conditions of an unsolidified-partially saturated column representing a slice through the thickness of an unidirectional laminate with fibres transverse to the x-z plane shown and subjected to mechanical loading at the top surface ............................. 139 Figure 5.44 Predicted distribution of pressure along the partially saturated column at different time for the case study shown in Figure 5.43 ........................................................................... 139 Figure 5.45 Predicted distribution of air volume fraction along the unsaturated column at different times for the case study shown in Figure 5.43 ........................................................... 140 Figure 5.46 Predicted distribution of air volume fraction along the unsaturated column at different time for the case study shown in Figure 5.43 ............................................................ 140 Figure 5.47 Predicted time-histories of the resin velocity for the case study shown in Figure 5.14 with permeable and free top surface (top right figure) ............................................ 141 Figure 5.48 Predicted time-histories of the air velocity for the case study shown in Figure 5.14 with permeable and free top surface (top right figure) ............................................................... 141 Figure 5.49 Predicted time-histories of the strain in vertical direction for the case study shown in Figure 5.14 with permeable and free top surface (top right figure) ....................................... 142 Figure 5.50 Predicted time-histories of the resin and air volume fractions for the example shown in Figure 5.14 with permeable and free top surface (top right figure) ....................................... 142 Figure 5.51 Predicted time-histories of the stress in horizontal direction for the example shown in Figure 5.14 with permeable and free top surface (top right figure) ....................................... 143 Figure 5.52 Predicted time-histories of the strain in the vertical direction for the example shown in Figure 5.14 with impermeable and free top surface (bottom right figure) ............................. 143 xix  Figure 5.53 Predicted time-histories of the resin and air volume fractions for the example shown in Figure 5.14 with impermeable and free top surface (bottom right figure) ............................. 144 Figure 5.54 Initial and final average volume fraction of a) each phase and b) fibre and combined resin and air for the example shown in Figure 5.33 .................................................................. 145 Figure 5.55 Final distribution of a) air and b) fibre volume fraction over the domain for 𝜑0𝑎 =0.3 for the example shown in Figure 5.33 ................................................................................ 146 Figure 5.56 Normal displacement for different values of initial air content a) key point A, b) key point C for the example shown in Figure 5.33 ......................................................................... 146  xx  Acknowledgements  I offer my sincere gratitude to my supervisors, Professor Reza Vaziri and Professor Anoush Poursartip, for providing valued insight, constructive criticism, and paternal support during my PhD studies. This thesis would not have been possible without their professional and patiently supervision.  I would like to present my enduring gratitude to my friend, Dr. Alireza Forghani, for his true friendship, positive encouragement, constructive supervision, and willingness to share his endless scientific expertise. His role in technical side of my research was nothing less than an official supervisor. I would also like to thank my colleagues in Composites Research Network (CRN), Kamyar Gordnian, Mina Shahbazi, Andrew Stewart, Amir Hossein Afrasiabi Garekani, Mehdi Haghshenas, Sardar Malek Mohammadi, Hamidreza Bakhtiarizadeh, and Casey Keulen for their helpful discussions and cooperation during my research. Many thanks are owed to my supportive friends, Ashkan Babaie, Amirhossein Hadi Hossein Abadi, Hamid Moradi Kashkouli, Reza Nickmanesh, and Moein Javadian who have always been there for me. For parents, nothing is harder than absence of their children. I would like to express my sincere appreciation to my exceptional parents and my brother for their patience during my absence and their unconditional love, continual support, and encouragement throughout my whole life. And last but not least, my special thanks go to my lovely wife whom this thesis, more than all, is indebted to. Words are unable to express the role of her kindness, courage, intelligence, and pure love from start to the end of my PhD.  Sina Amini Niaki The University of British Columbia September 2017  xxi  Dedication     To my pillar of strength,  my dad, Mansour peace of my life,  my mom, Simin hope of my life,  my brother, Aria and happiness of my life,  my wife, Naghmeh  1  Chapter 1: Introduction Fibre-reinforced polymeric composite materials have attracted increasingly interest over past decades. Traditional materials such as steel and aluminum are being replaced with lightweight composite materials even on critical structural components such as wings and fuselage of airplanes because of their distinct properties. Applications are generally included in advanced key economic sectors such as renewable energy, aerospace, and ground transportation. Composite materials show excellent strength- and stiffness-to-weight ratio as well as impressive durability. Also, manufacturing of composite structures enables production of complex-shaped large-scale structures from the raw materials in only one step. However, their advantageous manufacturing process involves considerable risk if it is not performed carefully along with scientific knowledge about their behaviour during production. In order to predict process-induced defects (e.g. residual stresses, shape distortions, voids, etc.) and the resulting variations in material properties, it is important to have predictive models that can capture the behaviour of the composite structure as it undergoes the complete sequence of events during manufacturing. An accurate process simulation tool, that incorporates the interactive multi-physics phenomena involved during processing, can help mitigate the manufacturing risks and associated costs. However, the complex behaviour of composite materials as a result of the evolving properties of their various constituents poses challenges to their predictive modelling during processing.  Different types of processing are used for manufacturing of composite parts depending on reinforcement architecture (e.g. unidirectional, woven fabric, etc.), type of matrix and reinforcement materials, industrial application of the composite material (e.g. aerospace, automation), and so forth. In some cases, pre-impregnated sheets of unidirectional fibres 2  (prepreg) are stacked at different orientations on a tool having the shape of final products. The lay-up process is typically followed by a debulking process at room temperature to remove the entrapped volatiles. In an alternative manufacture processing, dry mat or textile of reinforcement is impregnated with liquid resin flowing through the reinforcement. In either case, the manufacturing is followed by subjecting the whole part to a controlled cycle of temperature and pressure outside or inside an autoclave (Figure 1.1) that cures the polymeric resin and forms a solid composite material with all desired properties.   1.1 Motivations  In processing of composite materials, the presence and interaction of constituents with diverse properties and physical forms, such as solid, liquid and gas, makes their behaviour very complex to be fully understood and accurately predicted. Specifically, the presence of the gas phase and its highly compressible behaviour compared to the solid and liquid constituents, gives rise to difficulties in process modeling of composite materials. Also, extensive physical and chemical changes occurring simultaneously during the course of process add more complexity into the process simulation. For example, as the resin cures, it undergoes very diverse behavioural regimes such as viscous (low molecular weight monomers), rubber-like (continuous networks of molecules), and glassy (cross-linked macromolecular structures) behaviour at different stages of processing (Figure 1.2). Distinct formulations associated with each stage of material, liquid, gel, or solid, are applied at different times during composite materials processing. Also, these extensive transformations of the polymeric material does not occur uniformly in a composite part during processing which adds further complication to its modelling.   3  In composite materials, presence of voids and spaces that are filled with gas, represent an important manufacturing quality issue. Porosity, defined as volume fraction of voids, is considered to be a process-induced defect that can potentially lead to deterioration of the structural performance of the final part. Availability of modeling tools that can quantify the total amount and spatial distribution of porosity trapped in final composite structures is highly desirable in order to take required precautions to minimize it or limit its occurrence to locations that are far removed from critical regions of the structure. This requires a detailed understanding of composite materials and their constituents’ behaviour during such processes.  Accurate modelling of composite materials behaviour during different types of processing requires a unified model that takes into account all the diverse physical and chemical properties of constituents as they evolve. Such integrated process models provide the fundamental capability to model various types of composite materials processing.   1.2 Objectives and Outline of Thesis The main objective of this thesis is to develop a framework to integrate flow and stress development in processing of composite materials regarded as a three-phase material system. The three-phase system, involving liquid phase (e.g. resin), gas phase (e.g. air or water vapor), and solid phase (e.g. fibre), considers compressibility of all phases with no restriction on their level of compressibility.    Integrating flow and stress development in a single formulation needs to capture the behaviour of the material at the two extremes of processing, namely fluidic behaviour of the resin at early 4  stages of the process and the final solid behaviour of the composite material, as well as the intermediate behaviour between them to enhance the efficiency and accuracy of the predicted results. Also, having the capability to account for highly compressible materials is imperative where the very high compressibility of the gas phase plays a significant role in the system response. The framework for the 3-phase Integrated Flow-Stress (3IFS) model will be developed in this study such that it provides the capability of including further details and mechanisms in subsequent studies.   The thesis is organized as follows:  Chapter 2: Literature review - A brief literature review of composite process models is presented. The previous fundamental studies on flow through porous media and stress development in composite systems and their main formulations, such as effective stress and flow equations, are presented briefly. Some of the relevant studies on three-phase modelling of soils are reviewed. Finally, a brief review of literature on porosity investigations in composite processing is presented.  Chapter 3: Isotropic Flow-Stress Model (IFS-Isotropic) - The basis of the IFS model is presented in this chapter for composite materials with individual phases possessing isotropic properties. First, the two-phase (2IFS-Isotropic) model is presented through introduction of the governing conservation laws and constitutive equations when the fluid phase solidifies during processing. Then, the 2IFS model is extended to a three-phase (3IFS-Isotropic) model. The finite element (FE) formulation and the solution method is presented for both the two-phase and three-phase models. 5   Chapter 4: Orthotropic Flow-Stress Model (IFS-Orthotropic) – The IFS-Isotropic model presented in the previous chapter is extended to the more general case of orthotropic composite materials in this chapter. The IFS formulation is generalized for the condition when both the solid and solid-skeleton possess orthotropic properties. First, the orthotropic version of the key parameters in the IFS model are derived and then the governing equations are presented for both 2IFS-Orthotropic and 3IFS-Orthotropic cases.  Chapter 5: Numerical Applications – First, preliminary and simple examples are provided for the isotropic model that serve to verify the accuracy of the developed integrated model compared to previous independent modules. Then, the IFS-Orthotropic model is applied to process modeling of a configured L-shaped composite laminate as a representative case study for assuring the performance of different process models and results are compared qualitatively to previous experimental investigations.   Chapter 6: Summary, Conclusions, and Suggestions for Future Research – A summary and conclusion of this work as well as suggestions to improve the current 3IFS model and recommendations for future works are presented in this chapter.  6  1.3 Figures  Figure 1.1 Schematic of a controlled loop in an autoclave processing [1]                          (a)                                          (b)                                      (c)                                       (d)  Figure 1.2 Schematic illustration of resin polymerization during cure: a) pre-cure stage, b) start of curing as molecular size increases, c) continuous network at gelation, d) fully cured solid resin {{158 Berglunds, L. A. 1991; }}  ControllerPressureSourceVacuumPumpHeater( T )( )( )AutoclaveCompositein vacuum bagSensoroutputController output7  Chapter 2: Literature Review Processing of composite materials involves various phenomena such as flow of resin and gas through the fibre-bed, thermochemical changes and phase transformations during resin cure, and the resulting build-up of residual stress and dimensional variations in the final part. These phenomena are usually simulated through the application of independent sub-models. This approach, called sequential method or integrated sub-models approach, was first applied to composite processing by Loos and Springer [2] and later by Bogetti and Gillespie [3,4] and White and Hahn [5,6]. In the sequential method, several modules such as cure kinetics, thermal and heat transfer analysis, resin flow, stress development, etc., are carried out independently and sequentially where outcome at the end of each module is used as initial condition for the next module. Johnston et al. [7,8] developed a multi-physics 2D plane strain finite element code, called COMPRO, for polymeric composites. In COMPRO with a schematic of its flow chart shown in Figure 2.1, thermochemical, flow, and stress modules are used sequentially where resulting geometry, volume fraction of phases, degree of cure, and/or temperature (if heat transfer analysis is done) at the end of each module is used as initial state for subsequent step which is performed during the whole process cycle. Arafath [9] later introduced the COMPRO Component Architecture (CCA) into the commercial finite element software, Abaqus.  In the sequential method, the interaction between the resin flow and stress development is not captured despite the fact that these phenomena occur concurrently in processing of composite materials. Furthermore, the history of pressure during the flow regime of processing is ignored at the start of the stress development module which potentially results in inaccurate prediction of the final process-induced stresses and deformations in the composite material. Mapping of the 8  results from one state to another in this manner is a tedious and inefficient process from a computational standpoint. Moreover, in reality, not all regions of a large and complex composite part undergo the same curing process and, consequently, resin gelation does not occur uniformly over the whole spatial domain of the structure. In a non-interactive sub-model approach, however, resin gelation can only be considered to occur simultaneously for the whole composite structure at the start of the stress module and therefore local spatial and temporal variations in gelation cannot be considered in such an approach.  On the other hand, if flow and stress modules are integrated into a unified formulation, a non-uniform gelation (curing) of the resin can easily be accounted for when combined with a heat-transfer analysis.   In this chapter, a literature review on the previous works on flow and stress models is first presented. This is then followed by a discussion of integrated approach applied to processing of composite materials. The most recent studies on porosity and gas transport investigations in processing of composite materials are then reviewed briefly. Finally, research objectives for this work are stated at the end of the chapter.  2.1 Flow Models Flow models predict movement and migration of the fluid phase, liquid or gas, through the solid porous media. Flow models typically incorporate governing equations including mass conservation, equilibrium equation of the fluid phase(s), and equilibrium equation of the system. Some of these models predict flow during compaction of composite materials in consolidation modelling of composites processed outside an autoclave, i.e. in Out-Of-Autoclave (OOA) processing, or within an autoclave under high pressure and temperature. Other category of flow 9  models predict flow of the fluid phase in Liquid Composite Molding (LCM) processing where the liquid phase impregnates the dry fibre reinforcements (preforms). Although most of the developed compaction models deal primarily with saturated flows (single phase flow), some of the liquid molding models have succeeded in predicting voids formation during the resin infusion process.   2.1.1 Composite Materials Compaction Effective stress model was introduced in soil mechanics in 1940s for consolidation modeling of soils. Terzaghi [10] introduced the effective stress formulation for a fluid flowing through a porous structure of an incompressible solid phase: 𝜎𝑖𝑗 = 𝜎𝑖𝑗′ −  𝑃𝛿𝑖𝑗 (2-1) where the subscripts 𝑖 and 𝑗 denote coordinate directions aligned with 𝑥1, 𝑥2, and 𝑥3 axes, 𝛿𝑖𝑗 is the Kronecker delta (or identity tensor), 𝜎 is the total stress in the system (including both the fluid and solid phases), 𝑃 is the pore pressure of the fluid phase and 𝜎′ is so-called effective stress of the solid-skeleton. The spring-piston system shown in Figure 2.2 represents a mechanical analogy for this equation where the fibre-bed is assumed to be a spring made of an incompressible material in a cylinder full of resin. The total external load is shared between the fibre-bed, as stress, and resin, as pressure, each taking a fraction of the load depending on the flow condition of the resin and deformation of the fibre-bed [1].   In order to take into account the compressibility of the phases, the flow model developed by Biot et al. [11,12] and its extension by Zienkiewicz and Shiomi [13] can be considered. Biot et al. [11,12] developed equations to model consolidation of the porous media taking into account the 10  compressibility of the phases through the introduction of some physical constants for the porous media. Zienkiewicz and Shiomi [13] presented an alternative form of the Biot’s model extended to the dynamic behaviour of saturated porous media. They demonstrated the effectiveness of the model for slow compaction as well as shock excitation scenarios. In Biot’s effective stress formulation, a coefficient termed Biot’s coefficient takes on a value between zero and unity, determines the share of pore pressure contributing to the overall load carrying capacity, and represents deformability of the porous media. The Biot’s effective stress formulation for a system made up isotropic phases is given as: 𝜎𝑖𝑗 = 𝜎𝑆𝐾𝑖𝑗 − 𝑏 𝑃𝛿𝑖𝑗 (2-2) where subscript 𝑆𝐾 denotes solid-skeleton and the scalar parameter 𝑏 is the Biot’s coefficient. For isotropic materials, Biot’s coefficient can be derived in terms of the bulk moduli of the solid-skeleton and the solid phase, 𝐾𝑆𝐾 and 𝐾𝑆, as follows:  𝑏 = 1 − 𝐾𝑆𝐾/𝐾𝑆 (2-3) where the subscript 𝑆 denotes solid phase.   The work by Gutowski et al. [14] can be considered as one of the pioneering studies that applied such theories to consolidation of composite materials during curing. Using experimental measurements on compaction of fibre-beds, they applied Terzaghi’s effective stress formulation to the consolidation of composite materials.   Hubert et al. [1,15] developed a flow model for processing of composite materials with constituents that were assumed to be incompressible.  They used Terzaghi’s effective stress theory (Eq. (2-1)) for load sharing between the resin and the fibre-bed, so deformability of the 11  solid grains within the porous medium was not taken into account. For flow of fluid through porous media, Darcy’s law [16] was employed as the continuum (macroscale) representation of the physics of flow in porous media. Darcy’s law is valid for laminar flows through porous media which can be shown to be the case even for low viscosity fluids (such as gaseous volatiles) flowing through the miniscule space between the fibres in composite materials. Darcy’s law is the macroscopic or smeared equivalent of the microscale approach where the solid is explicitly represented at the pore-scale and the general equations of fluid mechanics (e.g. the Stokes’ equations) are solved directly subject to appropriate boundary conditions. In the continuum approach that forms the basis of Darcy’s law, each point of the domain contains both fluid and solid in an average sense and the average velocity and pressure quantities are governed by the following equation: 𝑃,𝑖 + 𝜇 𝑆𝑖𝑗−1 𝑣𝑗 = 0 (2-4) where comma denotes spatial derivative i.e. ( ),𝑖 =𝜕𝜕𝑥𝑖( ), 𝜇 is the fluid viscosity and 𝑆𝑖𝑗 is the permeability tensor of the solid-skeleton. Permeability describes the resistance of the porous solid-skeleton to the fluid flow in an average sense (or effective property). Considering Darcy’s equation as the equilibrium equation of the fluid phase as well as employing Terzaghi’s effective stress formulation in the equilibrium equation of the system along with mass conservation equation of the system, compaction of composite materials during cure was modelled in the flow module of a multi-physics 2D plane strain finite element code, called COMPRO [8].    Haghshenas [17] used the Darcy-Brinkman [18] equation for momentum balance of the fluid phase where the shear stress terms of the fluid phase are taken into account:  12  𝑃,𝑖 + 𝜇 𝑆𝑖𝑗−1 𝑣𝑗 + 𝜏𝑖𝑗,𝑗 = 0 (2-5) In Eq. (2-5), 𝜏 represents the shear stress in the fluid phase and is given by: 𝜏𝑖𝑗,𝑗 = [𝜇(𝑣𝑖,𝑗 + 𝜑?̇?𝑖,𝑗) + 𝜇(𝑣𝑗,𝑖 + 𝜑?̇?𝑗,𝑖)],𝑗 (2-6) where 𝜑 is the volume fraction of the fluid phase, 𝑢 is the displacement of the system (associated with the solid structure), and overdot denotes time derivative. In most applications, the shear stress can be ignored since it is much smaller than the drag force represented by 𝜇 𝑆𝑖𝑗−1 𝑣𝑗.  2.1.2 Liquid Composite Molding (LCM) Liquid composite molding involves the saturation of dry-fibre reinforcement with liquid resin by means of a pressure gradient. Fibre preform is placed in a mold, then the relatively viscous liquid resin is injected through porous preform until it is saturated. For simple molds, flow pattern can be predicted based on simple closed form equations. However, an advanced numerical flow simulation is often needed for flow simulation in complex molds. These numerical simulations involve prediction of resin flow front and pressure distribution during/after resin injection.   Wang et al. [19], Fracchia et al. [20] and Bruschke and Advani [21] developed a Control Volume-Finite Element Method (CV-FEM) method for Resin Transfer Molding (RTM) process. In the CV-FEM technique, the spatial domain is discretized with a FE mesh where several sub-volumes are created within each element through connecting the centroids of the elements to the midpoints of elements sides. The generated sub-volumes from adjacent elements form the control volumes (Figure 2.3). Each of the constructed CVs surrounds one of the nodes in the FE mesh. For this reason, FE nodes are sometimes called nodal control volumes as well. In this 13  technique, moving of the resin front is tracked by monitoring the state of each nodal control volume. A quantity taking values between zero and unity, called fill factor, represents the ratio of fluid volume in the CVs to the total volume of CVs.  Loos and MacRae [22] used CV-FEM method for two-dimensional resin flow simulation in RTM process of a single blade-stiffened panel. In their model, the textile preform is considered as a heterogeneous and anisotropic porous medium. The resin is assumed to be incompressible, and capillary and inertia effects are neglected (low Reynolds number flow). Darcy’s law and continuity equation are assumed for the flow within the FE domain. The FEM discretization is used to obtain pressure distribution in the saturated domain. After calculating the pressure distribution in the saturated domain, flow velocity, assumed to be constant within each element, is determined at the centroid of each element using Darcy’s law.  Simacek and Advani [23] and Kuentzer et al. [24] extended CV-FEM technique to be used for a dual-scale porous media. Darcy’s law was assumed for flow of resin through pores between the tows while a separate one-directional linear flow was considered for flow through the pores inside the tows. Combining continuity equation and Darcy’s law for flow at the macro-scale, the effect of flow at the micro-scale is introduced in the governing equation as a sink term representing flow into the tows. Using this method, they were able to consider the partially saturated region between the macroscopic flow front and the fully saturated region. This method was implemented as Liquid Injection Molding Simulation (LIMS) program and used successfully for industrial LCM applications.  14  Tan and Pillai [25] added more complexity to the previously developed dual-scale models. Instead of considering the micro-scale medium as only extra volume and solving the fluid movement within this medium as one-dimensional linear flow, a complex two-dimensional nonlinear flow was considered for the fluid movement within the tows. Flow inside the tows was modelled in parallel to the macro-flow through a two-dimensional micro-scale FE simulation while the solution within the two scales were coupled through the fluid pressure at the boundaries of the two-scales.  Trochu et al. [26,27] used non-conforming elements in CV-FEM in order to skip the need for separate CVs in their model. The domain was discretized with triangular elements with nodes at the middle of the element borders. This numerical scheme satisfies conservation of the fluid mass across the element boundaries. Therefore, the elements themselves could be considered as CVs. Computing pressure distribution from FE model, fluid velocity on the flow front was computed using Darcy’s law. As a result, volume of resin flowing across each element during a given time step could be calculated.   In an alternative approach to CV-FEM method, Mohan et al. [28,29] used a pure FE based method to simulate resin flow through porous media. In this approach, the degree of saturation (fill factor) was included in the FE formulation as a degree of freedom. Hence, no identification of CVs was required.   15  2.1.3 Three-Phase Flow Models Composite materials usually undergo a debulking process at room temperature that involves flow and removal of entrapped air and other gaseous volatiles. Furthermore, compaction of composite materials during the curing process requires consideration of the presence and flow of residual gas from previous steps, as well as gases produced during the in- or out-of-autoclave processing. For simulation of such processes it is essential to account for the gas phase and its flow in the formulation of the governing equations. The evolution of porosity during the curing cycle as influenced by gas flow, thermal effects, and resin hardening, necessitates consideration of a three-phase integrated flow-stress model. Despite extended experimental investigations in porosity and gas transport phenomena of composite materials, a multiphase modelling approach that explicitly considers the gas phase in composites processing is still not available and numerical modelling of voids formation in composites process modeling is limited to LCM processes. For soils on the other hand, models have been developed considering the gas phase in the three-phase solid-liquid-gas models. For gas transport simulation in composite materials, the same principles of multiphase fluid flow previously utilized in modelling partially saturated soils can be adopted for simulating uncured composite materials undergoing manufacturing process.  In a two-phase system (consisting of a solid and a single fluid phase), it is assumed that all the fluid pressure in the system is transferred directly to the solid phase as applied in Eqs. (2-1) and (2-2). In a three-phase system, the matrix is a mixture of two fluids, denoted by subscripts 𝐹 and ?̂?. Consequently, it is assumed that the pressure of the equivalent single phase matrix, denoted by subscripts 𝑚, is transferred to the solid grains. In an unsaturated system including multiple fluid phases (e.g. liquid and gas), Bishop [30] proposed that the pore pressure of a mechanically 16  equivalent single phase matrix be used as representative of the multiphase matrix in effective stress formula: 𝜎𝑖𝑗 = 𝜎𝑆𝐾𝑖𝑗 − 𝑏𝑃𝑚𝛿𝑖𝑗 (2-7) where 𝑃𝑚 is the effective pore pressure of the homogenized matrix and defined as: 𝑃𝑚 = 𝑃?̂? −ℬ(𝑃?̂? − 𝑃𝐹) (2-8) and 0 ≤ ℬ ≤ 1 is the so-called Bishop’s parameter. It is assumed that the Bishop’s parameter is solely dependent on the degree of saturation of the fluid 𝐹, which is defined as: 𝜓𝐹 =𝜑𝐹𝜑𝐹 + 𝜑?̂? (2-9) with 𝜓𝐹 + 𝜓?̂? = 1 [31] where 𝜑 represents volume fraction. Different experimental attempts have been made to find the dependency of the Bishop’s parameter on the degree of saturation (as reported in [31]). The simplest choice made by Schrefler [32] is ℬ = 𝜓𝐹. This choice of the Bishop’s parameter yields the matrix pore pressure as the volumetric weighted average of the pressure of its constituent fluids:  𝑃𝑚 = 𝑃?̂? − 𝜓𝐹(𝑃?̂? − 𝑃𝐹) (2-10)  In a macroscopic representation of the fluids flow in a three-phase system, the momentum exchange across the interface of the two fluid phases is assumed to be negligible compared to the momentum exchange between each fluid phase and solid phase. This is because from a practical standpoint, the significance of the flow coupling between the two fluid phases is questionable when the aim of the modelling is not to capture bubble movements [33]. Using the decoupled flow model for fluids, the averaged momentum balance equation, in the form of Darcy’s equation, can be written separately for each fluid phase as [33] 17  𝑃Φ,𝑖 + 𝜇Φ𝑆Φ𝑖𝑗−1𝑣Φ𝑗 = 0 (2-11) where Φ ≡ 𝐹 or ?̂?, and 𝑆Φ𝑖𝑗 is the effective permeability tensor of the solid-skeleton with respect to the fluid Φ (𝐹 or ?̂?). Effective permeabilities of the solid-skeleton with respect to each fluid phase in Eq. (2-11) are related to the absolute permeability of the skeleton as: 𝑆Φ𝑖𝑗 = 𝑆Φ̅𝑆𝑖𝑗 (2-12) where 𝑆 represents the absolute permeability tensor of the solid-skeleton, and 𝑆Φ̅ is a scalar which takes on values between zero and unity and is a measure of relative permeability to fluid Φ. Bear and Bachmat [34] showed that the relative permeabilities can be taken to be only dependent on degree of saturation of each phase, 𝜓𝐹 and 𝜓?̂?. Complicated functions for dependency of the relative permeabilities on degree of saturation can be incorporated as listed in [33]. However, relative permeabilities are assumed to vary linearly with respect to the degree of saturation for each phase i.e. 𝑆?̅? = 𝜓𝐹 and 𝑆?̂̅? = 𝜓?̂? [35]. In other words, in a macroscopic model, permeability of each fluid phase is taken to be proportional to the area available for them to flow through the solid-skeleton.   In processing of composite materials, pressure difference across the interface of two immiscible fluid phases shows significant changes as a result of extensive physical and chemical evolution in the constituents. A diverse set of formulations and values have been reported in the literature for capillary pressure between the liquid and gas phase in soils or composite materials. As an example, Leverett function [36,37] is a very appropriate function that considers dependency of capillary pressure on volume fraction of the matrix, absolute permeability, state of curing (i.e. surface tension of fluid phase and contact angle) and degree of saturation 18  𝑃𝑐 = √𝜑𝑆𝛾𝑓′(𝜃) 𝐽(𝜓𝐹) (2-13) where 𝛾 is surface tension between the two fluid phases, 𝑓′(𝜃) is a function that accounts for the  dependency of 𝑃𝑐 on contact angle 𝜃 for solid-fluid-fluid system where for straight capillary tubes cos 𝜃 can be used for the contact angle function [37,38]. 𝐽(𝜓𝐹) is also an empirical function called Leverett’s J-function and considers dependency of capillary pressure solely on degree of saturation 𝜓𝐹. Surface tension of different resins, as well as contact angle of different fibre-resin-air systems, have been widely investigated experimentally [39-43]. In some cases, their dependency on temperature, degree of cure, or time is reported [41-43]. Leverett’s J-function also is studied meticulously in [37] for different materials, matrix volume fractions, and fibre architectures.   Khoei and Mohammadnejad [44] presented a three-phase finite element model for seismic analysis of earth and rockfill dams. In their formulation the displacements of the solid phase, the pressure of the wetting phase and the capillary pressure (pressure difference between the two fluid phases) are taken as the primary unknowns while the degree of saturation (fraction of wetting fluid phase in the pore space) is calculated in their model by means of the capillary pressure-saturation relationship from the experiments. In their model, the Bishop’s effective stress is also considered to account for compressibility of soil grains.   Oettl et al. [45] used a three-phase deformable finite element model for dewatering simulation of fully or partially saturated soils using compressed air. Using Bishop’s effective stress formulation, equilibrium of the multiphase system and mass balance for each fluid phase, 19  combined with momentum balance equations, are taken as governing equations of the model. Displacement of the porous medium and pressures of the individual fluid phases are the unknowns of the FE model.   Nagel and Meschkle [46] presented a FE formulation for partially saturated deformable soils as a three-phase model where the liquid phase was assumed to be incompressible. The governing equations of the model were based on the mass conservation of fluids combined with equilibrium of the three-phase system along with the relationship between the capillary pressure and the degree of saturation, Darcy’s law for two-phase flow, and the effective stress formulation.   Schümann [47] also showed that the built-in module for two- or three-phase models in the commercial finite element code, Abaqus, can be extended through implementing a user element (UEL) subroutine for unsaturated soils based upon the three-phase model developed by Holler [48].  2.2 Stress Models Residual stress and process induced deformations are important parameters in processing of composite materials. Residual stresses can potentially lead to deterioration of the structural performance of the final part as a result of reduction in the ultimate strength and fatigue life. Induced deformations also have significant impact on the application of the final composite part as well as assembly of different composite parts into larger structures.  20  Bogetti and Gillespie [3,4] developed an anisotropic cure and stress model for two-dimensional analysis of thick thermoset composites. Accounting for thermal and chemical interactions associated with cure, temperature and degree of cure distributions are predicted and correlated against experimental measurements. A methodology for evolution of process-induced stresses as well as associated deformation in thick composite laminates was proposed. The developed heat transfer and cure analysis was combined and coupled with an incremental laminate plate theory.  Johnston et al. [7,8,49] developed the stress development module which was coupled with other modules, such as temperature, degree of cure, and flow, within COMPRO. The polymeric resin was modelled as a Cure-Hardening Instantaneously Linear Elastic (CHILE) material model which assumed the linear elastic moduli of resin increase instantaneously and monotonically with the progression of cure. Calculated elastic constants for the resin were used along with properties of the fibres in micromechanics models to determine the composite ply properties.   Zobeiry et al. [50-52] extended the material model used in COMPRO to a fully viscoelastic transversely isotropic material model. The results from fully viscoelastic model were compared against the special case of pseudo-viscoelastic material model for processing of composite with different cure cycles. The accuracy of the results from either of the models was investigated and the applicability of each model was reported.  The pseudo-viscoelastic material model was proven to yield the same results as a full viscoelastic model for standard cure cycles (one-hold or two-hold cure cycles without post-hardening) while for non-standard cure cycles (e.g. involving post cure) use of a full viscoelastic model was recommended.  21  Gigliotti [53] used a thermo-chemically-elastic model in a finite element framework to predict residual stresses developed after the gel point during processing of composite materials. Also, Lacoste et al. [54] incorporated a multi-scale approach for in-plane stress prediction as a result of thermal and cure shrinkage during processing of isotropic composite materials. Incorporating a two-step self-consistent scale transition for homogenization of composite materials, the stress at the microscopic level due to mismatch of properties between the fibres and the matrix was also computed for a simplified cure process.   2.3 Integrated Approaches  Incorporating the above mentioned main modules of sequential method into a unified model that captures the various phenomena involved during processing in an integrated and continuous process modelling would help overcome the foregoing drawbacks of the sequential approach. A schematic of a full integrated model flow chart is shown in Figure 2.4 (which integrates all different module shown in Figure 2.1). Haghshenas [17] developed a framework that integrates both resin and stress development into a unified process model for composite materials. In this IFS model, modifications were made to the classical flow through porous media models (e.g. Biot’s flow model) to enable a seamless connection to the regime of the process simulation when the composite material is in its fully cured (solid) state. This integrated model was successfully applied to the case studies with practical applications when typical boundary conditions are considered. However, although the above two-phase (resin-fibre) IFS model considers compressibility of its constituents, the accuracy of the model during the flow regime, especially for drained (permeable) condition, is sacrificed due to using micromechanics formulations for composite materials which inherently assume that the resin cannot leave the system due to flow. 22  For this reason, the model breaks down when it is applied to a composite system with constituents that have extreme bulk properties (e.g. very low bulk moduli or very high coefficients of thermal expansion). Also, for using such an approach in the transversely isotropic IFS model, a switch (instantaneous change in the formulation) had to be introduced in the model at the nominal point of resin gelation. As an undesirable consequence of having this switch, two different formulations were considered for each stage of processing and therefore a continuous formulation for transversely isotropic composites was not achieved. Furthermore, using this approach to integrate the flow and stress build-up, the fibre (or fibre-bed) free strain before gelation had to be ignored in the model which in some extreme cases led to unrealistic predictions of residual stresses at the end of the process.  2.4 Porosity Process-induced porosity in composite parts can affect their mechanical performance drastically. An approximate 20% reduction in interlaminar and flexural shear strength can be seen because of a 2% increase in the porosity content [55]. Porosity is affected by multiple phenomena involved in processing such as gas transport and resin thermochemical evolutions which makes producing high quality composite parts with a low porosity content highly challenging.   Multiple sources of porosity formation in composite materials as well as several alternatives to bring it under an acceptable level have been discussed in the literature. The work by Kardos et al. [56] can be considered as one of the pioneer work on porosity in prepreg composites. Bubble growth as a result of moisture dissolution was studied and the effect of temperature and pressure 23  was investigated to propose a pressure-temperature-humidity stability map to identify the condition for porosity formation during the cure cycle.  Boey et al. [57,58] studied the effect of high pressure in an autoclave as well as vacuum bagging in order to reduce the void content in thermoset composites. The bending strength and modulus of composite parts manufactured under different process conditions were also analyzed where an increase in the strength was reported as a result of post-curing.   Thorfinnson and Biermann [59,60] investigated processing parameters that affect the laminate void content and found that controlling the degree of resin impregnation was important as partial resin impregnation of the prepreg creates gas transport pathways within the laid-up laminates.  Arafath et al. [61] proposed a simple one-dimensional gas transport model for composite materials. Assuming Darcy’s law and ideal gas law, a closed form solution was developed to find the fraction of gas that is extracted from the composite during processing.   Farhang and Fernlund [62,63] performed a characterization study of void and porosity evolution in prepregs undergoing OOA processing. Using experimental measurement of gas permeability and porosity, an understanding of gas transport in OOA prepreg processing was developed by examining the time scales for gas transport by Darcy flow and molecular diffusion.  24  Wells [64] studied the resin voids and their growth mechanisms during processing of composites. The processing conditions affecting surface porosity were investigated and correlated to bulk porosity inside the composite product.   Kay et al. [65,66] investigated the gas transport mechanisms in processing of composite materials to develop a simple model for formation and transport of voids in the composite materials. Experimental measurements were conducted to determine the connection between parameters such as debulk time, moisture content and part size with porosity in laminates during the OOA process.   Roy [67] examined the relationship between the mechanisms driving void evolution. He parametrically investigated the effect of ply drops and caul sheets on the porosity content and final thickness profile of the composite parts through examining different configured composites. The lack of resin pressure as a result of pressure shielding due to relative rigidity of the caul sheet can lead to porosity in the composite material.  2.5 Summary  According to the brief literature review presented in this chapter, an accurate and efficient process model applicable to the whole process cycle of composite materials is not available. Process modelling of composite materials is typically performed within separate decoupled modules and results are shared between them in a tedious and inefficient process. In such methods, important phenomena and physics, such as pressure history and spatial variation of gelation in the resin, are ignored. The previously developed two-phase integrated model is not 25  sufficiently accurate, as it ignores some important material parameters in the transversely isotropic formulation of the integrated model (e.g. the fibres thermal expansion coefficient before gelation), and is not applicable for all ranges of constituents’ properties of the composite material (highly compressible phases). It is important to consider the gas phase separately in the formulation as it leads to prediction of porosity in the final composite part which is a very important quality issue in composite processing. Developing an accurate three-phase integrated flow-stress model is highly advantageous and provides the motivation for this thesis and its content presented in the subsequent chapters.  26  2.6 Figures   Figure 2.1 Schematic showing the sub-model structure used in COMPRO [8]  Autoclave Simulation ModuleThermochemical ModuleFlow ModuleStress ModuleOutputModuleState VariablesMaterial ModuleInput Module27    Figure 2.2 Spring-piston analogy for effective stress [68]    Figure 2.3 CV-FEM Discretization of the domain  NodeControl VolumeElementControl Surface28   Figure 2.4 Schematic showing an integrated model structures   State VariablesMaterial ModuleInput ModuleIntegrated Model29  Chapter 3: Isotropic Integrated Flow-Stress Model (IFS-Isotropic) In this chapter, a new methodology is presented to integrate the simulation of flow and stress development into a unified computational modelling framework for processing of isotropic composite materials. Isotropic model is applicable to uniformly distributed short fibre composites or the transverse plane of continuously reinforced unidirectional composites. The motivation for development of the isotropic formulation prior to the more general orthotropic formulation presented in Chapter 4, is to establish the groundwork and convey the basic concepts of the integrated flow-stress model without the complications that arise from material orthotropy. It should be noted that the integrated model developed previously by Haghshenas [17] provided  the foundation for constructing the current IFS model which is compatible with physics of the composite material behavior at different stages of processing and includes a broader range of constituents’ properties compared to the previous IFS model. First, the governing equations are developed for the general case of a two-phase composite material system that, as a consequence of curing, undergoes a transition from a fluid-like state into an elastic solid. The constitutive equations employed are such that they provide a continuous representation of the evolving material behaviour while maintaining consistency with the formulations that are typically used to represent the material at each of the two extremes. The two-phase formulation is capable of handling highly compressible phases, which is an important consideration when extending the model to a three-phase model that includes gas as a distinct phase. The 2D plane strain finite element implementation of the two-phase isotropic model within a 𝑢-𝑣-𝑃 framework is presented.  30  In subsequent sections, an extension of the developed two-phase integrated flow-stress process model for composites is presented to account for the presence of the gas phase within a three-phase composite system. The extended three-phase model incorporates the multi-physics phenomena involving flow of both the resin and gas phases through the fibre-bed as well as the previously considered thermochemical transformations during resin cure. Finally, the finite element implementation of the three-phase integrated flow-stress model is presented.   3.1 Formulation of Two-Phase Model (2IFS-Isotropic) The conservation equations of the 2IFS-Isotropic model is firstly presented in this section followed by the constructive equations that account for continuous evolution of the material properties as the polymeric resin solidifies during curing. Then, final form of the governing equations of the two-phase isotropic system is shown through combination of the constitutive and conservation equations. Implementation of governing equations in a 𝑢-𝑣-𝑃 FE formulation is presented in the subsequent subsection.  3.1.1 Conservation Equations This section presents the conservation equations (balance laws) for the poroelastic two-phase system. The two-phase system considers a compressible fluid phase (which can solidify during the process) flowing through a porous medium (fibre-bed) referred to as the solid-skeleton which is made of solid material (fibre). A schematic of the two-phase composite system and its constituents is shown in Figure 3.1. The conservation equations include the mass and momentum conservation of the matrix phase, as well as the momentum conservation (equilibrium) of the system.  31   3.1.1.1 Mass Conservation  Assuming small strains and ignoring spatial variation of densities and phase volume fractions, the mass conservation equations for a compressible solid and a compressible fluid phase can be written in indicial notation as: (1 − 𝜑𝐹)?̇?𝑖,𝑖 +(1 − 𝜑𝐹)𝜌𝑆?̇?𝑆 − ?̇?𝐹 = 0 (3-1) 𝑣𝐹𝑖,𝑖 + 𝜑𝐹?̇?𝑖,𝑖 +𝜑𝐹𝜌𝐹?̇?𝐹 + ?̇?𝐹 = 0 (3-2) respectively, where the subscripts 𝐹, and 𝑆 denote fluid (e.g. liquid or gas) and solid (e.g. solid grains in soil mechanics or fibres in composite materials) phases. Also, the dummy index 𝑖=1,…,3 denotes coordinate directions and as is customary in index notation, repeated index implies summation. Coma and overdot denote spatial and time derivatives, respectively, 𝑢 is the displacement of the system (associated with the solid structure), 𝜌 is density, 𝜑𝐹 is the volume fraction of the fluid phase, and 𝑣𝐹 is the volume-averaged relative velocity of the fluid phase (i.e. flow velocity).  The latter is defined in terms of the fluid velocity, ?̅?𝐹, and displacement of the system, 𝑢, as: 𝑣𝐹 = 𝜑𝐹(?̅?𝐹 − ?̇?) (3-3)  Summation of Eqs. (3-1) and (3-2) yields the following mass conservation equation for the two-phase system: ?̇?𝑖,𝑖 + 𝑣𝐹𝑖,𝑖 − ?̇? = 0 (3-4) where ?̇? represents the rate of change of the compressibility of the system given by 32  ?̇? = −𝜑𝐹𝜌𝐹?̇?𝐹 −(1 − 𝜑)𝜌𝑆?̇?𝑆 (3-5) When flow is inhibited, i.e. 𝑣𝐹𝑖,𝑖=0, this quantity represents the rate of volumetric strain of the system denoted by ?̇? = ?̇?𝑖,𝑖 in this work.  It is noted that the advection terms are ignored in the mass conservation equation, i.e.   1𝜌𝑆(𝜌𝑆(1 − 𝜑𝐹)),𝑖?̇?𝑖 ≅ 0 and 1𝜌𝐹(𝜌𝐹𝜑𝐹),𝑖(𝑣𝐹𝑖𝜑𝐹+ ?̇?𝑖) ≅ 0. This is because spatial variation of density and volume fraction are assumed to be small compared to other terms. These are assumptions made in previous formulations for both composite materials and soil mechanics [13,33] applications (among other examples of published work in the literature). Using this linearization assumption for the model to be solved numerically, time step sizes have to be chosen so that strain increments are also small in each time step.   3.1.1.2 Equilibrium of the Matrix In this work, Darcy’s law (see section 2.1.1) is used as equilibrium equation of the fluid phase: 𝑃𝐹,𝑖 + 𝜇𝐹 𝑆𝑖𝑗−1 𝑣𝐹𝑗 = 0 (3-6) where 𝑃𝐹 is the fluid pore pressure, 𝜇𝐹 is the fluid viscosity and 𝑆𝑖𝑗 is the permeability tensor of the solid-skeleton which is an average representation (or effective property) of the porous medium. In Figure 3.2, a macroscopic representation of fluid flow through the porous medium is presented representing both axial and transverse flow with respect to the fibre orientations. The equation governing the flow of the fluid phase is not affected when the fluidic resin undergoes cure.  However, the evolution of viscosity, 𝜇𝐹, during the cure process needs to be taken into account in the model. Darcy’s phenomenological equilibrium equation remains valid in the case 33  of a curing resin and therefore Eq. (3-6) for the unsolidified material applies identically to a solidifying material.  3.1.1.3 Equilibrium of the System Neglecting body forces and inertia effects, the equilibrium equation of the system, including both fluid and solid phases, can simply be written as  𝜎𝑖𝑗,𝑗 = 0 (3-7) where 𝜎 is the total stress in the system. As mentioned in Chapter 1, in soil mechanics for a two-phase system with compressible phases, the total stress of the system, 𝜎, is defined as: 𝜎𝑖𝑗 = 𝜎𝑆𝐾𝑖𝑗 − 𝑏 𝑃𝐹𝛿𝑖𝑗 (3-8) where the subscript 𝑆𝐾 denotes solid-skeleton (e.g. fibre-bed in this case) and 𝛿𝑖𝑗 is the Kronecker delta (or identity tensor). 𝜎𝑆𝐾 is identical to the so-called effective stress of the skeleton when Biot’s coefficient, 𝑏, is unity given by Eq. (2-3).  Assuming that the spatial variation of the volume fraction of the phases is negligible, we can ignore the spatial variation of the solid-skeleton properties, i.e. 𝐾𝑆𝐾,𝑖 ≈ 0 which leads to 𝑏,𝑖 ≈ 0. Combining Eqs. (3-7) and (3-8), the equilibrium equation of the system can be written as: 𝜎𝑆𝐾𝑖𝑗,𝑗 − 𝑏 𝑃𝐹,𝑖 = 0 (3-9)  3.1.2 Constitutive Equations The constitutive equations are constructed in such a way that at the intermediate stages of curing, they account for the varying composition of the constituents of the composite material (e.g. solid, 34  skeleton, and resin) as the resin properties become a mixture of the properties of a liquid phase and a solid phase with a ratio that varies according to the degree of solidification of the resin.  3.1.2.1 Bulk Response The rate of change of the density of the system, ?̇?, as appeared in the mass conservation equation (Eqs. (3-4) and (3-5)) is presented in more detail in this section for the two cases of (i) unsolidified matrix phase, and (ii) solidifying matrix phase.    Unsolidified Matrix Phase Based on the definition of bulk modulus and taking into account the effect of free strains, the following equation can be written for the fluid phase  1𝜌𝐹?̇?𝐹 =1𝐾𝐹?̇?𝐹 − 𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗 (3-10) where 𝐾𝐹 is the bulk modulus of the fluid and the superscript 𝑓 is used to denote the free component of strain (or strain rate) quantity to which it is attached.  The free strain of the curing fluid includes both the thermal and cure shrinkage effects as follows: 𝜀?̇?𝑓 = 𝛼𝐹𝑡ℎ?̇? + 𝛼𝐹𝑐𝑠?̇? (3-11) where 𝛼𝑡ℎ denotes the coefficient of thermal expansion, 𝛼𝑐𝑠 is coefficient of cure shrinkage, ?̇? is the rate of change of temperature, and ?̇? is the rate of change of the degree of cure of the resin. It is important to note that the matrix may cure without solidifying (i.e. it is not transforming from a liquid to a solid). Therefore, the free strains from cure shrinkage should be taken into account even for the case of unsolidified matrix.  35  Incorporating the effect of the free strains in the relation introduced by Lewis and Schrefler [69] for the solid phase, we can write 1𝜌𝑆?̇?𝑆 =1𝐾𝑆?̇?𝑆 −13(1 − 𝜑𝐹)𝐾𝑆?̇?𝑖𝑖′ − 𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗 (3-12) where 𝑃𝑆 is the pressure transferred to the solid phase and 𝜎′𝑖𝑗 is the effective stress tensor in the solid phase (Eq. (2-1)). Lewis and Schrefler [69] also presented the following constitutive equation for the rate of the first invariant (or trace) of the effective stress tensor ?̇?𝑖𝑖′ = 3𝐾𝑆𝐾 (?̇?𝑖,𝑖 +1𝐾𝑆?̇?𝑆) (3-13) According to Zienkiewicz and Shiomi [13], the rate of change of the skeleton volume (per unit volume of the system) due to free strain of the skeleton is given by: ?̇?𝑆𝐾𝑓 =13𝐾𝑆𝛿𝑖𝑗𝐷𝑆𝐾𝑖𝑗𝑘𝑙 𝜀?̇?𝐾𝑓𝑘𝑙 (3-14) where 𝐷 denote the elastic constitutive stiffness matrix of the material. The rate of free strain of the solid and solid-skeleton, 𝜀?̇?𝑓 and  𝜀?̇?𝐾𝑓, respectively, only includes the thermal effects which for the case of an unsolidified matrix phase is defined as: 𝜀Φ̇𝑓 = 𝛼Φ𝑡ℎ?̇? ,   Φ ≡ 𝑆 or 𝑆𝐾 (3-15)  Assuming that all the pressure in the fluid is transferred to the solid phase, 𝑃𝐹 = 𝑃𝑆, we can introduce Eqs. (3-10)-(3-14) into Eq.(3-5). Then, using the definition of Biot’s coefficient in Eq. (2-3) and including the effect of free-strains for the solid-skeleton (Eq. (3-14)), the total rate of change of the system density for the case of a liquid-like (unsolidified) fluid phase (that is representative of the behaviour of the resin in its pre-gelation regime, applicable to the early stages of the cure cycle) can be obtained as: 36  ?̇? = (1 − 𝑏)?̇?𝑖,𝑖 − (𝜑𝐹𝐾𝐹+𝑏 − 𝜑𝐹𝐾𝑆) ?̇?𝐹 + 𝜑𝐹𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗+ (1 − 𝜑𝐹)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗+13𝐾𝑆𝛿𝑖𝑗𝐷𝑆𝐾𝑖𝑗𝑘𝑙 𝜀?̇?𝐾𝑓𝑘𝑙 (3-16)  In Eq. (3-16), the term, (1 − 𝑏)?̇?𝑖,𝑖 , is storage due to volumetric expansion of the solid-skeleton while the term, (𝜑𝐹𝐾𝐹+𝑏−𝜑𝐹𝐾𝑆) ?̇?𝐹 , represents storage due to compressibility of the fluid and solid grains as a result of the pore pressure in the fluid phase.    Solidifying Matrix Phase In processing of composite materials, the matrix phase exhibits different behaviour at different stages of the process. For example, thermoset resins at elevated temperatures usually undergo three stages as a consequence of chemical reactions. These materials first show fluid-like viscous behaviour, followed by a gel- or rubber-like behaviour, and finally a solid- or glassy-like behaviour [7,15,50].   Eq. (3-16) applies to a composite material with a fluid phase that represents its unsolidified state. This equation is further extended here to the general case when the resin undergoes a continuous transformation into a fully solid material. To enable the integration within a unified IFS model and account for a continuously hardening (solidifying) fluid phase, a new scalar quantity, 𝜆, termed ‘solidification factor’ is introduced. This parameter is a measure of the solidification of the fluid phase (and in a complex interpretation, its shear load transfer capability between individual fibres within the composite material) and varies between zero and unity, 0 ≤ 𝜆 ≤ 1 37  (𝜆 = 0 corresponds to the unsolidified fluid resin and 𝜆 = 1 corresponds to the fully solid resin).  The solidification factor could be a function of any state variable used in process modelling, such as degree of cure (𝜒), temperature (𝑇), viscosity (𝜇), or even fibre volume fraction and the fibre-bed architecture, depending on the type of material, type of process, and so forth.  As shown schematically in Figure 3.3 for 𝜆 = 0, the fluid phase (resin or gas, as a special case of fluid), solid phase (fibres), and solid-skeleton (fibre-bed) are incorporated in the model according to the Biot’s poroelastic model [11-13]. At the other extreme when 𝜆 = 1, we have a solid composite material. In the latter case, there are no distinct solid and solid-skeleton phases as everything is considered homogenized (smeared) into a continuum, solid composite material consisting of fibres in a fully solid matrix. In the intermediate regime (when 0 < 𝜆 < 1), the fluid phase is partially solidified and has a combination of fluid- and solid-like behaviour. While the material undergoes curing, the fluid phase/solid grains gradually become more and more part of the composite material (through contributing to the solid-skeleton) and less and less as a pure fluid/solid. As a simple mathematical descriptor of this behaviour using the scalar parameter, 𝜆, the fluid phase is homogenized partially with the solid-skeleton. Therefore, all the properties of the solid-skeleton evolves from the fibre-bed properties at one extreme to the fully homogenized composite material at the other extreme as a function of 𝜆. For solid grain properties, this evolution starts from the fibre properties and terminates with the homogenized composite properties. In the general functional form, we can write the elastic constitutive stiffness matrix, the coefficient of thermal expansion (𝛼𝑡ℎ), and the coefficient of cure shrinkage (𝛼𝑐𝑠) as: 𝐃𝑆 = 𝑓(𝐃𝑓 , 𝐃𝑐, 𝜆)  ,  𝐃𝑆𝐾 = 𝑓′(𝐃𝑓𝑏, 𝐃𝑐, 𝜆)  (3-17) 38  𝛼𝑆𝐾𝑡ℎ = 𝑔(𝛼𝑓𝑏𝑡ℎ, 𝛼𝑐𝑡ℎ, 𝜆)  ,  𝛼𝑆𝐾𝑐𝑠 = 𝑔′(𝛼𝑐𝑐𝑠, 𝜆) where subscripts 𝑓, 𝑓𝑏, and 𝑐 denote fibre, fibre-bed, and composite quantities, respectively, and bold symbols represent matrix (tensor) quantities. It should be noted that these material properties depend on other variables as well. For example, in composite materials, the fibre-bed or skeleton stiffness properties, 𝐃𝑆𝐾, is dependent nonlinearly on the instantaneous volume fractions of the fibre even in the absence of solidification in the fluid. In Eq. (3-17), however, only their dependency on solidification factor is highlighted. Composite properties in Eq. (3-17), i.e. 𝐃𝑐, 𝛼𝑐𝑡ℎ, and 𝛼𝑐𝑐𝑠, at any instant during the simulation (i.e. at any stage of the processing as shown in Figure 3.3) can be derived from different sources e.g. any appropriate numerical or analytical micromechanics models. The Biot coefficient, 𝑏, evolves as both 𝐾𝑆 and 𝐾𝑆𝐾 undergo changes during the course of curing, according to Eq. (2-3).  In the unsolidified state of the fluid (i.e. when 𝜆 = 0), 𝐾𝑆 = 𝐾𝑓 and 𝐾𝑆𝐾 = 𝐾𝑓𝑏. In this case, 𝑏 is usually very close to unity as 𝐾𝑓𝑏 ≪𝐾𝑓 for composite materials. On the other hand, for a fully cured composite material (i.e. when 𝜆 = 1), 𝐾𝑆 = 𝐾𝑆𝐾 = 𝐾𝑐, and consequently 𝑏 = 0. Therefore, the variable 𝑏 varies continuously from a value very close to unity for an uncured material to zero for a fully cured material. The gradual transformation of the fluid phase from an unsolidified to a solidified state is introduced in the expression for the free volumetric strain rate by incorporating the solidification factor, 𝜆 as follows: ?̇?𝐹𝑓 = (1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗 (3-18) ?̇?𝑆𝑓 = (1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗 (3-19)  39  Making use of Eqs. (3-18) and (3-19) for the evolving free volumetric strain rate due to solidification, results in the following modified form of Eq. (3-16) for the constitutive bulk response: ?̇? = (1 − 𝑏)?̇?𝑖,𝑖 − (𝜑𝐹𝐾𝐹+𝑏 − 𝜑𝐹𝐾𝑆) ?̇?𝐹 + 𝜑𝐹(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗+ (1 − 𝜑𝐹)(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗+13𝐾𝑆𝛿𝑖𝑗𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝑓𝑘𝑙 (3-20)  Although the quantitative choice of solidification factor 𝜆 and its evolution are beyond the scope of this work, it should be remarked that its evolution must be synchronized with other fluid material properties evolution. For example, a fluid phase whose bulk and shear moduli are in the range of liquid materials may not have a solidification factor 𝜆 = 1 which is assigned to solid materials. Poor choice of the solidification factor might cause ill-conditioning of the matrices in the simulations that generally arise from using non-physical material properties in the model.  3.1.2.2 Effective Stress (Non-Bulk) Response It is customary to express the constitutive equations in differential or rate form when they are introduced in the system of equations used to solve nonlinear problems (e.g. through nonlinear finite element analysis procedures). In the current problem, nonlinearities can arise from at least two sources, one being the inherently nonlinear constitutive properties of the fibre-bed under compaction [15,70] and the other is the geometrical nonlinearities (large strains) that arise from possible high compressibility of the phases. Moreover, the constitutive response of the composite as it undergoes curing is in general viscoelastic in nature.  In their simplified form they can be assumed to be hypoelastic (e.g. cure hardening instantaneously linear elastic, as has successfully 40  been used in COMPRO). Therefore, in deriving the constitutive response that relates the effective stress of the system to the corresponding strains and hence displacements of the system, we will employ the differential form (rather than total form) of all stress quantities for both unsolidified and solidifying matrix.   Unsolidified Matrix Phase Writing the rate form of the total stress given by Eq. (2-1), we have: ?̇?𝑖𝑗 = ?̇?𝑆𝐾𝑖𝑗 − 𝑏 ?̇?𝐹𝛿𝑖𝑗 − ?̇? 𝑃𝐹𝛿𝑖𝑗 (3-21) where ?̇? is the rate of change of 𝑏 and ?̇?𝑆𝐾 is related to the rate of the mechanical strain of the solid-skeleton, 𝜀?̇?𝐾𝜎 , which is the same as the mechanical strain of the system, through its elastic constitutive stiffness tensor (which coincides with the stiffness properties of the fibre-bed when the matrix phase is unsolidified), as follows: ?̇?𝑆𝐾𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙 (3-22)  In unsolidified composite materials as the bulk modulus of the solid-skeleton phase, 𝐾𝑆𝐾, is negligible compared to the bulk modulus of the solid grains (fibres), 𝐾𝑆, the value of the Biot’s coefficient  remains very close to unity. It can be shown that ?̇? in an unsolidified material, which reflects the variation of 𝑏 solely as a result of the nonlinearity in solid or skeleton properties (e.g. due to compaction), is almost zero for composite materials. This is because in polymeric fibrous composite materials, as mentioned before, the magnitude of solid-skeleton (fibre-bed) bulk modulus is negligible compared to the solid (fibres) bulk modulus which leads to Biot’s coefficient of almost unity. Although, the solid-skeleton properties are affected by fibres volume 41  fraction change during the compaction, its bulk modulus typically remains negligible compared to the fibres even at high volume fractions. This results in a minor variation in 𝑏, i.e. ?̇? ≅ 0. Consequently, we can write: ?̇?𝑆𝐾𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙≅ 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙+ ?̇? 𝑃𝐹𝛿𝑖𝑗 (3-23) Combining Eqs. (3-21) and (3-23) yields  the total stress rate as: ?̇?𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙− 𝑏 ?̇?𝐹𝛿𝑖𝑗 (3-24) which is identical to the stress rate equation in classical porous media formulation used for soils where nonlinearity arising from the evolution of material properties of the porous structure is ignored (?̇? ≅ 0). In solid mechanics studies, the last term in Eq. (3-21), −?̇? 𝑃𝐹𝛿𝑖𝑗, is usually ignored and Eq. (3-22) is used for solid-skeleton stress which results in the same total stress constitutive equation as Eq. (3-24) used in this work.   Solidifying Matrix Phase In a solidifying material, both the solid and solid-skeleton properties evolve significantly from their initial values to the fully solid composite material (based on the methodology adopted here and explained in Section 3.1.2.1.2). This leads to a significant change in Biot’s coefficient value, 𝑏, from almost unity to zero.  Therefore, as resin solidifies, the magnitude of ?̇? can be significant as a result of variation of material properties due to solidification as described in Eq. (3-17) (i.e. changes in both 𝐾𝑆𝐾 and 𝐾𝑆). Conceptually, as the resin material undergoes solidification and hence contributes to load (stress) carrying capacity of the system, a portion of the pressure (?̇? 𝑃𝐹) in the unsolidified fluid should be incrementally added to the stress carried by the solid-skeleton. This additional stress can be regarded as a residual stress induced in the material due to 42  transformation of the pressure bearing fluid into a homogenized solid composite material. In this regard, ?̇?𝑆𝐾 is related to the skeleton strain through its elastic constitutive stiffness tensor, as follows: ?̇?𝑆𝐾𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙+ ?̇? 𝑃𝐹𝛿𝑖𝑗 (3-25)  In the sequential method, the history of the pressure is ignored at the start of the stress development module and this can affect the final stress state as demonstrated conceptually through the example shown in Figure 3.4. This figure shows a cube of neat resin in its fluid state (the presence of a solid phase is ignored for sake of simplicity). If this resin undergoes pressure (Figure 3.4 (a)) or experiences a temperature rise and then cures (Figure 3.4 (b)), the pressure induced in the fluid resin has to be transformed to the compressive stresses carried by the solid (cured) resin. In contrast to the sequential approach, the proposed methodology allows for the transformation of the pressure carried by the liquid resin to the solid composite material in a seamless fashion. Note that the form of the above equation is different than what is used in classical Biot’s model (for an unsolidified fluid flow through a porous medium) as shown in Eq. (3-22). Combining Eqs. (3-21) and (3-25) yields:  ?̇?𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙− 𝑏 ?̇?𝐹𝛿𝑖𝑗 (3-26)  In summary, the portion of ?̇? as a result of changes in 𝐾𝑆𝐾 due to compaction of the solid-skeleton is negligible and can be ignored in the model. However, the variation of 𝑏 due to fluid solidification is quite significant and has to be taken into account as some form of residual stress during the process that augments the stress in the solid-skeleton, 𝜎𝑆𝐾 (Eq. (3-25)).  43   3.1.3 Summary of the Governing Equations (2IFS-Isotropic) Introducing the rate of change of the density of the system, ?̇? as given by Eq. (3-20), into Eq. (3-4) yields the mass conservation equation for the IFS model. Within the framework of the IFS model, the balance laws or conservation equations consisting of the mass conservation of fluid, equilibrium of fluid, and equilibrium of the system written in terms of the basic variables, 𝑢, 𝑣𝐹 and 𝑃𝐹 can be summarized as follows:  Mass conservation: 𝑏?̇?𝑖,𝑖 + 𝑣𝐹𝑖,𝑖 + (𝜑𝐹𝐾𝐹+𝑏 − 𝜑𝐹𝐾𝑆) ?̇?𝐹 − 𝜑𝐹(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗− (1 − 𝜑𝐹)(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗−13𝐾𝑆𝛿𝑖𝑗𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝑓𝑘𝑙= 0 (3-27)  Darcy’s flow for the fluid: 𝑃𝐹,𝑖 + 𝜇𝐹𝑆𝐹𝑖𝑗−1𝑣𝐹𝑗 = 0 (3-28)  Equilibrium equation: 𝜎𝑆𝐾𝑖𝑗,𝑗 − 𝑏 𝑃𝐹,𝑖 = 0 (3-29) where, using the strain-displacement relationship, the rate of stress in the solid-skeleton can be expressed as: ?̇?𝑆𝐾𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙 (12(?̇?𝑘,𝑙 + ?̇?𝑙,𝑘) − 𝜀?̇?𝐾𝑓𝑘𝑙) + ?̇? 𝑃𝐹𝛿𝑖𝑗 (3-30) Also, 𝑏 = 1 − 𝐾𝑆𝐾/𝐾𝑆, 𝐃𝑆 = 𝑓(𝐃𝑓 , 𝐃𝑐, 𝜆) , 𝐃𝑆𝐾 = 𝑓(𝐃𝑓𝑏 , 𝐃𝑐, 𝜆), 𝛼𝑆𝐾𝑡ℎ = 𝑓(𝛼𝑓𝑏𝑡ℎ, 𝛼𝑐𝑡ℎ, 𝜆), and 𝛼𝑆𝐾𝑐𝑠 = 𝑓(𝛼𝑐𝑐𝑠, 𝜆). The elastic constitutive stiffness matrix, 𝐃, for a 3D isotropic solid material is given by: 44  𝐃 =[     𝐾 + 4𝐺/3 𝐾 − 2𝐺/3 𝐾 − 2𝐺/3 𝐾 + 4𝐺/3 𝐾 − 2𝐺/3  𝐾 + 4𝐺/30 0 00 0 00 0 0    𝑆𝑦𝑚𝑚.     𝐺 0 0 𝐺 0  𝐺]      (3-31) Eqs. (3-27)-(3-29), govern the response of the poroelastic system during the curing process and can be used to determine the flow of fluid and stress development in the composite material at any instant of time. In the following sections we will check for the consistency of the governing equations with those that apply under the two extreme conditions of uncured (unsolidified) and fully cured (solidified) materials.   3.1.3.1 Special Case: Unsolidified Material (𝝀 = 𝟎) When the fluid phase (e.g. resin) is unsolidified (corresponding to the far left schematic representation shown in Figure 3.3), 𝜆 = 0, 𝐃𝑆 = 𝐃𝑓, 𝐃𝑆𝐾 = 𝐃𝑓𝑏, 𝛼𝑆𝐾𝑡ℎ = 𝛼𝑓𝑏𝑡ℎ, and 𝛼𝑆𝐾𝑐𝑠 = 0. Therefore, the model will be identical to that of the classical Biot’s model for flow though porous media used widely in soil mechanics applications in which the following governing equations apply: 𝑏?̇?𝑖,𝑖 + 𝑣𝐹𝑖,𝑖 + (𝜑𝐹𝐾𝐹+𝑏 − 𝜑𝐹𝐾𝑓) ?̇?𝐹 − 𝜑𝐹𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗− (1 − 𝜑𝐹)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗−13𝐾𝑓𝛿𝑖𝑗𝐷𝑓𝑏𝑖𝑗𝑘𝑙𝜀?̇?𝑏𝑓𝑘𝑙= 0 𝑃𝐹,𝑖 + 𝜇𝐹 𝑆𝑖𝑗−1 𝑣𝐹𝑗 = 0 𝜎𝑓𝑏𝑖𝑗,𝑗 − 𝑏 𝑃𝐹,𝑖 = 0 (3-32) where 45  ?̇?𝑓𝑏𝑖𝑗 = 𝐷𝑓𝑏𝑖𝑗𝑘𝑙 (12(?̇?𝑘,𝑙 + ?̇?𝑙,𝑘) − 𝜀?̇?𝑏𝑓𝑘𝑙) + ?̇? 𝑃𝐹𝛿𝑖𝑗 (3-33) where ?̇? ≅ 0. In this case, the elastic constitutive stiffness matrices of the skeleton and the solid phase under 2D plane strain conditions, for isotropic materials, are given by 𝐃𝑆𝐾|𝜆=0 = 𝐃𝑓𝑏 = [𝐾𝑓𝑏 + 4𝐺𝑓𝑏 3⁄ 𝐾𝑓𝑏 − 2𝐺𝑓𝑏 3⁄ 0 𝐾𝑓𝑏 + 4𝐺𝑓𝑏 3⁄ 0𝑆𝑦𝑚𝑚.  𝐺𝑓𝑏] (3-34) 𝐃𝑆|𝜆=0 = 𝐃𝑓 = [𝐾𝑓 + 4𝐺𝑓 3⁄ 𝐾𝑓 − 2𝐺𝑓 3⁄ 0 𝐾𝑓 + 4𝐺𝑓 3⁄ 0𝑆𝑦𝑚𝑚.  𝐺𝑓] (3-35)  It should be noted that in the IFS model developed previously [17], the mass conservation equation was considered for both the fluid and solid phases as ?̇?𝑖,𝑖 + 𝑣𝐹𝑖,𝑖 − (?̇?σ + ?̇?f) = 0 (3-36) where the superscript, 𝜎,  denotes the mechanical component of the quantity to which it is attached. The assumption of lumping the rate of change of the density of the two phases into a single system (or composite) quantity, ?̇?𝑐 = ?̇?𝜎 + ?̇?𝑓, given in the above equation, enabled bridging of the fluid and solid regimes of the process simulation with the emphasis being placed on attaining the correct (consistent) micromechanics of the composite in its cured (solid) state. However, the undesirable consequence of this assumption is that it imposes a constraint on the load sharing between the two phases during the early stages of the process (when the resin is in a fluid-like state) leading to inaccurate response predictions. In fact, for cases where the fluid phase is highly compressible, the numerical solution using the previous IFS model breaks down. An eigenvalue analysis (on the 𝐈 −12(𝐾𝑐+𝐺𝑐)𝐃𝑠−𝑃 matrix in Eq. (4.135) of reference [17]) 46  indicates that the previous IFS model only applies to cases in which 𝐾𝑐 > 𝐾𝑓𝑏 + 𝐺𝑓𝑏/3. For low values of 𝐾𝑐 that conceivably arise when the fluid phase is highly compressible (e.g. gas or partially saturated resins), the above inequality may be violated leading to an ill-conditioned stiffness matrix and hence unstable numerical solution. In the newly formulated IFS model presented in this work, there is no requirement to invoke the solid micromechanics model when the fluid phase is unsolidified (see Eq.(3-27)).   3.1.3.2 Special Case: Fully Solidified Material (𝝀 = 𝟏) When the fluid phase (e.g. resin) is in a fully cured and solidified state (corresponding to the far right schematic representation shown in Figure 3.3), there is no distinction between the solid and skeleton components of the composite. Therefore, under this condition, 𝑏 = 0, 𝐃𝑆 = 𝐃𝑆𝐾 = 𝐃𝐜, 𝛼𝑆𝐾𝑡ℎ = 𝛼𝑐𝑡ℎ, 𝛼𝑆𝐾𝑐𝑠 = 𝛼𝑐𝑐𝑠, and Eqs. (3-27) and (3-28) become decoupled from Eq. (3-29). Essentially, the first two equations reduce to a coupled system of equations in terms of the fluid velocity and pressure (𝑣𝐹 , 𝑃𝐹) both of which have no physical meaning in this limiting case of  a fully solidified material. What matters at this stage are the displacements, 𝑢𝑖, and the stresses in the homogenized composite material which can be obtained solely from the Eq. (3-29) written in its following special form:   𝜎𝑐𝑖𝑗,𝑗 = 0 (3-37) where ?̇?𝑐𝑖𝑗 = 𝐷𝑐𝑖𝑗𝑘𝑙 (12(?̇?𝑘,𝑙 + ?̇?𝑙,𝑘) − 𝜀?̇?𝑓𝑘𝑙) (3-38) The 2D plane strain, elastic constitutive stiffness matrix of the composite, 𝐃𝑐 is given by  47  𝐃𝑐 = 𝐃𝑆𝐾|𝜆=1 = 𝐃𝑆|𝜆=1 = [𝐾𝑐 + 4𝐺𝑐 3⁄ 𝐾𝑐 − 2𝐺𝑐 3⁄ 0 𝐾𝑐 + 4𝐺𝑐 3⁄ 0𝑆𝑦𝑚𝑚.  𝐺𝑐] (3-39) where 𝐾𝑐 and 𝐺𝑐 are obtained from the adopted micromechanics model for the solid composite.   It should be remarked in passing that while the coupled system of equations (Eqs. (3-27)-(3-29)) are solved for the basic variables 𝑢, 𝑣𝐹 and 𝑃𝐹 for all values of 𝜆 between 0 and 1, in order to avoid computing redundant quantities and thereby increase the efficiency of the numerical model at the instant when full solidification is attained, only Eq. (3-37), which represents the equilibrium equation of the homogenized solid composite material is solved for the unknown displacement variable, 𝑢.  3.1.4 Finite Element Implementation (2IFS-Isotropic) The finite element implementation of the IFS methodology developed in the preceding section can be achieved by converting the governing differential equations of the system (Eqs. (3-27)-(3-29)) into their weak form using the Galerkin method. Within a 𝑢-𝑣-𝑃 FE formulation adopted here, the system displacement vector, 𝑢𝑖, relative velocity vector of the fluid phase, 𝑣𝐹𝑖, and fluid pore pressure, 𝑃𝐹, are discretized spatially.  Applying the backward Euler scheme for temporal integration combined with the Newton’s nonlinear solution method, results in an incremental form of a system of algebraic equations that are solved iteratively. The iterative matrix equation for the 𝑘𝑡ℎ iteration within the current time station 𝑡𝑛 (or the 𝑛𝑡ℎ time-step) of the solution can be written as 48  [Δ𝑡𝐊𝒖𝒖 𝟎 Δ𝑡𝐊𝒖𝑷𝐹𝟎 Δ𝑡𝐊𝒗𝐹𝒗𝐹 Δ𝑡𝐊𝒗𝑷𝐹𝐂𝑷𝐹𝒖 Δ𝑡𝐊𝑷𝐹𝒗𝐹 𝐂𝑷𝐹𝑷𝐹]𝑘−1{∆?̅?𝑘?̅?𝐹𝑘∆?̅?𝐹𝑘}={    Δ𝑡[𝐟𝒖 + 𝐟𝒖𝑓𝑆𝐾𝑘−1−𝐟𝑢𝑖𝑛𝑡𝑘−1]∆𝑡[𝐟𝒗𝐹 − 𝐊𝒗𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1]Δ𝑡𝐟𝑷𝐹 + 𝐂𝑷𝐹𝒖𝑘−1(?̅?0 − ?̅?𝑘−1) + 𝐂𝑷𝐹𝑷𝐹𝑘−1(?̅?𝐹0− ?̅?𝐹𝑘−1)}     (3-40) where Δ𝑡 = 𝑡𝑛 − 𝑡𝑛−1 is the time step, and ∆( )𝑘 = ( )𝑘 − ( )𝑘−1 where superscript zero indicates initial value at the start of time step (converged value at end of previous time step).  A low order bi-linear isoparametric element with 4 corner nodes for the system displacements and relative velocities of the fluid phase, ?̅? and ?̅?𝐹, respectively, as degrees of freedom and only one internal central node assigned to the fluid pressure ?̅?𝐹 (the so-called 4-1 element shown in Figure 3.5) is used to carry out the FE solution where a 2 × 2 Gauss quadrature is used for evaluating the integrals. All matrix and vector quantities in Eq. (3-40) are defined in Table 3.1. Details of derivation of Eq. (3-40) is presented in Appendix A.1.  A schematic of the algorithm adopted for implementation of the two-phase IFS model in the MATLAB® code is shown in Figure 3.6. In FE implementation, the free strain of the skeleton and fluid phase, 𝛆𝑆𝐾𝑓 and 𝛆𝐹𝑓, are additively decomposed into the thermal and cure shrinkage components as follows: Δ𝛆𝑓𝑘= Δ𝛆𝑡ℎ𝑘+ Δ𝛆𝑐𝑠𝑘 (3-41) where Δ𝛆𝑡ℎ𝑘= 𝛼𝑡ℎ Δ𝑇𝑘𝛅 (3-42) Δ𝛆𝑡ℎ𝑘= 𝛼𝑐𝑠 Δ𝜒𝑘𝛅 (3-43) 49  where 𝛅𝑇 = {1 1 0}.  As mentioned before, both coefficients of thermal expansion and cure shrinkage of the skeleton, 𝛼𝑆𝐾𝑡ℎ  and 𝛼𝑆𝐾𝑐𝑠 , are functions of the solidification factor, 𝜆. For the free strain of the solid phase, 𝛆𝑆𝑓, only the thermal effect is taken into account as the solid phase does not undergo solidification (solid phase does not cure).  Once Eq. (3-40) is solved for the new values of ?̅?, ?̅?𝐹, and ?̅?𝐹 at the 𝑛𝑡ℎ  time step and 𝑘𝑡ℎ iteration, the stress and the volume fraction of the fluid. Using Eq. (3-26) for stress computation, we can write:  𝛔𝑘 = 𝛔𝑘−1 + 𝐃𝑆𝐾𝑘−1𝐁𝑘Δ?̅?𝑘 − 𝐃𝑆𝐾𝑘−1Δ𝛆𝑆𝐾𝑓 𝑘−𝑏𝑘−1𝛅Δ?̅?𝐹𝑘 (3-44) where 𝐁 is the matrix of spatial derivatives of the displacement (or velocity) shape functions.   It is important to note that in the current IFS model which is intended to handle process simulation of composite materials with a highly compressible fluid phase, the deformations can be significant enough to warrant consideration of finite displacements in the formulations. Therefore, in order to capture nonlinearities due to large deformations, the geometry of the mesh (nodal coordinates) and all dependent finite element matrices and vectors are updated after each iteration. Provided that small enough time steps are used in the computations, this approach can be shown to be adequate for achieving numerical accuracy (convergence) of results.  3.1.4.1 Volume Fraction Calculation The volume of the fluid and solid phases at different locations within the domain, vary through time due to such factors as flow of the fluid phase, compressibility of each phase, and free 50  strains. In order to update the volume fraction of the fluid, Haghshenas [17] only considered flow of the fluid phase using the following equation 𝜑𝐹 =𝜑0𝐹 + ∆𝑒𝐹𝑣1 + ∆𝑒𝐹𝑣  (3-45) where 𝜑0𝐹 is the initial fluid volume fraction at the beginning of a time step, and ∆𝑒𝐹𝑣, representing the incremental volumetric strain of the system solely due to fluid flow, is given by: ∆𝑒𝐹𝑣 = − ∫ 𝑣𝐹𝑖,𝑖𝑑𝑡𝑡+Δ𝑡𝑡 (3-46)  Equation (3-45) neglects volume changes due to compressibility of the fluid and solid phases which is a plausible assumption when dealing with nearly incompressible materials such as fully saturated resins. However, when the fluid phase is highly compressible, e.g. gas or partially saturated resins with low degrees of saturation and therefore a prominent presence of a gas phase, compressibility should be taken into account in calculating and updating the volume fraction.  For this purpose, we may write the following equation for the current volume fraction of the fluid phase in terms of its initial volume fraction: 𝜑𝐹 =𝑉0𝐹 + ∆𝑉𝐹𝑉0 + ∆𝑉 (3-47) where 𝑉0𝐹 and 𝑉0 are the initial volumes of the fluid phase and the system, respectively, at the start of the time step, and ∆𝑉𝐹 and ∆𝑉 are the corresponding incremental changes. Eq. (3-47) can be re-written in the following form:  𝜑𝐹 =𝜑0𝐹 + ∆𝑉𝐹/𝑉01 + ∆𝑒 (3-48) 51  where 𝜑𝐹0 = 𝑉𝐹0/𝑉0 and  ∆𝑒 = ∆𝑉/𝑉0 represents the incremental volumetric strain of the system which can be obtained from the deformation of the system, 𝑢𝑖. The change in the volume of the fluid phase, ∆𝑉𝐹, can be expressed as:  ∆𝑉𝐹 = 𝑉0∆𝑒𝐹𝑣 + 𝑉𝐹0(∆𝑒𝐹𝜎 + ∆𝑒𝐹𝑓) (3-49) where ∆𝑒𝐹𝜎 and ∆𝑒𝐹𝑓 are the incremental volumetric mechanical (due to induced pressure) and free strains in the fluid phase, respectively. Using Eq. (3-49) in Eq. (3-48), the following equation is obtained for the calculation of the updated volume fraction at the end of the time step: 𝜑𝐹 =𝜑0𝐹(1 + ∆𝑒𝐹𝜎 + ∆𝑒𝐹𝑓) + ∆𝑒𝐹𝑣1 + ∆𝑒 (3-50) The equation for updating the volume fraction of the fluid as a result of the spatial discretization of Eq. (3-50) then becomes: 𝜑𝐹𝑘 =𝜑𝐹0 (exp (− ?̅?𝐹𝑘− ?̅?𝐹𝟎𝐾𝐹) + 𝛅𝑇Δ𝛆𝐹𝑓𝑘) − Δ𝑡𝛅𝑇𝐁𝑘?̅?𝐹𝑘1 + 𝛅𝑇𝐁𝑘(?̅?𝑘 − ?̅?0) (3-51)  In solid micromechanics models, changes in the constituents’ volumes are not considered in the calculation of homogenized composite material properties. Consequently, even though the volume of each constituent is changing throughout the process, Eq. (3-51) is used to update the volume fractions only up to the point when full solidification occurs (i.e. it is valid for λ < 1 ).  Therefore, the volume fraction of phases are kept constant in the model once 𝜆 reaches unity.   52  It should be noted that in the derivation of the volume fraction update equation, Eq. (3-50), no assumption is made on the magnitude of strains. In other words, this equation is valid for general and unrestricted magnitude of strains while the mass conservation equation is applicable to small strains (small during each time step). If the mass conservation equation of the fluid phase, Eq. (3-2), is used for updating the volume fractions, we will end up with the following equation:  𝜑𝐹𝑆 = 𝜑0𝐹(1 + ∆𝑒𝐹𝜎 + ∆𝑒𝐹𝑓 − ∆𝑒) + ∆𝑒𝐹𝑣 (3-52) where superscript 𝑆 stands for small. The volume fraction in Eq. (3-50) can be written in terms of 𝜑𝐹𝑆 as: 𝜑𝐹𝐿 =𝜑𝐹𝑆 + 𝜑0𝐹∆𝑒1 + ∆𝑒 (3-53) where superscript 𝐿 stands for large. For small incremental volumetric strains, i.e. ∆𝑒~0, we have 𝜑𝐹𝐿 = 𝜑𝐹𝑆. Although, both equations yield the same value for the updated volume fraction when the time steps are small (small strain during each time step), the examples in Figure 3.7 show the inaccuracy of Eq. (3-52) (and accuracy of the Eq. (3-53)) when we have very large strains. In all the three cases shown, the composite system experiencing a total volumetric strain of ∆𝑒 = 0.5. In Figure 3.7(a), the fluid phase expands due to a total mechanical and free strain of unity (∆𝑒𝐹𝜎 + ∆𝑒𝐹𝑓 = 1). The system volumetric strain in Figure 3.7(b) come from the solid phase expansion while the fluid volume remains constant. In Figure 3.7(c), flow of fluid phase into the system (∆𝑒𝐹𝑣 = 0.5) is the source of volumetric strain of the system where no mechanical and free strain is induced in the fluid phase.  53  3.2 Formulation of the Three-Phase Model (3IFS-Isotropic) The 2IFS-Isotropic model presented in the previous sections is extended in this part of Chapter 3 to include a second fluid phase as well. A schematic of the three-phase composite system is shown in Figure 3.8. In a similar format to the two-phase model, the conservation and constitutive equations are presented and then followed by a summary of governing equations and FE implementation of the 3IFS-Isotropic model.  3.2.1 Conservation Equations  The conversation equations of the two-phase system are extended here for a composite system including a solid phase and two fluid phases. In addition to the three conservation equations presented before for the two-phase system, an extra pressure equilibrium equation is also introduced to complete the number of governing equations required to describe the response of the three-phase system. As the two-phase system, these conservation equations apply at any stage of the processing whether the fluid phases are unsolidified or fully solidified to form a solid three-phase composite material.  3.2.1.1 Mass Conservation Similar to the two-phase system, we assume small strains and ignore spatial variation of densities and phase volume fractions. Therefore, the mass conservation equation for a system that comprises of a compressible solid phase and two compressible fluid phases, 𝐹 and ?̂?, can be written as: ?̇?𝑖,𝑖 + 𝑣𝐹𝑖,𝑖 + 𝑣?̂?𝑖,𝑖 − ?̇? = 0 (3-54) 54  The volume-averaged velocity for each fluid phase relative to the solid-skeleton (i.e. flow velocities 𝑣𝐹 and 𝑣?̂?) is defined in terms of fluid velocities, ?̃?𝐹 and ?̃??̂?, fluids volume fractions, 𝜑𝐹 and 𝜑?̂?, and displacement of the system and can be written as follows: 𝑣Φ = 𝜑Φ(?̃?Φ − ?̇?),   Φ ≡ 𝐹 or ?̂? (3-55) In Eq. (3-54), ?̇? represents the rate of change of the compressibility of the three-phase system given by: ?̇? = −𝜑𝐹𝜌𝐹?̇?𝐹 −𝜑?̂?𝜌?̂??̇??̂? −(1 − 𝜑)𝜌𝑆?̇?𝑆 (3-56) where 𝜌 denotes density and 𝜑 = 𝜑𝐹 + 𝜑?̂?.   3.2.1.2 Equilibrium of the Fluid Phases Using the macroscopic representation of the fluids flow as explained in Section 2.1.3, in a decoupled flow model which ignores momentum exchange across the interface between two fluids, the equilibrium equation for each fluid phase can be written as: 𝑃Φ,𝑖 + 𝜇Φ𝑆Φ𝑖𝑗−1𝑣Φ𝑗 = 0,   Φ ≡ 𝐹 or ?̂? (3-57) A macroscopic representation of two fluids flowing through the porous medium is presented in Figure 3.9 where the flow of the fluids are decoupled as the drag force (momentum exchange) between them is ignored.  As before, the equations governing the flow of the fluid phases are not affected when the fluidic resin undergoes curing. However, the model needs to account for the evolution of viscosity of each fluid phase during the cure process. Figure 3.10 shows a representation of the three-phase model for the flow of the fluids where the viscosity of fluid 𝐹 increases. It should be noted that 55  in this figure only the viscosity change of one of the two fluid phases and its effect on flow is highlighted.  3.2.1.3 Pressure Equilibrium In the three-phase system, the interfacial equilibrium condition between the two fluids can be written as: 𝑃?̂? − 𝑃𝐹 = 𝑃𝑐 (3-58) where 𝑃𝑐 is the macroscopic capillary pressure defined as the difference in pressure across the interface between the two fluids. Capillary pressure in composite materials is dependent on material properties such as degree of cure, volume fraction of phases, surface tension of fluids, permeability of porous media, etc. [36,37].  3.2.1.4 Equilibrium of the Three-Phase System In a macroscopic representation, the equilibrium equation of the system regardless of of number of phases involved, can be written in terms of the total stress, 𝜎, as: 𝜎𝑖𝑗,𝑗 = 0 (3-59) in which body forces and inertia effects are neglected.  Using Bishop’s formulation [30] for unsaturated flow through porous media as presented in Section 2.1.3, the total stress of the system can be written as: 𝜎𝑖𝑗 = 𝜎𝑆𝐾𝑖𝑗 − 𝑏𝑃𝑚𝛿𝑖𝑗 (3-60) where 𝑃𝑚 is the effective pore pressure of the mechanically equivalent single phase matrix, with the subscript 𝑚 representing the multiphase matrix. This effective pressure is defined as: 56  𝑃𝑚 = 𝑃?̂? − ℬ(𝑃?̂? − 𝑃𝐹) (3-61) where 0 ≤ ℬ ≤ 1 is the Bishop’s parameter.  3.2.2 Constitutive Equations The constitutive equations are constructed in such a way that at the intermediate stages of curing, they account for the varying composition of the constituents of the composite material (e.g. solid, resin, and gas) as the resin properties become a mixture of the properties of a liquid phase and a solid phase with a ratio that varies according to the degree of solidification of the resin.   3.2.2.1 Rate of Compressibility In a three-phase system, the same equation as that applied to the density of the single fluid phase in a two-phase system (Eq. (3-10)) is applied separately to each fluid phase. Also, for the solid phase also, regardless of the number of the fluid phases in the system, the same density equation as Eq. (3-12). For the pressure transferred to the solid phase, 𝑃𝑆, the  pressure of the equivalent single phase matrix, 𝑃𝑚 in Eq. (3-61), is considered based on the Bishop’s model, i.e. 𝑃𝑆 = 𝑃𝑚. For the Bishop’s parameter, the simplest choice of ℬ = 𝜓𝐹 is applied as discussed in Section 2.1.3 where 𝜓𝐹 = 𝜑𝐹/(𝜑𝐹 + 𝜑?̂?) is called the degree of saturation of the fluid phase, 𝐹.  Substituting Eqs. (3-10)-(3-14) in Eq. (3-56) and based on the same arguments used as those for two-phase system the solidifying matrix phase in the two-phase system, we obtain: 57  ?̇? = (1 − 𝑏)?̇?𝑖,𝑖 −𝑏 − 𝜑𝐾𝑆?̇?𝐹𝑃𝐹 −𝑏 − 𝜑𝐾𝑆?̇??̂?𝑃?̂? − (𝜑𝐹𝐾𝐹+𝜓𝐹(𝑏 − 𝜑)𝐾𝑆) ?̇?𝐹− (𝜑?̂?𝐾?̂?+𝜓?̂?(𝑏 − 𝜑)𝐾𝑆) ?̇??̂? + 𝜑𝐹(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗+ 𝜑?̂?𝛿𝑖𝑗𝜀?̂̇?𝑓𝑖𝑗+ (1 − 𝜑)(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗+13𝐾𝑆𝛿𝑖𝑗𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝑓𝑘𝑙= 0 (3-62) where 𝜓?̂? = 1 − 𝜓𝐹. It is noted that all evolving properties of the solid and solid-skeleton discussed in Section 3.1.2.1.2 and shown in Eq. (3-17) are also applicable to the three-phase system.   In deriving Eq. (3-62), in the current 3IFS model it is assumed that only one of the two fluid phases (i.e. resin) undergoes curing during the process while the other fluid phase (e.g. air) remains unsolidified. The gradual transformation of the resin fluid phase, 𝐹, from an unsolidified to a fully solidified state affects the free volumetric strain rate of the system. However, it is possible to apply this approach to a composite system with more than one solidifying fluid phase as well each having an independent solidification behaviour described by two distinct solidification factors, namely 𝜆𝐹 and 𝜆?̂?. However, as the objective of the current three-phase model is inclusion of the gas phase as one of the two fluid phases, which does not solidify, only one solidification factor is considered in the model.  At the final stages of processing the composite is a three-phase solid material and its properties in Eq. (3-17), i.e. 𝐃𝑐, 𝛼𝑐𝑡ℎ, and 𝛼𝑐𝑐𝑠 must be obtained considering all three phases in a solid composite material, e.g. through a three-phase analytical or numerical micromechanics formulation.  58   3.2.2.2 Effective Stress  Considering the degree of saturation of the fluid phase 𝐹 as Bishop’s parameter, i.e. ℬ = 𝜓𝐹 in Eq. (3-61), and using Eq. (3-60) yields the following Bishop’s generalized effective stress equation for the load-sharing of the three-phase system: 𝜎𝑖𝑗 = 𝜎𝑆𝐾𝑖𝑗 − 𝑏(𝑃?̂? − 𝜓𝐹(𝑃?̂? − 𝑃𝐹))𝛿𝑖𝑗 (3-63) To consider nonlinearity due to material or geometrical effects, the effective stress (Eq. (3-63)) needs to be expressed in its differential, or rate, form and then solved incrementally/iteratively within a nonlinear finite element analysis procedure: ?̇?𝑖𝑗 = ?̇?𝑆𝐾𝑖𝑗 − 𝑏 (?̇??̂? − 𝜓𝐹(?̇??̂? − ?̇?𝐹) − ?̇?𝐹(𝑃?̂? − 𝑃𝐹)) 𝛿𝑖𝑗− ?̇?(𝑃?̂? − 𝜓𝐹(𝑃?̂? − 𝑃𝐹))𝛿𝑖𝑗 (3-64) As shown schematically in Figure 3.11 for a three-phase system, as the resin material solidifies, the matrix (denoted by 𝑚 and representing a combination of 𝐹 and ?̂?) contributes to the load carrying capacity of the solid-skeleton to finally form a three-phase solid composite material. As the resin material undergoes solidification, a portion of the equivalent matrix pressure, ?̇?𝑃𝑚 (with 𝑃𝑚 defined in Eq. (3-61)), should be added incrementally to the stress carried by the solid-skeleton. This transformation of pressure to stress due to solidification of the matrix is identical to the two-phase system where the matrix phase included only a single fluid phase 𝐹 (see Section 3.1.2.2.2). As a result of this transformation in the three-phase system, ?̇?𝑆𝐾 is related to the mechanical strain in the solid-skeleton through its elastic constitutive stiffness tensor, 𝐷𝑆𝐾, as follows: ?̇?𝑆𝐾𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙+ ?̇? (𝑃?̂? − 𝜓𝐹(𝑃?̂? − 𝑃𝐹))𝛿𝑖𝑗 (3-65) 59  Combining Eqs. (3-64) and (3-65) yields the differential (rate) form of the total stress in terms of the skeleton strain and fluids’ pressures as ?̇?𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝜎𝑘𝑙− 𝑏 (?̇??̂? − 𝜓𝐹(?̇??̂? − ?̇?𝐹) − ?̇?𝐹(𝑃?̂? − 𝑃𝐹)) 𝛿𝑖𝑗 (3-66)  3.2.3 Summary of the Governing Equations (3IFS-Isotropic) Taking into account the constitutive equations and the conservation equations, the governing equations for the 3IFS model written in terms of the primary variables, 𝑢, 𝑣𝐹, 𝑣?̂?, 𝑃𝐹, and 𝑃?̂? can be summarized as follows:  Mass conservation: 𝑏?̇?𝑖,𝑖 + 𝑣𝐹𝑖,𝑖 + 𝑣?̂?𝑖,𝑖 +𝑏 − 𝜑𝐾𝑆?̇?𝐹𝑃𝐹 +𝑏 − 𝜑𝐾𝑆?̇??̂?𝑃?̂? + (𝜑𝐹𝐾𝐹+𝜓𝐹(𝑏 − 𝜑)𝐾𝑆) ?̇?𝐹+ (𝜑?̂?𝐾?̂?+𝜓?̂?(𝑏 − 𝜑)𝐾𝑆) ?̇??̂? − 𝜑𝐹(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗− 𝜑?̂?𝛿𝑖𝑗𝜀?̂̇?𝑓𝑖𝑗− (1− 𝜑)(1 − 𝜆)𝛿𝑖𝑗𝜀?̇?𝑓𝑖𝑗−13𝐾𝑆𝛿𝑖𝑗𝐷𝑆𝐾𝑖𝑗𝑘𝑙𝜀?̇?𝐾𝑓𝑘𝑙= 0 (3-67)  Darcy’s flow for fluid 𝐹: 𝑃𝐹,𝑖 + 𝜇𝐹𝑆𝐹𝑖𝑗−1𝑣𝐹𝑗 = 0 (3-68)  Darcy’s flow for fluid ?̂?: 𝑃?̂? ,𝑖 + 𝜇?̂?𝑆?̂?𝑖𝑗−1𝑣?̂?𝑗 = 0 (3-69)  Pressure equation: 𝑃?̂? − 𝑃𝐹 = 𝑃𝑐 (3-70)  Equilibrium equation: 𝜎𝑆𝐾𝑖𝑗,𝑗 − 𝑏 (𝑃?̂? ,𝑖 − 𝜓𝐹 (𝑃?̂? ,𝑖 − 𝑃𝐹,𝑖)) = 0 (3-71) 60  with the rate of stress in the solid-skeleton given by: ?̇?𝑆𝐾𝑖𝑗 = 𝐷𝑆𝐾𝑖𝑗𝑘𝑙 (12(?̇?𝑘,𝑙 + ?̇?𝑙,𝑘) − 𝜀?̇?𝐾𝑓𝑘𝑙) + ?̇? (𝑃?̂? − 𝜓𝐹(𝑃?̂? − 𝑃𝐹))𝛿𝑖𝑗 (3-72) where 𝑏 = 1 − 𝐾𝑆𝐾/𝐾𝑆, 𝐃𝑆 = 𝑓(𝐃𝑓 , 𝐃𝑐, 𝜆) , 𝐃𝑆𝐾 = 𝑓′(𝐃𝑓𝑏 , 𝐃𝑐, 𝜆), 𝛼𝑆𝐾𝑡ℎ = 𝑔(𝛼𝑓𝑏𝑡ℎ, 𝛼𝑐𝑡ℎ, 𝜆), and 𝛼𝑆𝐾𝑐𝑠 = 𝑔′(𝛼𝑐𝑐𝑠, 𝜆). In the special case when 𝜆 = 0, all solid and solid-skeleton properties are the same as those of the fibre and fibre-bed, respectively, while at the other extreme when 𝜆 = 1, they both have properties of the three-phase fully solid composite material for which 𝑏 = 0. This leads to a decoupled equilibrium equation of the system (Eq. (3-29) which is identical to Eq. (3-37) noting that 𝑏 = 0 and 𝜎𝑆𝐾 = 𝜎𝑐 in this case).   3.2.4 Finite Element Implementation (3IFS-Isotropic) The finite element implementation of the 3IFS methodology developed in the preceding section can be achieved by converting the governing differential equations (Eqs. (3-67)-(3-71)) into their weak form using the Galerkin method. Within a 𝑢-𝑣-𝑃 FE framework adopted here, the system displacement vector, 𝑢𝑖, relative velocity vectors of each fluid phase, 𝑣𝐹𝑖 and 𝑣?̂?𝑖, and fluid pore pressures, 𝑃𝐹 and 𝑃?̂?, are discretized spatially. As shown in Appendix A.2, applying the backward Euler scheme for temporal integration combined with the Newton’s nonlinear solution method, the matrix equation for the 𝑘𝑡ℎ iteration within the 𝑛𝑡ℎ time step can be written as: [      ∆𝑡𝐊𝒖𝒖 𝟎 𝟎 ∆𝑡𝐊𝒖𝑷𝐹 ∆𝑡𝐊𝒖𝑷?̂?𝟎 ∆𝑡𝐊𝒗𝐹𝒗𝐹 𝟎 ∆𝑡𝐊𝒗𝐹𝑷𝐹 𝟎𝟎 𝟎 ∆𝑡𝐊𝒗?̂?𝒗?̂? 𝟎 ∆𝑡𝐊𝒗?̂?𝑷?̂?𝟎 𝟎 𝟎 ∆𝑡𝐊𝑷𝐹𝑷𝐹 ∆𝑡𝐊𝑷𝐹𝑷?̂?𝐂𝑷?̂?𝒖 ∆𝑡𝐊𝑷?̂?𝒗𝐹 ∆𝑡𝐊𝑷?̂?𝒗?̂? ∆𝑡𝑲𝑷?̂?𝑷𝐹 + 𝑪𝑷?̂?𝑷𝐹 ∆𝑡𝑲𝑷?̂?𝑷?̂? + 𝑪𝑷?̂?𝑷?̂?]      𝑘−1{      ∆?̅?𝑘?̅?𝐹𝑘?̅??̂?𝑘∆?̅?𝐹𝑘∆?̅??̂?𝑘}       (3-73) 61  ={      Δ𝑡[𝐟𝒖 + 𝐟𝒖𝑓𝑆𝐾𝑘−𝐟𝒖𝑖𝑛𝑡𝑘−1−𝐟𝒖𝑏𝑃𝑘−1 (∆𝜓?̂?𝑘−1?̅??̂?𝑘−1+ ∆𝜓𝐹𝑘−1?̅?𝐹𝑘−1)]∆𝑡[𝐟𝒗𝐹 − 𝐊𝒗𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1]∆𝑡[𝐟𝒗?̂? − 𝐊𝒗?̂?𝑷?̂?𝑘−1?̅??̂?𝑘−1]∆𝑡[𝐟𝑷𝐹 − 𝐊𝑷𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1− 𝐊𝑷𝐹𝑷?̂?𝑘−1?̅??̂?𝑘−1]𝐅𝑷?̂? }       where 𝐅𝑷?̂? = ∆𝑡[𝐟𝑷?̂? − 𝐊𝑷?̂?𝑷𝐹𝑘−1?̅?𝐹𝑘−1− 𝐊𝑷?̂?𝑷?̂?𝑘−1?̅??̂?𝑘−1] + 𝐂𝑷?̂?𝒖𝑘−1(?̅?0 − ?̅?𝑘−1)+ 𝐂𝑷?̂?𝑷𝐹𝑘−1(?̅?𝐹0− ?̅?𝐹𝑘−1) + 𝐂𝑷?̂?𝑷?̂?𝑘−1(?̅??̂?0− ?̅??̂?𝑘−1) (3-74) and ∆( )𝑘 = ( )𝑘 − ( )𝑘−1. A schematic of the algorithm adopted for implementation of the three-phase IFS model in the MATLAB® code is shown in Figure 3.12.  As in the previous 2IFS case, we use the same low order bi-linear two-dimensional (plane) isoparametric element in which the degrees of freedom of the system displacements and relative velocities of the fluid phases, i.e. ?̅?, ?̅?𝐹, and ?̅??̂?, respectively, are assigned to the corner nodes and the fluid pressures, ?̅?𝐹 and ?̅??̂?, are assigned to the internal (central) node. All matrix and vector quantities in Eq. (3-73) are defined in Table 3.2.  Using the new values of ?̅?, ?̅?𝐹, ?̅??̂?, ?̅?𝐹, and ?̅??̂? from Eq. (3-73) in the form of Eq. (3-66), the stress is updated as follows: 𝛔𝑘 = 𝛔𝑘−1 + 𝐃𝑆𝐾𝑘−1𝐁𝑘∆?̅?𝑘− 𝐃𝑆𝐾𝑘−1Δ𝛆𝑆𝐾𝑓 𝑘−𝑏𝑘−1𝛅 [∆?̅??̂?𝑘− 𝜓𝐹𝑘−1 (∆?̅??̂?𝑘− ∆?̅?𝐹𝑘)− ∆𝜓𝐹𝑘−1(?̅??̂?𝑘− ?̅?𝐹𝑘)] (3-75) 62   3.2.4.1 Volume Fraction Calculation Using Eq. (3-50) for each fluid phase separately, the discretized form of the equation for computing the volume fraction can be written as:  𝜑Φ𝑘 =𝜓Φ0𝜑0(exp(− ?̅?Φ𝑘−?̅?Φ𝟎𝐾Φ)+𝛅𝑇Δ𝛆Φ𝑓 𝑘)−Δ𝑡𝛅𝑇𝐁𝑘?̅?Φ𝑘1+𝛅𝑇𝐁𝑘(?̅?𝑘−?̅?0),   Φ ≡ 𝐹 or ?̂?  (3-76)  It is noted that again the volume fractions of constituents are updated (based on Eq. (3-76)) only up to the solidification point after which they are kept constant according to classical solid mechanics micromechanics formulations wherein volume fraction of phases are independent of phase deformations. Also, once the solidification factor 𝜆 reaches one, the pressure of the gas phase is kept constant as no volume change for gas phase is considered after full solidification. Also, for the free strain of the solid phase, 𝛆𝑆𝑓, and the second fluid phase, 𝛆?̂?𝑓, the cure shrinkage component of the free strain are zero as these two phases do not undergo curing.  3.3 Summary  In this chapter, the formulation and FE implementation of the IFS-isotropic model is presented. First, the two-phase model formulation comprising of one fluid phase flowing through a porous solid medium, i.e. liquid which can cure/solidify or gas, is presented. Conservation equations are presented and then followed by constitutive equations developed in a way to integrate the flow and stress development regimes in a seamless fashion. The basis of IFS model and its compatibility with the fundamental physics in processing of polymeric composite materials, such 63  as transformation of pressure into residual stresses, is discussed in detail. Then, the two-phase model is extended to the three-phase IFS model including two fluid phases, i.e. liquid and gas, flowing simultaneously and independently while the liquid phase can solidify (e.g. as a result of curing). In the next chapter, the orthotropic extension of both the two- and three-phase IFS models are presented.  64  3.4 Tables  𝐊𝒖𝒖 = ∫𝐁T Ω𝐃𝑆𝐾𝐁𝑑Ω 𝐊𝒖𝑷 = −𝑏 ∫𝐁T Ω𝛅𝐍𝑃𝑑Ω 𝐊𝑷𝒗 = ∫𝐍𝑃𝑇 Ω𝛅𝑇𝐁𝑑Ω 𝐊𝒗𝒗 = 𝜇𝐹 ∫𝐍𝑑𝑇𝐒−1𝐍𝑑 Ω𝑑Ω 𝐊𝒗𝑷 = −∫𝐁𝑇 Ω𝛅𝐍𝑃𝑑Ω  𝐂𝑷𝒖 = 𝑏 ∫𝐍𝑃𝑇 Ω𝛅𝑇𝐁𝑑Ω 𝐂𝑷𝑷 = (𝜑𝐹𝐾𝐹+𝑏 − 𝜑𝐹𝐾𝑆) ∫𝐍𝑃𝑇 Ω𝐍𝑃𝑑Ω  𝐟𝒖 = − ∫𝐍𝒅𝑻 Γ𝜎𝐭 𝑑Γ 𝐟𝒖𝑖𝑛𝑡 = ∫𝐁𝑇𝛔𝑑Ω Ω 𝐟𝒖𝑓𝑆𝐾= ∫𝐁T Ω𝐃𝑠𝑘d𝛆SK𝑓 dΩ 𝐟𝒗 = − ∫𝐍𝑑𝑇 Γ𝑃𝑃𝑏𝑑Γ 𝐟𝑷 = ∫𝛅𝑇(𝜑𝐹(1 − 𝜆)?̇?𝐹𝑓 + (1 − 𝜑𝐹)(1 − 𝜆)?̇?𝑆𝑓 +13𝐾𝑆𝐃𝑆𝐾?̇?𝑆𝐾𝑓 ) Ω𝐍𝑃𝑇𝑑Ω  𝐍𝑑 = [𝐍 𝟎𝟎 𝐍], where 𝐍 is the matrix of displacement or velocity shape functions 𝐍𝑃 = Matrix of pressure shape functions 𝐁 = Matrix of spatial derivatives of the displacement or velocity shape functions 𝛅𝑇 = {1 1 0}  𝐭 = Vector of total tractions applied on the traction boundaries (Γ𝜎) of the system 𝑃𝑏 = Applied pressure to the fluid at permeable boundaries (Γ𝑃) of the system  Table 3.1 Definition of matrices and vectors used in the two-phase isotropic FE model         65  𝐊𝒖𝒖 = ∫𝐁T Ω𝐃𝑆𝐾𝐁𝑑Ω 𝐊𝒖𝑷Φ = −𝑏𝜓Φ ∫𝐁T Ω𝛅𝐍𝑃𝑑Ω 𝐊𝑷Φ𝒗Φ = ∫𝐍𝑃𝑇 Ω𝛅𝑇𝐁𝑑Ω 𝐊𝒗Φ𝒗Φ = 𝜇Φ ∫𝐍𝑑𝑇𝐒Φ−1𝐍𝑑 Ω𝑑Ω 𝐊𝒗Φ𝑷Φ = −∫𝐁𝑇 Ω𝛅𝐍𝑃𝑑Ω  𝐊𝑷𝐹𝑷𝐹 = −𝐊𝑷𝐹𝑷?̂? = −∫𝐍𝑃𝑇 Ω𝐍𝑃𝑑Ω 𝐊𝑷?̂?𝑷Φ =∆𝜓Φ (𝑏 − 𝜑)𝐾𝑆∫𝐍𝑃𝑇 Ω𝐍𝑃𝑑Ω 𝐂𝑷?̂?𝒖 = 𝑏 ∫𝐍𝑃𝑇 Ω𝛅𝑇𝐁𝑑Ω 𝐂𝑷?̂?𝑷Φ = (𝜑Φ𝐾Φ+𝜓Φ(𝑏 − 𝜑)𝐾𝑆) ∫𝐍𝑃𝑇 Ω𝐍𝑃𝑑Ω 𝐟𝒖 = − ∫𝐍𝒅𝑻 Γ𝜎𝐭 𝑑Γ𝜎 𝐟𝒖𝑖𝑛𝑡 = ∫𝐁𝑇𝛔𝑑Ω Ω 𝐟𝒖𝑓𝑆𝐾= ∫𝐁𝑇 Ω𝐃𝑠𝑘𝑑𝛆SK𝑓 dΩ 𝐟𝒖𝑏𝑷 = 𝑏 ∫𝐁𝐓 Ω𝛅𝐍𝑷𝑑Ω 𝐟𝑷𝐹 = ∫𝑃𝑐𝑑Ω  Ω 𝐟𝒗𝛽 = − ∫𝐍𝑑𝑇 Γ𝑃𝑃𝑏Φ𝑑Γ𝑃 𝐟𝑷 ?̂? = ∫𝛅𝑇(𝜑𝐹(1 − 𝜆)?̇?𝐹𝑓 + 𝜑?̂??̇??̂?𝑓 + (1 − 𝜑)(1 − 𝜆)?̇?𝑆𝑓 +13𝐾𝑆𝐃𝑆𝐾?̇?𝑆𝐾𝑓 ) Ω𝐍𝑃𝑇𝑑Ω Φ = 𝐹 or ?̂?   Table 3.2 Definition of matrices and vectors used in the three-phase isotropic FE model     66  3.5 Figures   Figure 3.1 A representation of the two-phase composite system and its constituents      Figure 3.2 A representation of axial and transverse flow of a single fluid through fibrous porous medium  Unsolidified Fluid ( )Plane of isotropy (transverse to the fibres) Longitudinal plane of the fibresPorous structure:Fibre-bed ( ) / Solid Skeleton ( )Fibres ( ) / Solid ( )axial flowtransverse flowFluidFibres67    Figure 3.3 Evolution of the properties of the solid/solid-skeleton and fluid phases with increasing solidity of the fluid phase  a)   b)  Figure 3.4 Conceptual example showing the effect of liquid resin pressure on the stress in the solidified resin after curing Biot’s Model Solid Mechanics, fibre ( ) , fibre-bed ( )increasing Homogenized composite material ( )composed of solid grains andsolidified fluidSolid grains ( )forming porous solid-skeleton ( )Unsolidified fluid ( )Liquid resin PressurefSolid resin StressfCuring (solidification)Liquid resin PressureLiquid resin Temp Curing (solidification)Solid resin Stress68    Figure 3.5 Schematic representation of 4-1 element                     and nodesnode69     Figure 3.6 Algorithm for the finite element implementation of the two-phase IFS model  Solidification factorCuring Fluid ( ) Degree of CureFluids ViscosityPermeability TensorsFluid properties, )Composites properties, , )Skeleton properties)Fibre-bed properties( )Fibre properties, , )Biot’s Coefficient)Convergence Criteria, Converged?Solid properties( )matrices for each elementand matricesfor systemRHS and LHS matrices in solution methodNew , , valuesPhase Volume FractionsTwo-phase MicromechanicsAssemblyInternal ForcesStress UpdateYesNoVolume FractionUpdateBoundary Conditions Updatevectors for each elementvectotr for systemAssemblyTemperature Update( )Initial Values(Temperature, Phase Volume Fractions, Mechanical Loading, Geometry)Is ?YesNoSimulation FinishedStrat of simulation70   a)  b)  c)  Figure 3.7 Comparison of volume fraction equations based on small or large strain assumptions: a) expansion of fluid phase b) expansion of solid phase c) flow of the fluid into the system   FluidFluid expands to twice of its volumeFluidSolid SolidFluidSolid expands to twice of its volumeFluidSolidSolidFluidFluid has twice of its original volume due to flowFluidSolid Solid71   Figure 3.8 A representation of three-phase composite system and its constituents       Figure 3.9 A representation of decoupled multiphase fluids axial and transverse flow through fibrous porous medium   Fluid ( )Porous structure:Fibre-bed ( ) / Solid Skeleton ( )Fibres ( ) / Solid ( )Fluid ( )Plane of isotropy (transverse to the fibres) Longitudinal plane of the fibresFluid axial flowFluid transverse flowFluid axial flowFluid transverse flowFluid FibresFluid 72    Figure 3.10 Representation of three-phase composite system and the fluids flow mechanism as fluid 𝑭 solidifies    Figure 3.11 Evolution of the properties of the solid/solid-skeleton (𝑺/𝑺𝑲) and multiphase matrix phase (combined of fluid phases 𝑭 and ?̂?) with increasing solidity of the fluid phase   Fluid ( )Solid grains ( )forming porous solid-skeleton ( )Unsolidified fluid ( )Fluid ( )evolution of Bshop’s Model Solid Mechanics, fibre ( ) , fibre-bed ( )Homogenized three-phase composite material ( )Multiphase matrix ( ) Multiphase matrix ( ) Multiphase matrix ( )Solid grains ( )increasing 73   Figure 3.12 Algorithm for the finite element implementation of the three-phase IFS model  Solidification factorCuring Fluid ( ) Degree of CureFluids ViscosityPermeability TensorsFluid properties, )Composites properties, , )Skeleton properties)Fibre-bed properties( )Fluid properties, )Fibre properties, , )Biot’s Coefficient)Convergence Criteria, Converged?Solid properties( )matrices for each elementand matricesfor systemRHS and LHS matrices in solution methodNew , , valuesPhase Volume FractionsThree-phase MicromechanicsAssemblyInternal ForcesStress UpdateYesNoVolume FractionUpdateBoundary Conditions Updatevectors for each elementvectotr for systemAssemblyTemperature Update( )Initial Values(Temperature, Phase Volume Fractions, Mechanical Loading, Geometry)Is ?YesNoSimulation FinishedStrat of simulationCapillary pressure)Surface tension74  Chapter 4: Orthotropic Integrated Flow-Stress Model (IFS-Orthotropic) Fibrous composite materials are intrinsically orthotropic. Although a unidirectional composite with perfectly aligned fibres is transversely isotropic, in reality due to waviness of fibres, the fibre-bed and consequently the composite material, show orthotropic characteristics.  Furthermore, the carbon fibres that are typically used in prepreg composites, i.e. the solid phase in the IFS model, also show transversely isotropic behaviour. Therefore, for a broader application of the model, the isotropic IFS model developed in the previous chapter is extended to the more general case when both the solid and solid-skeleton show orthotropic properties.   Extending the isotropic model requires a broader definition of material properties in the poroelasticity model as well as their definition in the presented IFS model. For instance, Biot’s coefficient discussed in previous chapters represents deformability of solid grains versus pores deformation in a porous medium under loading. In a system where the solid phase and the structure made out of it, solid-skeleton, are both isotropic, the direction of deformation does not matter. Therefore, a scalar can represent the relative deformability of pores and solid grains. However, for an orthotropic system, not only the deformability in different directions are different but also they interact as well. As a result, the scalar Biot’s coefficient used widely in the poroelasticity models needs to be extended to a 6 × 1 vector (or an equivalent symmetric 3 × 3 matrix) in a three-dimensional model. This vector would depend on the full elasticity stiffness tensor of the constituents, solid and solid-skeleton, as there is no explicit definition of bulk modulus in orthotropic materials. The vectorial representation of the Biot’s coefficient is: 75  𝐛 =[     𝑏11𝑏22𝑏33𝑏23𝑏31𝑏12]      (4-1) The first three components of this vector (or matrix if written as a square array) represents the Biot’s coefficient in each direction while the other three (off-diagonal terms in Biot’s coefficient matrix) represents the interplay influence coefficient between different directions. For the special case of an isotropic system, the Biot’s coefficient is the same in every direction while the three off-diagonal terms are zero. Consequently, a scalar coefficient is enough for representing the material system in the poroelasticity model: 𝐛 =[     𝑏𝑏𝑏000]     → 𝑏 (4-2) In a general orthotropic IFS model, an extension of the methodology used for evolution of material properties, and consequently Biot’s coefficient, as a result of solidification of the fluid phase is also required.   In this chapter, first, the orthotropic formulations of quantities and variables, such as Biot’ coefficient, rate of change of compressibility, etc., are developed. For each case, the developed formulation is evaluated for the special case of isotropic constituents and compared to equations presented in the previous chapter. Then the orthotropic version of the governing equations for the IFS model is presented briefly, and finally the matrix equations required for FE implementation of the model are derived.  76     4.1 Biot’s Coefficient for Orthotropic System In this section, the Biot’s coefficient for the general case when both the solid and solid-skeleton have orthotropic properties is obtained. First, the general Biot’s coefficient is derived and then it is evaluated for the special case of isotropic solid and solid-skeleton through comparison with the isotropic expressions presented before in Eq. (2-3) and used in the IFS-Isotropic model in Chapter 3. In order to simplify the form of the equations, we will adopt the contracted boldface notation for the tensorial (matrix and vector) quantities rather than the index notation used in the previous chapter.  In the general case of orthotropic solid and solid-skeleton constituents, the total stress 𝛔 for a composite system with liquid fluid phase is defined in terms of either 𝛔𝑆𝐾 (same as 𝛔′′ used in soil mechanics)  or 𝛔′ effective stresses as follows [13,69]: 𝛔 = 𝛔′ − 𝑃𝑆𝛅 (4-3) 𝛔 = 𝛔𝑆𝐾 − 𝑃𝑆𝐛 (4-4) where 𝛔′ and 𝛔𝑆𝐾 are effective stress vectors, 𝐛 is the vector of Biot’s coefficients, 𝑃𝑆 is the  hydrostatic pressure transmitted to the solid phase from the fluid phase, and 𝛅 = [1 1 1 0 0 0]𝐓 is the identity vector. Effective stresses 𝛔′ and 𝛔𝑆𝐾 can be expressed as:  𝛔′ = 𝐃𝑆𝐾(𝛆 − 𝛆𝑃𝑆) (4-5) 𝛔𝑆𝐾 = 𝐃𝑆𝐾𝛆 (4-6) 77  In Eqs. (4-5) and (4-6), 𝛆 is the total strain vector of the system (associated with the solid-skeleton) and 𝛆𝑃𝑆 is defined as the solid phase strain as a result of the hydrostatic pressure transferred to the solid, 𝑃𝑠, for the general case of an orthotropic solid phase (𝑆) written as: 𝛆𝑃𝑆 = −𝑃𝑆𝐃𝑆−1𝛅 (4-7) Using Eqs. (4-3)- (4-7), the vector of Biot’s coefficients for the general case of orthotropic materials can be written as: 𝐛 = 𝛅 − 𝐃𝑆𝐾𝐃𝑆−1𝛅 (4-8) The above equation is in agreement with the expression presented by Carroll [71] in index notation where homogeneous surface tractions were used to obtain the effective stress law of the elastic porous structure. Carroll [71] and Thompson and Willis [72] provided further insights in the anisotropic constants of poroelastic constitutive constants obtained twenty-eight years earlier by Biot [73] and Biot and Willis [12]. Further physical insight of the anisotropic constants of poroelastic theory and their engineering implication as well as a list of previous works in this area can be found in the work by Cheng [74].  4.1.1 Special Case of Isotropic Materials The general Biot’s coefficients presented in Eq. (4-8) is evaluated and checked for the special case of isotropic phases. When the solid phase has isotropic properties, the general form of the solid phase strain as result of hydrostatic pressure, Eq. (4-7), is reduced to 𝛆𝑃𝑆 = −𝑃𝑆𝐃𝑆−1𝛅 = −13𝐾𝑆𝑃𝑆𝛅 (4-9) Taking Eq. (4-9) into account along with assuming an isotropic solid-skeleton stiffness matrix, 𝐃𝑆𝐾, in Eq. (4-8) yields the Biot’s coefficient as 78  𝐛 = (1 −𝐾𝑆𝐾𝐾𝑆)𝛅 (4-10) which is in identical to Eq. (2-3) for isotropic materials.  4.2 Rate of Change of Compressibility, ?̇?, for Orthotopic System Density of a solid phase 𝑆 in a composite system varies as a result of two different volumetric strains namely, those associated with mechanical loading and free strains: 1𝜌𝑆?̇?𝑆 = −?̇?𝛔 − ?̇?𝑆𝑓 (4-11) where ?̇?𝛔 and ?̇?𝑆𝑓 are the rate of mechanical and free volumetric strain, respectively. For free volumetric strain, we simply have: ?̇?𝑆𝑓 = 𝛅𝐓?̇?𝑆𝑓 (4-12) The mechanical volumetric strain, ?̇?𝛔, can be decomposed into the hydrostatic pressure and effective stress 𝛔′ (as defined in Eq. (4-5)) components as: ?̇?𝜎 = ?̇?𝛔′+ ?̇?𝑃𝑆 (4-13) Taking Eq. (4-7) into account and assuming that in a two-phase system the total fluid pore pressure is transmitted to the solid phase, i.e. 𝑃𝑆 = 𝑃𝐹, we can write:  ?̇?𝑃𝑆 = −?̇?𝐹𝛅T𝐃𝑆−1𝛅 (4-14) Also, for ?̇?𝛔′ we have: ?̇?𝛔′=1(1 − 𝜑)𝛅T?̇?𝛔′ (4-15) 79  where (1 − 𝜑) is the volume fraction of the solid phase (in a two-phase system) and the strain due to 𝛔′ can be easily obtained from the stiffness properties of the orthotropic solid phase as follows: 𝛆𝛔′= 𝐃𝑆−1𝛔′ (4-16) Using Eqs. (4-14)-(4-16) in Eq. (4-13) and substituting the resulting expression for ?̇?𝜎 into Eq. (4-11) and making use of Eq. (4-12), we arrive at: 1𝜌𝑆?̇?𝑆 = −1(1 − 𝜑)𝛅T𝐃𝑆−1?̇?′ + ?̇?𝐹𝛅T𝐃𝑆−1𝛅 − 𝛅T?̇?𝑆𝑓 (4-17) Considering Eq. (4-5) for 𝛔′ along with Eq. (4-7) as well as the orthotropic Biot’s coefficients as presented in Eq. (4-8) and taking advantage of the symmetry of both 𝐃𝑆𝐾 and 𝐃𝑆, we obtain: 1𝜌𝑆?̇?𝑆 = −1(1 − 𝜑)𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾?̇? + ?̇?𝐹𝛅𝐓 (𝐃𝑆−1𝛅 −1(1 − 𝜑)𝐃𝑆−1(𝛅 − 𝐛)) − 𝛅𝐓?̇?𝑆𝑓 (4-18)  The total strain of the system, 𝛆, can be defined in terms of the displacement filed as: 𝛆 = 𝛁T𝐮 (4-19) where 𝐮 is the 3 × 1 displacement vector and 𝛁 is the gradient operator defined in the following matrix form: 𝛁 =[      𝜕𝜕𝑥0 0 0𝜕𝜕𝑧𝜕𝜕𝑦0𝜕𝜕𝑦0𝜕𝜕𝑧0𝜕𝜕𝑥0 0𝜕𝜕𝑧𝜕𝜕𝑦𝜕𝜕𝑥0]       (4-20) Simplifying Eq. (4-18) and using Eq. (4-19), we can write: (1 − 𝜑)𝜌𝑆?̇?𝑆 = 𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾𝛁T?̇? + ?̇?𝐹𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅) − (1 − 𝜑)𝛅𝐓?̇?𝑆𝑓 (4-21) 80   In the rate of change of compressibility, the effect of the solid-skeleton’s free strain should also be taken into account [13] which for the orthotropic solid-skeleton can be written as follows: ?̇?𝑆𝐾𝑓 = 𝛅𝐓𝐃𝑆−1𝛔𝑆𝐾𝑓  (4-22) where 𝛔𝑆𝐾𝑓 is the stress associated with free strains in  the skeleton (this would be the negative of the stresses that would develop if all the free strains were to be fully constrained): 𝛔𝑆𝐾𝑓 = 𝐃𝑆𝐾𝛆𝑆𝐾𝑓  (4-23) So, ?̇?𝑆𝐾𝑓 = 𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾𝛆𝑆𝐾𝑓  (4-24)  Using the vector form of the rate of change of the fluid phase density, Eq. (3-10), along with Eqs. (4-21) and (4-24) in the rate of change of compressibility of the system, Eq. (3-5),  the final form of the rate of change of compressibility of a two-phase system (with orthotropic solid and solid-skeleton) is given by: ?̇? = (𝛅 − 𝐛)T𝛁T?̇? − ?̇?𝐹(𝜑𝐾𝐹+ 𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)) + 𝜑𝛅𝐓?̇?𝐹𝑓 + (1 − 𝜑)𝛅𝐓?̇?𝑆𝑓+ 𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾𝛆𝑆𝐾𝑓 (4-25) For the case when the fluid phase remains fluid and does not evolve in to a solid.  For a three-phase system with two fluid phases, 𝐹 and ?̂?, this equation is extended to: 81  ?̇? = (𝛅 − 𝐛)T𝛁T?̇? − ?̇?𝐹𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)𝑃𝐹 − ?̇??̂?𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)𝑃?̂?− ?̇?𝐹 (𝜑𝐹𝐾𝐹+ 𝜓𝐹𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)) − ?̇??̂? (𝜑?̂?𝐾?̂?+ 𝜓?̂?𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅))+ 𝜑𝐹𝛅𝐓?̇?𝐹𝑓 +𝜑?̂?𝛅𝐓?̇??̂?𝑓 + (1 − 𝜑)𝛅𝐓?̇?𝑆𝑓 + 𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾𝛆𝑆𝐾𝑓 (4-26)  4.2.1 Special Case of ?̇? for Isotropic Materials Considering the fact that the Biot’s coefficient is a scalar for isotropic materials, we can write: 𝛅𝐓𝐃𝑆−1𝐛 = 𝑏𝛅𝐓𝐃𝑆−1𝛅 (4-27) (𝛅 − 𝐛)T𝛁T?̇? = (1 − 𝑏)𝛅𝐓𝛁T?̇? (4-28) For an isotropic solid phase, Eqs. (4-27) and (4-28) simplify to: 𝑏𝛅𝐓𝐃𝑆−1𝛅 =𝑏𝐾𝑆 (4-29) Also we can write 𝜑𝛅𝐓𝐃𝑆−1𝛅 =𝜑𝐾𝑆 (4-30) The volumetric strain in the solid-skeleton due to free strain, Eq. (4-22), can be written as follows when the solid phase is isotropic  ?̇?𝑆𝐾𝑓 =13𝐾𝑆𝛅𝐓𝛔𝑆𝐾𝑓  (4-31) where 𝛔𝑆𝐾𝑓 is defined in Eq. (4-23). Using Eqs. (4-29)-(4-31) in Eq. (4-25), gives the rate of  change of the compressibility for a composite system with isotropic solid and solid-skeleton constituents as: 82  ?̇? = (1 − 𝑏)𝛅𝐓𝛁T?̇? − (𝜑𝐾𝐹+(𝑏 − 𝜑)𝐾𝑆) ?̇?𝐹 + 𝜑𝛅𝐓?̇?𝐹𝑓 + (1 − 𝜑)𝛅𝐓?̇?𝑆𝑓+13𝐾𝑆𝛅𝐓𝐃𝑆𝐾  ?̇?𝑆𝐾𝑓  (4-32) which in index notation is identical to the expression introduced by Zienkiewicz and Shiomi [13] (see Eq. (3-16)). The last term in Eq. (4-32) is shown in terms of 𝐃𝑆𝐾 to maintain its similarity with the work by Zienkiewicz and Shiomi [13] while it can also be written in a more simplified form of 𝐾𝑆𝐾𝐾𝑆𝛅𝐓 ?̇?𝑆𝐾𝑓 for an isotropic skeleton. For the three-phase system, following the same steps, it can be easily shown that Eq. (4-26) transforms to Eq. (3-62) for an isotropic system.  4.2.2 Special Case of Solid Density Equation for Isotropic Materials The special case of Eq. (4-17) for isotropic material is examined here. For an isotropic solid phase, we can write: ?̇?𝑆𝛅T𝐃𝑆−1𝛅 =1𝐾𝑆?̇?𝑆 (4-33) as well as: 𝛅𝐓𝐃𝑆−1?̇?′ =13𝐾𝑆𝛅𝐓?̇?′ (4-34) Substituting Eqs. (4-33) and (4-34) into Eq. (4-17) yields: 1𝜌𝑆?̇?𝑆 =1𝐾𝑆?̇?𝑆 −13(1 − 𝜑)𝐾𝑆𝛅𝐓?̇?′ − 𝛅𝐓?̇?𝑆𝑓 (4-35) which, in index notation, is identical to the equation introduced by Lewis and Schrefler [69] (see Eq.(3-12)).  83  4.3 First Invariant of 𝛔′ The first invariant of the effective stress tensor 𝛔′ was introduced by Lewis and Schrefler [69] (see Eq. (3-13)) for of the isotropic model). Here, this quantity is derived for the more general case of orthotropic constituents and subsequently checked to verify that the special case of the equations for isotropic system can be recovered.  Considering Eq. (4-3) and (4-6), the effective stress 𝛔′ can be written as: 𝛔′ = 𝐃𝑆𝐾𝛆 + 𝑃𝑆(𝛅 − 𝐛) (4-36) where its first invariant in terms of the displacement field is given by 𝛅𝐓𝛔′ = 𝛅𝐓𝐃𝑆𝐾𝛁T?̇? + 𝑃𝑆(3 − 𝛅𝐓𝐛) (4-37) Using the orthotropic Biot’s coefficients (Eq. (4-8)) we arrive at: 𝛅𝐓𝛔′ = 𝛅𝐓𝐃𝑆𝐾𝛁T?̇? + 𝑃𝑆𝛅𝐓𝐃𝑆𝐾𝐃𝑆−1𝛅 (4-38)  4.3.1 Special Case for Isotropic Materials When the solid and solid-skeleton possess isotropic properties, we can write: 𝛅𝐓𝐃𝑆𝐾𝛁T?̇? = 3𝐾𝑆𝑘𝛅𝐓𝛁T?̇? (4-39) 𝛅𝐓𝐃𝑆𝐾𝐃𝑆−1𝛅 =𝐾𝑆𝑘𝐾𝑆 (4-40) Using Eqs. (4-39) and (4-40) in (4-38) results in: 𝛅T𝛔′ = 3𝐾𝑆𝑘(𝛅𝐓𝛁T?̇? +1𝐾𝑆𝑃𝑆) (4-41) This equation was introduced by Lewis and Schrefler [69] for the case of isotropic materials and written in index notation as Eq. (3-13) in Chapter 3. 84   4.4 Evolution of Material Properties  For evolution of the mechanical properties of the material in the IFS-Isotropic model (shown in general form in Eq. (3-17)), only two independent elastic constants of the solid and solid-skeleton, say the bulk and shear moduli, are evolving as the resin material cures. For example, 𝐾𝑆𝐾 transforms from the fibre-bed bulk modulus (𝐾𝑓𝑏) for the unsolidified fluid (𝜆 = 0) to the bulk modulus of the fully cured composite material (𝐾𝑐) for the fully solid resin (𝜆 = 1). For an orthotropic material system, the stiffness matrix (shown in Eq. (3-31) for isotropic materials) is defined as follows: 𝐃 =[           1 − 𝜈23𝜈32𝐸2𝐸3Ξ𝜈21 + 𝜈23𝜈31𝐸2𝐸3Ξ𝜈31 + 𝜈21𝜈32𝐸2𝐸3Ξ 1 − 𝜈13𝜈31𝐸1𝐸3Ξ𝜈32 + 𝜈31𝜈12𝐸1𝐸3Ξ  1 − 𝜈12𝜈21𝐸1𝐸2Ξ0 0 00 0 00 0 0    𝑆𝑦𝑚𝑚.     𝐺23 0 0 𝐺13 0  𝐺12]            (4-42) where  Ξ =1 − 𝜈12𝜈21 − 𝜈23𝜈32 − 𝜈13𝜈31 − 2𝜈21𝜈32𝜈13𝐸1𝐸2𝐸3 (4-43) and the following reciprocity relation applies. 𝜈𝑖𝑗𝐸𝑗 = 𝜈𝑗𝑖𝐸𝑖 (4-44) where repeated indices do not imply summation.   85  For orthotropic solid and solid-skeleton, all the quantities in the stiffness matrix of the solid and solid-skeleton are evolving. In other words, each component of the stiffness matrix changes from the corresponding component in the fibre/fibre-bed stiffness matrix to the corresponding component in the composite stiffness matrix during the course of curing. These evolving properties result in evolution of Biot’s coefficient vector (for the orthotropic system) from the following form for the unsolidified fluid phase at one extreme 𝐛 = 𝛅 − 𝐃𝑓𝑏𝐃𝑓−1𝛅 (4-45) to a 6 × 1 vector of zero values for the fully solidified system at the other extreme.  4.5 Summary of Governing Equations Incorporating the orthotropic from of the quantities presented in the previous sections of this chapter in the governing equation of the integrated model presented in Chapter 3, the governing equations of the two- and three-phase IFS-Orthotropic model can be derived.   4.5.1 Two-Phase Model The equation involving the mass conservation, equilibrium of the matrix, and equilibrium of the system for the 2IFS-Orthotropic model can be presented in vector form as follows:  Mass conservation: (𝛁𝐛)T?̇? + (𝛁𝛅)T𝐯𝐹 + ?̇?𝐹 (𝜑𝐹𝐾𝐹+ 𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝐹𝛅)) − 𝜑𝐹(1 − 𝜆)𝛅𝐓?̇?𝐹𝑓− (1 − 𝜑𝐹)(1 − 𝜆)𝛅𝐓?̇?𝑆𝑓 − 𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾𝛆𝑆𝐾𝑓 = 0 (4-46)  Darcy’s flow for the fluid: 86  (𝛁(𝑃𝐹𝐈))𝛅 + 𝜇𝐹𝐒−1 𝐯𝐹 = 𝟎 (4-47)  Equilibrium equation: 𝛁𝛔𝑆𝐾 − (𝛁(𝑃𝐹𝐈))𝐛 = 𝟎 (4-48) where 𝐛 = 𝛅 − 𝐃𝑆𝐾𝐃𝑆−1𝛅 and the rate of stress in the solid-skeleton is given by: ?̇?𝑆𝐾 = 𝐃𝑆𝐾?̇?𝑆𝐾𝜎 + 𝑃𝐹?̇? (4-49)  The iterative solution for the orthotropic model is the similar to that used for the isotropic model (see Eq. (3-40)). Some of the matrix quantities are different from those listed in Table 3.1. These slightly different components are presented in Table 4.1.   4.5.2 Three-Phase Model The governing equations for the 3IFS-Orthotropic model are:  Mass conservation: (𝛁𝐛)T?̇? + (𝛁𝛅)T𝐯𝐹 + (𝛁𝛅)T𝐯?̂? + ?̇?𝐹𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)𝑃𝐹 + ?̇??̂?𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)𝑃?̂?+ ?̇?𝐹 (𝜑𝐹𝐾𝐹+ 𝜓𝐹𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅)) + ?̇??̂? (𝜑?̂?𝐾?̂?+ 𝜓?̂?𝛅𝐓𝐃𝑆−1(𝐛 − 𝜑𝛅))− 𝜑𝐹(1 − 𝜆)𝛅𝐓?̇?𝐹𝑓 − 𝜑?̂?𝛅𝐓?̇??̂?𝑓 − (1 − 𝜑𝐹)(1 − 𝜆)𝛅𝐓?̇?𝑆𝑓− 𝛅𝐓𝐃𝑆−1𝐃𝑆𝐾𝛆𝑆𝐾𝑓 = 0 (4-50)  Darcy’s flow for fluid 𝐹: (𝛁(𝑃𝐹𝐈))𝛅 + 𝜇𝐹𝐒𝐹−1 𝐯𝐹 = 𝟎 (4-51)  Darcy’s flow for fluid ?̂?: (𝛁(𝑃?̂?𝐈))𝛅 + 𝜇?̂?𝐒?̂?−1 𝐯?̂? = 𝟎 (4-52) 87   Pressure equation: 𝑃?̂? − 𝑃𝐹 = 𝑃𝑐 (4-53)  Equilibrium equation: 𝛁𝛔𝑆𝐾 − (𝛁(𝑃𝐹𝐈) − 𝜓𝐹(𝛁(𝑃?̂?𝐈) − 𝛁(𝑃𝐹𝐈)))𝐛 = 𝟎 (4-54) where: ?̇?𝑆𝐾 = 𝐃𝑆𝐾?̇?𝑆𝐾𝜎 + (𝑃?̂? − 𝜓𝐹(𝑃?̂? − 𝑃𝐹))?̇? (4-55)  The iterative solution scheme for the discretized form of these equations takes on a similar form to that used for the 3IFS-Isotropic model (Eq. (3-73) and Table 3.2), with notable differences in some of the matrices and vectors that are listed in Table 4.2.  4.6 Summary  The more general form of the IFS model is presented in this chapter for both the two- and three-phase composite systems with orthotropic solid (fibre) and solid-skeleton (fibre-bed). First, the orthotropic form of the quantities and equations in the IFS model are derived where their special case of isotropic form is also examined and compared to those that are presented in the previous chapter. In the next chapter, numerical case studies are presented to show performance and effectiveness of the both isotropic (Chapter 3) and orthotropic (Chapter 4) IFS models for both two- and three-phase composite systems.   88  4.7 Tables 𝐊𝒖𝑷 = −∫𝐁𝑇 Ω𝐛𝐍𝑃𝑑Ω   𝑪𝑷𝒖 = ∫𝐍𝑃𝑇 Ω𝐛𝑇𝐁𝑑Ω 𝑪𝑷𝑷 = ∫𝐍𝑃𝑇 Ω(𝜑𝐾𝐹+ 𝛅𝑇𝐃𝑆−1(𝐛 − 𝜑𝛅))𝐍𝑃𝑑Ω  Table 4.1 Definition of specific matrices and vectors used in the two-phase orthotropic FE model that are different from those defined for the two-phase isotropic model (listed in Table 3.1)   𝐊𝒖𝑷Φ = −𝜓Φ ∫𝐁T Ω𝐛𝐍𝑃𝑑Ω 𝐊𝑷?̂?𝑷Φ = ∆𝜓Φ ∫𝐍𝑃𝑇 Ω𝛅𝑇𝐃𝑆−1(𝐛 − 𝜑𝛅)𝐍𝑃𝑑Ω 𝐂𝑷?̂?𝒖 = ∫𝐍𝑃𝑇 Ω𝐛𝑇𝐁𝑑Ω 𝐂𝑷?̂?𝑷Φ = ∫𝐍𝑃𝑇 Ω(𝜑Φ𝐾Φ+ 𝜓Φ𝛅𝑇𝐃𝑆−1(𝐛 − 𝜑𝛅))𝐍𝑃𝑑Ω 𝐟𝒖𝑏𝑷 = ∫𝐁𝐓 Ω𝐛𝐍𝑷𝑑Ω   Table 4.2 Definition of specific matrices and vectors used in the three-phase orthotropic FE model that are different from those defined for the two-phase isotropic model (listed in Table 3.2)89  Chapter 5: Numerical Applications The FE matrix equations have been coded into a MATLAB® software developed in-house. Numerical case studies are presented in this chapter to show the accuracy and effectiveness of the developed IFS model for composite materials under various loading and boundary conditions experienced at different stages of processing, and including constituents with diverse material properties. Initial examples in each section of this chapter serve to verify the accuracy of the model and its implementation. For this purpose, examples are presented at both unsolidified and fully solidified extremes of processing. In some curing examples where the resin material solidifies, the significance of considering flow of resin throughout the process by adopting the IFS method is demonstrated through comparison of predicted results with those obtained from the stress model alone developed by Zobeiry [50].   5.1 Assumptions In all examples presented, 2D plane strain conditions are assumed to prevail. In keeping with the plane strain assumptions, the change in in-plane strains due to suppression of all out-of-plane strains is taken into account [8]. Also, a uniform temperature over the domain of the problems is assumed in all examples in order to avoid running the heat transfer analysis. All the simulations are carried out quasi-statically as dynamic (inertia) effects are ignored in the current model.  In all examples, the composite material reinforcement is AS4 carbon fibre. The matrix material  is either 3501-6 epoxy resin, air, or a mixture of both which is the case in the three-phase model. Table 5.1 shows the material properties for each constituent. For the AS4 fibre, direction 1 refers to the fibre direction and directions 2 and 3 are transverse to the fibres. The thermal properties of 90  the fibre-bed are assumed to be identical to those of the fibres. Applying the ideal gas assumptions to air, its bulk modulus is taken to be equal to its absolute pressure. Also, a very low value (almost zero) is assumed for the shear modulus of air. The cure kinetics and viscosity models for the 3501-6 resin are adopted from Lee et al. [75] and Hubert et al. [15] and presented in Table 5.2. For the fibre-bed, Table 5.3 presents the permeability and elasticity models that are based on the formulation by Cai and Gutaowski [70] with the corresponding model parameters adopted from Hubert [1]. Also, the time-dependent mechanical properties of the resin are assumed to follow the pseudo-viscoelastic (PVE) formulation [51] summarized in Table 5.4 with relaxation time and weight factors taken from [76] and shown in Table 5.5.   In those examples where the resin cures and solidifies, the solidification factor is assumed to be a piece-wise linear (tri-linear) function of the degree of cure as follows: 𝜆 = 0          for 𝜒 < 𝜒𝑎 𝜆 =𝜒−𝜒𝑎𝜒𝑏−𝜒𝑎            for 𝜒𝑎 ≤ 𝜒 ≤ 𝜒𝑏 𝜆 = 1         for 𝜒𝑏 < 𝜒 (5-1) Figure 5.1 shows the demarcation points 𝜒𝑎 and 𝜒𝑏 that are chosen to be 0.7 and 0.9, respectively. Choosing these values, the solidification factor starts to rise in synchronization with evolution of resin properties, i.e. 𝐾𝑟 and 𝐺𝑟, predicted based on the PVE formulation. The solid and solid-skeleton properties are also assumed to be linearly dependent on the solidification factor evolving from fibre and fibre-bed properties, respectively, to those of the composite as follows: 𝐃𝑆 = 𝐃𝑓 + (𝐃𝑐 − 𝐃𝑓)𝜆 (5-2) 91  𝐃𝑆𝐾 = 𝐃𝑓𝑏 + (𝐃𝑐 − 𝐃𝑓𝑏)𝜆  𝛼𝑆𝐾𝑡ℎ = 𝛼𝑓𝑏𝑡ℎ + (𝛼𝑐𝑡ℎ − 𝛼𝑓𝑏𝑡ℎ)𝜆 𝛼𝑆𝐾𝑐𝑠 = 𝛼𝑐𝑐𝑠𝜆 which is shown in Figure 5.2 for bulk and shear moduli as well as thermal expansion and cure shrinkage coefficients.  For the two-phase system, the mechanical properties of the composite material, i.e. 𝛼𝑐𝑡ℎ, 𝛼𝑐𝑐𝑠, and components of 𝐃𝑐, are calculated during the simulation using a micromechanics model adopted by Johnston et al. [49] and Bogetti and Gillespie [4] using resin and fibres properties (shown schematically in Figure 5.3(a)). For the three-phase system, these properties are computed using a three-phase micromechanics model outlined below and shown schematically in Figure 5.3(b). First, mechanical properties of the matrix, which consists of air and resin, are calculated using the Generalized Self Consistent (GSC) micromechanics model by Christensen [77] while the micromechanics model by Levin [78,79] is adopted to compute the effective coefficients of thermal expansion (𝛼𝑡ℎ) and cure shrinkage (𝛼𝑐𝑠). The effective properties of the two-phase matrix are then combined with the fibre properties through the same micromechanics model used for the two-phase system developed by Bogetti and Gillespie [4].   In the forthcoming sections, first, an example is presented that shows the applicability of the IFS model to fully solid phases (the condition at the end of the process). Then examples are categorized into two- and three-phase models where both isotropic and orthotropic cases are discussed for fully unsolidified fluid phase(s) as well as solidifying resin material as the 92  composite system undergoes curing. In the two-phase examples, the importance of using the current IFS method for modelling the flow of a highly compressible fluid phase is highlighted. Also, in all the isotropic examples considered, the composite material is assumed to be a unidirectional laminate with fibre orientations perpendicular to the domain of interest.  5.2 Fully Solidified Material Example To investigate the performance of the IFS model for the fully solidified material (𝜆 = 1), a single element with isotropic properties is subjected to a combination of mechanical and thermal loadings as shown in Figure 5.4. To highlight the effect of nonlinear geometry due to large deformations, extremely high values of applied loadings are used thus inducing very large strains in the element. A sensitivity analysis is performed to obtain the required size of the time step for convergence of the result. Based on this sensitivity study, a small time step, ∆𝑡 = 0.1 s, is used. Also, both mechanical and thermal loadings are applied as ramp functions through the 40 seconds of the simulation time. All other required material properties are shown in Figure 5.4 that apply to either two- or three-phase fully solidified material. The predicted time-history of developed strains and stresses in the element are presented in Figure 5.5 and Figure 5.6, respectively. For comparison, the equivalent results obtained using the 4-noded bilinear plane strain quadrilateral element (CPE4) in Abaqus are also shown in the figure. The excellent agreement between the results from the IFS model and the Abaqus stress analysis serves as a verification of the accuracy of the code (model) for the fully solidified material.  It should be noted that the restriction that applies to the time step size of the IFS model (because of infinitesimal strain assumption) is not applied to the Abaqus model due to utilization of the 93  updated Lagrangian method (the flag NLGEOM was activated) which can accommodate significantly larger time steps. Fully solid isotropic example is considered and shown in Figure 5.7 to show the effect of time step size in the IFS model where the load is applied gradually over 460 s. It should be noted that a very large value of the force has been used intentionally to exaggerate and highlight the effect of nonlinear geometry. The load-displacement curve for the top surface is shown in Figure 5.8 to highlight the effect of nonlinear geometry. It can be seen that Abaqus yields the same results independent of the size of the time step as it uses the updated Lagrangian method. However, in our current IFS model, the time step has to be chosen to be small enough to capture the effect of nonlinear geometry using the infinitesimal strain assumptions made in the formulation.   5.3 Two-Phase Examples 5.3.1 Unsolidified Isotropic Material - Mechanical Loading As mentioned before, incorporating the capability to handle highly compressible materials in process modelling of composite materials is one of the primary impetus for developing the current IFS model. In this example, a single element of unsolidified composite material (𝜆 = 0) is tested under a uniform pressure at a constant temperature as shown in Figure 5.9. An applied pressure, 𝑃𝑏 = 0, is imposed at the top permeable surface. Two distinct materials are considered for the fluid phase of the composite material: (i) 3501-6 resin, which is nearly incompressible, and (ii) atmospheric air, which is highly compressible. In this example, all resin properties correspond to its relaxed conditions and therefore, for the purposes of simulations, the temperature is assumed to be above the resin’s glass transition temperature, 𝑇𝑔 [76]. The air’s bulk modulus is assumed to be the sum of its pressure and the atmospheric pressure. Therefore, 94  all pressures in the model are relative to the atmosphere, e.g. the pressure imposed at the top surface is the atmospheric pressure. Using a constant temperature of 20℃ along with a zero degree of cure for the resin (𝜒 =0), results in a constant viscosity of  1.177 × 104 Pa.s for the resin and 1.813 × 10−5 Pa.s for air.  The bulk modulus of the fibre-bed (assumed to be isotropic in this example) is obtained using the expression, 𝐾 = 𝐷22 − 4𝐺23 3⁄  where 𝐷22 is the transverse stiffness modulus (second diagonal component of the stiffness tensor, see Table 5.3). To compare the results obtained for resin and air as selected fluid material for the composite system, the same distributed load, 𝑞 = 50 kPa, is applied instantaneously to both cases and maintained for a duration of 4000 min for the resin (using a time step of  ∆𝑡 = 30 s) and 0.002 min for air (using a time step of ∆𝑡 = 0.001 s). The numerical results for the time-histories of the developed fluid pore pressure, fluid velocity at the top surface, strain, and fluid volume fraction are presented in Figure 5.10 to Figure 5.13 for both cases of the resin and air as separate fluid phases of the sample.  It can be seen from Figure 5.10 that the resin pressure remains constant and equal to the applied mechanical load for the impermeable boundary condition. For permeable top surface, it drops slightly after a very long time. On the other hand, pressure of air drops to zero in less than a second for permeable top surface as a considerable amount of air evacuates quickly and the load is carried entirely by the fibre-bed (note that inertia effects are ignored in the model). For the impermeable boundary condition, the compressed air takes a major share of the applied load which is seen in its constant pressure (~20 kPa) while the rest of the load (~30 kPa) is carried by the fibre-bed.   95  For the permeable case, the resin velocity shown in Figure 5.11 is very low due to its high viscosity compared to air velocity. The resin velocity decays with time because of the decrease in resin pressure seen in the previous figure. Same as pressure, the air velocity drops to zero very quickly for permeable top surface as a result of the drop in air pressure due to the load being taken by the fibre-bed.  In Figure 5.12, again the difference seen in the total vertical strains for the two cases (resin and air) under both permeable and impermeable conditions show the effect of high compressibility of air compared to resin. For example, for the impermeable boundary condition, air compresses instantaneously after the load application and a significant strain is obtained while there is almost no strain for the case of resin as it is nearly incompressible.  The evolution of fluid volume fraction is shown in Figure 5.13 where it can be seen that even for the impermeable condition a significant change in volume fraction of air is obtained due to its high compressibility. This is a direct result of considering compressibility in the volume fraction calculations using Eqs. (3-50) and (3-51) which could not be predicted if only the flow of fluid was taken into account in the calculation of volume fraction (Eq. (3-45)).  It is noted that the constant total compressive stress (−50 kPa) in the vertical direction, which is equal to the applied load, was observed throughout the simulations for all cases considered which further serves as a verification of the accuracy of the IFS model and its FE implementation in the MATLAB® code.    96  5.3.2 Curing Isotropic Composite Material In this example, an isotropic composite system containing resin as its matrix phase undergoes curing under different boundary conditions imposed on the fluid velocity and system displacement as shown in Figure 5.14. The boundary conditions considered for the presented single element test case represent the multitude of conditions that can possibly prevail at different locations (local elements of the FE mesh) within more complex and practical geometries such as that shown in Figure 5.15. For instance, the top left boundary condition in Figure 5.14 represents the conditions due to placement of a bleeder layer on the surface of the part being manufactured while the displacements are restricted by a dam (element in the red circled zone #1 in Figure 5.15). The top right boundary condition in Figure 5.14 is representative of the conditions that apply to the inner elements of a large practical example where both inter-element displacement and flow are permitted (element in the circled zone #2 in Figure 5.15). The bottom left condition in Figure 5.14 corresponds to that of an element adjacent to a mould where both flow and displacement are restricted (element in the circled zone #3 in Figure 5.15). The bottom right condition in Figure 5.14 is representative of the case when an impermeable bag, rather than a bleeder, is used on the surface of the part (element in the circled zone #4 in Figure 5.15).   The applied temperature cycle is shown in Figure 5.16 along with the predicted time history of viscosity (𝜇), degree of cure (𝜒) (using equations in Table 5.2), and solidification factor (𝜆) (using Eq. (5-1)). As the resin is less compressible and has a higher viscosity than air (considered as the matrix in the next example), the strains that would be induced in the sample are relatively 97  small and therefore convergence (accuracy) of the numerical results can be achieved with a relatively coarse time step, ∆𝑡 = 30 s. Same sample with some of the boundary conditions shown in Figure 5.14 was also investigated by Haghshenas [17] using the previous IFS model. For isotropic composite system with nearly incompressible phases (which is the case in this example), the new IFS model developed in this work predicts nearly identical results as the previous IFS model [17].   Resin pressure and velocity are only shown for the first 60 min of the cure cycle in Figure 5.17 and Figure 5.18 when the solidification factor 𝜆 < 1 and there is a fluid phase in the system with corresponding pressure and velocity. Furthermore, as shown in Figure 5.19 for resin volume fraction, the stress model has no means of capturing the continuous changes in volume fraction of the resin for all boundary conditions considered. This is because the stress model does not account for flow.  The most increase in pressure is seen for the sample with impermeable and constrained conditions at the top surface in Figure 5.17. For this case, the increase of pressure due to expansion of resin cannot be relieved neither from its flow out of the sample nor expansion of the sample. In the sample with permeable and constrained conditions, an increase in pressure is seen as the resin is initially extremely viscous and its flow through the top surface is very limited. Shortly afterwards, pressure drops as viscosity decreases and the resin can flow through the top permeable surface and relieve the induced pressure. For samples with free condition at 98  the top surface, no significant pressure would be induced in the resin as the samples freely expand due to applied temperature.  Velocity of resin through the top surface is shown for both samples with permeable top surface in Figure 5.18. For the constrained sample, the resin velocity trend is straightforward and agrees with its pressure-time history discussed in the previous paragraph. As the movement of the top surface is constrained, any volumetric change in the system leads to an instantaneous flow of resin through the top surface. For the sample with free top surface, the condition is more complex as the top surface displaces and the presented resin velocities are relative to this moving surface. For this condition, the resin expands initially while its viscosity is still quite high. Because of composite systems constituents expansion, the top surface moves upward (sample expands). Consequently, resin velocity is negligible initially. After a while, resin viscosity drops allowing it to readily flow out of the top permeable surface without imposing any strain on the sample. At this stage, the top surface moves downward (as will be discussed subsequently) and the induced strain in the sample recovers. Resin velocity at this stage is relatively high as its viscosity is relatively low and the top surface moves downward. The combined effect of these two phenomena result in a higher relative velocity of the resin (relative to the top surface). The negative flow seen for both conditions just before the velocity goes to zero (because of the increase in resin viscosity) is an important phenomenon that is well captured by our IFS model. This negative velocity occurs as cure shrinkage starts just before the viscosity increases and material fully solidifies causing inward flow of resin into the sample. This is a real phenomenon which is a very important source of porosity in composite materials. It is worth noting that 99  although the magnitude of velocities are very small, their general trend under different conditions considered here should be evaluated in detail as was done in this example. In this way, the capability of the IFS model is assessed allowing it to be used with confidence in the future for cases where the velocities are higher, e.g. larger composite structures, low viscosity fluids (such as what will be observed for air in subsequent sections), or in other types of processing such as infusion where resin velocity is significantly higher.  The stresses predicted by the two-phase IFS model and the stress development model (Zobeiry et al. [51]) for two selected permeable and impermeable boundary conditions but fully constrained displacement are shown in Figure 5.20. As expected, the stresses obtained from both the IFS and stress models are in perfect agreement for the impermeable boundary conditions (red lines) when flow is inhibited. For permeable but fully constrained displacement boundary condition (blue lines), at the early stages of the process, the resin flows out of the system thus relieving some of the developed pressure as it undergoes thermal expansion due to temperature rise. While the resin is still in its pre-gelation state and flows out, no stresses develop in the system. This behaviour is well captured by the IFS model. Stresses only start to develop after gelation of the resin. This observation explains the notable difference in the final predictions of residual stresses in the system between the IFS model and the stress model for permeable top surface. However, the stress model yields the same results for both permeable and impermeable cases which further highlights the errors involved when flow is ignored in the process simulation. Stress evolution due to resin cure shrinkage as well as cool-down at the end of processing are also evident from these curves. 100   The predicted time history of the vertical strain for the case when the top surface is free is shown in Figure 5.21. Again the stress model erroneously predicts the same strains for both the permeable and impermeable top surface conditions. In contrast, the IFS model predicts a completely different time history of strain for the case of permeable top surface. For impermeable and free conditions at the top surface more strain is induced in the sample as resin cannot leave the system thus imposing tensile strain in the composite system due to its expansion. For permeable and free condition, initially, the same behavior is observed as resin cannot flow through the top surface due to its high viscosity even though this surface is permeable. After a while, the viscosity decreases and resin starts to flow out of the top permeable surface. Consequently, the top surface experiences less drag force and recovers some of the strain that was previously induced in the sample. For both samples, cure shrinkage effect is observed after about 50 min followed by thermal contraction during the cool-down stage of the temperature cycle.  5.3.3 Composite with Highly Compressible Fluid Phase under Thermal Loading In order to check the performance of the IFS model for highly compressible materials under thermal loading, the same temperature cycle as that used in the previous example is applied to an isotropic composite system containing air as its fluid phase instead of resin. The geometry and boundary conditions are the same as in the previous example (Figure 5.14). The temperature dependent air viscosity, which is obtained from the model presented in Table 5.2, is shown in Figure 5.22. Since air does not solidify, 𝜆 = 0 throughout the process. Dealing with a more compressible material and consequently higher values for strain, it was found that a smaller time 101  step, ∆𝑡 = 15 s, compared to that used in the previous example was required to achieve convergence. The time histories of air pressure, air velocity through the top permeable surface, air volume fraction, stresses, and vertical strain are shown in Figure 5.23 to Figure 5.27. The predicted lower pressures and stresses combined with higher velocities and strains compared to the previous example where the fluid phase was a pre-gelled resin demonstrates clearly the effect that the high compressibility of air and its low viscosity have on the system response which is well captured by the current IFS model. For both impermeable samples, an increase in pressure is observed in Figure 5.23 as a result of the temperature ramp. For permeable samples, however, no pressure is induced as air readily flows out and relieves its pressure. The air viscosity is only dependent on the temperature (Table 5.1 and Figure 5.22) and there is no cure shrinkage involved when considering air. Therefore, all increased pressure will recover after the cool-down. For velocity shown in Figure 5.24, air with its low viscosity readily flows through the top permeable surface after expansion due to applied temperature. As temperature increases (with constant rate) viscosity of air increases proportionately (as shown in Figure 5.22). Consequently, the flow velocity decreases. During the temperature hold, no expansion and no flow through the top surface is expected as seen in the figure. The opposite inward flow is seen during cool-down as a result of thermal contraction of air. It is noted that due to very low viscosity of air, no significant strain is expected in the sample with free top surface (low drag force on the top surface during flow of air). As a result, the top surface is not moving (as opposed to the cases with resin matrix) and consequently no considerable difference is seen for air velocity for constrained and free samples.  102  The appreciable changes in the volume fraction of air (Figure 5.25) for the case involving impermeable and free boundary conditions is a direct consequence of using volume fraction calculations that explicitly account for compressibility of phases (Eq. (3-50) or (3-51)) and clearly shows its significance for highly compressible materials. Comparison of the volume fraction variation for this example with the previous example where a nearly incompressible resin was considered as the fluid phase (Figure 5.19) also highlights the accurate performance of the model in capturing the behaviour of composite materials with a highly compressible fluid phase. As expected, the most change in volume fraction is seen for the case with impermeable and free conditions. The highly compressible air expands considerably while it cannot leave the system. This induces a significant amount of expansion in the sample while the volume of fibers do not change significantly compared to air. Therefore, the air volume fraction in the sample increases. For permeable cases, as noted before, the air flows out without imposing a considerable expansion in the sample. Therefore, the total volume of both phases (air and fibre) remains almost constant. It is interesting that only a very minor decrease in the air volume fraction is observed in the sample with impermeable and constrained conditions. Temperature increase leads to a very minor increase in the volume of the fibres. Although air wants to expand due to the same temperature increase, because of its high compressibility it cannot apply enough load on the stiff fibres. The result is a minor expansion of the fibres and a minor contraction of the air resulting in a minor decrease in air volume fraction (note that the total volume of the sample is constant in this case).   No significant stress is expected for samples with free top surface (stress is the same as applied mechanical load which is zero in this example). For impermeable and constrained condition 103  shown in Figure 5.26, as the flow of air out of the sample as well as the expansion of the sample (due to considerable expansion of air) as a whole is restricted, a significant compressive stress is induced. On the other hand, the stress induced in the sample with permeable condition is minimal as air flows out of the top surface and the minor compressive stress is only induced because of the minor expansion of the fibres in both directions. The significant difference between the induced strains in impermeable and permeable cases shown in Figure 5.27 can easily be explained based on the same arguments presented in this section.  5.3.4 Curing Orthotropic Composite Material  In this example, the implementation of the IFS-Orthotropic model and anisotropic features of the code is evaluated. A composite sample undergoes curing as shown in Figure 5.28 where both the top and right surfaces of the sample are assumed to be free and permeable. Fibres are in-plane and three different values of 𝛽 =0°, 45°, and 90° for their orientation in the 𝑥 − 𝑧 plane of interest are considered. The sample undergoes the same temperature cycle as that in the previous examples (see Figure 5.16).  The predicted deformed shape of the sample by both the IFS and stress model are shown in Figure 5.29 for all three orientations (note that the deformations are scaled by a factor of 10 for ease of illustration). Similarity of the deformed shapes for 𝛽 = 0° and 90° cases, confirms the validity of the coordinate transformations implemented in the orthotropic IFS model. The expected shear deformation for the case where 𝛽 = 45° can also be observed. As it is seen, the stress model erroneously predicts an expansion in the transverse direction to the fibres for 𝛽 =0° and 90° which stems from ignoring flow of the resin in that formulation.  104   The time-history of the resin velocities through the top and right surfaces of the sample are presented in Figure 5.30. Note the consistency in the predicted resin velocity for 𝛽 = 0° and 90°, i. e. 𝑣𝑥|0° = 𝑣𝑧|90°. The difference in velocity between the longitudinal and transverse directions to the fibre represents the mismatch between the effect of permeabilities along and perpendicular to the fibre.  The evolution of normal strain in the transverse direction with respect to the fibres, 𝜀𝑥|𝛽=90° =𝜀𝑧|𝛽=0°, as well as shear strain for 𝛽 = 45°, 𝛾𝑥𝑧|𝛽=45°, are presented in Figure 5.31. At the early stages of the process, the resin undergoes thermal expansion due to temperature rise which leads to the expansion of the sample in both transverse and longitudinal directions with respect to the fibres. Because of the very high stiffness properties of the composite sample in the fibre direction, the strain in the transverse direction (i.e. 𝜀𝑥|𝛽=90° = 𝜀𝑧|𝛽=0°) is much greater than the strain in the longitudinal direction (i.e. 𝜀𝑥|𝛽=0° = 𝜀𝑧|𝛽=90°), which are almost zero and not shown in the figure. Furthermore, permeability in the transverse direction is lower than that in the longitudinal direction and the highly viscous resin experiences more resistance while flowing in the transverse direction to the fibres (through the top surface for the case of 𝛽 = 0° and through the right surface for the case of 𝛽 = 90°) which leads to higher strain in the transverse directions (𝜀𝑥|𝛽=90° = 𝜀𝑧|𝛽=0°) compared to the longitudinal direction (i.e. 𝜀𝑥|𝛽=0° = 𝜀𝑧|𝛽=90°). Consequently, it induces more strain in the sample at this stage. As a result of high viscosity of the resin, the strain predictions of the IFS model are close to those of the stress model initially. But afterwards the resin viscosity drops and it can readily flow out without imposing any tensile 105  strain on the composite system. Therefore, the strain in the sample drops. This behaviour cannot be captured by the stress model. After gelation of the resin, the transverse strain starts to develop once again due to cure shrinkage with a trend very similar to that trend from the stress model. The initial difference between the predicted strains from the IFS and stress models leads to a noticeable difference between the predicted final strains both in the magnitude and sign of strain. This example illustrates the significance of using the IFS model for process simulation of composite materials.   The stress prediction for 𝛽 = 45° is shown in Figure 5.32 where again a considerable difference between the results from the two models is apparent. It is noted that as expected, due to the absence of mechanical loading on the system, all other shear and normal stresses for various cases, i.e. 𝜎𝑥|𝛽=0°, 𝜎𝑧|𝛽=0°, 𝜎𝑥|𝛽=90°, 𝜎𝑧|𝛽=90°, 𝜏𝑥𝑧|𝛽=0°, and 𝜏𝑥𝑧|𝛽=90°, are negligible.  5.3.5 Unidirectional [0°] Angle Laminate Undergoing Compaction and Cure In this example, we consider a unidirectional [0°] laminate conformed and fully bonded to an L-shaped convex tool as shown in Figure 5.33. Only half of the laminate is modeled due to symmetry and three key points A, B, and C corresponding to the corner, transition between the curved and flat part, and middle of the flat part, respectively, are labeled. The top and side surfaces are permeable and a 540 kPa distributed load is applied on the top surface. The laminate undergoes the temperature-time history shown in Figure 5.34 where the resulting viscosity, resin degree of cure, and solidification factor are also presented. The evolution of the resin bulk and shear moduli predicted from the pseudo-viscoelastic model [51] are also presented in Figure 5.35. Hubert [1] performed a sensitivity analysis and showed the effect of the fibre-bed 106  shear modulus, i.e. 𝐺𝑓𝑏12 = 𝐺𝑓𝑏13, on the final shape of the laminate. Comparing the numerical simulation of the deformed shape of the laminates with different fibre-bed shear moduli to the experimentally measured deformed shape, he chose 𝐺𝑓𝑏12 = 𝐺𝑓𝑏13 = 500 kPa which is used in the current example. The problem is modelled with different levels of discretization as shown in Figure 5.36 to ensure convergence of the results. Because of higher gradient at the corner of the laminate, a finer mesh is used there while the mesh gradually becomes coarser in the flat part of the laminate. The 36 × 60 mesh, which guarantees the convergence of results and at the same time yields a practical solution time, is chosen to report the results.   The final shape of the laminate is shown in Figure 5.37. It can be seen that the corner thickening as well as the expected lateral deformation of the edge for a [0°] laminate processed on a convex tool are qualitatively captured by the IFS model. These two observed patterns in deformation, are in agreement with the previous experimental and numerical investigations [1,15]. Hubert and Poursartip [1,80] tested different laminated composite materials on both convex and concave tools. For [0°] laminate of the same composite material (AS4-3501-6) on a convex tool, for permeable conditions, they observed corner thickening after the curing process as shown in Figure 5.38 . The predicted deformed shape of the laminate at the end of compaction using their flow model (Figure 5.39) also shows same lateral deformation at the edge as observed from the current IFS model predictions. The same corner thickening effect and edge lateral deformation seen after simulating the sample using the current IFS model further confirms the capability of the IFS model for processing of composite materials. The normal displacement of key points A, B, and C as well as the deformation of the free edge as a function of time are presented in 107  Figure 5.40 using both IFS and stress models. It can be seen from this figure that the full compaction is achieved in almost 60 minutes after which no further deformation (compaction) of the laminate takes place. The full compaction corresponds to the condition where all the applied load is carried by the fibre-bed and little pressure is carried by the resin as a result of having permeable boundary conditions. Also, the incorrect tensile deformation of all the points (thickening of the laminate everywhere compared to its initial thickness) predicted by the stress model is evident from the figure.  Capturing continuous change of resin/fibre volume fraction distribution is one of the notable advantages of the IFS model over the stress model as the stress model has no means of accommodating for changes in the volume fraction calculation of constituents (phases). The predicted distribution of fibre volume fraction, (1 − 𝜑), at different instances of time is shown in Figure 5.41. At the end of the process, a uniform compaction in the flat part of the laminate is evident while there is less compaction at the corner, associated with lower fibre volume fraction (a higher resin volume fraction). This is in qualitative agreement with the expected compaction behaviour of [0°] laminate on a convex tool [15].   Figure 5.42 presents the distribution of resultant axial force and bending moment along the mid-plane of the laminate predicted by the IFS and stress models at three different times; namely, onset of second ramp (𝑡 = 137 min), onset of second hold (𝑡 = 285 min), and at the end of the cycle (𝑡 = 340 min). Bending moment diagrams are essential in predicting the final deformed shape of the laminate after its removal from the tool. The curved geometry of the laminate in the corner region plays an important role in the development of residual stress and eventually 108  defining the shape of the axial force and bending moment diagrams.  Also, higher values of bending moment at the corner of the laminate potentially results in high spring-in or spring-out after removal from the mold.  5.4 Three-Phase Examples Some numerical case studies that are designed to highlight the various features of the model formulation are presented here to assess the performance of the developed three-phase model for different initial volume fractions of the resin and air (𝜑0𝑟 and 𝜑0𝑎). For the sake of simplicity in analyzing the results and evaluating the performance of the code, the surface tension at the interface between the resin and air phases is ignored. Also, in all the examples considered here, the effective permeability of each fluid phase, e.g. 𝐒𝐹 and 𝐒?̂? in Eq. (4-50), are assumed to vary linearly with respect to the degree of saturation of that fluid phase, i.e. 𝐒𝐹 = 𝜓𝐹𝐒 and 𝐒?̂? = 𝜓?̂?𝐒 as discussed in Section 2.1.3 where 𝐒 is the absolute permeability of the fibre-bed.   5.4.1 Debulking of a Column of Partially Saturated Composite A column of partially saturated (𝜑0𝑎 ≠ 0) and unsolidified (𝜆 = 0) composite system (i.e. a composite with an uncured or liquid resin and initial air content) is subjected to a mechanical loading of intensity 𝑞 = 50 kPa at the top permeable surface as shown in Figure 5.43. A zero pressure condition applied at the top surface on both fluid phases (resin and air), 𝑃𝑏𝑟 = 0 and 𝑃𝑏𝑎 = 0 as defined in Table 3.1 corresponds to the atmospheric pressure (as bulk modulus of air is 𝐾𝑎 = 𝑃𝑎 + 𝑃𝑎𝑡𝑚). The simulation is performed at a constant temperature of 20℃ (constant viscosity of  1.177 × 104 Pa.s for the resin and 1.813 × 10−5 Pa.s for the air). The composite 109  sample is discretized using 20 elements along its height and the initial fibre volume fraction is taken to be 0.58 (𝜑0 = 𝜑0𝑟 + 𝜑0𝑎 = 0.42). To guarantee convergence (accuracy) of the numerical results, a very small time step (∆𝑡 = 0.001 s) is used to capture the flow of air.  Figure 5.44 shows the predicted pressure for three different values of the initial air content, 𝜑0𝑎, and at different times after the application of the load. As inertia effects are neglected, a drop in pressure is observed instantaneously after the load is applied. As expected, at any instant of time the pressure drops as we approach the top permeable surface (𝑧/𝐻 = 1) so that the extrapolated pressure at the top surface vanishes due to the applied boundary condition. Also, for all initial air contents, the pressure decreases with time since the combined flow of resin and air out of the system leads to a greater share of the load being carried by the fibre-bed. This behaviour is a consequence of the permeable condition at the top surface. The lower pressure at higher initial air contents (𝜑0𝑎) also reflects the effect of the relatively high compressibility of air compared to the resin which results in a greater portion of the load being carried by the fibre-bed.  The air content (𝜑𝑎) along the height of the column at different times is also shown in Figure 5.45. It can be seen that 𝜑𝑎 decreases with time as more air flows out of the column due to its very low viscosity. Most of the air evacuates from the column and also compresses very quickly after the application of the load. This shows the influence of the applied load on porosity in a composite system when gas flow is taken into account.   110  The displacement of the top surface, i.e. total consolidation of the column, versus time is shown in Figure 5.46. Once again, for samples with higher initial air content, the consolidation becomes more pronounced due to the higher compressibility of the material. Also, the two main steps in the displacement-time graph can be distinguished. The first step is mainly due to compression and occurs soon after the load is applied to the sample. This sudden displacement is then followed by a very slow consolidation due to flow of fluids (primarily air, because of its very low viscosity) through the top permeable surface.  5.4.2 Curing Isotropic Composite Material – Permeable B.C. This example considers a three-phase composite system with permeable and unconstrained (free) top surface (top right example in Figure 5.14) that undergoes a curing temperature cycle as shown in Figure 5.16.   The resin and air flow velocities through the top permeable boundary of the system for different initial volume fractions of resin and air are shown in Figure 5.47 and Figure 5.48. Comparison of these figures reveals a significant difference in the resin and air velocity as expected due to their marked difference in viscosity.  Figure 5.49 shows the effect of resin cure shrinkage on the transverse strain component, 𝜀𝑧, for samples that contain resin (i.e. fully or partially saturated sample with 𝜑0𝑎 < 42%). However, at the other extreme when the composite sample contains no resin and the matrix material only consists of air (𝜑0𝑎 = 42%), all the strain in the system is finally recovered as thermal expansion/contraction is the sole source of strain. It is also interesting to note that for the fully 111  saturated sample (𝜑0𝑎 = 0%), a considerable amount of strain is induced during the heat-up stage (𝑡 < 52 min). As the resin is highly viscous, it imposes some tensile strain on the system while it expands and flows out of the sample. However, in the presence of air in the system (𝜑0𝑎 > 0%), air evacuates promptly and, in so doing, induces very little strain in the system as it flows out of the top permeable surface due to its very low viscosity. The resin expands and fills the space left behind in the system in the wake of air evacuation. The minor strain observed during heat up for these cases (𝜀𝑧 ≅ 0.0025) stems from the expansion of the fibre-bed only. Same argument can be used to explain the considerable difference in resin flow velocity through the top permeable surface for fully saturated (𝜑0𝑎 = 0%) and partially saturated (𝜑0𝑎 = 10% and 𝜑0𝑎 = 20%) conditions observed in Figure 5.47.    The time history of the volume fractions of air and resin shown in Figure 5.50 confirms this phenomenon where it can be seen that the air volume fraction decreases due to its evacuation from the system while that of the resin increases when it fills the available space left behind. As observed for the case where the matrix material solely consists of air (𝜑0𝑎 = 42%), the volume fraction of the air remains almost constant as air flows out of the system through the top surface and cannot induce any significant strain in the system due to its low viscosity.  Finally as shown in Figure 5.51, the lower stiffness of samples with more air content, results in insignificant values of stress in the constrained 𝑥-direction. As the amount of gas in the composite system increases, the three-phase solid composite material will be less stiff after 112  solidification. This lower stiffness results in lower values of induced tensile stress in the constrained 𝑥-direction during the final cool-down stage of processing.  5.4.3 Curing Isotropic Composite Material – Impermeable B.C.  This example is identical to the previous example with the exception that the boundary condition at the top surface is impermeable (bottom right sample in Figure 5.14). As seen in Figure 5.52, and compared to the previous example (Figure 5.49), the impermeable condition induces a considerable amount of strain in the system during the heat-up regime of the cure cycle. Impermeability inhibits the flow of fluid out of the system while thermal expansion results in transverse tensile strain in the sample. With introduction of more air content in the system, the induced strain increases due to the higher expansion of the air. Again, when the matrix material solely consists of air (𝜑0𝑎 = 42%), since there is no cure shrinkage effect and the only source of strain is thermal, which is reversible, a full recovery of the strain is achieved. For samples that include resin, cure shrinkage effect is easily detectable in the strain versus time curve starting at around 50 min.  For this example, both the resin and air volume fractions, 𝜑𝑟 and 𝜑𝑎 shown in Figure 5.53, increase at early stages of the heat-up phase, however, the rate of increase of resin volume fraction is less due to its lower expansion. Furthermore, with advancement of cure in the resin, the increase in the resin volume fraction is negated by the contraction due to cure shrinkage.   113  5.4.4 Unsaturated Unidirectional Angle Laminate Undergoing Compaction and Cure In this example, the same L-shaped unidirectional [0°] laminate in Figure 5.33 undergoes curing using 3IFS-Othotropic model. The laminate initially contains different amounts of air, 𝜑0𝑎 =0, 0.1, 0.2, and 0.3, distributed uniformly over the domain. For all cases considered, initial volume fraction of the fibres is the same as other previous examples, i.e. 𝜑0𝑓 = 0.58.  The initial and final air, resin, and fibre volume fractions averaged over the laminate surface are shown in Figure 5.54. For all cases, the final volume fraction of air (blue color) is less than its initial value due to its evacuation as a result of the applied load. For 𝜑0𝑎 = 0.1, all the air in the laminate is removed as a result of the applied force. For the fibre (green color), regardless of the initial air content, an almost equivalent final fibre volume fraction is achieved corresponding to full compaction of the laminate where the applied load is entirely resisted by the fibre-bed. For an initially fully saturated laminate (𝜑0𝑎 = 0), the resin volume fraction (red color) decreases after compaction. However, when there is initial air content in the system, the resin fills the area provided as a result of air evacuation. Consequently, in each case the final resin volume fraction rises above its initial value. In other words, as the applied load fully compacts the laminate before resin gelation, the combined resin and air/gas volume fractions for all cases are almost the same as shown in Figure 5.54(b). The final distribution of air in the cured laminate is also shown in Figure 5.55 for the case where higher air content at the corner of laminate is observed. Hubert and Poursartip [1,80] obtained levels of voids in the flat and corner regions of the cured laminate from a visual inspection. They graded the level of void in the final cured part from no voids (0) to many voids (3). For [0°] laminate on the convex tool with permeable boundaries, voids level 114  ranked at highest level (3) at the corner while no void (level 0) was observed for the flat part. This observation is in total qualitatively agreement of the three-phase IFS model prediction shown in Figure 5.55.  Figure 5.56 shows the normal displacement of key points A and C (as shown in Figure 5.33-b) versus time for different initial air content. For all cases, the load is entirely taken by the fibre-bed as indicated by the plateau region of the graphs. For the cases containing air (𝜑0𝑎 = 0.1, 0.2, and 0.3) an instantaneous displacement after application of the load is observed which is due to air evacuation from the system. For higher air content, i.e. 𝜑0𝑎 =0.2 and 0.3, the full compaction is achieved only because of air removal. However, for 𝜑0𝑎 = 0.1, the less instantaneous displacement is followed by further displacement within 50 min because of resin flow through permeable boundaries after its viscosity drops at higher temperature. No instantaneous displacement is seen for the case with no air, 𝜑0𝑎 = 0. In this condition only a gradual displacement is observed as a result of resin flow at higher temperatures.   It is noted that for samples with final porosity (void), residual stress and the consequential wrapage after tool removal may not matter as from a manufacturing standpoint the composite product is considered to have failed (rejected). For samples without final porosity, the final residual stress would be the same as those predicted using the two-phase model.   115  5.5 Tables 3501-6 resin (White and Kim [81], Zobeiry [82]) AS4 fibre (White and Kim [81]) Unrelaxed bulk modulus, 𝐾𝑢 (GPa) 3.556 Poisson’s ratio in-plane, 𝜈12 0.2 Relaxed bulk modulus, 𝐾𝑟 (GPa) 0.8 Poisson’s ratio out-of-plane, 𝜈23 0.3 Unrelaxed shear modulus, 𝐺𝑢 (GPa) 1.185 Longitudinal elastic modulus, 𝐸1 (GPa) 207 Relaxed shear modulus, 𝐺𝑟 (GPa) 1 × 10−6 Transverse elastic modulus, 𝐸2 (GPa) 20.7 Glassy 𝛼𝑡ℎ (T < Tg), (/℃) 5.76 × 10−5 In-plane shear modulus, 𝐺12 (GPa) 27.6 Rubbery 𝛼𝑡ℎ (T ≥ Tg), (/℃) 1.854 × 10−4 Longitudinal 𝛼𝑡ℎ (/°C) −0.9 × 10−6 Volumetric 𝛼𝑐𝑠  −0.055 Transverse 𝛼𝑡ℎ (/°C) 7.2 × 10−6 Atmospheric air (Dixon [83])    Bulk Modulus (GPa) Absolute Pressure of air   Shear Modulus (GPa) 1 × 10−11   Volumetric 𝛼𝑡ℎ  (/℃)  1/𝑇   Table 5.1 Mechanical properties of the fibre, resin and atmospheric air   Resin’s evolution of degree of cure, 𝝌 (Lee et al. [75]) Resin viscosity, 𝝁 (Pa.s) (Hubert et al. [15]) 𝑑𝜒𝑑𝑡= {(𝐶1 + 𝐶2𝜒)(1 − 𝜒)(𝐵 − 𝜒)        𝜒 ≤ 0.3𝐶3(1 − 𝜒)                                        𝜒 > 0.3 𝜇 = 𝜇∞𝑒𝑈 𝑅𝑇⁄ 𝑒𝜅𝜒 𝐶𝐼 = 𝐴𝐼𝑒−∆𝐸𝐼 𝑅𝑇⁄  𝐵 = 0.47  𝐴1 = 3.5017 × 107 𝑠−1 𝐴2 = −3.3567 × 107 𝑠−1 𝐴3 = 3.2667 × 103 𝑠−1   ∆𝐸1 = 80700 𝐽 𝑚𝑜𝑙⁄  ∆𝐸2 = 77800 𝐽 𝑚𝑜𝑙⁄  ∆𝐸3 = 56600 𝐽 𝑚𝑜𝑙⁄   𝜇∞ = 4.6 × 10−17 𝑃𝑎. 𝑠 𝑈 = 114477 𝐽 𝑚𝑜𝑙⁄  𝜅 = 14.8 Air viscosity, 𝝁 (Pa.s) (Dixon [83])  𝜇 =1.458 × 10−6 𝑇1.5110.4 + 𝑇  Table 5.2 Cure kinetics and viscosity models for 3501-6 resin and atmospheric air    116  Fibre-bed permeability (Gutowski et al. [14]) Fibre-bed elastic properties (Haghshenas [17], Hubert [15], Cai and Gutowski [70])  𝑆1 =𝑟𝑓24𝑘(1 − 𝑉𝑓)3𝑉𝑓2  ,  𝑆2 =𝑟𝑓24𝑘′(√𝑉𝑎′ 𝑉𝑓⁄ − 1)3𝑉𝑎′ 𝑉𝑓⁄ + 1 {      𝜎1𝜎2𝜎3𝜏23𝜏13𝜏12}      =[      𝐸1 0 0 𝐷22 𝐷22 − 2𝐺23  𝐷220 0 00 0 00 0 0    𝑆𝑦𝑚𝑚.     𝐺23 0 0 𝐺12 0  𝐺12]      {      𝜀1𝜀2𝜀3𝛾23𝛾13𝛾12}       𝑉𝑓 = 1 − 𝜑 𝑉𝑎′ = 0.81 𝑘 = 0.7 𝑘′ = 0.2 𝑟𝑓 = 4 × 10−6 m 𝐸1 = 100 GPa 𝐺12 = 𝐺23 = 0.001  MPa 𝐷22(MPa) =3𝜋𝐸𝑓𝑙𝑒𝑥𝜔4(5𝑣𝑚𝑎𝑥 − 4𝑣0𝑣𝑚𝑎𝑥 − 1)(𝑣𝑚𝑎𝑥 − 1)5 𝐸𝑓𝑙𝑒𝑥 = 230 GPa  𝑣𝑚𝑎𝑥 = √𝑉𝑓𝑚𝑎𝑥 𝑉𝑓 ,     𝑣0 = √𝑉0 𝑉𝑓 ,    𝑉0 = 0.45 ,   𝑉𝑓𝑚𝑎𝑥 = 0.81 𝜔 = 350 𝐷22    : Transverse component of the constitutive stiffness matrix of the fibre-bed, 𝐾𝑓𝑏 = 𝐷22 − 4𝐺23/3  𝐸𝑓𝑙𝑒𝑥 : Flexural modulus of the fibre 𝜔       : Waviness ratio of fibre 𝑉𝑓𝑚𝑎𝑥 : Maximum fibre volume fraction achievable 𝑉0      : Fibre volume fraction for the uncompacted fibre-bed   Table 5.3 Fibre-bed properties for AS4/3501-6 CFRP                     117  𝐾 = 𝐾𝑟 + (𝐾𝑢 − 𝐾𝑟)∑𝑊?̅?𝑛𝜔=1exp (−𝑡𝑒𝑎𝑇(𝜒)𝜏?̅?(𝜒)) 𝐺 = 𝐺𝑟 + (𝐺𝑢 − 𝐺𝑟)∑𝑊?̅?𝑛𝜔=1exp (−𝑡𝑒𝑎𝑇(𝜒)𝜏?̅?(𝜒)) log(𝜏?̅?(𝜒)) = log(𝜏?̅?(𝜒0)) + (𝑓′(𝜒) − (𝜒 − 𝜒0) log(𝜆?̅?′ )) 𝑓′(𝜒) = −9.3694 + 0.6089𝜒 + 9.1347𝜒2 𝜆?̅?′ =𝜏𝑝(𝜒0)𝜏?̅?(𝜒0) 𝑡𝑒 = −log(𝑒)𝑐1(𝜒𝑓)𝑚                                 for 𝑡 > 𝑡𝑓  𝑡𝑒 = −log(𝑒)𝑐1(𝜒𝑓)𝑚+ (𝑡𝑓 − 𝑡)𝑎𝑇(𝜒)𝑎𝑇(𝜒𝑓)     for 𝑡 ≤ 𝑡𝑓 log(𝑎𝑇) = 𝑐1(𝜒)𝑇 + 𝑐2(𝜒) 𝑐1(𝜒) = −𝑎1 exp (1𝜒−1) − 𝑎2  𝑐2(𝜒) = −𝑇0𝑐1(𝜒) 𝜏𝑝 : peak relaxation time (7.94× 109 min [76]) 𝑡𝑒 : effective time 𝑡𝑓 : time at the onset of cool-down 𝜒0: reference degree of cure (0.98 [76]) 𝜒𝑓 : degree of cure at the onset of cool-down 𝑚 : rate of cool-down 𝑎𝑇(𝜒) : shift factor 𝑇0 : reference temperature (30℃  [76]) 𝑎1 = 1.4 /℃ [76] 𝑎2 = 0.0712/℃  [76]  Table 5.4 Fibre-bed Pseudo-Viscoelastic (PVE) formulation used for the 3501-6 resin (Zobeiry et al. [51])  ?̅? 𝒕𝝎 (min) 𝑾𝝎 1 2 3 4 5 6 7 8 9 2.92E+1 2.92E+3 1.82E+5 1.10E+7 2.83E+8 7.94E+9 1.95E+11 3.32E+12 4.92E+14 0.059 0.066 0.083 0.112 0.154 0.262 0.184 0.049 0.025  Table 5.5 Relaxation times and weight factor for 3501-6 resin for 𝝌𝟎 = 𝟎. 𝟗𝟖 (Kim and White [76])  118  5.6 Figures  Figure 5.1 Assumed variation of the solidification factor with the degree of cure    Figure 5.2 Assumed variation of solid and solid-skeleton properties as a function of the solidification factor (shown for bulk and shear moduli as example)   for for for 119  (a)  (b)   Figure 5.3 Schematic of the miromechanics used in this study (a) two-phase (b) three-phase        Figure 5.4 A solid material under mechanical and thermal loading Fibre        ( )Two-phase Composite( )Resin ( )Micromechanics IAir ( )Resin ( )MatrixFibre        ( )Three-phase Composite( )Matrix  ( )Micromechanics IMicromechanics IIH = 4 mmL = 12 mmzxFz = 8400 NFzFx = 700 NFx120     Figure 5.5 Comparison of the strains obtained using the current IFS model and Abaqus 4-noded plane strain element for fully solidified material   Figure 5.6 Comparison of the stresses obtained using the current IFS model and Abaqus 4-noded plane strain element for fully solidified material -0.07-0.06-0.05-0.04-0.03-0.02-0.01000.050.10.150.20.250.30.350.40 3.7 7.4 11.1 14.8 18.5 22.2 25.9 29.6 33.3 37Strain_xTime (s)IFS - Strain_xABAQUS - Strain_xIFS - Strain_zABAQUS - Strain_z-1200-1000-800-600-400-20000 3.7 7.4 11.1 14.8 18.5 22.2 25.9 29.6 33.3 37Strain_xTime (s)IFS - Strain_xABAQUS - Strain_xIFS - Strain_zABAQUS - Strain_z(MPa)121   Figure 5.7 A solid material under mechanical and thermal loading to investigate the effect of time step on the the IFS model predictions    Figure 5.8 Load-displacement for top surface of the example shown in Figure 5.7                       (a)                                        (b) MPaMPaf=2600 MPaH=4 mmLoad is applied as a ramp in ~460 sL=4 mm-3000-2500-2000-1500-1000-5000-2.5 -2 -1.5 -1 -0.5 0Load (MPa)Displacement of Top Surface (mm)IFS Model - LinearIFS Model - Nonlinear - dt=1sIFS Model- Nonlinear - dt=42.6 sIFS Model - Nonlinear - delt=85.2 sIFS Model - Nonlinear - dt=231 sABAQUS - LinearABAQUS - Nonlinear - dt=1 sABAQUS - Nonlinear - dt=42.6 sABAQUS - Nonlinear - dt=85.2 sABAQUS - Nonlinear - dt=231 sIFS model predictions approach ABAQUS results (with nonlinear geometry option) as smaller time steps are usedIFS and ABAQUS results coincide for linear analysisssssssss122   Figure 5.9 An unsolidified composite system subjected to uniaxial strain loading (a) permeable boundary condition at the top surface only, and (b) impermeable boundary condition at all surfaces      Figure 5.10 Predicted time-histories of the resin and air pressures for both examples shown in Figure 5.9(a) (permeable top surface) and Figure 5.9(b) (impermeable top surface)  H = 4 mmL = 1 mmzxqHLzxq q = 50 kPa0 0.0005 0.001 0.0015 0.00201020304050600 1000 2000 3000 4000Pressure (kPa)Time (min)Impermeable Top - ResinPermebale Top - ResinImpermeable Top - AirPermeable Top - Air123   Figure 5.11 Predicted time-histories of the resin and air velocities for both examples shown in Figure 5.9(a) (permeable top surface)   Figure 5.12 Predicted time-histories of the sample vertical strain containing resin and air as matrix for both examples shown in Figure 5.9(a) (permeable top surface) and Figure 5.9(b) (impermeable top surface) 0 0.0005 0.001 0.0015 0.00205101520250.00E+002.00E-084.00E-086.00E-088.00E-081.00E-071.20E-071.40E-071.60E-071.80E-070 1000 2000 3000 4000Resin Velocity (mm/s)Time (min)Permebale Top - ResinPermeable Top - Air AirVelocity (mm/s)0 0.0005 0.001 0.0015 0.002-0.1-0.09-0.08-0.07-0.06-0.05-0.04-0.03-0.02-0.0100 1000 2000 3000 4000Pressure (kPa)Time (min)Impermeable Top - ResinPermebale Top - ResinImpermeable Top - AirPermeable Top - Air124   Figure 5.13 Predicted time-histories of the resin and air volume fractions for both examples shown in Figure 5.9(a) (permeable top surface) and Figure 5.9(b) (impermeable top surface)  Figure 5.14 Geometry and different boundary conditions of an isotropic (e.g. a unidirectional composite with fibre directions normal to the plane being analyzed) sample undergoing curing 0 0.0005 0.001 0.0015 0.0020.360.370.380.390.40.410.420.430 1000 2000 3000 4000Pressure (kPa)Time (min)Impermeable Top - ResinPermebale Top - ResinImpermeable Top - AirPermeable Top - AirHLzxH = 4 mmL=12mmzxHLzxHLzxPermeable - Constrained Permeable - FreeImpermeable - Constrained Impermeable - Free125    Figure 5.15 Possible location of boundary conditions for the examples shown in Figure 5.14 within a representative practical geometry  (a)                                                                         (b)   Figure 5.16 (a) Time-history of applied temperature and predicted viscosity, (b) predicted degree of cure, and solidification factor of material for the case studies shown in Figure 5.14  4231zxzxzxzx1.00E-011.00E+001.00E+011.00E+021.00E+031.00E+04040801201602000 50 100 150 200Viscosity (Pa.s)Time (min)TemperatureResin Viscosity(Pa.s)Temperature ()00.20.40.60.811.2040801201602000 50 100 150 200Viscosity (Pa.s)Time (min)TemperatureDOClambda (Pa.s) and    Temperature (℃)126   Figure 5.17 Predicted time-histories of the resin pressure for all examples shown in Figure 5.14      Figure 5.18 Predicted time-histories of the resin velocity in the vertical direction at the top surface for examples shown in Figure 5.14 with permeable top surface (two top cases)   -10010203040500 10 20 30 40 50 60Pressure (MPa)Time (min)Permeable     – ConstrainedPermeable     – FreeImpermeable – ConstrainedImpermeable – Free-8.00E-05-4.00E-050.00E+004.00E-058.00E-051.20E-041.60E-040 20 40 60Velocity (mm/s)Time (min)Permeable     – ConstrainedPermeable     – Free127      Figure 5.19 Predicted time-histories of the resin volume fraction for all examples shown in Figure 5.14   Figure 5.20 Predicted time-histories of stresses for the examples shown in Figure 5.14 with resin matrix and constrained top surface (two left cases)   0.410.420.430.440 50 100 150 200 250Matrix Volume FractionTime (min)Permeable     – ConstrainedPermeable     – FreeImpermeable – ConstrainedImpermeable – FreeIntegrated ModelStress Model-80-40040801200 50 100 150 200 250Stress_X (MPa)Time (min)(MPa)both stress model graphsand impermeable – constrained from IFS modelPermeable     – ConstrainedImpermeable – ConstrainedIntegrated ModelStress Model128   Figure 5.21 Predicted time-histories of vertical strain for the examples shown in Figure 5.14 with resin matrix and free top surface (two right cases)     Figure 5.22 Time-history of applied temperature and predicted viscosity of air     -0.0200.020.040 50 100 150 200 250Strain_Z Time (min)both stress model graphsand impermeable – free from IFS modelPermeable     – FreeImpermeable – FreeIntegrated ModelStress Model1.50E-051.70E-051.90E-052.10E-052.30E-052.50E-052.70E-052.90E-05040801201602000 50 100 150 200Viscosity (Pa.s)Time (min)TemperatureAir Viscosity (Pa.s)Temperature (℃)129   Figure 5.23 Predicted time-histories of the air pressure for all examples shown in Figure 5.14    Figure 5.24 Predicted time-histories of the air velocity in the vertical direction at the top surface for examples shown in Figure 5.14 with permeable top surface (two top cases)   -0.0100.010.020.030.040.050.060 50 100 150 200 250Pressure (MPa)Time (min)Permeable     – ConstrainedPermeable     – FreeImpermeable – ConstrainedImpermeable – FreeIntegrated ModelStress Model-4.00E-04-3.00E-04-2.00E-04-1.00E-040.00E+001.00E-042.00E-043.00E-044.00E-040 50 100 150 200 250Velocity (mm/s)Time (min)Permeable     – ConstrainedPermeable     – Free130   Figure 5.25 Predicted time-histories of the air volume fraction for all examples shown in Figure 5.14   Figure 5.26 Predicted time-histories of stresses for examples shown in Figure 5.14 with air matrix and constrained top surface (two left cases)    0.410.420.430.440.450.460.470.480 50 100 150 200 250Matrix Volume FractionTime (min)Permeable     – ConstrainedPermeable     – FreeImpermeable – ConstrainedImpermeable – FreeIntegrated ModelStress Model-0.06-0.05-0.04-0.03-0.02-0.0100.010 50 100 150 200 250Stress_X (MPa)Time (min)(MPa)Permeable     – ConstrainedImpermeable – ConstrainedIntegrated ModelStress Model131    Figure 5.27 Predicted time-histories of vertical strain for free samples at top surface for examples shown in Figure 5.14 with air matrix and free top surface (two right cases)    Figure 5.28 Geometry and boundary condition of orthotropic sample undergoing curing  -0.0200.020.040.060.080.10.120 50 100 150 200 250Strain_Z Time (min)Permeable     – FreeImpermeable – FreeIntegrated ModelStress ModelH= 4 mmzxL= 4 mm132   Figure 5.29 Final deformed shape of the orthotropic sample shown in Figure 5.28 for three different fibre angles (deformations are scaled by a factor of 10 for illustration purpose)    (a)                                                                         (b)   Figure 5.30 Predicted time-histories of the resin velocity through the (a) right vertical surface (b) top horizontal surface for the example shown in Figure 5.28    040 4040 4040 4Deformed Shape – Stress ModelInitial ShapeDeformed Shape – IFS Model-4.00E-050.00E+004.00E-058.00E-051.20E-041.60E-040 20 40 60Velocity XTime (min) =0° =45° =90°(mm/s) -4.00E-050.00E+004.00E-058.00E-051.20E-041.60E-040 10 20 30 40 50 60Velocity YTime (min) =90° =45° =0°(mm/s)133    (a)                                                                         (b)   Figure 5.31 Predicted time-histories of the (a) normal strain for 𝜷 = 𝟎° and 𝜷 = 𝟗𝟎° (b) shear strain for 𝜷 = 𝟒𝟓° for the example shown in Figure 5.28       (a)                                                                         (b)   Figure 5.32 Predicted time-histories of the stress for 𝜷 = 𝟒𝟓° (a) normal stress (b) shear stress for the example shown in Figure 5.28    -0.03-0.02-0.0100.010.020.030.040 50 100 150 200 250Strain XTime (min)Stress ModelIFS Model-0.005-0.004-0.003-0.002-0.00100.0010.0020.0030 50 100 150 200 250Strain XYTime (min)Stress ModelIFS Model-0.1-0.08-0.06-0.04-0.0200.020 50 100 150 200 250Strain XTime (min)IFS ModelStress Model(MPa)-100-80-60-40-200204060801000 50 100 150 200 250Strain XYTime (min)(MPa)Stress ModelIFS Model134   (a)                                                                             (b)  Figure 5.33 (a) Geometry and boundary conditions of the convex angle laminate; and (b) location of key points and the centre-line of the angle laminate  (a)                                                                         (b)   Figure 5.34 (a) Time history of autoclave and part temperature and predicted resin viscosity (b) predicted degree of cure and solidification factor for the case study shown in Figure 5.33   L = 63 mmW = 4.96 mmR = 4.57 mmf = 540 kPa= 0.42WR1ABC1.00E-011.00E+001.00E+011.00E+021.00E+031.00E+040204060801001201401601800 100 200 300Viscosity (Pa.s)Time (min)Autoclave Temp.Part Temp.Series3(Pa.s)Temperature ()00.20.40.60.811.20204060801001201401601800 100 200 300Viscosity (Pa.s)Time (min)Part Temp.Series2lambdaand Temperature ()135   Figure 5.35 Predicted time history of resin bulk and shear moduli based on pseudo-viscoelastic model for the case study shown in Figure 5.33  Figure 5.36 FE meshes used in the analysis of angle laminates  00.20.40.60.811.200.511.522.533.540 50 100 150 200 250 300 350Viscosity (Pa.s)Time (min)and (GPa)-12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44-12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44-12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44-12 -8 -4 0 4 8 12 16 20 24 28 32 36 40 44136   Figure 5.37 Predicted deformed shape of the laminate for the case study shown in Figure 5.33 (deformations are scaled up by a factor of 2 for illustration purposes)      Figure 5.38 Photograph of the cross-section of [0°] laminate on a convex tool with permeable boundary (bleeding system) showing the corner thickening (with permission from Hubert [1])           Initial shapeDeformed shapecorner thickeningedge lateral deformation137      Figure 5.39 Predicted deformed shape of a [0°] laminate on a convex tool with permeable boundary (bleeding system) at end of compaction using flow model (with permission from Hubert [1])            (a)                                                                            (b)   Figure 5.40 (a) Normal displacement of critical points A, B, and C through the curing process (b) tangential displacement (in 𝝃 direction in Figure 5.37) for different through thickness position (in 𝜼 direction in Figure 5.37) of the laminate edge at different times for the case study shown in Figure 5.33 -0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.100.10.20 50 100 150 200 250 300 350Normal Displacement (mm)Time (min)ABCstress model for all the three points (A, B, C)01234560 0.1 0.2 0.3 0.4Vertical Position - (mm)Tangential Displacement -  (mm)time = 20 min40506070340138   Figure 5.41 Predicted fibre volume fraction, (𝟏 − 𝝋), distribution at different instances of time for case study shown in Figure 5.33  (a)                                                                         (b)   Figure 5.42 Predicted distribution of the resultant (a) axial force (b) bending moment, along the mid-plane of the laminate at three different times for the case study shown in Figure 5.33  time = 0 mintime = 20 mintime = 40 mintime = 60 mintime = 50 mintime = 340 min-180-160-140-120-100-80-60-40-200200 10 20 30 40 50 60 70 80Axial Force (N/mm)Distance from corner (mm)137 min285 mint = 340 minStress ModelIFS Model0501001502002500 10 20 30 40 50 60 70 80Moment (N.mm/mm)Distance from corner (mm)137 min285 mint = 340 minStress ModelIFS Model139   Figure 5.43 Geometry and boundary conditions of an unsolidified-partially saturated column representing a slice through the thickness of an unidirectional laminate with fibres transverse to the x-z plane shown and subjected to mechanical loading at the top surface  Figure 5.44 Predicted distribution of pressure along the partially saturated column at different time for the case study shown in Figure 5.43 H= 25 mmL = 2.5 mmq = 50 kPazx0510152025303540450 0.2 0.4 0.6 0.8 1Pressure (kPa)z/H, , , Time = 0.0001 minTime = 0.0005 minTime = 0.001   minTime = 0.002   min140   Figure 5.45 Predicted distribution of air volume fraction along the unsaturated column at different times for the case study shown in Figure 5.43   Figure 5.46 Predicted distribution of air volume fraction along the unsaturated column at different time for the case study shown in Figure 5.43 00.10.20.30.40.50 0.2 0.4 0.6 0.8 1Phi_gz/HTime = 0          minTime = 0.0001 minTime = 0.0005 minTime = 0.001   minTime = 0.002   min-1.6-1.4-1.2-1-0.8-0.6-0.4-0.200 0.0005 0.001 0.0015 0.002Displacement (mm)Time (min)Compression and flow of airOnly flow141    Figure 5.47 Predicted time-histories of the resin velocity for the case study shown in Figure 5.14 with permeable and free top surface (top right figure)    Figure 5.48 Predicted time-histories of the air velocity for the case study shown in Figure 5.14 with permeable and free top surface (top right figure) -1.00E-080.00E+001.00E-082.00E-083.00E-084.00E-080 10 20 30 40 50 60Resin Velocity (mm/s)Time (min)=10%20%0.00E+001.00E-042.00E-040 20 40 600%-1.00E-040.00E+001.00E-042.00E-043.00E-04-5 5 15 25 35 45 55 65Air Velocity (mm/s)Time (min)=42%10%20%142   Figure 5.49 Predicted time-histories of the strain in vertical direction for the case study shown in Figure 5.14 with permeable and free top surface (top right figure)  Figure 5.50 Predicted time-histories of the resin and air volume fractions for the example shown in Figure 5.14 with permeable and free top surface (top right figure) -0.02-0.015-0.01-0.00500.0050.010.0150 50 100 150 200 250StrainTime (min)00.050.10.150.20.250.30.350.40.450 50 100 150 200 250Vr and VgTime (min)and  143   Figure 5.51 Predicted time-histories of the stress in horizontal direction for the example shown in Figure 5.14 with permeable and free top surface (top right figure)  Figure 5.52 Predicted time-histories of the strain in the vertical direction for the example shown in Figure 5.14 with impermeable and free top surface (bottom right figure)  010203040500 50 100 150 200 250(MPa)Time (min)=0%42%20%10%00.020.040.060.080.10 50 100 150 200 250StrainTime (min)144   Figure 5.53 Predicted time-histories of the resin and air volume fractions for the example shown in Figure 5.14 with impermeable and free top surface (bottom right figure)                        00.050.10.150.20.250.30.350.40.450.50 50 100 150 200 250Vr and VgTime (min)and  145  (a)                                                                           (b)  Figure 5.54 Initial and final average volume fraction of a) each phase and b) fibre and combined resin and air for the example shown in Figure 5.33 00.10.20.30.40.50.60.70.80 0.1 0.2 0.3Averahe PhiInitial phigAverage volume fraction ()finalinitial00.10.20.30.40.50.60.70.80 0.1 0.2 0.3Averahe PhiInitial phigfinalinitialAverage volume fraction ()146  (a)                                                                         (b)         Figure 5.55 Final distribution of a) air and b) fibre volume fraction over the domain for 𝝋𝟎𝒂= 𝟎. 𝟑 for the example shown in Figure 5.33   (a)                                                                          (b)  Figure 5.56 Normal displacement for different values of initial air content a) key point A, b) key point C for the example shown in Figure 5.33 -0.7-0.6-0.5-0.4-0.3-0.2-0.100 50 100 150 200 250 300 350Normal Displacement (mm)Time (min)=00.10.2 and 0.3-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.100 50 100 150 200 250 300 350Normal Displacement (mm)Time (min)=00.10.2 and 0.3147  Chapter 6: Summary, Conclusions, and Suggestions for Future Works 6.1 Summary and Conclusions A methodology is presented in this thesis to integrate the simulation of flow of multiple fluids and stress development into a unified modelling framework for processing of composite materials. The model first developed for a material system in which phases are assumed to possess isotropic properties in the plane of interest. The isotropic model first presented for a two-phase system comprising of a solid-skeleton phase (fibre-bed) and a fluid phase (e.g. resin or air). Introducing a state variable to account for solidification of the fluid phase, the governing equations are developed for the general case when the fluid phase (namely the resin phase) evolves during curing from an unsolidified fluid to a fully solid material. In this way, the consistency of the governing equations is maintained with those used for each of the two extremes of processing: i) the fluid flow through porous media formulation applicable to early stage of processing when the fluid phase behaves as a liquid, and ii) the solid mechanics formulation that is used for a solid composite material when the liquid phase is fully solidified while accounting for all principal phenomena that occur during such transformation. Incorporating a second fluid phase in the flow formulation as well as using a three-phase micromechanics for the solid end of processing, the two-phase model is extended to a three-phase integrated flow-stress model that accounts for the presence of multiphase matrix. The isotropic models, 2IFS-Isotropic and 3IFS-Isotropic, are then generalized to a model in which the constituents can possess orthotropic properties, and where the isotropic models can be regarded as a special case of the 3IFS-Orthotropic model.   148  The weak form of the governing equations have been developed within a 2D plane strain finite element modeling framework with solid-skeleton displacement (𝑢), fluid velocities (𝑣 ) and pressures (𝑃), as the primary degrees of freedom where two separate velocity and pressure degrees of freedom are assigned to each fluid phase; namely, gas and liquid, in the three-phase model. The finite element implementation of the model is presented in detail for each of the four cases, i.e. two- or three-phase, isotropic and orthotropic model. The FE equations have been implemented in a MATLAB® code developed in-house that use a backward Euler (implicit) time integration scheme combined with the Newton’s iterative solution to account for both material and geometrical nonlinearities.  By considering air, resin, and carbon fibre as the constituents (phases) of composite materials, the model has been applied to various numerical test cases that include fully solidified composite material under thermal and mechanical loading, compaction of composite material at room temperature with either air or resin as its matrix, debulking of a partially saturated uncured three-phase composite material at room temperature, two- or three-phase composite material that undergoes curing, etc. Where possible, the results are compared to those obtained from the decoupled stress model as well as previous IFS model to show the significance of using the current way of integrated methodology in processing of composite materials. The demonstrated capability of the developed integrated flow-stress model in capturing the essence of the composite material behaviour under processing conditions instills confidence in accuracy of the developed framework to conduct process simulation of complex composite structures.  149  Using the developed framework, more accurate and efficient simulation tool that is based on the physics involved during processing of composite materials will be available. The drawbacks and issues with previous separated flow and stress modules, as well as the previous integrated model are resolved through the new methodology adopted in this work. The model is applicable to composite materials with highly compressible constituents and accounts for multiple phase matrix material. Important process-induced variables such as porosity, volume fraction of phases, residual stresses, pressure, etc. can be predicted using this model. Moreover, the framework can be applied to various types of processing such as debulking (in autoclave processing), infusion, among others. In summary, contribution of this work in process modelling of composite materials can be categorized as:  Integration of flow and stress development into a unified framework with applicability to highly compressible phases  Development of a three-phase model for processing of composite materials  6.2 Suggestions for Future Works The proposed framework is developed in a way that more physics and phenomena can be incorporated as our knowledge and understanding about composite materials during their processing will be enhanced in the future. These phenomena will extend applications of the model to broader types of processing. Some of the complexities and future studies that can be carried out for this purpose are listed in this section.  150  6.2.1 Incorporation of Fluids Phase Exchange In processing of composite materials at high temperature and due to lack of local pressure at some regions of the composite part, the moisture and other volatiles solved in the resin can come out of the solution in their gas state. These dissolved volatiles would show up as porosity in the final composite product. This phenomena can be incorporated into the IFS model as a mass/volume exchange between two fluid phases. A new sub-model needs to control this mass exchange which can be included into the volume fraction update formulation provided in this work as a new exchange term.   6.2.2 Full Form of Governing Equations  As mentioned before, the mass conservation equation used in the current model is based on the small strain assumption. Also, the spatial derivatives of densities and volume fractions are ignored in the mass conservation equations as well as in the Darcy’s law used in this work governing flow of the fluid phase. The body forces and inertia effects are also assumed negligible in governing equations developments. Inclusion of all these terms can be easily included in the equation which makes the framework applicable to broader applications.  6.2.3 Adding More Phases to the Framework There are multiple types of porosity in the composite materials as classified in [62,63], such as resin voids, gaps between different layers of a laminate, channels of gas between fibre bundles in the fibre-bed, etc. The state of gas phase in each of these locations is different and in a more detailed macro model should not be considered as a single gas phase. The current 3-phase model can be extended to a more complicated multiple-phase model where separate phases in the model 151  are assigned to each type of voids mentioned above. Using specific flow and load sharing mechanisms to each of new phases based on the physics that governs their behaviour, provides a more detailed and accurate prediction tool for porosity in composite materials.  6.2.4 Use of Coupled Flow Models As mentioned, decoupled equations are used in the current 3IFS model for flow of each fluid phase. Adopting a more complicated coupled flow mechanism for multiphase matrix, interaction of fluid phases, while they are flowing through porous media, can be considered. In other words, the drag force between two fluid phases would also be taken into account and the gas bubble movements as result of resin flow can be a very useful consequence of bringing in such complexity into the framework.  6.2.5 Full Viscoelastic Material Models In the current 3IFS model, a pseudo viscoelastic material model is used for the resin. Using a full viscoelastic material model can make the developed framework applicable to non-standard complicated cure cycles as well.   6.2.6 𝒖-𝑷 Formualtion  The model can be transformed into a 𝑢-𝑃 formulation which currently is available in some commercial software such as Abaqus. This can be done through elimination of the velocity variable, 𝑣 in the governing equations as result of its expression in terms of pressure using flow equations. Although the adopted 𝑢-𝑣-𝑃 formulation in the work has advantages over the 𝑢-𝑃 152  formulation, a new 𝑢-𝑃 model can be readily implemented in commercial software. Also, it might show higher numerical stability as a results of fewer number of degrees of freedom.  6.2.7 Large Strain Formulation Although the effect of nonlinear geometry is considered in the current model through updating the mesh geometry in each iteration within the nonlinear solution scheme, the model has been developed on the basis of small strains assumption. Using infinitesimal strains along with small time steps, this approach yields results that are just as accurate as the updated Lagrangian method which can accommodate significantly larger time steps. In other words, linearization of the strains (using infinitesimal strains) is valid so long as sufficiently small time steps are used. Using updated Lagrangian method which would include nonlinear terms in the strain-displacement formulation allows coarser time steps to be used. Also, incorporation of nonlinear stress-strain relation along with considering effect of large rotations in definition of strains makes the model beneficial for cases in which these nonlinear terms are indispensable, such as wrinkling simulations.  6.2.8 Other Type of Elements The developed 3IFS framework can be implemented in other kinds of elements than the so-called 4-1 element used in this work. For example, the model can be implemented in a multi-layer element where a layer(s) of shell element is (are) also included. Adding bending stiffness to the elements in an IFS framework could be very useful for forming and wrinkling predictions applications where all of the flow of liquid and gas phases, their compressibility, bending 153  stiffness of laminates, and evolution of properties during processing are essential for an accurate simulation.    6.2.9 Solidification Factor Investigation In this work, a robust framework is devoloped to be used in different kinds of applications and for different material systems. For each application the appropriate input parameters are required. The connection of the fluid flow through porous media and classical solid mechanics end of processing required introduction of a state variable that accounts for the physical state of the fluid phase to determine its solidity, i.e. solidification factor. This requires a wide range of experimental investigations on the solidification factor and its more accurate dependency on the other state variables and parameters involved in the processing such as degree of cure, temperature, volume fractions of phases, architecture of the reinforcement, etc.     154  Bibliography [1] P Hubert, Aspects of Flow and Compaction of Laminated Composite Shapes During Cure, PhD Thesis, The University of British Columbia. (1996).  [2] AC Loos, GS Springer. Curing of epoxy matrix composites, Journal of Composite Materials. 17 (1983) 135-169.  [3] TA Bogetti, JW Gillespie. Two-dimensional cure simulation of thick thermosetting composites, Journal of Composite Materials. 25 (1991) 239-273.  [4] TA Bogetti, JW Gillespie. Process-induced stress and deformation in thick-section thermoset composite laminates, Journal of Composite Materials. 26 (1992) 626-660.  [5] SR White, HT Hahn. Process modeling of composite materials: residual stress development during curin. Part I. model formulation, Journal of Composite Materials. 26 (1992) 2402-2422.  [6] SR White, HT Hahn. Process modeling of composite materials: residual stress development during cure. Part II. Experimental validation, Journal of Composite Materials. 26 (1992) 2423-2453.  [7] A Johnston, An Integrated Model of The Development of Process-Induced Deformation in Autoclave Processing Of Composite Structures, PhD Thesis, The University of British Columbia. (1997).  [8] A Johnston, P Hubert, G Fernlund, R Vaziri, A Poursartip. Process modeling of composite structures employing a virtual autoclave concept, Science and Engineering of Composite Materials. 5 (1996) 235-252.  [9] ARA Arafath, Efficient Numerical Techniques for Predicting Process-Induced Stresses and Deformations in Composite Structures, PhD Thesis, The University of British Columbia. (2007).  [10] K Terzaghi, Theoretical Soil Mechanics, John Wiley and Sons, New York, 1943.  [11] MA Biot. General Theory of Three‐Dimensional Consolidation, Journal of Applied Physics. 12 (1941) 155-164.  [12] MA Biot, DG Willis. The elastic coefficients of the theory of consolidation, Journal of Applied Mechanics. 24 (1957) 594-601.  [13] OC Zienkiewicz, T Shiomi. Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution, International journal for numerical and analytical methods in geomehanics. 8 (1984) 71-96.  155  [14] TG Gutowski, Z Cai, S Bauer, D Boucher, J Kingery, S Wineman. Consolidation experiments for laminate composites, Journal of Composite Materials. 21 (1987) 650-669.  [15] P Hubert, R Vaziri, A Poursartip. A two-dimensional flow model for the process simulation of complex shape composite laminates, International Journal for Numerical Methods in Engineering. 44 (1999) 1-26.  [16] H Darcy, Les Fontaines Publiques De La Ville De Dijon: Exposition et Application, Victor Dalmont, 1856.  [17] SM Haghshenas, Integrating Resin Flow and Stress Development in Process Modeling of Thermoset Composites, PhD Thesis, The University of British Columbia (UBC). (2012).  [18] HC Brinkman. Calculation of viscous force exerted by flowing fluid on dense swarm of particles, Applied Scientific Research. A (1949) 27-34.  [19] VW Wang, CA Hieber, KK Wang. Dynamic simulation and graphics for the injection molding of three-dimensional thin parts, Journal of Polymer Engineering. 7 (1986) 21-45.  [20] CA Fracchia, J Castro, CL Tucker, A finite element/control volume simulation of resin transfer mold filling, Proceedings of American Society for Composites - 4th Technical Conference, Lancaster, PA, (1989) 157-166.  [21] MV Bruschke, SG Advani. A finite element/control volume approach to mold filling in anisotropic porous media, Polymer Composites. 11 (1990) 398-405.  [22] AC Loos, JD MacRae. A process simulation model for the manufacture of a blade-stiffened panel by the resin film infusion process, Composites Science and Technology. 56 (1996) 273-289.  [23] P Simacek, SG Advani. A numerical model to predict fiber tow saturation during liquid composite molding, Composites Science and Technology. 63 (2003) 1725-1736.  [24] N Kuentzer, P Simacek, SG Advani, S Walsh. Correlation of void distribution to VARTM manufacturing techniques, Composites Part A: Applied Science and Manufacturing. 38 (2007) 802-813.  [25] H Tan, KM Pillai. Multiscale modeling of unsaturated flow in dual-scale fiber preforms of liquid composite molding I: Isothermal flows, Composites Part A: Applied Science and Manufacturing. 43 (2012) 1-13.  [26] F Trochu, R Gauvin, Z Zhang, Simulation of Mold Filling in Resin Transfer Molding by Non-Conforming Finite Elements, in: Computer Aided Design in Composite Material Technology III, Springer Netherlands, 1992, pp. 109-120.  156  [27] F Trochu, R Gauvin, D Gao. Numerical analysis of the resin transfer molding process by the finite element method, Advances in Polymer Technology. 12 (1993) 329-342.  [28] RV Mohan, ND Ngo, KK Tamma. On a pure finite‐element‐based methodology for resin transfer molds filling simulations, Polymer Engineering and Science. 39 (1999) 26-43.  [29] RV Mohan, KK Ngo, RV Tamma. Three-dimensional resin transfer molding: isothermal process modeling and implicit tracking of moving fronts for thick, geometrically complex composites manufacturing applications - Part 2, Numerical Heat Transfer, Part A: Applications. 35 (1999) 839-858.  [30] AW Bishop. The effective stress principle, Teknisk Ukeblad. 39 (1959) 859-863.  [31] M Nuth, L Laloui . Effective stress concept in unsaturated soils: clarification and validation of a unified framework, International Journal for Numerical and Analytical Methods in Geomechanics. 32 (2008) 771-801.  [32] BA Schrefler, The Finite Element Method in Soil Consolidation (with applications to surface subsidence), PhD Thesis, University College of Swansea. (1984).  [33] J Bear, AHD Cheng, Modeling Groundwater Flow and Contaminant Transport, Springer Science & Business Media 2010.  [34] J Bear, Y Bachmat, Introduction to Modeling of Transport Phenomena in Porous Media, Kluwer, Dordrecht, 1990.  [35] N Patel, LJ Lee. Modeling of void formation and removal in liquid composite molding. Part II: Model development and implementation, Polymer Composites. 17 (1996) 104-114.  [36] MC Leverett. Capillary behaviour in porous media, Transactions of Society of Petroleum Engineering AIME. 142 (1941) 152-269.  [37] N Patel, LJ Lee. Modeling of void formation and removal in liquid composite molding. Part I: Wettability analysis, Polymer Composites. 17 (1996) 96-103.  [38] J Bear, Dynamics of Fluids in Porous Media, Courier Corporation 2013.  [39] JDH Hughes. The carbon fibre/epoxy interface - A review, Composites Science and Technology. 41 (1991) 13-45.  [40] HL Needles, KW Alger, R Okamoto. Adhesion at the interface in cured graphite fiber epoxy-amine resin composites, Journal of Reinforced Plastics and Composites. 6 (1987) 357-366.  157  [41] SA Page, R Mezzenga, L Boogh, JC Berg, JE Månson. Surface Energetics Evolution during Processing of Epoxy Resins, Journal of Colloid Interface Science. 222 (2000) 55-62.  [42] SA Page, JC Berg, JE Månson. Characterization of epoxy resin surface energetics, Journal of Adhesion Science and Technology. 15 (2001) 153-170.  [43] M Li, C Yuan, SK Wang, YZ Gu, K Potter, ZG Zhang. Evolution of the wettability between carbon fiber and epoxy as a function of temperature and resin curing, Journal of Applied Polymer Science. 128 (2013) 4095-4101.  [44] AR Khoei, T Mohammadnejad. Numerical modeling of multiphase fluid flow in deforming porous media: A comparison between two- and three-phase models for seismic analysis of earth and rockfill dams, Computers and Geotechnics. 38 (2011) 142-166.  [45] G Oettl, RF Stark, G Hofstetter. Numerical simulation of geotechnical problems based on a multi-phase finite element approach, Computers and Geotechnics. 31 (2004) 643-664.  [46] F Nagel, G Meschke . An elasto‐plastic three phase model for partially saturated soil for the finite element simulation of compressed air support in tunnelling, International Journal for Numerical and Analytical Methods in Geomechanics. 34 (2010) 605-625.  [47] B Schümann, Modeling of soils as multiphase materials with Abaqus, Proceedings of SIMULIA Customer Conference 2010 at Providence, Rhode Island/USA (2010) 348-399.  [48] S Holler, K Meskouris, Dynamisches Mehrphasenmodell mit hypoplastischer Materialformulierung der Feststoffphase, PhD Thesis, Fakultät für Bauingenieurwesen. (2006).  [49] A Johnston, R Vaziri, A Poursartip. A plane strain model for process-induced deformation of laminated composite structures, Journal of Composite Materials. 35 (2001) 1435-1469.  [50] N Zobeiry, R Vaziri, A Poursartip. Differential implementation of the viscoelastic response of a curing thermoset matrix for composites processing, Journal of Engineering Materials and Technology. 128 (2005) 90-95.  [51] N Zobeiry, R Vaziri, A Poursartip. Computationally efficient pseudo-viscoelastic models for evaluation of residual stresses in thermoset polymer composites during cure, Composites Part A: Applied Science and Manufacturing. 41 (2010) 247-256.  [52] N Zobeiry, S Malek, R Vaziri, A Poursartip. A differential approach to finite element modelling of isotropic and transversely isotropic viscoelastic materials, Mechanics of Materials. 97 (2016) 76-91.  [53] M Gigliotti, F Jacquemin, J Molimard, A Vautrin. Transient and cyclical hygrothermoelastic stress in laminated composite plates: Modelling and experimental assessment, Mechanics of materials. 39 (2007) 729-745.  158  [54] E Lacoste, K Szymanska, S Terekhina, S Fréour, F Jacquemin, M Salvia. A multi-scale analysis of local stresses development during the cure of a composite tooling material, International Journal of Material Forming. 6 (2013) 467-482.  [55] SR Ghiorse. Effect of void content on the mechanical properties of carbon/epoxy laminates, SAMPE Quarterly. 24 (1993) 54-59.  [56] J Kardos, M Duduković, R Dave. Void growth and resin transport during processing of thermosetting—matrix composites, Epoxy Resins and Composites IV. (1986) 101-123.  [57] FYC Boey, SW Lye. Effects of vacuum and pressure in an autoclave curing process for a thermosetting fibre-reinforced composite, Journal of Materials Processing Technology. 23 (1990) 121-131.  [58] FYC Boey, SW Lye. Void reduction in autoclave processing of thermoset composites: Part 1: High pressure effects on void reduction, Composites. 23 (1992) 261-265.  [59] B Thorfinnson, TF Biermann, Production of void free composite parts without debulking, 31st  International SAMPE Symposium, Las Vegas, Nevada, 1986 (1986) 480-490.  [60] B Thorfinnson, TF Biermann. Degree of Impregnation of Prepregs--Effects on Porosity, Advanced Materials Technology. 87 (1987) 1500-1509.  [61] ARA Arafath, G Fernlund, A Poursartip, Gas transport in prepregs: model and permeability experiments, Proceedings of 17th International Conference on Composite Materials (ICCM), Edinburg, United Kingdom (2009).  [62] L Farhang, Void Evolution During Processing of Out-Of-Autoclave Prepreg Laminates. PhD Thesis, The University of British Columbia. (2014).  [63] L Farhang, G Fernlund. Void and porosity characterization of uncured and partially cured prepregs. Journal of Composite Materials. 50 (2015) 937-948.  [64] J Wells, Behaviour of Resin Voids in Out-Of-Autoclave Prepreg Processing. Master Thesis, The University of British Columbia. (2015).  [65] J Kay, L Farhang, K Hsiao, G Fernlund, Effect of process conditions on porosity in out-of-autoclave prepreg laminates, (2011).  [66] J Kay, Gas Transport and Void Evolution in Out-of-autoclave Composite Prepregs, PhD Thesis (under preparation), The University of British Columbia (2017).  [67] M Roy, Porosity in Configured Structures: Effect of Ply Drops and Caul Sheets in the Processing of Composite Parts, Master Thesis, The University of British Columbia  (2015).  159  [68] R Dave, JL Kardos, MP Duduković. A model for resin flow during composite processing: Part 1—general mathematical development, Polymer Composites. 8 (1987) 29-38.  [69] RW Lewis, BA Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, 2nd ed., Wiley, New York, 1998.  [70] Z Cai, TG Gutowski. The 3-D deformation behavior of a lubricated fiber bundle, Journal of Composite Materials. 26 (1992) 1207-1237.  [71] MM Carroll. An effective stress law for anisotropic elastic deformation, Journal of Geophysical Research: Solid Earth. 84 (1979) 7510-7512.  [72] M Thompson, JR Willis. A reformation of the equations of anisotropic poroelasticity, Journal of Applied Mechanics. 58 (1991) 612-616.  [73] MA Biot. Theory of elasticity and consolidation for a porous anisotropic solid, Journal of Applied Physics. 26 (1955) 182-185.  [74] A Cheng. Material coefficients of anisotropic poroelasticity, International Journal of Rock Mechanics and Mining Sciences. 34 (1997) 199-205.  [75] WI Lee, AC Loos, GS Springer. Heat of reaction, degree of cure, and viscosity of Hercules 3501-6 resin, Journal of Composite Materials. 16 (1982) 510-520.  [76] YK Kim, SR White. Stress relaxation behavior of 3501-6 epoxy resin during cure, Polymer Engineering & Science. 36 (1996) 2852-2862.  [77] RM Christensen. A critical evaluation for a class of micro-mechanics models, Journal of Mechanics of Physics and Solids. 38 (1990) 379-404.  [78] VM Levin. Thermal expansion coefficient of heterogeneous materials, Mechanics of Solids. 2 (1967) 58-61.  [79] S Nie, C Basaran. A micromechanical model for effective elastic properties of particulate composites with imperfect interfacial bonds, International Journal of Solids and Structures. 42 (2005) 4179-4191.  [80] P Hubert, A Poursartip. Aspects of the compaction of composite angle laminates: an experimental investigation, Journal of Composite Materials. 35 (2001) 2-26.  [81] SR White, YK Kim. Process-induced residual stress analysis of AS4/3501-6 composite material, Mechanics of Composite Materials and Structures an International Journal. 5 (1998) 153-186.  160  [82] N Zobeiry, Viscoelastic Constitutive Models for Evaluation of Residual Stresses in Thermoset Composites During Cure, PhD Thesis, The University of British Columbia. (2006).  [83] JC Dixon, The Shock Absorber Handbook, 2nd ed., John Wiley & Sons, Chichester, 2007.      161  Appendices Appendix A  Iterative Solution for System of Governing Equations A.1 Two-Phase System To reach to iterative solution of the two-phase system, Eq. (3-40), two different method is applied for first two governing equation, mass conservation (Eq. (3-27)) and equilibrium of the matrix (Eq.(3-28)), compared to the last equation, equilibrium of the system (Eq.(3-29)). The first and second governing differential equations of the system, can be presented in following matrix form: [𝟎 𝟎 0𝟎 𝐊𝒗𝐹𝒗𝐹 𝐊𝒗𝑷𝐹𝟎 𝐊𝑷𝐹𝒗𝐹 𝟎] {?̅??̅?𝐹?̅?𝐹} + [𝟎 𝟎 𝟎𝟎 𝟎 𝟎𝐂𝑷𝐹𝒖 𝟎 𝐂𝑷𝐹𝑷𝐹] {?̇̅??̇̅?𝐹?̇̅?𝐹} = {𝟎𝐟𝒗𝐹𝐟𝑷𝐹} (A-1) where the last row corresponds to the mass conservation equation (Eq. (3-27)). All matrix and vector components shown in Table 3.1 are obtained through applying the Galerking method combined with integration by parts as presented in [17]. Eq. (A-1) can be written in short form as: 𝐊𝐗 + 𝐂?̇? = 𝐅 (A-2) where: 𝐊 = [𝟎 𝟎 0𝟎 𝐊𝒗𝐹𝒗𝐹 𝐊𝒗𝑷𝐹𝟎 𝐊𝑷𝐹𝒗𝐹 𝟎] (A-3) 𝐂 = [𝟎 𝟎 𝟎𝟎 𝟎 𝟎𝐂𝑷𝐹𝒖 𝟎 𝐂𝑷𝐹𝑷𝐹] (A-4) 𝐅 = {𝟎𝐟𝒗𝐹𝐟𝑷𝐹} (A-5) 162  To solve the first order differential Eq. (A-2) using backward Euler method, the following assumption is made for solution at 𝑛th time step: 𝐗𝑛 = 𝐗𝑛−1 + Δ𝑡?̇?𝑛 (A-6) Substituting ?̇?𝑛 from Eq. (A-2) and rearranging of Eq. (A-6) gives: 𝐗𝑛−1 + Δ𝑡 (𝐅𝑛 − 𝐊𝑛𝐗𝑛𝑪𝑛) − 𝐗𝑛 = 0 (A-7) To simplify mathematic manipulations, the above equation can be regarded as a function in terms of variable 𝐗𝑛 as follows where all other quantities are regarded as constant values: 𝐹(𝐗𝑛) = 0 (A-8) Based on Newton’s method, the iterative solution of the above equation at 𝑘th iteration is given as: 𝐗𝑛𝑘 − 𝐗𝑛𝑘−1 = −((𝝏𝐹𝝏 𝐗𝑛(𝐗𝑛))𝑘)−1(𝐹(𝐗𝑛))𝑘 (A-9) where: (𝝏𝐹𝝏 𝐗𝑛(𝐗𝑛))𝑘= −Δ𝑡𝐊𝑛𝑘−1𝐂𝑛𝑘−1 − 1 (A-10) (𝐹(𝐗𝑛))𝑘= 𝐗𝑛−1 + Δ𝑡 (𝐅𝑛𝑘−1 − 𝐊𝑛𝑘−1𝐗𝑛𝑘−1𝑪𝑛𝑘−1 ) − 𝐗𝑛𝑘  (A-11) Substituting Eqs. (A-10) and (A-11) in Eq. (A-9), we arrive at:  𝐗𝑛𝑘 − 𝐗𝑛𝑘−1 = −(−Δ𝑡𝐊𝑛𝑘−1𝐂𝑛𝑘−1 − 1)−1(𝐗𝑛−1 + Δ𝑡 (𝐅𝑛𝑘−1 −𝐊𝑛𝑘−1𝐗𝑛𝑘−1𝑪𝑛𝑘−1 ) − 𝐗𝑛𝑘) (A-12) Rearranging the above equation results in: (𝐂𝑛𝑘−1 + Δ𝑡𝐊𝑛𝑘−1)𝐗𝑛𝑘  − 𝐂𝑛𝑘−1𝐗𝑛−1 = Δ𝑡𝐅𝑛𝑘−1 (A-13) 163  Substituting quantities in Eq. (A-3)-(A-5) in Eq. (A-13) yields: [𝟎 𝟎 𝟎𝟎 Δ𝑡𝐊𝒗𝐹𝒗𝐹 Δ𝑡𝐊𝒗𝐹𝑷𝐹𝐂𝑷𝐹𝒖 Δ𝑡𝐊𝑷𝐹𝒗𝐹 𝐂𝑷𝐹𝑷𝐹]𝑛𝑘−1{?̅?𝑛𝑘 − ?̅?𝑛𝑘−1?̅?𝐹𝑛𝑘?̅?𝐹𝑛𝑘− ?̅?𝐹𝑛𝑘−1}= {𝟎∆𝑡[𝐟𝒗𝐹𝑛 − 𝐊𝒗𝐹𝑷𝑛𝑘?̅?𝐹𝑛𝑘−1]Δ𝑡𝐟𝑷𝐹𝑛 + 𝐂𝑷𝐹𝒖𝑛𝑘−1(?̅?𝑛−1 − ?̅?𝑛𝑘−1) + 𝐂𝑷𝐹𝑷𝐹𝑛𝑘(?̅?𝐹𝑛−1 − ?̅?𝐹𝑛𝑘−1)} (A-14) Removing all 𝑛 subscripts, defining (  )𝑛−1 = (  )0 as well as ∆( )𝑘 = ( )𝑘 − ( )𝑘−1 in order to simplify appearance of the equations we arrive at: [𝟎 𝟎 𝟎𝟎 Δ𝑡𝐊𝒗𝐹𝒗𝐹 Δ𝑡𝐊𝒗𝐹𝑷𝐹𝐂𝑷𝐹𝒖 Δ𝑡𝐊𝑷𝐹𝒗𝐹 𝐂𝑷𝐹𝑷𝐹]𝑘−1{∆?̅?𝑘?̅?𝐹𝑘∆?̅?𝐹𝑘}= {𝟎∆𝑡[𝐟𝒗𝐹 − 𝐊𝒗𝐹𝑷𝑘?̅?𝐹𝑘−1]Δ𝑡𝐟𝑷𝐹0 + 𝐂𝑷𝐹𝒖𝑘−1(?̅?0 − ?̅?𝑘−1) + 𝐂𝑷𝐹𝑷𝐹𝑛𝑘(?̅?𝐹0 − ?̅?𝐹𝑘−1)} (A-15)  To obtain the first row of iterative solution (3-40), following the Galerking method, the equilibrium equation of the system, Eq. (3-7), over the domain Ω uis integrated using some weight function 𝑤𝑝: ∫𝜎𝑖𝑗,𝑗 Ω𝑤𝑝𝑑Ω = 0 (A-16) Using integration by parts, we arrive at:  −∫𝜎𝑖𝑗 Ω𝑤𝑝,𝑗𝑑Ω + ∫𝜎𝑖𝑗𝑛𝑗𝑤𝑝𝑑Γ Γ𝜎= 0 (A-17) 164  where Γ and Γ𝜎 denote total boundary of the domain and traction boundary of the domain, respectively, and 𝑛𝑗  represents the component of the unit normal vector to the Γ𝜎. Writing the above equation at 𝑘th iteration in current time step of 𝑛 gives: −∫(𝜎𝑖𝑗𝑤𝑝,𝑗)𝑘 Ω𝑑Ω + ∫(𝜎𝑖𝑗𝑛𝑗𝑤𝑝)𝑘𝑑Γ Γ𝜎= 0 (A-18) The temporal discretization of the stress is written as: 𝜎𝑖𝑗𝑘 = 𝜎𝑖𝑗𝑘−1 + 𝑑𝜎𝑖𝑗𝑘 (A-19) where based on Eq. (3-24) we can write: 𝑑𝜎𝑖𝑗𝑘 = 𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝜎 𝑘+1 − 𝑏𝑘−1 𝑑𝑃𝑘𝛿𝑖𝑗 (A-20) Using definition of mechanical strain on deformation and free strains in Eq. (A-20) combined with Eq. (A-19): 𝜎𝑖𝑗𝑘 = 𝜎𝑖𝑗𝑘−1 +12𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1[(𝑢𝑞,𝑟𝑘 + 𝑢𝑟,𝑞𝑘) − (𝑢𝑞,𝑟𝑘−1 + 𝑢𝑟,𝑞𝑘−1)]− 𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝑓𝑆𝐾𝑘− 𝑏𝑘−1(𝑃𝑘 − 𝑃𝑘−1)𝛿𝑖𝑗 (A-21) Substituting Eq. (A-21) in integral Eq. (A-18) yields: 12∫𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑛𝑘−1[(𝑢𝑞,𝑟𝑘 + 𝑢𝑟,𝑞𝑘) − (𝑢𝑞,𝑟𝑘−1 + 𝑢𝑟,𝑞𝑘−1)]𝑤𝑝,𝑗𝑑Ω Ω− ∫𝑏𝑘−1(𝑃𝐹𝑘 − 𝑃𝐹𝑘−1)𝑤𝑝,𝑖𝑑Ω Ω= ∫𝑡𝑖𝑘𝑑Γ Γσ+ ∫𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝑓𝑆𝐾𝑘 Ω𝑤𝑝,𝑗𝑑Ω − ∫𝜎𝑖𝑗𝑘−1 Ω𝑤𝑝,𝑗𝑑Ω (A-22) 165  where 𝑡𝑖 represent component of vector of total tractions applied on Γσ. The above equation can be written in matrix form as: [∆𝑡𝐊𝒖𝒖 𝟎 ∆𝑡𝐊𝒖𝑷𝐹𝟎 𝟎 𝟎𝟎 𝟎 𝟎]𝑘−1{∆?̅?𝑘?̅?𝑘∆?̅?𝑘} = {Δ𝑡[𝐟𝒖 + 𝐟𝒖𝑓𝑆𝐾𝑘−𝐟𝒖𝑖𝑛𝑡𝑘−1]𝟎𝟎} (A-23) where its components are defined in Table 3.1. Combining Eq. (A-14) and (A-23), the final form of iterative matrix equation for the 𝑘th iteration within the current time station 𝑡𝑛 (or the 𝑛th time-step) is written as: [∆𝑡𝐊𝒖𝒖 𝟎 ∆𝑡𝐊𝒖𝑷𝐹𝟎 Δ𝑡𝐊𝒗𝐹𝒗𝐹 Δ𝑡𝐊𝒗𝐹𝑷𝐹𝐂𝑷𝐹𝒖 Δ𝑡𝐊𝑷𝐹𝒗𝐹 𝐂𝑷𝐹𝑷𝐹]𝑘−1{∆?̅?𝑘?̅?𝐹𝑘∆?̅?𝐹𝑘}={    Δ𝑡[𝐟𝒖 + 𝐟𝒖𝑓𝑆𝐾𝑘−𝐟𝒖𝑖𝑛𝑡𝑘−1]∆𝑡[𝐟𝒗𝐹 − 𝐊𝒗𝐹𝑷𝑘?̅?𝐹𝑘−1]Δ𝑡𝐟𝑷𝐹 + 𝐂𝑷𝐹𝒖𝑘−1(?̅?0 − ?̅?𝑘−1) + 𝐂𝑷𝐹𝑷𝐹𝑛𝑘(?̅?𝐹0 − ?̅?𝐹𝑘−1)}     (A-24)  A.2 Three-Phase System The matrix form of first four equations in governing equation of three-phase system, (Eqs. (3-67)-(3-70)), is written as: 166  [      𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝐊𝒗𝐹𝒗𝐹 𝟎 𝐊𝒗𝐹𝑷𝐹 𝟎𝟎 𝟎 𝐊𝒗?̂?𝒗?̂? 𝟎 𝐊𝒗?̂?𝑷?̂?𝟎 𝟎 𝟎 𝐊𝑷𝐹𝑷𝐹 𝐊𝑷𝐹𝑷?̂?𝟎 𝐊𝑷?̂?𝒗𝐹 𝐊𝑷?̂?𝒗?̂? 𝐊𝑷?̂?𝑷𝐹 𝑲𝑷?̂?𝑷?̂?]      {      ?̅??̅?𝐹?̅??̂??̅?𝐹?̅??̂?}      +[      𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝟎 𝟎𝐂𝑷?̂?𝒖 𝟎 𝟎 𝐂𝑷?̂?𝑷𝐹 𝐂𝑷?̂?𝑷?̂?]      {      ?̇̅??̇̅?𝐹?̇̅??̂??̇̅?𝐹?̇̅??̂?}      ={      𝟎𝐟𝒗𝐹𝐟𝒗?̂?𝐟𝑷𝐹𝐟𝑷?̂?}       (A-25) Using the iterative solution (A-13) which is based on backward Euler and Newton solution method and rearranging gives: [      𝟎 𝟎 𝟎 𝟎 𝟎𝟎 ∆𝑡𝐊𝒗𝐹𝒗𝐹 𝟎 ∆𝑡𝐊𝒗𝐹𝑷𝐹 𝟎𝟎 𝟎 ∆𝑡𝐊𝒗?̂?𝒗?̂? 𝟎 ∆𝑡𝐊𝒗?̂?𝑷?̂?𝟎 𝟎 𝟎 ∆𝑡𝐊𝑷𝐹𝑷𝐹 ∆𝑡𝐊𝑷𝐹𝑷?̂?𝐂𝑷?̂?𝒖 ∆𝑡𝐊𝑷?̂?𝒗𝐹 ∆𝑡𝐊𝑷?̂?𝒗?̂? ∆𝑡𝑲𝑷?̂?𝑷𝐹 + 𝑪𝑷?̂?𝑷𝐹 ∆𝑡𝑲𝑷?̂?𝑷?̂? + 𝑪𝑷?̂?𝑷?̂?]      𝑘−1{      ∆?̅?𝑛𝑘?̅?𝐹𝑘?̅??̂?𝑛𝑘∆?̅?𝐹𝑛𝑘∆?̅??̂?𝑛𝑘}       ={      𝟎∆𝑡[𝐟𝒗𝐹 −𝐊𝒗𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1]∆𝑡[𝐟𝒗?̂? − 𝐊𝒗?̂?𝑷?̂?𝑘−1?̅??̂?𝑘−1]∆𝑡[𝐟𝑷𝐹 − 𝐊𝑷𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1− 𝐊𝑷𝐹𝑷?̂?𝑘−1?̅??̂?𝑘−1]𝐅𝑷?̂? }       (A-26) where 𝐅𝑷?̂? = ∆𝑡[𝐟𝑷?̂? − 𝐊𝑷?̂?𝑷𝐹𝑘−1?̅?𝐹𝑘−1− 𝐊𝑷?̂?𝑷?̂?𝑘−1?̅??̂?𝑘−1] + 𝐂𝑷?̂?𝒖𝑘−1(?̅?0 − ?̅?𝑘−1)+ 𝐂𝑷?̂?𝑷𝐹𝑘−1(?̅?𝐹0− ?̅?𝐹𝑘−1) + 𝐂𝑷?̂?𝑷?̂?𝑘−1(?̅??̂?0− ?̅??̂?𝑘−1) (6-1)  167  To obtain the first row corresponds to equilibrium equation of the system, regarding Eq.  (3-66) we can write: 𝑑𝜎𝑖𝑗𝑘 = 𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝜎 𝑘+1− 𝑏𝑛𝑘−1 (𝑑𝑃?̂?𝑘 − 𝜓𝐹𝑘(𝑑𝑃?̂?𝑘 − 𝑑𝑃𝐹𝑘) − (𝜓𝐹𝑘 − 𝜓𝐹𝑘−1)(𝑃?̂?𝑘− 𝑃𝐹𝑘))𝛿𝑖𝑗 (A-27) Substituting Eq. (A-27) in (A-19) and using definition of mechanical strain on deformation and free strains same as two-phase model results in: 𝜎𝑖𝑗𝑘 = 𝜎𝑖𝑗𝑘−1 +12𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1[(𝑢𝑞,𝑟𝑘 + 𝑢𝑟,𝑞𝑘) − (𝑢𝑞,𝑟𝑘−1 + 𝑢𝑟,𝑞𝑘−1)]− 𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝑓𝑆𝐾𝑘− 𝑏𝑘−1 (𝑑𝑃?̂?𝑘 − 𝜓𝐹𝑘(𝑑𝑃?̂?𝑘 − 𝑑𝑃𝐹𝑘) − (𝜓𝐹𝑘 − 𝜓𝐹𝑘−1)(𝑃?̂?𝑘− 𝑃𝐹𝑘))𝛿𝑖𝑗 (A-28) Substituting Eq. (A-28) in integral Eq. (A-18): 168  12∫𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1[(𝑢𝑞,𝑟𝑘 + 𝑢𝑟,𝑞𝑘) − (𝑢𝑞,𝑟𝑘−1 + 𝑢𝑟,𝑞𝑘−1)]𝑤𝑝,𝑗𝑑Ω Ω− ∫𝑏𝑘−1 ((1 − 𝜓𝐹𝑘−1)(𝑃?̂?𝑘 − 𝑃?̂?𝑘−1))𝑤𝑝,𝑖𝑑Ω Ω− ∫𝑏𝑘−1 (𝜓𝐹𝑘−1(𝑃𝐹𝑘 − 𝑃𝐹𝑘−1))𝑤𝑝,𝑖𝑑Ω Ω= ∫𝑡𝑖𝑘𝑑Γ Γσ+ ∫𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝑓𝑆𝐾𝑘 Ω𝑤𝑝,𝑗𝑑Ω − ∫𝜎𝑖𝑗𝑘−1 Ω𝑤𝑝,𝑗𝑑Ω+ ∫𝑏𝑘−1(𝜓𝐹𝑘 − 𝜓𝐹𝑘−1)(𝑃?̂?𝑘−1 − 𝑃𝐹𝑘−1) Ω𝑤𝑝,𝑗𝑑Ω (A-29) If we assume  (𝜓𝐹𝑘 − 𝜓𝐹𝑘−1) ≈ (𝜓𝐹𝑘−1 − 𝜓𝐹𝑘−2) (A-30) , we arrive at: 169  12∫𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1[(𝑢𝑞,𝑟𝑘 + 𝑢𝑟,𝑞𝑘) − (𝑢𝑞,𝑟𝑘−1 + 𝑢𝑟,𝑞𝑘−1)]𝑤𝑝,𝑗𝑑Ω Ω− ∫𝑏𝑘−1 ((1 − 𝜓𝐹𝑘−1)(𝑃?̂?𝑘 − 𝑃?̂?𝑘−1))𝑤𝑝,𝑖𝑑Ω Ω− ∫𝑏𝑘−1 (𝜓𝐹𝑘−1(𝑃𝐹𝑘 − 𝑃𝐹𝑘−1))𝑤𝑝,𝑖𝑑Ω Ω= ∫𝑡𝑖𝑘𝑑Γ Γσ+ ∫𝐷𝑆𝐾𝑖𝑗𝑞𝑟𝑘−1𝑑𝜀𝑞𝑟𝑓𝑆𝐾𝑘 Ω𝑤𝑝,𝑗𝑑Ω − ∫𝜎𝑖𝑗𝑘−1 Ω𝑤𝑝,𝑗𝑑Ω+ ∫𝑏𝑘−1(𝜓𝐹𝑘−1 − 𝜓𝐹𝑘−2)(𝑃?̂?𝑘−1 − 𝑃𝐹𝑘−1) Ω𝑤𝑝,𝑗𝑑Ω (A-31) The above equation can be written in matrix form as: [      ∆𝑡𝐊𝒖𝒖 𝟎 𝟎 ∆𝑡𝐊𝒖𝑷𝐹 ∆𝑡𝐊𝒖𝑷?̂?𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝟎 𝟎𝟎 𝟎 𝟎 𝟎 𝟎 ]      𝑘−1{      ∆?̅?𝑘?̅?𝐹𝑘?̅??̂?𝑘∆?̅?𝐹𝑘∆?̅??̂?𝑘}      ={      Δ𝑡[𝐟𝒖 + 𝐟𝒖𝑓𝑆𝐾𝑘−𝐟𝒖𝑖𝑛𝑡𝑘−1−𝐟𝒖𝑏𝑃𝑘−1 (∆𝜓?̂?𝑘−1?̅??̂?𝑛𝑘−1+ ∆𝜓𝐹𝑘−1?̅?𝐹𝑛𝑘−1)]𝟎𝟎𝟎𝟎 }       (A-32) where combined with Eq. (A-26) gives the full matrix form of iterative solution for three-phase system as: 170  [      ∆𝑡𝐊𝒖𝒖 𝟎 𝟎 ∆𝑡𝐊𝒖𝑷𝐹 ∆𝑡𝐊𝒖𝑷?̂?𝟎 ∆𝑡𝐊𝒗𝐹𝒗𝐹 𝟎 ∆𝑡𝐊𝒗𝐹𝑷𝐹 𝟎𝟎 𝟎 ∆𝑡𝐊𝒗?̂?𝒗?̂? 𝟎 ∆𝑡𝐊𝒗?̂?𝑷?̂?𝟎 𝟎 𝟎 ∆𝑡𝐊𝑷𝐹𝑷𝐹 ∆𝑡𝐊𝑷𝐹𝑷?̂?𝐂𝑷?̂?𝒖 ∆𝑡𝐊𝑷?̂?𝒗𝐹 ∆𝑡𝐊𝑷?̂?𝒗?̂? ∆𝑡𝑲𝑷?̂?𝑷𝐹 + 𝑪𝑷?̂?𝑷𝐹 ∆𝑡𝑲𝑷?̂?𝑷?̂? + 𝑪𝑷?̂?𝑷?̂?]      𝑘−1{      ∆?̅?𝑘?̅?𝐹𝑘?̅??̂?𝑘∆?̅?𝐹𝑘∆?̅??̂?𝑘}       ={      Δ𝑡[𝐟𝒖 + 𝐟𝒖𝑓𝑆𝐾𝑘−𝐟𝒖𝑖𝑛𝑡𝑘−1−𝐟𝒖𝑏𝑃𝑘−1 (∆𝜓?̂?𝑘−1?̅??̂?𝑛𝑘−1+ ∆𝜓𝐹𝑘−1?̅?𝐹𝑛𝑘−1)]∆𝑡[𝐟𝒗𝐹 −𝐊𝒗𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1]∆𝑡[𝐟𝒗?̂? − 𝐊𝒗?̂?𝑷?̂?𝑘−1?̅??̂?𝑘−1]∆𝑡[𝐟𝑷𝐹 − 𝐊𝑷𝐹𝑷𝐹𝑘−1?̅?𝐹𝑘−1− 𝐊𝑷𝐹𝑷?̂?𝑘−1?̅??̂?𝑘−1]𝐅𝑷?̂? }       (A-33) where 𝐅𝑷?̂? = ∆𝑡[𝐟𝑷?̂? − 𝐊𝑷?̂?𝑷𝐹𝑘−1?̅?𝐹𝑘−1− 𝐊𝑷?̂?𝑷?̂?𝑘−1?̅??̂?𝑘−1] + 𝐂𝑷?̂?𝒖𝑘−1(?̅?0 − ?̅?𝑘−1)+ 𝐂𝑷?̂?𝑷𝐹𝑘−1(?̅?𝐹0− ?̅?𝐹𝑘−1) + 𝐂𝑷?̂?𝑷?̂?𝑘−1(?̅??̂?0− ?̅??̂?𝑘−1) (6-2) and all components are defined in Table 3.2.    

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-0355735/manifest

Comment

Related Items