@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Civil Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Haghshenas, Seyed Mehdi"@en ; dcterms:issued "2013-01-02T15:42:23Z"@en, "2012"@en ; vivo:relatedDegree "Doctor of Philosophy - PhD"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms:description """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 two-phase 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 pseudo-viscoelastic 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."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/43758?expand=metadata"@en ; skos:note " 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 ii 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 two- phase 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 pseudo- viscoelastic 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. iii 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 Table of Contents iv 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 Table of Contents v 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 Table of Contents vi 6.3. Further Work ..................................................................................................................... 207 Bibliography .................................................................................................................. 209 vii 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 viii 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 List of Figures ix 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 well- pronounced ................................................................................................................ 75 List of Figures x 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 m 2 ..................................................................................................... 78 Figure 3.26: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10 -8 m 2 ..................................................................................................... 78 Figure 3.27: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10 -10 m 2 .................................................................................................... 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 List of Figures xi 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 List of Figures xii 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 List of Figures xiii 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 List of Figures xiv 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 List of Figures xv 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 List of Figures xvi 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 List of Figures xvii 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 xviii List of Symbols Ta 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 id f Vector of drag force G12 In-plane shear modulus G23 Transverse shear modulus Gc Shear modulus of the composite material List of Symbols xix Gfb Shear modulus of the fibre-bed Gm Shear modulus of resin G r Relaxed shear modulus of resin G u 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 K r Relaxed bulk modulus of resin K u Un-relaxed bulk modulus of resin m Cool-down rate in a cure cycle iN Displacement shape function for the i th node p iN Pressure shape function for the i th pressure node  iN Shape function of fluid volume fraction for the i th Gauss point P Resin pressure p Vector of pressure DOFs te Variable time, used in the pseudo-viscoelastic model T Temperature iu Displacement vector of the system u Vector of displacement DOFs List of Symbols xx iv Vector of resin relative velocity v Vector of relative velocity DOFs pw 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 sv  Volumetric strain of the system sv  Rate of volumetric strain of the system mech sv  Rate of mechanical volumetric strain of the system free sv  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 List of Symbols xxi ζ  , ij  Biot effective stress tensor  Vector of stress rates ijs  Total stress tensor of the system ijps  „s-p‟ effective stress tensor  Dummy time variable (for time integration) ij Shear stress tensor ijm  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 xxii 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. 1 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. Chapter 1: Introduction 2 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, thermo- chemical changes and curing of resin, and development of residual stresses in the composite material. This complexity has led many researchers to use an „integrated sub- model‟ 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 Chapter 1: Introduction 3 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 Chapter 1: Introduction 4 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 Chapter 1: Introduction 5 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 fibre- reinforced 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 fibre- reinforced composites. Liquid composite moulding (LCM) processes such as resin Chapter 1: Introduction 6 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 Chapter 1: Introduction 7 composite matrix resin, and micromechanics models were used to determine composite ply mechanical properties and behaviour, including thermal expansion and cure- shrinkage. 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. Chapter 1: Introduction 8 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 pseudo- viscoelastic 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. Chapter 1: Introduction 9 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 two- phase media and an alternative FE approach in u-v-p formulation. The additional Chapter 1: Introduction 10 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. Chapter 1: Introduction 11 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. Chapter 1: Introduction 12 1.6. Figures (a) (b) Figure 1.1: Autoclave processing steps, (a) lay-up, (b) cure with a controller loop (adapted from Hubert, 1996) Tool Vacuum plug Breather BleederVacuum bag LaminateDam Pressure source Heater Autoclave Controller output Controller Vac. pump Laminate-tool in vacuum bag Sensors output (Pb) () (T) (T) () (Pb) Time Cure cycle Chapter 1: Introduction 13 Figure 1.2: Integrated sub-model approach in process modeling (adapted from Johnston et al., 2001) Figure 1.3: Different regimes of resin behaviour during a typical autoclave processing history State Variables Database Output module Stress module Autoclave controller module Thermo-chemical module Flow module Input module 0 20 40 60 80 100 120 140 160 180 0 50 100 150 200 250 300 350 Time (min) T em p er a tu re (° C ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 D eg r ee o f cu re Autoclave Temp. Part Temp. Degree of cure Flow regime Stress development regime T em p er a tu re (° C ) D eg r ee o f cu re Chapter 1: Introduction 14 Figure 1.4: Time-history of changes in resin viscosity during a typical autoclave processing history Figure 1.5: Schematic representation of resin flow through the fibre-bed 0.1 1 10 100 1000 0 20 40 60 80 100 120 140 160 180 0 50 100 150 200 250 300 350 R e s in v is c o s it y ( P a .s ) T e m p e ra tu re ( °C ) time (min) Autoclave Temp. Part Temp. Resin viscosity 1 2, 3 Fibre-bed Compaction 15 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)]. Chapter 2: Governing Equations and FE Model of Two-phase Media 16 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. Chapter 2: Governing Equations and FE Model of Two-phase Media 17 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 mixed- penalty 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 Chapter 2: Governing Equations and FE Model of Two-phase Media 18 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 Darcy- Brinkman 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-v- p 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     phase- in the liet doesn' if0 phase- in the lies if1 )(    x x xX (2.1) Chapter 2: Governing Equations and FE Model of Two-phase Media 19 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  V BdV V B 1 (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   VV dVXB V dVB V B   11 (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     V V V dVX dVXB dVB V B        1 (2.4) Volume fraction of the α-phase is defined as dVX VV V V     1 (2.5) The phase and intrinsic phase averages for an arbitrary α-phase are related by    BB (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)]    S dSB V BB n 1 (2.7) Chapter 2: Governing Equations and FE Model of Two-phase Media 20 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)] . 1     S dS V nAAA (2.8) 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 0 ,0   f m v v (2.9) 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 0 1   mfS mfmm dS V nvv (2.10) where mfS is the interfacial surface between the fluid and solid phase, and mfn is the unit vector normal to mfS 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-hand- side of the above equation is equal to the rate of resin volume fraction, i.e. t dS V m S mfm mf     nv 1 (2.11) Chapter 2: Governing Equations and FE Model of Two-phase Media 21 Therefore, the macroscopic form of the mass conservation equation of the fluid phase may be written as 0   m m t v  (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 0   f f t v  (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 1 mf  (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 0 mf vv (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 0ζ 0ζ   f m , (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 0nζζ   mfS mfmm dS V 1 (2.17) The total stress tensor of the fluid phase may be decomposed into a hydrostatic pressure term Pm and a deviatoric stress tensor ηm Chapter 2: Governing Equations and FE Model of Two-phase Media 22 Iηζ mmm P (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 m m mmmmmm PPP   ηηζ (2.19) Substituting (2.19) into (2.17) leads to 0nζη   mfS mfmm m mm m mm dS V PP 1  (2.20) If the vector of drag forces between the fluid phase and the solid phase is defined as   mfS mfmm m md dS V P nζf 1  (2.21) the equilibrium equation of the fluid phase takes the final form of 0fη  dm m mm P (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 0nζζ   mfS fmff dS V 1 (2.23) If we assume that the surface tension between the fluid and solid particles may be neglected, then we may write [Tucker and Dessenberger (1994)] mfmfmfmf Son 0nζnζ  (2.24) which leads to   mfmf S mfm S fmf dS V dS V nζnζ 11 (2.25) Considering the definition of drag force in (2.21, the above equation may be rewritten as Chapter 2: Governing Equations and FE Model of Two-phase Media 23 m m md S fmf PdS V mf  fnζ 1 (2.26) Substituting (2.26) into (2.23) leads to the following form of the equilibrium equation of the solid phase 0fζ  m m mdf P  (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 ff η Iηζ mm f f f f P  (2.28) where ff η represents a stress component dependent on the deformation of the fibre- bed 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 Iηζ mmf f fff P   (2.29) Substituting (2.14) in (2.29), we arrive at  Iηζ m m mff P  1 (2.30) which leads to      m m mm m mf m m mff PP P     1 1 η ηζ (2.31) Substituting (2.31) into (2.27) leads to   0fη  m m mdf P 1 (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 two- phase system Chapter 2: Governing Equations and FE Model of Two-phase Media 24         0fη 0fη vv m m mdf dm m mm mf P P   1 0 (2.33) The constitutive behaviour of the fibre-bed could be expressed as   Tfffff  uuCη 2 1 (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  Tmmm )( vvη   (2.35) where μ is the viscosity of the fluid phase. The volume-averaged form of the above equation is  Tmmm  vvη  (2.36) in which, for  mv , we may rewrite the volume averaging theorem in (2.7) to have dS V mfS mfmmm  nvvv 1 (2.37) It is commonly assumed that the slip velocity at the interface is zero [Tucker and Dessenberger (1994)], and therefore mffm Son vv  (2.38) Thus, the integral term on the right-hand-side of (2.37), may be rewritten as dS V mfS mffmm  nvvv 1 (2.39) Arguing that the changes in the solid phase velocity fv are small in comparison to the changes in fluid phase velocity mv , we assume that fv on the interface can be Chapter 2: Governing Equations and FE Model of Two-phase Media 25 approximated with f f v . Incorporating this assumption into the integral term in (2.39), we have           dSV dS V dS V mfmfmf S mf f f S mf f f S mfm nvnvnv 111 (2.40) Applying the averaging theorem in (2.7) to the phase function Xm leads to  mfS mfmmm dSX V XX n 1 (2.41) Within the fluid phase, Xm is equal to unity and mX equals zero. Also, it is readily deduced from (2.3) that  mX is simply the fluid volume fraction m . Therefore, (2.41) may be simplified to [Tucker and Dessenberger (1994)] dS V mfS mfm  n 1  (2.42) Substituting (2.42) into (2.40) results in m f f S mfm dS V mf  vnv 1 (2.43) which in turn, if substituted in (2.39) leads to m f fmm  vvv (2.44) Substituting the above into (2.36) and neglecting the approximations, we get   ffmmffTmmm  vvvvη  (2.45) Rewriting the above equation leads to     Tmffmmffmm   vvvvη (2.46) 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)] Chapter 2: Governing Equations and FE Model of Two-phase Media 26  ffmmmd   vvSf 12 (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         0vvSuuC   mmmffmmmTffff P  1 2 1 12 (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         0vvS vvvv    f f m mm T m f fmm f fm m mm P 12   (2.49) Thus, the set of equations (2.33) may be expressed in the form                              0vvSuuC 0vvS vvvv vv m m m f f m mm Tf f f f f f m mm T m f fmm f fm m mm mf P P    1 2 1 - 0 12 12 (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 ffc  vv (2.51) and the volume-averaged relative velocity (or flow velocity) of the resin as  ffmmmflow  vvv  (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 Chapter 2: Governing Equations and FE Model of Two-phase Media 27   cmflowm mcmflow m m cff c f f vvv vvv vv vv        (2.53) Using the above equations, one may easily arrive at the governing equations with respect to the new velocity variables                        0vSuuC 0vSvvvv vv m mmflowm T cc flowm T cmflowcmflow m mm flowc P P   1 2 1 - 0 1 1 (2.54) where cc uv  (2.55) and uc is the composite displacement field. For the sake of convenience, variables uc, vflow, m mP  , and m are replaced with u, v, P, and  respectively, i.e.      m m m flow c PP vv uu (2.56) Using index notations, the set of governing equations (2.54) may be written in the form 0,,  iiii vu (2.57)       01 ,,,,,,   jijjijijjijii vSuvuvP   (2.58)      01 2 1 , 1 ,,,   ijijjkllkijkl PvSuuC  (2.59) Adding up the two equilibrium equations of (2.58) and (2.59), leads to the total equilibrium equation of the domain Chapter 2: Governing Equations and FE Model of Two-phase Media 28          0 2 1 ,,,,,,,,,  ijijijjijijkllkijkl PuvuvuuC   (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                          0 2 1 0 0 ,,,,,,,,, 1 ,,,,,, ,, ijijijjijijkllkijkl jijjijijjijii iiii PuvuvuuC vSuvuvP vu      (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.              0 2 1 0 0 ,,,, 1 , ,, ijkllkijkl jiji iiii PuuC vSP vu   (2.62) Chapter 2: Governing Equations and FE Model of Two-phase Media 29 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 2 nd 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. Chapter 2: Governing Equations and FE Model of Two-phase Media 30 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 0)( ,,   dwvu piiii (2.63) 0)()( )( 1 ,,, ,,,,,            dwnPdwvSdwuu dwvvdPwdPw pjijijpjijjpijji jpijjipiip  (2.64) 0)()( )()( ,,,, ,,,,,,2 1          dwnPdPwdwuu dwvvdwuuC pjijijijipjpijji jpijjijpkllkijkl  (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 )( 2 1 ,, kllkijklij uuC  (2.66) and the components of shear stress tensor in the resin are expressed as     ijijjijiij uvuv ,,,,    (2.67) The term jijijij nP )(   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 jijij nP )(   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. Chapter 2: Governing Equations and FE Model of Two-phase Media 31 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 139                                                       0 f f p v u CCC CCC CCC p v u KKK KKK KKK v u pppvpu vpvvvu upuvuu pppvpu vpvvvu upuvuu t (2.68) where 118u , 118v , and 13p 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   jijijijni nPt  )( (2.69) we may rewrite (2.65) to have          dwtdPwdwuu dwvvdwuuC p n iipjpijji jpijjijpkllkijkl )( ,,,, ,,,,,,2 1 )( )()(  (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) 1818   dT T BDBKuu (2.71) 1818 1    d T BBKuv  (2.72) 318   d pTδNBKup (2.73) 1818 1    d T BBCuu  (2.74) 1818 0Cuv (2.75) Chapter 2: Governing Equations and FE Model of Two-phase Media 32 318 0Cup (2.76) 118 )(     dnTd tNfu (2.77) Where t is the total traction vector on the boundary of the medium with components in x- and y-directions in the form of          )( )( )( n y n xn t t t (2.78) and Nd is defined as        N0 0N Nd (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  921 ... NNNN (2.80) and assuming that the element has three internal nodes for pressure interpolation, we may introduce the matrix of pressure shape functions as  pppp NNN 321N (2.81) The displacement, resin relative velocity, and the pressure fields are interpolated on the domain of the element using the following equations Chapter 2: Governing Equations and FE Model of Two-phase Media 33 i i i iy i iyy ix i ixx iy i iyy ix i ixx pNp vNv vNv uNu uNu                9 1 9 1 9 1 9 1 9 1 pN vN vN uN uN (2.82) The displacement shape functions may be written as follows                           229 2 8 2 7 2 6 2 5 4 3 2 1 11 11 2 1 11 2 1 11 2 1 11 2 1 11 4 1 11 4 1 11 4 1 11 4 1                   N N N N N N N N N (2.83) The definition presented in (2.79) is relevant if the vector of main variables is defined as 139                 p v v u u y x y x (2.84) Chapter 2: Governing Equations and FE Model of Two-phase Media 34 while we are interested to present the vectors of displacement and relative resin velocity degrees of freedom in the form                                               y x y x y x y x y x y x v v v v v v u u u u u u 9 9 2 2 1 1 9 9 2 2 1 1 ,  vu (2.85) which leads to rewriting (2.79) so that        921 921 000 000 NNN NNN d   N (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 183,9,9,2,2,1,1 ,9,2,1 ,9,2,1 000 000             xyxyxy yyy xxx NNNNNN NNN NNN    B (2.87) B1 is another matrix containing the spatial derivatives of the shape functions defined in (2.83), and is defined as 183,9,9,2,2,1,1 ,9,2,1 ,9,2,1 1 202020 020202             xyxyxy yyy xxx NNNNNN NNN NNN    B (2.88) δ is defined by  T011δ (2.89) Chapter 2: Governing Equations and FE Model of Two-phase Media 35 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                        yx yx J (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 1 9 1 y , i ii i ii yNxNx (2.91) which upon substitution in (2.90) yields        2221 1211 JJ JJ J (2.92) where . , , , 9 1 ,12 9 1 ,21 9 1 ,12 9 1 ,11       i ii i ii i ii i ii yNJxNJ yNJxNJ   (2.93) Using the inverse of Jacobian matrix, we may define a matrix consisting of the spatial derivative of the shape functions as                 , ,1 , , N N J N N B y x (2.94) where the inverse of Jacobian can be easily written as Chapter 2: Governing Equations and FE Model of Two-phase Media 36           1121 1222 21122211 1 1 JJ JJ JJJJ J (2.95) B facilitates the definition of B based on the components of B as 183192912221121 292221 191211 000 000             BBBBBB BBB BBB    B (2.96) To determine the pressure shape functions N p, let‟s assume that the coordinates of the three pressure nodes are      333222111 , ,, ,,  PPP (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 0 23 23 22 23 23                      (2.98) We may set the first shape function for pressure to be                        23 23 22 23 23 1      aN p (2.99) and the parameter a may be calculated by setting   1, 111 N . Substituting the calculated value of a into (2.99) leads to                               23 23 22 23 23 133132232112 23 1       pN (2.100) The other two pressure shape functions of the 9-3 element can be derived in a similar approach Chapter 2: Governing Equations and FE Model of Two-phase Media 37                               31 31 33 31 31 133132232112 31 2       pN (2.101)                               12 12 11 12 12 133132232112 12 3       pN (2.102) We choose the location of the pressure nodes to be      31,0 ,31,31 ,31,31 321 PPP  which leads to the simplification of the equations for the pressure shape functions as follows        3 1 2 4 3 1  pN (2.103)        3 1 2 4 3 2  pN (2.104)        3 1 2 3 3  pN (2.105) 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                     22100 01 01 211     E TD (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   jijiji n m nPt   )( (2.107) Chapter 2: Governing Equations and FE Model of Two-phase Media 38 we may rewrite (2.64) to have            dwtdwuudPw dPwdwvSdwvv pi n mjpijjipi ippjijjpijji )( ,,,, , 1 ,,, )( )(  (2.108) 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) 1818 0Kvu (2.109)   1818 1 1 1 1        ddd d T d T d T d T NSNBBNSNBBK vv  (2.110)   318 22 1 22 1    ddd pTTpTpT δNBBδNBδNBK vp  (2.111) 1818 1    d T BBCvu  (2.112) 1818 0Cvv (2.113) 318 0Cvp (2.114) 118 )(     dnm T d tNfv (2.115) Where B2 is a matrix consisting of the displacement shape functions and the spatial derivatives and the resin volume fraction  in the form of 1839,9,2,2,1,1, 9,2,1, 9,2,1, 2 202020 020202             NNNNNN NNN NNN xyxyxy yyy xxx       B (2.116) Chapter 2: Governing Equations and FE Model of Two-phase Media 39 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 mc mm c m VV VV V V    0 0 0 0 0 ,  (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 0c m v V V  (2.118) leads to the following equation for the volume fraction of resin v v       1 0 (2.119) Substituting the relationship between the volumetric strain and the total displacement field results in nn mm u u , ,0 1     (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, 515 ) as Chapter 2: Governing Equations and FE Model of Two-phase Media 40                           229 2 8 2 7 2 6 2 5 4 3 2 1 5353 9 25 53515 18 25 51553 18 25 53515 18 25 51553 18 25 515515 36 25 515515 36 25 515515 36 25 515515 36 25                            N N N N N N N N N (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 i i iN      9 1 θN (2.122) where θ is the vector of the values of fluid volume fraction at Gauss points. Spatial gradients of  may be obtained by θNθN yyxx ,,,, ,    (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: 0,,    dwudwv piipii  (2.124) which then leads to defining the components on the last row of the matrices in (2.68) 183 0CK pvpu (2.125) Chapter 2: Governing Equations and FE Model of Two-phase Media 41 183   d TTp BδNCK pupv (2.126) 33 0CK pppp (2.127) Hence, the matrix equation (2.68) may be written in the following simplified form 139                                                       0 f f p v u 00C 00C 00C p v u 0K0 KK0 KKK v u pu vu uu pv vpvv upuvuu t (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  4321 NNNNN (2.129) where the kinematic shape functions are defined as                    11 4 1 11 4 1 11 4 1 11 4 1 4 3 2 1 N N N N (2.130) Chapter 2: Governing Equations and FE Model of Two-phase Media 42 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. 1pN (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. 117                 p v v u u y x y x (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 ( 33 ) as                        3333 4 3 3333 4 3 3333 4 3 3333 4 3 4 3 2 1 N N N N (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, Chapter 2: Governing Equations and FE Model of Two-phase Media 43 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: . )( appl n i ft  .defii vv   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: . )( appl n i ft  mappli n m ft . )(  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 .. applmappl ff   .  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: .defii uu  .defii vv  Chapter 2: Governing Equations and FE Model of Two-phase Media 44  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: .defii uu  mappli n m ft . )(  Since in this case, the only plausible applied traction is due to ambient air pressure on the boundary, we may write .. applmappl ff   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 FXCKX   (2.134) where            p v u X (2.135) 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   tnnn   XXX 1 (2.136)   11   nnn XXX  (2.137) Chapter 2: Governing Equations and FE Model of Two-phase Media 45 where 10  . A few well-known special cases of the above scheme are Forward Euler ( 0 ), Backward Euler ( 1 ), and Crank-Nicholson scheme ( 5.0 ). Writing equation (2.134) at time tn+θ leads to    nnnnn FXCXK  (2.138) Substituting (2.136) and (2.137) into the above equation, we have          nnnnnnn t FXXCXXK 111 (2.139) which after some manipulations, leads to a system of equations for the vector of main variables at the current time tn+1          nnnnn ttt FXKCXKC 11 (2.140) or in an expanded form as                                                                                 nnn nn t t tt ttt t tt ttt 0 f f p v u 0KC KKC KKKC p v u 0KC KKC KKKC v u pvpu vpvvvu upuvuuuu pvpu vpvvvu upuvuuuu 1 11 111 1 (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. Chapter 2: Governing Equations and FE Model of Two-phase Media 46     nn uKuu (2.142) is replaced by the step-wise description        nnnnn uuKff uuuu    intint (2.143)   nKu f is an accumulated load vector at step n from previous time-steps obtained by the following recursive relationship        nnnnn uuKff uuuu   11int1int (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             uuuupuvuuu fuCpKvKuuKf   11111111int nnnnnnnnnn  (2.145) At every time step, an iterative solution is performed. At the k th iteration, the matrix equation of the system is                                                          k nn k nnn nn k n k n k n k n k n t tt ttt 1 11 int1 1 1 1 1 1 1 1 1 uuC uuCf ff p v uu 0KC KKC KKKC pu vuv uu pvpu vpvvvu upuvuuuu (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. Chapter 2: Governing Equations and FE Model of Two-phase Media 47 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         dPwdPwdwPn dwPdwPndwP piippi ippipi ,, ,,     (2.147) We should remind the reader that the assumption for performing integration by parts in the general form of     b a b a b a vduuvudv (2.148) 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       Rd dL xxxxP xxxxP xP )( )( )( 2 1 (2.149) Chapter 2: Governing Equations and FE Model of Two-phase Media 48 with a jump equal to α at the interface of the two elements. Differentiating the pressure function along the length of the domain leads to          Rd dd dL xxxxP xxxx xxxxP xP )( )( )( )( 2 1  (2.150) At this point, it is helpful to review the following formula regarding integration of Dirac delta function )()()( dd xfdxxxxf     (2.151) Therefore, the actual value of the above-mentioned integral term obtained without integration by parts will take the form     d d R d d L R L x x d x x x x x x dxxxxdxxPxdxxPxdxPI )()()()()()( 22111  (2.152) Taking advantage of (2.151), and considering the symmetry of δ function, we may conclude that  )()( 2 )()()()( 2122111 dd x x x x xxdxxPxdxxPxI R d d L      (2.153) We now calculate the value of the same integral term that the current FE approach with the integration by parts would result     R d d L R L R L x x x x x x x x PdxPdxPdxPI 2 (2.154) Expanding the above relationship leads to     R d d L R L x x x x x x dxxPxdxxPxPI )()()()( 22112  (2.155) and applying integration by parts on the integral terms in the above, we have Chapter 2: Governing Equations and FE Model of Two-phase Media 49         R d d L R d d L R L x x x x x x x x x x dxxPxdxxPxPPPI )()()()( 221122112  (2.156) After simplifications, we arrive at )()()()()()()()( 221122112 dddd x x x x xPxxPxdxxPxdxxPxI R d d L    (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 )()( 12 dd xPxP  and after simplifications, we arrive at the correction term for the 1D case       dxx avgdddd PxxxPxPII    )()( )()( 2 1 212121 (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      d dwnP piavg  (2.159) 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 Chapter 2: Governing Equations and FE Model of Two-phase Media 50 multiplied by pressure degrees of freedom, they contribute to the vpK components of the stiffness matrix. Assuming a general interface between the elements m1 and m2, we may rewrite (2.159) in the form          dd dwnPPdwnPP pimmmmpimm 212121 2 1 2 1  (2.160) The corrections needed in vpK (or additional vpK matrices) for the two involved elements may then be defined as               318 318 2221 2 1121 1       d d d d m p m Tn mmmextra m p m Tn mmmextra NNK NNK vp vp   (2.161) The above integrals are estimated using 3 Gauss points along the specified boundary located at 515,0or  . nN represents the product of the shape functions by the normal vector to the boundary of each element, and is written as  992211 NnNnNnNnNnNn yxyxyxn N (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. Chapter 2: Governing Equations and FE Model of Two-phase Media 51 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)]:     2 2 1 1     ijij  (2.163) Assuming that 1 represent the viscous fluid and therefore 11  , one may remove the subscript from 2 and rewrite (2.163) to have            2,,221 1111 ijjiijijij vv   (2.164) To apply the above condition to the FE approach, a virtual traction force equal to     jijji tvv 2,,11   (2.165) is applied to the fluid phase at the boundary between the two media. jt 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 vvK matrix for element m1 as follows         1818 13 1111 11    d d mm T mmextra BBK vv  (2.166) where 3B is defined by Chapter 2: Governing Equations and FE Model of Two-phase Media 52 183992211 921 921 3 000 000             NnNnNnNnNnNn NnNnNn NnNnNn xyxyxy yyy xxx    B (2.167) Chapter 2: Governing Equations and FE Model of Two-phase Media 53 2.7. Figures Figure 2.1: Schematic representation of Q2P-1 (9-3) element Figure 2.2: Schematic representation of Q1P0 (4-1) element   x y Pressure nodes u / v nodes x y Pressure node u / v nodes   Chapter 2: Governing Equations and FE Model of Two-phase Media 54 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 x P,φ α β P φ P1(x) P2(x) φ1(x) φ2(x) xL xd xR 55 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 Chapter 3: Verification of the Two-phase Model 56 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×10 3 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×10 3 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                                t h ESj h yj j f tyP y j j 2 22 1 1 0 4 12 exp 2 12 cos 12 14 ,    (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 Chapter 3: Verification of the Two-phase Model 57 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 (2 nd equation in (2.62)) into the mass conservation(1 st equation in (2.62)); leading to the following set of two equations:           0 2 1 0 ,,,, ,,, ijkllkijkl ijijii PuuC PSu  (3.2) As a result, the pressure undergoes a second-order differentiation in space, and inter- element continuity of pressure becomes necessary. This is not the case in the present u-v- p formulation where all three of the mass conservation equation, fluid equilibrium equation, , and the total equilibrium equation (3 rd 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 Chapter 3: Verification of the Two-phase Model 58 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, Chapter 3: Verification of the Two-phase Model 59 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             0 0 0 ,0 22 2 12 2 zhv Sz v hz z v    (3.3) 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)          00 1 zz z v z v vv  (3.4) 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                                                                                       0 coshsinh sinh 0 , coshsinh cosh 1 2 212 20 11 212 2 0 zh S h S h S h hz S V hzhz S h S h S h S h S V v        (3.5) Chapter 3: Verification of the Two-phase Model 60 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                   0 0 , 22 2 12 2 zh x P v Sz v hz x P z v    (3.6) 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                                                                                                                                            0 1sinh coshsinh 2 1 0 , coshsinh cosh12sinh 2 1 22 212 2 1 1 1 212 221 1 2 22 zhehz S S h S h S h S h e S he x PS hz S h S h S h S hS S hh zhz x P v hz S S h S h                (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 Beavers- Joseph problem. We present numerical results for both of the above problems. Chapter 3: Verification of the Two-phase Model 61 It is assumed that mmhh 421  , 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         01 ,,,,,,   jijjijijjijii vSuvuvP   (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 Chapter 3: Verification of the Two-phase Model 62 (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 Chapter 3: Verification of the Two-phase Model 63 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 Chapter 3: Verification of the Two-phase Model 64 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×10 7 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. Chapter 3: Verification of the Two-phase Model 65 3.3. Tables Table 3.1: Resin and fibre-bed properties for the AS4/3501-6 angle laminate Table 3.2: Resin and fibre-bed properties for the AS4/8552 angle laminate Fibre-bed elastic properties Resin viscosity and fibre-bed permeability Resin degree of cure     exp/exp RTU 8.14 /114477 .106.4 17       molJU sPa             0.3 1 0.3 1 3 21     C dt d BCC dt d RTE II IeAC  47.0 56600 77800 80700 102667.3 103567.3 105017.3 3 2 1 13 3 17 2 17 1           B molJE molJE molJE sA sA sA              MPaSymm E GPa D 5.0. 0 00100 33  1fV   2 32 1 1 4 f ff V V k r S      1 1 4 3 2 2     fa faf VV VV k r S 2.0k mrf 4 81.0aV               03 2.003.41 16.02.015.15 1.016.033.3 01.067.0 3 3 3 3 3 3      MPa MPa MPa MPa MPa E Resin degree of cure Resin viscosity Fibre-bed elastic propertiesFibre-bed permeability              MPaSymm E GPa D 5.0. 0 00100 33  1fV   2 32 1 1 4 f ff V V k r S      1 1 4 3 2 2     fa faf VV VV k r S 2.0k mrf 4 68.0aV     TC nmRTE CTC a e Ae dt d        01 1 C C n m sA molJE CT C a          /19.0 1048.5 1.43 74.2 813.0 1053.1 66500 3 0 15   47.0 5.2 8.3 /76536 .1045.3 10       g B A molJE sPaA                  06.6 1.01.62 079.01.022 035.0079.065.8 0035.008.2 3 3 3 3 3 3      MPa MPa MPa MPa MPa E        BA g g RTEA            exp Chapter 3: Verification of the Two-phase Model 66 3.4. Figures 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 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 b h f0 x y h = 25 mm b = 2.5 mm υ0 = 0.5 E = 70 MPa ν = 0 f0 = 700 kPa Sy = 4.5×10 -14 m2 b f0 2h 2h x y h = 25 mm b = 2.5 mm υ0 = 0.5 E = 70 MPa ν = 0 f0 = 700 kPa Sy = 4.5×10 -14 m2 Chapter 3: Verification of the Two-phase Model 67 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 Chapter 3: Verification of the Two-phase Model 68 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 Chapter 3: Verification of the Two-phase Model 69 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 Chapter 3: Verification of the Two-phase Model 70 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 Chapter 3: Verification of the Two-phase Model 71 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 Chapter 3: Verification of the Two-phase Model 72 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 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 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 9.E+05 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 y/h P ( P a ) closed-form solution FE: 10 9-3 elements FE: 20 9-4 elements 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 9.E+05 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 y/h P ( P a ) closed-form solution FE: 10 9-3 elements FE: 20 9-4 elements Chapter 3: Verification of the Two-phase Model 73 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 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 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 y/h P ( P a ) closed-form solution FE: 10 9-3 elements FE: 20 9-4 elements 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 y/h P ( P a ) closed-form solution FE: 10 9-3 elements FE: 20 9-4 elements Chapter 3: Verification of the Two-phase Model 74 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 Figure 3.18: Pressure surfaces of the two end elements at a very early stage of consolidation of the column 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 y/h P ( P a ) closed-form solution FE: 10 9-3 elements FE: 20 9-4 elements 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 9.E+05 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 y/h P ( P a ) closed-form solution FE: 9-3 elements pressure surface of element Element #9 Element #10 Chapter 3: Verification of the Two-phase Model 75 Figure 3.19: Pressure surfaces of the two end elements at a time that oscillation is well- pronounced 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 0.E+00 1.E+05 2.E+05 3.E+05 4.E+05 5.E+05 6.E+05 7.E+05 8.E+05 9.E+05 1.E+06 0.8 0.82 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 y/h P ( P a ) closed-form solution FE: 9-3 elements pressure surface of element Element #9 Element #10 Chapter 3: Verification of the Two-phase Model 76 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 Chapter 3: Verification of the Two-phase Model 77 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 Figure 3.24: Schematic diagram of two unidirectional flow problems over and through a porous wall; (a) Taylor problem, (b) Beavers-Joseph problem Porous wall x z h1 h2 V0 0   x P (a) (b) Porous wall x z h1 h2 0   x P Chapter 3: Verification of the Two-phase Model 78 Figure 3.25: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10 -6 m 2 Figure 3.26: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10 -8 m 2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 1.2 Velocity (mm/s) z/ h Analytical FE: 4 elements FE: 8 elements -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -0.1 0.1 0.3 0.5 0.7 0.9 1.1 Velocity (mm/s) z/ h Analytical FE: 4 elements FE: 8 elements Chapter 3: Verification of the Two-phase Model 79 Figure 3.27: Profile of volume-averaged velocity for the Taylor problem with μ=5 Pa.s, and S=4.5×10 -10 m 2 Figure 3.28: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10 -6 m 2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -0.1 0.1 0.3 0.5 0.7 0.9 1.1 Velocity (mm/s) z/ h Analytical FE: 4 elements FE: 8 elements -1 -0.8 -0.6 -0.4 -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) z/ h Analytical FE: 4 elements FE: 8 elements Chapter 3: Verification of the Two-phase Model 80 Figure 3.29: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10 -8 m 2 Figure 3.30: Profile of volume-averaged velocity for the Beavers-Joseph problem with μ=5 Pa.s, and S=4.5×10 -10 m 2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 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) z/ h Analytical FE: 4 elements FE: 8 elements -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 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) z/ h Analytical FE: 4 elements FE: 8 elements Chapter 3: Verification of the Two-phase Model 81 Figure 3.31: Geometry, BC, and processing cycle used for compaction of a unidirectional angle laminate on a convex tool Figure 3.32: (a) 8×4 mesh and (b) 16×6 mesh of the angle laminate using symmetry H R L f f Permeable B.C. Permeable B.C. x z H = 4.96 mm for AS4/3501-6 4.80 mm for AS4/8552 L = 63 mm R = 4.57 mm T & P time Bag pressure: 0 kPa Autoclave pressure: 540 kPa 2°C/min 2°C/min 122°C 70 min 163°C 2.5 hrs 25°C 2.5°C/minT & P (b)A B (a)A B Chapter 3: Verification of the Two-phase Model 82 Figure 3.33: Time-history of normal displacement for unidirectional AS4/3501-6 angle with [0°] fibres on a convex tool Figure 3.34: Time-history of normal displacement for unidirectional AS4/8552 angle with [0°] fibres on a convex tool -1 -0.8 -0.6 -0.4 -0.2 0 0 10 20 30 40 50 60 70 80 Time (min) u n ( m m ) @ A - Hubert @ B - Hubert @ A - mesh (a) @ B - mesh (a) @ A - mesh (b) @ B - mesh (b) -0.25 -0.2 -0.15 -0.1 -0.05 0 0 20 40 60 80 100 120 140 160 Time (min) u n ( m m ) @ A - Hubert @ B - Hubert @ A - mesh (a) @ B - mesh (a) @ A - mesh (b) @ B - mesh (b) Chapter 3: Verification of the Two-phase Model 83 Figure 3.35: Time-history of normal displacement for unidirectional AS4/3501-6 angle with [90°] fibres on a convex tool Figure 3.36: Time-history of normal displacement for unidirectional AS4/8552 angle with [90°] fibres on a convex tool -1 -0.8 -0.6 -0.4 -0.2 0 0 10 20 30 40 50 60 70 80 Time (min) u n ( m m ) @ A - Hubert @ B - Hubert @ A - mesh (a) @ B - mesh (a) @ A - mesh (b) @ B - mesh (b) -0.25 -0.2 -0.15 -0.1 -0.05 0 0 20 40 60 80 100 120 140 160 Time (min) u n ( m m ) @ A - Hubert @ B - Hubert @ A - mesh (a) @ B - mesh (a) @ A - mesh (b) @ B - mesh (b) 84 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:       d d d TtEt t e 0 )(),(,)( (4.1) where T is temperature, α is the degree of cure, and et is an effective value of time that is defined as       f fT T f f e f f e t ta a t mc e t t mc e t         )( )( )( log )( log 1 1 (4.2) Chapter 4: Integration of Modeling from Fluid to Cured Resin 85 where ft 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 Ta is the shift factor that relates the reduced time, ξ, to the time variable by the following equation dt a t T  0 1  (4.3) The shift factor, Ta , is given by       21log cTcaT  (4.4) where the coefficients 1c and 2c are defined as   211 1 1 exp aac           (4.5)     1 0 2 cTc  (4.6) In the case of 3501-6 epoxy resin; Ca  /4.11 , Ca  /0712.02 , and CT  30 0 [Kim and White (1996)]. 0T 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                  T e n rur m a t WGGGG exp 1 (4.7)                  T e n rur m a t WKKKK exp 1 (4.8) 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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 86 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)(log 00    f (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   21347.96089.03694.9  f (4.10)  is given by    0 0      p  (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 KG KG εεζ  2 (4.12) in term of stress and strain tensors, ζ and ε, or KGKG KG   (4.13) Chapter 4: Integration of Modeling from Fluid to Cured Resin 87 in vector form where a single underscore denotes a vector quantity. Assuming pseudo- viscoelastic 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 GG G   (4.14) KK K   (4.15) Summation of the above two equations leads to KG KG    (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                         G GSymm G GK GKGK GKGKGK D 0. 00 00034 0003234 000323234 (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 23122 23122 GGK GGK εεεεεζ   (4.19) or Chapter 4: Integration of Modeling from Fluid to Cured Resin 88    5 1 23122 p pGGK ζζζζζζζ  (4.20) In the expanded form , we may rewrite (4.19) in the following form                                                             002 002 220 20 20 000 00 00 000 00 00 00 000 000 00 13 12 1312 12 332223 23332223 3322 33222 11 11 332211            GG Kζ (4.21) where 2K 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 1 2 12 23 2 414 1 EGE K    (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  2112 2 212    E K (4.23) where 32 EEE  is the transverse modulus, and 23  is the Poisson‟s ratio in the plane of isotropy. 12G is the axial shear modulus, and 23G is the transverse shear modulus (shear modulus in the plane of isotropy). The two new moduli  and  are defined respectively as 2 1221 4  KE  (4.24) 1222 K (4.25) In the case of pseudo-viscoelastic formulation, the material behaviour at any time-step may be written in the form of Chapter 4: Integration of Modeling from Fluid to Cured Resin 89 1212 2323 22 12 23 2 GG GG KK G G K                  (4.26) Summing the equations in (4.26) and taking into account (4.20), we may write 12232 12232 GGK GGK      (4.27) which could also be written in the form   D (4.28) where                        12 12 23 232 232232 0. 00 000 000 000 G GSymm G GK GKGK D  (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 two- phase 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. Chapter 4: Integration of Modeling from Fluid to Cured Resin 90 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)]   01 1        c f f tt v    (4.30) and 0      cflow m m tt vv     (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 volume- averaged 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 0 1        tt m m f f flowc      vv (4.32) which may be rewritten in the indicial form of 0 1 ,,        tt vu m m f f iiii       (4.33) Chapter 4: Integration of Modeling from Fluid to Cured Resin 91 where the following notation is used for simplicity iflowiici vvvu  , (4.34) Lewis and Schrefler (1998) introduce the following relationships for the mechanical components of the last two terms of the above equation t P Kt m m m m      11   (4.35) and     t tr Kt P Kt fii f f f f          ζ 1 111    (4.36) where mK and fK 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 Iζζ P (4.37) where it is assumed that fm PPP  . 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)] Iζζ bP (4.38) where b is Biot coefficient and is defined by f fb K K b 1 (4.39) where fbK is the bulk modulus of the fibre-bed (solid skeleton in general). Lewis and Schrefler (1998) introduce the following constitutive equation Chapter 4: Integration of Modeling from Fluid to Cured Resin 92                t P K uK t tr f iifbii 1 ,  ζ (4.40) which is then substituted in (4.36) to result in t Pb K u b t P K u b t P Kt f ii f ii f f f                                       1 1 1 11 1 111 ,,  (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: 0,,               t P KK b vub mf iiii   (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 tt m m f f sv             1  (4.43) which enables us to write the mass conservation equation in (4.33) in the following compact form 0,,  sviiii vu  (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 0,,  free sv mech sviiii vu   (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 mech cv  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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 93 general micromechanics formulation may be written in the form of the following two equations  ,, fmGc GGfG  (4.46)  ,, fmKc KKfK  (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                    3 1 1 s cii s c mech sv tr tK tr tK ζζ   (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                        bP tr tKK iimf mech sv    ζ 1  (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 poro- mechanics 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 iitr ζ is not equal to cK . For 2D plane strain, the constitutive equation for the rate of the volumetric strain is Chapter 4: Integration of Modeling from Fluid to Cured Resin 94                      2 3 1 3 1 s ccii s cc mech sv tr tGK tr tGK ζζ   (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 0 ,,  idjijmi fP  (4.51) where id f are the components of the drag force that, if defined as follows, leads to the Darcy‟s law. jijid vSf 1  (4.52) and ijm  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 0,  idi fP (4.53) while the resin equilibrium equation in (4.51) that we intend to use here, is the Darcy- Brinkman (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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 95 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   01 ,,  iidjij Pfσ  (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 0,,,  ijijmjij Pσ  (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 0,,,  ijijmjij bPσ  (4.56) We may now revisit the governing equations in(2.61) and rewrite them in the form         0 0 0 ,,, ,, ,, ijijmjij idjijmi sviiii bPσ fP vu    (4.57) As was mentioned earlier in this section, in most cases the term involving the deviatoric stress of the matrix (i.e. ijm  ) are negligible compared to id f 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 ijijmijijs bPσ   (4.58) Let us define the „s-p‟ stress as ijmijijijsijps σbPσ   (4.59) 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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 96 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          0 0 0 ,, , ,, ijijps idi sviiii bPσ fP vu   (4.60) Assuming elastic behaviour for the fibre-bed, a constitutive relationship for Biot effective stress tensor may be written as klijklfbij C   (4.61) where ijklfb C 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 thermo- chemical changes from fluid to a solid cross-linked state. Conceptually, this new stress tensor ijps  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 ijps  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, ijps  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 ijps  except for its hydrostatic pressure (bulk) component. Further development will show that ijps  plays an essential role in linking the modeling of the porous flow-compaction and the fully solid behaviour using micromechanics. Chapter 4: Integration of Modeling from Fluid to Cured Resin 97 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 ijps  may be assumed in the form of fbps KK  (4.62) for the bulk modulus, and fbcps GGG  (4.63) for the shear modulus. The index fb represents the fibre-bed or the solid skeleton of the system, and cG 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, ijps  . 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 KpsGpsKpsGpsps KG    (4.64) The fibre-bed is assumed to be elastic. The constitutive relations for the two stress components may be written as   GfbcGps GG    (4.65) and KfbKps K    (4.66) Summation of the above two equations leads to   KfbGfbcps KGG    (4.67) Chapter 4: Integration of Modeling from Fluid to Cured Resin 98 which could also be written in the form   psps D   (4.68) where                                             fbc fbc fbc fbcfb fbcfbfbcfb fbcfbfbcfbfbcfb ps GG GGSymm GG GGK GGKGGK GGKGGKGGK D 0. 00 000 3 4 000 3 2 3 4 000 3 2 3 2 3 4 (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 ijijpsijs bPσ    (4.70) Chapter 4: Integration of Modeling from Fluid to Cured Resin 99 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. bP G GKGK GKGK c cfbcfb cfbcfb                                             0 1 1 00 03432 03234 12 2 1 12 2 1       (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  32 cc smech sv GK tr   ζ  (4.72) which upon considering the components of ζs in (4.71) results in  32 21 cc mech sv GK      (4.73) In the absence of flow and free strains, the mass conservation equation then becomes  32 21 21 cc GK      (4.74) Combining the mass conservation in the above equation with (4.71) we have       bPGKGK cfbcc 23232 2121   (4.75) which after simplification leads to   21   cfb KKbP (4.76) Chapter 4: Integration of Modeling from Fluid to Cured Resin 100 If the above equation for the pressure is substituted into (4.71), we will arrive at                                  12 2 1 12 2 1 00 03432 03234       c cccc cccc G GKGK GKGK (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  32 2 2121 cc f GK      (4.78) where f is the free strain in either 1 or 2 directions. Combining the above with (4.71) we eventually arrive at      fcccfb GKKKbP  3221  (4.79) Again, if the above equation is substituted into (4.71) we will have                                    12 2 1 12 2 1 00 03432 03234       f f c cccc cccc G GKGK GKGK (4.80) Chapter 4: Integration of Modeling from Fluid to Cured Resin 101 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 ijIijpsijs Pb    (4.81) where Ib 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)]            b b b1 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 f fb K K b 2 2 1 (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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 102 P b b b pss                                                                0 0 0 1 12 13 23 33 22 11 12 13 23 33 22 11             (4.84) for transversely isotropic materials. With a slight change from the 3 rd equation in (4.60), the total equilibrium equation of the system is now written as   0,,  iIjijps Pb (4.85) Taking a similar approach to the isotropic case, we decompose the „s-p‟ stress to have 23122 G psGpsKpspspsps   ζζζζζζ  (4.86) In vector form, we may write the foregoing equation as 23122 23122 23122 GpsGpsKpspsps GpsGpsKpspspsps GGK            (4.87) A reasonable representation for the elastic moduli related to „s-p‟ stress may be assumed in the form fbps KK 22  (4.88) fbps   (4.89) for the moduli dealing with the bulk behaviour of the material, and   cps (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. fbcps GGG 121212  (4.91) Chapter 4: Integration of Modeling from Fluid to Cured Resin 103 fbcps GGG 232323  (4.92) As it is seen here, the concept is the same as in the isotropic case. Assuming pseudo- viscoelastic behaviour for the resin, and therefore the composite, and that the fibre-bed is elastic, we write 22 2 KfbKps K    (4.93)       cps (4.94)    fbps  (4.95)   1212 1212 GfbcGps GG    (4.96)   2323 2323 GfbcGps GG    (4.97) Summation of the five components of the rate of „s-p‟ stress vector leads to   psps D   (4.98) where                                   fbc fbc fbc fbcfb fbcfbfbcfb fbfbc ps GG GGSymm GG GGK GGKGGK D 1212 1212 2323 23232 2323223232 0. 00 000 000 000 (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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 104               fb fbfbfb fbc ps G GKD 12 232 00 0 0   (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               c cfbfb fbc ps G GKD 12 232 00 0 0   (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 time- step 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. Chapter 4: Integration of Modeling from Fluid to Cured Resin 105                                12 2 1 12 232 12 2 1 00 0 0       c ccc cc G GK  (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. Pb b G GK c cfbfb fbc                                             000 0 0 1 12 2 1 12 232 12 2 1         (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           22321223212 211211   ccccfbfb ccfbc GKbPGK Pb   (4.104) The above may be rewritten as           2221 211   cfbcfb cfb KKbP Pb   (4.105) As the above should be valid for any combination of strain components, the following must be always correct       cfb cfb cfb KKb b 22 1          (4.106) Chapter 4: Integration of Modeling from Fluid to Cured Resin 106 that leads to the values for b1 and  that are consistent with the two constitutive equations being equivalent, i.e. b KK b cfb cfb 22 1     (4.107)   cfb cfb KK 22 2      (4.108) Substituting the above results into (4.103), the constitutive equation may be written as         bP KK G GK KK cfbcfb c cfbfb fbcfbcfbc                                              0 1 00 0 0 22 12 2 1 12 232 22 2 12 2 1          (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                                        2 1232 2 2322 1 1 2322 1 1       cc ccc ccccccc cc GK GKGK     (4.110) Adding the two normal components of strain, the volumetric strain of the system is obtained as       212322 232 21 1    ccccc cccc v GK GK      (4.111) We substitute the stress components in the above equation with the relationships obtained by our approach in (4.109) to arrive at Chapter 4: Integration of Modeling from Fluid to Cured Resin 107                                                 bPGK bP KKKK GK GK cfbfbcc cfb cfb fb cfb cfb cccc cccc 22321 22 21 22 2 232 2 232 21 1           (4.112) which after extensive manipulations leads to     2221  cfbcfb KKbP   (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 f 11   and f 22   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 mechv f vv   (4.114) to the constitutive equation of two-phase systems in the form of Pb b G GK f f c cfbfb fbc                                               000 0 0 1 12 22 11 12 232 12 2 1         (4.115) leads to an approach equivalent to using the common constitutive equation of transversely isotropic materials which follows                                  12 22 11 12 232 12 2 1 00 0 0       f f c ccc cc G GK  (4.116) Chapter 4: Integration of Modeling from Fluid to Cured Resin 108 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                   f cc f ccfbfb f c f cfbc GKbPGK Pb 2223211223212 22111211     (4.117) that may be rewritten as              f cc f ccfbcfb f c f ccfb GKKKbP Pb 223212221 21211     (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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 109         cc c c c cfb cfb cfb GKKKb b 23222 1              (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                         100 0 0 100 0 0 CS SC CosSin SinCos   n (4.120) The transformation between the stress tensor in the laminate coordinates, ζ̂ , and the stress tensor in the material coordinate systems takes the form ζnnζ Tˆ (4.121) or in indicial notations klljkiij nn  ˆ (4.122) Using the Voigt notation, the stress vector in the laminate coordinates 321  may be written as Chapter 4: Integration of Modeling from Fluid to Cured Resin 110                                           12 13 23 33 22 11 1 12 13 23 33 22 11 ˆ ˆ ˆ ˆ ˆ ˆ              T (4.123) where the transformation matrix is defined by                          22 22 22 1 000 0000 0000 000100 2000 2000 SCCSCS CS SC CSCS CSSC T  (4.124) To transform the stress vector back to the material coordinate system we need to have                         22 22 22 000 0000 0000 000100 2000 2000 SCCSCS CS SC CSCS CSSC T  (4.125) The stress vector in the laminate coordinates 321  may be written as                                           12 13 23 33 22 11 1 12 13 23 33 22 11 ˆ ˆ ˆ ˆ ˆ ˆ              T (4.126) where Chapter 4: Integration of Modeling from Fluid to Cured Resin 111                         22 22 22 00022 0000 0000 000100 000 000 SCCSCS CS SC CSCS CSSC T  (4.127) The matrix of elastic moduli in the laminate coordinate system may be written in the form  TDTD 1ˆ  (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   P bbCS b bCbS bSbC P b b b TT pspss                                                                    1 2 1 2 2 1 2 1 1 0 0 0 0 0   (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. Chapter 4: Integration of Modeling from Fluid to Cured Resin 112 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                            fbc fbcfb fbcfbfbcfb ps GGSymm GGK GGKGGK D . 0 3 4 0 3 2 3 4 (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 0,,    dwdwdwudwv p free svp mech svpiipii   (4.131) where mech sv  can be obtained from (4.50). Substituting (4.59) into (4.50) leads to             bP tr tGK ps cc mech sv 2 3 1 ζ  (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 183 0CK pvpu (4.133) 183   d TTp BδNKpv (4.134)   183183 32 1       d GK d ps TTp cc TTp BDδNBδNCpu (4.135)   33 3     d GK b pTp cc NNCpp (4.136) 33 0Kpp (4.137) Chapter 4: Integration of Modeling from Fluid to Cured Resin 113 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 13 1      dd TpfreeTfree psc TTp NεδεDDδNfp  (4.138) 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 118 )(    dd freeTnT d ζBtNfu  (4.139) Combining the above with the implementation of the other two governing equations leads to the matrix equations of an element in the form of 139                                                       p v u pppu vu uu pv vpvv upuvuu f f f p v u C0C 00C 00C p v u 0K0 KK0 KKK t (4.140) Transversely isotropic materials: In 2D, the global x-y plane is equivalent to 31  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.)                   12 2 23 2 232 232 22 12 22 232 4224 . 0 042 GCGSSymm GK GKSCGSCGKSSCC D  (4.141) In special cases, we may be able to simplify the above equation, e.g.            12 232 . 0 0 0 GSymm GKD   (4.142) Chapter 4: Integration of Modeling from Fluid to Cured Resin 114              23 232 232232 . 0 0 90 GSymm GK GKGK D (4.143) In the x-y plane, (4.129) is simplified to Pb bSbC psxy y x sxy y x                                   0 2 1 2       (4.144) As the plane strain conditions are assumed (i.e. 0 ;0  free z mech z  ), the volumetric strain of the composite and its rate are respectively expressed by mech y mech x mech sv   (4.145) mech y mech x mech sv    (4.146) Using the expanded form of the constitutive equation of the composite material, we may write         yscxscyscxsc mech sv DDDD   1 22 1 12 1 12 1 11   (4.147) We expand the first two rows of (4.144) to have           bP PbSbC ypsys xpsxs       21 2 (4.148) and then the above is substituted into (4.147) to result in            PDDbDDbSbC DDDD cccc ypsccxpscc mech sv   1 22 1 12 1 12 1 11 2 1 2 1 22 1 12 1 12 1 11         (4.149) where  freeklklijklps mech klijklpsijps DDσ    (4.150) Chapter 4: Integration of Modeling from Fluid to Cured Resin 115 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   183 1      dpsc TTp BDDBδNCpu (4.151)       33 1 22 1 12 1 12 1 11 2 1 2     dDDbDDbSbC pTp cccc NNCpp (4.152) 183 2 1 2 0                 db bSbC pT NBK up (4.153) 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                                                          n k nnn k nnn nn k n k n k n k n k n t tt ttt pCuuCf uuCf ff p v uu CKC KKC KKKC pppup vuv uu pppvpu vpvvvu upuvuuuu 11 11 int1 1 1 1 1 1 1 1 1 (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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 116 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 1 st equation of (2.62) may be rewritten in the form flowvviiii vu    or ,, (4.155) Integrating the above equation through time leads to  t iiflowvv dtv 0 ,  (4.156) 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 free sv mech svii t iiflowv dtudtv 0 , 0 ,   (4.157) and substituted into the following equation (similar to (2.119)) for updating the fluid volume fraction flowv flowv       1 0 (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 Chapter 4: Integration of Modeling from Fluid to Cured Resin 117 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 thermo- mechanical 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. Chapter 4: Integration of Modeling from Fluid to Cured Resin 118 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 119 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 cKfb KCoeffK  (5.1) where CoeffK provides a simple representation of the value of fibre-bed bulk modulus. Chapter 5: Numerical Applications 120 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 4.0 ,3.0 1 ,10 )4.0(or 4.0 ,6.0     mf mf mf GPaEGPaE    kPaf mmLmmh 700 12 ,4 0   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 GPa GG GG G mffm fm c 7837.0    GPa KK KK K mffm fm c 205.3    The matrix of elastic moduli is then expressed as GPa SymmGSymm GK GKGK D c cc cccc                               7837.0. 02501.4 06827.22501.4 . 0 3 4 0 3 2 3 4 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                     y x GPa kPa   0 2501.46827.2 6827.22501.4 700 and solve for the unknowns in stress and strain fields Chapter 5: Numerical Applications 121       kPax y 84.441 10647.1 4   The strain energy of the system may be calculated by      mmmmmMPaVolE yys 112410647.144184.0 2 1 . 2 1 4  mNEs .002767.0 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 .1000 ,400 21 kPafkPaf  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 Chapter 5: Numerical Applications 122 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 stsPa 18 ,.101 7  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. Chapter 5: Numerical Applications 123 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 .100 ,700 21 kPafkPaf  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). Chapter 5: Numerical Applications 124 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 9- 4 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 Chapter 5: Numerical Applications 125 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 Chapter 5: Numerical Applications 126 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 Chapter 5: Numerical Applications 127 cKfb KCoeffK 22  (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 398.0 6198.7 7262.2 4055.6 068.4 253.0 4.127 23 2 23 2 12 12 1          GPaE GPaG GPaK GPaG GPaE consisting of fibres with the following properties GPaKGPaG GPaE GPaE 57.11 25.0 2.0 6.27 2.17 210 2 23 12 12 2 1                 Other specifications of the problem include 4.0 ,6.0  mf  and kPaf mmLmmh 700 12 ,4 0   Chapter 5: Numerical Applications 128 Substituting the properties of the composite material into (4.142) leads to the matrix of elastic moduli of the problem GPa Symm D            068.4. 01317.9 02412.304.129 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                     y x GPa kPa   0 1317.92412.3 2412.304.129 700 and solve for the unknowns in stress and strain fields       kPax y 46.248 10666.7 5   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 fb K2 . 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 .100 ,700 21 kPafkPaf  Chapter 5: Numerical Applications 129 Three different situations for the direction of the fibres are assumed;  90 and ,45 ,0 . 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 Chapter 5: Numerical Applications 130 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/3501- 6. 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 Chapter 5: Numerical Applications 131 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 Chapter 5: Numerical Applications 132 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 Chapter 5: Numerical Applications 133 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 Chapter 5: Numerical Applications 134 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 non- existent 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 Chapter 5: Numerical Applications 135 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 Chapter 5: Numerical Applications 136 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. Chapter 5: Numerical Applications 137 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 Chapter 5: Numerical Applications 138 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 time- steps 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 non- uniform 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. Chapter 5: Numerical Applications 139 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 Chapter 5: Numerical Applications 140 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 Chapter 5: Numerical Applications 141 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. Chapter 5: Numerical Applications 142 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 time- history 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. Chapter 5: Numerical Applications 143 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. Chapter 5: Numerical Applications 144 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 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 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 6×4 9×6 Mesh Resin velocity conditions 2×1 3×2 CoeffK = 1 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 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 24×16 1.6404 48×32 1.646 96×64 1.6482 192×128 1.649 12×8 Mesh uA(×10 -7 m) ABAQUS 1.5749 9×6 uA (×10 -7 m) 3×2 6×4 Resin velocity conditions Chapter 5: Numerical Applications 145 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 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 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 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 Strain energy (N.mm)Resin velocity conditions 3×2 12×8 6×4 9×6 Mesh CoeffK=1 0.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 24×16 1.6404 48×32 1.646 96×64 1.6482 192×128 1.649 Mesh uA(×10 -7 m) ABAQUS Boundary Conditions 12×8 6×4 1.5749 uA (×10 -7 m) 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 Strain energy (N.mm) 12×8 Mesh Boundary Conditions 6×4 Chapter 5: Numerical Applications 146 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 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 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 192×128 θ = 0° The displacement values are in 10 -8 m, and the strain energy values are in 10 -4 N.m. 12×8 Mesh ABAQUS (disp.) 9.7984 9×6 CoeffK 3×2 6×4 Output type 1 0.1 0.01 0.001 0.0001 0.00001 displacement 0.8639 1.3198 1.4509 1.4673 1.469 1.4692 strain energy 4.4961 4.6509 4.6917 4.6967 4.6972 4.6972 displacement 1.0823 1.1804 1.2037 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 192×128 1.0485 θ = 90° The displacement values are in 10 -8 m, and the strain energy values are in 10 -4 N.m. ABAQUS (disp.) 3×2 12×8 6×4 9×6 Mesh Output type CoeffK Chapter 5: Numerical Applications 147 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 Table 5.9: Resin properties in Examples 1 & 2, where Poisson‟s ratio is assumed constant Property 3501-6 resin E u (GPa) 3.2 E r (GPa) 0.031 ν 0.35 Table 5.10: Mechanical properties of fibre and resin in the examples of Section 5.2 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 96×64 9.0092 192×128 θ = 45° The displacement values are in 10 -8 m, and the strain energy values are in 10 -4 N.m. 9×6 12×8 Mesh ABAQUS (disp.) Output type 6×4 8.9136 3×2 CoeffK Property AS4 fibre Property 3501-6 resin 12 0.2 K u (GPa) 3.556 23 0.3 K r (GPa) 0.8 E1 (GPa) 207 G u (GPa) 1.185 E2 (GPa) 20.7 G r (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 Chapter 5: Numerical Applications 148 5.4. Figures Figure 5.1: The geometry and specifications of the composite sample under constant pressure load (uni-axial strain condition) Figure 5.2: The geometry and specifications of the composite sample under linear pressure load f0 L h x y x y h f1 L f2 Chapter 5: Numerical Applications 149 Figure 5.3: The profile of the vertical displacement at the top of the constrained sample under linear distribution of pressure load Figure 5.4: Composite sample under pressure gradient -1.E-03 -9.E-04 -8.E-04 -7.E-04 -6.E-04 -5.E-04 -4.E-04 -3.E-04 -2.E-04 0 2 4 6 8 10 12 14 x (mm) u y ( m m ) 1×2 mesh 3×2 mesh 6×4 mesh f1 L 2h f2 A f1 = 700 kPa f2 = 100 kPa h = 4 mm L = 12 mm 4.0 205.3 7837.0     GPaK GPaG c c z x Chapter 5: Numerical Applications 150 Figure 5.5: Shape of the displacement profiles of the sample under pressure gradient 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 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 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) y ( m m ) Left side's displacement profile Right side's displacement profile 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 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) y ( m m ) 3×2 mesh/ 9-3 elements/ Coeff.K=1 3×2 mesh/ 9-3 elements/ Coeff.K=0.1 3×2 mesh/ 9-3 elements/ Coeff.K=0.01 3×2 mesh/ 9-3 elements/ Coeff.K=0.0001 Chapter 5: Numerical Applications 151 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 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 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 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) y ( m m ) 6×4 mesh/ 9-3 elements/ Coeff.K=1 6×4 mesh/ 9-3 elements/ Coeff.K=0.1 6×4 mesh/ 9-3 elements/ Coeff.K=0.01 6×4 mesh/ 9-3 elements/ Coeff.K=0.0001 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 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) y ( m m ) 12×8 mesh/ 9-3 elements/ Coeff.K=1 12×8 mesh/ 9-3 elements/ Coeff.K=0.1 12×8 mesh/ 9-3 elements/ Coeff.K=0.01 12×8 mesh/ 9-3 elements/ Coeff.K=0.0001 Chapter 5: Numerical Applications 152 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 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 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 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) y ( m m ) 9-3 elements/ 3×2 mesh 9-3 elements/ 6×4 mesh 9-3 elements/ 12×8 mesh 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 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) y ( m m ) 9-3 elements/ 3×2 mesh 9-3 elements/ 6×4 mesh 9-3 elements/ 12×8 mesh Chapter 5: Numerical Applications 153 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 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 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.0E+00 5.0E-08 1.0E-07 1.5E-07 2.0E-07 2.5E-07 ux (m) y ( m m ) 9-3 elements 9-4 elements/Permeable B.C. 9-4 elements/Impermeable B.C. 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 0.0E+00 5.0E-08 1.0E-07 1.5E-07 2.0E-07 2.5E-07 ux (m) y ( m m ) 9-3 elements 9-4 elements/Permeable B.C. 9-4 elements/Impermeable B.C. Chapter 5: Numerical Applications 154 Figure 5.13: Pressure distribution of the sample shown in Figure 5.4, using a 6×4 mesh of 9-3 elements 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 00.0020.0040.006 0.0080.010.012 0 1 2 3 4 x 10 -3 0 1 2 3 4 5 6 7 x 10 5 x (m)y (m) P ( P a ) 00.0020.0040.0060.0080.010.012 0 1 2 3 4 x 10 -3 0 1 2 3 4 5 6 7 8 x 10 5 x (m)y (m) P ( P a ) Chapter 5: Numerical Applications 155 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 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 00.0020.0040.006 0.0080.010.012 0 1 2 3 4 x 10 -3 -1 0 1 2 3 4 5 6 7 x 10 5 x (m)y (m) P ( P a ) 00.0020.0040.006 0.0080.010.012 0 1 2 3 4 x 10 -3 -1 0 1 2 3 4 5 6 7 x 10 5 x (m)y (m) P ( P a ) Chapter 5: Numerical Applications 156 Figure 5.17: Time-history of the applied temperature and the resulting degree of cure in Examples 1-7 Figure 5.18: Time-history of the applied temperature and the resulting resin viscosity in Examples 1-7 0 20 40 60 80 100 120 140 160 180 200 0 50 100 150 200 250 time (min) T e m p e ra tu re ( °C ) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 D e g re e o f c u re Temperature Degree of cure 0 20 40 60 80 100 120 140 160 180 200 0 50 100 150 200 250 time (min) T e m p e ra tu re ( °C ) 1E-1 1E+0 1E+1 1E+2 1E+3 V is c o s it y ( P a .s ) Temperature Resin viscosity Chapter 5: Numerical Applications 157 Figure 5.19: Time-history of the development of longitudinal stress in Example 1. The predictions are essentially identical. Figure 5.20: Time-history of the development of transverse stress in Examples 1 and 2. The predictions are essentially identical. 0 5 10 15 20 25 30 35 40 45 50 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model (Zobeiry) -10 0 10 20 30 40 50 60 70 80 0 50 100 150 200 250 time (min) σ z ( M P a ) Integrated model Stress model (Zobeiry) Chapter 5: Numerical Applications 158 Figure 5.21: Time-history of the development of x in Example 2. The predictions are essentially identical. Figure 5.22: Time-history of the development of x in Example 3. The predictions are essentially identical. 0 10 20 30 40 50 60 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model (Zobeiry) -20 -10 0 10 20 30 40 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model Chapter 5: Numerical Applications 159 Figure 5.23: Time-history of the development of z in Examples 3 and 4. The predictions are essentially identical. Figure 5.24: Time-history of the development of resin pressure in Examples 3 and 4 -60 -40 -20 0 20 40 60 0 50 100 150 200 250 time (min) σ z ( M P a ) Integrated model Stress model -40 -30 -20 -10 0 10 20 30 40 50 60 0 50 100 150 200 250 time (min) P ( M P a ) Chapter 5: Numerical Applications 160 Figure 5.25: Time-history of the development of x in Example 4. The predictions are essentially identical. Figure 5.26: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 5 -30 -20 -10 0 10 20 30 40 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model 0.E+00 1.E-05 2.E-05 3.E-05 4.E-05 5.E-05 6.E-05 7.E-05 8.E-05 0 50 100 150 200 250 time (min) R e s in v e lo c it y ( m m /s ) Chapter 5: Numerical Applications 161 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 Figure 5.28: Time-history of the changes in resin volume fraction in Example 5 0 25 50 75 100 125 150 175 200 15 20 25 30 35 40 45 50 55 time (min) T e m p e ra tu re ( °C ) 0.0E+0 1.0E-4 2.0E-4 3.0E-4 4.0E-4 5.0E-4 6.0E-4 7.0E-4 8.0E-4 R a te o f fr e e v o lu m e tr ic s tr a in ( /m in ) Temperature Rate of free vol. strain 0.36 0.38 0.4 0.42 0.44 0 50 100 150 200 250 time (min) φ Integrated model Stress model Chapter 5: Numerical Applications 162 Figure 5.29: Time-history of the development of resin pressure in Example 5 Figure 5.30: Time-history of the development of x in Example 5 -80 -70 -60 -50 -40 -30 -20 -10 0 10 0 50 100 150 200 250 time (min) P ( M P a ) -20 -10 0 10 20 30 40 50 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model Chapter 5: Numerical Applications 163 Figure 5.31: Time-history of the development of z in Example 5 Figure 5.32: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 6 -60 -40 -20 0 20 40 60 80 100 0 50 100 150 200 250 time (min) σ z ( M P a ) Integrated model Stress model 0.0E+00 2.0E-05 4.0E-05 6.0E-05 8.0E-05 1.0E-04 1.2E-04 1.4E-04 1.6E-04 0 50 100 150 200 250 time (min) R e s in v e lo c it y ( m m /s ) Chapter 5: Numerical Applications 164 Figure 5.33: Time-history of the changes in resin volume fraction in Example 6 Figure 5.34: Time-history of the development of resin pressure in Example 6 0.36 0.38 0.4 0.42 0.44 0 50 100 150 200 250 time (min) φ Integrated model Stress model -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 0 50 100 150 200 250 time (min) P ( M P a ) Chapter 5: Numerical Applications 165 Figure 5.35: Time-history of the transverse strain in Example 6 Figure 5.36: Time-history of the development of x in Example 6 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0 50 100 150 200 250 time (min) ε z Integrated model Stress model -5 0 5 10 15 20 25 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model Chapter 5: Numerical Applications 166 Figure 5.37: Time-history of the development of resin pressure in Example 7 Figure 5.38: Time-history of the changes in resin velocity in the vertical direction at the top surface of Example 7 0 0.1 0.2 0.3 0.4 0.5 0.6 0 50 100 150 200 250 time (min) P ( M P a ) 0.0E+00 2.0E-04 4.0E-04 6.0E-04 8.0E-04 1.0E-03 1.2E-03 1.4E-03 0 50 100 150 200 250 time (min) R e s in v e lo c it y ( m m /s ) Chapter 5: Numerical Applications 167 Figure 5.39: Time-history of the changes in resin volume fraction in Example 7 Figure 5.40: Time-history of the transverse strain in Example 7 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0 50 100 150 200 250 time (min) φ Integrated model Stress model Stress model - final φ -0.2 -0.16 -0.12 -0.08 -0.04 0 0.04 0 50 100 150 200 250 time (min) ε z Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 168 Figure 5.41: Time-history of the development of x in Example 7 Figure 5.42: Geometry and BC of the flat laminate undergoing cure -10 -5 0 5 10 15 20 25 0 50 100 150 200 250 time (min) σ x ( M P a ) Integrated model Stress model Stress model - final φ x z 63 mm 4.96 mm Permeable BC f A B C D Chapter 5: Numerical Applications 169 Figure 5.43: Time-history of the autoclave and part temperatures and the resulting degree of cure in the flat and angle laminate examples 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 0 20 40 60 80 100 120 140 160 180 200 0 50 100 150 200 250 300 350 time (min) T e m p e ra tu re ( °C ) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 D e g re e o f c u re Autoclave Temp. Part Temp. Degree of cure 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -5 0 5 10 15 20 25 30 35 σx (MPa) z ( m m ) [0°] Stress model [0°] ABAQUS [90°] Stress Model [90°] ABAQUS Chapter 5: Numerical Applications 170 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 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 0 10 20 30 40 50 60 0 10 20 30 40 50 60 70 x (mm) σ x ( M P a ) Stress model ABAQUS 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.2 0.4 0.6 0.8 1 φ z ( m m ) [0°] lay-up [45°] lay-up [90°] lay-up Initial φ Chapter 5: Numerical Applications 171 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 Chapter 5: Numerical Applications 172 Figure 5.50: Final profile of σx at section AB of the [0°] flat laminate Figure 5.51: Final profile of σx at section AB of the [45°] flat laminate 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 -10 -5 0 5 10 15 20 σx (MPa) z ( m m ) Integrated model Stress model Stress model - final φ 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 5 10 15 20 25 σx (MPa) z ( m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 173 Figure 5.52: Final profile of σx at section AB of the [90°] flat laminate Figure 5.53: Final distribution of axial force along the length of the [0°] flat laminate 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5 10 15 20 25 30 35 σx (MPa) z ( m m ) Integrated model Stress model Stress model - final φ -20 -15 -10 -5 0 5 10 15 0 10 20 30 40 50 60 70 x (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 174 Figure 5.54: Final distribution of bending moment along the length of the [0°] flat laminate Figure 5.55: Final distribution of axial force along the length of the [45°] flat laminate -10 -8 -6 -4 -2 0 2 4 6 0 10 20 30 40 50 60 70 x (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ 0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 x (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 175 Figure 5.56: Final distribution of bending moment along the length of the [45°] flat laminate Figure 5.57: Final distribution of axial force along the length of the [90°] flat laminate -30 -25 -20 -15 -10 -5 0 0 10 20 30 40 50 60 70 x (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ 0 20 40 60 80 100 120 140 160 0 10 20 30 40 50 60 70 x (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 176 Figure 5.58: Final distribution of bending moment along the length of the [90°] flat laminate Figure 5.59: Geometry and boundary conditions of the convex angle laminate -80 -70 -60 -50 -40 -30 -20 -10 0 10 0 10 20 30 40 50 60 70 x (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ H R L f f Permeable B.C. Permeable B.C. x z H = 4.96 mm L = 63 mm R = 4.57 mm υ = 0.442 f = 540 kPa Chapter 5: Numerical Applications 177 Figure 5.60: Geometry and boundary conditions of the concave angle laminate Figure 5.61: Location of key points on the geometry of the angle laminate, and transformation of geometry into the ξη coordinates H R f f Permeable B.C. L x z H = 4.96 mm L = 63 mm R = 1.39 mm υ = 0.442 f = 540 kPa A D C F B E ξ η Chapter 5: Numerical Applications 178 Figure 5.62: FE meshes used in the analysis of angle laminates 24×16b 24×16a 24×16 20×12 16×8 Chapter 5: Numerical Applications 179 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 0 1 2 3 4 5 0.3 0.305 0.31 0.315 0.32 0.325 0.33 φ η ( m m ) 16×8 mesh 20×12 mesh 24×16 mesh 24×16a mesh 24×16b mesh Chapter 5: Numerical Applications 180 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 Chapter 5: Numerical Applications 181 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 Chapter 5: Numerical Applications 182 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 Chapter 5: Numerical Applications 183 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 Chapter 5: Numerical Applications 184 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 Chapter 5: Numerical Applications 185 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 Chapter 5: Numerical Applications 186 Figure 5.77: Time-history of the development of tangential stress at different locations of the convex [0°] angle laminate Figure 5.78: Time-history of the development of tangential stress at point C on the convex [0°] angle laminate -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) @ C @ D @ A @ B -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 187 Figure 5.79: Time-history of the development of tangential stress at point D on the convex [0°] angle laminate Figure 5.80: Time-history of the development of tangential stress at point A on the convex [0°] angle laminate -4 -2 0 2 4 6 8 10 12 14 16 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ -160 -140 -120 -100 -80 -60 -40 -20 0 20 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 188 Figure 5.81: Time-history of the development of tangential stress at point B on the convex [0°] angle laminate Figure 5.82: Time-history of the development of tangential stress at different locations of the convex [45°] angle laminate -30 -25 -20 -15 -10 -5 0 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) @ C @ D @ A @ B Chapter 5: Numerical Applications 189 Figure 5.83: Time-history of the development of tangential stress at point C on the convex [45°] angle laminate Figure 5.84: Time-history of the development of tangential stress at point D on the convex [45°] angle laminate -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ -5 0 5 10 15 20 25 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 190 Figure 5.85: Time-history of the development of tangential stress at point A on the convex [45°] angle laminate Figure 5.86: Time-history of the development of tangential stress at point B on the convex [45°] angle laminate -80 -70 -60 -50 -40 -30 -20 -10 0 10 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ -10 -5 0 5 10 15 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 191 Figure 5.87: Time-history of the development of tangential stress at different locations of the convex [90°] angle laminate Figure 5.88: Time-history of the development of tangential stress at point C on the convex [90°] angle laminate -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) @ C @ D @ A @ B -5 0 5 10 15 20 25 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 192 Figure 5.89: Time-history of the development of tangential stress at point D on the convex [90°] angle laminate Figure 5.90: Time-history of the development of tangential stress at point A on the convex [90°] angle laminate -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ -2 0 2 4 6 8 10 12 14 16 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 193 Figure 5.91: Time-history of the development of tangential stress at point B on the convex [90°] angle laminate Figure 5.92: Time-history of the development of tangential stress at point F on the convex [90°] angle laminate -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ ξ (M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 194 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 Chapter 5: Numerical Applications 195 Figure 5.95: Final distribution of tangential stress for the concave [90°] angle laminate obtained by the integrated model Figure 5.96: Time-history of the development of tangential stress at different locations of the concave [0°] angle laminate -50 0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) @ C @ D @ A @ B Chapter 5: Numerical Applications 196 Figure 5.97: Time-history of the development of tangential stress at point C on the concave [0°] angle laminate Figure 5.98: Time-history of the development of tangential stress at point D on the concave [0°] angle laminate -50 0 50 100 150 200 250 300 350 400 450 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ -2 0 2 4 6 8 10 12 14 16 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 197 Figure 5.99: Time-history of the development of tangential stress at point A on the concave [0°] angle laminate Figure 5.100: Time-history of the development of tangential stress at point B on the concave [0°] angle laminate -5 0 5 10 15 20 25 30 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ -4 -2 0 2 4 6 8 10 12 14 16 0 50 100 150 200 250 300 350 time (min) σ 1 1 ( M P a ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 198 Figure 5.101: Final distribution of axial force along the length of the convex [0°] angle laminate Figure 5.102: Final distribution of bending moment along the length of the convex [0°] angle laminate -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 0 10 20 30 40 50 60 70 distance from corner (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ -250 -200 -150 -100 -50 0 50 0 10 20 30 40 50 60 70 distance from corner (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 199 Figure 5.103: Final distribution of axial force along the length of the convex [45°] angle laminate Figure 5.104: Final distribution of bending moment along the length of the convex [45°] angle laminate -60 -40 -20 0 20 40 60 80 0 10 20 30 40 50 60 70 distance from corner (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ -160 -140 -120 -100 -80 -60 -40 -20 0 0 10 20 30 40 50 60 70 distance from corner (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 200 Figure 5.105: Final distribution of axial force along the length of the convex [90°] angle laminate Figure 5.106: Final distribution of bending moment along the length of the convex [90°] angle laminate 0 20 40 60 80 100 120 140 160 0 10 20 30 40 50 60 70 distance from corner (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ -60 -50 -40 -30 -20 -10 0 10 0 10 20 30 40 50 60 70 distance from corner (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 201 Figure 5.107: Final distribution of axial force along the length of the concave [0°] angle laminate Figure 5.108: Final distribution of bending moment along the length of the concave [0°] angle laminate -50 0 50 100 150 200 250 300 0 10 20 30 40 50 60 70 distance from corner (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ -350 -300 -250 -200 -150 -100 -50 0 50 0 10 20 30 40 50 60 70 distance from corner (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 202 Figure 5.109: Final distribution of axial force along the length of the concave [45°] angle laminate Figure 5.110: Final distribution of bending moment along the length of the concave [45°] angle laminate 0 50 100 150 200 250 300 0 10 20 30 40 50 60 70 distance from corner (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ -200 -150 -100 -50 0 50 0 10 20 30 40 50 60 70 distance from corner (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ Chapter 5: Numerical Applications 203 Figure 5.111: Final distribution of axial force along the length of the concave [90°] angle laminate Figure 5.112: Final distribution of bending moment along the length of the concave [90°] angle laminate 0 50 100 150 200 250 300 0 10 20 30 40 50 60 70 distance from corner (mm) A x ia l fo rc e ( N /m m ) Integrated model Stress model Stress model - final φ -80 -60 -40 -20 0 20 40 60 80 0 10 20 30 40 50 60 70 distance from corner (mm) M o m e n t (N .m m /m m ) Integrated model Stress model Stress model - final φ 204 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. Chapter 6: Conclusions and Further Work 205 4. The proposed two-phase approach was adapted so that it could successfully represent the stress development in a curing composite based on the pseudo- viscoelastic 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 Chapter 6: Conclusions and Further Work 206 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 flow- compaction 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. Chapter 6: Conclusions and Further Work 207 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)]. Chapter 6: Conclusions and Further Work 208  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 3- phase 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. 209 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 process- induced 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. 1106- 1117, 7, 2008. Arafath, A.R.A., Vaziri, R., and Poursartip, A., \"Closed-form solution for process- induced 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. Bibliography 210 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 Thick- Section 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. Bibliography 211 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. Bibliography 212 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. Bibliography 213 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. 2 nd 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. Bibliography 214 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 one- dimensional models.\" International Journal for Numerical and Analytical Methods in Geomechanics, vol. 10, pp. 461, 1986a. Bibliography 215 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. Bibliography 216 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. Bibliography 217 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."@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "2013-05"@en ; edm:isShownAt "10.14288/1.0071859"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Civil Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "Attribution-NonCommercial-NoDerivatives 4.0 International"@en ; ns0:rightsURI "http://creativecommons.org/licenses/by-nc-nd/4.0/"@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Integrating resin flow and stress development in process modeling of thermoset composites"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/43758"@en .