Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Integrating resin flow and stress development in process modeling of thermoset composites Haghshenas, Seyed Mehdi 2012

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

Item Metadata

Download

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

Full Text

INTEGRATING RESIN FLOW AND STRESS DEVELOPMENT IN PROCESS MODELING OF THERMOSET COMPOSITES by SEYED MEHDI HAGHSHENAS B.Sc. (Civil Engineering), University of Tehran, 2000 M.Sc. (Civil Engineering/Structural Engineering), University of Tehran, 2002  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Civil Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (VANCOUVER) December 2012 © Seyed Mehdi Haghshenas, 2012  Abstract The usual approach to the process modeling of thermoset matrix composites is to divide the analysis into two distinct and sequential steps, first of flow-deformation behaviour and then of stress-deformation. In the current processing models, each of these two aspects is dealt with in a separate sub-model, typically called the flow module and stress module respectively. The flow module is relevant to the pre-gelation behaviour of resin, while the stress module is valid for the post-gelation composite material. In this thesis, the framework to integrate the flow and the stress modules into a unified module in finite element processing models is presented. The work is based on a twophase model for analysis of resin flow and its resulting deformations in the composite material. Special measures are introduced to provide for additional capability of this model to account for the development of stresses in the curing composite material. These modifications are needed to ensure the accuracy of the model in both of resin flow and stress development regimes, and include the introduction of consistent compressibility in the mass conservation equation of the two-phase system, and a special decomposition of stresses of the system. The formulation is implemented for a pseudo-viscoelastic stress model in a 2D plane strain FE code in MATLAB. The approach may readily be extended to fully viscoelastic models. Various examples from single-element problems dealing with the development of residual stresses throughout a single-hold cure cycle to more geometrically complex composite laminates undergoing standard cure cycles are modeled by the integrated model and comparisons are made in one extreme to the flow-compaction behaviour by the standard flow models, and in the other extreme to the results obtained by the pseudoviscoelastic approach. The model developed here is a promising tool for simulating processing of large-scale composite structures continuously from the very early stages of the process when the resin behaves in a fluid-like manner all the way to the final stage when it behaves as a 3D solid. ii  Table of Contents  Abstract .............................................................................................................................. ii Table of Contents ............................................................................................................. iii List of Tables ................................................................................................................... vii List of Figures ................................................................................................................. viii List of Symbols ............................................................................................................. xviii Acknowledgements ....................................................................................................... xxii Chapter 1 : Introduction .................................................................................................. 1 1.1. Background ............................................................................................................................. 1 1.2. Resin Flow during Processing ............................................................................................... 4 1.3. Stress Development during Cure .......................................................................................... 6 1.4. Motivation ............................................................................................................................... 8 1.5. Research Objectives and Thesis Outline .............................................................................. 9 1.6. Figures ................................................................................................................................... 12  Chapter 2 : Governing Equations and FE Model of Two-phase Media .................... 15 2.1. Introduction .......................................................................................................................... 15 2.2. Volume Averaging Method.................................................................................................. 18  iii  Table of Contents 2.3. Governing Equations............................................................................................................ 20 2.4. Finite Element Implementation........................................................................................... 29 2.4.1. Q2P-1 element ..................................................................................................................................29 2.4.2. Galerkin finite element equations ...................................................................................................30 2.4.3. Q1P0 element ...................................................................................................................................41 2.4.4. Boundary conditions .......................................................................................................................42 2.4.5. Time integration..............................................................................................................................44 2.4.6. Non-linear solution scheme ............................................................................................................45  2.5. Treatment of Discontinuity in Volume Fraction ............................................................... 46 2.6. Interfaces with Purely Viscous Fluids ................................................................................ 50 2.7. Figures ................................................................................................................................... 53  Chapter 3 : Verification of the Two-phase Model ....................................................... 55 3.1. General Two-phase Problems ............................................................................................. 55 3.1.1. Column of saturated porous medium under loading at the top .......................................................55 3.1.2. Column of saturated porous medium under loading at mid-height .................................................57 3.1.3. Taylor and Beavers-Joseph problems .............................................................................................58  3.2. Flow-compaction in Processing of Composite Angle Laminates ...................................... 62 3.3. Tables..................................................................................................................................... 65 3.4. Figures ................................................................................................................................... 66  Chapter 4 : Integration of Modeling from Fluid to Cured Resin ............................... 84 4.1. Pseudo-viscoelastic Stress Model ........................................................................................ 84 4.2. Constitutive Equations ......................................................................................................... 86 4.2.1. Isotropic materials ..........................................................................................................................86 4.2.2. Transversely isotropic materials .....................................................................................................87  4.3. Two-phase Model for Stress Development in Cured Materials ....................................... 89 4.3.1. Consistent compressibility in mass conservation equation .............................................................90 iv  Table of Contents 4.3.2. The concept of modified effective stress („s-p‟ stress) ...................................................................94 4.3.3. Elastic moduli for isotropic materials .............................................................................................97 4.3.4. Consistency of constitutive equations for isotropic materials .........................................................98 4.3.5. Elastic moduli for transversely isotropic materials .......................................................................101 4.3.6. Consistency of constitutive equations for transversely isotropic materials...................................104 4.3.7. Decomposition of stresses in 3D...................................................................................................109  4.4. Modifications in the Numerical Technique ...................................................................... 111 4.4.1. FE implementation in 2D plane strain ..........................................................................................112 4.4.2. Modifications in the non-linear solution scheme ..........................................................................115 4.4.3. Calculation of volume fractions ....................................................................................................115  4.5. Resin Properties during Cure............................................................................................ 117 4.6. Tables................................................................................................................................... 118  Chapter 5 : Numerical Applications ........................................................................... 119 5.1. Cured Materials .................................................................................................................. 119 5.1.1. Cured composite sample under uni-axial strain conditions ..........................................................120 5.1.2. Cured composite sample under linear distribution of pressure .....................................................121 5.1.3. Cured Composite sample under pressure gradient ........................................................................123 5.1.4. Transversely isotropic samples .....................................................................................................126  5.2. Materials Undergoing Cure ............................................................................................... 129 5.2.1. One-element examples ..................................................................................................................129 5.2.2. Flat unidirectional laminate undergoing cure ...............................................................................134 5.2.3. Flow-compaction and stress development in angle laminates ......................................................137  5.3. Tables................................................................................................................................... 144 5.4. Figures ................................................................................................................................. 148  Chapter 6 : Conclusions and Further Work .............................................................. 204 6.1. Summary and Conclusions ................................................................................................ 204 6.2. Contributions ...................................................................................................................... 205  v  Table of Contents 6.3. Further Work ..................................................................................................................... 207  Bibliography .................................................................................................................. 209  vi  List of Tables  Table 3.1: Resin and fibre-bed properties for the AS4/3501-6 angle laminate ................ 65 Table 3.2: Resin and fibre-bed properties for the AS4/8552 angle laminate.................... 65 Table 4.1: Relaxation times and weight factors for 3501-6 resin for α0 = 0.98 .............. 118 Table 5.1: The comparison and convergence of strain energy values (in N.mm) for the isotropic sample shown in Figure 5.2, using 9-3 elements ...................................... 144 Table 5.2: The comparison and convergence of uA (×10-7 m) for the isotropic sample shown in Figure 5.4, using 9-3 elements ................................................................. 144 Table 5.3: The comparison and convergence of strain energy values (in N.mm) for the isotropic sample shown in Figure 5.4, using 9-3 elements ...................................... 145 Table 5.4: The comparison and convergence of uA (×10-7 m) for the isotropic sample shown in Figure 5.4, using 9-4 elements ................................................................. 145 Table 5.5: The comparison and convergence of strain energy values (in N.mm) for the isotropic sample shown in Figure 5.4, using 9-4 elements ...................................... 145 Table 5.6: The comparison and convergence of strain energy and uA (×10-7 m) for the transversely isotropic sample (   0 ) shown in Figure 5.4, using 9-3 elements .. 146 Table 5.7: The comparison and convergence of strain energy and uA (×10-7 m) for the transversely isotropic sample (   90 ) shown in Figure 5.4, using 9-3 elements 146 Table 5.8: The comparison and convergence of strain energy and uA (×10-7 m) for the transversely isotropic sample (   0 ) shown in Figure 5.4, using 9-3 elements .. 147 Table 5.9: Resin properties in Examples 1 & 2, where Poisson‟s ratio is assumed constant ................................................................................................................................. 147 Table 5.10: Mechanical properties of fibre and resin in the examples of Section 5.2 .... 147  vii  List of Figures  Figure 1.1: Autoclave processing steps, (a) lay-up, (b) cure with a controller loop (adapted from Hubert, 1996) ..................................................................................... 12 Figure 1.2: Integrated sub-model approach in process modeling (adapted from Johnston et al., 2001) ................................................................................................................ 13 Figure 1.3: Different regimes of resin behaviour during a typical autoclave processing history ........................................................................................................................ 13 Figure 1.4: Time-history of changes in resin viscosity during a typical autoclave processing history ...................................................................................................... 14 Figure 1.5: Schematic representation of resin flow through the fibre-bed ....................... 14 Figure 2.1: Schematic representation of Q2P-1 (9-3) element ........................................... 53 Figure 2.2: Schematic representation of Q1P0 (4-1) element ............................................ 53 Figure 2.3: The schematic diagram of the pressure and volume fraction distribution for a 1-D system consisting of 2 elements with a jump at their interface .......................... 54 Figure 3.1: A column of saturated porous medium under uniform loading at the top, and meshes of 10 9-3 elements and 20 9-4 elements ....................................................... 66 Figure 3.2: A column of saturated porous medium with impermeable BC all around with uniform loading at mid-height, and mesh of 10 9-3 elements ................................... 66 Figure 3.3: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=5×103Pa.s ................................................................................................................................... 67 Figure 3.4: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=5×103 Pa.s ................................................................................................................................... 67 Figure 3.5: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=5×102 Pa.s ............................................................................................................................. 68  viii  List of Figures  Figure 3.6: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=5×102 Pa.s ................................................................................................................................... 68 Figure 3.7: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=50 Pa.s 69 Figure 3.8: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=50 Pa.s ... 69 Figure 3.9: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=5 Pa.s .. 70 Figure 3.10: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=5 Pa.s ..... 70 Figure 3.11: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=0.5 Pa.s 71 Figure 3.12: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=0.5 Pa.s .. 71 Figure 3.13: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=5×103 Pa.s .............................................................. 72 Figure 3.14: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=5×102 Pa.s .............................................................. 72 Figure 3.15: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=50 Pa.s.................................................................... 73 Figure 3.16: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=5 Pa.s...................................................................... 73 Figure 3.17: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=0.5 Pa.s................................................................... 74 Figure 3.18: Pressure surfaces of the two end elements at a very early stage of consolidation of the column ...................................................................................... 74 Figure 3.19: Pressure surfaces of the two end elements at a time that oscillation is wellpronounced ................................................................................................................ 75  ix  List of Figures  Figure 3.20: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained without the correction for discontinuity of volume fraction. μ=5×10-2 Pa.s ....................................................... 75 Figure 3.21: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained without the correction for discontinuity of volume fraction. μ=5×10-3 Pa.s ....................................................... 76 Figure 3.22: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained with the correction for discontinuity of volume fraction. μ=5×10-2 Pa.s ....................................................... 76 Figure 3.23: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained with the correction for discontinuity of volume fraction. μ=5×10-3 Pa.s ....................................................... 77 Figure 3.24: Schematic diagram of two unidirectional flow problems over and through a porous wall; (a) Taylor problem, (b) Beavers-Joseph problem ................................. 77 Figure 3.25: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10-6 m2 ..................................................................................................... 78 Figure 3.26: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10-8 m2 ..................................................................................................... 78 Figure 3.27: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10-10 m2 .................................................................................................... 79 Figure 3.28: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10-6 m2 ..................................................................................... 79 Figure 3.29: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10-8 m2 ..................................................................................... 80 Figure 3.30: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10-10 m2 .................................................................................... 80 Figure 3.31: Geometry, BC, and processing cycle used for compaction of a unidirectional angle laminate on a convex tool ................................................................................ 81 Figure 3.32: (a) 8×4 mesh and (b) 16×6 mesh of the angle laminate using symmetry .... 81 Figure 3.33: Time-history of normal displacement for unidirectional AS4/3501-6 angle with [0°] fibres on a convex tool ............................................................................... 82 x  List of Figures  Figure 3.34: Time-history of normal displacement for unidirectional AS4/8552 angle with [0°] fibres on a convex tool ............................................................................... 82 Figure 3.35: Time-history of normal displacement for unidirectional AS4/3501-6 angle with [90°] fibres on a convex tool ............................................................................. 83 Figure 3.36: Time-history of normal displacement for unidirectional AS4/8552 angle with [90°] fibres on a convex tool ............................................................................. 83 Figure 5.1: The geometry and specifications of the composite sample under constant pressure load (uni-axial strain condition) ................................................................ 148 Figure 5.2: The geometry and specifications of the composite sample under linear pressure load ............................................................................................................ 148 Figure 5.3: The profile of the vertical displacement at the top of the constrained sample under linear distribution of pressure load ................................................................ 149 Figure 5.4: Composite sample under pressure gradient .................................................. 149 Figure 5.5: Shape of the displacement profiles of the sample under pressure gradient.. 150 Figure 5.6: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for different values of CoeffK, obtained using a 3×2 mesh of 9-3 elements ..... 150 Figure 5.7: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for different values of CoeffK, obtained using a 6×4 mesh of 9-3 elements ..... 151 Figure 5.8: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for different values of CoeffK, obtained using a 12×8 mesh of 9-3 elements ... 151 Figure 5.9: Comparison of displacement profiles at the right side of the sample in Figure 5.4 using different meshes of 9-3 elements, in the case where CoeffK=0.01 .......... 152 Figure 5.10: Comparison of displacement profiles at the right side of the sample in Figure 5.4 using different meshes of 9-3 elements, in the case where CoeffK=1 ............... 152 Figure 5.11: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for 9-3 and 9-4 elements, obtained using a 6×4 mesh ....................................... 153 Figure 5.12: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for 9-3 and 9-4 elements, obtained using a 12×8 mesh ..................................... 153 Figure 5.13: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-3 elements ............................................................................................................. 154  xi  List of Figures  Figure 5.14: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-4 elements with permeable BC on the sides ......................................................... 154 Figure 5.15: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-4 elements with impermeable BC on the sides .................................................... 155 Figure 5.16: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-4 elements based on the u-v-p formulation .......................................................... 155 Figure 5.17: Time-history of the applied temperature and the resulting degree of cure in Examples 1-7 ........................................................................................................... 156 Figure 5.18: Time-history of the applied temperature and the resulting resin viscosity in Examples 1-7 ........................................................................................................... 156 Figure 5.19: Time-history of the development of longitudinal stress in Example 1. The predictions are essentially identical. ........................................................................ 157 Figure 5.20: Time-history of the development of transverse stress in Examples 1 and 2. The predictions are essentially identical. ................................................................. 157 Figure 5.21: Time-history of the development of  x in Example 2. The predictions are essentially identical.................................................................................................. 158 Figure 5.22: Time-history of the development of  x in Example 3. The predictions are essentially identical.................................................................................................. 158 Figure 5.23: Time-history of the development of  z in Examples 3 and 4. The predictions are essentially identical. ........................................................................ 159 Figure 5.24: Time-history of the development of resin pressure in Examples 3 and 4 .. 159 Figure 5.25: Time-history of the development of  x in Example 4. The predictions are essentially identical.................................................................................................. 160 Figure 5.26: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 5 ......................................................................................... 160 Figure 5.27: Time-history of the applied temperature and the changes in the rate of volumetric strain during t=15-55 min in Example 5 ............................................... 161 Figure 5.28: Time-history of the changes in resin volume fraction in Example 5 ......... 161 Figure 5.29: Time-history of the development of resin pressure in Example 5 ............. 162 Figure 5.30: Time-history of the development of  x in Example 5............................... 162  xii  List of Figures  Figure 5.31: Time-history of the development of  z in Example 5............................... 163 Figure 5.32: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 6 ......................................................................................... 163 Figure 5.33: Time-history of the changes in resin volume fraction in Example 6 ......... 164 Figure 5.34: Time-history of the development of resin pressure in Example 6 ............. 164 Figure 5.35: Time-history of the transverse strain in Example 6 ................................... 165 Figure 5.36: Time-history of the development of  x in Example 6............................... 165 Figure 5.37: Time-history of the development of resin pressure in Example 7 ............. 166 Figure 5.38: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 7 ......................................................................................... 166 Figure 5.39: Time-history of the changes in resin volume fraction in Example 7 ......... 167 Figure 5.40: Time-history of the transverse strain in Example 7 ................................... 167 Figure 5.41: Time-history of the development of  x in Example 7............................... 168 Figure 5.42: Geometry and BC of the flat laminate undergoing cure ............................ 168 Figure 5.43: Time-history of the autoclave and part temperatures and the resulting degree of cure in the flat and angle laminate examples ...................................................... 169 Figure 5.44: Final profiles of σx predicted by the stress model at section AB of the [0°] and [90°] flat laminates, and their comparison with response obtained from ABAQUS ................................................................................................................. 169 Figure 5.45: Final distribution of σx predicted by the stress model along the bottom surface of the [90°] flat laminates in comparison with response obtained from ABAQUS ................................................................................................................. 170 Figure 5.46: Comparison of final profiles of resin volume fraction predicted by the integrated model at AB and CD sections for different lay-ups of the flat laminate 170 Figure 5.47: Final distribution of σx for the [0°] flat laminate obtained by the integrated model ....................................................................................................................... 171 Figure 5.48: Final distribution of σx for the [45°] flat laminate obtained by the integrated model ....................................................................................................................... 171 Figure 5.49: Final distribution of σx for the [90°] flat laminate obtained by the integrated model ....................................................................................................................... 171 Figure 5.50: Final profile of σx at section AB of the [0°] flat laminate .......................... 172 xiii  List of Figures  Figure 5.51: Final profile of σx at section AB of the [45°] flat laminate ........................ 172 Figure 5.52: Final profile of σx at section AB of the [90°] flat laminate ........................ 173 Figure 5.53: Final distribution of axial force along the length of the [0°] flat laminate. 173 Figure 5.54: Final distribution of bending moment along the length of the [0°] flat laminate ................................................................................................................... 174 Figure 5.55: Final distribution of axial force along the length of the [45°] flat laminate174 Figure 5.56: Final distribution of bending moment along the length of the [45°] flat laminate ................................................................................................................... 175 Figure 5.57: Final distribution of axial force along the length of the [90°] flat laminate175 Figure 5.58: Final distribution of bending moment along the length of the [90°] flat laminate ................................................................................................................... 176 Figure 5.59: Geometry and boundary conditions of the convex angle laminate ............ 176 Figure 5.60: Geometry and boundary conditions of the concave angle laminate ........... 177 Figure 5.61: Location of key points on the geometry of the angle laminate, and transformation of geometry into the ξη coordinates ................................................ 177 Figure 5.62: FE meshes used in the analysis of angle laminates .................................... 178 Figure 5.63: Final profile of resin volume fraction for the convex [0°] angle laminate at section EF obtained by the integrated model with different meshes ....................... 179 Figure 5.64: Distribution of resin volume fraction for the convex [0°] angle laminate obtained by the integrated model at t = 14 min ....................................................... 179 Figure 5.65: Distribution of resin volume fraction for the convex [0°] angle laminate obtained by the integrated model at t = 44 min ....................................................... 180 Figure 5.66: Final distribution of resin volume fraction for the convex [0°] angle laminate obtained by the integrated model............................................................................. 180 Figure 5.67: Final distribution of resin volume fraction for the convex [45°] angle laminate obtained by the integrated model .............................................................. 181 Figure 5.68: Final distribution of resin volume fraction for the convex [90°] angle laminate obtained by the integrated model .............................................................. 181 Figure 5.69: Final distribution of tangential stress for the convex [0°] angle laminate obtained by the integrated model............................................................................. 182  xiv  List of Figures  Figure 5.70: Final distribution of tangential stress for the convex [45°] angle laminate obtained by the integrated model............................................................................. 182 Figure 5.71: Final distribution of tangential stress for the convex [45°] angle laminate obtained by the integrated model............................................................................. 183 Figure 5.72: Final distribution of tangential stress for the convex [0°] angle laminate obtained by the stress model .................................................................................... 183 Figure 5.73: Final distribution of tangential stress for the convex [45°] angle laminate obtained by the stress model .................................................................................... 184 Figure 5.74: Final distribution of tangential stress for the convex [90°] angle laminate obtained by the stress model .................................................................................... 184 Figure 5.75: Distribution of tangential stress for the convex [0°] angle laminate obtained by the integrated model at t = 144 min .................................................................... 185 Figure 5.76: Distribution of tangential stress for the convex [0°] angle laminate obtained by the stress model at t = 144 min ........................................................................... 185 Figure 5.77: Time-history of the development of tangential stress at different locations of the convex [0°] angle laminate ................................................................................ 186 Figure 5.78: Time-history of the development of tangential stress at point C on the convex [0°] angle laminate ...................................................................................... 186 Figure 5.79: Time-history of the development of tangential stress at point D on the convex [0°] angle laminate ...................................................................................... 187 Figure 5.80: Time-history of the development of tangential stress at point A on the convex [0°] angle laminate ...................................................................................... 187 Figure 5.81: Time-history of the development of tangential stress at point B on the convex [0°] angle laminate ...................................................................................... 188 Figure 5.82: Time-history of the development of tangential stress at different locations of the convex [45°] angle laminate .............................................................................. 188 Figure 5.83: Time-history of the development of tangential stress at point C on the convex [45°] angle laminate .................................................................................... 189 Figure 5.84: Time-history of the development of tangential stress at point D on the convex [45°] angle laminate .................................................................................... 189  xv  List of Figures  Figure 5.85: Time-history of the development of tangential stress at point A on the convex [45°] angle laminate .................................................................................... 190 Figure 5.86: Time-history of the development of tangential stress at point B on the convex [45°] angle laminate .................................................................................... 190 Figure 5.87: Time-history of the development of tangential stress at different locations of the convex [90°] angle laminate .............................................................................. 191 Figure 5.88: Time-history of the development of tangential stress at point C on the convex [90°] angle laminate .................................................................................... 191 Figure 5.89: Time-history of the development of tangential stress at point D on the convex [90°] angle laminate .................................................................................... 192 Figure 5.90: Time-history of the development of tangential stress at point A on the convex [90°] angle laminate .................................................................................... 192 Figure 5.91: Time-history of the development of tangential stress at point B on the convex [90°] angle laminate .................................................................................... 193 Figure 5.92: Time-history of the development of tangential stress at point F on the convex [90°] angle laminate .................................................................................... 193 Figure 5.93: Final distribution of tangential stress for the concave [0°] angle laminate obtained by the integrated model............................................................................. 194 Figure 5.94: Final distribution of tangential stress for the concave [45°] angle laminate obtained by the integrated model............................................................................. 194 Figure 5.95: Final distribution of tangential stress for the concave [90°] angle laminate obtained by the integrated model............................................................................. 195 Figure 5.96: Time-history of the development of tangential stress at different locations of the concave [0°] angle laminate............................................................................... 195 Figure 5.97: Time-history of the development of tangential stress at point C on the concave [0°] angle laminate .................................................................................... 196 Figure 5.98: Time-history of the development of tangential stress at point D on the concave [0°] angle laminate .................................................................................... 196 Figure 5.99: Time-history of the development of tangential stress at point A on the concave [0°] angle laminate .................................................................................... 197  xvi  List of Figures  Figure 5.100: Time-history of the development of tangential stress at point B on the concave [0°] angle laminate .................................................................................... 197 Figure 5.101: Final distribution of axial force along the length of the convex [0°] angle laminate ................................................................................................................... 198 Figure 5.102: Final distribution of bending moment along the length of the convex [0°] angle laminate .......................................................................................................... 198 Figure 5.103: Final distribution of axial force along the length of the convex [45°] angle laminate ................................................................................................................... 199 Figure 5.104: Final distribution of bending moment along the length of the convex [45°] angle laminate .......................................................................................................... 199 Figure 5.105: Final distribution of axial force along the length of the convex [90°] angle laminate ................................................................................................................... 200 Figure 5.106: Final distribution of bending moment along the length of the convex [90°] angle laminate .......................................................................................................... 200 Figure 5.107: Final distribution of axial force along the length of the concave [0°] angle laminate ................................................................................................................... 201 Figure 5.108: Final distribution of bending moment along the length of the concave [0°] angle laminate .......................................................................................................... 201 Figure 5.109: Final distribution of axial force along the length of the concave [45°] angle laminate ................................................................................................................... 202 Figure 5.110: Final distribution of bending moment along the length of the concave [45°] angle laminate .......................................................................................................... 202 Figure 5.111: Final distribution of axial force along the length of the concave [90°] angle laminate ................................................................................................................... 203 Figure 5.112: Final distribution of bending moment along the length of the concave [90°] angle laminate .......................................................................................................... 203  xvii  List of Symbols  aT  Shift factor  b  Biot coefficient  b1  Biot coefficient in the fibre direction  C  Matrix of coefficients of main variables‟ rate in FE  Cijkl  Stiffness tensor  CTE  Coefficient of thermal expansion  CTE1  Longitudinal coefficient of thermal expansion  CTE2  Transverse coefficient of thermal expansion  CTEf  Fibre coefficient of thermal expansion  CTEg  Glassy coefficient of thermal expansion  CTEm  Resin coefficient of thermal expansion  CTEr  Rubbery coefficient of thermal expansion  CSC  Volumetric cure shrinkage coefficient  D  Material stiffness matrix  E  Young‟s modulus/ Transverse modulus in transversely isotropic media  E1  Longitudinal modulus in transversely isotropic media  F  Vector of applied loads  fd i  Vector of drag force  G12  In-plane shear modulus  G23  Transverse shear modulus  Gc  Shear modulus of the composite material xviii  List of Symbols  Gfb  Shear modulus of the fibre-bed  Gm  Shear modulus of resin  Gr  Relaxed shear modulus of resin  Gu  Un-relaxed shear modulus of resin  I  Unit tensor  J  Jacobian matrix  K  Matrix of coefficients of main variables in FE  K2  In-plane bulk modulus in the plane of isotropy  Kc  Bulk modulus of the composite material  Kf  Bulk modulus of fibres  Kfb  Bulk modulus of fibre-bed  Km  Bulk modulus of resin  Kr  Relaxed bulk modulus of resin  Ku  Un-relaxed bulk modulus of resin  m  Cool-down rate in a cure cycle  Ni  Displacement shape function for the ith node  N ip  Pressure shape function for the ith pressure node  N i  Shape function of fluid volume fraction for the ith Gauss point  P  Resin pressure  p  Vector of pressure DOFs  te  Variable time, used in the pseudo-viscoelastic model  T  Temperature  ui  Displacement vector of the system  u  Vector of displacement DOFs  xix  List of Symbols  vi  Vector of resin relative velocity  v  Vector of relative velocity DOFs  wp  Weight functions in FE  W  Weight factor related to the ωth Maxwell element    Resin‟s degree of cure  0  Resin‟s reference degree of cure    Boundary of the integration domain   ij  Kronecker delta  t  Time-step size    Vector of strain rates  vs  Volumetric strain of the system  v s  Rate of volumetric strain of the system  v s mech  Rate of mechanical volumetric strain of the system  v s free  Rate of free volumetric strain of the system    A material property in transversely isotropic media    Viscosity    Poisson‟s ratio    Reduced time  f  Density of fibres  m  Density of resin  ζ ,  ij  Stress tensor  ζ ,  ij  Effective stress tensor  xx  List of Symbols  ζ ,  ij  Biot effective stress tensor    Vector of stress rates   s ij  Total stress tensor of the system   s p ij  „s-p‟ effective stress tensor    Dummy time variable (for time integration)   ij  Shear stress tensor   m ij  Shear stress tensor of the resin  p  Peak relaxation time    Relaxation time of the ωth Maxwell element    Volume fraction of resin  0  Initial volume fraction of resin    Volume fraction of α-phase    Number representing any of the Maxwell elements    Domain of integration    A material property in transversely isotropic media  B  Gradient of B  A  Divergence of vector A  B  Spatial average of variable B   B   Phase average of variable B over α-phase  B    Intrinsic phase average of variable B over α-phase  xxi  Acknowledgements I would like to take this opportunity to acknowledge the help of those without whom this work would have not come to fruition. Firstly, I would like to thank my supervisors Dr. Reza Vaziri and Dr. Anoush Poursartip for their invaluable guidance, expertise, and patience. Also, the advice and encouragement generously offered by Dr. Göran Fernlund and Dr. Fernand Ellyin is truly acknowledged. Many thanks to several past and present members of UBC Composites Group for their friendship and support. Especially, I would like to thank Dr. Abdul R. Arafath for his invaluable help at different stages of my work. I am grateful to Mr. Robert Courdji for his great help throughout this work. Many fruitful discussions with Dr. Amir Osooly, Dr. Alireza Forghani, Dr. Nima Zobeiry, and Dr. Navid Zobeiry are greatly appreciated. I truly appreciate the companionship of my close friends, and I appreciate their enthusiasm, support and encouragement through the course of my degree. Financial support from the University of British Columbia and Killam Trusts is truly appreciated. My special gratitude is due to my mother, sister and two brothers for their continuous love, encouragement, and support. I am deeply grateful to Samaneh for her patience, encouragement, and sacrifices throughout the past five years. This dissertation is dedicated to my mother and the memory of my father both of whom encouraged, inspired, and celebrated their children‟s achievements with exceptional enthusiasm.  xxii  Chapter 1: Introduction  1.1. Background Fibre-reinforced plastic composites are replacing traditional materials in engineering applications. The growing popularity of these materials is due to advantages such as high strength to weight ratio, high stiffness to weight ratio, and impressive durability. Manufacturing of composite structures is distinct from that of metallic structures in that a complex and large-scale structure is processed from the raw materials in one step in order to create otherwise unachievable geometries and to reduce manufacturing cost. The focus of this work is on the high-performance structural components made of advanced thermoset matrix composites that are mainly used in the aerospace industry. The making of such structures is commonly performed by autoclave processing. This process involves stacking of pre-impregnated sheets of unidirectional fibres (prepreg) at specially designed orientations over a tool of desired shape. The whole assembly is then subjected to a controlled cycle of temperature and pressure change inside an autoclave. The end result of this process is the compaction and curing of the composite part. Finally, the cured composite part is removed from the tool However, development of residual stresses during the curing process leads to difficulty and uncertainty in prediction of the precise shape and dimensions of the part after tool removal. It is essential to have a good understanding of the phenomena that the material undergoes throughout the processing cycle. In the past three decades as the popularity of composite materials has increased, efforts to numerically simulate the processing cycle of composite materials have increasingly gained more traction in order to curb the costs of production. A bagging schematic for a typical autoclave processing system is shown in Figure 1.1 (a). The prepreg plies are placed on a hard tool to form the desired shape of the laminate. A typical prepreg ply has a thickness of 0.2 mm and a fibre volume fraction of 0.5-0.6. The plies can be oriented in different directions to achieve the target mechanical properties.  1  Chapter 1: Introduction  Typical tooling materials include aluminum, steel, invar, and even composite. Dams are placed around the perimeter of the laminate to control resin flow from the edges of the plies. In many cases, a release film is used on top of the laminate to allow for easy removal of adjacent bagging material. For laminates that need resin bleed, an absorbing cloth forming the bleeder is placed on top and around the laminate. Inserts and honeycomb cores may be included in the laminate for moulding purposes or structural requirements. A breather cloth covers the assembly to provide a path for air flow. The whole assembly consisting of the tool, laminate, bleeder, and breather is placed in a vacuum bag and is sealed. A vacuum plug connects the interior of the bag to an external vacuum pump. The cycle of curing the part is shown in Figure 1.1 (b). At this stage, the tool-laminate assembly is put in an autoclave and the sealed bag is attached to the vacuum system. The vacuum is applied in order to debulk the composite laminate, partially motivate the flow of excess resin out of the laminate, and remove any excess air that is trapped between the plies. In order to cure the part, pressure and temperature are applied in a pre-determined cure cycle. The temperature cycle is crucial to trigger the reaction of resin polymerization. The pressure is applied to conform the laminate to the tool surface, and to compact the laminate to the target value for fibre volume fraction and collapse any voids that may develop and grow during resin cure. At the end of the processing, the cured part is debagged and removed from the rest of the assembly and the tool and is ready for any finishing processes. In an ideal curing process, the cure cycle, the tool, and the bagging procedure should be designed in a way that leads to a fully cured, void free and undistorted final part in the shortest time and most economical fashion. The processing of composite materials is a complex event during which various physical and chemical changes occur simultaneously or at different stages. These phenomena include melting of resin, flow of resin through the fibre-bed, heat transfer, thermochemical changes and curing of resin, and development of residual stresses in the composite material. This complexity has led many researchers to use an „integrated submodel‟ to model the processing of composites. Based on this approach, the complex process model is divided into several sub-models that may be studied quite independently. The sub-model approach is shown schematically in Figure 1.2. As shown 2  Chapter 1: Introduction  in the schematics of this approach, the modeling procedure includes several modules each responsible for one aspect of the material behaviour during cure. The modules that are of interest in this work are the stress module and the flow module. In chronological order some of the notable work on the application of this method to autoclave process modeling were done by Loos and Springer (1983), Bogetti and Gillespie (1991, 1992), White and Hahn (1992a, 1992b), Johnston et al. (2001), and Zhu et al. (2001). Zobeiry (2006) presents a more comprehensive review of the relevant literature. The processing models range from very simple one-dimensional elastic analyses to very complicated 3D viscoelastic models. One of these models was developed at UBC in the form of a multi-physics, 2D finite element code, COMPRO, to analyze industrial autoclave processing of composite materials of intermediate size and complexity [Hubert (1996), Hubert et al. (1999), Johnston (1997), Johnston et al. (2001)]. The model predicts a number of essential processing parameters and the development of residual stresses and deformations. It also accounts for the effects of tool/part interaction in the development of residual stresses. COMPRO has been used in the past two decades to solve many practical problems in autoclave processing faced by the industry. To expand on the capabilities of COMPRO to different material constitutive behaviour and more importantly 3D modeling, Arafath (2007) introduced the concept of COMPRO Component Architecture (CCA) and incorporated it into ABAQUS. CCA is a modular approach that consists of all the different material subroutines in COMPRO that can be incorporated into many commercial finite element codes using a specific interface module. In this work, we refer to the CCA approach as COMPRO 3D. It is very common that the literature on modeling of autoclave processing focus either on resin flow during processing or on stress development throughout cure and residual deformations in the final part. The same was true in the development of COMPRO as Hubert (1996) developed the flow module while Johnston‟s (1997) focus was on the stress development side. Even the current version of COMPRO 3D approaches the modeling of autoclave processing in two distinct steps. Initially, the flow of resin and the compaction of composite laminate are modeled using the flow module in COMPRO. Once the flow analysis is performed, the geometry and fibre volume fraction of the 3  Chapter 1: Introduction  laminate are updated. Then the laminate is re-meshed in FE and the new mesh is used in the prediction of development of stresses during the processing cycle. The reason for this separation could be that resin flow is typically formulated by Darcy‟s law and is modeled by two-phase elements in a FE code, while modeling of stress development requires the use of common solid finite elements.  1.2. Resin Flow during Processing Curing of a polymer resin may be defined as the gradual formation of a 3D network made out of the polymer molecules. During the processing, resin transforms from a vitrified fluid to a viscous fluid (upon heating up from the room temperature) and eventually to a viscoelastic material after gelation. Pre-gelation resin may be considered a viscous fluid with its viscosity a function of temperature and degree of cure. Therefore flow phenomena play a major role in defining the behaviour of resin and the composite material at early stages of processing (See Figure 1.3). Figure 1.4 shows how the resin viscosity changes in typical cure cycle. It is evident form the figure that lower viscosity values of resin occur during the first isothermal hold in the autoclave temperature history. In fact the first hold in temperature is typically designed to provide enough time for fluid resin to flow and also allow for air transport as a result. Figure 1.5 shows a schematic presentation of flow of resin through the porous structure consisting of the fibres and its resulting compaction in the composite laminate. Hubert and Poursartip (1998) did a comprehensive review of the modelling approaches in the literature regarding flow and compaction in the processing of thermoset matrix laminates. Springer (1982) studied the relationship between the applied pressure and resin flow and observed that the layers consolidated in a wavelike manner. He developed the sequential compaction model based on those experimental results. Loos and Springer (1983) developed a one-dimensional resin flow model for the curing process. The resin velocity was related to the pressure gradient, fibre-bed permeability, and resin viscosity through Darcy‟s law. Gutowski et al. (1987) presented 3D flow and 1D consolidation models for the processing of composites, where the resin flow was modeled using 4  Chapter 1: Introduction  Darcy‟s law for an anisotropic porous medium. Their consolidation model was based on the application of the effective stress theory of soil mechanics [Terzaghi(1943), Biot(1941)] to the compaction of fibre-beds. Davé (1990) derived the general flow equation applicable to all composite processes which accounted for the presence of both resin and air in the fibre-bed. Cai and Gutowski (1992) developed a general 3D deformation model of lubricated fibre bundle subjected to a multidirectional stress state. Smith and Poursartip (1992) showed that the sequential compaction model is an special case of the effective stress formulation. Gutowski and Dillon (1997) reviewed the stresses of an aligned fibre bundle and compared the response for transverse compression coupled with axial extension in a 3D fibre deformation model. Hubert (1996) and Hubert et al. (1999) developed a 2D plane strain flow-compaction FE model as part of COMPRO to simulate resin flow in autoclave processing of fibrereinforced composite laminates. They adapted an incremental, quasi-linear elastic model assuming infinitesimal strains. Hubert and Poursartip (2001) presented a method the fibre-bed compaction curve directly from composite prepregs. Li and Tucker (2002) developed a two-phase continuum based process model assuming hyperelasticity for the fibre-bed stress, which is valid for large deformations. The advantage of their work over the model of Hubert et al. (1999) is that the large deformation formulation allows them to update the mesh geometry and fibre orientation during the course of process. This is especially important for composites with lower viscosity resins such as AS4/3501-6. Li and Tucker (2002) stated that their formulation has another advantage in that it accounts for spatial variation of fluid volume fraction while that of Hubert(1996) and Hubert et al. (1999) does not. This is not correct as Hubert‟s (1996) simplification of mass conservation equation includes an error that has lead him to unnecessarily assume that the spatial variation of volume fraction be negligible. As a matter of fact, the two models practically use the same mass conservation equation, and therefore both account for spatial variations in volume fraction. Larsson et al. (2004) presented a biphasic continuum model with large deformation capabilities for resin flow and deformation in a family of forming processes for fibrereinforced composites. Liquid composite moulding (LCM) processes such as resin 5  Chapter 1: Introduction  transfer moulding (RTM), vacuum assisted RTM (VARTM), and injection/compression liquid composite moulding (I/C-LCM) also involve resin flow phenomena. These include infusion of resin through the dry fibrous reinforcement, and also flow-compaction of the system under vacuum or applied pressure. Some relevant work on modeling the resin flow phenomena in LCM processes include Bruschke and Advani (1990), Gebart (1992), Trochu et al. (1993, 2006), Tucker and Dessenberger (1994), Pillai and Advani (1998), Pillai et al. (2000,2001), Pillai (2002), Bréard et al. (2003a, 2003b), and Tan and Pillai (2012a, 2012b) in chronological order.  1.3. Stress Development during Cure At gelation, a 3D network of molecules has formed all over the body of thermoset resin. This network causes the resin to attain elastic properties, and therefore from this point on it can develop residual stresses. By the progress of cure, the network of molecules gets more intricate and strong and the elastic moduli of the resin grow exponentially. This does not contradict the general understanding that the mechanical behaviour of polymer resins at any stage of cure is best represented by viscoelastic behaviour as viscoelastic formulation includes both viscous and elastic behaviour. In early process models such as the works done by Loos and Springer (1983) and Nelson and Cairns (1989) very simple elastic relationship were used to represent the behaviour of the composite material. The next generation of stress models made the assumption that at every time-step of the simulation the material behaviour may be considered elastic, but the moduli increase with the progress of cure. Some examples include the works of White and Hahn (1992), Bogetti and Gillespie (1992), Johnston et al (2001), Fernlund et al. (2002a, 2002b, 2003), and Antonucci et al. (2006). These models are referred to as Cure Hardening Instantaneously Linear Elastic or CHILE models, where the modulus of elasticity changes as a function of the instantaneous temperature and degree of cure. As a part of COMPRO, Johnston et al. (2001) presented a plane strain FE model for simulation of the development of process-induced deformation during autoclave processing of composites. A CHILE model to represent the mechanical behaviour of the 6  Chapter 1: Introduction  composite matrix resin, and micromechanics models were used to determine composite ply mechanical properties and behaviour, including thermal expansion and cureshrinkage. The effect of tooling on the final shape of the composite were also considered through simulation of tool-part interfaces and post-processing tool removal. Twigg et al. (2004a, 2004b) performed an experimental investigation into the effects of tool-part interaction on the final shape of cured composite parts, and presented analytical and numerical simulations to validate the experimentally obtained relationships. Osooly (2008) developed large deformation and contact capabilities applicable to the numerical simulation of stress development and deformation of curing composites, and numerically modeled the experiments done by Twigg et al. (2004a, 2004b) on tool-part interaction. Arafath et al. (2008, 2009) presented a cost-effective analytical solution for the process induced stresses and deformation in composite parts with flat and curved geometries cured on a solid tool. According to Zobeiry (2006) and Zobeiry et al. (2010), since the polymer spends considerable time in the viscoelastic regime during cure, it is important to check the validity and generality of a CHILE constitutive model relative to a full viscoelastic model. Using fundamental principles of viscoelasticity Zobeiry et al. (2010) showed that in a large majority of cases CHILE model is quite accurate and efficient in predicting the behaviour of the composite during processing leading to accurate predictions of residual stresses and/or dimensional distortions of the composite after tool removal. Using viscoelastic formulations to model the stress development during the cure of composites has gained more popularity during the past two decades. Some examples of these investigations include White and Kim (1998), Adolf and Martin (1996), Prasatya et al (2001), and Zhu et al. (2001), and Melo and Radford (2004). For these models, viscoelastic characterization of the behaviour of thermosetting polymers is crucial. Kim and White (1996) used dynamic mechanical analyzer (DMA) to characterize the relaxation behaviour of 3501-6 resin at different degrees of cure ranging from 0.57 to 0.98. Prasatya et al. (2001) presented a similar viscoelastic model for the bulk behaviour of Hexcel 8551-7 resin.  7  Chapter 1: Introduction  Zobeiry (2006) and Zobeiry et al. (2006) presented a differential implementation of the viscoelastic response of curing thermoset composites to predict the residual stresses and the final geometry of the composite part. A viscoelastic solid behaviour based on generalized Maxwell elements were assumed for the resin, and a solid micromechanics scheme was applied in the Laplace space. The results were then transferred back to time domain to obtain the coefficients of a generalized Maxwell model for the composite material. Zobeiry (2006) and Zobeiry et al. (2010) also presented a pseudo-viscoelastic formulation as an efficient approximation to fully viscoelastic models for the behaviour of thermoset resins during typical autoclave cycles. Zobeiry et al. (2010) showed that CHILE model is merely one form of the pseudo-viscoelastic approach where the modulus of resin is assumed to be instantaneously elastic, but the value of the modulus is set equal to the relaxation modulus of resin at a variable time, or its storage modulus at a certain frequency. In the current work, we will use the „variable time‟ version of the pseudoviscoelastic approach.  1.4. Motivation As it was mentioned before, the current approach to process modeling is to first run the flow module until the gelation of resin, then the geometry and volume fractions of the system are updated and the laminate is remeshed if needed. The final step is to run the stress model from the beginning of the processing cycle to capture the development of stresses in the composite system. There are a few drawbacks in the above approach; one is that there is an overlap in the simulation process as both flow model and stress model need to simulate the processing before the gelation of resin leading to an inherent inefficiency in the above method. Another disadvantage is that considering cases (e.g. sandwich panels) where at a certain time during the processing resin has gelled at one side of the composite structure but it is still a viscous fluid at the other side the identification of a gelation point for the system may prove to be difficult. More importantly, using the two-step approach does not capture the interaction of resin flow and its effects on how the stresses develop in the composite laminate.  8  Chapter 1: Introduction  If the flow model and stress model can be integrated into a unique model that simultaneously captures and simulates both the resin flow and stress development, the new integrated model would not have any of the disadvantages mentioned above for the separate application of the two models. This integrated approach would result in a more realistic representation of the complex phenomena which occur during the processing of composite. Not only does this method capture the development of the state variables relevant to flow and stress at the same time, but it will also predict the interactions of these phenomena during the processing.  1.5. Research Objectives and Thesis Outline The main objective of this thesis is to integrate the resin flow and stress development modules in COMPRO into a unique module that can capture both of these effects at the same time and in seamless fashion. In order to make this integration happen, we need to have a set of governing equations that encompasses both the two-phase formulation for flow through porous media (governing equations for flow) and the equilibrium equation (governing equation for stress development). To that end, we base our framework on a two-phase system for resin flow through the fibre-bed and identify and implement the necessary adjustments so that it also accurately predicts the stress development response of the composite material during cure. The most common approach to FE analysis of flow in porous media is the so-called „u-p‟ formulation, where the displacement of the solid skeleton and the pressure of the fluid phase are assumed to be main variables. In the u-p formulation pressure values are required to be pre-defined as boundary conditions at any permeable boundary. This may cause problems in our efforts to utilize the two-phase formulation for stress development in cured composites, as the fluid pressure values in that context are in fact a part of the total stress of the system. Also as will be discussed later, the finite elements based on u-p formulation are susceptible to pressure oscillations at the start of analysis which is likely to cause convergence problems whenever a nonlinear solution scheme is needed. These issues has led us to investigate the derivation of the governing equations for the twophase media and an alternative FE approach in u-v-p formulation. The additional 9  Chapter 1: Introduction  variable, v, represents the velocity of the fluid phase. The u-v-p formulation is superior to u-p formulation in the above-mentioned issues. It also accommodates the inclusion of shear stress components for the fluid phase that enables the approach to also model purely viscous fluids and their boundary layer with a two-phase system. In this work, we present the framework of integration in modeling of flow and stress for the pseudo-viscoelastic formulation presented by Zobeiry (2010). However, the framework can readily be adapted to the more comprehensive differential viscoelastic formulation introduced by Zobeiry (2006) and Zobeiry et al. (2006). Based on the above discussion, the thesis is organized as follows: Chapter 2 - Governing Equations and FE Model of Two-phase Media Derivation of the equations governing the behaviour of two-phase media involving distinct phases of fluid and solid constituents is reviewed. Specialized FE elements are introduced and formulated to model the two-phase behaviour based on a u-v-p formulation. Chapter 3 - Validation and Verification of Two-phase model Through examples varying from general porous flow problems to problems pertaining to flow of resin in the processing of composite laminates , the response of the developed FE model is validated. Also, the capability of the proposed formulation at capturing the flow behaviour at a boundary between a viscous fluid and a porous wall is demonstrated. Chapter 4 - Integration of Modeling from Fluid to Cured Resin The necessary modifications in the two-phase model in order to capture stress development as well as resin flow are identified and developed. These include but are not limited to modifications in the mass conservation equation to account for a compressibility consistent with the stress models for cured composites, introduction of a modified concept of effective stress, and modifications to the numerical solution technique.  10  Chapter 1: Introduction  Chapter 5 - Numerical Applications In the first section of this chapter, the application of the current two-phase model to stress behaviour in cured elastic composite materials is verified. The second part involves various examples from single-element problems dealing with the development of residual stresses in curing composite materials throughout a single-hold cure cycle to more geometrically complex composite laminates undergoing standard cure cycles. The response of the current integrated approach is compared with those obtained by common pseudo-viscoelastic approach. Chapter 6 - Conclusions and Further Work A discussion of the significance and contributions of the present work is provided. Recommendations for improvements and further development of the current framework are also presented.  11  Chapter 1: Introduction  1.6. Figures  Vacuum bag  Breather  Dam  Bleeder  Vacuum plug  Laminate  Tool (a)  (T) ()  Controller output  (Pb)  Heater (T)  Time  Cure cycle  () Laminate-tool in vacuum bag (Pb) Controller  Sensors output  Pressure source  Vac. pump Autoclave  (b) Figure 1.1: Autoclave processing steps, (a) lay-up, (b) cure with a controller loop (adapted from Hubert, 1996)  12  Chapter 1: Introduction  Input module  Autoclave controller module  Output module  State Variables Database Stress module  Thermo-chemical module  Flow module  Figure 1.2: Integrated sub-model approach in process modeling (adapted from Johnston et al., 2001)  180  1  160  0.9 0.8 0.7  120  0.6 100  Stress development regime  80  0.5 0.4  60  0.3  Flow regime  40  Autoclave Temp. Part Temp.  20  Degree of cure  0 0  50  100  Degree of cure  Temperature(°C)  140  150  200  250  300  0.2 0.1 0 350  Time (min)  Figure 1.3: Different regimes of resin behaviour during a typical autoclave processing history 13  Chapter 1: Introduction  180  1000  160 100  120 100 10 80 60  Autoclave Temp. 40  1  Part Temp.  Resin viscosity (Pa.s)  Temperature (°C)  140  Resin viscosity  20 0  0.1 0  50  100  150  200  250  300  350  time (min) Figure 1.4: Time-history of changes in resin viscosity during a typical autoclave processing history  Compaction  2, 3 1  Fibre-bed Figure 1.5: Schematic representation of resin flow through the fibre-bed  14  Chapter 2: Governing Equations and FE Model of Two-phase Media  2.1. Introduction Deformation and flow analysis of porous media using the finite element method has been carried out extensively during the past four decades. The pioneering works on this subject were due to Sandhu and Wilson (1969) on a linear displacement-pressure (u-p) formulation for seepage analysis, and Ghaboussi and Wilson (1972) on a variational formulation for the dynamic response of elastic porous solids. Both approaches dealt with elastic solid structures with the latter formulation leading to a displacement-relative velocity (u-v) FE treatment. Since then, the common agreement in the work of various authors has been to use the u-p representation in quasi-static problems such as consolidation (of porous media, specifically soils), and the u-v formulation in analysing the dynamic response of porous media such as liquefaction of soils under earthquake loading. This viewpoint, even though representing the common trend among the relevant articles published so far, does not seem to have entered the notable textbooks on the finite element method, as the u-p formulation seems to be the only FE representation of porous media discussed and documented in the FE textbooks for either quasi-static or dynamic analyses. Prevost (1983) presented displacement-pressure FE formulation for the linear and nonlinear analysis of consolidation. He also developed displacement-velocity formulations for the transient response of saturated porous media [Prevost(1982,1985)]. To the best of the author‟s knowledge, he was the first to use a penalty parameter to eliminate the pressure for the case of incompressible fluids. He also suggested that the fluid contribution to the equilibrium equations be always treated implicitly in order to remove the stringent time-step size restriction associated with the presence of a stiff fluid [Prevost(1985)].  15  Chapter 2: Governing Equations and FE Model of Two-phase Media  Simon et al. (1986a) evaluated and compared u-p and u-v formulations in modeling the dynamic response of saturated porous media in 1D problems. In a dynamic analysis, they concluded that the u-v method is slightly more accurate than the u-p formulation. In the derivation of the u-p formulation, one needs to rewrite the momentum balance equation of the fluid phase to express the relative velocity (or relative displacement) as a function of the other two major variables (u & p) to be substituted into the mass conservation equation. This procedure leads to the elimination of the relative velocity from the governing equations of the system. In order to rewrite the fluid phase equilibrium equation in the above fashion, it is commonly assumed that in the generalized Darcy‟s law the inertial terms containing the time derivative of fluid relative velocity are negligible. This should play a major role in the u-p formulation‟s slightly inferior results compared to those of u-v in the dynamic behaviour of flow in porous media. Simon et al. (1986b) reviewed and compared the u-v, u-v-p, u-v-p-σ formulations among others for the dynamic analysis of saturated porous media. They concluded that compared to u-v formulation, u-v-p and u-v-p-σ approaches result in slightly more accurate results for fluid pressure and total stress with a minimal loss in efficiency. In their representation of the sub-matrices for the u-v-p approach, however, it seems that there is a typographical error: the vp component of the stiffness matrix includes the divergence of the pressure shape functions which is not consistent with their relevant force term as it includes the applied pressure on the boundary of the system. As we will see in the formulation presented here, to avoid C0 continuity for pressure degrees of freedom between elements, we use integration by parts to eliminate any spatial differentiation of the pressure shape functions, and it is only then that the applied pressure on the boundary of the system appears in the relevant force term. To derive the u-v formulation, one needs to substitute the mass conservation equation into the equilibrium equations of the system. Pressure is solved with respect to the other variables and is then substituted into the generalized Darcy‟s law and the total equilibrium equation. If the system is incompressible, a u-v formulation cannot be achieved since there is no pressure term in the mass conservation equation. In order to derive u-v formulations for incompressible cases, some researchers have introduced a large penalty parameter into the mass conservation of the system as the bulk modulus of the system, and thus incorporating pressure in the continuity equation. 16  Chapter 2: Governing Equations and FE Model of Two-phase Media  Depending on the elimination of pressure at the differential equations level or after FE implementation, these approaches have been called penalty and mixed-penalty methods respectively by Spilker and Maxian (1990). They developed a so-called mixed-penalty FE formulation for linear biphasic response of soft tissues and compared it with a penalty form previously proposed. They found the mixed-penalty method to be superior over the penalty method in locking issues and sensitivity to mesh irregularities. However, they also noted that the two methods produce the same equations for four-node quadrilaterals. In a two-part paper, Almeida and Spilker (1997,1998) used the so-called mixed penalty method and also velocity-pressure (similar and equivalent to u-p) FE formulations to solve non-linear biphasic equations including strain-dependent permeability and a hyperelastic behaviour for the solid phase. They found their velocity-pressure method to be more efficient computationally due to the lower number of DOFs involved. However, they observed that the velocity-pressure formulation produces fluid and relative velocity fields that are not as accurate as the ones from the mixed-penalty method. They also reported some difficulties with the velocity-pressure formulation associated with oscillations of the stress and strain fields at early times and near the loaded surface that prevented the convergence of the nonlinear solution. However, considering the computational time, they gave the advantage to the velocity-pressure formulation at least for quasi-static problems. They considered dynamic problems the turf in which mixedpenalty formulation would fare stronger than the velocity-pressure formulation. In dynamic problems such as wave propagation, the extension of the velocity-pressure formulation would mean that either the fluid inertia terms be neglected or the fluid momentum equation be integrated in time. They argued that the former option is acceptable in soil mechanics with porosities around 0.3, but would not be very suitable for soft tissues with porosities in the range of 0.75 or more. Chan et al. (2000) presented a mixed-penalty finite element model for the behaviour of articular cartilage in the biomechanics of diarthrodial joints. Using mixture theory, they added a viscous shear term to the hydrostatic pressure to have a complete stress tensor for the fluid phase. They argue that this change makes it possible to model a single-phase continuum as a limiting case of biphasic material if the solid or fluid volume fraction is 17  Chapter 2: Governing Equations and FE Model of Two-phase Media  set to zero. This characteristic is helpful in modeling the interfacial response and interaction of porous flow and viscous flow. By the help of special assembly considerations on the interface, and using their FE model, they simulated Taylor‟s problem of Couette flow over a rigid (or deformable) porous layer. Here, volume averaging technique is used to derive the complete set of governing equations for the response of a two-phase system including a solid structure and fluid matrix. The shear stress terms of the fluid phase are not neglected leading to the DarcyBrinkman equation for the momentum balance of the fluid phase. Neither fluid pressure nor relative velocity of fluid are eliminated from the formulation to have a complete u-vp representation of the behaviour.  2.2. Volume Averaging Method In the general case of a two-phase medium where a fluid interacts with a porous structure, one needs to know how the fluid is transported through the pores from one point to another in the domain of the problem. The direct analysis of this problem in terms of transport equations that are valid within the pores is practically impossible due to the complex structure of the typical porous medium. Rather than tackle this problem in the form of equations that are valid in the pores, the equations relevant to each phase can be spatially averaged to produce equations that are valid within the domain of the two-phase medium. The method of volume averaging is a technique developed to rigorously generalize continuum equations to multiphase media [Whitaker (1998)]. Tucker and Dessenberger (1994) applied this method to the composite processing method of resin transfer moulding (RTM). They developed the governing equations for flow and heat transfer through stationary fibre-beds. Considering the solid phase to be stationary, is a very relevant assumption in RTM. By defining a phase function, it is possible to have a formal description of the microscopic geometry of the porous medium. The phase function for the  phase in a multi-phase medium is if x lies in the  - phase 1 X  ( x)   0 if x doesn' t lie in the  - phase  (2.1)  18  Chapter 2: Governing Equations and FE Model of Two-phase Media  Three distinct average values for a variable in a multi-phase medium may be defined. The spatial average is the average value of a parameter in all phases within a representative volume V   B   1 BdV V V  (2.2)  Phase average is defined as the average value of a parameter at points that lie within a single phase α, but still over the entire volume V  B    1 1 B dV   B X  dV  V V VV  (2.3)  Intrinsic phase average is defined as the average value of a parameter at points that lie within a single phase α over the volume occupied by that phase. It can be represented as  1  B   V    B dV    B X  dV V  V   X  dV  (2.4)  V  Volume fraction of the α-phase is defined as     V 1  X  dV V V V  (2.5)  The phase and intrinsic phase averages for an arbitrary α-phase are related by  B     B    (2.6)  In the following, we review some useful theorems of volume averaging that are very useful in deriving the averaged governing equations of a two-phase system. If B is continuous within the α-phase, the first theorem of volume averaging states that the average of gradient of B in the α-phase may be expressed as [Gray and Lee (1977), Whitaker (1998)] B    B    1 V   B n dS  (2.7)  S  19  Chapter 2: Governing Equations and FE Model of Two-phase Media  where  denotes the other phase of the two-phase medium. S is the interfacial surface between α-phase and β-phase, and n  is the unit vector normal to S directed from αto β-phase. A similar equation can be written for the average of divergence of a tensor as follows [Whitaker (1998)]   A       A     1 V   A  n dS.  (2.8)  S  2.3. Governing Equations Assuming a two-phase system consisting of a fluid phase and a solid porous structure where each phase is incompressible, the microscopic equations for mass conservation of the fluid phase and the solid phase may be written as   v m  0,  (2.9)  v f  0  where v is the vector of velocity, and f and m represent the solid phase (fibres) and the fluid phase (resin matrix) respectively. Taking the phase average of the microscopic equation of mass conservation of the fluid phase in (2.9) and incorporating the theorem presented in (2.8) leads to   vm    1 V  v  m   n mf dS  0  (2.10)  S mf  where S mf is the interfacial surface between the fluid and solid phase, and n mf is the unit vector normal to S mf directed from the fluid phase to the solid phase. Using the transport theorem, Tucker and Dessenberger (1994) showed that the second term on the left-handside of the above equation is equal to the rate of resin volume fraction, i.e. 1 V  v S mf  m   n mf dS    m t  (2.11)  20  Chapter 2: Governing Equations and FE Model of Two-phase Media  Therefore, the macroscopic form of the mass conservation equation of the fluid phase may be written as   m    vm   0 t  (2.12)  The same procedure may be performed to obtain the macroscopic form of the mass conservation equation of the solid phase in the form of   f     v f   0  t  (2.13)  Assuming that the pores are fully saturated by the fluid phase, the volume fractions of the fibre-bed and the matrix are related by   f  m  1  (2.14)  Adding the two mass conservation equations in (2.12) and (2.13) together and combining the result with (2.14) leads to the mass conservation equation of the two-phase system   v f     vm   0  (2.15)  Assuming that the inertia is not of consequence and the effect of body forces including gravity is negligible, the microscopic equations for equilibrium equation of the fluid phase and the solid phase take the simple form of   ζ m  0, ζ f  0  (2.16)  where ζ represent the stress tensor. Taking the phase average of the microscopic equilibrium equation of the fluid phase in (2.16) and incorporating the theorem presented in (2.8) leads to   ζ m    1 V  ζ  m   n mf dS  0  (2.17)  S mf  The total stress tensor of the fluid phase may be decomposed into a hydrostatic pressure term Pm and a deviatoric stress tensor ηm 21  Chapter 2: Governing Equations and FE Model of Two-phase Media  ζ m  η m  Pm I  (2.18)  where I is the second-order unit tensor. Thus, the first term on the left-hand-side of (2.17) may be written in the form of   ζ m      η m    Pm      η m   m Pm  m   Pm  m m  (2.19)  Substituting (2.19) into (2.17) leads to 1 V    m  Pm  m     η m    Pm  m  m   ζ  m   n mf dS  0  (2.20)  S mf  If the vector of drag forces between the fluid phase and the solid phase is defined as f d   Pm  m  m   1 V  ζ  m   n mf dS  (2.21)  S mf  the equilibrium equation of the fluid phase takes the final form of  m  Pm  m     η m   f d  0  (2.22)  Taking the phase average of the microscopic equilibrium equation of the solid phase in (2.16) and incorporating the theorem presented in (2.8) leads to   ζ f    1 V  ζ  f   n fm dS  0  (2.23)  S mf  If we assume that the surface tension between the fluid and solid particles may be neglected, then we may write [Tucker and Dessenberger (1994)] ζ f  n fm  ζ m  n mf  0  on Smf  (2.24)  which leads to 1 V  ζ S mf  f   n fm dS    1 V  ζ  m   n mf dS  (2.25)  Smf  Considering the definition of drag force in (2.21, the above equation may be rewritten as  22  Chapter 2: Governing Equations and FE Model of Two-phase Media  1 V  ζ  f   n fm dS  f d   Pm  m  m  (2.26)  S mf  Substituting (2.26) into (2.23) leads to the following form of the equilibrium equation of the solid phase   ζ f   f d   Pm  m m  0  (2.27)  Assuming that the absolute value of pressure is directly transferred to the solid particles one may decompose the intrinsic phase average of the fibre-bed stress tensor into hydrostatic pressure of the fluid (which is transferred to the fibres) and η f  f ζ f  f   η f  f   Pm  m I  (2.28)  where η f  f represents a stress component dependent on the deformation of the fibrebed and may be related to the effective stress of the fibre-bed,  η f  , by (2.6). Applying (2.6) to the left-hand-side of the above equation leads to ζ f    f  η f  f   f  Pm  m I  (2.29)  Substituting (2.14) in (2.29), we arrive at ζ f    η f    Pm  m 1  m I  (2.30)  which leads to    ζ f      η f    Pm  m 1   m      η f    Pm  m 1   m    Pm  m  m  (2.31)  Substituting (2.31) into (2.27) leads to    η f   f d  Pm  m 1  m   0  (2.32)  The mass conservation equation in (2.15) combined with the equilibrium equations of the two phases in (2.22) and (2.32) leads to the system of governing equations of the twophase system  23  Chapter 2: Governing Equations and FE Model of Two-phase Media     v f     v m   0  m    m  Pm      η m   f d  0    η   f   P  m 1     0 f d m m   (2.33)  The constitutive behaviour of the fibre-bed could be expressed as      1  η f   C u f  f  u f  f 2   T  (2.34)  where u is the displacement vector, and C is the material stiffness tensor of the fibre-bed. For the shear stress in a viscous fluid, we may write    η m   v m  (v m )T    (2.35)  where μ is the viscosity of the fluid phase. The volume-averaged form of the above equation is     η m    v m   v m  T    (2.36)  in which, for v m  , we may rewrite the volume averaging theorem in (2.7) to have v m    v m    1 V  v  m  n mf dS  (2.37)  S mf  It is commonly assumed that the slip velocity at the interface is zero [Tucker and Dessenberger (1994)], and therefore vm  v f  on Smf  (2.38)  Thus, the integral term on the right-hand-side of (2.37), may be rewritten as v m    v m    1 V  v  f  n mf dS  (2.39)  Smf  Arguing that the changes in the solid phase velocity v f are small in comparison to the changes in fluid phase velocity v m , we assume that v f on the interface can be  24  Chapter 2: Governing Equations and FE Model of Two-phase Media  approximated with v f  f . Incorporating this assumption into the integral term in (2.39), we have    1 1 f f 1  v n dS   v  n dS   v  n dS m mf f mf f mf  V S  V Smf V Smf  mf   (2.40)  Applying the averaging theorem in (2.7) to the phase function Xm leads to 1 V  X m    X m    X  m  (2.41)  n mf dS  S mf  Within the fluid phase, Xm is equal to unity and X m equals zero. Also, it is readily deduced from (2.3) that  X m  is simply the fluid volume fraction m . Therefore, (2.41) may be simplified to [Tucker and Dessenberger (1994)]  m    1 V  n  mf  (2.42)  dS  S mf  Substituting (2.42) into (2.40) results in 1 V  v  m  n mf dS   v f  f m  (2.43)  S mf  which in turn, if substituted in (2.39) leads to v m    v m    v f  f m  (2.44)  Substituting the above into (2.36) and neglecting the approximations, we get        (2.45)       (2.46)   η m     v m    v m  T   v f  f m  m  v f  f  Rewriting the above equation leads to     η m     v m    v f  f m   v m    v f  f m  T  The drag force between the matrix and fibre-bed may be written in the following form that is consistent with Darcy‟s law [Whitaker (1998), Tucker and Dessenberger (1994)] 25  Chapter 2: Governing Equations and FE Model of Two-phase Media    f d  m2 S 1  v m  m   v f  f    (2.47)  where S is the matrix of permeability of the fibre-bed. Substituting (2.34) and (2.47) into (2.32), we arrive at the final form of the solid phase equilibrium equation      1   C u f  f  u f  f 2     T  2 m  S 1  v m  m   v f  f    Pm  m 1  m   0  (2.48)  Substituting (2.46) and (2.47) into the second equation of (2.33) leads to the final form of the fluid matrix equilibrium equation        v    0    m  Pm  m      v m    v f  f  m   v m    v f  f  m      m2 S 1  v m  m    T  f  (2.49)  f  Thus, the set of equations (2.33) may be expressed in the form     v f     v m   0  T m f f    m  Pm       v m    v f   m   v m    v f   m  -  m2 S 1  v m  m   v f  f  0  1 f f T   m2 S 1  v m  m   v f  f   Pm  m 1   m   0    C u f   u f  2                          (2.50)  In order to simplify the above equation, we define the velocity vector of the composite system (which is essentially equal to the velocity vector of the fibre-bed) as vc   v f  f  (2.51)  and the volume-averaged relative velocity (or flow velocity) of the resin as    v flow  m  v m  m   v f  f    (2.52)  The above definition of flow velocity is commonly used in the literature in FE formulation of flow in porous media [e.g. Lewis and Schrefler (1998)]. Based on (2.51) and (2.52) we may write the following relationships between the different velocity fields  26  Chapter 2: Governing Equations and FE Model of Two-phase Media   v f  f  vc  v f    f vc   v m  m  v flow  m v c  m  (2.53)   v m   v flow  m v c Using the above equations, one may easily arrive at the governing equations with respect to the new velocity variables     v c    v flow  0   T m 1   m  Pm      v flow   m v c   v flow   m v c  -  m S v flow  0  1 T   C u c  u c    m S 1 v flow  1   m  Pm  m  0  2           (2.54)  where v c  u c  (2.55)  and uc is the composite displacement field. For the sake of convenience, variables uc, vflow,  Pm  m , and m are replaced with u, v, P, and  respectively, i.e.  uc  u v flow  v  (2.56)   Pm  m  P  m   Using index notations, the set of governing equations (2.54) may be written in the form u i ,i  vi ,i  0    (2.57)     P,i  vi , j  ui , j   v j ,i  u j ,i  , j   S ij1v j  0      1 Cijkl u k ,l  ul ,k  , j   S ij1v j  1   P,i  0 2  (2.58)  (2.59)  Adding up the two equilibrium equations of (2.58) and (2.59), leads to the total equilibrium equation of the domain 27  Chapter 2: Governing Equations and FE Model of Two-phase Media         1 Cijkl u k ,l  ul ,k  , j   vi , j  u i , j   v j ,i  u j ,i  , j  P,i  0 2  (2.60)  The mass conservation equation in (2.57), along with the equilibrium equation of the matrix in (2.58), and the total equilibrium equation (2.60) constitute the set of governing equations    u i ,i  vi ,i  0   P,i   vi , j  u i , j   v j ,i  u j ,i  , j   S ij1v j  0  1  Cijkl u k ,l  u l ,k  , j   vi , j  u i , j   v j ,i  u j ,i  , j  P,i  0 2           (2.61)    that will be used in formulating the finite element discretization equations. The complete form of the governing equations presented above is especially useful in modeling the interaction at the boundary of a two-phase system with a neighbouring viscous fluid. As the shear stress of the fluid is not neglected, a few adjustment in the element leads to a correct representation of a purely viscous fluid. This will be discussed more in a later section. Also, the complete presentation of the stress components of the fluid phase helps in establishing the theoretical framework to extend the two-phase formulation to the cured composite material in which the resin phase is capable of carrying shear stress. In many cases of real-life porous media, the terms involving the shear stresses of the fluid phase are usually negligible compared to the terms involving the permeability of the system [Tucker and Dessenberger (1994)]. Therefore, the above set of governing equations in (2.61) may be commonly found in the literature in the following simplified form where the shear stress contributions of the fluid phase are neglected.   u i ,i  vi ,i  0    P,i   S ij1v j  0  1  Cijkl u k ,l  u l ,k  , j  P,i  0 2    (2.62)    28  Chapter 2: Governing Equations and FE Model of Two-phase Media  2.4. Finite Element Implementation In this section, the general set of governing equations in (2.61) is implemented in the finite element method. Two different elements are introduced and formulated based on a u-v-p formulation where the main degrees of freedom are displacement of the system, relative fluid velocity, and the pressure of the fluid phase.  2.4.1. Q2P-1 element A 2D Q2P-1 finite element is developed considering the set of governing equations in (2.61) compared to the regular finite elements based on the u-p formulation utilized for modeling flow in porous media in the commercial FE codes (such as ABAQUS). In the FE modeling of incompressible Navier-Stokes equations, the terminology of QmP-n is related to 2D quadrilateral and 3D hexahedral elements. It indicates that the velocity is approximated by a continuous piecewise polynomial of degree m in each direction and a pressure approximation that is a discontinuous piecewise polynomial of degree n (not of degree n in each direction, for 2D quadrilaterals as if the pressure is represented by a triangle within the quadrilateral). This element passes the LBB condition for stability and is considered to be one of the most accurate 2D elements in incompressible viscous flow [Gresho and Sani (2000)]. The Q2P-1 element used in FE representation of incompressible Navier-Stokes flow is the inspiration for the namesake element developed in this work where the kinematic degrees of freedom (system displacement and fluid velocity) are approximated by 2nd degree polynomials in each direction and the fluid pressure is approximated by linear distribution over the element. The Q2P-1 element presented here is a 2D bi-quadratic isoparametric element with 9 nodes for the system displacement and relative velocity of the fluid phase. As depicted in Figure 2.1, three internal nodes are assigned to pressure of the fluid phase, therefore enabling every element to represent its internal pressure distribution as a linear surface. For simplicity, we will also refer to this element as the 9-3 element.  29  Chapter 2: Governing Equations and FE Model of Two-phase Media  2.4.2. Galerkin finite element equations Integrating the equations in (2.61) over a domain Ω using some weight functions wp, and using integration by parts on occasion we arrive at   (u  i ,i   vi ,i ) wp d  0  (2.63)     Pw    p ,i  d   ,i Pw p d    (vi , j  v j ,i ) w p , j d        (u i , j  u j ,i ) w p , j d    S v j w p d   ( ij  P ij )n j w p d  0 1 ij      (2.64)     12  Cijkl (u k ,l  u l ,k ) w p , j d    (vi , j  v j ,i ) w p , j d        (u i , j  u j ,i ) w p , j d   Pw p ,i d   ( ij   ij  P ij )n j w p d  0     (2.65)    where Γ denotes the boundary of the domain and Γσ is the portion of the boundary where tractions are specified. The components of effective stress tensor in the fibre-bed are defined by  1 2   ij  Cijkl (u k ,l  ul ,k )  (2.66)  and the components of shear stress tensor in the resin are expressed as   ij  vi , j  ui , j   v j ,i  u j ,i   (2.67)  The term ( ij   ij  P ij )n j in (2.65) appears as the result of the application of integration by parts to all the three terms in the last equation of (2.61) and represents the vector of total traction force applied at the boundary of the system. The term ( ij  P ij )n j in (2.64) appears as the result of the application of integration by parts to  the first and second terms in the second equation in (2.61) and represents the fluid-phase share of the vector of total traction force applied at the boundary of the system.  30  Chapter 2: Governing Equations and FE Model of Two-phase Media  Writing Equations (2.65), (2.64), and (2.63) respectively in matrix form for the 9-3 finite element depicted in Figure 2.1, we arrive at  K uu  K vu K pu   K uv K vv K pv  K up  u  C uu    K vp  v   C vu K pp  p  C pu  C up   u   f u      C vp   v    f v  t C pp   p   0  391  C uv C vv C pv  (2.68)  where u181 , v 181 , and p31 are the vectors of the degrees of freedom relating to the total displacement field, relative resin velocity, and resin pressure respectively. Defining the traction vector on the boundary of the two-phase medium as t i( n )   ij   ij  P ij n j  (2.69)  we may rewrite (2.65) to have 1 2  C    ijkl  (u k ,l  u l ,k ) w p , j d    (vi , j  v j ,i ) w p , j d      (u i , j  u j ,i ) w p , j d   Pw p ,i d   t i( n ) w p d     (2.70)    which leads to the equation in the first row of (2.68). Hence, we may arrive at the definitions of some of the components of the matrices in (2.68) K uu   BT DT Bd   K uv   B T B1d   (2.71)  1818  (2.72) 1818  K up    B T δN p d   C uu    B T B1d   Cuv  01818  (2.73)  183  (2.74)  1818  (2.75)  31  Chapter 2: Governing Equations and FE Model of Two-phase Media  Cup  0183  (2.76)  f u   N Td t ( n ) d   (2.77) 181  Where t is the total traction vector on the boundary of the medium with components in xand y-directions in the form of  t  (n)  (n)  t x     (n)   t y    (2.78)  and Nd is defined as N 0  Nd     0 N  (2.79)  where N is the matrix of shape functions which interpolate the displacement and relative velocity fields on the element. For a bi-quadratic 9-node element, N is expressed in the form N  N1  N2  ... N9   (2.80)  and assuming that the element has three internal nodes for pressure interpolation, we may introduce the matrix of pressure shape functions as    N p  N1p  N 2p  N3p    (2.81)  The displacement, resin relative velocity, and the pressure fields are interpolated on the domain of the element using the following equations  32  Chapter 2: Governing Equations and FE Model of Two-phase Media 9  u x  Nu x   N iu x i i 1 9  u y  Nu y   N iu y i i 1 9  v x  Nv x   N i v x i  (2.82)  i 1 9  v y  Nv y   N i v y i i 1  9  p  N p   N i pi i 1  The displacement shape functions may be written as follows 1 N 1   1   1    4 1 N 2    1   1    4 1 N 3   1   1    4 1 N 4    1   1    4 1 N 5    1   2 1    2 1 N 6   1    1   2 2 1 N 7   1   2 1    2 1 N 8    1    1   2 2 N9  1   2 1  2                (2.83)        The definition presented in (2.79) is relevant if the vector of main variables is defined as  ux     uy  v   x vy  p   391  (2.84)  33  Chapter 2: Governing Equations and FE Model of Two-phase Media  while we are interested to present the vectors of displacement and relative resin velocity degrees of freedom in the form   u1x     u1 y  u   2x  u   u2 y ,       u9 x     u9 y    v1x     v1 y  v   2x  v   v2 y        v9 x     v9 y   (2.85)  which leads to rewriting (2.79) so that N Nd   1 0  0 N 9    N9  0  N2  0  N1  0  N2   0  (2.86)  B is a matrix containing the spatial derivatives of the shape functions defined in (2.83), and may be expressed in the form   N1, x  B 0  N1, y    N 9, x  0  N 2, x  0  N1, y  0  N 2, y    N1, x  N 2, y  N 2, x   N 9, y  0  0   N 9, y  N 9, x  318  (2.87)  B1 is another matrix containing the spatial derivatives of the shape functions defined in (2.83), and is defined as  2 N1, x  B1   0  N1, y   0  2 N 2, x  2 N1, y  0  N1, x  N 2, y   2 N 9, x  0  2 N 2, y    N 2, x  0 N 9, y  0   2 N 9, y  N 9, x  318  (2.88)  δ is defined by δ  1 1 0  T  (2.89)  34  Chapter 2: Governing Equations and FE Model of Two-phase Media  To calculate the spatial derivatives of the shape functions, we need to establish a relationship between the general (x,y) coordinates and the local (ξ,η) coordinates through the Jacobian matrix J as follows  x   J  x    y     y     (2.90)  As the element is isoparametric, the geometry is interpolated over the element by the same shape functions as displacement and velocity fields, thus leading to 9  9  i 1  i 1  x   N i xi , y   N i yi  (2.91)  which upon substitution in (2.90) yields J J   11  J 21  J12  J 22   (2.92)  where 9  J11   N i , xi , i 1 9  J 21   N i , xi , i 1  9  J12   N i , yi , i 1 9  J12   N i , yi .  (2.93)  i 1  Using the inverse of Jacobian matrix, we may define a matrix consisting of the spatial derivative of the shape functions as N, x  N ,  B     J 1   N ,  N , y   (2.94)  where the inverse of Jacobian can be easily written as  35  Chapter 2: Governing Equations and FE Model of Two-phase Media  J 1    J 12  J 11    J 22 1 J 11 J 22  J 12 J 21  J 21  (2.95)  B facilitates the definition of B based on the components of B as   B11 0 B12  B   0 B21 0  B21 B11 B22   0   B19  B22   0  B12  B29  0   B29  B19  318  (2.96)  To determine the pressure shape functions Np, let‟s assume that the coordinates of the three pressure nodes are P1 1,1 , P2 2 ,2 , P3 3 ,3   (2.97)  To find the shape function for the first pressure node, we write the equation of the straight line that passes through nodes 2 and 3 to have  3  2     2       2   2  3   0  3   2   3   2   (2.98)  We may set the first shape function for pressure to be     2      2       2   2  3  N1p  a  3       3 2 3 2      (2.99)  and the parameter a may be calculated by setting N1 1 ,1   1 . Substituting the calculated value of a into (2.99) leads to  N1p    3  2       3   2     2   2  3 2    21  12   32   23   13  31   3   2    3   2   (2.100) The other two pressure shape functions of the 9-3 element can be derived in a similar approach  36  Chapter 2: Governing Equations and FE Model of Two-phase Media  N 2p    1  3       1  3     3  3  1 3  (2.101)   21  12   32   23   13  31   1  3   1  3   N 3p    2  1        2  1     1  1  2 1  (2.102)   21  12   32   23   13  31    2  1    2  1   We choose the location of the pressure nodes to be    P1  1  3 ,1     3 , P2 1  3 , 1     3 , P3 0,1  3    which leads to the simplification of the equations for the pressure shape functions as follows N1p   N 2p   3 1    2     4  3  (2.103)  3 1   2     4  3  (2.104)  3 1     2  3  (2.105)  N 3p   In equation (2.71) DT is the matrix of tangent moduli of the fibre-bed, which in the case of a 2D isotropic plane strain problem may be described as   1  E  DT   1  1  1  2    0 0    0  1  2  2 0  (2.106)  with E being the Young‟s modulus and  the Poisson‟s ratio of the fibre-bed. Defining the resin‟s share of the applied traction vector on the boundary of the system as t m( n) i   ij  P ij n j  (2.107)  37  Chapter 2: Governing Equations and FE Model of Two-phase Media  we may rewrite (2.64) to have    (v    i, j   v j ,i ) w p , j d    S ij1v j w p d   Pw p ,i d     (2.108)    ,i Pw p d    (u i , j  u j ,i ) w p , j d   t m( n ) i w p d       Equation (2.108) then leads to the system of equations related to the second row of matrices in Eqn. (2.68). Thus, we may arrive at the definition of the components on the second row of the matrices in (2.68) K vu  01818  (2.109)      K vv   B T B1d    N Td S 1N d d    B T B1  N Td S 1N d d       (2.110)  1818      K vp    BT δN p d   12 BT2 δN p d    BT  12 BT2 δN p d  (2.111)  C vu    B T B1d  (2.112)          183  1818  Cvv  01818  (2.113)  Cvp  0183  (2.114)  f v   N Td t (mn ) d   (2.115) 181  Where B2 is a matrix consisting of the displacement shape functions and the spatial derivatives and the resin volume fraction  in the form of  2, x N1 0 2, x N 2  B2   0 2, y N1 0  , y N1 , x N1 , y N 2   0   2, x N 9  2, y N 2   , x N 2    0  , y N 9    2, y N 9  , x N 9  318 0  (2.116)  38  Chapter 2: Governing Equations and FE Model of Two-phase Media  With the assumption of incompressibility of the phases, we may write the following equations for the initial and current values of the resin volume fraction  0   Vm 0 Vc 0    ,  Vm 0  Vm Vc 0  Vm  (2.117)  where 0 is the initial value of resin volume fraction, Vm0 is the initial volume of the resin, and ΔVm is the change in the volume of resin due to flow which is the only source of any change in the volume of the system. Vc0 is the initial total volume of the composite system. Dividing the numerator and denominator of the second equation in (2.117) and taking into account that that the volumetric strain of the system is defined as  v   Vm Vc 0  (2.118)  leads to the following equation for the volume fraction of resin    0   v 1 v  (2.119)  Substituting the relationship between the volumetric strain and the total displacement field results in    0  u m , m 1 un ,n  (2.120)  To obtain the spatial differentiation of the fluid phase volume fraction on the domain of an element from the values of  at the Gauss points, we introduce a special set of shape functions with base points located at the 3×3 set of Gauss points (0,  15 5 ) as  39  Chapter 2: Governing Equations and FE Model of Two-phase Media                     25  15 5   15 5   36 25 N 2    15 5   15 5   36 25  N 3   15 5   15 5   36 25  N 4    15 5   15 5   36 25  N 5    3 5   2 15 5   18 25  N 6   15 5   3 5   2 18 25  N 7   3 5   2 15 5   18 25  N 8    15 5   3 5   2 18 25  N9  3 5   2 3 5  2 9 N 1                      (2.121)          Using the above set of shape functions we may write an equation for the field of fluid volume fraction in the form of 9    N θ   N ii   (2.122)  i 1  where θ is the vector of the values of fluid volume fraction at Gauss points. Spatial gradients of  may be obtained by  , x  N , x θ,  , y  N , y θ  (2.123)  The system of equations in the last row of the matrix equation (2.68) represents the mass conservation equation (2.63) that may be rewritten as:  v  i ,i    wp d   ui ,i wp d  0  (2.124)    which then leads to defining the components on the last row of the matrices in (2.68) K pu  Cpv  0 318  (2.125)  40  Chapter 2: Governing Equations and FE Model of Two-phase Media  K pv  Cpu   N p δ T Bd T    (2.126) 318  K pp  Cpp  0 33  (2.127)  Hence, the matrix equation (2.68) may be written in the following simplified form  K uu   0  0   K uv K vv K pv  K up  u  C uu    K vp  v   C vu 0  p  C pu  0 0  u   f u      0 0  v    f v  t 0 0  p   0  391  (2.128)  2.4.3. Q1P0 element A Q1P0 element is also developed based on the same set of governing equations as the previous element. As the naming system suggests, the kinematic degrees of freedom are approximated by a linear polynomial in each direction and the pressure is approximated by a constant value over the surface of the element. This element is a bi-linear isoparametric element with 4 nodes for the system displacement and relative velocity of the fluid phase, and only one central node assigned to pressure of the fluid phase. This element is depicted schematically in Figure 2.2. This element will also be referred to as the 4-1 element. The matrix of shape functions, N, is expressed in the form N  N1  N2  N3  N4   (2.129)  where the kinematic shape functions are defined as  1 1   1    4 1 N 2  1   1    4 1 N 3  1   1    4 1 N 4  1   1    4 N1   (2.130)  41  Chapter 2: Governing Equations and FE Model of Two-phase Media  For this element, there is only one shape function for pressure values which is defined as unity since the pressure node is located at the origin of ξ and η axes. N p 1  (2.131)  The vector of main variables for this element consists of 8 displacement variables, 8 relative velocity values, and only one variable representing pressure of the fluid phase.  ux     uy  v   x vy  p  171  (2.132)  Similar to the case of 9-3 element, we introduce a special set of shape functions with base points located at the 2×2 set of Gauss points (  3 3 ) as  3 4 3  N2  4 3  N3  4 3  N4  4 N 1         3 3    3 3    3 3    3 3   3 3   3 3  3 3  3 3   (2.133)  The rest of the matrices and vectors relevant to the 4-1 element, are calculated in a similar fashion to those of the 9-3 element. One would find the only difference in the size of these matrices, and therefore we avoid repeating the formulas here.  2.4.4. Boundary conditions Due to the inclusion of shear stress components in the stress tensor of the fluid phase in the current formulation, we need to specify a complete set of boundary conditions for the fluid phase as well as the whole medium. The BC‟s relevant to the fluid phase are essentially similar to the BC‟s used in the FE analysis of 2D incompressible viscous flow,  42  Chapter 2: Governing Equations and FE Model of Two-phase Media  and the BC‟s relevant to the whole medium are the same as the ones introduced in FE treatment of 2D elastic solids. Effectively, the two-phase problem is treated as a superposition of two problems; the whole medium with displacement/traction force boundary conditions, and the fluid phase with relative velocity/traction force boundary conditions. Distinct combinations of BC may be assumed, including:   Free and impermeable BC, where the total traction vector on the system and the relative velocity of the fluid are pre-determined: t i( n )  f appl. vi  vi def .    Free and permeable BC, where the total traction vector on the system and also the traction vector on the fluid phase are set to the applied external forces on the system and the matrix phase respectively: t i( n )  f appl. t m( n ) i  f appl. m  If the only applied force on the boundary is due to air pressure, then the share of the matrix phase of the applied force on the boundary is proportional to the area occupied by the matrix phase on the relevant boundary surface. In the case of isotropic materials, this ratio is equal to the volume fraction of the matrix phase, and therefore we have f appl. m    f appl. .   Constrained and impermeable BC, where both the displacement vector of the system and the relative velocity of the matrix phase are set to pre-defined values: ui  ui def . vi  vi def .  43  Chapter 2: Governing Equations and FE Model of Two-phase Media    Constrained and permeable BC, where the displacement vector of the system is set to pre-defined values. Also, the traction vector applied to the matrix phase is set to the share of that phase of the total traction applied to the boundary: ui  ui def . t m( n ) i  f appl. m  Since in this case, the only plausible applied traction is due to ambient air pressure on the boundary, we may write  f appl. m    f appl. for isotropic two-phase  materials.  2.4.5. Time integration To discuss the numerical solution of (2.128), we rewrite the equation in a short form to have  F KX  CX  (2.134)  u   X  v p    (2.135)  where  is the generalized vector of main variables. Using the generalized mid-point integration rule, we may write the values of X and its time derivative at an arbitrary time between tn and tn+1 in the form  X n   X n 1  X n  t  (2.136)  Xn  1   Xn  Xn1  (2.137)  44  Chapter 2: Governing Equations and FE Model of Two-phase Media  where 0    1 . A few well-known special cases of the above scheme are Forward Euler (   0 ), Backward Euler (   1 ), and Crank-Nicholson scheme (   0.5 ). Writing equation (2.134) at time tn+θ leads to  K n Xn  Cn X n   Fn   (2.138)  Substituting (2.136) and (2.137) into the above equation, we have K n 1   Xn  Xn1   Cn Xn1  Xn  t  Fn  (2.139)  which after some manipulations, leads to a system of equations for the vector of main variables at the current time tn+1  C  tKn Xn1  C  1   tKn Xn  tFn  (2.140)  or in an expanded form as  Cuu  tK uu tK uv tK up    C vu tK vv tK vp    Cpu tK pv 0   Cuu  1   tK uu  C vu   Cpu    1   tK uv   1   tK vv  1   tK pv  n   u    v  p   n 1   1   tK up    1   tK vp   0   n   u  fu       v   t  f v  p 0  n   n   (2.141)  Due to irregularities in the coefficient matrices K and C, such as zero values on the diagonals and huge numerical difference between the terms, care needs to be exercised in the selection of the parameter θ for the solution to be stable in time. In all the numerical examples discussed in this work, the backward Euler approach was used and led to stable results in all the cases.  2.4.6. Non-linear solution scheme To implement the step-wise approach for the solution of the displacement vector, the component of the first term in (2.138) that includes Kuu i.e.  45  Chapter 2: Governing Equations and FE Model of Two-phase Media  K uu n un  (2.142)  is replaced by the step-wise description  f u int n  f u int n  K uu n un  un  fu K n  (2.143)  is an accumulated load vector at step n from previous time-steps obtained by the  following recursive relationship  f u int n1  f u int n  K uu n1 un1  un   (2.144)  Let us choose the backward Euler scheme (   1 ) as the method of solution in time domain. As a result, the first row in (2.128) at t=tn+1 may be written in the form  f u int n  K uu n1 un1  un   K uv n1 v n1  K up n1 pn1  Cuu n1 u n1  f u  (2.145)  At every time step, an iterative solution is performed. At the kth iteration, the matrix equation of the system is  C uu  tK uu  C vu   C pu   tK uv tK vv tK pv  tK up   u nk11  u nk1   f u n 1  f u int n     tK vp   v kn 11    f v n1  C vu u n  u nk1 0   p nk11   C pu u n  u nk1 n 1 k           (2.146)      2.5. Treatment of Discontinuity in Volume Fraction In some cases, there may be a jump in the distribution of fluid volume fractions. This could be due to the placement of two layers of porous structure with different compaction amounts beside each other. In the composite processing, such cases may arise when resin flows out of the fibre-bed and pools in a corner leading to the presence of a purely viscous resin adjacent to the composite system consisting of fibre-bed and resin. As we will observe in an example in the next chapter, the discontinuity could also arise as a result of fluid flow motivated by s specific loading situation. Here, we introduce the measures that are required in the FE formulation of such discontinuities. 46  Chapter 2: Governing Equations and FE Model of Two-phase Media  In the derivation of the weak form of the fluid phase equilibrium in the second equation of (2.61), integration by parts (or Gauss theorem) is applied on the pressure term as follows   P,i w p d    Pni w p d   Pw p ,i d          Pni w p d   Pw p ,i d   ,i Pw p d     (2.147)    We should remind the reader that the assumption for performing integration by parts in the general form of b  b   udv  uv   vdu b a  a  (2.148)  a  is that both u and v be continuously differentiable in the relevant domain. Equation (2.147) holds for any smooth distribution of volume fraction over the domain, but it does not hold in cases such as layered media with different volume fractions, where a sudden jump occurs in the volume fraction. Here, we assume that any jump in  values happens only at the boundary between elements, and calculate a general formula for the amount of this error in case of a jump in volume fraction. To this end, we assume a 1D case involving two elements with a jump in the distribution of volume fraction at their interface. The schematic diagram of the pressure and volume fraction functions in this example is depicted in Figure 2.3. Let us assume that the volume fraction functions depicted in the diagram are in fact equal to the multiplication of the actual volume fraction distribution by the FE shape functions wp. We calculate the value of the left-hand-side integral in (2.147) over the length of each of the two elements, and add them up to arrive at the actual value of the above-mentioned integral over the domain obtained without using integration by parts. The distribution of pressure along the domain is assumed to take the following form  P ( x) P( x)   1 P2 ( x)  xL  x  xd xd  x  xR  (2.149)  47  Chapter 2: Governing Equations and FE Model of Two-phase Media  with a jump equal to α at the interface of the two elements. Differentiating the pressure function along the length of the domain leads to   P1( x)  P( x)   ( x  xd ) P( x)  2  x L  x  xd x  xd  (2.150)  xd  x  x R  At this point, it is helpful to review the following formula regarding integration of Dirac delta function    f ( x) ( x  x  d  )dx  f ( xd )  (2.151)    Therefore, the actual value of the above-mentioned integral term obtained without integration by parts will take the form xR  xd  xR  xd  xL  xL  xd  x d  I1   Pdx   1 ( x) P1( x)dx   2 ( x) P2( x)dx    ( x) ( x  xd )dx  (2.152)  Taking advantage of (2.151), and considering the symmetry of δ function, we may conclude that xd  xR  xL  xd  I1   1 ( x) P1( x)dx   2 ( x) P2( x)dx    2  1 ( xd )  2 ( xd )  (2.153)  We now calculate the value of the same integral term that the current FE approach with the integration by parts would result xR  xd  xR  xL  xd  I 2   Pdx  PxRL    Pdx    Pdx x  xL  (2.154)  Expanding the above relationship leads to xd  xR  xL  xd  I 2  P   1( x) P1 ( x)dx   2 ( x) P2 ( x)dx xR xL  (2.155)  and applying integration by parts on the integral terms in the above, we have 48  Chapter 2: Governing Equations and FE Model of Two-phase Media xd  xR  xL  xd  I 2  P   P    P    1 ( x) P1( x)dx   2 ( x) P2( x)dx xR xL  xd 1 1 xL  xR 2 2 xd  (2.156)  After simplifications, we arrive at xd  xR  xL  xd  I 2   1 ( x) P1( x)dx   2 ( x) P2( x)dx  1 ( xd ) P1 ( xd )  2 ( xd ) P2 ( xd )  (2.157)  The difference between I1 and I2 represents the term that should be added to the FE formulation to account for any discontinuity of volume fraction distribution between elements. Considering that   P2 ( xd )  P1 ( xd ) and after simplifications, we arrive at the correction term for the 1D case  I1  I 2   1 P1 ( xd )  P2 ( xd )1 ( xd )  2 ( xd )  Pavg    x xd 2  (2.158)  where [[ ]] represents the jump of a function. Taking into account the point that in the above derivation of the correction term,  represents the multiplication of fluid volume fraction by shape functions wp, one may argue that the correction term for a general 2D or 3D problem takes the form of   P  n w d avg  i  p  (2.159)  d  for all inter-element boundaries, d where fluid volume fraction incurs a discontinuity. As a matter of fact, the above correction term is a force term that should be applied on any inter-element boundary where a jump in volume fraction occurs in order to avoid any discontinuity in the pressure distribution over the domain of the problem. If efficiency is of no concern, one may also implement this term into the code so that it is applied on all the inter-element boundaries and it would work just the same since if there is no jump in the volume fraction the above term will be equal to zero. As it is obvious from (2.159), this force term is a function of fluid pressure, a main degree of freedom of any two-phase problem, and therefore it is best to introduce this term not as an applied force but as some extra coefficients in the stiffness matrices of the relevant elements. Since these extra coefficients are related to the equilibrium equation of the fluid phase and are to be 49  Chapter 2: Governing Equations and FE Model of Two-phase Media  multiplied by pressure degrees of freedom, they contribute to the K vp components of the stiffness matrix. Assuming a general interface between the elements m1 and m2, we may rewrite (2.159) in the form            1 1 Pm1  Pm2  ni wp d   Pm1  Pm2 m1  m2 ni wp d  2 d 2 d  (2.160)  The corrections needed in K vp (or additional K vp matrices) for the two involved elements may then be defined as  K K    vpextra m 1    vpextra m 2       m1    d        N  d N  N  d    m2 N n  T  p  m1  m1  183  m1    m2  nT  (2.161)  p  m2  m2  d  183  The above integrals are estimated using 3 Gauss points along the specified boundary located at  or   0, 15 5 . N n represents the product of the shape functions by the normal vector to the boundary of each element, and is written as    N n  nx N1  ny N1  nx N 2  ny N 2  nx N9  n y N9    (2.162)  The normal vector on the boundary is in fact the outward normal to m1 element. The effect of this correction will be demonstrated here through numerical examples.  2.6. Interfaces with Purely Viscous Fluids The need to model the interface of a two-phase porous system with a purely viscous fluid arises quite often in the field of biomechanics, e.g. the lubrication mechanics of joints. Some of the work done on modeling such interfaces in biomechanics include Mow et al. (1980), Hou et al. (1989), Spilker and Maxian (1990), Almeida and Spilker (1997,1998), Ateshian and Wang (1998), and Chan et al. (2000). As was mentioned in the previous section, in the processing of composite laminates with curved geometry such interfaces may arise when resin flows out of the fibre-bed and pools in a space provided by the complex geometry. 50  Chapter 2: Governing Equations and FE Model of Two-phase Media  To model a viscous medium using the current two-phase approach, we assume fluid volume fraction of unity, and very small values for the moduli of the solid structure. At the boundary between a viscous fluid and two-phase medium, the shear force, which is carried only by the viscous fluid on one side, is transferred partly to the viscous fluid flowing through the porous medium on the other side and partly to the porous solid structure. The actual (phase-averaged) value of shear stress carried by the fluid on both sides is the same, but there is a jump in the volume-averaged values of the fluid shear stress at the boundary. Assuming that the ratio of the areas of the two phases on the boundary is the same as the ratio of their volumes (which is the case for isotropic materials), the boundary condition between the two media with distinct volume fractions may be described as [Hou et al. (1989)]:      ij 1  1    ij 2  (2.163)  2  Assuming that 1 represent the viscous fluid and therefore 1  1, one may remove the subscript from 2 and rewrite (2.163) to have        1   1   1   1 v ij 1  ij 2  ij 2  i, j   v j ,i 2  (2.164)  To apply the above condition to the FE approach, a virtual traction force equal to  1   1 vi , j  v j ,i 2 t j  (2.165)  is applied to the fluid phase at the boundary between the two media. t j represent the components of the tangent vector of the inter-media boundary. We assume a general interface between the elements m1 and m2, with the latter being the one representing a purely viscous fluid. Once implemented in the current FE approach, the virtual force in (2.165) contributes to the components of K vv matrix for element m1 as follows  K vv extra m  1     1   d  m1     B    1  B3  T  m1  1 m1  d  (2.166) 1818  where B 3 is defined by  51  Chapter 2: Governing Equations and FE Model of Two-phase Media   nx N1 0 nx N 2  B3   0 n y N1 0 n y N1 nx N1 n y N 2   0   nx N 9  ny N 2   0  nx N 2  n y N 9  0   ny N9  nx N 9  318  (2.167)  52  Chapter 2: Governing Equations and FE Model of Two-phase Media  2.7. Figures u / v nodes  Pressure nodes      y  x Figure 2.1: Schematic representation of Q2P-1 (9-3) element  u / v nodes  Pressure node      y  x Figure 2.2: Schematic representation of Q1P0 (4-1) element  53  Chapter 2: Governing Equations and FE Model of Two-phase Media  P,φ P1(x)  P  P2(x) α  φ1(x)  φ  φ2(x) β  xL  xd  xR  x  Figure 2.3: The schematic diagram of the pressure and volume fraction distribution for a 1-D system consisting of 2 elements with a jump at their interface  54  Chapter 3: Verification of the Two-phase Model  The FE approach presented in the previous chapter is implemented into a MATLAB code using 9-3 and 4-1 elements. These two elements are based on the u-v-p formulation. To validate the two-phase model, a few examples are solved numerically with the code and the results are compared to analytical solutions (if available) and other numerical predictions cited in the literature. Also, for the sake of more in-depth comparison, the u-p formulation is implemented in a 9-4 element which is rated to behave in a very similar fashion to the 8-4 standard element in ABAQUS for multi-phase porous materials. The naming criterion is based on the number of nodes used to discretize the kinematic degrees of freedom and the number of nodes assigned to the pressure of the fluid phase. In the u-p formulation, continuity of pressure must be satisfied. The examples chosen for the purpose of verification of the flow aspect of analysis include: (1) a column of porous material undergoing a uniform load at the top with permeable BC at the top surface, (2) a column of porous material undergoing a uniform load at and mid-height with impermeable BC all around, (3) flow in a channel with a porous wall, and (4) flow-compaction in processing of composite angle laminates.  3.1. General Two-phase Problems 3.1.1. Column of saturated porous medium under loading at the top Here, we consider an example to verify the predictions of the developed two-phase model and also demonstrate the phenomenon of initial oscillations in the pressure distribution. Figure 3.1 shows the schematic representation of a column of saturated porous medium with permeable BC at the top surface under vertical load at the top (left figure). The dimensions of the problem, boundary conditions, material properties and the loading conditions are depicted in Figure 3.1. The total load of f0 is applied from t = 0. Two FE meshes are assumed; one with 10 9-3 elements and another with 20 9-4 elements along 55  Chapter 3: Verification of the Two-phase Model  the height of the column. Note that 10 9-3 elements lead to 20 independent pressure DOFs along the height while 20 9-4 elements have 21 independent pressure DOFs along the height of the column. We will be comparing the pressure response in this example, and that is why the number of the 9-4 elements is chosen to be twice that of the 9-3 element. Figure 3.3 and Figure 3.4 depict the time-history of pressure distribution along the height of the column for a fluid viscosity equal to 5×103 Pa.s predicted respectively by 10 9-3 elements and 20 9-4 elements. Figure 3.5 to Figure 3.12 present the same data but for lower viscosities that gradually decrease from 5×103 Pa.s to 0.5 Pa.s by increments of one order of magnitude. By changing the viscosity of the fluid, we change the time scale in that higher values for viscosity represent the behaviour at earlier stages of compaction and vice versa. It is evident that in the case of 9-4 elements the oscillations happen at the very start of the analysis while the 9-3 element initially does not show any significant oscillations but they develop to a maximum value at a critical time. The values of the observed oscillations are much smaller in the case of 9-3 elements. There is an analytical model for the case of 1D consolidation that applies to the current example. Terzaghi (1943) presented a closed form solution for the distribution of fluid pressure along the height of the column as a function of time. That formula may be written in terms of the material parameters introduced in the current work as  P y, t    4 f0       1 j 1 cos  2 j  1y  exp  2 j  1  2 S y E t    2 j 1 j 1  2    2h       4h 2     (3.1)  where Sy is the permeability of the porous material along the height of the column and E represents the Young‟s modulus of porous structure. Figure 3.13 shows the pressure distribution at an early time during the consolidation of the column predicted by both the 9-3 and 9-4 elements, comparing them with the analytical solution in (3.1). Superiority of the 9-3 element‟s prediction is due to the fact that using this element there is no need to enforce the value of pressure (to zero here) at the permeable boundary. Figure 3.14 to Figure 3.17 present the pressure distribution at the same time but for cases with gradually decreasing viscosity respectively. These diagrams represent later stages of consolidation of the porous column. The oscillations 56  Chapter 3: Verification of the Two-phase Model  related to the 9-4 elements gradually vanish as we move further from the initial stages of compaction. In all the cases a very good match is observed between the distribution predicted by the FE code and the closed form solution. Let us focus on the vicinity of the permeable BC to investigate how the oscillations of pressure develops in the case of 9-3 elements. Figure 3.18 shows the pressure surface of the top two 9-3 elements at an early stage of compaction. Figure 3.19 presents the same data at a time that the oscillation is more pronounced. It is evident that when the point of maximum curvature of the analytical solution is located between the pressure nodes of the top element, the oscillations are more pronounced in the 9-3 elements. The u-p approach (implemented in the 9-4 element) in FE formulation of the governing equations of two-phase media is to substitute the equilibrium equation of the fluid phase (2nd equation in (2.62)) into the mass conservation(1st equation in (2.62)); leading to the following set of two equations: u i ,i  S ij P, j  ,i  0   1 Cijkl u k ,l  u l ,k  , j  P,i  0  2      (3.2)  As a result, the pressure undergoes a second-order differentiation in space, and interelement continuity of pressure becomes necessary. This is not the case in the present u-vp formulation where all three of the mass conservation equation, fluid equilibrium equation, , and the total equilibrium equation (3rd equation in (2.62)) are formulated into the finite element technique. Discontinuity of pressure gives the 9-3 element some advantages over the common elements in terms of stability and observed initial oscillations in pressure profile.  3.1.2. Column of saturated porous medium under loading at mid-height To demonstrate the effect of the treatment of jumps in the volume fraction, an example will be presented here in the form of a column of saturated incompressible porous medium under uniform loading at its mid-height with impermeable B.C. all around. The 57  Chapter 3: Verification of the Two-phase Model  dimensions of the problem, boundary conditions, material properties and the loading conditions are depicted in Figure 3.2. The total load of f0 is applied from t = 0. Resin flows from the bottom half to the top half of the column, therefore gradually causing a discontinuity in the volume fraction of the porous medium. Figure 3.20 and Figure 3.21 depict the prediction of the 9-3 element for the time-history of pressure distribution along the height of the column for μ= 5×10-2 Pa.s and μ= 5×10-3 Pa.s, respectively. The gradual formation of an unrealistic jump in the pressure profile as the flow of the fluid occurs is evident. Figure 3.21 shows a clear jump in the pressure distribution at the final state of equilibrium. Figure 3.22 and Figure 3.23 present the same predictions where no discontinuities are observed after taking the corrective measure discussed in the previous chapter regarding the jump in volume fraction.  3.1.3. Taylor and Beavers-Joseph problems Here we consider the unidirectional flow of a viscous fluid over and through a rigid porous wall due to a constant velocity applied to the top surface of the fluid. The constant velocity V0 is applied to the fluid via a moving rigid wall at the top boundary of the fluid layer as it is depicted in Figure 3.24(a). This problem has indeed a Couette-type behaviour as the only distinction between the two is the porous wall underneath the viscous fluid in the current configuration, and one would imagine if the porous wall happens to be impermeable enough, then unidirectional Couette flow is a limiting case of the problem at hand. Another problem of interest may be described by the unidirectional flow of a viscous fluid over and through a porous wall sustained by pressure gradient in a 2D channel (Figure 3.24(b)). This case may be considered as a Poiseuille-type flow, as its limiting case when the permeability of the porous structure diminishes to zero is in fact the Poiseuille flow. Hou et al. (1989) refer to the two above-mentioned problems as the Taylor problem [Taylor (1971)] and the Beavers-Joseph problem [Beavers and Joseph (1967)], respectively. They established a set of boundary conditions for such interface problems in their work. They also presented an analytical solution for fluid velocity profile for the Beavers-Joseph problem based on their choice of interfacial boundary conditions. Also,  58  Chapter 3: Verification of the Two-phase Model  the solution of Taylor problem was presented via graphs depicting the flow profile. However, their original differential equations were based on the fluid‟s phase-averaged velocity (absolute value of velocity) compared to volume-averaged fluid velocity in this work, and also drag coefficient for the porous medium as opposed to the permeability coefficient in this work. Here, we present the analytical solution for the two problems using the nomenclature of the current work. Considering the assumptions involved in the Taylor problem, the equilibrium equation of the fluid phase for the case of this unidirectional flow problem may be written as   2v   0,   z 2  2    v   v  0  S  z 2  0  z  h1  (3.3)  h2  z  0  where v is the volume-averaged velocity of the fluid in the x direction. The boundary conditions may be readily identified as v(h1)=V0, and v(-h2)=0. As for the interfacial boundary conditions between the viscous fluid at the top and the underlying two-phase system, we use the relationships derived by Hou et al. (1989)  v  v 1 v  z   z 0  v z  (3.4) z 0  Solving the set of differential equations in (3.3) with their assigned boundary conditions, and the specified interfacial conditions in (3.4) leads to the analytical solution for the volume-averaged velocity of fluid along z axis         cosh h2    S S     z  h1 , 0  z  h1 V0 1     sinh  h    h  cosh h     2 S 1 S  2 S        v      z  h2  V0 sinh    S    h2  z  0           h1  cosh h2   sinh  h2   S S S       (3.5)  59  Chapter 3: Verification of the Two-phase Model  In a very similar fashion to the Taylor problem, the equilibrium equations of the fluid phase for the case of Beavers-Joseph problem may be introduced as   2 v P   ,   z 2 x  2    v   v  P  S x  z 2  0  z  h1  (3.6)  h2  z  0  where v is the volume-averaged velocity of the fluid in the x direction. The boundary conditions may be identified as v(h1)=v(-h2)=0. The interfacial boundary conditions between the viscous fluid at the top and the underlying two-phase system are the same as the ones introduced in (3.4), since the system geometry of the two problems are identical. Solving the set of differential equations in (3.6) with the appropriate boundary and interfacial conditions leads to the analytical response for the volume-averaged velocity of fluid along z axis in the Beavers-Joseph problem             h2  h1 sinh  h2  2  S 1  cosh     S  S    1 P     z  h1  z   ,         2 x     h1   sinh  h2 cosh h2      S S S      v    2    h2   h2 S h1  S    1  e  h1 e    S 2S       S  P    sinh  z  h2   1  e    x          S      sinh h  h cosh h 2 1 2      S S S       0  z  h1      z  h2   S       h2  z  0  (3.7) The above result matches the analytical solution obtained by Hou et al. (1989) for the profile of the phase-averaged velocity of the viscous fluid. The closed-form solutions in (3.5) and (3.7) will be used here for verification of the approach presented in this work. Chan et al. (2000) compared the results of their FE approach with the analytical solutions for the Taylor problem, but they did not report any numerical results for the BeaversJoseph problem. We present numerical results for both of the above problems.  60  Chapter 3: Verification of the Two-phase Model  It is assumed that h1  h2  4mm , and viscosity of the fluid, μ, is equal to 5 Pa.s. Three different values for the permeability of the porous layer are assumed to capture the limits of velocity profile that could be obtained in these two problems. Two uniform FE meshes are considered; one with four 9-3 elements through the thickness, and another with eight. As these problems are 1D and we have constrained the degrees of freedom in the z direction, the number of elements in the x direction is irrelevant. Figure 3.25 to Figure 3.27 present the velocity profiles (for decreasing values of permeability of the porous wall) obtained by the current approach for the Taylor problem and compares them with the analytical solution in (3.5). The convergence of the results even with these relatively small number of elements is clear, and an essentially perfect match is observed between the numerical and analytical solutions. Figure 3.28 to Figure 3.30 show the velocity profiles (for decreasing values of permeability of the porous wall) obtained by the current approach for the Beavers-Joseph problem and compares them with the close form solution in (3.7). In a similar fashion to the Taylor problem, convergence of the response is evident, and a very good match is observed between the numerical and analytical response. Chan et al. (2000) added a viscous shear term to the hydrostatic pressure to have a complete stress tensor for the fluid phase. Here, we compare their formulation with the current study. A very important distinction is in the differential equation pertaining to the equilibrium equation of the fluid phase. Rewriting Chan et al.‟s version of the fluid equilibrium equation using the notation of the current work for the sake of comparison, we arrive at a form equivalent to       P ,i  vi , j  ui , j   v j ,i  u j ,i  , j   S ij1v j  0  (3.8)  The major difference between the above and the fluid equilibrium equation in the current study in (2.58) is in the first term which includes a form of pressure gradient. Equation (2.58) used in the current approach is consistent with the governing equations of porous media presented in the prominent literature in this field (e.g. Lewis and Schreffler, 1998) while equation (3.8) is not. Furthermore, if one assumes that the velocity terms vanish in  61  Chapter 3: Verification of the Two-phase Model  (3.8), it would lead to an equilibrium state with variable pressure for a case with variable volume fractions throughout the domain, an unrealistic and undesirable outcome. Another difference between the work by Chan et al. and the current approach is that they based their set of differential equations on the phase-averaged velocity of fluid while the volume-averaged relative velocity of the fluid is the velocity parameter chosen in this work. Not only does the latter choice lead to a very straight-forward presentation of the mass conservation equation as in (2.57), but also it allows for the continuity of the fluid velocity degrees of freedom at the interface of a purely viscous fluid and a two-phase medium. However, the use of phase-averaged fluid velocity in the formulation presented by Chan et al. (2000) leads a discontinuity in the velocity field at the interface which requires special assembly considerations that are not desirable in a robust code.  3.2. Flow-compaction in Processing of Composite Angle Laminates Hubert (1996) and Hubert et al. (1999) performed simulations involving resin flow and compaction of angle shaped laminates. Hubert et al. developed a 4-4 element in COMPRO based on the u-p formulation. Of two-phase media. They carried out a comprehensive parametric study of the numerical response and the sensitivity of the effect of various constitutive properties on the compaction behaviour of the laminates. Hubert (1996) also performed an experimental study on the compaction of unidirectional angle laminates on convex and concave tools. The samples were made of two different materials including AS4/3501-6 and AS4/8552. He numerically simulated examples closely representative of the geometry of the actual specimens made for the experiments. The change in the thickness of the final products at the corner and along the flat part were studied and compared with the results predicted by the simulations. The numerical predictions of the final thickness changes for the specimens made of AS4/3501-6 compared much better than AS4/8552 specimens with those obtained from the corresponding experiments. Here, the 9-3 element is used to model the same problem  62  Chapter 3: Verification of the Two-phase Model  with the same material properties assumed in Hubert‟s work to compare the results from the two approaches (See Table 3.1 and Table 3.2). Figure 3.31 shows the geometry and boundary conditions of the problem assuming a convex tool and the time-history of autoclave conditions. Since the current code does not deal with heat transfer analysis, we have assumed that the temperature of the laminate changes uniformly throughout the process and follows the temperature history of point A obtained from Hubert‟s numerical run for the same problem. Two different FE meshes are presented in Figure 3.32 to ensure convergence of the results. Point A is located on the top corner of the laminate, and point B is located on the top surface of the laminate at the mid-length of the flat section. The normal displacements at A and B are compared to those results reported by Hubert. Figure 3.33 shows the time-history of the normal displacement at points A and B for the case of the AS4/3501-6 [0°] unidirectional angle laminate on a convex tool and compares the response with those reported by Hubert. Figure 3.34 presents the same comparison for the case of AS4/8552 [0°] unidirectional angle laminate. The convergence of the response is verified and a good match is observed between the current predictions and those reported by Hubert for both materials. Any observed difference in the history of normal displacements especially in the case of AS4/3501-6 material is attributed to the fact that in the present work the temperature is assumed to be uniform over the domain. It was observed that for cases involving [90°] unidirectional angle laminates, the current predictions did not match the response reported by Hubert. Investigating further into this matter made it clear that in the runs for [90°] cases, Hubert‟s code had erroneously set the permeability matrix equal to those of a [0°] laminate. In a modified version of COMPRO, we fixed that error and reanalyzed the [90°] cases. The response histories that will be presented in Figure 3.35 and Figure 3.36 as Hubert‟s are in fact the results obtained by the updated version of COMPRO. Figure 3.35 presents the time-history of the normal displacement at points A and B for the case of the AS4/3501-6 [90°] unidirectional angle laminate on a convex tool and compares the response with those reported by our version of Hubert‟s approach. Figure  63  Chapter 3: Verification of the Two-phase Model  3.36 shows the same comparison for the case of convex AS4/8552 [90°] unidirectional angle laminate. In both of the above figures, the convergence of the response is clear and a good match is obtained between the displacement history predicted by the 9-3 element and the response obtained by the corrected version of Hubert‟s approach. The 4-4 element developed by Hubert et al. (1999) approximates the displacement and pressure degrees of freedom in the same fashion by continuous linear polynomials in each direction. Due to having the same level of discretization for displacement and pressure, the 4-4 element sustains instability issues and also its accuracy is conditional. Vermeer and Verruijt (1981) presented an accuracy condition for such finite elements. The condition sets a minimum size for the time steps or a maximum for mesh size to ensure accuracy and avoid excessive pressure oscillations. This issue led Hubert et al. to run the COMPRO flow module once the viscosity of resin fell below 1000 Pa.s. In the present work, there are no such limitations but since the goal here was to compare the results of the 9-3 with the predictions of Hubert et al., for viscosity values above 1000 Pa.s the resin viscosity was artificially increased to 1×107 Pa.s to prevent any resin flow. This explains the plateau observed at the beginning of the displacement time-history predictions in Figure 3.33 to Figure 3.36.  64  Chapter 3: Verification of the Two-phase Model  3.3. Tables Table 3.1: Resin and fibre-bed properties for the AS4/3501-6 angle laminate  Resin viscosity and fibre-bed permeability  Resin degree of cure  d  dt  C1  C2 1   B     d   C3 1     dt  CI  AI e  EI    0.3  Fibre-bed elastic properties     exp U / RT exp   100GPa D    Symm.    4.6 10 17 Pa.s    0.3  U  114477 J / mol    14.8  RT  S1   A1  3.5017  10 7 s 1 A2  3.3567  10 7 s 1  S2   A3  3.2667  103 s 1 E1  80700 J mol  rf  rf  2  3  2  Vf   V V a  f  a  f  4k   V  V   0.67 MPa  3.33MPa   E3  15.15MPa 41.03MPa    3MPa    3  1   1  Va  0.81  E3  56600 J mol  0  f  Vf  1  E2  77800 J mol     0.5MPa  0  1  V   4k 2  0  E3  3    0.1   3  0  0.16   3  0.1  0.2   3  0.16   3  0.2 3  0  k   0.2 rf  4m  B  0.47  Table 3.2: Resin and fibre-bed properties for the AS4/8552 angle laminate  Resin degree of cure  Resin viscosity  Fibre-bed permeability  A B d Ae  Ea RT  m 1     g     dt 1  eC   C 0  CT T    A exp E RT   g     n    Ea  66500 J mol  A  3.45 10 10 Pa.s  m  0.813  E  76536 J / mol  n  2.74  A  3.8 B  2.5   C 0  5.48 10  CT  0.19 / C  3   g  0.47  S2   2  rf  2  4k   1  V   3  f  4k    A  1.53  105 s 1  C  43.1  rf  S1     Fibre-bed elastic properties  Vf  2    Va V f  1  V  V a  Vf  1  Va  0.68  k   0.2 rf  4m  f   1  3  100GPa D    Symm. 2.08MPa 8.65MPa   E3   22 MPa 62.1MPa    6.6 MPa  0  E3  3      0.5MPa   0 0   0.035   3  0  0.079   3  0.035  0.1   3  0.079   3  0.1 3  0  65  Chapter 3: Verification of the Two-phase Model  3.4. Figures y  f0 h = 25 mm b = 2.5 mm υ0 = 0.5 E = 70 MPa ν=0 f0 = 700 kPa Sy = 4.5×10-14 m2  h  x b  Figure 3.1: A column of saturated porous medium under uniform loading at the top, and meshes of 10 9-3 elements and 20 9-4 elements  y  h = 25 mm b = 2.5 mm υ0 = 0.5 E = 70 MPa ν=0 f0 = 700 kPa Sy = 4.5×10-14 m2  h2  f0  h2  x b  Figure 3.2: A column of saturated porous medium with impermeable BC all around with uniform loading at mid-height, and mesh of 10 9-3 elements 66  Chapter 3: Verification of the Two-phase Model  Figure 3.3: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=5×103Pa.s  Figure 3.4: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=5×103 Pa.s 67  Chapter 3: Verification of the Two-phase Model  Figure 3.5: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=5×102 Pa.s  Figure 3.6: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=5×102 Pa.s 68  Chapter 3: Verification of the Two-phase Model  Figure 3.7: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=50 Pa.s  Figure 3.8: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=50 Pa.s 69  Chapter 3: Verification of the Two-phase Model  Figure 3.9: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=5 Pa.s  Figure 3.10: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=5 Pa.s 70  Chapter 3: Verification of the Two-phase Model  Figure 3.11: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-3 (u-v-p) elements. μ=0.5 Pa.s  Figure 3.12: Time-history of pressure distribution along the column with top permeable BC and under applied load at the top, obtained by 9-4 (u-p) elements. μ=0.5 Pa.s 71  Chapter 3: Verification of the Two-phase Model  9.E+05  8.E+05  7.E+05  P (Pa)  6.E+05  5.E+05 closed-form solution  4.E+05  FE: 10 9-3 elements  3.E+05  FE: 20 9-4 elements  2.E+05  1.E+05  0.E+00 0.5  0.55  0.6  0.65  0.7  0.75  0.8  0.85  0.9  0.95  1  y/h  Figure 3.13: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=5×103 Pa.s  9.E+05  8.E+05  7.E+05  P (Pa)  6.E+05  5.E+05 closed-form solution  4.E+05  FE: 10 9-3 elements  3.E+05  FE: 20 9-4 elements  2.E+05  1.E+05  0.E+00 0.5  0.55  0.6  0.65  0.7  0.75  0.8  0.85  0.9  0.95  1  y/h  Figure 3.14: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=5×102 Pa.s 72  Chapter 3: Verification of the Two-phase Model  8.E+05  7.E+05  6.E+05  P (Pa)  5.E+05  4.E+05 closed-form solution FE: 10 9-3 elements  3.E+05  FE: 20 9-4 elements  2.E+05  1.E+05  0.E+00 0.5  0.55  0.6  0.65  0.7  0.75  0.8  0.85  0.9  0.95  1  y/h  Figure 3.15: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=50 Pa.s  8.E+05  7.E+05  6.E+05  P (Pa)  5.E+05  4.E+05 closed-form solution FE: 10 9-3 elements  3.E+05  FE: 20 9-4 elements  2.E+05  1.E+05  0.E+00 0.5  0.55  0.6  0.65  0.7  0.75  0.8  0.85  0.9  0.95  1  y/h  Figure 3.16: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=5 Pa.s 73  Chapter 3: Verification of the Two-phase Model  8.E+05  7.E+05  6.E+05  P (Pa)  5.E+05  4.E+05  closed-form solution  3.E+05  FE: 10 9-3 elements  2.E+05  FE: 20 9-4 elements  1.E+05  0.E+00 0.5  0.55  0.6  0.65  0.7  0.75  0.8  0.85  0.9  0.95  1  y/h  Figure 3.17: Pressure distribution at t=3 s, along the column with top permeable BC and under applied load at the top. μ=0.5 Pa.s  9.E+05  Element #9  8.E+05  Element #10  7.E+05  P (Pa)  6.E+05  5.E+05 closed-form solution  4.E+05  FE: 9-3 elements  3.E+05  pressure surface of element  2.E+05  1.E+05  0.E+00 0.8  0.82  0.84  0.86  0.88  0.9  0.92  0.94  0.96  0.98  1  y/h  Figure 3.18: Pressure surfaces of the two end elements at a very early stage of consolidation of the column 74  Chapter 3: Verification of the Two-phase Model  1.E+06  Element #9  Element #10  9.E+05 8.E+05 7.E+05  P (Pa)  6.E+05 5.E+05 closed-form solution  4.E+05  FE: 9-3 elements pressure surface of element  3.E+05 2.E+05 1.E+05 0.E+00 0.8  0.82  0.84  0.86  0.88  0.9  0.92  0.94  0.96  0.98  1  y/h  Figure 3.19: Pressure surfaces of the two end elements at a time that oscillation is wellpronounced  Figure 3.20: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained without the correction for discontinuity of volume fraction. μ=5×10-2 Pa.s 75  Chapter 3: Verification of the Two-phase Model  Figure 3.21: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained without the correction for discontinuity of volume fraction. μ=5×10-3 Pa.s  Figure 3.22: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained with the correction for discontinuity of volume fraction. μ=5×10-2 Pa.s  76  Chapter 3: Verification of the Two-phase Model  Figure 3.23: Time-history of pressure distribution along the column with impermeable BC under applied load at the mid-height, obtained with the correction for discontinuity of volume fraction. μ=5×10-3 Pa.s  z  V0 P 0 x  h1  (a) h2  x Porous wall  z  P 0 x  h1  (b) h2  x Porous wall  Figure 3.24: Schematic diagram of two unidirectional flow problems over and through a porous wall; (a) Taylor problem, (b) Beavers-Joseph problem  77  Chapter 3: Verification of the Two-phase Model  1 0.8 0.6 0.4  z/h  0.2 Analytical  0  FE: 4 elements FE: 8 elements  -0.2 -0.4 -0.6 -0.8 -1 0  0.2  0.4  0.6  0.8  1  1.2  Velocity (mm/s)  Figure 3.25: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10-6 m2  1 0.8 0.6 0.4  z/h  0.2 0 -0.2 Analytical  -0.4  FE: 4 elements FE: 8 elements  -0.6 -0.8 -1 -0.1  0.1  0.3  0.5  0.7  0.9  1.1  Velocity (mm/s)  Figure 3.26: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10-8 m2 78  Chapter 3: Verification of the Two-phase Model  1 0.8 0.6 0.4  z/h  0.2 0 -0.2 Analytical  -0.4  FE: 4 elements FE: 8 elements  -0.6 -0.8 -1 -0.1  0.1  0.3  0.5  0.7  0.9  1.1  Velocity (mm/s)  Figure 3.27: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10-10 m2  1 Analytical FE: 4 elements FE: 8 elements  0.8 0.6 0.4  z/h  0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0  0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2  Velocity (mm/s)  Figure 3.28: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10-6 m2 79  Chapter 3: Verification of the Two-phase Model  1 0.8 0.6 0.4  z/h  0.2 0 -0.2 Analytical  -0.4  FE: 4 elements FE: 8 elements  -0.6 -0.8 -1 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  Velocity (mm/s)  Figure 3.29: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10-8 m2  1 0.8 0.6 0.4  z/h  0.2 0 -0.2 Analytical  -0.4  FE: 4 elements FE: 8 elements  -0.6 -0.8 -1 -0.1  0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  Velocity (mm/s)  Figure 3.30: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10-10 m2 80  Chapter 3: Verification of the Two-phase Model  f H = 4.96 mm for AS4/3501-6 4.80 mm for AS4/8552 Permeable B.C.  L = 63 mm R = 4.57 mm  f L 163°C 2.5 hrs  R  2°C/min  T&P  Permeable B.C.  122°C 70 min  2.5°C/min  2°C/min  z  Autoclave pressure: 540 kPa  H  Bag pressure: 0 kPa  x  25°C  time  Figure 3.31: Geometry, BC, and processing cycle used for compaction of a unidirectional angle laminate on a convex tool  B  A  B  (a)  A  (b)  Figure 3.32: (a) 8×4 mesh and (b) 16×6 mesh of the angle laminate using symmetry  81  Chapter 3: Verification of the Two-phase Model  0  un (mm)  -0.2  -0.4  @ A - Hubert @ B - Hubert @ A - mesh (a)  -0.6  @ B - mesh (a) @ A - mesh (b) @ B - mesh (b)  -0.8  -1 0  10  20  30  40 Time (min)  50  60  70  80  Figure 3.33: Time-history of normal displacement for unidirectional AS4/3501-6 angle with [0°] fibres on a convex tool  0  un (mm)  -0.05  -0.1  @ A - Hubert @ B - Hubert @ A - mesh (a)  -0.15  @ B - mesh (a) @ A - mesh (b) @ B - mesh (b)  -0.2  -0.25 0  20  40  60  80 100 Time (min)  120  140  160  Figure 3.34: Time-history of normal displacement for unidirectional AS4/8552 angle with [0°] fibres on a convex tool 82  Chapter 3: Verification of the Two-phase Model  0  un (mm)  -0.2  -0.4  @ A - Hubert @ B - Hubert @ A - mesh (a)  -0.6  @ B - mesh (a) @ A - mesh (b) @ B - mesh (b)  -0.8  -1 0  10  20  30  40 Time (min)  50  60  70  80  Figure 3.35: Time-history of normal displacement for unidirectional AS4/3501-6 angle with [90°] fibres on a convex tool  0  un (mm)  -0.05  -0.1  @ A - Hubert @ B - Hubert @ A - mesh (a)  -0.15  @ B - mesh (a) @ A - mesh (b) @ B - mesh (b)  -0.2  -0.25 0  20  40  60  80 100 Time (min)  120  140  160  Figure 3.36: Time-history of normal displacement for unidirectional AS4/8552 angle with [90°] fibres on a convex tool 83  Chapter 4: Integration of Modeling from Fluid to Cured Resin  Having verified our FE approach in modeling porous flow in the previous chapter, we embark on identification and implementation of necessary modifications in the two-phase model in order to capture stress development as well as resin flow. To that end, we review the pseudo-viscoelastic stress model developed by Zobeiry et al. (2010) and present the formulation for adapting that into the current approach.  4.1. Pseudo-viscoelastic Stress Model In the pseudo-viscoelastic or PVE formulation [Zobeiry (2006), Zobeiry et al. (2010)], it is assumed that the moduli of resin are constant at any instant of time, but changing as a function of temperature and degree of cure. The moduli however, are derived so that the response will be approximately equivalent to a fully viscoelastic approach. In order to do that, Zobeiry et al. (2010) compared the CHILE (cure hardening instantaneously linear elastic) model [Johnston et al. (2001)]with a viscoelastic model. In their work, Zobeiry et al. (2010) presented the following mathematical formulation for stress development in 1D: t   (t )   E t e , T ( ),  ( )  0  d d d  (4.1)  where T is temperature, α is the degree of cure, and t e is an effective value of time that is defined as  log e  c1 ( f )m a ( ) log e  te    t f    T c1 ( f )m aT (t f )  te      tf   tf  (4.2)  84  Chapter 4: Integration of Modeling from Fluid to Cured Resin  where t f is the time at the onset of cool-down,  f is the degree of cure at the onset of cool-down , m is the rate of cool-down, and aT is the shift factor that relates the reduced time, ξ, to the time variable by the following equation t    0  1 dt aT  (4.3)  The shift factor, aT , is given by  logaT   c1  T  c2    (4.4)  where the coefficients c1 and c 2 are defined as   1  c1    a1 exp    a2   1  (4.5)  c2    T 0 c1    (4.6)  In the case of 3501-6 epoxy resin; a1  1.4 / C , a2  0.0712 / C , and T 0  30C [Kim and White (1996)]. T 0 is the reference temperature. The instantaneous elastic moduli for the resin (namely the shear and bulk moduli, G and K) are obtained from the following equations    W      te   exp   aT      (4.7)    W    te   exp   aT      (4.8)  Gm  G r  G u  G r  Km  K r  K u  K r  n   1  n   1    where superscripts r and u denote the fully relaxed and un-relaxed moduli of resin, respectively.   and W are the stress relaxation times and their associated weight factors for the  th Maxwell element in a generalized Maxwell model consisting of n elements. The relaxed and un-relaxed values of moduli for the 3501-6 resin are assumed to be 85  Chapter 4: Integration of Modeling from Fluid to Cured Resin  independent of degree of cure and temperature. For this material, the weight factors are also assumed to be constant during cure [Kim and White (1996)]. The stress relaxation times change with the degree of cure according to the following relationship       log  ( )  log   ( 0 )  f ( )  (   0 ) log(  )    (4.9)  where  0 is the reference degree of cure (equal to 0.98 in the case of 3501-6 resin). Also,  f   for this resin is defined as f    9.3694  0.6089  9.1347 2  (4.10)   is given by      p  0      0   (4.11)  where  p is the peak relaxation time. The relaxation times and their associated weight factors for 3501-6 resin are obtained from Kim and White (1996) and presented in Table 4.1.  4.2. Constitutive Equations In this section, we will review the elastic constitutive equations for isotropic and transversely isotropic materials.  4.2.1. Isotropic materials The elastic constitutive equation for 2D isotropic materials is written in the form of ζ  2Gε G  Kε K  (4.12)  in term of stress and strain tensors, ζ and ε, or     G   K  G G  K  K  (4.13)  86  Chapter 4: Integration of Modeling from Fluid to Cured Resin  in vector form where a single underscore denotes a vector quantity. Assuming pseudoviscoelastic formulation for the behaviour of resin, leads to a PVE representation for the composite material as well. Thus, at any time-step the behaviour of the composite is elastic and the rates of the two components of stress may be written as   G  G G  (4.14)   K  K  K  (4.15)  Summation of the above two equations leads to    GG  K  K  (4.16)  which could also be written in the form    D  (4.17)  where D is the elastic constitutive stiffness matrix (note that double underscore denotes a matrix quantity) given by   K  4G 3 K  2G 3 K  2G 3 0 0 0   K  4G 3 K  2G 3 0 0 0    K  4G 3 0 0 0  D  G 0 0   Symm. G 0   G    (4.18)  4.2.2. Transversely isotropic materials We review the formulation presented by Zobeiry (2006) for the constitutive behaviour of transversely isotropic materials. Based on a system of decoupled moduli presented by Hashin (1972) for use in micromechanics, the constitutive equation is written as ζ  ε  ε   K 2ε K2  G12εG12  G23εG23  (4.19)  or 87  Chapter 4: Integration of Modeling from Fluid to Cured Resin 5  ζ  ζ  ζ   ζ K 2  ζ G12  ζ G23   ζ p  (4.20)  p 1  In the expanded form , we may rewrite (4.19) in the following form 0 0  11 0 0  22   33 0 0  0      ζ    0 0 0    0 11 0   K 2 0  22   33 0   0 0 0  0 0 0 11  0  22   33  0 0 212 213  0   0     G23 0  22   33 2 23 0 0    G12 212 0 213 2 23   22   33  0 0   (4.21)  where K 2 is the in-plane bulk modulus of the material in its plane of isotropy (2-3 plane). It may be calculated by the following relationship  K2   1 4 1 4 12   E G23 E1  2  (4.22)  It can also be written in a format similar to the relationship for the in-plane bulk modulus of an isotropic material as  K2   E 21   2 12 21   (4.23)  where E  E2  E3 is the transverse modulus, and    23 is the Poisson‟s ratio in the plane of isotropy. G12 is the axial shear modulus, and G23 is the transverse shear modulus (shear modulus in the plane of isotropy). The two new moduli  and  are defined respectively as    E1  4K 2 122  (4.24)    2K 2 12  (4.25)  In the case of pseudo-viscoelastic formulation, the material behaviour at any time-step may be written in the form of  88  Chapter 4: Integration of Modeling from Fluid to Cured Resin             K  K 2  K 2   G  G23  G 23  (4.26)  2  23   G  G12  G 12  12  Summing the equations in (4.26) and taking into account (4.20), we may write         K 2  K  G23  G  G12  G 2  23  12  (4.27)  which could also be written in the form    D  (4.28)  where     D          0  0  K 2  G23  K 2  G23  0  0  K 2  G23  0  0  G23  0  Symm.  G12  0  0  0   0  0   G12   (4.29)  In this work, the material properties of the composite material are obtained from the equations of micromechanics adapted from Bogetti and Gillespie (1992).  4.3. Two-phase Model for Stress Development in Cured Materials An important aspect of this work is to implement the modifications required in the twophase element so that it can model the different stages of the cure of the material from a fluid resin to a completely cured polymer matrix. It is naturally expected of the element to successfully model the initial stages of the process consisting of resin flow and compaction of the sample. This was verified in the previous chapter, and also will be discussed in examples later in this work.  89  Chapter 4: Integration of Modeling from Fluid to Cured Resin  4.3.1. Consistent compressibility in mass conservation equation In order for the two-phase element to successfully model the cured composite, a few critical implementations are required. Many times in the modeling of flow, satisfactory results are obtained by assuming that the solid particles and the fluid phase (or one of them) are incompressible. As a matter of fact, Hubert et al. (1999) made both of the above-mentioned assumptions in his flow-compaction model. However, in the analysis of solid materials such assumptions would introduce considerable errors in obtaining the strain and stress fields. Therefore, compressibility of the phases needs to be considered in the governing equations implemented in the finite element scheme. Assuming small strains, the equations of mass conservation for compressible fibres and a compressible resin matrix may be written as [Lewis and Schrefler (1998)]    1    f   1     v c  0 t  f t  (4.30)  and     m     v flow    v c  0 t  m t  (4.31)  respectively, where ρf and ρm denote the density of fibres and resin matrix. vc is the velocity vector of the composite system defined in (2.51) and vflow is the volumeaveraged relative velocity (or flow velocity) of the resin given by (2.52). In the derivation of the above, it is assumed that the spatial gradients of the density of the phases are negligible. The summation of (4.30) and (4.31) results in the mass conservation equation of the system   v c    v flow   1    f   m  0  f t  m t  (4.32)  which may be rewritten in the indicial form of u i ,i  vi ,i   1    f   m  0  f t  m t  (4.33)  90  Chapter 4: Integration of Modeling from Fluid to Cured Resin  where the following notation is used for simplicity ui  vc i ,  vi  v flowi  (4.34)  Lewis and Schrefler (1998) introduce the following relationships for the mechanical components of the last two terms of the above equation  1  m 1 Pm   m t K m t  (4.35)  1  f 1 Pf 1 tr ζ     f t K f t  ii 1   K f t  (4.36)  and  where K m and K f are bulk moduli of resin and fibre respectively. ζ , the effective stress tensor related to the skeleton of the porous media (in this case; the fibre-bed), is defined by ζ  ζ  PI  (4.37)  where it is assumed that P  Pm  Pf . Here, it is worthwhile to note that if the solid grains (in this case; individual fibres) are compressible, we introduce “Biot effective” stress tensor as [Biot (1941), Terzaghi (1943)] ζ  ζ  bPI  (4.38)  where b is Biot coefficient and is defined by b  1  K fb Kf  (4.39)  where K fb is the bulk modulus of the fibre-bed (solid skeleton in general). Lewis and Schrefler (1998) introduce the following constitutive equation  91  Chapter 4: Integration of Modeling from Fluid to Cured Resin   tr ζ  1 P    ii K fb  u i ,i   t K f t    (4.40)  which is then substituted in (4.36) to result in 1  f 1 P 1  b  1 P  1 b 1   u i ,i   u i ,i     f t K f t 1    K f t  1 Kf   b    P    1    t  (4.41)  Incorporating (4.35) and (4.41) into (4.33) leads to the following mass conservation equation used in the treatment of porous soils:  b    P bu i ,i  vi ,i    0  K K m  t  f  (4.42)  In this work, we aim to introduce the micromechanics of the interaction of the phases and their load sharing settlement at the latest stage possible. To that end, we introduce the rate of volumetric strain of the two-phase system as  v s    1    f   m   f t  m t  (4.43)  which enables us to write the mass conservation equation in (4.33) in the following compact form ui ,i  vi ,i  v s  0  (4.44)  In the presence of thermochemical strains, the rate of volumetric strain is divided into mechanical and free components and the above equation is rewritten as mech free ui ,i  vi ,i  v s  v s 0  (4.45)  Having augmented the contributions of the phases into a single value for the volumetric strain rate of the composite, we are then able to attribute a constitutive equation to v c mech that is consistent with the micromechanical formulation that we opt for at the extreme end of cured composite. Assuming that the phases and the composite material are isotropic, a  92  Chapter 4: Integration of Modeling from Fluid to Cured Resin  general micromechanics formulation may be written in the form of the following two equations Gc  f G Gm , G f ,    (4.46)  K c  f K K m , K f ,    (4.47)  where the index c refers to the properties of the composite material. The constitutive equation for the bulk behaviour of the composite in 3D would be simply written in the form  v s mech   1   tr ζ s  1   tr ζ s      Kc t   ii  Kc t  3   (4.48)  where ζ s is the total stress of the two-phase system, and basically equal to the total stress  ζ in (4.38). In a format more familiar to the mass conservation equations frequent in the formulation of porous soils, one may write   1       tr ζ      bP    K  K m  t   ii   f  v s mech    (4.49)  The above equation is mentioned here just for the sake of comparison between the consistent treatment of the rate of volumetric strain presented here and the usual poromechanics formulations. The very small discrepancy that the consistent (with the solid micromechanics) approach would introduce into the model‟s poro-elasticity aspect of modeling, is a reasonable penalty to pay considering the gains in the consistency of the results for the analysis of solid composite with cured resin. In 2D plane strain however, the constitutive equation for the rate of volumetric strain needs to be expressed in a slightly different way than in (4.48). The reason is that in 2D plane strain (or plane stress for that matter!), the stiffness coefficient relating the volumetric strain and tr ζ  ii is not equal to K c . For 2D plane strain, the constitutive equation for the rate of the volumetric strain is  93  Chapter 4: Integration of Modeling from Fluid to Cured Resin  v s mech   1   tr ζ s  1   tr ζ s       K c  Gc 3 t   ii  Kc  Gc 3 t  2   (4.50)  4.3.2. The concept of modified effective stress (‘s-p’ stress) In the treatment of flow through porous media, an equilibrium equation is written for each of the solid and fluid phases that are related together with a common drag force. The equilibrium equation for the fluid resin, which was presented in (2.22), may be written in indicial form as  P,i   m ij , j  f d i  0  (4.51)  where f d i are the components of the drag force that, if defined as follows, leads to the Darcy‟s law. f d i   S ij1v j  (4.52)  and  m ij represent the deviatoric stress components of the resin matrix. Here, it should be reminded that, in the current notations, Darcy‟s law may be written as  P,i  f d i  0  (4.53)  while the resin equilibrium equation in (4.51) that we intend to use here, is the DarcyBrinkman (or Brinkman) equation [Brinkman (1949)]. The second term in (4.51) (i.e. the divergence of the deviatoric stress of the resin matrix) is usually much smaller than the drag force, so in most applications, it is enough to incorporate (4.53) in the governing equations [Tucker and Dessenberger (1994)]. However, in this work we aim for the porous flow equations to be able to model the behaviour of the cured composite at one extreme. Having that in mind, we choose to use (4.51) as the resin equilibrium equation. This will lead to a total equilibrium equation with all components of stress for the resin phase, which is consistent with the behaviour of the cured composite. We will later see that once any micromechanics scheme is incorporated in the total equilibrium equations  94  Chapter 4: Integration of Modeling from Fluid to Cured Resin  of the system, the term containing the deviatoric stresses of resin will not appear explicitly anymore, and it will be implicitly considered in a decomposition of total stress. The equilibrium equation for the solid skeleton (here; fibre-bed) may be written in the form σ ij , j  f d i  1   P,i  0  (4.54)  where ζ is the effective stress tensor of the fibre-bed. Summation of the two equilibrium equations in (4.51) and (4.54) leads to the total equilibrium equation of the system σ ij , j   m ij , j  P,i  0  (4.55)  which could also be written in terms of Biot effective stress, to yield a more general form covering for the compressibility of solid particles σ ij, j   m ij , j  bP,i  0  (4.56)  We may now revisit the governing equations in(2.61) and rewrite them in the form  u i ,i  vi ,i  v s  0   P,i   m ij , j  f d i  0  σ    m ij , j  bP,i  0  ij , j  (4.57)  As was mentioned earlier in this section, in most cases the term involving the deviatoric stress of the matrix (i.e.  m ij ) are negligible compared to f d i in the matrix phase equilibrium equation (second equation of the above). Based on the total equilibrium equation in (4.56) the total stress of the system may be defined as   s ij  σij   mij  bP ij  (4.58)   s  p ij  σ s ij  bP ij  σij   mij  (4.59)  Let us define the „s-p‟ stress as  Combining the deviatoric components of matrix stress with the stress components of the porous structure paves the way toward a formulation for the porous two-phase system 95  Chapter 4: Integration of Modeling from Fluid to Cured Resin  that is consistent with the micromechanics representation of choice for the final solid composite material. Considering the above points, we may rewrite (4.57) as ui ,i  vi ,i  v  0 s    P  f  ,i di  0 σ  bP,i  0  s  p ij , j  (4.60)  Assuming elastic behaviour for the fibre-bed, a constitutive relationship for Biot effective stress tensor may be written as  ij  C fb ijkl  kl  (4.61)  where C fb ijkl represent the stiffness moduli of the fibre-bed. Equation (4.56) should be representative of the equilibrium state of the composite material throughout the different stages of processing during which the resin matrix cures and goes through thermochemical changes from fluid to a solid cross-linked state. Conceptually, this new stress tensor  s  p ij can be looked at from two different viewpoints. Chronologically, the first viewpoint should be based on the porous flow side of the behaviour, and is explained by considering  s  p ij as an “improved” fibre-bed stress since the shear stress contribution from the resin is added to Biot effective stress. Therefore, one could observe that a very viscous resin would naturally contribute to “non-bulk” components of the modified fibre-bed stiffness. The second viewpoint is from the other extreme of the process, i.e. the solid micromechanics of the composite. Here,  s  p ij is viewed as (as the symbol suggests presently!) the components of the composite stress excluding a share of the total pressure. This suggests that any micromechanical model that we would use for the regular solid composite could be introduced for the components of  s  p ij except for its hydrostatic pressure (bulk) component. Further development will show that  s  p ij plays an essential role in linking the modeling of the porous flow-compaction and the fully solid behaviour using micromechanics. 96  Chapter 4: Integration of Modeling from Fluid to Cured Resin  4.3.3. Elastic moduli for isotropic materials Let‟s assume that the fibre-bed (solid skeleton would make more sense here!) is isotropic for the time being. A reasonable representation for the equivalent elastic moduli relevant to  s  p ij may be assumed in the form of K s  p  K fb  (4.62)  Gs  p  Gc  G fb  (4.63)  for the bulk modulus, and  for the shear modulus. The index fb represents the fibre-bed or the solid skeleton of the system, and Gc is the shear modulus of the composite obtained from a micromechanics scheme. To further elaborate, the new bulk modulus is equal to the bulk modulus of the fibre-bed since there is no contribution from the resin pressure in the relevant stress tensor,  s  p ij . However, in shear we have a combination of two load sharing mechanisms; one is the interaction of resin with individual fibres, and the other could be considered as a load sharing system between the previous mechanism and the fibre-bed. Taking a similar approach to (4.13), we write the constitutive equation in vector form   s  p   s  p G   s  p K  Gs  p  G  K s  p  K  (4.64)  The fibre-bed is assumed to be elastic. The constitutive relations for the two stress components may be written as   s  p G  Gc  G fb  G  (4.65)   s  p K  K fb  K  (4.66)  and  Summation of the above two equations leads to   s  p  Gc  G fb G  K fb  K  (4.67) 97  Chapter 4: Integration of Modeling from Fluid to Cured Resin  which could also be written in the form   s  p  D s  p   (4.68)  where  D s p  4 2 2   K fb  3 Gc  G fb  K fb  3 Gc  G fb  K fb  3 Gc  G fb   4 2  K fb  Gc  G fb  K fb  Gc  G fb  3 3  4  K fb  Gc  G fb   3  Gc  Symm .     0  0  0  0  0  0   G fb  0 Gc  G fb      0   0   0  0  Gc  G fb  0  (4.69)  When the resin is a viscous fluid any elastic contributions from the resin to the shear modulus of the composite is negligible. As a result Gc tends to be much smaller than Gfb, and therefore Gs-p obtained from (4.63) presents a very good approximation of Gfb (Gs-p ≈ Gfb). This effectively means that for pre-gelation composite, the mechanical behaviour of the „s-p‟ system is a very good approximation of the behaviour of the actual fibre-bed and therefore the flow-compaction behaviour can be captured by this formulation. At the other extreme when the composite is cured, Gc is much larger than Gfb which, once considered in (4.63), leads to Gs-p≈ Gc. The question now is whether the total stress of the system obtained by such a two-phase model is equivalent to the stress response predicted by a regular stress model. This will be referred to as the consistency of the approach and will be addressed in the next section.  4.3.4. Consistency of constitutive equations for isotropic materials Disregarding free strains: Rewriting the decomposition of the stresses of the composite system in (4.59) leads to σ s ij   s  p ij  bP ij  (4.70)  98  Chapter 4: Integration of Modeling from Fluid to Cured Resin  Let us investigate whether the following constitutive equation, obtained by expanding (4.70) for 2D plane strain, is equal to the constitutive equation of a solid composite once the mass conservation equation is also applied.   1   K fb  4Gc 3 K fb  2Gc 3 0    1  1         2    K fb  2Gc 3 K fb  4Gc 3 0    2   1bP    0 0 Gc   12  0  12    (4.71)  Note that we are considering elastic constitutive equation as the PVE response is elastic during every time-step. Stress and strain components and resin pressure in the above equation are in fact representative of the change of these parameters during an arbitrary time-step. Therefore, during this arbitrary time-step equation (4.50) may be written without the time differentiation as  v s mech   tr ζ s 2K c  Gc 3  (4.72)  which upon considering the components of ζs in (4.71) results in   v s mech   1   2  2K c  Gc 3  (4.73)  In the absence of flow and free strains, the mass conservation equation then becomes  1   2   1   2  2K c  Gc 3  (4.74)  Combining the mass conservation in the above equation with (4.71) we have 2K c  Gc 3 1   2   2K fb  Gc 3 1   2   2bP  (4.75)  which after simplification leads to bP  K fb  K c  1   2   (4.76)  99  Chapter 4: Integration of Modeling from Fluid to Cured Resin  If the above equation for the pressure is substituted into (4.71), we will arrive at   1   K c  4Gc 3 K c  2Gc 3 0    1        2    K c  2Gc 3 K c  4Gc 3 0    2     0 0 Gc   12   12    (4.77)  which is essentially the constitutive equation of elastic solid composites in 2D plane strain. This proves the consistency of our approach with the elastic behaviour of solid isotropic composites. Including free strains: In the presence of free strains, the above derivations are slightly different. We assume that the porous skeleton does not undergo any thermochemical strains and therefore equation (4.71) remains the same. This assumption makes sense as the majority of the thermochemical effects occur in the resin, and using the total strains in (4.71) facilitates the flow of resin through the porous skeleton due to its volume change. Equations (4.72) and (4.73) remain the same, but the mass conservation equation is written as  1   2  2 f   1   2  2Kc  Gc 3  (4.78)  where  f is the free strain in either 1 or 2 directions. Combining the above with (4.71) we eventually arrive at bP  K fb  Kc 1   2   2Kc  Gc 3 f  (4.79)  Again, if the above equation is substituted into (4.71) we will have f  1   K c  4Gc 3 K c  2Gc 3 0  1        f   2    K c  2Gc 3 K c  4Gc 3 0   2       0 0 Gc    12   12    (4.80)  100  Chapter 4: Integration of Modeling from Fluid to Cured Resin  which is the constitutive equation of elastic solid composites in 2D plane strain including free strains. This proves the consistency of our approach with the elastic behaviour of solid isotropic composites in the presence of free strains.  4.3.5. Elastic moduli for transversely isotropic materials For anisotropic materials, we need to present the components of the composite stress tensor in a slightly different format from (4.59) to have  s ij   s  p ij  bI P ij  (4.81)  where bI represent the components of the vector of Biot coefficients, b, associated with the anisotropic porous structure. The index I takes the same value as i, and is presented in upper-case form to indicate that the summation rule in index algebra is not to be applied. In the case of transversely isotropic materials, b can be presented as follows [Coussy (2004)]  b1    b  b  b     (4.82)  where b is the Biot coefficient in the plane of isotropy, and b1 is the Biot coefficient in the direction of the fibres. In this work, we assume that b may be defined by b  1  K 2 fb K2 f  (4.83)  Due to the very high stiffness of the fibre-bed in the longitudinal direction, the value of b1 has no discernable effect on the phenomena related to the porous flow of the resin through the system and the resulting compaction. However, as we will see later, the value of b1 plays an important role in making the approach consistent with the regular constitutive model for cured composite materials. If the tensors of stress and strain are represented as vectors in Voigt notation, (4.81) may be expressed as  101  Chapter 4: Integration of Modeling from Fluid to Cured Resin   11   11  b1      b   22   22     33   33   b         P  23   23  0  13   13  0        12  s  12  s  p  0   (4.84)  for transversely isotropic materials. With a slight change from the 3rd equation in (4.60), the total equilibrium equation of the system is now written as     s  p ij , j   bI P,i  0  (4.85)  Taking a similar approach to the isotropic case, we decompose the „s-p‟ stress to have ζ s  p  ζ s  p  ζ s  p   ζ s  p K  ζ s  p G  ζ s  p G 2  12  (4.86)  23  In vector form, we may write the foregoing equation as   s p   s p   s p    s p K   s p G   s p G 2  12  23    s p     s p    K 2 s p  K2  G12s p  G12  G23s p  G23  (4.87)  A reasonable representation for the elastic moduli related to „s-p‟ stress may be assumed in the form K 2 s p  K 2 fb  (4.88)   s p   fb  (4.89)  for the moduli dealing with the bulk behaviour of the material, and   s  p   c    (4.90)  as the bulk response has a small role in the value of η evident from (4.24). Finally, the shear moduli are introduced similar to the formulation of shear modulus in the isotropic system. G12s p  G12c  G12 fb  (4.91)  102  Chapter 4: Integration of Modeling from Fluid to Cured Resin  G23s p  G23c  G23 fb  (4.92)  As it is seen here, the concept is the same as in the isotropic case. Assuming pseudoviscoelastic behaviour for the resin, and therefore the composite, and that the fibre-bed is elastic, we write  s  p K  K 2 fb  K 2  (4.93)  2   s  p   c     (4.94)   s  p    fb    (4.95)   s  p G  G12c  G12 fb G  (4.96)  12  12   s  p G  G23c  G23 fb  G 23  (4.97)  23  Summation of the five components of the rate of „s-p‟ stress vector leads to  s  p  D s  p   (4.98)  where  D s p   c               fb  K 2 fb  G 23c  G 23 fb     fb  K 2 fb Symm.    G  K 2 fb  G 23c  G 23 fb 23 c   G 23 fb     0  0  0  0  0  0  G 23c  G 23 fb  0    0   0  0   0   G12 fb  0  G12 c  G12 fb G12 c  (4.99) Similar to the arguments made in the case of isotropic materials, when the resin is a viscous fluid the shear moduli of the „s-p‟ system obtained from (4.91) and (4.92) are good approximations of the corresponding shear moduli of the fibre-bed. Therefore, the 2D plane strain version of (4.99) for a pre-gelation unidirectional [0°] composite system may be written as  103  Chapter 4: Integration of Modeling from Fluid to Cured Resin  D s p   c       fb  0    fb K 2 fb  G23 fb 0  0   0  G12 fb   (4.100)  The only difference between the material stiffness matrix of the „s-p‟ system in the above equation and its counterpart for the fibre-bed is in the 11 component where ηfb is replaced by ηc+Δη. However, both of these values are typically in the order of 100 GPa, and much larger than the other component of the material stiffness matrix. Therefore, the above difference will not have any discernable effect on the compaction response of the system. Similar to the isotropic case, this means that for pre-gelation composite the mechanical behaviour of the „s-p‟ system is a very good approximation of the behaviour of the actual fibre-bed and hence the flow-compaction behaviour is correctly represented. At the other extreme of cured composite, using similar arguments made for isotropic materials, the shear moduli of the „s-p‟ system are very good approximations of the corresponding shear moduli of the composite obtained by micromechanics. Therefore, the 2D plane strain version of (4.99) for a cured unidirectional [0°] composite system may be written as  D s p   c       fb  0   fb K 2 fb  G23c 0  0   0  G12 c   (4.101)  The consistency of the approach using (4.101) for predicting the stress response of transversely isotropic materials will be addressed in the next section.  4.3.6. Consistency of constitutive equations for transversely isotropic materials Disregarding free strains: Similar to the isotropic case, let us focus on an arbitrary timestep and investigate the consistency of the approach for the elastic behaviour during that time-step. The following is the constitutive equation of a unidirectional solid composite with [0°] fibres in 2D plane strain.  104  Chapter 4: Integration of Modeling from Fluid to Cured Resin   1  c     2    c    0  12    c K 2 c  G23c 0  0   1    0   2  G12c   12   (4.102)  Let us investigate whether the following constitutive equation obtained by expanding (4.81) could be equal to the above relationship once the mass conservation equation is also applied.   1   c       2     fb    0  12     fb K 2 fb  G23c 0  0    1  b1      0    2    b P G12c   12   0   (4.103)  The change in the longitudinal modulus  is considered in a general form as  is defined by (4.24) where the second term is a function of both K2 and  . As both of these values are changed, it is not easily intuitive (as we will see later) to obtain the consistent change in the longitudinal modulus. Therefore, we have opted to apply an arbitrary change to the value of longitudinal modulus. We set the stresses from the two matrix equations in (4.102) and (4.103) to be equal and obtain   1   c    1   fb  2  b1 P   c  1   c  2    2   fb  1  K 2 fb  G23c  2  bP   c  1  K 2 c  G23c  2      (4.104)  The above may be rewritten as b1 P   1   fb   c  2   bP   fb   c  1  K 2 fb  K 2 c  2      (4.105)  As the above should be valid for any combination of strain components, the following must be always correct   fb   c  b1    b  fb   c  K 2 fb  K 2 c      (4.106)  105  Chapter 4: Integration of Modeling from Fluid to Cured Resin  that leads to the values for b1 and  that are consistent with the two constitutive equations being equivalent, i.e. b1    fb   c K 2 fb  K 2 c       (4.107)  b   c   2  fb  (4.108)  K 2 fb  K 2 c  Substituting the above results into (4.103), the constitutive equation may be written as    2  1  c   fb   c  K 2 fb  K 2 c     fb  2       0  12       fb K 2 fb  G23c 0    0   1   fb   c  K 2 fb  K 2 c    0   2    1     G12c   12   0     bP    (4.109) If we can demonstrate that applying the mass conservation equation to the above will lead to the same definition for bP that was obtained in (4.105), we have proved that the two approaches are equivalent. The normal components of strain are calculated as 1  c    2   c  1    1   K 2 c  G23c 1    2  K 2 c  G23c   2  c K 2 c  G23c    c    c c    c   1    c   2   (4.110)  Adding the two normal components of strain, the volumetric strain of the system is obtained as   v  1   2   1 K 2c  G23c   c 1  c   c  2  c K 2 c  G23c    c 2  (4.111)  We substitute the stress components in the above equation with the relationships obtained by our approach in (4.109) to arrive at  106  Chapter 4: Integration of Modeling from Fluid to Cured Resin  1   2   1 c K 2 c  G23c    c 2       c 2       fb   c bP  K 2 c  G23c   c c  fb 1 fb 2 K 2 fb  K 2 c  K 2 fb  K 2 c                 K  G   bP   c c fb 1 2 fb 23c 2          (4.112) which after extensive manipulations leads to      bP   fb   c 1  K 2 fb  K 2 c  2  (4.113)  therefore proving that the two constitutive equations in (4.102) and (4.103) are in fact equivalent. Including free strains: Assuming the presence of free strains, the above proof is still valid providing that we substitute  1 and  2 with  1   1f and  2   2f respectively. Using an identical approach to the above, it is readily proved that applying the mass conservation equation in its simplified form in the absence of resin flow as below   v   vf   vmech  (4.114)  to the constitutive equation of two-phase systems in the form of   1   c       2     fb    0  12     fb K 2 fb  G23c 0  0   1   1f  b1      0   2   2f    b P G12c    12   0   (4.115)  leads to an approach equivalent to using the common constitutive equation of transversely isotropic materials which follows   1   c     2    c    0  12    c K 2 c  G23c 0  0   1   1f    0   2   2f  G12 c    12   (4.116)  107  Chapter 4: Integration of Modeling from Fluid to Cured Resin  This ensures that in the absence of flow, the predictions of the current approach for the development of residual stresses of the composite system are equivalent to the response observed using the regular approach. However, in modeling resin flow, using (4.115) practically means that we have assumed that the free strain of the system (predominantly due to thermal expansion of resin in the stages that flow occurs) occurs in the fibre-bed. Therefore using this approach, the thermal free strains will not provide a motivation for flow of resin as the effective fibre-bed does not undergo any stress due to the system thermal expansion and hence does not put the resin under pressure to flow. This is clearly in contrast with the micromechanical behaviour of the system at the early stages of processing when the resin can move relative to the fibre-bed and the thermal strains of the fibre-bed are practically negligible compared to those of the resin. In such conditions, it is expected that the thermal expansion of the resin will cause the fibre-bed to expand as well initially. This leads to tensile stresses in the fibre-bed, which in turn puts the resin under pressure that leads to flow of resin out of the system until the pressure is vanished and the whole system reaches equilibrium. One should base the two-phase approach on (4.103) to capture the flow of resin correctly since the stresses of the effective fibre-bed are independent of the free strains of the system. To correctly predict the development of stresses due to cure, applying the mass conservation equation to (4.103) must lead to (4.116). If we set the stresses from the two matrix equations in (4.103) and (4.116) to be equal, we have            1   c    1   fb  2  b1 P   c  1   1f   c  2   2f   f f   2   fb  1  K 2 fb  G23c  2  bP   c  1   1  K 2 c  G23c   2   2              (4.117)  that may be rewritten as  b1 P   1   fb   c  2   c  1f   c  2f   f f  bP   fb   c  1  K 2 fb  K 2 c  2   c  1  K 2 c  G23c  2      (4.118)  As was argued before, the above should be valid for any combination of strain components and free strain components, and therefore the following must be always correct 108  Chapter 4: Integration of Modeling from Fluid to Cured Resin   fb   c   c c b1      b  fb   c  K 2 fb  K 2 c  c K 2 c  G23c       (4.119)  which leads to 4 inconsistent equations to obtain the two unknowns of b1 and  . This shows that (4.103) and (4.116) cannot be consistent with each other. The inconsistency demonstrated above for transversely isotropic composites has led us to an approach including a switch at a nominal point of gelation for the resin system. We start the analysis based on (4.103) and past the gelation point we include the free strains in the behaviour of the effective fibre-bed, and therefore use (4.115) from that point forward. This ensures consistency of the approach at both stages of uncured fluid resin and cross-linked resin with the actual behaviour of the composite system at those stages respectively.  4.3.7. Decomposition of stresses in 3D The fibres could be placed in any direction in the plane of the composite laminate. If the laminate coordinate system subtends an angle,  , with the material coordinate system in the 1-2 plane of the laminate, the transformation tensor between the two coordinate systems is then expressed by  Cos n   Sin  0   Sin Cos 0  0 C  S 0 0   S C 0 1  0 0 1  (4.120)  The transformation between the stress tensor in the laminate coordinates, ζˆ , and the stress tensor in the material coordinate systems takes the form ζˆ  n T ζn  (4.121)  ˆ ij  nki nlj kl  (4.122)  or in indicial notations  Using the Voigt notation, the stress vector in the laminate coordinates 1  2  3 may be written as 109  Chapter 4: Integration of Modeling from Fluid to Cured Resin  ˆ 11   11  ˆ     22   22   33  ˆ 33  1     T   ˆ 23   23  ˆ 13   13      ˆ 12   12   (4.123)  where the transformation matrix is defined by  T  1   C2  2  S  0   0  0   CS  S2  0  0  0  C2  0  0  0  0  1  0  0  0  0 C  S  0  0  S  C  CS  0  0  0  2CS    2CS  0   0  0   C 2  S 2   (4.124)  To transform the stress vector back to the material coordinate system we need to have C 2  2 S 0 T   0 0  CS  S2  0  0  2  0  0  0  1  0  0  0  C  0  0 S  C   CS  0  0   2CS   0 2CS  0 0   S 0  C 0   0 C 2  S 2  0  (4.125)  The stress vector in the laminate coordinates 1  2  3 may be written as ˆ11   11  ˆ     22   22  ˆ33   33  1    T   ˆ 23   23  ˆ13   13       12  ˆ12   (4.126)  where  110  Chapter 4: Integration of Modeling from Fluid to Cured Resin   C2  2 S  0 T    0  0  2CS  S2  0  2  0  C  0  1  0  0  0  0   2CS  0   CS   0 0 CS  0 0 0   C S 0  S C 0   0 0 C 2  S 2  0  0  (4.127)  The matrix of elastic moduli in the laminate coordinate system may be written in the form 1 Dˆ  T  DT   (4.128)  The decomposition of stress to „s-p‟ stress and resin pressure is applied in the material coordinate system. Therefore, to obtain the total stress of the composite, one needs to transfer the „s-p‟ stress into the material coordinate system to be added to resin pressure components, and then transfer the result back to the laminate coordinate system. As a result, (4.84) may be written in the form   C 2 b1  S 2 b  b1      2 b  2    S b1  C b        b    b     s  T  1  T   s  p    P    s  p   P  0 0     0    0           0   CS  b1  b     (4.129)  4.4. Modifications in the Numerical Technique The additional implementations in order to capture the integration of flow and stress in the FE formulation are discussed in this section. The FE equations for the isotropic and transversely isotropic conditions are presented. Also, the required modifications in the non-linear solution scheme and the updating procedure of the resin volume fraction are discussed.  111  Chapter 4: Integration of Modeling from Fluid to Cured Resin  4.4.1. FE implementation in 2D plane strain Isotropic materials: Assuming 2D plane strain conditions, the „s-p‟ matrix of elastic moduli is expressed as  D s p  4 2   0  K fb  3 Gc  G fb  K fb  3 Gc  G fb     4   K fb  Gc  G fb  0 3   Symm. Gc  G fb      (4.130)  The most important changes need to be applied in the FE implementation of the mass conservation equation. Integrating (4.45) over the domain Ω using some weight functions wp we arrive at  v  i ,i    wp d   ui ,i wp d   v s   mech    wp d   v s  free  wp d  0  (4.131)    where v s mech can be obtained from (4.50). Substituting (4.59) into (4.50) leads to  1   tr ζ s  p   bP  K c  Gc 3 t  2   v s mech   (4.132)  Incorporating (4.132) in (4.131) and discretization of the result for a 9-3 element gives the components of the last row of matrices of coefficients as K pu  Cpv  0 318  (4.133)  K pv   N p δ T Bd  (4.134)  T    Cpu   N p δT Bd T    C pp    318  318  T 1 N p δT Ds  p Bd  2Kc  Gc 3  318  T b N p N p d  K c  Gc 3  33  K pp  0 33  (4.135)  (4.136) (4.137)  112  Chapter 4: Integration of Modeling from Fluid to Cured Resin  There is a new force term defined as a function of the rate of free strains is introduced in the FE form of the mass conservation equation  f p    N p δT D c D s  p ε free d   δ T ε free N p d 1  T  T      (4.138) 31  Also in the FE formulation of the third equation in (4.60), an additional term is added to the force term on the right-hand side of the matrix equations  f u   N Td t ( n ) d   B T ζ free d     (4.139)  181  Combining the above with the implementation of the other two governing equations leads to the matrix equations of an element in the form of  K uu   0  0   K uv K vv K pv  K up  u  Cuu 0 0   u   fu         K vp  v   C vu 0 0   v    f v  t 0  p  Cpu 0 Cpp   p   fp  391  (4.140)  Transversely isotropic materials: In 2D, the global x-y plane is equivalent to 1  3 plane in the 3D laminate coordinate system. Thus, we may obtain the full 3D matrix of elastic moduli from (4.128), and extract the relevant components to have (note that we have dropped the “hat” symbol for simplicity in this section.)  C 4  2C 2 S 2  S 4 K 2  G23   4C 2 S 2G12 C 2  S 2 K 2  G23   0   D K 2  G23 0  2 2  Symm. S G23  C G12   (4.141) In special cases, we may be able to simplify the above equation, e.g.    0    0     D K 2  G23 0   Symm. G12   (4.142)  113  Chapter 4: Integration of Modeling from Fluid to Cured Resin    90    K 2  G23 D    Symm.  K 2  G23  0  0  G23   K 2  G23  (4.143)  In the x-y plane, (4.129) is simplified to C 2 b1  S 2 b  x   x        b  y    y    P       0  xy  s  xy  s  p    As the plane strain conditions are assumed (i.e.  z  (4.144)   0;  z  mech  free   0 ), the volumetric  strain of the composite and its rate are respectively expressed by  v s mech   x mech   y mech  (4.145)  v s mech  x mech  y mech  (4.146)  Using the expanded form of the constitutive equation of the composite material, we may write 1  s x  Dc121  s y  Dc121  s x  Dc 221  s y v s mech  Dc11  (4.147)  We expand the first two rows of (4.144) to have   s x   s  p x  C 2 b1  S 2 bP  s  y   s  p y  bP  (4.148)  and then the above is substituted into (4.147) to result in         bD    1 1 1 1 v s mech  Dc11  Dc12  s  p x  Dc12  Dc 22  s  p y      1  1  1 c12   C 2b1  S 2b Dc11  Dc12    1  Dc 22 P  (4.149)  where    σ s  p ij  Ds  p ijkl  klmech  Ds  p ijkl  kl   klfree    (4.150)  114  Chapter 4: Integration of Modeling from Fluid to Cured Resin  Substituting the above into (4.131) and discretization of the result for a 9-3 element gives the modified form of some of the components in of the FE matrices compared to the isotropic case as      Cpu   N p δ T B  D c D s  p B d T  1           (4.151) 318    Cpp  C 2 b1  S 2 b Dc 11  Dc 12  b Dc 12  Dc 22   N p N p d 1  1  1  1  T    K up  C 2 b1  S 2 b   p   BT  b N d    0    (4.152) 33  (4.153) 318  to be utilized in the matrix equations of the finite element presented in (4.140).  4.4.2. Modifications in the non-linear solution scheme With the new matrix equations in (4.140), we need to revisit the nonlinear solution scheme. The iterative matrix equation for the solution of the system is therefore slightly modified from (2.146) to be presented as  C uu  tK uu  C vu   C pu   tK uv tK vv tK pv   f u n1  f u int n tK up   u nk11  u nk1        k 1 k tK vp   v n 1    f v n 1  C vu u n  u n 1      k  1 k  C pp   p n 1   f p n1  C pu u n  u n 1  C pp p n  n 1 k          (4.154)  4.4.3. Calculation of volume fractions With the assumption of compressible phases and also considering the free strains, it is essential to revisit our method for obtaining the values of fluid volume fraction over the domain. The formula used to update the fluid volume fraction in (2.119) and (2.120) is  115  Chapter 4: Integration of Modeling from Fluid to Cured Resin  based on the assumption that the phases are incompressible and therefore the volumetric strain of the system is solely dependent on the volume of fluid flowing into or out of the domain. The mass conservation equation in the 1st equation of (2.62) may be rewritten in the form ui ,i  vi ,i or v  v flow  (4.155)  Integrating the above equation through time leads to t   v   v flow    vi ,i dt  (4.156)  0  As it is evident from (4.156), if the phases are assumed to be incompressible, the volumetric strain of the system is equal to the component of volumetric strain that happens as a result of fluid flow. Therefore, using the total volumetric strain in (2.119) and (2.120) is correct. Here, the mass conservation equation has the general form as in (4.45), and as a result, the flow-dependent component of volumetric strain may be obtained from t  t  0  0       v flow    vi ,i dt   u i ,i  v s mech  v s free dt  (4.157)  and substituted into the following equation (similar to (2.119)) for updating the fluid volume fraction    0   v flow 1   v flow  (4.158)  In the models dealing with cured material and stress development, the change of volume in the phases is assumed to have no effect on the volume fractions of the composite material. For our approach to be consistent with that, we need to obtain the component of volumetric strain of the system that is only due to flow of the fluid phase, and incorporate that into updating the volume fractions. That is the reason why in (4.158), we have opted for an equation that assumes that the phases are incompressible. If one decides to consider the effect of volume change of the phases in updating the volume fraction 116  Chapter 4: Integration of Modeling from Fluid to Cured Resin  values, even when no flow occurs, a non-zero change in the volume fraction of the phases will be predicted which is not consistent with the models only dealing with stress development in curing composites.  4.5. Resin Properties during Cure Even though quite a lot of work has been done in the literature to characterize the properties of resin during cure, a material model with consistent viscosity and viscoelastic properties of resin was nonexistent at the start of this work. Typically, the thermomechanical characterization of resin is focused either on finding the relationship between degree of cure and viscosity or developing formulas for resin properties as a viscoelastic or cure-hardening elastic solid. Therefore in this work we have opted to use separate set of data for the viscosity and the moduli of resin. Ideally, comprehensive rheological and mechanical experiments should be performed on the resin at various states of temperature and cure to capture the full range of behaviour of resin from a viscous fluid to a viscoelastic solid. Having full-range data, ensures that the evolution of viscosity during cure is consistent with how the viscoelastic properties of the resin change throughout processing. Very close to the completion of this thesis, Thorpe (2012) presented a consistent material model that predicts the viscoelastic liquid and viscoelastic solid behaviour of a commercially available thermoset resin, MTM45-1 epoxy, for the full range of degree of cure. This development is very promising as this material characterization is designed to be directly fed into the current integrated approach for resin flow and stress modeling.  117  Chapter 4: Integration of Modeling from Fluid to Cured Resin  4.6. Tables  Table 4.1: Relaxation times and weight factors for 3501-6 resin for α0 = 0.98  ω  τω (min)  Wω  1  2.92E+01  0.059  2  2.92E+03  0.066  3  1.82E+05  0.083  4  1.10E+07  0.112  5  2.83E+08  0.154  6  7.94E+09  0.262  7  1.95E+11  0.184  8  3.32E+12  0.049  9  4.92E+14  0.025  118  Chapter 5: Numerical Applications  To verify the current approach toward integrating the numerical modeling throughout the processing of composite materials a few examples are presented. These examples will verify some or all the modeling aspects (such as compaction behaviour, resin flow, residual stress development) that are of consequence in processing problems.  5.1. Cured Materials In this section, we focus our attention towards analysing fully cured, linear elastic, isotropic composites using the approach developed in this thesis and comparing the results with the classical approach for elastic solids. As benchmark solution, we will decouple the solution of the problem in the MATLAB code so that it will just solve the elastic solid problem using the classical approach, and then compare the results with the response obtained by running the code in the coupled multi-physics mode for the same example. We will also use the commercial code ABAQUS to verify the results, and compare the convergence rate of the two approaches. Other than using 9-3 elements, we will also apply the above approach to regular bi-quadratic 9-4 element based on the u-p formulation and compare their response with the 9-3 element. All the examples are assumed to be under 2D plane strain conditions. We will start with a very simple example, whose exact solution can be obtained theoretically, and move along to further complicated problems. We assume the following relationship between the fibre-bed and composite bulk moduli K fb  Coeff K  Kc  (5.1)  where CoeffK provides a simple representation of the value of fibre-bed bulk modulus.  119  Chapter 5: Numerical Applications  5.1.1. Cured composite sample under uni-axial strain conditions A rectangular sample of an isotropic two-phase composite material, as shown in Figure 5.1, is constrained on 3 sides, and is put under a constant distribution of normal pressure load. The input properties for this problem are   f  0.6,   m  0.4  (or   0.4) E f  10GPa,  E m  1GPa   f  0.3, h  4mm,   m  0.4 L  12mm  f 0  700kPa For simplicity, we assume that the phases act in series together in both shear and volumetric responses. Thus, the composite shear and bulk moduli may be obtained by Gc   Kc   Gm G f G m f  G f  m   0.7837GPa  Km K f K m f  K f  m   3.205GPa  The matrix of elastic moduli is then expressed as 4   K c  3 Gc  D   Symm.   2 K c  Gc 3 4 K c  Gc 3   0 4.2501 2.6827 0     0 4.2501 0 GPa  0.7837 Gc   Symm.   There is a state of uniform stress in the sample under uni-axial strain conditions, so we just apply the boundary conditions to the constitutive equation to have 0  x   4.2501 2.6827 GPa        700kPa 2.6827 4.2501  y   and solve for the unknowns in stress and strain fields 120  Chapter 5: Numerical Applications   y  1.647 10 4    x  441.84kPa  The strain energy of the system may be calculated by      1 1 E s   y  yVol .   0.44184MPa   1.647  10 4 4mm12mm1m 2 2 Es  0.002767 N .m  This example was run using the developed 9-3 elements for flow through deforming porous media. Due to the simplicity of the stress and strain fields in this specific problem, the exact results are obtained using any number of elements. Different meshes were chosen including: 1 element, 3×2 elements, 6×4 elements, and 9×6 elements all giving the exact results in terms of stress and strain components and also the calculated value for the strain energy. It is also very important to note that the elastic response of the composite system using this approach is independent of the value of the bulk modulus of the fibre-bed Kfb. This point numerically confirms the capability of the integrated approach to model the stress response of the post-gelation composite, which was proved theoretically in the section on the consistency for isotropic materials in Chapter 4.  5.1.2. Cured composite sample under linear distribution of pressure A sample with the same geometry and constitutive properties as the previous example is under a linear distribution of pressure load as shown in Figure 5.2. The new loading condition is defined by  f1  400kPa,  f 2  1000kPa.  The linear distribution of the load introduced significantly more complexity into the stress and strain fields of the sample compared to the previous problem. Shear deformation occurs in the whole domain, and the stress and strain fields are functions of the location of the point on the domain. Figure 5.3 shows the profile of the vertical displacement at the top surface of the sample for different meshes of 9-3 elements. In Table 5.1, the value of strain energy for the system is reported for the combination of different values of CoeffK, mesh densities, and the conditions on relative velocity of 121  Chapter 5: Numerical Applications  resin. The convergence of the results with refining the mesh and decreasing the value of CoeffK is very fast in this case. It is noteworthy to remind the reader that assigning unity to CoeffK in 9-3 element leads to the regular approach of solid mechanics. Since the total bulk modulus of the two-phase medium is considered in the total equilibrium equation, there is no motivation for flow of resin. In other words, setting CoeffK=1 leads to a very stiff fibre-bed (with a precise stiffness) that would carry all the hydrostatic pressure in the system, leading to the pressure of the resin to take zero values on the domain and thus eliminating the cause of flow. Here, since we aim to solve the solid composite problem using the formulation of flow in deforming porous media, the relative velocity of resin should be zero. For efficiency of the solution and also to converge to the exact response of the solid composite, one may constrain all the DOFs corresponding to the resin relative velocity. However, this would mean that in process modeling we would have to change the essential boundary conditions midway of the analysis at gelation point. We expect that introduction of a high value for the “viscosity” of solid resin leads to practically zero values for relative velocity of resin, and therefore omitting any specific need for constraining the velocities at gelation. Table 5.1 confirms the above deduction as the values of strain energy are very close in the two cases of constrained and free relative velocities. It should be noted that even with high values of viscosity, given enough time resin will eventually flow out of the system if the relative velocity DOFs are not constrained. This would be considered undesirable when we intend to model a solid composite problem. In the numerical examples in this work with free conditions on relative resin velocity, we have assumed the following for viscosity of resin and the end-time of analysis    1 10 7 Pa.s,  t  18s  An interesting observation is that refining the mesh not only does not reduce the very small difference between free and constrained conditions, but also increases it. One may explain this using the argument that by refining the mesh the number of velocity DOFs increases, leading to a softer response on part of the relative flow of resin. 122  Chapter 5: Numerical Applications  5.1.3. Cured Composite sample under pressure gradient A rectangular sample of an isotropic two-phase composite material, as shown in Figure 5.4, is constrained at top and bottom surfaces, and is put under a pressure gradient in the horizontal direction. The properties of the problem are the same as previous examples except for  f1  700kPa,  f 2  100kPa.  Using axial symmetry half of the problem will be modeled using meshes of 9-3 and 9-4 elements. To have a schematic shape of the deformation of the sample under these loading conditions, Figure 5.5 shows the displacement profiles of the two vertical sides of the specimen. The combination of the large difference in the pressure values on the two sides and the geometry of the problem force all points of the domain to move to the right. The displacement values on the right side are much smaller than the left. These smaller displacements resulted in a slower convergence of the numerical results on the right side. Consequently, we opt to focus our attention on the displacement profile and the maximum displacement of the right side for comparison and discussion. Table 5.2 reports the obtained values of the maximum horizontal displacement (uA at point A shown in Figure 5.4) using ABAQUS and also the 9-3 element under various conditions. In terms of comparing the free and constrained conditions on relative resin velocity, the observations made in the previous example are applicable again. In short, the results of the free and constrained conditions are practically the same given that the viscosity of the resin is assumed high enough for the time-scale of modeling. Assuming values for the viscosity of the solid resin that are too high could possibly introduce numerical instabilities due to ill-conditioning of the matrices of coefficients. The ABAQUS results are obtained using the 2D plane strain 8-noded element. This element is specialized for solid problems, and takes the regular approach of discretizing the total equilibrium equation. As it is seen in Table 5.2, the convergence of the displacement values with refining the mesh happens quite slowly (using a uniform mesh). 123  Chapter 5: Numerical Applications  As was mentioned earlier, with CoeffK=1 the results of the porous media approach to the solid problem at hand, are comparable to the regular approach (ABAQUS results for instance), since in this case only the total equilibrium equation (among the three governing equations) is effectively considered. In the case of the 6×4 mesh, if we compare the displacement values obtained from ABAQUS 8-noded element and 9-3 element with CoeffK=1, while both values are in the same range, the result from the latter is closer to the converged value (~1.649 using ABAQUS). This can be attributed to the advantage of 9-noded over 8-noded element. Due to limitations of the developed MATLAB code the finest mesh that we have modeled is 12×8. Figure 5.6 to Figure 5.8 show the convergence of the displacement profiles with decreasing values of CoeffK using various meshes. Figure 5.9 shows the convergence of the displacement profiles with refining the FE mesh for the case of CoeffK=0.01. Figure 5.10 represents the same information but for CoeffK=1. Comparing the two graphs, it is evident that the convergence of the displacement response as a whole is very quick for CoeffK=0.01 compared to CoeffK=1 (which effectively represents the regular solid mechanics approach as was mentioned before). Reducing CoeffK from unity introduces tangible improvement in the rate of convergence of the displacement results by refining the mesh. For the composite materials that we deal with, CoeffK~0.0001-0.01 could be a realistic range. Table 5.3 offers similar comparisons to the ones made in Table 5.2, but in terms of total strain energy of the system, where similar observations can be made. Table 5.4 and Table 5.5 report the obtained values of the maximum horizontal displacement and total strain energy respectively, using the 9-4 element under various conditions. Figure 5.11 and Figure 5.12 compare the displacement profiles for 9-3 and 94 elements for the case of CoeffK=0.01. It is evident that the results obtained using 9-4 elements assuming permeable B.C. on the two sides, are not satisfactory. Yet, they show improvement with mesh refinement. In the impermeable case, the displacement and energy results are satisfactory but slightly inferior to 9-3 results in terms of rate of convergence both with refining the mesh and reducing CoeffK. The common property of  124  Chapter 5: Numerical Applications  the 9-3 and 9-4 elements is the 9-noded representation for the solid structure. As a result, we should expect the two elements to produce identical results in terms of displacement profiles and strain energy values for the case of CoeffK=1, which is exactly the case here (assuming impermeable B.C. for 9-4 elements). Discussion: So far, not only we are able to model the cured (isotropic, for the time being) composite using the framework that we use for pre-gelation composite with fluid resin, but also we achieve improvements in terms of the convergence rate at least for modeling challenging examples, such as the last example above. As for the reason for the above improvement, It may be attributed to the fact that using the integrated approach we introduce an extra condition in the form of mass conservation equation. This extra equation, enforces a kinematic constraint on the volumetric behaviour of the composite material, which tends to be dominant over the equilibrium equation if activated through using CoeffK<1. Having this kinematic constraint explicitly in one of the governing equations seems to be a stronger form than its regular implicit introduction as a part of constitutive behaviour in total equilibrium equation. This may serve as a qualitative explanation for this behaviour at the present time. Regarding the 9-4 element and the significant difference between the predictions of the results from the case with permeable B.C. and impermeable B.C., the explanation is that in the permeable case, we need to force the pressure values at the two boundaries to predefined values of pressure loading, and this causes the irregularities that we observe in the last example. To examine the above theory about the effect of the constraining conditions of pressure degrees of freedom on the general response of the problem, we look at the distribution of pressure on the surface of the example in Figure 5.4 for various cases. Figure 5.13 presents the pressure distribution using a 6×4 mesh of 9-3 elements. Figure 5.14 refers to the same data and the same mesh but using 9-4 elements assuming permeable B.C. on the sides. It is obvious from the comparison of the two surfaces that the permeability assumption has deteriorated the pressure surface by forcing the values to be equal to the applied pressure on the boundaries. Figure 5.15 shows the pressure distribution for the 125  Chapter 5: Numerical Applications  same mesh of 9-4 elements but with impermeable B.C. on the sides. Without the constraints (i.e. permeable B.C.), it is observed that the pressure surface is much closer to the one in Figure 5.13 except for the pressure value at the node right on the top-left corner of the sample. This improvement in the general distribution of pressure explains the improvement that was observed in the displacement and strain energy results previously. Figure 5.16 is obtained using a proposed 9-4 element based on the u-v-p formulation of the 9-3 element presented in this work rather than the regular u-p formulation in FE for porous media. It is clearly seen that the graph in Figure 5.16 is practically the same as the one in Figure 5.15. In the 9-3 element‟s formulation, there is no essential BC on pressure DOF‟s and that is the reason why the last graph here, is the same as the one for the case of 9-4 elements with impermeable BC. Also, the agreement between the responses for the two formulations shows that both approaches give the same response as expected, and that the negative value in pressure at the top left corner is not a problem of the formulation, but is a result of the placement of a pressure node at a singular point with a considerable shear deformation.  5.1.4. Transversely isotropic samples In this section, we use the presented approach to analyse fully cured linear elastic transversely isotropic composite samples and compare the results with the classical approach used in elastic solids. Having compared the results obtained from 9-3 elements with different constraining conditions, and also 9-4 elements having different B.C., we opt to focus our attention just on 9-3 elements with constrained relative velocity of resin. The major goal here is to examine the results of the specific approach presented for the treatment of transversely isotropic composites. We will also use the commercial code ABAQUS to verify the results, and compare the convergence rate of the two approaches. All the examples are assumed to be under 2D plane strain conditions. We will start with a very simple example, whose exact solution can be obtained theoretically, and progress to more complicated problems. Similar to the isotropic case, we assume the following relationship between the fibre-bed and composite in-plane bulk moduli  126  Chapter 5: Numerical Applications  K2 fb  Coeff K  K 2c  (5.2)  where CoeffK provides a simple representation of the value of fibre-bed in-plane bulk modulus. Constrained sample under uni-axial strain condition: A rectangular sample of a transversely isotropic two-phase composite material, as shown in Figure 5.1, is constrained at 3 sides, and is put under a constant distribution of normal pressure load. The fibres are placed so that they are along x axis (   0 ). The properties of the composite material are  E1  127.4GPa   12  0.253 G12  4.068GPa K 2  6.4055GPa G23  2.7262GPa E2  7.6198GPa   23  0.398 consisting of fibres with the following properties E1  210GPa  E 2  17.2GPa    G12  27.6GPa   K 2  11.57GPa  12  0.2    23  0.25   Other specifications of the problem include   f  0.6,  m  0.4  and  h  4mm,  L  12mm  f 0  700kPa  127  Chapter 5: Numerical Applications  Substituting the properties of the composite material into (4.142) leads to the matrix of elastic moduli of the problem  0  129.04 3.2412  D 9.1317 0 GPa  Symm. 4.068 There is a uniform state of stress all over the sample under uni-axial strain conditions, so we can apply the boundary conditions to the constitutive equation to have 0  x  129.04 3.2412 GPa        700kPa 3.2412 9.1317  y   and solve for the unknowns in stress and strain fields  y  7.666  105    x  248.46kPa  This example was run using the developed 9-3 elements for flow through deforming porous media. Due to the simplicity of the stress and strain fields in this specific problem, the exact results are obtained using any number of elements. A mesh of 6×4 elements was analysed giving the exact results in terms of stress and strain components. It is also very important to note that the elastic response of the composite system using this approach is independent of the value of the in-plane bulk modulus of the fibre-bed K 2 fb . Similar to the isotropic case, this can be considered a numerical confirmation of the consistency of the integrated approach to regular stress modeling that was proved in the previous chapter. Constrained sample under pressure gradient: A rectangular sample of a transversely isotropic two-phase composite material, as shown in Figure 5.4, is constrained at top and bottom surfaces, and is put under a pressure gradient in the horizontal direction. The properties of the problem are the same as previous examples except for  f1  700kPa,  f 2  100kPa.  128  Chapter 5: Numerical Applications  Three different situations for the direction of the fibres are assumed;   0, 45, and 90 . Using axial symmetry, half of the problem will be modeled using meshes of 9-3 elements. The shape of the deformation profile of the sample under these loading conditions was presented in Figure 5.5 for the isotropic case. Table 5.6 to Table 5.8 report the obtained values of the maximum horizontal displacement (uA) using  ABAQUS and the 9-3 element under various conditions and different directions of fibres. Also, the values of the strain energy of the system are obtained and reported using 9-3 elements. The ABAQUS results are obtained using the 3D brick elements. This element is specialized for solid problems, and takes the regular approach of discretizing the total equilibrium equation. The plane strain conditions are forced on the structure using appropriate boundary conditions. As was the case for the isotropic material, the convergence of the displacement values with refining the mesh happens quite slowly (using a uniform mesh). In all cases, there is a considerable improvement in the rate of convergence of results with refining the mesh once values smaller than unity are assigned to CoeffK. This observation is in accordance with what we have already seen in the isotropic case of this problem.  5.2. Materials Undergoing Cure In this section, various examples with ascending degree of complexity will be presented to investigate the response of the integrated model in cases where resin undergoes cure and transforms from a fluid to a completely cured and solid form. The interaction of the flow behaviour of resin and the stress response of the composites will also be discussed.  5.2.1. One-element examples Example 1.  Fully constrained [0°] uni-directional laminate, assuming constant  Poisson’s ratio for resin  129  Chapter 5: Numerical Applications  In this example, the stress development response of a fully constrained composite material undergoing cure is studied. The sample is a unidirectional [0°] laminate with an assumed length of 1 mm and a height of 4 mm. The material is assumed to be AS4/35016. The temperature history and the resulting degree of cure of resin are presented in Figure 5.17. Figure 5.18 shows the changes in the viscosity of the resin during the cure cycle. The Poisson‟s ratio for the 3501-6 epoxy is assumed constant (see Table 5.9) so that the response can be compared with those obtained by Zobeiry (2006). The problem is modeled by a single 4-1 element (see Figure 2.2) and the time-step size chosen for the numerical solution is 30 seconds. All boundaries are assumed to be impermeable, and since all the kinematic DOFs are located on the boundary of the single element resin flow is completely constrained. Figure 5.19 depicts the time-history of stress in the direction of the fibres using the integrated model and compares the response with those reported by Zobeiry. Figure 5.20 presents the comparison of stress response along the z-axis (through-thickness direction) with the results reported by Zobeiry (2006). In both of the above graphs a perfect match is obtained. Example 2.  Fully constrained [30°] uni-directional laminate, assuming constant  Poisson’s ratio for resin The only difference that may be observed between this example and Example 1 is in the orientation of the fibres at [30°]. Figure 5.21 presents the time-history of σx using the integrated model and compares the response with those reported by Zobeiry (2006). Figure 5.20 presents the comparison of stress response along the z-axis with the results reported by Zobeiry. Again, in both graphs a perfect match is observed. Since the system is fully constrained, the z-axis stress response in this example is the same as the one in Example 3. Example 3.  Fully constrained [0°] uni-directional laminate  This example is basically the same as Example 1 with the exception that here, the Poisson‟s ratio of the 3501-6 is assumed to be evolving based on the pseudo-viscoelastic  130  Chapter 5: Numerical Applications  formulation that was presented in the previous chapter. We assume that both G and K evolve with changes in temperature and degree of cure throughout the processing. Figure 5.22 depicts the time-history of σx using the integrated model and compares the response with those obtained from the stress model. Figure 5.23 presents the comparison of stress response along the z-axis with the results obtained from the stress model. In both graphs a perfect match is observed. The time-history of developed resin pressure predicted by the integrated model is also presented in Figure 5.24. Example 4.  Fully constrained [30°] uni-directional laminate  Similar to the previous problem, the only difference here with Example 2 is that the Poisson‟s ratio is assumed to be evolving similar to Example 3. Figure 5.25 depicts the time-history of σx using the integrated model and compares the response with those obtained from the stress model. Figure 5.23 presents the comparison of stress response along the z-axis with the results obtained from the stress model. In both directions a perfect match is observed. Since the system is fully constrained, the z-axis stress response in this example is the same as the one in Example 3. The time-history of developed resin pressure predicted by the integrated model is also presented in Figure 5.24. Example 5.  Fully constrained [0°] uni-directional laminate, with permeable BC  In the previous single-element examples, the resin velocity was constrained to ensure the accuracy of the stress aspect of the integrated model. In this example, all the properties match those of Example 3 with the exception that the top boundary is assumed to be permeable, so the resin can flow out of the sample. Figure 5.26 depicts the changes in the vertical velocity of resin at the top surface of the sample through the process. It is evident that after the resin gels, the relative velocity of resin approaches zero. At around t=45 min, there is a kink observed in the time-history of resin velocity. Since the displacement of the composite system are constrained, any volumetric changes in the system instantaneously leads to the flow of resin out of the system. Therefore, the velocity of the resin in this case is a direct function of the rate of volumetric strain of the system. Figure 5.27 shows the changes in the rate of volumetric 131  Chapter 5: Numerical Applications  strain of the system for t=15-55 min. The pattern of changes in the rate of volumetric strain are the same as that of the resin velocity depicted in Figure 5.26. The rate of the volumetric strain is a function of the rate of cure, and the model used in this work for rate of cure (see Table 3.1) has a discontinuity at α=0.3 which causes the above anomaly. The time-history of the changes in the resin volume fraction obtained by the integrated model is presented in Figure 5.28, and compared with the initial (and constant) volume fraction used in the stress model. Figure 5.29 presents the time-history of the developed resin pressure in the system. Figure 5.30 and Figure 5.31 depict the time-history of stress response along x and z axes respectively, and compare the results with those obtained by the stress model. At the early stages of the process, the resin undergoes expansion due to temperature rise and, in the presence of a permeable boundary, resin flows out of the system to relieve some of the developed pressure. In the first few minutes however, the resin does not have the time and opportunity yet to leave the system and the stress response of the integrated model more or less follows the response obtained by the stress model. But after those initial stages are passed more resin can flow out and no stresses develop in the system while resin has not gelled yet and is able to flow. After gelation of resin, the stresses start to develop once again with a trend very similar to the one observed by the stress model (as the volume fractions in the two cases are quite close in value). The above observation explains the notable difference in the final predictions of residual stresses in the system between the integrated model and stress model. Example 6.  Uni-axially constrained [0°] laminate, with permeable BC on top  In this example, the stress development response of a uni-axially constrained composite material undergoing cure is studied. The geometry and B.C. of the problem are the same as Figure 5.1 (if y-direction is replaced by z-direction) for the case of f0=0. The sample is a unidirectional [0°] laminate with an assumed width of 1 mm and a height of 4 mm. The temperature history and the resulting degree of cure of resin are the same as the previous examples and are presented in Figure 5.17. The material is assumed to be AS4/3501-6  132  Chapter 5: Numerical Applications  with the same flow-compaction properties that are listed in Chapter 3. The only exception is that in the nonlinear stress-strain behaviour of the fibre-bed, an exponential curve is fitted to the steps above the strain value of 0.1. Also, the tensile modulus of the fibre-bed is assumed to be equal to the initial compressive modulus to avoid a sudden change of modulus. The problem is modeled by a single 4-1 element. The top surface of the sample is assumed to be permeable. Figure 5.32 depicts the changes in the vertical velocity of resin at the top surface of the sample through the process. The time-history of the changes in the resin volume fraction obtained by the integrated model is presented in Figure 5.33, and contrasted with the constant volume fraction used the stress model. Figure 5.34 presents the time-history of the developed resin pressure in the system. The changes in the transverse strain of the system are presented in Figure 5.35 and compared with the results of the stress model. Figure 5.36 depicts the time-history of stresses along x-axis, and compares it with those obtained by the stress model. The observations made in the case of the previous example are also valid here. The initial agreement between the transverse strain predicted by the stress model and the integrated model is more pronounced here and continues for around ten minutes. The reason is that the system is not constrained in the z direction and the sample can initially expand until enough pressure is developed to motivate the flow of resin out of the system. Example 7.  Uni-axially constrained [0°] laminate under pressure loading applied on  the permeable boundary This example is basically the same as Example 6 with the exception that there is a uniform pressure load applied to the top surface of the sample. The geometry and B.C. of the problem can be represented Figure 5.1 (if y-direction is replaced by z-direction) with f0=540 kPa. Figure 5.37 presents the time-history of the developed resin pressure in the system. Figure 5.38 depicts the changes in the vertical velocity of resin at the top surface of the sample through the process. The time-history of the changes in the resin volume fraction  133  Chapter 5: Numerical Applications  obtained by the integrated model is presented in Figure 5.39, and compared with the two runs of stress model with initial and final volume fractions. The changes in the transverse strain of the system are presented in Figure 5.40 and compared with the results of the stress model assuming both initial and final volume fractions. Figure 5.41 depicts the time-history of stress along the x-axis, and compares the results with those obtained by the stress model assuming both initial and final volume fractions. The applied pressure provides a significant motivation for the resin flow and leads to a maximum resin velocity one order of magnitude larger than that of the previous example. This could be a reason why the initial effects observed in the last two examples are nonexistent here. The large amount of resin flow, leads to a large difference between the initial and final values of resin volume fraction. This causes the trend of the post-gelation behaviour of the system obtained by the integrated model to be different from the response predicted by the stress model with the assumption of initial volume fractions. This point, which is clearly observed in Figure 5.41, has led us to include a run of the stress model with the final value of resin volume fraction for comparison.  5.2.2. Flat unidirectional laminate undergoing cure Here, we consider a flat laminate fully constrained on (or bounded to) a rigid tool. All material properties remain the same as the previous example. The material is assumed to be AS4/3501-6, and the geometry and BC of the problem are depicted in Figure 5.42. Due to symmetry, only half of the length of the laminate is shown. The mechanical properties of the fibres and resin are presented in Table 5.10. CTEg and CTEr are the glassy and rubbery coefficients of thermal expansion, respectively. CSC is the volumetric cure shrinkage coefficient. Figure 5.43 shows the time-history of autoclave temperature. It also presents the assumed temperature time-history of the part and the resulting degree of cure of resin throughout the processing. The predictions obtained by three approaches will be analysed and compared here; the integrated model, the stress model using initial distribution of resin volume fraction, and finally the stress model using the final distribution of resin volume fraction (obtained from the integrated model). We investigate 3 different directions for the placement of the fibres: [0°], [45°], and [90°]. The applied  134  Chapter 5: Numerical Applications  load on the top and side surfaces of the sample is 540 kPa. The problem is modeled by 12×6, 16×8, and 20×16 meshes of 4-1 elements to ensure the convergence of results. Time-step sizes of 60, 12, and 6 seconds are used and the convergence of results with refining the time-steps was confirmed. The 20×16 mesh and time-step size of 6 seconds are chosen to report the results and comparisons. To verify the stresses predicted by the stress model, the beta version of the commercially available COMPRO CCA in ABAQUS [Arafath (2012)] was used incorporating the variable time PVE formulation of Zobeiry (2006) and Zobeiry et al. (2010). A 3D equivalent of our 2D plane strain was modeled using a mesh of 20×16 elements for [0°] and [90°] lay-ups. Figure 5.44 shows the comparison of profiles of σx at section AB for [0°] and [90°] lay-ups, and a very good match between the results is observed. Figure 5.45 presents the distribution of σx at the bottom surface of the [90°] laminate obtained by the stress model and compares it with the one predicted by ABAQUS. Again, a very good match between the two predictions is observed. Note that in the latter figure, oscillation is observed in the distribution of stress in the neighbourhood of the bottom right corner. These oscillations occur as a result of the inherent singularity at the bottom right corner of the laminate where the fully constrained BC at the bottom meets the free BC on the side. Due to their simple geometry, in all the lay-ups compaction of the laminate occurs in a very uniform fashion. Figure 5.46 shows the original distribution of resin volume fraction at sections AB and CD (constant at 0.442) along with the final distributions of volume fractions at those sections. Due to uniformity of the compaction, no discernible differences are observed between the profiles at sections AB and CD. Figure 5.47 to Figure 5.49 present the final distribution of σx over the laminate for [0°], [45°], and [90°] lay-ups respectively. It is evident that the oscillation phenomenon that was discussed earlier becomes more and more pronounced if we move from [0°] lay-up to the case of [90°] lay-up. Figure 5.50 to Figure 5.52 show the profile of σx at section AB of the laminate obtained by the three approaches for [0°], [45°], and [90°] lay-ups respectively. In Figure 5.50, it is observed that the two stress model approaches lead to a  135  Chapter 5: Numerical Applications  rather large change in the amount of developed stress in the x direction in the neighbourhood of the constrained boundary. These sudden changes are accompanied by some oscillation that dies down quite rapidly when we move away from the constrained boundary. The sudden development of stress in the vicinity of the boundary occurs because the constraint on the bottom surface effectively constrains the fibres located at the bottom of the laminate from any movement in the x direction. This combined with the thermal expansion of resin prior to gelation leads to significant development of stress in that region. The response obtained from the integrated model is different in that this boundary layer effect is not observed. The predictions of the stress model with the final volume fractions agrees very well with the one obtained from the integrated model if we move away from the constrained boundary. This shows that the difference between the two results (which is located in the vicinity of the constrained bottom surface) is in fact due to stresses that are predicted by the stress model to develop before the nominal gelation of resin while in the integrated model no significant stresses develop since excess amount of resin can flow out of the system when it undergoes thermal expansion. Similar observations can be made for the case of the [45°] layup as the high stiffness of the system along the x-axis leads to a very similar constrained situation to the [0°] layup (see Figure 5.51). Figure 5.52 does not show a significant difference in the predicted stresses responses by the three modeling approaches applied to the [90°] layup. There are no signs of sudden changes in σx using the two stress models as the fibre directions are transverse to the x-axis. Figure 5.53 to Figure 5.58 present the final distribution of axial forces and bending moments along the length of the laminate for the [0°], [45°] and [90°] layups respectively. As a general rule, the results predicted by the stress model with final volume fractions are closer to those predicted by the integrated model than the regular stress model approach. This is especially more visible in the axial force diagrams. In the case of [90°] layup, there is not much difference between the axial force and bending moment diagrams predicted by the three approaches. The only difference between the results of the integrated model and the two stress models is in the shape of the oscillations in the stress distribution (seeFigure 5.57 and Figure 5.58). This difference could be attributed to the flow of resin through and out of the system in the integrated approach. 136  Chapter 5: Numerical Applications  5.2.3. Flow-compaction and stress development in angle laminates Here, we basically extend upon the previous study of flow-compaction response of angle laminates in Chapter 4, and investigate the interaction of flow-compaction with stress development in the composite to predict the final stress distribution of the angle laminate undergoing cure. The material is assumed to be AS4/3501-6, and the geometry and BC of the problem are presented in Figure 5.59 and Figure 5.60 for the convex and concave cases respectively. The mechanical properties of the fibres and resin are presented in Table 5.10. Figure 5.43 shows the time-history of autoclave temperature. It also presents the assumed time-history of the part and the resulting degree of cure of resin throughout the processing. Similar to the previous examples of this section, the response obtained by three approaches will be analysed and compared here; the integrated model, the stress model using initial distribution of resin volume fraction, and finally the stress model using the final distribution of resin volume fraction obtained from the integrated model. Figure 5.61 shows the location of key points that are considered as critical points where the response of the system is studied and compared for the three different approaches. This figure also presents the transformation of the physical coordinates into a rectilinear ξη coordinate system in which the geometry can be presented in simplified form as a rectangle. Five different meshes of 4-1 elements are used to model every case in this example using the integrated approach. From coarse to fine as presented in Figure 5.62, these meshes are 16×8, 20×12, 24×16, 24×16a, and 24×16b. Due to the limitations of the MATLAB code, it is not possible to model meshes much finer than 24×16 elements at this time. Since, as we will see later, the top and bottom boundaries of the laminate undergo considerable amounts of deformation and stress gradients we opted for refining the 24×16 mesh in the vicinity of the top and bottom boundaries. 16×8, 20×12, 24×16 are all uniformly meshed through the thickness of the laminate. 24×16a mesh is only different from 24×16 mesh in the through thickness direction as the first 2 rows of elements on the top and bottom are 40% the thickness of the elements in the uniform 24×16 mesh. 24×16b is even finer at the top and bottom boundaries and the above ratio is reduced to 20% the thickness of the  137  Chapter 5: Numerical Applications  elements in 24×16 mesh. In all the cases, the results obtained from these five meshes were compared to ensure the convergence of the numerical response. Time-step sizes of 60, 12, and 6 seconds were used and the convergence of results with refining the timesteps was confirmed. The time-step size of 6 seconds was chosen to report the results and comparisons. Figure 5.63 compares the final distribution of resin volume fraction along the EF cross section located at the interface of the curved part and the flat part of the convex [0°] laminate. The convergence of the distribution of resin volume fraction, which is a good representation of the deformation field, in this case is quite evident. Similar comparisons were made for other layups and also for the concave cases and the convergence was verified. The 24×16a mesh was chosen to report the various results from all the three approaches. Figure 5.64 to Figure 5.66 show the distribution of resin volume fraction, υ, over the convex [0°] angle laminate at 3 different times. The combination of these surfaces indicates how the flow of resin through the fibre bed has occurred during the processing. After 14 minutes, only the elements located at the very top of the specimen have seen some resin flow out leading to a decrease in resin volume fraction. At 44 minutes, Figure 5.65 presents a different picture where two flow patterns may be observed; one is the flow of resin in the normal direction to the tool and out of the laminate, and the second is flow of resin from the flat part of the laminate to the curved part and out of the system through that region. The reason for the latter pattern of flow is that the placement of very stiff fibres at the curved part leads to quite a large stiffness normal to the direction of the curved part as any normal displacement would induce an axial strain in the fibres. The final distribution of resin volume fraction is presented in Figure 5.66 where it is evident that the flat part has achieved a uniform compaction while the curved part shows a nonuniform compaction with much higher amounts of resin volume fraction. The final distribution in this example is obtained at the end of the analysis (344 minutes), but after about 70 minutes most of the compaction is completed, and the volume fractions remain almost constant except for small changes due to flow of resin caused by thermal effects before gelation.  138  Chapter 5: Numerical Applications  Figure 5.67 shows the final distribution of resin volume fractions for the case of convex [45°] angle laminate. This surface is very similar to the one in Figure 5.66 because as far as flow is concerned [0°] and [45°] layups are quite similar in 2D plane strain as porous structures with one direction much stiffer than the other direction. Figure 5.68 presents the same distribution but for the convex [90°] angle laminate. In this case, we are practically dealing with an isotropic fibre bed and therefore the laminate compacts quite uniformly and the resin volume fraction reaches values very close to 0.3 at all locations. Figure 5.69, Figure 5.70, and Figure 5.71 present the final distribution of tangential stress (the component of stress in the direction of the laminate) for the convex example with [0°], [45°], and [90°] layups respectively. Stress distributions follow a similar pattern for [0°] and [45°] layups but the values are quite different. As expected, the [90°] layup shows much less variance and lower absolute values for stresses than the other two layups. Figure 5.72, Figure 5.73, and Figure 5.74 show the final distribution of tangential stress obtained by the stress model for the convex example with [0°], [45°], and [90°] layups respectively. Comparing Figure 5.69 with its counterpart, Figure 5.72, the first difference that is observed is the sudden jump in the stress values in the ultimate row of elements at the bottom of the laminate by the tool side in Figure 5.72 that is completely eliminated in Figure 5.69. In the results obtained from the stress model, initially due to the increase in temperature the resin undergoes thermal expansion. As the composite material is considered a solid material as a whole throughout the analysis, it undergoes thermal expansion considering the thermal volumetric changes in the fibres are negligible compared to those of the resin. At the constrained boundary between the composite and the (assumed rigid) tool, the situation becomes very similar to that in Example 6. In the direction of the fibres the system is practically constrained but any strain in the transverse direction is allowed, so the composite material expands in the normal direction to the boundary. The Poisson effect would have led to a small contraction in the longitudinal direction had it not been for the effect of the constrained boundary. Therefore tensile stresses form in the direction of the fibres at a small boundary layer immediate to the constrained boundary which is exactly what is observed in Figure 5.72. This is not the 139  Chapter 5: Numerical Applications  case for the integrated model, as early on resin is free to flow with respect to the fibre-bed and the system does not see any thermal strains due to the thermal expansion of resin prior to its gelation. Figure 5.69 attests to the above argument as there is no sign of a boundary layer at the interface with the tool. The second difference between the two figures is observed in yet another boundary layer effect this time occurring in the stresses (obtained from the integrated model) at the top permeable boundary. Figure 5.69 shows considerably larger values for the stress than those obtained by the stress model at a thin boundary layer on top of the laminate. This effect is not observed on the flat end of the composite and it becomes more pronounced when we move toward and into the curved part. This is attributed to the effect of normal compaction of the laminate due to resin flow and its effect on the fibres in the form of negative longitudinal strains as they move to smaller radii. In the case of [45°] layup, comparing Figure 5.70 and Figure 5.73 leads us to similar differences as were mentioned above for the [0°] layup. the reasoning is similar when we consider that [45°] layup is still extremely anisotropic but with a slightly smaller stiffness in the tangential direction compared to the [0°] layup. As a result, the magnitudes of the boundary layer effects on the top and bottom boundaries are less than the [0°] case which is in harmony with the whole distribution of stresses in [45°] layup if compared with the [0°] layup. For the [90°] layup, comparing Figure 5.71 and Figure 5.74 leads us to no major differences in the distribution of stresses. This makes sense intuitively as all the arguments explaining the above-mentioned boundary layer effects are not valid here as they were inherently based on the anisotropy caused by the directions of the fibres while [90°] layup is fully isotropic in 2D plane strain. In order to have a better understanding of these two boundary layer effects, let‟s look at the distribution of developed stresses just before the degree of cure reaches 0.5 (nominal gelation point) which in this case occurs between 145 and 146 minutes. Figure 5.75 presents the distribution of tangential stresses in the convex [0°] laminate at 144 minutes obtained by the integrated model. The dip in stress in the close proximity of the top 140  Chapter 5: Numerical Applications  permeable boundary is very clear as no significant stress has developed so far in the domain. The amount of compressive stress is clearly reduced and dies down as we move farther from the curved part. Figure 5.76 shows the equivalent results obtained from the stress model. Here, the jump in stress in the first layer of elements located at the interface of the laminate with the rigid tool is clearly highlighted as again there is no significant development of stress if we move away from this boundary. In both of these cases, the localization of the stresses is attributed to the extreme anisotropy in the behaviour of the material due to the high stiffness of fibres in one direction and insignificant shear modulus to carry the stress deeper into the thickness of the composite laminate. Figure 5.77 shows the history of development of tangential stresses at 4 different points on the domain of the convex [0°] laminate. These points, as shown in Figure 5.61, are located at the top and bottom of the symmetric line at the corner and the top and bottom of the centre of the flat part of the laminate. The reported values of stress in this graph and similar figures are from the nearest Gauss point on the right-hand-side of these points. Figure 5.78 compares the time-history of development of tangential stress at point C obtained by three approaches; integrated model, stress model, and stress model using the final distribution of volume fractions obtained from the integrated model. The two sets of results obtained from the stress model show significant amounts of stress development during the pre-gelation period while the integrated model does not. After gelation at around 145 minutes, the stress time history of the integrated model behaves as an offset of the stress model with final volume fractions. Figure 5.79 shows the same comparison as above but for point D. Here, the initial and final volume fractions lead to drastically different trends in the development of stresses especially during the cool-down period. Using this graph, it is very easy to show that the post gelation trend of stress development predicted by the integrated model follows very closely the behaviour observed by the stress model using the final volume fractions. Figure 5.80 and Figure 5.81 present the same comparison as above but for points A and B respectively. As these points are located at the top permeable boundary of the laminate, the integrated model predicts consistent growth of compressive stresses during the time that resin flow is significant. The stress model on the other hand predicts no discernible development of tangential stresses before gelation. 141  Chapter 5: Numerical Applications  Figure 5.82 presents the history of development of tangential stresses at 4 different points on the domain of the convex [45°] laminate. Figure 5.83 to Figure 5.86 compare the results obtained by the three approaches for the time-history of development of tangential stresses at points C, D, A, and B respectively. The comparisons here are very similar to the case of [0°] layup. Figure 5.87 presents the history of development of tangential stresses at 4 different points on the domain of the convex [90°] laminate. Figure 5.88 to Figure 5.91 compare the results obtained by the three approaches for the time-history of development of tangential stresses at points C, D, A, and B respectively. It is evident that in the case of [90°], there is no substantial difference between the responses obtained by the 3 different approaches. For this case, we also show the same comparison in Figure 5.92 for point F (which represents the maximum amount of difference between the obtained results) located at the interface of the flat part and curved part of the laminate. Even at this point, the predictions are indeed very close. Figure 5.93 to Figure 5.95 show the final distribution of the tangential stress in the domain of the concave laminate for [0°],[45°], and [90°] layups respectively. The stress concentration in the case of concave laminates occurs at point C which is locateed at the inside of the corner and adjacent to the permeable boundary. The values of the stress concentration are much higher for all the layups compared to their convex counterparts. Even in the case of [90°] layup, where in the convex sample no concentration of stress was observed, the concave case shows a considerable amount of increase in the final value of tangential stress at point C. Figure 5.96 presents the time-history of development of tangential stresses at points C,D,A and B for the concave [0°] laminate. Due to its sheer magnitude, the development of stress at point C overshadows the history of stress development at the other points. Figure 5.97 to Figure 5.100 show and compare the timehistory of tangential stress development obtained by the three approaches at points C, D, A, and B respectively. In terms of differences between the predicted results, similar observations to the ones made for the convex case can be made here. The only difference here, is that the top and bottom boundary conditions are switched in the concave example as the bottom surface is permeable and the top is rigidly constrained.  142  Chapter 5: Numerical Applications  Figure 5.101 to Figure 5.112 present the final distribution of axial forces and bending moments along the length of the laminate for the 6 combination of convex/concave geometry and [0°]/[45°]/[90°] layup. Bending moment diagrams are essential in the prediction of the final deformed shape of the laminate after it is removed from the tool. In all the cases, the curved geometry of the laminate at the vicinity of the corner plays an integral role in the development of residual stress and eventually defining how the axial force and bending moment diagrams look like. As a general rule, the results predicted by the stress model with final volume fractions are closer to the results predicted by the integrated model than the regular stress model approach. This is especially more visible in the axial force diagrams. Also in the concave samples the above point is more noticeable as compared to the convex cases, the stress concentration at point C overshadows the boundary layer differences between the results obtained by the integrated model and the stress model. The fully constrained boundary at the tool-side leads to some local oscillations in the numerical predictions of tangential stresses at the end of flat part of the laminate. This effect is especially more visible in the case of [45°] layups, and considerably more so for [90°] cases. These oscillations can be seen for instance in Figure 5.74, Figure 5.95, and Figure 5.106. The pattern of these oscillations are similar for the two responses predicted by the stress model but different for the response obtained by the integrated model. This point is showcased more clearly for the case of [90°] layup in Figure 5.105, Figure 5.106, Figure 5.111, and Figure 5.112.  143  Chapter 5: Numerical Applications  5.3. Tables Table 5.1: The comparison and convergence of strain energy values (in N.mm) for the isotropic sample shown in Figure 5.2, using 9-3 elements  Mesh 2×1 3×2 6×4 9×6  Resin velocity conditions  Strain energy (N.mm) CoeffK = 1  0.1  0.01  0.001  0.0001  constrained  2.9627  2.9653  2.9656  2.9657  2.9657  free  2.9627  2.9654  2.9658  2.9658  2.9658  constrained  2.9647  2.9652  2.9652  2.9652  2.9652  free  2.9647  2.9654  2.9655  2.9655  2.9655  constrained  2.9651  2.9651  2.9651  2.9651  2.9651  free  2.9651  2.9657  2.9657  2.9657  2.9657  constrained  2.9651  2.9651  2.9651  2.9651  2.9651  free  2.9651  2.9659  2.966  2.966  2.966  Table 5.2: The comparison and convergence of uA (×10-7 m) for the isotropic sample shown in Figure 5.4, using 9-3 elements  Mesh  0.1  0.01  0.001  0.0001  0.00001  constrained  1.4784  1.6404  1.6835  1.6887  1.6893  1.6893  free  1.4784  1.6404  1.6835  1.6887  1.6893  1.6893  constrained  1.5956  1.6381  1.6484  1.6497  1.6498  1.6498  free  1.5957  1.6382  1.6485  1.6498  1.6499  1.6499  constrained  1.6229  1.6428  1.6477  1.6483  1.6484  1.6484  3×2 6×4  1.5749  9×6 12×8 24×16  uA (×10-7 m)  uA(×10-7m) Resin velocity conditions CoeffK = 1 ABAQUS  free  1.6229  1.6428  1.6477  1.6483  1.6484  1.6484  constrained  1.6326  1.6449  1.648  1.6484  1.6484  1.6484  free  1.6328  1.645  1.6481  1.6485  1.6486  1.6486  1.6404  48×32  1.646  96×64  1.6482  192×128  1.649  144  Chapter 5: Numerical Applications Table 5.3: The comparison and convergence of strain energy values (in N.mm) for the isotropic sample shown in Figure 5.4, using 9-3 elements  Mesh 3×2 6×4 9×6 12×8  Resin velocity conditions  Strain energy (N.mm) CoeffK = 1  0.1  0.01  0.001  0.0001  0.00001  constrained  1.209  1.2409  1.2496  1.2507  1.2508  1.2508  free  1.209  1.241  1.2496  1.2507  1.2508  1.2508  constrained  1.2366  1.25  1.2535  1.2539  1.254  1.254  free  1.2365  1.2499  1.2535  1.2539  1.2539  1.254  constrained  1.2447  1.2525  1.2545  1.2548  1.2548  1.2548  free  1.2447  1.2526  1.2547  1.2549  1.2549  1.2549  constrained  1.2483  1.2536  1.255  1.2551  1.2552  1.2552  free  1.2481  1.2535  1.2549  1.2551  1.2551  1.2551  Table 5.4: The comparison and convergence of uA (×10-7 m) for the isotropic sample shown in Figure 5.4, using 9-4 elements  Mesh 6×4  uA(×10-7m) Boundary ABAQUS Conditions  uA (×10-7 m) CoeffK=1  0.01  0.001  0.0001  1×10-5  1×10-6  1×10-7  Impermeable  1.5956  1.623  1.6276  1.6281  1.6282  1.6282  1.6282  1.6282  Permeable  1.685  2.1044  2.3458  2.3832  2.3871  2.3875  2.3876  2.3876  Impermeable  1.6326  1.6404  1.6419  1.6421  1.6421  1.6421  1.6421  1.6421  Permeable  1.677  1.8907  1.9894  2.0032  2.0047  2.0048  2.0048  2.0048  1.5749  12×8 24×16  0.1  1.6404  48×32  1.646  96×64  1.6482  192×128  1.649  Table 5.5: The comparison and convergence of strain energy values (in N.mm) for the isotropic sample shown in Figure 5.4, using 9-4 elements  Mesh  Boundary Conditions  Strain energy (N.mm) CoeffK= 1  0.1  0.01  0.001  0.0001  1×10-5  1×10-6  1×10-7  Impermeable  1.2366  1.2452  1.2469  1.2471  1.2472  1.2472  1.2472  1.2472  Permeable  1.159  1.1765  1.1931  1.1959  1.1962  1.1962  1.1962  1.1962  Impermeable  1.2483  1.2517  1.2524  1.2524  1.2525  1.2525  1.2525  1.2525  Permeable  1.2085  1.2204  1.2279  1.2291  1.2292  1.2292  1.2292  1.2292  6×4 12×8  145  Chapter 5: Numerical Applications Table 5.6: The comparison and convergence of strain energy and uA (×10-7 m) for the transversely isotropic sample (   0 ) shown in Figure 5.4, using 9-3 elements  Mesh  ABAQUS (disp.)  3×2 6×4  9.7984  9×6 12×8  CoeffK  Output type  1  0.1  0.01  0.001  0.0001  0.00001  displacement  9.6715  9.9541  10.0153  10.0223  10.023  10.0231  strain energy  1.1715  1.1956  1.201  1.2016  1.2017  1.2017  displacement  9.8174  9.8784  9.8915  9.893  9.8932  9.8932  strain energy  1.1877  1.1952  1.1968  1.197  1.197  1.197  displacement  9.8434  9.868  9.8732  9.8738  9.8739  9.8739  strain energy  1.1912  1.1948  1.1956  1.1957  1.1957  1.1957  displacement  9.8518  9.8651  9.8679  9.8682  9.8683  9.8683  strain energy  1.1925  1.1947  1.1952  1.1952  1.1952  1.1952  24×16 48×32 96×64  9.8614  θ = 0° The displacement values are in 10-8 m, and the strain energy values are in 10-4 N.m.  192×128  Table 5.7: The comparison and convergence of strain energy and uA (×10-7 m) for the transversely isotropic sample (   90 ) shown in Figure 5.4, using 9-3 elements  Mesh  ABAQUS (disp.)  3×2 6×4  1.0485  9×6 12×8  CoeffK  Output type  1  0.1  0.01  0.001  displacement  0.8639  1.3198  1.4509  strain energy  4.4961  4.6509  4.6917  displacement  1.0823  1.1804  1.2037  0.0001  0.00001  1.4673  1.469  1.4692  4.6967  4.6972  4.6972  1.2065  1.2068  1.2069  strain energy  4.6062  4.6673  4.6831  4.685  4.6852  4.6852  displacement  1.131  1.1737  1.1839  1.1851  1.1852  1.1853  strain energy  4.6367  4.6712  4.6801  4.6811  4.6813  4.6813  displacement  1.1481  1.1729  1.179  1.1797  1.1798  1.1798  strain energy  4.6499  4.6727  4.6786  4.6793  4.6794  4.6794  24×16 48×32 96×64  1.1749  θ = 90° The displacement values are in 10-8 m, and the strain energy values are in 10-4 N.m.  192×128  146  Chapter 5: Numerical Applications Table 5.8: The comparison and convergence of strain energy and uA (×10-7 m) for the transversely isotropic sample (   0 ) shown in Figure 5.4, using 9-3 elements  Mesh  ABAQUS (disp.)  3×2 8.9136  6×4 9×6 12×8  CoeffK  Output type  1  0.1  0.01  0.001  0.0001  0.00001  displacement  8.7336  9.1036  9.1848  9.1941  9.195  9.1951  strain energy  1.9287  1.9672  1.976  1.977  1.9771  1.9771  displacement  8.942  9.0248  9.0428  9.0448  9.045  9.0451  strain energy  1.9574  1.9713  1.9744  1.9748  1.9748  1.9748  displacement  8.9808  9.0145  9.0218  9.0227  9.0227  9.0227  strain energy  1.9645  1.9718  1.9735  1.9736  1.9737  1.9737  displacement  8.9935  9.012  9.016  9.0164  9.0165  9.0165  strain energy  1.9674  1.972  1.973  1.9731  1.9731  1.9731  24×16 48×32 9.0092  96×64  θ = 45° The displacement values are in 10-8 m, and the strain energy values are in 10-4 N.m.  192×128  Table 5.9: Resin properties in Examples 1 & 2, where Poisson‟s ratio is assumed constant  Property  3501-6 resin  u  3.2  r  E (GPa)  0.031  ν  0.35  E (GPa)  Table 5.10: Mechanical properties of fibre and resin in the examples of Section 5.2  Property  AS4 fibre  Property  3501-6 resin   12  0.2  Ku (GPa)  3.556   23  0.3  Kr (GPa)  0.8  E1 (GPa)  207  Gu (GPa)  1.185  E2 (GPa)  20.7  Gr (kPa)  1  G12 (GPa)  27.6  CTEg (/°C)  5.76×10-5  CTE1 (/°C)  -0.9×10-6  CTEr (/°C)  1.854×10-4  CTE2 (/°C)  7.2×10-6  CSC (volumetric)  -0.05488  147  Chapter 5: Numerical Applications  5.4. Figures  y f0  h  x L Figure 5.1: The geometry and specifications of the composite sample under constant pressure load (uni-axial strain condition)  y f2  f1  h  x L Figure 5.2: The geometry and specifications of the composite sample under linear pressure load 148  Chapter 5: Numerical Applications  -2.E-04  -3.E-04  -4.E-04  uy (mm)  -5.E-04  1×2 mesh 3×2 mesh  -6.E-04  6×4 mesh -7.E-04  -8.E-04  -9.E-04  -1.E-03 0  2  4  6  8  10  12  14  x (mm)  Figure 5.3: The profile of the vertical displacement at the top of the constrained sample under linear distribution of pressure load  z f2  f1 f1 = 700 kPa f2 = 100 kPa 2h h = 4 mm L = 12 mm  Gc  0.7837GPa  A K  3.205GPa c  x    0.4  L Figure 5.4: Composite sample under pressure gradient  149  Chapter 5: Numerical Applications  4.5  4 Left side's displacement profile 3.5  Right side's displacement profile  y (mm)  3  2.5  2  1.5  1  0.5  0 0.00E+00  2.00E-07  4.00E-07  6.00E-07  8.00E-07  1.00E-06  1.20E-06  1.40E-06  ux (m)  Figure 5.5: Shape of the displacement profiles of the sample under pressure gradient  4.5  4  3×2 mesh/ 9-3 elements/ Coeff.K=1 3.5  3×2 mesh/ 9-3 elements/ Coeff.K=0.1 3×2 mesh/ 9-3 elements/ Coeff.K=0.01  y (mm)  3  3×2 mesh/ 9-3 elements/ Coeff.K=0.0001  2.5  2  1.5  1  0.5  0 0.00E+00  2.00E-08  4.00E-08  6.00E-08  8.00E-08  1.00E-07  1.20E-07  1.40E-07  1.60E-07  1.80E-07  ux (m)  Figure 5.6: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for different values of CoeffK, obtained using a 3×2 mesh of 9-3 elements 150  Chapter 5: Numerical Applications  4.5  4  6×4 mesh/ 9-3 elements/ Coeff.K=1 3.5  6×4 mesh/ 9-3 elements/ Coeff.K=0.1 6×4 mesh/ 9-3 elements/ Coeff.K=0.01  y (mm)  3  6×4 mesh/ 9-3 elements/ Coeff.K=0.0001  2.5  2  1.5  1  0.5  0 0.00E+00  2.00E-08  4.00E-08  6.00E-08  8.00E-08  1.00E-07  1.20E-07  1.40E-07  1.60E-07  1.80E-07  ux (m)  Figure 5.7: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for different values of CoeffK, obtained using a 6×4 mesh of 9-3 elements  4.5  4  12×8 mesh/ 9-3 elements/ Coeff.K=1 3.5  12×8 mesh/ 9-3 elements/ Coeff.K=0.1 12×8 mesh/ 9-3 elements/ Coeff.K=0.01  y (mm)  3  12×8 mesh/ 9-3 elements/ Coeff.K=0.0001  2.5  2  1.5  1  0.5  0 0.00E+00  2.00E-08  4.00E-08  6.00E-08  8.00E-08  1.00E-07  1.20E-07  1.40E-07  1.60E-07  1.80E-07  ux (m)  Figure 5.8: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for different values of CoeffK, obtained using a 12×8 mesh of 9-3 elements 151  Chapter 5: Numerical Applications  4.5  4  3.5 9-3 elements/ 3×2 mesh 9-3 elements/ 6×4 mesh  y (mm)  3  9-3 elements/ 12×8 mesh  2.5  2  1.5  1  0.5  0 0.0E+00  2.0E-08  4.0E-08  6.0E-08  8.0E-08  1.0E-07  1.2E-07  1.4E-07  1.6E-07  1.8E-07  ux (m)  Figure 5.9: Comparison of displacement profiles at the right side of the sample in Figure 5.4 using different meshes of 9-3 elements, in the case where CoeffK=0.01  4.5  4  3.5 9-3 elements/ 3×2 mesh 9-3 elements/ 6×4 mesh  y (mm)  3  9-3 elements/ 12×8 mesh  2.5  2  1.5  1  0.5  0 0.0E+00  2.0E-08  4.0E-08  6.0E-08  8.0E-08  1.0E-07  1.2E-07  1.4E-07  1.6E-07  1.8E-07  ux (m)  Figure 5.10: Comparison of displacement profiles at the right side of the sample in Figure 5.4 using different meshes of 9-3 elements, in the case where CoeffK=1 152  Chapter 5: Numerical Applications  4.5  4 9-3 elements  3.5  9-4 elements/Permeable B.C. 9-4 elements/Impermeable B.C.  y (mm)  3  2.5  2  1.5  1  0.5  0 0.0E+00  5.0E-08  1.0E-07  1.5E-07  2.0E-07  2.5E-07  ux (m)  Figure 5.11: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for 9-3 and 9-4 elements, obtained using a 6×4 mesh  4.5  4 9-3 elements  3.5  9-4 elements/Permeable B.C. 9-4 elements/Impermeable B.C.  y (mm)  3  2.5  2  1.5  1  0.5  0 0.0E+00  5.0E-08  1.0E-07  1.5E-07  2.0E-07  2.5E-07  ux (m)  Figure 5.12: Comparison of displacement profiles at the right side of the sample in Figure 5.4 for 9-3 and 9-4 elements, obtained using a 12×8 mesh 153  Chapter 5: Numerical Applications  5  x 10 7  6  5  P (Pa)  4  3  2  1  0 0 1 2 3  -3  x 10  4  0.012  0.01  0.008  0.006  0.004  0.002  0  x (m)  y (m)  Figure 5.13: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-3 elements  5  x 10 8 7 6  P (Pa)  5 4 3 2 1 0 0 1 2 3  -3  x 10  4 y (m)  0.012  0.01  0.008  0.006  0.004  0.002  0  x (m)  Figure 5.14: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-4 elements with permeable BC on the sides 154  Chapter 5: Numerical Applications  5  x 10 7 6 5  P (Pa)  4 3 2 1 0 -1 0 1 2 3  -3  x 10  4  0.012  0.01  0.008  0.006  0.004  0.002  0  x (m)  y (m)  Figure 5.15: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-4 elements with impermeable BC on the sides  5  x 10 7 6 5  P (Pa)  4 3 2 1 0 -1 0 1 2 3  -3  x 10  4 y (m)  0.012  0.01  0.008  0.006  0.004  0.002  0  x (m)  Figure 5.16: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-4 elements based on the u-v-p formulation  155  200  1  180  0.9  160  0.8  140  0.7  120  0.6  100  0.5  Temperature 80  0.4  Degree of cure  60  0.3  40  0.2  20  0.1  0 0  50  100  150  200  Degree of cure  Temperature (°C)  Chapter 5: Numerical Applications  0 250  time (min)  Figure 5.17: Time-history of the applied temperature and the resulting degree of cure in Examples 1-7  200  1E+3  180 160 1E+2  120 100  1E+1  Temperature Resin viscosity  80  Viscosity (Pa.s)  Temperature (°C)  140  60 1E+0 40 20 0 0  50  100  150  200  1E-1 250  time (min)  Figure 5.18: Time-history of the applied temperature and the resulting resin viscosity in Examples 1-7 156  Chapter 5: Numerical Applications  50 45  σx (MPa)  40 35  Integrated model  30  Stress model (Zobeiry)  25 20 15 10 5 0 0  50  100  150  200  250  time (min)  Figure 5.19: Time-history of the development of longitudinal stress in Example 1. The predictions are essentially identical.  80 70 60 50  σz (MPa)  Integrated model  40  Stress model (Zobeiry)  30 20 10 0 0  50  100  150  200  250  -10 time (min)  Figure 5.20: Time-history of the development of transverse stress in Examples 1 and 2. The predictions are essentially identical. 157  Chapter 5: Numerical Applications  60  50  σx (MPa)  40 Integrated model Stress model (Zobeiry)  30  20  10  0 0  50  100  150  200  250  time (min)  Figure 5.21: Time-history of the development of  x in Example 2. The predictions are essentially identical.  40  30  20  σx (MPa)  Integrated model Stress model  10  0 0  50  100  150  200  250  -10  -20 time (min)  Figure 5.22: Time-history of the development of  x in Example 3. The predictions are essentially identical. 158  Chapter 5: Numerical Applications  60  40 Integrated model  σz (MPa)  20  Stress model  0 0  50  100  150  200  250  -20  -40  -60 time (min)  Figure 5.23: Time-history of the development of  z in Examples 3 and 4. The predictions are essentially identical.  60 50 40 30  P (MPa)  20 10 0 0  50  100  150  200  250  -10 -20 -30 -40 time (min)  Figure 5.24: Time-history of the development of resin pressure in Examples 3 and 4  159  Chapter 5: Numerical Applications  40  30  Integrated model  20  σx (MPa)  Stress model 10  0 0  50  100  150  200  250  -10  -20  -30 time (min)  Figure 5.25: Time-history of the development of  x in Example 4. The predictions are essentially identical.  8.E-05 7.E-05  Resin velocity (mm/s)  6.E-05 5.E-05 4.E-05 3.E-05 2.E-05 1.E-05 0.E+00 0  50  100  150  200  250  time (min)  Figure 5.26: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 5 160  200  8.0E-4  175  7.0E-4  150  6.0E-4  125  5.0E-4  100  4.0E-4  75  3.0E-4  50  2.0E-4  Temperature Rate of free vol. strain  25  Rate of free volumetric strain (/min)  Temperature (°C)  Chapter 5: Numerical Applications  1.0E-4  0  0.0E+0 15  20  25  30  35  40  45  50  55  time (min)  Figure 5.27: Time-history of the applied temperature and the changes in the rate of volumetric strain during t=15-55 min in Example 5  0.44  φ  0.42  0.4  Integrated model  0.38  Stress model  0.36 0  50  100  150  200  250  time (min)  Figure 5.28: Time-history of the changes in resin volume fraction in Example 5  161  Chapter 5: Numerical Applications  10 0 0  50  100  150  200  250  -10  P (MPa)  -20 -30 -40 -50 -60 -70 -80 time (min)  Figure 5.29: Time-history of the development of resin pressure in Example 5  50  40  30  Integrated model  σx (MPa)  Stress model  20  10  0 0  50  100  150  200  250  -10  -20 time (min)  Figure 5.30: Time-history of the development of  x in Example 5  162  Chapter 5: Numerical Applications  100 80 Integrated model  60  Stress model  σz (MPa)  40 20 0 0  50  100  150  200  250  -20 -40 -60 time (min)  Figure 5.31: Time-history of the development of  z in Example 5  1.6E-04 1.4E-04  Resin velocity (mm/s)  1.2E-04 1.0E-04 8.0E-05 6.0E-05 4.0E-05 2.0E-05 0.0E+00 0  50  100  150  200  250  time (min)  Figure 5.32: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 6 163  Chapter 5: Numerical Applications  0.44  φ  0.42  0.4  Integrated model  0.38  Stress model  0.36 0  50  100  150  200  250  time (min)  Figure 5.33: Time-history of the changes in resin volume fraction in Example 6  0.03 0.025 0.02 0.015  P (MPa)  0.01 0.005 0 0  50  100  150  200  250  -0.005 -0.01 -0.015 -0.02 time (min)  Figure 5.34: Time-history of the development of resin pressure in Example 6  164  Chapter 5: Numerical Applications  0.04  0.03  Integrated model Stress model  0.02  εz  0.01  0 0  50  100  150  200  250  -0.01  -0.02  -0.03 time (min)  Figure 5.35: Time-history of the transverse strain in Example 6  25  20  15  σx (MPa)  Integrated model Stress model  10  5  0 0  50  100  150  200  250  -5 time (min)  Figure 5.36: Time-history of the development of  x in Example 6  165  Chapter 5: Numerical Applications  0.6  0.5  P (MPa)  0.4  0.3  0.2  0.1  0 0  50  100  150  200  250  time (min)  Figure 5.37: Time-history of the development of resin pressure in Example 7  1.4E-03  Resin velocity (mm/s)  1.2E-03 1.0E-03 8.0E-04 6.0E-04 4.0E-04 2.0E-04 0.0E+00 0  50  100  150  200  250  time (min)  Figure 5.38: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 7  166  Chapter 5: Numerical Applications  0.44 0.42 0.4 0.38 Integrated model  φ  0.36  Stress model  0.34  Stress model - final φ  0.32 0.3 0.28 0.26 0  50  100  150  200  250  time (min)  Figure 5.39: Time-history of the changes in resin volume fraction in Example 7  0.04  0 0  50  100  150  200  250  εz  -0.04  -0.08  Integrated model Stress model Stress model - final φ  -0.12  -0.16  -0.2 time (min)  Figure 5.40: Time-history of the transverse strain in Example 7  167  Chapter 5: Numerical Applications  25  20  15 Integrated model  σx (MPa)  10  Stress model Stress model - final φ  5  0 0  50  100  150  200  250  -5  -10 time (min)  Figure 5.41: Time-history of the development of  x in Example 7  Permeable BC  f  4.96 mm  z 63 mm  x B A  D C  Figure 5.42: Geometry and BC of the flat laminate undergoing cure  168  200  1.0  180  0.9  160  0.8  140  0.7  120  0.6  100  0.5  80  0.4  60  0.3  Autoclave Temp. Part Temp.  40  Degree of cure  Temperature (°C)  Chapter 5: Numerical Applications  0.2  Degree of cure  20  0.1  0 0  50  100  150  200  250  0.0 350  300  time (min)  Figure 5.43: Time-history of the autoclave and part temperatures and the resulting degree of cure in the flat and angle laminate examples  5 4.5 4 3.5  [0°] Stress model [0°] ABAQUS  z (mm)  3  [90°] Stress Model  2.5  [90°] ABAQUS  2 1.5 1 0.5 0 -5  0  5  10  15 σx (MPa)  20  25  30  35  Figure 5.44: Final profiles of σx predicted by the stress model at section AB of the [0°] and [90°] flat laminates, and their comparison with response obtained from ABAQUS 169  Chapter 5: Numerical Applications  60 Stress model  50  ABAQUS  σx (MPa)  40  30  20  10  0 0  10  20  30  40  50  60  70  x (mm)  Figure 5.45: Final distribution of σx predicted by the stress model along the bottom surface of the [90°] flat laminates in comparison with response obtained from ABAQUS  5 4.5 4 [0°] lay-up  z (mm)  3.5  [45°] lay-up  3  [90°] lay-up  2.5  Initial φ  2 1.5 1 0.5 0 0  0.2  0.4  0.6  0.8  1  φ  Figure 5.46: Comparison of final profiles of resin volume fraction predicted by the integrated model at AB and CD sections for different lay-ups of the flat laminate 170  Chapter 5: Numerical Applications  Figure 5.47: Final distribution of σx for the [0°] flat laminate obtained by the integrated model  Figure 5.48: Final distribution of σx for the [45°] flat laminate obtained by the integrated model  Figure 5.49: Final distribution of σx for the [90°] flat laminate obtained by the integrated model 171  Chapter 5: Numerical Applications  5 4.5 4 3.5  Integrated model  z (mm)  3  Stress model Stress model - final φ  2.5 2 1.5 1 0.5 0 -10  -5  0  5  10  15  20  σx (MPa)  Figure 5.50: Final profile of σx at section AB of the [0°] flat laminate  5 4.5 4  z (mm)  3.5  Integrated model Stress model  3  Stress model - final φ  2.5 2 1.5 1 0.5 0 0  5  10  15  20  25  σx (MPa)  Figure 5.51: Final profile of σx at section AB of the [45°] flat laminate  172  Chapter 5: Numerical Applications  5 4.5 4  z (mm)  3.5  Integrated model  3  Stress model Stress model - final φ  2.5 2 1.5 1 0.5 0 5  10  15  20  25  30  35  σx (MPa)  Figure 5.52: Final profile of σx at section AB of the [90°] flat laminate  15 10  Axial force (N/mm)  Integrated model  5  Stress model Stress model - final φ  0 0  10  20  30  40  50  60  70  -5 -10 -15 -20 x (mm)  Figure 5.53: Final distribution of axial force along the length of the [0°] flat laminate  173  Chapter 5: Numerical Applications  6 Integrated model  4  Stress model Stress model - final φ  Moment (N.mm/mm)  2 0 0  10  20  30  40  50  60  70  -2 -4 -6 -8 -10 x (mm)  Figure 5.54: Final distribution of bending moment along the length of the [0°] flat laminate  90 80  Axial force (N/mm)  70 60 50 40  Integrated model  30  Stress model Stress model - final φ  20 10 0 0  10  20  30  40  50  60  70  x (mm)  Figure 5.55: Final distribution of axial force along the length of the [45°] flat laminate  174  Chapter 5: Numerical Applications  0 0  10  20  30  40  50  60  70  Moment (N.mm/mm)  -5  -10  -15  Integrated model Stress model Stress model - final φ  -20  -25  -30 x (mm)  Figure 5.56: Final distribution of bending moment along the length of the [45°] flat laminate  160 140  Axial force (N/mm)  120 100 Integrated model  80  Stress model  60  Stress model - final φ  40 20 0 0  10  20  30  40  50  60  70  x (mm)  Figure 5.57: Final distribution of axial force along the length of the [90°] flat laminate  175  Chapter 5: Numerical Applications  10 0 0  10  20  30  40  50  60  70  Moment (N.mm/mm)  -10 -20 -30 Integrated model  -40  Stress model Stress model - final φ  -50 -60 -70 -80  x (mm)  Figure 5.58: Final distribution of bending moment along the length of the [90°] flat laminate  f  Permeable B.C.  f L Permeable B.C.  R H  z  H = 4.96 mm L = 63 mm R = 4.57 mm υ = 0.442 f = 540 kPa  x  Figure 5.59: Geometry and boundary conditions of the convex angle laminate  176  Chapter 5: Numerical Applications  f L  R  z  H  x  f  Permeable B.C.  H = 4.96 mm L = 63 mm R = 1.39 mm υ = 0.442 f = 540 kPa  Figure 5.60: Geometry and boundary conditions of the concave angle laminate  B D E A  F C  η  ξ Figure 5.61: Location of key points on the geometry of the angle laminate, and transformation of geometry into the ξη coordinates 177  Chapter 5: Numerical Applications  16×8  20×12  24×16  24×16a  24×16b  Figure 5.62: FE meshes used in the analysis of angle laminates 178  Chapter 5: Numerical Applications  5  4  η (mm)  16×8 mesh  3  20×12 mesh 24×16 mesh 24×16a mesh  2  24×16b mesh  1  0 0.3  0.305  0.31  0.315 φ  0.32  0.325  0.33  Figure 5.63: Final profile of resin volume fraction for the convex [0°] angle laminate at section EF obtained by the integrated model with different meshes  Figure 5.64: Distribution of resin volume fraction for the convex [0°] angle laminate obtained by the integrated model at t = 14 min  179  Chapter 5: Numerical Applications  Figure 5.65: Distribution of resin volume fraction for the convex [0°] angle laminate obtained by the integrated model at t = 44 min  Figure 5.66: Final distribution of resin volume fraction for the convex [0°] angle laminate obtained by the integrated model  180  Chapter 5: Numerical Applications  Figure 5.67: Final distribution of resin volume fraction for the convex [45°] angle laminate obtained by the integrated model  Figure 5.68: Final distribution of resin volume fraction for the convex [90°] angle laminate obtained by the integrated model  181  Chapter 5: Numerical Applications  Figure 5.69: Final distribution of tangential stress for the convex [0°] angle laminate obtained by the integrated model  Figure 5.70: Final distribution of tangential stress for the convex [45°] angle laminate obtained by the integrated model  182  Chapter 5: Numerical Applications  Figure 5.71: Final distribution of tangential stress for the convex [45°] angle laminate obtained by the integrated model  Figure 5.72: Final distribution of tangential stress for the convex [0°] angle laminate obtained by the stress model  183  Chapter 5: Numerical Applications  Figure 5.73: Final distribution of tangential stress for the convex [45°] angle laminate obtained by the stress model  Figure 5.74: Final distribution of tangential stress for the convex [90°] angle laminate obtained by the stress model  184  Chapter 5: Numerical Applications  Figure 5.75: Distribution of tangential stress for the convex [0°] angle laminate obtained by the integrated model at t = 144 min  Figure 5.76: Distribution of tangential stress for the convex [0°] angle laminate obtained by the stress model at t = 144 min  185  Chapter 5: Numerical Applications  40 20 0  σ11 (MPa)  -20 -40 -60 @C  -80  @D  -100  @A  -120  @B  -140 -160 0  50  100  150  200  250  300  350  time (min)  Figure 5.77: Time-history of the development of tangential stress at different locations of the convex [0°] angle laminate  30 25  Integrated model Stress model  σ11 (MPa)  20  Stress model - final φ  15 10 5 0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.78: Time-history of the development of tangential stress at point C on the convex [0°] angle laminate 186  Chapter 5: Numerical Applications  16 14 12  σ11 (MPa)  10 8  Integrated model  6  Stress model Stress model - final φ  4 2 0 -2  0  50  100  150  200  250  300  350  -4 time (min)  Figure 5.79: Time-history of the development of tangential stress at point D on the convex [0°] angle laminate  20 0 0  50  100  150  200  250  300  350  -20  σ11 (MPa)  -40 -60 -80 -100  Integrated model Stress model  -120  Stress model - final φ  -140 -160 time (min)  Figure 5.80: Time-history of the development of tangential stress at point A on the convex [0°] angle laminate 187  Chapter 5: Numerical Applications  0 0  50  100  150  200  250  300  350  -5  σ11 (MPa)  -10  -15  -20  Integrated model Stress model Stress model - final φ  -25  -30 time (min)  Figure 5.81: Time-history of the development of tangential stress at point B on the convex [0°] angle laminate  30 20 10 0  σξ (MPa)  -10 -20 -30 @C  -40  @D  -50  @A  -60  @B  -70 -80 0  50  100  150  200  250  300  350  time (min)  Figure 5.82: Time-history of the development of tangential stress at different locations of the convex [45°] angle laminate 188  Chapter 5: Numerical Applications  30 25  Integrated model Stress model  σξ (MPa)  20  Stress model - final φ  15 10 5 0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.83: Time-history of the development of tangential stress at point C on the convex [45°] angle laminate  25  20  Integrated model Stress model Stress model - final φ  σξ (MPa)  15  10  5  0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.84: Time-history of the development of tangential stress at point D on the convex [45°] angle laminate 189  Chapter 5: Numerical Applications  10 0 0  50  100  150  200  250  300  350  -10  σξ (MPa)  -20 -30 -40 Integrated model  -50  Stress model  -60  Stress model - final φ  -70 -80 time (min)  Figure 5.85: Time-history of the development of tangential stress at point A on the convex [45°] angle laminate  15  Integrated model  10  Stress model  σξ (MPa)  Stress model - final φ  5  0 0  50  100  150  200  250  300  350  -5  -10 time (min)  Figure 5.86: Time-history of the development of tangential stress at point B on the convex [45°] angle laminate 190  Chapter 5: Numerical Applications  30 25 @C  σξ (MPa)  20  @D @A  15  @B  10 5 0 -5 0  50  100  150  200  250  300  350  time (min)  Figure 5.87: Time-history of the development of tangential stress at different locations of the convex [90°] angle laminate  25  20 Integrated model  σξ (MPa)  15  Stress model Stress model - final φ  10  5  0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.88: Time-history of the development of tangential stress at point C on the convex [90°] angle laminate 191  Chapter 5: Numerical Applications  30 25 Integrated model  20  Stress model  σξ (MPa)  Stress model - final φ  15 10 5 0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.89: Time-history of the development of tangential stress at point D on the convex [90°] angle laminate  16 14 12 Integrated model  10  σξ (MPa)  Stress model  8  Stress model - final φ  6 4 2 0 0  50  100  150  200  250  300  350  -2 time (min)  Figure 5.90: Time-history of the development of tangential stress at point A on the convex [90°] angle laminate 192  Chapter 5: Numerical Applications  30 25 20  Integrated model  σξ (MPa)  Stress model Stress model - final φ  15 10 5 0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.91: Time-history of the development of tangential stress at point B on the convex [90°] angle laminate  30 25 20  Integrated model  σξ (MPa)  Stress model Stress model - final φ  15 10 5 0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.92: Time-history of the development of tangential stress at point F on the convex [90°] angle laminate 193  Chapter 5: Numerical Applications  Figure 5.93: Final distribution of tangential stress for the concave [0°] angle laminate obtained by the integrated model  Figure 5.94: Final distribution of tangential stress for the concave [45°] angle laminate obtained by the integrated model  194  Chapter 5: Numerical Applications  Figure 5.95: Final distribution of tangential stress for the concave [90°] angle laminate obtained by the integrated model  350 300  σ11 (MPa)  250 200  @C @D  150  @A @B  100 50 0 -50 0  50  100  150  200  250  300  350  time (min)  Figure 5.96: Time-history of the development of tangential stress at different locations of the concave [0°] angle laminate  195  Chapter 5: Numerical Applications  450 400 350 Integrated model  300  σ11 (MPa)  Stress model  250  Stress model - final φ  200 150 100 50 0 -50  0  50  100  150  200  250  300  350  time (min)  Figure 5.97: Time-history of the development of tangential stress at point C on the concave [0°] angle laminate  16 14 12 Integrated model  σ11 (MPa)  10  Stress model  8  Stress model - final φ  6 4 2 0 0  50  100  150  200  250  300  350  -2 time (min)  Figure 5.98: Time-history of the development of tangential stress at point D on the concave [0°] angle laminate 196  Chapter 5: Numerical Applications  30 25  Integrated model Stress model  σ11 (MPa)  20  Stress model - final φ  15 10 5 0 0  50  100  150  200  250  300  350  -5 time (min)  Figure 5.99: Time-history of the development of tangential stress at point A on the concave [0°] angle laminate  16 14 12  σ11 (MPa)  10 8  Integrated model  6  Stress model Stress model - final φ  4 2 0 -2  0  50  100  150  200  250  300  350  -4 time (min)  Figure 5.100: Time-history of the development of tangential stress at point B on the concave [0°] angle laminate 197  Chapter 5: Numerical Applications  0 0  10  20  30  40  50  60  70  -20 -40  Axial force (N/mm)  -60 -80 -100  Integrated model Stress model  -120  Stress model - final φ  -140 -160 -180 -200 distance from corner (mm)  Figure 5.101: Final distribution of axial force along the length of the convex [0°] angle laminate  50  0  Moment (N.mm/mm)  0  10  20  30  40  50  60  70  -50  -100  Integrated model Stress model  -150  Stress model - final φ  -200  -250 distance from corner (mm)  Figure 5.102: Final distribution of bending moment along the length of the convex [0°] angle laminate  198  Chapter 5: Numerical Applications  80  60  Axial force (N/mm)  40  20  0 0  10  20  30  40  50  60  70  -20 Integrated model Stress model  -40  Stress model - final φ  -60 distance from corner (mm)  Figure 5.103: Final distribution of axial force along the length of the convex [45°] angle laminate  0 0  10  20  30  40  50  60  70  -20  Moment (N.mm/mm)  -40 -60 -80  Integrated model Stress model  -100  Stress model - final φ  -120 -140 -160 distance from corner (mm)  Figure 5.104: Final distribution of bending moment along the length of the convex [45°] angle laminate  199  Chapter 5: Numerical Applications  160 140  Axial force (N/mm)  120 100 80  Integrated model Stress model  60  Stress model - final φ  40 20 0 0  10  20  30  40  50  60  70  distance from corner (mm)  Figure 5.105: Final distribution of axial force along the length of the convex [90°] angle laminate  10  0 0  10  20  30  40  50  60  70  Moment (N.mm/mm)  -10  -20  -30  Integrated model Stress model  -40  Stress model - final φ  -50  -60 distance from corner (mm)  Figure 5.106: Final distribution of bending moment along the length of the convex [90°] angle laminate  200  Chapter 5: Numerical Applications  300  250  Axial force (N/mm)  200  Integrated model Stress model Stress model - final φ  150  100  50  0 0  10  20  30  40  50  60  70  -50 distance from corner (mm)  Figure 5.107: Final distribution of axial force along the length of the concave [0°] angle laminate  50 0 0  10  20  30  40  50  60  70  Moment (N.mm/mm)  -50 -100 -150  Integrated model Stress model  -200  Stress model - final φ  -250 -300 -350 distance from corner (mm)  Figure 5.108: Final distribution of bending moment along the length of the concave [0°] angle laminate  201  Chapter 5: Numerical Applications  300  Axial force (N/mm)  250  Integrated model  200  Stress model Stress model - final φ  150  100  50  0 0  10  20  30  40  50  60  70  distance from corner (mm)  Figure 5.109: Final distribution of axial force along the length of the concave [45°] angle laminate  50  0  Moment (N.mm/mm)  0  10  20  30  40  50  60  70  -50 Integrated model  -100  Stress model Stress model - final φ  -150  -200 distance from corner (mm)  Figure 5.110: Final distribution of bending moment along the length of the concave [45°] angle laminate 202  Chapter 5: Numerical Applications  300  250  Integrated model  Axial force (N/mm)  Stress model Stress model - final φ  200  150  100  50  0 0  10  20  30  40  50  60  70  distance from corner (mm)  Figure 5.111: Final distribution of axial force along the length of the concave [90°] angle laminate  80 60  Integrated model Stress model Stress model - final φ  Moment (N.mm/mm)  40 20 0 0  10  20  30  40  50  60  70  -20 -40 -60 -80 distance from corner (mm)  Figure 5.112: Final distribution of bending moment along the length of the concave [90°] angle laminate 203  Chapter 6: Conclusions and Further Work  6.1. Summary and Conclusions The main objective of this dissertation was to provide a framework to integrate two aspects of process modeling, namely, resin flow through the fibre-bed and development of residual stresses and deformations during cure of composite materials into one unique numerical model. To achieve that goal, a two-phase model was chosen to be the basis of the FE simulations, and an investigation was performed into measures that were necessary in order for the model to also be able to simulate the development of stresses in a post-gelled curing composite material. The following conclusions may be drawn from this work: 1. The volume averaging method was used to rigorously derive the governing equations for a general two-phase system consisting of a solid and fluid phase. The effect of shear stresses in the fluid phase were considered along with its pressure. This leads to a u-v-p formulation that is shown to be an ideal formulation for the ultimate goal of integration of flow and stress modeling. 2. Two finite elements, namely 9-3 and 4-1, were introduced and formulated using the obtained governing equations. Special measures were introduced in order to add the capability of modeling jumps in the distribution of volume fractions. Additional considerations were also made for the model to successfully represent the boundary between a viscous fluid and its adjacent porous structure. 3. The finite element formulation was implemented in an in-house MATLAB code. Other elements such as a 9-4 element based on the more common u-p approach were also implemented in a separate MATLAB code for the sake of comparison. The current approach was validated in a few numerical examples, and the response of the 9-3 and 9-4 elements were compared extensively.  204  Chapter 6: Conclusions and Further Work  4. The proposed two-phase approach was adapted so that it could successfully represent the stress development in a curing composite based on the pseudoviscoelastic model presented by Zobeiry et al. (2006). To this end, modifications such as the introduction of consistent compressibility in the mass conservation equation, and a special decomposition of stresses were introduced into the formulation of the two-phase model. 5. Various numerical examples were presented to validate the use of two-phase approach to obtain the predictions for the development of stress in cured and curing composite materials. The stress model was also separately implemented in MATLAB to compare its response with the results reported by the integrated approach. ABAQUS runs were also performed separately, to verify the response of the implemented stress models. In most of the examples involving curing composite, the stress response predicted by the integrated model was compared against those reported by the stress model assuming two distinct situations; one assuming the initial distribution of volume fractions, and the other with the final distribution of volume fractions obtained from the integrated model. Using these examples, the effect of resin flow on the development of residual stresses in the composite material is analyzed and discussed.  6.2. Contributions The contributions of the current work include:   Use of a rigorous approach in the form of volume averaging technique to derive a set of governing equations for the flow-deformation behaviour of two-phase systems.    Including shear stress contributions of the fluid phase in a two-phase fluid-solid system that enables the FE model of two-phase medium to be also capable of representing the response of viscous fluids as an individual problem and in the neighbourhood of two-phase systems. Comparing the formulation and the  205  Chapter 6: Conclusions and Further Work  approach with the one presented by Chan et al. (2000), and clarifying the differences of that work with the current study.   Implementation of special measures in the FE technique in order to achieve the capability of modeling the correct response in the presence of discontinuities in the distribution of the fluid volume fraction. To model the interface between a viscous fluid and a two-phase system, extra considerations were also included.    The most important contribution of this work is the integration of modeling two very different aspects of processing composite materials, namely resin flowcompaction and stress development in a seamless model. In order to achieve that, a few steps were necessary that will be discussed here: i. Addition of compressibility of the phases compared to the state-of-the-art models of the flow in thermosetting composite materials where the individual phases are assumed incompressible. This plays a major role in the integration of modeling flow and stress effects in the processing of thermosetting composite materials as consideration of compressibility is an integral part of the stress analysis of solid materials. ii. The  mass  conservation  equation  for  two-phase  materials  with  compressible phases is modified here for the bulk behaviour of the system to be consistent with the relevant micromechanics that one decides to be appropriate for the cured system. iii. In this approach, the correct representation of the bulk behaviour of the two-phase system is guaranteed by the consistent mass conservation equation. This enables one to use the bulk modulus of the porous structure throughout the analysis of cure cycle and still model the correct response for the volumetric response of the cured two-phase material.   Using this approach, the effect of parameters important in the resin flow and deformation of the porous fibre-bed on the final shape and the amount of residual stresses and strains of the laminate may be readily studied. 206  Chapter 6: Conclusions and Further Work  6.3. Further Work The approach introduced in this thesis provides the framework for the integration of modeling resin flow and stress development. It represents a major step toward an integrated and comprehensive model that could tackle complex industrial problems. To that end, the following recommendations for further work are introduced here:  In order for the approach to be tested along with other CCA modules, the 3D form of the current integrated model should be incorporated into COMPRO 3D in ABAQUS. This would provide for the current model the capability to analyze industrial problems of complex geometry and boundary conditions.1  Adaptation of this approach to the differential viscoelastic model of is very straightforward. The implementation of the integrated model into the model presented by Zobeiry et al. (2006) should be another step toward a more comprehensive model.  Due to the lack of experimental data, two separate sets of data are used for the behaviour of resin in terms of viscosity and moduli. Devising experimental measures for viscoelastic characterization of thermosetting resins during the full range of degrees of cure is crucial to the appropriate application of the current work. 2  Using the model to simulate the processing of real composite structures in order to predict their final geometry and residual stresses. This would facilitate the identification of the capabilities of the approach and also possibly further improvements.  1  This step has in fact been successfully accomplished by Dr. Abdul R. Arafath and colleagues at Convergent Manufacturing Technologies (CMT) [Arafath (2012)]. 2  Done by Mr. Ryan Thorpe as a part of his M.Sc. thesis at the UBC Composites Group [Thorpe (2012)].  207  Chapter 6: Conclusions and Further Work   Introduction of large deformation formulation to this approach is another recommended improvement. With this addition, the model would be able to update the geometry of the composite structure as resin flow progresses and produce a more realistic representation of the deformations and the stress state in the structure.  Introduction of the capability of predicting voids formation and transport via a 3phase model including solid, fluid and a gas phase. This development would make predictions of the integrated approach more realistic as presence of voids in the final product, at least in small percentages, occurs quite often and tends to affect the properties of the composite part.  208  Bibliography  Adolf, D. and Martin, J.E., "Calculation of Stresses in Crosslinking Polymers." Journal of Composite Materials, vol. 30, pp. 13-34, 1996. Almeida, E.S. and Spilker, R.L., "Mixed and Penalty Finite Element Models for the Nonlinear Behavior of Biphasic Soft Tissues in Finite Deformation: Part I - Alternate Formulations." Computer Methods in Biomechanics and Biomedical Engineering, vol. 1, pp. 25, 1997. Almeida, E.S. and Spilker, R.L., "Mixed and Penalty Finite Element Models for the Nonlinear Behavior of Biphasic Soft Tissues in Finite Deformation: Part II - Nonlinear Examples." Computer Methods in Biomechanics and Biomedical Engineering, vol. 1, pp. 151, 1998. Antonucci, V., Cusano, A., Giordano, M., Nasser, J., and Nicolais, L., "Cure-induced residual strain build-up in a thermoset resin," Composites Part A: Applied Science and Manufacturing, vol. 37, pp. 592-601, 4, 2006. Arafath, A.R.A., "Efficient numerical techniques for predicting process-induced stresses and deformations in composite structures," PhD thesis, UBC, Vancouver, 2007. Arafath, A.R.A., Vaziri, R., and Poursartip, A., "Closed-form solution for processinduced stresses and deformation of a composite part cured on a solid tool: Part I – Flat geometries," Composites Part A: Applied Science and Manufacturing, vol. 39, pp. 11061117, 7, 2008. Arafath, A.R.A., Vaziri, R., and Poursartip, A., "Closed-form solution for processinduced stresses and deformation of a composite part cured on a solid tool: Part II – Curved geometries," Composites Part A: Applied Science and Manufacturing, vol. 40, pp. 1545-1557, 10, 2009. 209  Bibliography  Arafath, A.R.A., Convergent Manufacturing Technologies Inc., Private communications, 2012. Ateshian, G.A., Wang, X., "Sliding tractions on a deformable porous layer," Journal of Tribology, vol. 120, pp. 89, 1998. Beavers, G.S. and Joseph, D.D., "Boundary conditions at a naturally permeable wall," Journal of Fluid Mechanics, vol. 30, pp. 197, 1967. Biot, M.A., "General Theory of Three-Dimensional Consolidation," Journal of Applied Physics, vol. 12, pp. 155-164, 1941. Bogetti, T.A. and Gillespie, J.W., "Two-Dimensional Cure Simulation of Thick Thermosetting Composites," Journal of Composite Materials, vol. 25, pp. 239-273, March 01, 1991. Bogetti, T.A. and Gillespie, J.W., "Process-Induced Stress and Deformation in ThickSection Thermoset Composite Laminates," Journal of Composite Materials, vol. 26, pp. 626, 1992. Bréard, J., Henzel, Y., Trochu, F., Gauvin, R., "Analysis of dynamic flows through porous media. Part II: Deformation of a double-scale fibrous reinforcement," Polymer Composites, vol. 24, pp. 409-421, 2003. Bréard, J., Henzel, Y., Trochu, F., Gauvin, R., "Analysis of dynamic flows through porous media. Part I: Comparison between saturated and unsaturated flows in fibrous reinforcements," Polymer Composites, vol. 24, pp. 391-408, 2003. Brinkman, H.C., "Calculation of viscous force exerted by flowing fluid on dense swarm of particles," Applied Scientific Research, vol. A 1, pp. 27-34, 1949. Bruschke, M.V., Advani, S.G., "A finite element/control volume approach to mold filling in anisotropic porous media," Polymer Composites, vol. 11, pp. 398-405, 1990.  210  Bibliography  Cai, Z. and Gutowski, T., "3-D deformation behavior of a lubricated fiber bundle," Journal of Composite Materials, vol. 26, pp. 1207, 1992. Chan, B., Donzelli, P.S., Spilker, R.L., "A Mixed-Penalty Biphasic Finite Element Formulation Incorporating Viscous Fluids and Material Interfaces," Annals of Biomedical Engineering, vol. 28, pp. 589-597, 2000. Coussy, O., Poromechanics. Wiley, 2004. Davé, R., "A Unified Approach to Modeling Resin Flow During Composite Processing," Journal of Composite Materials, vol. 24, pp. 22-41, January 01, 1990. Fernlund, G, Griffith, J., Courdji, R. and Poursartip, A., "Experimental and numerical study of the effect of caul-sheets on corner thinning of composite laminates," Composites Part A: Applied Science and Manufacturing, vol. 33, pp. 411-426, 3, 2002. Fernlund, G, Rahman, N., Courdji, R., Bresslauer, M., Poursartip, A., Willden, K. and Nelson, K., "Experimental and numerical study of the effect of cure cycle, tool surface, geometry, and lay-up on the dimensional fidelity of autoclave-processed composite parts," Composites Part A: Applied Science and Manufacturing, vol. 33, pp. 341-351, 3, 2002. Fernlund, G., Osooly, A., Poursartip, A., Vaziri, R., Courdji, R., Nelson, K., George, P., Hendrickson, L. and Griffith, J., "Finite element based prediction of process-induced deformation of autoclaved composite structures using 2D process analysis and 3D structural analysis," Composite Structures, vol. 62, pp. 223-234, 11, 2003. Gebart, B.R., "Permeability of Unidirectional Reinforcements for RTM," Journal of Composite Materials, vol. 26, pp. 1100-1133, August 01, 1992. Ghaboussi, J. and Wilson E.L., "Variational Formulation of Dynamics of Fluid-Saturated Porous Elastic Solids," Journal of Engineering Mechanics, vol. 98, pp. 947, 1972. Gray, W.G. and Lee, P.C.Y., "On the theorems for local volume averaging of multiphase systems," International Journal of Multiphase Flow, vol. 3, pp. 333-340, 1977. 211  Bibliography  Gresho, P.M. and Sani, R.L., Incompressible Flow and the Finite Element Method. Wiley, 2000. Gutowski, T.G., Cai, Z., Bauer, S., Boucher, D., Kingery, J., Wineman, S., "Consolidation Experiments for Laminate Composites," Journal of Composite Materials, vol. 21, pp. 650-669, 1987. Gutowski, T.G. and Dillon, G., "The elastic deformation of fiber bundles," in Advanced Composites Manufacturing, New York: Wiley, 1997. Hashin, Z., “Theory of Fiber Reinforced Materials”, NASA-CR-1974, 1972. Hou, J.S., Holmes, M.H., Lai, W.M., and Mow, V.C., "Boundary conditions at the cartilage-synovial fluid interface for joint lubrication and theoretical verifications," Journal of Biomechanical Engineering, vol. 111, pp. 78-87, 1989. Hubert, P., "Aspects of flow and compaction of laminated composite shapes during cure," PhD thesis, UBC, Vancouver, 1996. Hubert, P. and Poursartip, A., "Review of flow and compaction modelling relevant to thermoset matrix laminate processing," J Reinf Plast Compos, vol. 17, pp. 286-318, 1998. Hubert, P., Vaziri, R., and Poursartip, A., "Two-dimensional flow model for the process simulation of complex shape composite laminates," International Journal for Numerical Methods in Engineering, vol. 44, pp. 1-26, 1999. Hubert, P. and Poursartip, A., "A method for the direct measurement of the fibre bed compaction curve of composite prepregs," Composites. Part B, Engineering, vol. 32, pp. 179, 2001. Johnston, A.A., "An integrated model of the development of process-induced deformation in autoclave processing of composite structures," PhD thesis, UBC, Vancouver, 1997.  212  Bibliography  Johnston, A., Vaziri, R. and Poursartip, A., "A Plane Strain Model for Process-Induced Deformation of Laminated Composite Structures," Journal of Composite Materials, vol. 35, pp. 1435-1469, August 01, 2001. Kim, Y.K., White, S.R., "Stress relaxation behavior of 3501-6 epoxy resin during cure," Polymer Engineering & Science, vol. 36, pp. 2852-2862, 1996. R. Larsson, "Process-modeling of composites using two-phase porous media theory," European Journal of Mechanics. A, Solids, vol. 23, pp. 15-36, 2004. Lewis, R.W. and Schrefler, B.A., The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. 2nd Edition, Wiley, 1998. Li, M. and Tucker III, C.L., "Modeling and simulation of two-dimensional consolidation for thermoset matrix composites," Composites. Part B, Engineering, vol. 33, pp. 877, 2002. Li, M., Zhu, Q., Geubelle, P.H., and Tucker III, C.L., "Optimal curing for thermoset matrix composites: Thermochemical considerations," Polymer Composites, vol. 22, pp. 118-131, 2001. Melo, J.D.D., Radford, D.W., "Time and temperature dependence of the viscoelastic properties of PEEK/IM7," Journal of Composite Materials, vol. 38, pp. 1815-1830, 2004. Mow, V.C., Kuei, S.C., Lai, W.M., and Armstrong, C.G., "Biphasic creep and stress relaxation of articular cartilage in compression? Theory and experiments," J. Biomech. Eng., vol. 102, pp. 73-84, Feb, 1980. O'Brien, D.J., Mather, P.T., and White S.R., "Viscoelastic Properties of an Epoxy Resin during Cure," Journal of Composite Materials, vol. 35, pp. 883, 2001. Pillai, K.M. and Advani, S.G., "A Model for Unsaturated Flow in Woven Fiber Preforms during Mold Filling in Resin Transfer Molding," Journal of Composite Materials, vol. 32, pp. 1753-1783, October 01, 1998.  213  Bibliography  Pillai, K.M., Tucker III, C.L. and Phelan Jr., F.R., "Numerical simulation of injection/compression liquid composite molding. Part 1: Mesh generation," Composites Part A: Applied Science and Manufacturing, vol. 31, pp. 87-94, 1, 2000. Pillai, K.M., Tucker III, C.L. and Phelan Jr., F.R., "Numerical simulation of injection/compression liquid composite molding. Part 2: Preform compression," Composites Part A: Applied Science and Manufacturing, vol. 32, pp. 207-220, 2, 2001. Pillai, K.M., "Governing equations for unsaturated flow through woven fiber mats. Part 1. Isothermal flows," Composites. Part B, Engineering, vol. 33, pp. 1007-1019, 2002. Prasatya, "A Viscoelastic Model for Predicting Isotropic Residual Stresses in Thermosetting Materials: Effects of Processing Parameters," Journal of Composite Materials, vol. 35, pp. 826, 2001. Prevost, J.H., "Nonlinear transient phenomena in saturated porous media," Computer Methods in Applied Mechanics and Engineering, vol. 20, pp. 3, 1982. Prevost, J.H., "Implicit-explicit schemes for nonlinear consolidation," Computer Methods in Applied Mechanics and Engineering, vol. 39, pp. 225, 1983. Prevost, J.H., "Wave propagation in fluid-saturated porous media: an efficient finite element procedure, " International Journal of Soil Dynamics and Earthquake Engineering, vol. 4, pp. 183, 1985. Sandhu, R.S. and Wilson, E.L., "Finite-element analysis of seepage in elastic media," American Society of Civil Engineers, Journal of the Engineering Mechanics Division, vol. 95, pp. 641-652, 1969. Simon, B.R., Wu, J.S.S., and Zienkievicz, O.C., "Evaluation of u-w and u- Pi finite element methods for the dynamic response of saturated porous media using onedimensional models." International Journal for Numerical and Analytical Methods in Geomechanics, vol. 10, pp. 461, 1986a.  214  Bibliography  Simon, B.R., Wu, J.S.S., and Zienkievicz, O.C., "Evaluation of higher order, mixed and Hermitean finite element procedures for dynamic analysis of saturated prous media using one-dimensional models." International Journal for Numerical and Analytical Methods in Geomechanics, vol. 10, pp. 483, 1986b. Simon, S.L., McKenna, G.B., and Sindt, O., "Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation," Journal of Applied Polymer Science, vol. 76, pp. 495, 2000. Smith, G.D. and Poursartip, A., "A Comparison of Two Resin Flow Models for Laminate Processing," Journal of Composite Materials, vol. 27, pp. 1695, 1993. Spilker, R.L. and Maxian, T.A., "A mixed-penalty finite element formulation of the linear biphasic theory for soft tissues," Int J Numer Methods Eng, vol. 30, pp. 1063-1082, 1990. Tan, H. and Pillai, K.M., "Multiscale modeling of unsaturated flow in dual-scale fiber preforms of liquid composite molding I: Isothermal flows," Composites Part A: Applied Science and Manufacturing, vol. 43, pp. 1-13, 1, 2012. Tan, H. and Pillai, K.M., "Multiscale modeling of unsaturated flow of dual-scale fiber preform in liquid composite molding II: Non-isothermal flows," Composites Part A: Applied Science and Manufacturing, vol. 43, pp. 14-28, 1, 2012. Taylor, G.I., "A model for the boundary condition of a porous material. Part 1," Journal of Fluid Mechanics, vol. 49, pp. 319, 1971. Terzaghi, K., Theoretical Soil Mechanics. Wiley, New York, 1943. Thorpe, R., "Experimantal characterization of the viscoelastic behaviour of a curing epoxy matrix composite from pre-gelation to full cure" M.A.Sc. thesis, UBC, Vancouver, 2012. Trochu, F., Gauvin, R., Gao, D.M., "Numerical analysis of the resin transfer molding process by the finite element method," Advances in Polymer Technology, vol. 12, pp. 329-342, 1993. 215  Bibliography  Trochu, F., Ruiz, E., Achim, V. and Soukane, S., "Advanced numerical simulation of liquid composite molding for process analysis and optimization," Composites Part A: Applied Science and Manufacturing, vol. 37, pp. 890-902, 6, 2006. Tucker, C.L. and Dessenberger, R.B., "Governing equations for flow and heat transfer in stationary fiber beds," in Flow and Rheology in Polymer Composites Manufacturing, S. G. Advani, Ed. Amsterdam: Elsevier, 1994, pp. 257-323. Twigg, G., Poursartip, A. and Fernlund, G., "Tool–part interaction in composites processing. Part I: experimental investigation and analytical model," Composites Part A: Applied Science and Manufacturing, vol. 35, pp. 121-133, 1, 2004. Twigg, G., Poursartip, A. and Fernlund, G., "Tool–part interaction in composites processing. Part II: numerical modeling," Composites Part A: Applied Science and Manufacturing, vol. 35, pp. 135-141, 1, 2004. Vermeer, P.A. and Verruijt, A., "Accuracy condition for consolidation by finite elements," Int. J. Numer. Anal. Methods Geomech., vol. 5, pp. 1-14, 1981. Whitaker, S., The Method of Volume Averaging. Springer, 1998. White, S.R. and Hahn, H.T., "Process Modeling of Composite Materials: Residual Stress Development during Cure. Part I. Model Formulation," Journal of Composite Materials, vol. 26, pp. 2402-2422, January 01, 1992. White, S.R. and Kim, Y.K., "Process-induced residual stress analysis of AS4/3501-6 composite material," Mechanics of Advanced Materials and Structures, vol. 5, pp. 153, 1998. Zhu, Q., Geubelle, P.H., Li, M. and Tucker III, C.L., "Dimensional accuracy of thermoset composites: Simulation of process-induced residual stresses," J. Composite Mater., vol. 35, pp. 2171-2205, 2001. Zobeiry, N., "Viscoelastic constitutive models for evaluation of residual stresses in thermoset composites during cure," PhD thesis, UBC, Vancouver, 2006. 216  Bibliography  Zobeiry, N., Vaziri, R., and Poursartip, A., "Differential Implementation of the Viscoelastic Response of a Curing Thermoset Matrix for Composites Processing," Journal of Engineering Materials and Technology, vol. 128, pp. 90, 2006. Zobeiry, N., Vaziri, R. and Poursartip, A., "Computationally efficient pseudo-viscoelastic models for evaluation of residual stresses in thermoset polymer composites during cure," Composites Part A: Applied Science and Manufacturing, vol. 41, pp. 247-256, 2, 2010.  217  

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

Comment

Related Items