Analysis of Finite Length, Orthotropic Composite Cylinders Including Through-Thickness Shear Effects Byron A. Mansfield A dissertation submitted to the Faculty of Engineering, University of the Witwatersrand, Johannesburg, in fulfilment of the requirements for the degree of Master of Science in Engineering. Johannesburg, South Africa September 2006 Declaration I declare that this dissertation is my own, unaided work. It is being submitted for the Degree of Master of Science in Engineering in the University of the Witwatersrand, Johannesburg. It has not been submitted before for any degree or examination in any other University. Signed on this the day of B.A. Mansfield i Abstract Thick composite cylinders are important structural elements which cannot be analysed by traditional techniques due to through-thickness effects. This work presents analyses for thick composite tubes of finite length including through-thickness shear. A numerical thermal analysis is implemented for the determination of the transient through-thickness behaviour of tubes. A me- chanical analysis, based on the Rayleigh-Ritz technique, is presented which analyses finite length, composite tubes under a variety of loadings. The anal- yses are shown to be accurate and efficient and are validated against existing results. Results are presented for two ring-stiffened tubes under pressure and thermal loading and also for the transient thermal behaviour of these tubes. It was found that both the through-thickness and transient effects are of im- portance as the stress variation through the thickness and with time was significant for both tubes. ii Contents 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Effects of Geometry . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Mismatch of Poisson?s Ratios . . . . . . . . . . . . . . . . . . 4 1.1.3 Failure of Composite Materials . . . . . . . . . . . . . . . . . 6 1.1.4 Practical Significance . . . . . . . . . . . . . . . . . . . . . . . 8 1.2 Previous Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.1 Analysis Techniques . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.2 Use of Rayleigh-Ritz Technique . . . . . . . . . . . . . . . . . 13 1.2.3 Through-Thickness Shear Effects . . . . . . . . . . . . . . . . 14 1.2.4 Additional Cure Effects in Composites . . . . . . . . . . . . . 16 1.2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.3 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2 Objectives 20 2.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2 Project Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3 Development of Mechanical Analysis Technique 21 3.1 Elasticity Relationships . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2 Strain-Displacement Relationship . . . . . . . . . . . . . . . . . . . . 23 3.3 Selection of Displacement Functions . . . . . . . . . . . . . . . . . . . 24 3.3.1 Normalisation of Displacement Functions . . . . . . . . . . . . 25 3.4 Application of the Rayleigh-Ritz Technique . . . . . . . . . . . . . . . 27 iii 3.5 Implementation of Constraints . . . . . . . . . . . . . . . . . . . . . . 30 3.5.1 Constraint Matrix Approach . . . . . . . . . . . . . . . . . . . 31 3.6 Formulation of Constraints . . . . . . . . . . . . . . . . . . . . . . . . 33 3.6.1 Symmetry Conditions at the Tube Midplane . . . . . . . . . . 33 3.6.2 Displacement Compatibility between Axial Regions . . . . . . 35 3.6.3 Displacement Compatibility between Radial Regions . . . . . 37 3.7 Formulation of Pressure Loading . . . . . . . . . . . . . . . . . . . . . 39 3.8 Determination of Material Properties . . . . . . . . . . . . . . . . . . 40 3.9 Notes on Parametric Analyses . . . . . . . . . . . . . . . . . . . . . . 43 4 Development of Thermal Analysis Technique 44 4.1 Governing Heat Transfer Processes . . . . . . . . . . . . . . . . . . . 44 4.1.1 Conduction Heat Transfer . . . . . . . . . . . . . . . . . . . . 45 4.1.2 Convection Heat Transfer . . . . . . . . . . . . . . . . . . . . 46 4.2 Resistor-Capacitor Formulation . . . . . . . . . . . . . . . . . . . . . 46 5 Validation of Analysis Techniques 53 5.1 Description of FEM Models . . . . . . . . . . . . . . . . . . . . . . . 53 5.2 Thermal Model Validation - Heat Transfer of an Orthotropic Tube . . 54 5.3 Mechanical Model Validation - Orthotropic Tube with Temperature Profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5.3.1 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5.4 Rayleigh-Ritz Analysis with the Constant Axial Strain Assumption . 60 5.5 Results for a Finite Length, Orthotropic Tube with Regard to the Equilibrium Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.5.1 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 iv 5.6 Thick-Walled Isotropic Pressure Vessel Calculations . . . . . . . . . . 65 5.6.1 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.7 Free End Stresses Attained in an Unrestrained Multi-Material Tube under Internal Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.7.1 Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.8 Convergence Assessment of the Rayleigh-Ritz and FEM approaches . 71 5.9 On The Use of Rayleigh-Ritz vs. FEM . . . . . . . . . . . . . . . . . 75 6 Results 78 6.1 Analysis of Ring Stiffened Cylinders . . . . . . . . . . . . . . . . . . . 78 6.2 Predicting the Response of Ring Stiffeners . . . . . . . . . . . . . . . 78 6.3 Modelling Ring Stiffened Cylinders . . . . . . . . . . . . . . . . . . . 79 6.4 Material Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 6.5 Tube Geometric Properties and Loading Conditions . . . . . . . . . . 82 6.5.1 Filament Wound (FW) Tube . . . . . . . . . . . . . . . . . . . 83 6.5.2 Chopped Strand Mat (CSM) Tube . . . . . . . . . . . . . . . 83 6.6 Results for the FW Tube . . . . . . . . . . . . . . . . . . . . . . . . . 84 6.6.1 Thermal Response of FW Tube . . . . . . . . . . . . . . . . . 84 6.6.2 Mechanical Results of FW Tube . . . . . . . . . . . . . . . . . 85 6.7 Results for the CSM Tube . . . . . . . . . . . . . . . . . . . . . . . . 92 6.7.1 Thermal Response of CSM Tube . . . . . . . . . . . . . . . . 92 6.7.2 Mechanical Results of CSM Tube . . . . . . . . . . . . . . . . 93 6.8 Transient Thermal Stresses in Composite Tubes . . . . . . . . . . . . 101 6.8.1 Tube Geometric Parameters . . . . . . . . . . . . . . . . . . . 101 6.8.2 Loading Parameters . . . . . . . . . . . . . . . . . . . . . . . 102 v 6.8.3 Transient Stress Variation Results . . . . . . . . . . . . . . . . 102 7 Discussion 107 7.1 Ring Stiffened Tubes - Thermal Results . . . . . . . . . . . . . . . . . 107 7.2 Ring Stiffened Tubes - Mechanical Results . . . . . . . . . . . . . . . 107 7.3 Stress Discontinuities . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 7.4 Influence of Stress Discontinuities on other Stress Components . . . . 110 7.5 Transient Effects in Composite Tubes . . . . . . . . . . . . . . . . . . 111 7.6 Enforcing Stress Equilibrium Conditions by Means of Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.7 The Possible Benefits of this Analysis for Isotropic Systems . . . . . . 113 8 Conclusions 114 8.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 8.2 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 8.3 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 A Detail Derivation of Mechanical Model 122 A.1 Formulation of Strain Relationships . . . . . . . . . . . . . . . . . . . 122 A.2 Formulation of Energy Potential Relationship . . . . . . . . . . . . . 123 A.3 Formulation of the Loading matrix . . . . . . . . . . . . . . . . . . . 124 A.4 Formulation of Stiffness Matrix . . . . . . . . . . . . . . . . . . . . . 126 A.4.1 Integrations of Stiffness Matrix . . . . . . . . . . . . . . . . . 127 B Estimation of Secondary Poisson?s Ratio, ?23 130 C Sample Calc for Woven Roving 132 vi C.1 Uni-directional Properties for the Woven Roving Laminate . . . . . . 132 C.2 Woven Roving Laminate Properties . . . . . . . . . . . . . . . . . . . 134 vii List of Figures 1 -1 Deformations in a Curved Vessel Associated with a Change in Thickness 3 1 -2 Inter-Laminar Shear arising from Thermal Expansion of a Composite Laminate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1 -3 Manifestation of Inter-Laminar Shear Stress in Laminate . . . . . . . 6 1 -4 Free-Body of the Top Layer in Equilibrium Showing Inter-Laminar Stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1 -5 Sketches of Composite Failure Mechanisms . . . . . . . . . . . . . . . 9 1 -6 Springforward Effect in a Channel Section . . . . . . . . . . . . . . . 17 3 -1 Cylindrical Coordinate System with Associated Displacements . . . . 22 3 -2 Variations Unsuitable for Polynomial Approximations . . . . . . . . . 26 3 -3 Tube Cross-Section with Associated Dimensions . . . . . . . . . . . . 27 3 -4 Tube with Discrete External Pressure . . . . . . . . . . . . . . . . . . 39 3 -5 Laminate Coordinate System . . . . . . . . . . . . . . . . . . . . . . 42 4 -1 Discretization Scheme for Thermal Analysis . . . . . . . . . . . . . . 47 4 -2 Idealised Element System for Thermal Analysis . . . . . . . . . . . . 48 5 -1 Sketch of Loading Conditions on Orthotropic Tube . . . . . . . . . . 55 5 -2 Temperature Distribution Through the Thickness of an Orthotropic Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5 -3 Orthotropic Tube Stress State (?t = 0.25) . . . . . . . . . . . . . . . . 57 5 -4 Orthotropic Tube Stress State (?t = 0.5) . . . . . . . . . . . . . . . . . 57 5 -5 Orthotropic Tube Stress State (?t = 10) . . . . . . . . . . . . . . . . . 58 5 -6 Transient Through-Thickness Radial Displacements of Orthotropic Tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 5 -7 Comparison between Constant Strain and Finite Length Models at the Tube Midplane . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 viii 5 -8 Displacement Variation Along Tube Length (R = 0.5) . . . . . . . . . 63 5 -9 Axial Variation of Stress (R = 0.5) . . . . . . . . . . . . . . . . . . . 63 5 -10 Radial Variation of Stress at Free End (z = L) . . . . . . . . . . . . . 64 5 -11 Isotropic Cylinder With Internal Pressure . . . . . . . . . . . . . . . 65 5 -12 Axial Stress Variation in Multi-Layer Tube under Internal Pressure . 68 5 -13 Axial Stress Variation in Multi-Layer Tube under Internal Pressure (Detail 1) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5 -14 Axial Stress Variation in Multi-Layer Tube under Internal Pressure (Detail 2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5 -15 Displacement Variation Along Length of Tube at the Material Inter- face, (r = 0.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 5 -16 Displacement Variation Through the Tube Thickness at the Free End, (z = L) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5 -17 Stress Residuals for the Axial Stress at the Tube Free End . . . . . . 72 5 -18 Stress Residuals for the Shear Stress at the Tube Free End . . . . . . 73 5 -19 Stress Residuals for the Radial Stress at the Tube Inner Wall . . . . . 73 5 -20 Stress Residuals for the Shear Stress at the Tube Inner Wall . . . . . 74 5 -21 Stress Residuals for the Radial Stress at the Tube Outer Wall . . . . 74 5 -22 Stress Residuals for the Shear Stress at the Tube Outer Wall . . . . . 75 6 -1 Long Tube with Supporting Frames at Regular Intervals Along the Length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 6 -2 Transformation from Representative Region into Analysis Region . . 80 6 -3 Schematic of Ring Stiffened Tube . . . . . . . . . . . . . . . . . . . . 82 6 -4 Make-Up of Filament Wound Tube . . . . . . . . . . . . . . . . . . . 83 6 -5 Make-Up of Chopped Strand Mat Tube . . . . . . . . . . . . . . . . . 83 6 -6 Through Thickness Temperature Distribution - FW Tube . . . . . . . 84 6 -7 Displacement Variation Along Tube Length - FW Tube (r = 104mm) 86 ix 6 -8 Displacement Variation Along Tube Length - FW Tube, Detail 1 (r = 104mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6 -9 Displacement Variation Along Tube Length - FW Tube, Detail 2 (r = 104mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 6 -10 Stress Variation Along Tube Length - FW Tube, (r = 104mm) . . . . 88 6 -11 Stress Variation Along Tube Length - FW Tube, Detail 1 (r = 104mm) 89 6 -12 Stress Variation Along Tube Length - FW Tube, Detail 2 (r = 104mm) 89 6 -13 Stress Variation Through the Tube Thickness - FW Tube (z = 495mm) 91 6 -14 Stress Variation Through the Tube Thickness - FW Tube, Showing only Radial and Shear Stress, (z = 495mm) . . . . . . . . . . . . . . 91 6 -15 Strain Variation Through the Tube Thickness - FW Tube, (z = 495mm) 92 6 -16 Through Thickness Temperature Distribution - CSM Tube . . . . . . 93 6 -17 Displacement Variation Along Tube Length-CSM Tube,(r = 104.5mm) 94 6 -18 Displacement Variation Along Tube Length - CSM Tube, Detail 1 . . 94 6 -19 Displacement Variation Along Tube Length - CSM Tube, Detail 2 . . 95 6 -20 Stress Variation Along Tube Length - CSM Tube, (r = 100mm) . . . 96 6 -21 Stress Variation Along Tube Length at the Outer Material Interface - CSM Tube, (r = 104.5mm) . . . . . . . . . . . . . . . . . . . . . . . 97 6 -22 Stress Variation Along Tube Length at the Outer Wall - CSM Tube, (r = 108mm) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6 -23 Stress Variation Along Tube Length - CSM Tube, Detail 1 (r = 108mm) 98 6 -24 Stress Variation Through the Thickness - CSM Tube, (z = 495mm) . 99 6 -25 Stress Variation Through the Thickness - CSM Tube, Showing only Radial and Shear Stress, (z = 495mm) . . . . . . . . . . . . . . . . . 100 6 -26 Stress Variation Through the Thickness - CSM Tube, (z = 0mm) . . 100 6 -27 Stress Variation Through the Thickness - CSM Tube, Showing only Radial and Shear Stress, (z = 0mm) . . . . . . . . . . . . . . . . . . 101 x 6 -28 Transient Radial Stress Variation at the Axial Midplane for the CSM and FW Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6 -29 Transient Radial Stress Variation Along the Tube Length at the Outer Material Interface, FW Tube . . . . . . . . . . . . . . . . . . . . . . . 104 6 -30 Transient Radial Stress Variation Along the Tube Length at the Outer Material Interface, CSM Tube . . . . . . . . . . . . . . . . . . . . . . 104 6 -31 Transient Axial Stress Variation at the Axial Midplane for the CSM and FW Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6 -32 Transient Circumferential Stress Variation at the Axial Midplane for the CSM and FW Tubes . . . . . . . . . . . . . . . . . . . . . . . . . 106 B -1 Cross Section of a Composite Laminate . . . . . . . . . . . . . . . . . 130 List of Tables 5 -1 Material Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 5 -2 Material Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 5 -3 0? Material Layer Properties . . . . . . . . . . . . . . . . . . . . . . . 66 6 -1 Material Properties for Ring Stiffened Cylinder Simulation . . . . . . 81 C -1 Laminate Constituent Material Properties . . . . . . . . . . . . . . . 132 xi Nomenclature ? Coefficient of thermal expansion, (?C)?1 ? Strain Tensor ? Normal strain component ? Shear strain component ? Thermal diffusivity ? Constraint Forces ? Potential Energy, J ? Density, kg.m?3 ? Stress Tensor, MPa ? Normal stress component, MPa ? Shear stress componenet, MPa c Specific heat, J.kg?1 c Constraint matrix E Young?s modulus, GPa hc Heat transfer coefficient, W/m2K k Conduction Coefficient, W/m ?K kst Stefan-Boltzmann coefficient K Stiffness matrix ` Axial length, mm L Loading matrix p Applied pressure, MPa q Heat flux, W.m?2 Q Heat transfer, J.s?1 r, ?, z Coordinates in cylindrical reference frame ?R, ??, ?Z Surface forces in cylindrical coordinates R,?, Z Body Forces S Compliance Matrix t time, s xii T Temperature, ?C Tr Transformation matrix u, v, w Displacements in cylindrical coordinates, mm U,? Internal Energy, J ? External Work, J ?e External work, J V Volume, m3 W Displacement coefficients ?W Work rate, J.s?1 xiii 1 Introduction 1.1 Background The development of composite materialsi presents the engineer with the ability to tailor the material to meet exacting requirements. Coupled with the fact that com- posite materials tend to possess high stiffness-to-weight and strength-to-weight ra- tios when compared to metallic structures, they have the potential for the design of more efficient structures. To date however, composite materials are still found in limited use, primarily in the high performance industries (automotive, aeronautical and marine)1?3. This is due to a lack of knowledge and full understanding of the mechanisms governing composite materials. This is extremely true when dealing with thick composite parts as the majority of applications of composite materials are with thin laminates (and hence, the majority of analysis is centered around this). In metallic structures, increases in thickness do not generally possess any serious implications, apart from part cost. In contrast, the increase in thickness of composite parts has serious implications which need to be considered. To compound this, these effects are not accounted for by the traditional methods of analysis for anisotropic materials. Traditionally, composite analysis techniques tend to be based on the assumption of plane stress, and are therefore only applicable to thin laminates. Possibly the most well known of these techniques is Classical Lamination Plate Theory (CLPT)4. This theory is based upon the Kirchoff-Love5 assumptions for plate bending: ? Normals to the midplane remain straight and normal to the deformed midplane after deformation. ? Normals to the midplane do not change. By virtue of these assumptions, the through-thickness shear strains and normal stress must be zero. The implication is, therefore, that this theory is only applicable to thin laminates. Provided that this is the case, the results obtained from CLPT are accurate, which has led to its widespread use. In a large number of applications, this condition is sufficient (typically in the aero- nautical industry where very thin materials are utilised). For laminates such as these it can be assumed, in much the same way for isotropic structures, that a state of iThe term ?composite material? is very general and refers to a material which is comprised of multiple constituents. In the literature, materials such as metal-matrix materials are often referred to as composite materials. Following the loose definition above, even steel reinforced concrete can be described as a composite material. For the purposes of this work it is pertinent to limit the definition of composite materials to fibre reinforced plastics. 1 plane stress exists. However, as the thickness of the laminate increases, the through- thickness stresses become more prevalent and CLPT becomes unable to accurately analyse the laminate. Additional techniques must then be employed. Of course, this behaviour, which occurs with increasing thicknesses, is not peculiar only to composite materials. An example where such a situation occurs frequently in isotropic structures is the analysis of pressure vessels 6. Provided the vessel thickness is small compared with the radius, it is safe to utilise thin-walled pressure vessel theory (which is based upon a state of membrane stress only). However, as the vessel thickness increases the assumptions in the thin-walled theory no longer hold, and it becomes necessary to utilise another analysis technique (in this case, the Lame? equations6). However, this behaviour is of even greater importance in composite materials due to their anisotropic nature 2,7,8. While it is good design practice to avoid through-thickness loading when using com- posite materials, due to secondary effects which can contribute toward this loading, they can not be completely avoided. And while these effects are often not considered in the initial design, they can be of significant magnitude to promote failure2,7,9,10. There are two main secondary effects which cause through-thickness stresses, namely, geometry effects and effects of anisotropy. The effects of anisotropy are commonly referred to as the mismatch of Poisson?s ratios. These mechanisms are discussed in further detail in ?1.1.1 and ?1.1.2. 1.1.1 Effects of Geometry One of the problems associated with curved geometry is the coupling effect between the displacement components. In some cases this can lead to a self-imposed restraint on thermal expansion which in turn leads to the formation of a stress state. This effect is discussed in this section. Consider the thermal expansion of a flat, homogeneous, orthotropic plate. In the micro and meso scales a distinction is made between the fibre and the matrix. As these elements are constituted of different materials (the fibre and matrix) they will possess different coefficients of thermal expansion (CTE?s). As there is some chemical bond between the matrix and fibre (the adhesive effect which allows the composite material to function) it is not possible for the fibres and matrix to expand freelyii. There will therefore be a localised stress in the region of the fibre arising from the restraint of the desired thermal expansion of both components. In the macro analysis there is no distinction between matrix and fibre; therefore the material is treated as homogeneous. The aforementioned local effects are therefore not considered in this analysis. With this in mind, in conjunction with the fact that the plate is not subject to any form of external restraint, it can be concluded that thermal expansion of the plate will result in a null global stress state. iiassuming that the fibre and matrix remain perfectly bonded throughout the process 2 A very different scenario is created by the thermal expansion of a cylindrical com- posite vessel though. Initially, consider only the through-thickness expansion of this vessel. As the thickness expands, a change is induced in the curvature of the vessel. This change in curvature manifests itself as an in-plane strain; simply, the change in curvature has resulted in a change in circumference of the vessel. This effect is shown in Fig. 1 -1 and is an illustration of the coupling between displacement components within the curved vessel. For an isotropic vessel this is not a problem as the corresponding in-plane strains due to thermal expansion will compensate the induced in-plane strain such that it is negatediii. For an anisotropic vessel however, the in-plane expansions will not complement this induced in-plane strain. To proceed with the discussion it is assumed that the in- plane properties of the vessel are superior to the through-thickness properties (as is frequently the case). Thus the attempt at modifying the curvature of the vessel by the through-thickness strains, will therefore be, to a large part, resisted, as the in-plane stiffness is larger than that out-of-plane. Evidently, a through thickness stress state will then be formed in the vessel as a result of a self-imposed restraint. In the above argument, as the thickness of the part increases with respect to the nominal part dimensions, this effect becomes more pronounced and consequently, the stresses become more significant. This begins to illustrate some of the possible problems with producing thick-walled composite parts. Figure 1 -1: Deformations in a Curved Vessel Associated with a Change in Thickness iiiprimarily as a result of the isotropy. 3 1.1.2 Mismatch of Poisson?s Ratios In the preceding section the only mechanism whereby thermal strains resulted in a stress state was the coupling effects associated with geometry. However, there are other causes of such stresses; one particular mechanism arises from the differ- ing anisotropic properties between layers within a laminate. This effect is often termed the mismatch of Poisson?s ratios, however, effects such as differing CTE?s can produce similar results. The mechanisms by which the through-thickness stress state is created can be readily derived from simple equilibrium considerations. Consider the laminated flat plate shown in Fig. 1 -2 with a layup of 0?/90?iv. The orientation of the local plate coordinate system is defined in Fig. 1 -2 and the plate possesses a dimension of h in z direction and w in the y direction while possessing unit length in the x direction. It is assumed that the thickness of the individual plies are equal to half the total thickness, h/2. The coordinate system is located such that the width has a range of [?w/2;w/2] and the height, [?h/2;h/2]. The mechanism by which these stresses are generated is most easily described by considering a constant uni-directional applied strain, ?x, as indicated in Fig. 1 -2. Figure 1 -2: Inter-Laminar Shear arising from Thermal Expansion of a Composite Laminate ivWhile this investigation is conducted with composite materials in mind, it is interesting to note that this situation is exactly analogous to the thermal expansion of a bi-metallic strip. This is, of course, a problem which has been studied intensely. 4 Initially, consider the result if the plies within the laminate are not connected, as shown in Fig. 1 -2(a). As the layers are free, they will assume their desired displaced shapes independently. In this case of the 0? layer, the strain is applied in the primary fibre direction. Thus, as this layer is extremely stiff in this direction, the transverse Poisson contraction will be fairly large. In contrast, the strain is applied in the secondary fibre direction of the 90? layer. The associated transverse contraction of this this layer will be minimal. Due to the different anisotropic properties in the layers, the transverse contractions will differ in magnitude considerably. By considering the fact that the layers are bonded together, it is obvious that the laminate must attain some compatibility condition for the transverse contraction. Consequently, a transverse stress, ?y, must be present to realise this situation as can be seen in Fig. 1 -2(b). This stress can be characterised by some simple equilibrium considerations. For equilibrium throughout the laminate it is required that the resultant stress in each layer be equal in magnitude and opposite in sign. This can be expressed as: ? h/2 ?h/2 ?y dz = 0 (i) Additionally, the stress at the edge of the laminate must be zero as a consequence of the boundary conditions. The fact that this stress varies with the y coordinate implies that the shear stress, ?yz, is present. These stresses are sketched in Fig. 1 -2 and Fig. 1 -3. The existence of the this shear stress can also be confirmed by considering the mechanism by which the stress ?y is transmitted between the two layers. i.e. the only means by which the laminate can be forced to yield the displaced shape shown in Fig. 1 -2(b) is by means of shear stress along the interface. The inter-laminar stress state is examined more closely in Fig. 1 -3. By considering the free body of only the top layer it can be seen that the stress resultants of ?yz and ?y do not share a line of action. Of course in Fig. 1 -3 the stresses are only idealised as resultants without any form of calculation, but by considering the stresses in their general form, the previous statement will hold. Therefore an additional stress must be present to preserve the moment balance. This stress manifests as a interlaminar normal stress, ?z as can be seen in Fig. 1 -4. As can be seen from Fig. 1 -4 that the layer will be in equilibrium provided the following relationships hold: ? h/2 0 ?y dz = ? w/2 ?w/2 ?yz dy (ii) ? h/2 0 ?yz dz = ? w/2 ?w/2 ?zy dy (iii) 5 Figure 1 -3: Manifestation of Inter-Laminar Shear Stress in Laminate Figure 1 -4: Free-Body of the Top Layer in Equilibrium Showing Inter-Laminar Stresses One of the elementary results is that at some point along the interface this normal stress must reverse it?s sign in order to preserve vertical equilibrium, i.e.: ? w/2 ?w/2 ?z dy = 0 (iv) An important phenomenon arising from the interlaminar stresses discussed above, is that at the material interface at a free-edge, the interlaminar normal stress are predicted as stress singularities11,12. Simply, as the free edges are approached along the material interface, the interlaminar normal stress approaches infinity. 1.1.3 Failure of Composite Materials While ?1.1.1 and ?1.1.2 explain the existence of through-thickness stresses, they do not explain their significance. It may be argued that these are essentially secondary effects and can therefore be neglected in the initial analysis. As stated previously, while this may be safe design practice for thin laminates, it is not conservative design practice for thick laminates10. In order to assess the significance of these effects it is necessary, first, to consider the failure mechanisms for composite materials. Compared to metallic structures, failure of composite materials is not a trivial concept. The actual theories governing failure of composite materials are many and, at the time of writing, there is no single theory which adequately predicts the failure of composite materials for all conditions. Therefore it is at the discretion of the designer as to which theory best describes a given situation. For a comprehensive discussion on these failure theories 6 the reader is directed to Herakovich4. In metallic structures, static failure is governed by the material ultimate strengths and stability failure is governed by the material stiffness. Due to the isotropic nature of metals, both of the aforementioned properties are independent of the material orientation. However, in addition to anisotropic properties, there are a number of additional failure modes which must be accounted for in composite materials. In order to assess the failure of a composite laminate it is necessary first to investigate the failure mechanisms of the composite constituents. The two most obvious failure mechanisms are rupture of the fibre and matrix components, which is a particular concern under tensile loading. As in metals, this failure is governed by the material strengths (in the case of glass fibres this will be a constant value as the fibres are isotropic, while carbon fibres exhibit anisotropic strengths). While it is possible for the matrix to fail in compression, the fibres under compression are generally dominated by stability failure, often termed microbuckling (illustrated in Fig. 1 -5). The other failure mechanisms are dependent upon the adhesive bond between the matrix and the fibre. When this adhesive bond is overcome, the failure is termed debonding. When debonding is combined with rupture of the fibres, the failure is termed fibre pull-outv. Both of these failure mechanisms are also illustrated in Fig. 1 -5. An important aspect in composite material failure is the contributions of the con- stituent components to the final failure of the laminate. The main purpose of this discussion is to investigate which failures are predominantly due to the fibre (fibre dominated) and those which are governed by the matrix (matrix dominated). For the purposes of this discussion it is assumed that debonding and pull-out failures do not occur as both these failures are governed by the adhesive strength between the resin and fibre. Consider in-plane tension loading on a composite laminate. According to the afore- mentioned failure modes, the laminate is most likely to fail due to rupture of either the fibres or the matrix. Thus the laminate strength may be governed by either the fibre or matrix ultimate tensile strength. As the fibres are the main structural components in the laminate, good design practice would dictate that the desired failure mode would be fibre rupture. In the majority of cases, except where the fibre volume fraction is extremely low, this is typically the case. If the in-plane loading is compressive rather than tensile, the failure will be dominated by microbuckling of vNote that as composite materials generally consist of a laminate of many plies, failure in a ply may not lead to catastrophic failure of the entire laminate. The definition of catastrophic failure of a laminate therefore is not altogether clear as with metals. Some design codes specify that laminate failure has been deemed to occur when the first ply fails (FPF). Others are less conservative and allow failure to be determined by last ply failure (LPF). It is at the discretion of the designer to select the appropriate failure code. 7 the fibresvi. Therefore the in-plane properties are predominantly fibre dominated. Again, this makes sense as the fibres are typically the main structural constituents. As there are no fibres in the through-thickness direction, the main structural com- ponent in the laminate in this direction is the matrix. Therefore the out-of-plane properties, and in turn out-of-plane failures, are dominated by the matrix. Note, that as the matrix mechanical properties are not high, debonding failures can be prevalent between plies. This effect is known as delamination and is a particularly important failure mechanism. Comparing the matrix and fibre properties it can be seen that the fibre mechani- cal properties are far superior to those of the matrix. Therefore, the out-of-plane strength is typically far lower than that exhibited in-planevii. This low through thickness strength can now be evaluated in conjunction with the stresses arising from secondary effects. The obvious conclusion is that although the secondary stresses may be of lower magnitude than the design stresses, they can still give rise to a critical situation in light of the lower through-thickness strength. They cannot, therefore, simply be considered as non-critical and neglected. However, in typical design scenarios this is still commonly the case. 1.1.4 Practical Significance Up to this point it has been established that there are problems associated with thick composite laminates arising from secondary effects. The question must however be asked, of what practical significance are these problems? This can be answered firstly by considering the significance these stresses have on material failure, and secondly, by examining where such situations are likely to be found in reality. In addition to constant elastic properties, isotropic materials are assumed to possess material strengths which are independent of the material orientationviii. Thus, as secondary loading is typically of smaller magnitude than the design loading, it is frequently ignored. To some point, this procedure is safe as the stresses generated from secondary effects are typically smaller than the design loads, and the material is more than capable of supporting them. viAs the fibres are essentially extremely slender columns, they are highly prone to buckling under compression. And while they are supported laterally by the matrix, the buckling loads are typically still far below the resin compressive strengths. viiAs an aside, there is much research devoted to improving this situation by various means. In some cases, fibres are stitched through the laminate to attempt to include some fibres to carry through thickness loading13,14 thereby increasing the through-thickness strength. viiiNote that due to manufacturing processes (rolling, forging, etc) metals tend to possess a slight grain structure. However, the anisotropy due to this grain structure tends to be minimal and they can, for the most part, be safely treated as isotropic. 8 Figure 1 -5: Sketches of Composite Failure Mechanisms As discussed in the previous section, composite materials exhibit vastly different material strengths in different directions. For some laminates the through thickness strength can be 10% of that in-plane8. Therefore through thickness stresses which are of an order of magnitude smaller than the in-plane stresses may still be sufficient to promote failure in the material. Real world scenarios in which the aforementioned stresses can be attained are dis- cussed below. The stresses can be attained in two general cases: manufacture and in-service use. The stresses induced in manufacture occur primarily during the post-cure operation. During the normal curing process for a composite part, the gelation and polymeri- sation of the resin is not able to complete sufficiently. Thus it is common practice to subject a part to an elevated temperature post-cureix. This post-cure operation allows these processes to progress to a more acceptable state and thereby, improves the mechanical properties of the part. During this process it is widely assumed that the part attains a stress free state while at the elevated temperature. Now, upon cooling, residual stresses are created by a combination of two effects: thermal con- traction and resin shrinkage. The thermal contraction is obvious as the part cools to room temperature from the elevated temperature, stress free state. However, the ixThe exact details of the process are highly dependent on the nature of the part and the material makeup and are usually conducted along the material manufacturer?s guidelines. 9 elevated temperature allows the cross-linking of the resin to progress, thereby the resin exhibits a form of shrinkage due to the increased polymerisation. This effect also adds to the loading experienced by the part. In the case of tubes, these stresses can be of significant magnitude to promote failure in the part 7,15. Apart from the advantageous mechanical properties exhibited by composite ma- terials there are also advantages in other areas, for example, the ability of glass reinforced plastics (GRP) to resist chemical attack 3,7. To fully assess the signifi- cance of this attribute of GRP, it must be remembered that this is an area where metallic structures are particularly susceptible. This makes GRP a primary candi- date for use in transporting chemically aggressive fluid media (for example, in pipe networks and pressure vessels). Typically however, the mechanical properties of GRP are generally inferior than steel traditionally utilised for piping networks. In order to account for this, GRP pipes are comparatively thicker than their metallic counterparts. Apart from the manufacture stresses described above, through-thickness stresses can also be obtained from in-service loading. As an example, it is plausible to envisage the above pipes transporting fluid which is both at an elevated temperature while being pressurized. Thus the tubes will be subject to both thermal and mechanical loading. This loading, coupled with the thick-walled nature of the vessel, results in stresses via both the mechanisms of ?1.1.1 and ?1.1.2 being induced. It is highly likely that composite cylinders will be subject to both forms of the above loadings, i.e. manufacture and in-service. It has also been shown that these stresses are a significant factor in the failure of the vessels. These through-thickness effects are therefore of prime importance in the design and analysis of such vessels. 10 1.2 Previous Research In the preceding section it was shown that the likelihood of composite vessels ex- periencing through-thickness stresses from cure and thermal effects was high. The importance of these effects was also established; due to this, these effects have been the subject of much investigation. Some relevant studies are assessed in this section. 1.2.1 Analysis Techniques In general, mechanical analyses are governed by the theory of elasticity. For an in depth and detailed development of this theory, the reader is directed to the classical work by Love5. For a more practical approach, Timoshenko and Goodier 16 apply the theory of elasticity to solve a number of problems. The theory of elasticity in essence, describes the relationships which are required in order to attain equilibrium for a given solid body under the action of body and surface stress. The relationships are written in differential form and the solutions are attained by solving the governing differential equations. Elasticity solutions are useful as they predict the complete stress state throughout a body including the boundary effects. The results of the elasticity solutions are also exactx. In selected cases, the solution to the differential equations are easily attain- able. However, in most cases, the solutions are mathematically complex, usually involving special functions (Bessel functions for example). These special functions are often written in the form of an infinite series, thus they are not particularly easy to utilise. But aside from the mathematical complexity, possibly the largest disadvantage of elasticity solutions is their sensitivity to the boundary conditions of the problem. Even small changes in the boundary conditions require the problem to be refor- mulated. Additionally, the form of the solution can change radically with different boundary conditions. Thus elasticity solutions are usually limited to the selected analysis case. They are not, therefore, useful in analyses where the boundary condi- tions are expected to change frequently, and/or different load cases are to be anal- ysed. In situations such as these, the analysis is typically conducted by an adapted applied mechanics technique. A number of candidate techniques are discussed below. A popular technique for analysing cylinders (or any vessels possessing curvature) is by means of shell analysis17,18. As implied in the title, this technique assumes that the structure is a shell, hence the thickness is small, and through-thickness effects can be neglected. The advantage of this technique is the ease with which it can be applied to a large variety of geometries. The technique is also established and has been well developed19,20. It is often utilised for stability analyses. xAs with most engineering cases, the term exact should actually read: exact within the limits of the underlying analysis assumptions. 11 As the technique is primarily for use in thin structures, its ability to predict be- haviour in thick structures has to be proven. Although in some cases, some success was obtained by employing this technique for thicker structures (for example Chan- drashekhara and Gopalakrishnan21, where some bounds of limitations of shell analy- ses are presented), it is not particularly suitable for analysing thick structures17, and thus alternative techniques should be employed. Some alternatives are numerical solution approaches and elasticity solutions. Examples of each are discussed below. A prediction of the residual stresses obtained during the post-cure process on com- posite tubes was performed by Stone et al.7. In this work, the problem of post-cure was analysed in the form of a thermal expansion problem by replacing the thermal strain terms, ??T , by material shrinkages. These material shrinkages essentially characterise the laminate strains which would be expected upon the unrestrained post-cure of the laminate, were it flat. The stresses therefore arise from the restraint which is placed upon the tube which prevent the mechanical strains from attaining the desired free shrinkage values. The results of this work were attained numerically by use of the Rayleigh-Ritz technique. It was shown that this technique was indeed a viable one to predict these stresses. Due to the assumptions made during the analysis formulation, the model it is limited to predicting the stresses in long tubes at points away from the ends. This arises from the assumption that the axial strain is constant and independent of the axial coordinate. In addition to this assumption, axisymmetry was also assumed. The conclusion of this work showed that the stresses obtained during post-cure were of sufficient magnitude to warrant consideration. A particular area of concern was the presence of tensile stresses upon the inner wall of the tube. If the tube is transporting chemically aggressive media, the combination of the chemical attack and the tensile stresses could promote Environmentally Assisted Cracking (EAC) and hence, failure would be attained at levels far below what would be expected from the mechanical loading7. As the analysis of Stone et al.7 replaces the thermal strains by material shrinkages, the formulation is exactly the same as that of a thermal expansion problem with constant through-thickness temperature. It can however be expected that for a large majority of problems, the temperature profile through the tube wall will not be constant. This is especially true if the transient thermal response of a tube is considered. The response of an orthotropic tube to transient thermal behaviour was presented by Kardomateas8. As transient thermal behaviour is considered by Kardomateas8, allowance is made for a non-uniform temperature profile through the thickness of the tube. However, in a similar manner to Stone et al.7, the axial strain is assumed to be constant which follows from the assumption of a long tube, and axisymmetry is also assumed. In the work of Kardomateas8 both the thermal and the mechanical problems were solved by analytical means, with the resulting models being valid for a homoge- neous, orthotropic tube. The thermal loading is considered as a tube with fixed inner wall temperature with a free convection condition on the outer wall. The 12 model presented is a linear, uncoupled model and, as such, temperature and stress independent properties were assumed throughout. While the work of Kardomateas8 is concerned more with the thermal expansion problem than post-cure, the results presented indicate that stresses of significant or- der were also attained by this mechanism. Additionally, it is important to note that the results presented indicate that the stress states for a non-uniform temperature distribution can be of higher magnitude than those for a uniform through-thickness temperature. The major limitation of both of the above works lies in the assumption of constant axial strain. This assumption has the implication that both models are only capable of predicting the stress state far away from the ends in long tubes. In practice, not all tubes can be considered long and some account must be made for the length of the tube. Additionally, the effects of through-thickness shear can be important in thick tubes; neither of the above works (Stone et al. 7 and Kardomateas8) can account for these effects. 1.2.2 Use of Rayleigh-Ritz Technique As the Rayleigh-Ritz technique is hardly a new development, an analysis as to its suitability for problems of the type in this work can be accomplished by examining previous work of a similar nature in which the technique was utilised. A collection of such cases are discussed below. An examination of the literature yields many cases in which the Rayleigh-Ritz tech- nique has been utilised to predict the response of composite materials. A large num- ber of cases utilise the Rayleigh-Ritz technique to describe the vibrational problem of a system. However the basic principles behind the technique remain the same and there are only minor modifications to the problem formulation between the vibrational and mechanical systems. An example of such an analysis is presented by Dasgupta and Huang 1. In this work a combined Rayleigh-Ritz or Finite Element Method (FEM) hybrid model is utilised to predict the vibration of thick laminated composite spherical panels. Use is made within this work of the layer-wise approach. In this way, each layer of composite material within the laminate is treated as homogeneous. This effectively allows each layer within the laminate to be treated separately, i.e. possessing its own stiffness matrix. This approach deals with the composite material without having to explicitly deal with the anisotropic coupling effects. The final result is therefore very convenient requiring little computation in the formulation to account for the effects of anisotropy. Continuing with vibration problems, Ip et al. 22 uses the Rayleigh-Ritz technique to model vibrations in thin orthotropic tubes. In this work it is shown that for the selected cases, the Rayleigh-Ritz results compared well with test results (within 13 3%). The basis of this model is the Love shell theory and is therefore incapable of modelling tubes which are thickxi. With this aside, this work does show the applicability of the Rayleigh-Ritz technique in modelling orthotropic cylinders. In the above works there has been little consideration given to the through-thickness properties of composite materials. However the Rayleigh-Ritz technique has been applied successfully to model such effectsxii. Zhu and Lam23 presented a model based on the Rayleigh-Ritz technique which predicts local stresses in composite plates. The work is primarily concerned with predicting the through-thickness stress state in laminated plates. Although the analysis is conducted for flat plates, it is sufficient to be able to draw conclusions about the analysis technique for curvilinear geometry. This analysis utilises a layer-wise formulation to model the composite materials. The analysis is also formulated in such a manner so as to ensure that the compatibility conditions are explicitly satisfied by the displacement functions. The results are compared to existing results for composite plates and results obtained from a FEM analysis. Good correlation is shown in all cases. A small discrepancy was noted at the edge of a plate when compared to the FEM result although this was a localised effect; for the same plate the results away from this area compared excellently. Thus it was concluded that the approach chosen was successful in possessing the ability to model effects in laminated composite plates. An important conclusion from the work of both Zhu and Lam 23 and Dasgupta and Huang1 was that both works showed clearly that, by utilising the Rayleigh-Ritz approach, a significant reduction in the required Degrees of Freedom (DOF?s) was obtained for a similar accuracy when compared to the FEM method. 1.2.3 Through-Thickness Shear Effects Of the classical assumptions in the analysis of composite materials, one of the most common is the neglection of through-thickness effects, particularity through- thickness shear. However inter-laminar shear is a governing mechanism by which composite materials can fail. Additionally, through-thickness shear effects are im- portant when considering the required equilibrium conditions for finite length tubes. And while, for thin laminates this effect is not dominant, it can be critical in thicker laminates. A variety of methods have been utilised to account for this component to date. Some of the most popular techniques are the shear deformation theories24. These techniques originated with the first order shear deformation theory (FSDT). Subse- quently this theory has been extended to quadratic order and subsequently to higher order (HSDT). There is an abundance of literature for this technique, the reader is referred to Reddy and Wang24, Noor and Burton25 and Kant and Swaminathan26 which present reviews of some literature regarding shear deformation theories. Se- xiIp et al.22 uses a ratio of ri/t ? 20. xiithe work of Stone et al.7 falls under this category and has been discussed earlier. 14 lected implementations of the technique can be seen in Dube et al.27, Ferreira et al.28 and Ganapathi et al.29. A brief outline of this technique is presented below. The approach has been developed as an extension to classical plate theory. The Kirchoff-Love equations for flat plates neglect all forms of through-thickness effects, however, FSDT reformulates the problem by considering the displacement field24: u(x, y, z, t) = z?x(x, y, t) v(x, y, z, t) = z?y(x, y, t) (v) w(x, y, z, t) = w0(x, y, t) Where ?x and ??y are the rotations about the x and y axes. It can be divined from the above relationships that the through-thickness shear strain will be con- stant. Thus it cannot satisfy the equilibrium conditions as, in the absence of body forces, the shear strain on the boundaries of the plate must be zero. In FSDT additional factors are utilised known as shear correction factors. Essentially these are utilised to attempt to improve the approximate results obtained from FSDT. However these factors are highly dependent upon the problem parameters and are not easily determined in a general form. HSDT modifies the above relationships by including higher order terms of the through-thickness displacement expansion (i.e. the coordinate z). In particular, a third order theory is able to satisfy the required equilibrium conditions and it becomes unnecessary to include shear correction factors 24. There is however a sig- nificant increase in computational expense when transferring from FSDT to HSDT though. The principle of the shear deformation theories is to present a simple and compu- tationally efficient method of including through-thickness effects. An alternative is to simply consider the continuum problem, i.e. the fact that the through-thickness shear is directly dependent on the translational displacement field. In Cartesian coordinates this is written as: ?xy = ?u ?y + ?v ?x (vi) Where u, v are the translational displacements related to x, y respectively. Therefore, by simply allowing the displacements to vary with all coordinates, the shear strain can be calculated directly. There is no need then to utilise an as- sumption which relates the displacements to the through-thickness shear as in the shear deformation theories. This simplifies the problem formulation considerably. In essence, this is the technique followed by standard finite element analyses of through-thickness effects. 15 1.2.4 Additional Cure Effects in Composites Whilst most of the discussion thus far has been focused on the stresses obtained from thermal processes, also of importance are the deformations which occur during the curing and post-cure processes of composite parts. Although these effects are not the primary concern of this work, similar mechanisms are responsible for both cure deformations and residual stresses. The analysis techniques applied for determining cure deformations also have the potential to be applied to determine stress states. Additionally, some of the results obtained from analysis of cure deformation could have repercussions for the analysis contained in this work. Therefore it is important to investigate this area of research. Deformations occurring during cure and post-cure increase the difficulty of accu- rately manufacturing composite parts. A good example of such a situation is the manufacture of composite pipe flanges. As these flanges form part of a bolted joint (primarily under tension), they are often loaded through the thickness. This has the effect of producing relatively thick flanges (by virtue of the low through-thickness strength of the laminates). In order for the joint to function optimally, it is required that the mating flanges are flush. Due to the springback/springforward which is evident after curing, the flanges cannot simply be manufactured at right angles to the longitudinal axis of the tube. Therefore some allowance must be made to ac- count for the springback to yield a final product which is perpendicular to the tube. The alternative exists to manufacture the tube and machine the faces of the flange. However, this is not desirable as it incurs additional cost in machining as well as the additional time for manufacture. Importantly, the machining cuts the fibres and exposes them from the laminate, this is not ideal from many points of view (fatigue, environmental degradation, etc.). It is not a trivial matter to predict the springback of such parts, and therefore it becomes difficult to account for this during manufacture. To date, this process is somewhat of a black art and depends to a large degree on the manufacturer?s experience. Therefore some amount of trial and error is required in order to fine tune the manufacture process. This of course has negative connotations on part lead time and cost. Wiersma et al.9 investigate several techniques for predicting springforward effects. In this work a linear stress model is developed and validated. The displacements within this model are assumed to vary only with the radial coordinate of the ?L? shaped curve. As such, it is assumed that the part length is long and edge effects can be ignored. As this model is linear it cannot adequately describe irreversible processes such as those observed during a cure cycle. To account for these effects, the linear model is adapted to a visco-elastic model which incorporates nonlinear effects. The final model is solved by the Finite Element Method. It is shown that the prediction of the springforward effects were improved by including visco-elastic effects. However, some disagreement between the predicted results and measured test data was still found. 16 Interestingly, as this analysis9 assumes that the only variation in displacements occurs along the radial coordinate, the analysis is essentially conducted over a cylin- drical segment. The effect of cure induced deformations was also investigated by Jain and Mai2 for channel-type sections. In this work the analysis was conducted by means of modi- fied shell theory. While this model was shown to be accurate and computationally efficient, it is limited by the underlying assumptions. Throughout the derivation it is assumed that the through-thickness stress is zero. However, research has shown that this stress can be significant in thick vessels 7,8. In addition, the transverse shear strains are neglected. Therefore the model cannot satisfy the required equilib- rium conditions for a finite length and is therefore limited to the treatment of long sections. Again, it was noted by Jain and Mai2 that the majority of the deformation occurs across the corner fillet. The cause of this effect has already been shown in this work (see ?1.1.1). Thus, the situation is one where the straight flanges primarily exhibit rigid body rotation as shown in Fig. 1 -6. In order to account for this, Jain and Mai2 assume that the ?L? shaped section can be represented by a cylindrical section representing the corner fillet. This obviously reduces the amount of analysis required significantly. It is interesting to note that although the sections analysed by Wiersma et al.9 and Jain and Mai2 are not in themselves cylindrical, both analyses simplify the situation until the representative system is a cylindrical section. Thus it can be seen that cylindrical systems are highly important components of use even when the geometry is not purely cylindrical. Figure 1 -6: Springforward Effect in a Channel Section 17 1.2.5 Summary From the works discussed in this section it can be seen that thermal effects in com- posites materials are an area of high interest. However, in terms of thick laminates such as composite tubes, few works are able to consider a finite length composite tube, with most works on cylinders reverting to the assumption of constant axial strain which implies infinite length. Additionally, few works adequately incorporate through-thickness effects. A large portion of analyses are conducted on shell theory, which completely ignores through-thickness effects. Others use modified shell the- ory in an attempt to include through-thickness effects, however, these techniques all have some limitations. 1.3 Motivation The effect of through thickness stress fields is an important aspect in the design of thick composite parts; particularly in thick composite cylinders as they form important structural components30. Previous work7,8,15 has shown that stresses obtained from thermal loading as well as the residual stress obtained from post-cure can be large enough to promote failure in such vessels at levels much lower than would have been predicted in the design stage. The work completed to date tends to focus on tubes with infinite length and/or tubes in which the thickness is not considerable. These assumptions will not be able to cover all cases where this type of analysis is required. Additionally, important effects such as through thickness shear are not accounted for in these analyses. There exists therefore the opportunity to investigate finite length tubes. The validity of this analysis will be improved, and therefore will possess the ability to analyse a larger number of real life tubes without the need to revert to simplifying assumptions. It has been found that many works exist on this subject, however the above limita- tions notwithstanding, they tend to concentrate on different effects. Additionally, a variety of analysis techniques are employed to solve the problems. There is therefore an opportunity to combine these analyses into a single analysis such that various loadings can be analysed easily. While this in itself is not a tremendous benefit, it will of course allow combinations of loadings to be applied in the analysis. Thus the effects of these loading combinations can be analysed with ease. While it to is possible to analyse such cases by means of a commercial FEM software package, this approach has it?s disadvantages. The time associated with the creation, runtime and analysis of such models can be significant. Additionally these packages do not lend themselves to parametric analyses, thus if the parameters of the vessel change it can be expected that the model will have to be modified which also requires significant time. It is proposed that these are conducted in the design stage of such vessels. In this environment it can be expected that the vessel will undergo a number of modifications as the design progresses through iterations. Thus the total time invested in the FEM modelling can become significant. 18 There is therefore merit in employing alternate techniques to analyse these cases. The Rayleigh-Ritz technique shares some similarities with the finite element method in that they are both displacement based numerical techniques. However whereas the FEM is based upon the use of a number of general elements, it is possible to formulate the Rayleigh-Ritz model for a specific problem. Thus to an extent this model can be tailored for the analysis at hand. In this way it should be possible to generate the model so as to require the user to supply a number of standard values in order to perform the analysis. Additionally, past work 1,23 has shown that implementing the Rayleigh-Ritz technique in this way yields a model which is more computationally efficient than the FEM. In this way a model could be developed which would act as a tool in the design of such vessels by allowing these secondary effects to be accounted for. 19 2 Objectives 2.1 Problem Statement Currently the only technique available to analyse finite length composite tubes under general loading is by means of the Finite Element Method (FEM). While the FEM is a versatile technique, the time invested in creating, running and analysing FEM models can be significant. Therefore a need exists for a computationally efficient tool which allows parametric analysis of finite length composite tubes subject to a various loading conditions. 2.2 Project Objectives ? A method of analysis must be created which can accurately predict the tran- sient thermal response of a composite tube. ? An analysis method should be developed which predicts the stress state in a finite length, composite tube. ? In doing so it is necessary for allowance to be made for non-constant axial strain and for the existence of through-thickness shear effects. ? This analysis method must be able to account for any combinations of the following loadings: - Thermal loading from a temperature profile across the tube thickness. - Loading due to pressure (internal and external) over the entire tube length. - Loading due to pressure (internal and external) over a discrete length of the tube. - Loading arising from material shrinkages. ? The effect of the various loadings on composite tubes must be investigated. 20 3 Development of Mechanical Analysis Technique The analysis is based upon the Rayleigh-Ritz technique which is a numerical, vari- ational technique based upon the principle of minimisation of the system potential energy. A brief outline of the method is presented in this section, a more detailed derivation can be found in Appendix A. 3.1 Elasticity Relationships From the classical elasticity relations in conjunction with Hooke?s Law, the stress- strain relationship can be expressed by the constitutive equation: ? = C? (1) Where ? is the stress tensor, C is the material stiffness matrix, and ? is the strain tensor. The geometry of the tube suggests that it is most appropriate to perform the analysis in a cylindrical coordinate system. The exact definition of the coordinate system and associated displacement components is presented in Fig. 3 -1. From the problem statement it can be inferred that all loading and material properties are invariant with the angular coordinate. Thus the problem is completely described by the axial and radial coordinates, and in doing so fulfills the requirements of axisymmetry. This simplifies the necessary manipulations significantly. With this in mind, Eqn. 1 can be expanded for an orthotropic material in cylindrical coordinates to: ? ??? ?zz ??? ?rr ?rz ? ??? = ? ??? C11 C12 C13 0 C21 C22 C23 0 C31 C32 C33 0 0 0 0 C55 ? ??? ? ??? ?zz ? ?z?T ??? ? ???T ?rr ? ?r?T ?rz ? ??? (2) As the problem is assumed to be axisymmetric, it follows that the materials under consideration cannot be fully anisotropic and must be able to be described as or- thotropicxiii. Eqn. 2 has been written for thermal loading; this can of course be extended to other loadings by replacing the ??T terms with the relevant strains. From elasticity relations, the matrix C can be calculated for an orthotropic material xiiiThis statement of course considers that isotropic and transversely isotropic materials are special cases of an orthotropic material, just as an orthotropic material is a special case of an anisotropic material. 21 Figure 3 -1: Cylindrical Coordinate System with Associated Displacements simply by calculating the inverse of the material compliance matrix as: C = S?1 = ? ??? 1/Ez ??z?/Ez ??zr/Ez 0 1/E? ???r/E? 0 1/Er 0 sym 1/Gzr ? ??? ?1 (3) Where S is the material compliance matrix. By considering the equations of motion, it is possible to construct the equilibrium equations of a solid body in terms of the applied stresses. Eqn. 4 represents the equilibrium conditions for a solid in cylindrical coordinates. The complete derivation of the result in Eqn. 4 is presented in Love 5. ??rr ?r + 1 r ??r? ?? + ??rz ?z + ?rr ? ??? r +R = 0 ??rz ?r + 1 r ???z ?? + ??zz ?z + ?rz r + Z = 0 (4) ??r? ?r + 1 r ???? ?? + ???z ?z + 2?r? r + ? = 0 While it is not true for the general case, the body forces are neglected in this analysis. However, due to the axisymmetry assumptions, all loading in the circumferential direction must be zero (i.e. ? = 0). Additionally, the shear stress components ?r? and ??z must be zero as must any terms which are differentiated with respect to ?. Under the assumptions of this work, Eqn. 4 can be condensed to the form in Eqn. 5. ??rr ?r + ??rz ?z + ?rr ? ??? r +R = 0 22 ??rz ?r + ??zz ?z + ?rz r + Z = 0 (5) If the case of the tube without body forces is examined, Eqn. 5 describes the equi- librium conditions which must be satisfied throughout a solid body. At the surfaces, the stress must react the applied surface loads (normally denoted: ?R, ??, ?Z)16. Thus, force equilibrium must be satisfied. If there are no applied stresses, this dictates that the normalxiv surface stresses must be zero. This is logical as it dictates that the stresses must be self equilibrating (free tube). Eqn. 5 also indicates the importance of including the shear stress ?rz in the formulation of the model in order to satisfy the equilibrium conditions. The results presented in this section govern the stress-strain relationship as well as the potential energy and equilibrium of the system. However, as the Rayleigh-Ritz technique is a displacement based solution, it is necessary to relate the strain and the displacement. This relationship is presented in the following section. 3.2 Strain-Displacement Relationship The strain-displacement relationship can be determined by considering the effect of small displacements on a differential element. For the axisymmetric analysis in cylindrical coordinates the result is given in Eqn. 6. A more detailed and involved discussion about these results can be found in Love 5. ?z = ?u ?z , ?? = w r (6) ?r = ?w ?r , ?zr = ?u ?r + ?w ?z Note the use of the partial derivative as the displacements themselves can be written as functions of the spatial coordinates as: u = u(z, r), w = w(z, r) (7) v = 0 It can be seen that the displacements and strains are independent of the circumferen- tial coordinate, ? due to the assumptions of axisymmetry. Simply, the axisymmetry reduces the 3-Dimensional problem into a 2-Dimensional one as there are only two independent coordinates. The strain-displacement relationship is summarised in matrix form in Eqn. 8. xivThis can be claimed due to the orthogonal nature of the coordinate system. 23 ? = B.D (8) B = ? ?? ??z 0 0 ??r 0 1 r ? ?r ? ?z ? ?? T (9) Where D is the matrix of independent displacements, [u w]T . 3.3 Selection of Displacement Functions In order to implement the Rayleigh-Ritz technique it is necessary to describe the displacements as a general analytical function (or series thereof). The choice of the actual displacement functions is actually fairly arbitrary. Obviously though, the closer the displacement functions are able to describe the system behaviour, the better convergence will be and the overall accuracy of the solution will be higher. In the literature many forms of displacement functions have been utilised to varying degrees of success. Perhaps the most common form of the displacement functions is in the form of trigonometric functions. The great advantage of these trigonometric series is the orthogonality of the basis functions. This leads to stiffness matrices which possess terms only on the main diagonal. This is plainly a useful result as the required manipulations of the stiffness matrices become trivial. However, it becomes problematic with series of this form to satisfy arbitrary boundary conditions. In some cases, it has even been shown that literature which has utilised series of this form inherently violated the natural boundary conditions 31. As discussed before, this prevents the solution from satisfying the governing differential equations. Various polynomial forms have been utilised (Legendre 32, orthogonal, etc). These polynomials tend to be formulated such that the boundary conditions are inher- ently satisfied. Whilst this is convenient it does require considerably more effort to generate the functions. Additionally, modifying the boundary conditions require the polynomials to be reformulated thereby implying a reformulation of the entire problem. An alternative technique which avoids the above disadvantages is to apply general polynomial expansions. Due to their general nature these functions do not inherently satisfy any boundary conditions; they can be made to satisfy any required boundary condition as will be discussed later. Thus this technique has been employed for this work. The form of the displacement functions follows. 24 u(z, r) = M? m=0 N? n=0 Wmnzmrn (10) w(z, r) = P? p=0 Q? q=0 Wpqrpzq (11) There are a few disadvantages and/or limitations which do, however, need to be considered with general polynomial expansions. As can be seen from Eqn.?s 10 and 11, as the order of the polynomials increases, so too does the descriptive ability of the displacement function. In theory, as the displacement orders tend toward infinity, the displacement functions would be able to describe the exact displacements within any continuous region. However, the reality is that computational accuracy places limits on the displacement function orders; if these limits are exceeded the solution becomes numerically unstable. The other limitation which must considered is concerned with the ability of the poly- nomial functions to capture the required behaviour. Polynomials are best suited to describing functions which vary smoothly and by relatively small amounts. Con- versely, polynomials find difficulty in describing functions which possess discontinu- ous behaviour, large change over a small local area and sudden changes in curvature. Therefore, the polynomial series will be adequate provided the problem displace- ments are smooth and fairly moderate in variation. Examples of variations in forms in which polynomials may not be suitable are presented in Fig. 3 -2 3.3.1 Normalisation of Displacement Functions In the above forms, the displacement functions are described by the dimensional coordinates. Non-dimensionalising these coordinates offers a few advantages. If the tube physical dimensions are of extremely different orders of magnitude it is possible to obtain numerical instability when performing the analysis. By implementing nor- malised dimensions the effect of this problem is reduced. To illustrate this problem consider a tube with a large length compared with the radius. In this case, if the co- ordinate z is raised to a high power in the polynomial expansion the result will be an extremely large number. However, the relatively small radial dimension will result in a number which is very small in comparison to the length. Therefore, round-off error can lead to erroroneous results. It also aids significantly in the mathematical manipulations which are required. Thus the coordinates were normalised according to the transformation mapping: 0 ? z ? L 0 ? Z ? 1 ? z 7?? Z ? L (12) 25 Figure 3 -2: Variations Unsuitable for Polynomial Approximations ai ? r ? bi 0 ? R ? 1 ? r 7?? Rti + ai (13) The dimensions ai and bi are indicated in Fig. 3 -3 It follows that the differentials can be written as: ?z = L?Z (14) ?r = ti ?R (15) Incorporating these normalised variables into the displacement functions yields: u(Z,R) = M? m=0 N? n=0 WmnZmRn (16) w(Z,R) = P? p=0 Q? q=0 WpqRpZq (17) 26 Figure 3 -3: Tube Cross-Section with Associated Dimensions As the strain is explicitly dependent on the displacements, it is also necessary to incorporate the normalised coordinates into the strain relationship. This is a simple procedure conducted by an application of the chain rule for differentiation. As as example, the axial strain, ?z is treated below. The procedure can easily be extended to the other components of the strain tensor. ?z = ?u ?z = ?u ?Z ?Z ?z = 1 L ?u ?Z (18) 3.4 Application of the Rayleigh-Ritz Technique In general it is possible to express the potential energy of a system as: ? = ? ? ? (19) Where ? represents the system internal energy, and ? is the external work term. The stability of a system is characterised by the position of minimum potential energy. The Rayleigh-Ritz technique utilises the chosen displacement function and minimises the potential energy with respect to these functions. Thus, the greater the ability of the displacement functions to capture the actual behaviour of the system, the more accurate the solutionxv. The position of minimum potential energy may xvan interesting result is that if the displacement functions are able to capture the system be- haviour exactly, the Rayleigh-Ritz solution will generate the exact result. 27 therefore be determined as: ?? ?W = 0 = ?? ?W ? ?? ?W (20) Which can be rearranged to yield the familiar form: ?? ?W = ?? ?W (21) For the case of thermal loading, the system potential energy can be expressed as: ? = 1 2 ? ? V ? ?TC? dV (22) Where ? is the resultant strain of the system (i.e. system strains less the thermal strains for thermal loading). Eqn. 20 can then be written as: ?? ?W = 0 = 1 2 ? ? V ? ? ?W ( ?TC? ) dV (23) For convenience, the strain tensor can be rewritten as the system strains less the thermal strains (or other relevant strains) in the form: ? = ? ??? ?z ? ?z?T ?? ? ???T ?r ? ?r?T ?rz ? ??? = ?1 ? ?2 (24) Where: ?1 = ? ??? ?z ?? ?r ?rz ? ??? , ?2 = ? ??? ?z ?? ?r 0 ? ????T (25) By performing some simple mathematical manipulation it is possible to reduce Eqn. 23 to the familiar form: KW = L (26) 28 Where: L = ? ? V ? BTC?2 dV (27) K = ? ? V ? BTCB dV (28) A more detailed description of the necessary manipulations is presented in Appendix A. An important factor in Eqn. 28 is that the stiffness matrix does not in any way depend on the applied loading. This is one of the major advantages of the Rayleigh- Ritz technique. It is only necessary to calculate the stiffness matrix for a particular system once; it is then possible to analyse any loading on the system by simply modifying the loading matrix. At this point, some clarification of the results in Eqn. 28 is required. These equations imply that the loading matrix is independent of the system strains. This may be misleading as it would seem to imply that system strains in no way create stress states. It must be remembered that the loading matrix is, in essence, related to the external work terms. This external work can be of many forms, for example, strain energy due to mechanical loading. It must therefore be remembered that the above derivation is based on the assumption that the only form of loading is due to thermal expansion. It therefore makes sense that in the absence of a temperature differential, ?T , there will be no stress state. The inclusion of other forms of external work is dealt with later in this work. In the case of a multi-layered tube, the loading and stiffness matrix can be calculated for each layer. The individual matrices can then be compiled into a single global matrix equation in order to describe the full tube. KT = ? ???? K1 0 ? ? ? 0 0 K2 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? Ki ? ???? , LT = ? ???? L1 L2 . . . Li ? ???? (29) Where the subscript T denotes tube. Up to this point, reference has been made to multiple layers in the radial direction with the inference that only a single polynomial describes the displacement along the length of the tube. As it is expected that the displacements will vary by considerable amount along the length of the tube, it is possible that a single polynomial will not be adequate to accurately capture the behaviour. Thus a discretisation scheme was implemented to allow the length of the tube to be discretised into a number of regions. The process of implementing this scheme is exactly analogous to the 29 method for accounting for multiple radial layers as described above. Again, for each axial region a stiffness and loading matrix is calculated and in turn these matrices are compiled into Eqn. 29. By virtue of the normalised displacements it must be recognised that each region is transformed into a unit square region. Thus the essential form of the stiffness matrix for each region is identical. This has significant advantage in the formulation and implementation of the model as it requires the minimum of computation. 3.5 Implementation of Constraints As formulated thus far, Eqn. 29 will lead to a singular matrix problem. This is of course due to the fact that general polynomial displacement functions are utilised and hence, no boundary conditions have been applied to the problem. As can be seen in Eqn.?s 10 and 11, provision is made to allow non-zero initial displacements, and thus rigid body motion is possible. In order to perform the Rayleigh-Ritz analysis, it is required that the geometric boundary conditions are fulfilled for C0 compatibilityxvi. Importantly it should be noted that although the natural boundary conditions need not be satisfied, it is imperative that they are not inherently violated either. If this is the case, the governing differential equations will not be satisfied. The general nature of Eqn.?s 10 and 11 assures that the natural boundary conditions will not be inherently violated. Suppose that the constraints can be written in the form: cW = 0 (30) Then it is possible to augment the stiffness matrix with the matrix c in the form: [ KT cT c 0 ] [ W ? ] = [ LT 0 ] (31) Applying the constraints in this way is known as the Lagrange Multiplier Technique. The Lagrange Multiplier technique has the advantages that by implementing Eqn. 31 the constraint forces are immediately solved. Thus (in this case) the interface forces between regions can be determined without further calculation. Additionally formulating the constraints in this way is a simple and quick process. xviThe Cm notation is common when discussing variational techniques. Simply stated, C0 com- patibility implies that a function, ?, is continuous within a set space, s, but the derivatives, ???s , need not be. For this problem, the displacement within the system must be continuous, but the slopes between regions need not be. Therefore this implies that only displacement compatibility must be ensured33. 30 An alternative method for implementing the constraints was presented by Baharlou and Leissa34. This technique decomposes the global matrix problem into two smaller matrix problems. This is an advantage over the Lagrange Multiplier method in that the size of the matrix to be inverted is reduced, whereas the Lagrange Multiplier technique increases the size of the matrix which could lead to numerical instability. Of course, the corollary is that now two matrix problems must be solved which could require more computation time. However, it was felt that the increase in stability out-weighed the disadvantage of extra runtime and the method of Baharlou and Leissa34 was implemented in this work. This technique is developed in the following section. 3.5.1 Constraint Matrix Approach The method of Baharlou and Leissa34 is developed below. Extension is made to the original presentation of this method to allow non-zero constraints. Suppose the constraints can be written as: cW = d (32) As the model in this work is fairly complex, it is probable that there will be duplicate constraints applied simply by enforcing the required compatibility. Needless to say, this obviously leads to a singular stiffness matrix. Thus the initial step in dealing with the constraint matrix is to perform Gauss-Jordan elimination to obtain the matrix in its row reduced echelon form. At this point there will a number of zero filled rows evident by the fact the rank of the matrix, c, will be less than the number of rows. These rows must then be removed from the constraint matrix. By inference, Eqn. 32 suggests some interdependency between the elements of W . This is expected as the problem is statically indeterminant. The interdependence of the displacement coefficients can be simplified by considering a group of inde- pendent coefficients and a corresponding group of coefficients which are dependent on the independent coefficients. Now, except in trivial cases, the independent and dependent displacement coefficients are not arranged sequentially within the matrix W . Therefore it is necessary to perform some transformation to separate the dis- placement coefficient groups. Consider writing the displacement coefficients in the form: W = Tr [ ?W ?W ] (33) Written in the above form, ?W are considered to be independent coefficients and ?W are considered to then be dependent coefficients, i.e. coefficients which are dependent 31 on ?W . Tr is a transformation matrix. The exact form of this inter-dependence of the displacement coefficients is illustrated in the following derivation. Substituting into Eqn. 32 yields: cTr [ ?W ?W ] = d (34) Following basic linear matrix properties, the matrix product cTr can be written as: cTr = [ a : b ] (35) Such that: a ?W + b ?W = d ?W = b?1 ( ?a ?W + d ) (36) The above relationships clearly indicate the nature of the relationship between ?W and ?W and the dependence of ?W on ?W . Continuing, the displacement coefficients can be written as: W = Tr [ ?W b?1 ( ?a ?W + d ) ] W = Tr [ I ?b?1a ] ?W + Tr [ 0 b?1d ] W = Tr G ?W + Tr F (37) Following the model derivation in Appendix A, the matrix problem can be redefined as: KC = GTTrTKTTr G LC = GTTrTL?GTTrTKTTrF KC ?W = LC (38) As can now be seen, the matrix KC has smaller dimensions than the original stiffness matrix, KT . Thus the size of the matrix problem has been reduced. However, it is required to solve Eqn. 37 as well as Eqn. 38 to fully specify the problem. 32 It is necessary to elaborate on Eqn. 33 in order to explain how to separate the dependent and independent coefficients. Evidently, for problems of the magnitude solved in this work it is not possible (and certainly not feasible) to perform this oper- ation by inspection. However, looking at Eqn. 36 shows that the product cTr must produce an invertible matrix b. This can be obtained by arranging the columns and rows of the constraint matrix to obtain a square identity matrix at the end of matrix. This same column permutation is repeated on an identity matrix which possesses the same dimensions as the stiffness matrix. The resulting constraint matrix is in the form of cTr whilst the modified identity matrix becomes the transformation matrix, Tr. Whilst the implementation of the constraints has been treated in this section, it is still necessary to obtain the required constraint matrix. The form of this matrix is discussed in the following section. 3.6 Formulation of Constraints As mentioned previously it is only necessary to satisfy the geometric boundary conditions of the problem, these conditions can be summarized for the following cases: - Symmetry - Displacement Continuity - Forced Displacements 3.6.1 Symmetry Conditions at the Tube Midplane In a large number of cases the problem can be assumed to be symmetrical about the mid line of the length of the tube, thus only half the length need be modelled. In this case it is necessary to implement symmetry boundary constraints at the mid line. This can be written in the form: u|Z=0 = 0 ? R (39) M? m=0 N? n=0 Wmn0mRn = 0 (40) Now by definition, 0m = { 1 m = 0 0 m 6= 0 (41) 33 Therefore Eqn. 40 can be written in matrix form as: [ 1 R R2 ? ? ? RN 0 ? ? ? 0 ] ? ??????????????? W00 W01 W02 . . . W0N . . . WM0 . . . WMN ? ??????????????? = 0 (42) An inspection of Eqn. 42 indicates that the constraint is both a function of R and of Wmn. However, since R is not constant but varies through the thickness, the simplest, general solution to Eqn. 42 is given by: ? ?????? W00 W01 W02 . . . W0N ? ?????? = ? ?????? 0 0 0 . . . 0 ? ?????? (43) For purposes of implementation, the above equation can be considered as: ? ?????? 1 0 0 ? ? ? 0 0 1 0 ? ? ? 0 0 0 1 ? ? ? 0 . . . . . . . . . . . . . . . 0 0 0 ? ? ? 1 ? ?????? ? ?????? W00 W01 W02 . . . W0N ? ?????? = ? ?????? 0 0 0 . . . 0 ? ?????? IN ? ?????? W00 W01 W02 . . . W0N ? ?????? = ? ?????? 0 0 0 . . . 0 ? ?????? (44) Where IN is the identity matrix of order N + 1. IN can then be included in the constraint matrix c. The same procedure is followed for the formulation of the additional constraint matrices. The forms of the constraints for compatibility between axial and radial regions are treated in the ensuing sections. 34 3.6.2 Displacement Compatibility between Axial Regions Axial Displacement: u(1)(1, R) = u(2)(0, R) (45)( M? m=0 N? n=0 Wmn1mRn ) (1) ? ( M? m=0 N? n=0 Wmn0mRn ) (2) = 0 (46) Leading to: ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? W00 W01 . . . W0N ? ???? (1) + ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? W10 W11 . . . W1N ? ???? (1) + ... + ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? WM0 WM1 . . . WMN ? ???? (1) ? ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? W00 W01 . . . W0N ? ???? (2) = ? ???? 0 0 . . . 0 ? ???? INW (1)0n + INW (1) 1n + ? ? ? + INW (1) Mn ? INW (2) 0n = 0 [ IN IN ? ? ? IN ] W (1)mn ? [ IN 0 ? ? ? 0 ] W (2)mn = 0 (47) 35 Radial Displacement: w(1)(1, R) = w(2)(0, R) (48)( P? p=0 Q? q=0 Wpq1pRq ) (1) ? ( P? p=0 Q? q=0 Wpq0pRq ) (2) = 0 (49) Leading to: ? ???? 1 1 ? ? ? 1 0 0 ? ? ? 0 0 0 ? ? ? 0 0 0 ? ? ? 0 1 1 ? ? ? 1 0 0 ? ? ? 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ? ? ? 0 0 0 ? ? ? 0 1 1 ? ? ? 1 ? ???? ? ???????????????????????? W00 W01 . . . W0Q W10 W11 . . . W1Q . . . WP0 WP1 . . . WPQ ? ???????????????????????? (1) ? ? ???? 1 0 ? ? ? 0 0 0 ? ? ? 0 0 0 ? ? ? 0 0 0 ? ? ? 0 1 0 ? ? ? 0 0 0 ? ? ? 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ? ? ? 0 0 0 ? ? ? 0 1 0 ? ? ? 0 ? ???? ? ???????????????????????? W00 W01 . . . W0Q W10 W11 . . . W1Q . . . WP0 WP1 . . . WPQ ? ???????????????????????? (2) = ? ???????????????????????? 0 0 . . . 0 0 0 . . . 0 . . . 0 0 . . . 0 ? ???????????????????????? (50) 36 3.6.3 Displacement Compatibility between Radial Regions Axial Displacement: u(1)(Z, 1) = u(2)(Z, 0) (51)( M? m=0 N? n=0 WmnZm1n ) (1) ? ( M? m=0 N? n=0 WmnZm0n ) (2) = 0 (52) Leading to: ? ???? 1 1 ? ? ? 1 0 0 ? ? ? 0 0 0 ? ? ? 0 0 0 ? ? ? 0 1 1 ? ? ? 1 0 0 ? ? ? 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ? ? ? 0 0 0 ? ? ? 0 1 1 ? ? ? 1 ? ???? ? ???????????????????????? W00 W01 . . . W0N W10 W11 . . . W1N . . . WM0 WM1 . . . WMN ? ???????????????????????? (1) ? ? ???? 1 0 ? ? ? 0 0 0 ? ? ? 0 0 0 ? ? ? 0 0 0 ? ? ? 0 1 0 ? ? ? 0 0 0 ? ? ? 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 ? ? ? 0 0 0 ? ? ? 0 1 0 ? ? ? 0 ? ???? ? ???????????????????????? W00 W01 . . . W0N W10 W11 . . . W1N . . . WM0 WM1 . . . WMN ? ???????????????????????? (2) = ? ???????????????????????? 0 0 . . . 0 0 0 . . . 0 . . . 0 0 . . . 0 ? ???????????????????????? (53) 37 Radial Displacement: w(1)(Z, 1) = w(2)(Z, 0) (54)( P? p=0 Q? q=0 WpqZp1q ) (1) ? ( P? p=0 Q? q=0 WpqZp0q ) (2) = 0 (55) Leading to: ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? W00 W01 . . . W0Q ? ???? (1) + ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? W10 W11 . . . W1Q ? ???? (1) + ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? WP0 WP1 . . . WPQ ? ???? (1) ? ? ???? 1 0 ? ? ? 0 0 1 ? ? ? 0 . . . . . . . . . . . . 0 0 ? ? ? 1 ? ???? ? ???? W00 W01 . . . W0Q ? ???? (2) = ? ???? 0 0 . . . 0 ? ???? IQW (1)0n + IQW (1) 1q + ? ? ? + IQW (1) Pq ? IQW (2) 0q = 0 [ IQ IQ ? ? ? IQ ] W (1)pq ? [ IQ 0 ? ? ? 0 ] W (2)pq = 0 (56) Note that in each corner (i.e. Z = R = 1 say) for a particular region the axial displacement is dictated by the neighbouring axial region (Eqn. 46) in addition to being dictated by the adjacent radial region (Eqn. 52). Therefore, as alluded to in the previous section, the constraints will be duplicated. It is interesting to note from the above relationships an additional advantage to non- dimensionalising the displacements is evident (see ?3.3.1). In the form above, the constraint matrix for each of the regions follows the same core form. Thus, for example, a general matrix can be determined for the axial compatibility between regions. All that is required to be done then is to position this matrix such that it references the correct set of displacement coefficients. This reduces the computation significantly as the constraints need not be recalculated. 38 3.7 Formulation of Pressure Loading The effect of pressure can be incorporated into the loading matrix by considering the work done by the pressure load. In this derivation, the pressure load is assumed to be applied along the full length of a single axial region. This assumption is not strictly necessary although it simplifies both the derivation and the implementation of the loading. In this way provision is made to allow for a discrete pressurised region or a tube under pressure along its entire length. The derivation for external pressure over a discrete length is presented below. The idealised loading situation is shown in Fig. 3 -4. The approach followed below can readily be adapted for internal pressure. Figure 3 -4: Tube with Discrete External Pressure According to Fig. 3 -4 the differential radial force due to the pressure can be written as: dF = p ro dz d? (57) The work done by this force can be calculated from first principles as: d?e = ? dF ? d~r (58) Where ~r is a vector in the direction of the resulting displacement. By recognising that the dependence of the force on the circumferential coordinate can be neglected, the following relationship can be attained: d?e = ? p ro dz ? d~r (59) 39 An inspection of Eqn. 59 reveals that the integrand is independent from the inte- gration variable which yields the trivial result: d?e = p ro dz w(Z, 1) d?e = p ro dz P? p=0 Q? q=0 WpqZq1p (60) Thus the total work can be obtained by integrating Eqn. 60 with respect to the axial coordinate. The result is given in Eqn. 61. ?e = ? 1 0 p ro` P? p=0 Q? q=0 WpqZq1p dZ ?e = p ro` P? p=0 Q? q=0 Wpq 1 q + 1 (61) The above result can then be incorporated into the external work term, ?, in the mechanical analysis derivation. 3.8 Determination of Material Properties The methods utilised to estimate the material properties of composite laminates are briefly discussed in this section. It should be pointed out that all the meth- ods described here are largely approximate and test data should be used wherever possible. The variation in the makeup of composite materials contributes largely to the vari- ability in material properties. In metals the elastic properties are generally inherent to the metal in question. It is thus possible to define standards for these materials. Composite materials differ in that the process by which the part is manufactured and cured, can change the properties of the final material drastically. These processes are not repeatable to a very high degree and possess a large number of variables, examples are the fibre volume fractions, void content, degree of cure, etc. Thus it is not possiblexvii to obtain standards or a database of properties for composite materials. To overcome this there are methods whereby the properties for the laminate can be approximated from readily available data. As the constituents of the composite xviiAt least at the time of writing. 40 are fairly regular, most data is readily attainable for the fibre and the matrix. Mi- cromechanics is a common approach which utilises the matrix and fibre properties to predict the final laminate properties. This is accomplished by rule-of-mixtures approaches and other simple techniques4. In traditional composite material analyses, micromechanics provides sufficient in- formation to obtain all required information, essentially all the in-plane properties. However in this work, the primary concern is with through-thickness effects which implies that some knowledge as to the out-of-plane properties is required. As these effects are neglected in traditional analyses, there exists some difficulty in finding data and methods to determine these properties. Usually the through-thickness stiffness of a ply, E3, is assumed to be equal to the transverse stiffness of a uni-directional ply with the same fibre volume fraction, (E2)uni. This assumption relies on the fact that the structure in the through- thickness direction is very similar to the transverse direction in a uni-directional laminate. Methods for determining through-thickness Poisson?s ratios are presented by Her- akovich4; as are methods of determining the through-thickness thermal expansion coefficients. To determine the Poisson?s ratio, ?23, an alternate approach, based on micromechanics, is presented in Appendix B. Other data is calculated as mentioned in Stone et al.7. The shrinkages and CTE?s can be calculated by the approach of Schapery35. Typically, the properties are determined for a ply in terms of a coordinate system orientated as shown in Fig. 3 -5. Thus the axes are aligned in the fibre direction. However, plies are often incorporated into the laminate at other angles. Therefore it is convenient to perform the analysis in a coordinate system which is aligned with the laminate (which is typically the same coordinate system as the loading and boundary conditions). The ply constitutive properties can be described by the compliance matrix, S. As the constitutive relationships are tensor based, the properties for the ply (in this case, S) can be transformed into an alternative coordinate system by means of a rotation matrix, in this case denoted as T . Thus the properties for off-axis plies can be determined by rotating the compliance matrix through the desired angle by the equation: ?S = RT?1R?1ST (62) 41 The exact definition of the rotation matrix can be determined from tensor algebra as: T = ? ??????? c2 s2 0 0 0 2cs s2 c2 0 0 0 ?2cs 0 0 1 0 0 0 0 0 0 c ?s 0 0 0 0 s c 0 ?cs cs 0 0 0 c2 ? s2 ? ??????? (63) Where c and s are defined as cos(?) and sin(?) respectively, and ? is the angle of rotation about the through-thickness axis, axis 3 in Fig. 3 -5. R is the Reuter?s matrix defined as: R = ? ??????? 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 ? ??????? (64) The Reuter?s matrix is necessary for the conversion from tensor strain to engineer- ing strains. The transformed stiffness matrix can then be found by inverting the compliance matrix thus, ?C = ?S?1 (65) Figure 3 -5: Laminate Coordinate System 42 3.9 Notes on Parametric Analyses The point was raised in the introduction that the Rayleigh-Ritz analysis had the advantage over the FEM in that the problem could be formulated in such a way to assist parametric analyses. To validate this point it is only necessary to glance at the mechanical analysis derivation. It can be noted that the entire analysis derivation can be conducted by supplying the tube dimensions: L, ri, t, and by supplying the necessary tube make- up and material properties. In addition to these parameters, it is still necessary to supply the solution parameters: N,M,P,Q,Nl,Nm. Where Nl and Nm are the number of discretised regions along the length and throughout the thickness of the tube respectively. Of course, the input effort required for the material data is the same as for FEM. However, the effort required for the tube make-up and dimensions is far less, and is also far simpler. The discretisation scheme makes allowance for axial regions of differing length. Thus it is possible to decrease the size of the region in areas where there is large displace- ment variation expected. The process for determining the length of the axial regions can be done manually or alternatively, it can be automated to some extent. One of the simplest means for incorporating this is to consider the axial region length as a geometric series, i.e.: Lk+1 = 1 ? Lk (66) In this way, the length of the regions can be allowed to decrease towards the free end of the tube say (where it is expected that the most variation will occur). This approach is very simple as the relationships governing geometric series are well known and can be found in any work on basic mathematics 36. 43 4 Development of Thermal Analysis Technique In line with the development of the mechanical analysis, the heat transfer problem has been treated as axisymmetric. The primary direction of heat transfer is assumed to be in the radial direction, thus heat transfer across the end faces is ignored. The problem therefore becomes one dimensional. The development and basis for the analysis technique utilised in this work is described in the following sections. 4.1 Governing Heat Transfer Processes The law of conservation of energy can be expresses mathematically as presented in Eqn. 67. dU dt = Q+ ?W (67) In reality the problem considered in this work is a thermoelastic system. However as the only work performed by the system will be the strain energy associated with the loading (which is small), it is assumed that the system can be decoupled. Hence, the external work rate of the system, ?W , is zero. The change in internal energy is then characterised by the familiar equation: dU dt = ?V c dT dt (68) Where ? is the material density, V is the volume and c is the material specific heat. It should be noted that Eqn. 68 is written in the general case with c replaced with either cp or cv. This is due to the fact that for fluid media there is a large difference in the specific heat for processes under constant pressure (i.e cp) and processes under constant volume (i.e cv). However this work is concerned with the analysis of solids. To aid analysis, it has been assumed that the solids in question can be treated as incompressible. In this case it is appropriate to apply the assumption: cp ? cv = c (69) The term Q refers to the heat transfer across the system boundary. In general this will be due to one or a combination of the processes: - Conduction - Convection (forced and/or free) 44 - Radiation - Internal heat generation For the purposes of this work, all forms of internal heat generation are neglected. Possibly in some processes (particularly post-cure) where chemical reactions are expected, this assumption will not hold. However in situations where it is significant, a great deal of work will centre around quantifying this process. Radiation is governed by the Stefan-Boltzman Law expressed in the form of radiant heat exchange between two bodies in Eqn. 70 (It should be noted that the tem- peratures specified in Eqn. 70 are in fact absolute temperatures in Kelvin). As can be seen from the relationship in Eqn. 70, the dependence of the heat trans- fer on temperature is to the fourth power, thus making heat transfer to radiation inherently nonlinear. It must also be considered that in order to perform an accu- rate radiation analysis requires good knowledge of the environment of the system and the surroundings. As the heat transfer processes considered in this work are occurring across a fairly low absolute temperature, it was felt that the amount of radiant heat transfer would not be significant and therefore efforts can concentrated on conduction and convection. Qrad = kstA ( T 41 ? T 42 ) (70) 4.1.1 Conduction Heat Transfer Throughout the thickness of the tube, heat is transferred by means of conduction. It is assumed that there is a perfect bond between material layers, hence no contact resistance (although not included here, contact resistance can easily be incorporated into the formulation). Conduction heat transfer is governed by Fourier?s Law (See Eqn. 71). qcond = ?k?T (71) Incorporating Eqn. 68 and Eqn. 71 yields the thermal diffusivity equation presented in Eqn. 72. ?2T = ?ck ?T ?t (72) Following the definition for the Laplacian operator, the left hand side of Eqn. 72 can be expressed in cylindrical coordinates as: 45 ?2T = 1 r ? ?r ( r ?T ?r ) + 1 r2 ?2T ??2 + ?2T ?z2 (73) By virtue of the assumptions made ( ??? = ? ?z = 0), Eqn. 73 can be reduced simply to: 1 r ? ?r ( r ?T ?r ) = ?c k ?T ?t (74) 4.1.2 Convection Heat Transfer While the conduction heat transfer developed in ?4.1.1 is well described, convection heat transfer is more complex. For the basis of this work it is sufficient to assume that the convection is limited to free convection and can be adequately described by Newton?s law of cooling (Eqn. 75). It is of course possible to describe the process in more detail, but that is beyond the scope of this work. Qconv = hcA(T ? T?) (75) Where hc is the convection coefficient, A is the exposed area and T refers to the boundary temperature. In the form of Eqn. 75, hc is assumed to be constant. In reality it can be heavily temperature dependent. The determination of this constant is an involved process on its own made more difficult in the fact that there is an extremely large variance in order of magnitude for the values. 4.2 Resistor-Capacitor Formulation It is common to make an analogy between thermal and electrical systems. By relating the two systems, the problem can often be more easily visualised. The resistor-capacitor formulation37 is very convenient as it allows a problem to be de- composed into smaller building blocks. Thus the implementation of the formulation can be very flexible in terms of boundary conditions and/or loading conditions. The formulation also lends itself to numerical simulation. Mills37 presents an algorithm for applying this formulation to the transient heat transfer problem. The algorithm is based upon the relationships presented in the previous section, a detailed derivation of which can also be found in Mills37. The algorithm is essentially a finite difference, time marching solution. In this way, the nodal temperatures at any time step i + 1 are calculated from the temperatures at time i. In this form, the solution is explicit (i.e. the values obtained in one time step are entirely dependent on the prior time step). It is also possible to formulate 46 the problem implicitly. Each formulation has advantages. The implicit solution allows the temperatures to be determined at any point in time without the need to march to the solution via a number of time increments. As can be inferred from this statement, the implicit solution is inherently stable. The explicit solution is, however, more memory efficient than the implicit solu- tionxviii. While the implicit solution requires either the solution of a matrix problem (possibly by iterative means), the explicit solution makes use of very simple rela- tionships to calculate the temperatures with no need for iteration. From the point of view of including effects such as material nonlinearity xix the ex- plicit version is preferred as this becomes essentially a trivial exercise to implement; whilst the implicit solution requires to be reformulated. The explicit solution was chosen for the reason that the problems investigated in this work are not particularly complex (and hence will not require undue computation time). From a point of view of model flexibility, the explicit solution was preferred as it allows different effects and boundary conditions to be accounted for with relative ease. The discretisation scheme which was implemented for the thermal model can be seen in Fig. 4 -1. The idealised elements are shown in Fig. 4 -2. In particular, Fig. 4 -2 indicates the means by which the volumes are idealised to be connected with a thermal resistor. Figure 4 -1: Discretization Scheme for Thermal Analysis xviiiThe explicit solution requires solving simple equations at every point for a number of time steps. The implicit solution rather solves a matrix problem. As the time length of the solution increases, the computational cost of the explicit solution will overtake the implicit solution (whose computational cost is constant). xixFor the thermal analysis this implies material properties which vary with temperature. 47 Figure 4 -2: Idealised Element System for Thermal Analysis The relationships utilised in this analysis are presented in the ensuing part of this section. The analysis has been conducted on the assumption that a fixed temper- ature is applied at the inner wall of the tube, and the outer wall has been allowed to dissipate heat by means of free convection to an ambient temperature. It is im- portant to note that throughout the formulation of the analysis presented here it is not necessary to assume that the material properties are temperature independent, although this has been assumed to preserve the linearity of the analysis. The basic logic of this analysis is that the temperature at time, i + 1, is equal to the temperature at the current time, i, with the addition of a differential, ?T . This statement is represented mathematically in Eqn. 76. T i+1m = T im + ?T (76) In accordance with the governing relationships of heat transfer presented earlier, the temperature differential can occur due to a combination of change in internal energy and heat transfer by means of conduction or convection. Following Eqn. 68, the internal energy change can be written as: dU dt ? ?U ?t = ?cV?T ?t (77) Thus the change in internal energy at point m is: ?Um = ?mcm?Vm ( T i+1m ? T im ) (78) Which can be written as: ?Um = Cm ( T i+1m ? T im ) (79) 48 Where Cm is the capacitance of the element, in essence this is the ability of the element to store internal thermal energy, given by: Cm = ?mcm?Vm (80) ?Vm is the volume of the element which may be calculated in cylindrical coordinates as: ?Vm = rm?r???z (81) The change in temperature due to conduction at the element boundary can be determined as: qcond = ?k?T Qcond = ?kA?T Where qcond is the heat flux associated with conduction and Qcond is the associated heat transfer due to conduction, given by: Qcond = qcondA. Continuing: Qcond = ?k(rm + ?r/2)???z ( ?T ?r ) Qcond = ?k ( T im ? T in ) (rm + ?r/2))???z ?r (82) Note that in the above relationships, the term: (rm + ?r/2) is the radial distance to the boundary between the elements (see Fig. 4 -2 for clarification). The above relationship can be written in the form: Qcond = ( T in ? T im Rmn ) (83) R refers to the thermal resistance which the element is subject to along its bound- aries. In effect, this describes the heat transfer to other elements by means of conduction. The mathematical implementation is presented in Eqn. 84. The con- vention for determining R+m and R?m is shown in Fig. 4 -2. R+m = ?r (rm + ?r/2)???zk , R ? m = ?r (rm ??r/2)???zk (84) 49 The above relationship makes sense as the heat transfer is proportional to the tem- perature difference between point n and m and is inversely proportional to the thermal resistance between the points. The general expression for the temperature at any point of interest in the tube can, thus, be written as: T i+1m = T im ( 1? ?t Cm ? n 1 Rmn ) + ?t Cm ? n T in Rmn (85) Continuing with the derivation, as the problem is axisymmetric all dependence on the coordinate ? can be ignored (i.e. consider a unit change for ??). As the temper- ature has been assumed to be invariant along the length of the tube, a unit change for the axial coordinate is also employed. Expanding Eqn. 85 yields: T i+1m = T im [ 1? ?t ?mcmrm?r ( (rm + ?r/2)k + (rm ??r/2)k ?r )] ? ? ? + ?t ?mcmrm?r [ T im?1(rm ??r/2)k + T im+1(rm + ?r/2)k ?r ] T i+1m = T im [ 1? k?t ?mcm ( 2 ?r2 )] + k?t ?mcm [ T im?1(rm ??r/2)k + T im+1(rm + ?r/2)k rm?r2 ] (86) The above relationships are applicable to a point within any given layer of material in the interior of the tube wall. For a point on the outer surface of the tube the relationships become: Cm = 1 2 ?mcmrm?r???z (87) Rconv = 1 hcrm???z (88) T i+1m = T im [ 1? 2k?t(rm ??r/2) ?mcmrm?r2 ? 2?t ?r hc ?mcm ] + 2k?tT im?1(rm ??r/2) ?mcmrm?r2 + 2?tT? ?r hc ?mcm(89) At the boundary between two different materials the temperature can be determined by considering a point which consists equally of both materials. 50 Thus, Cmn = 1 2 (?mcm?rm???z + ?ncn?rn???z) (90) T i+1mn = T imn [ 1? ?tkm(rm ??rm/2) Cmn?rm ? ?tkn(rm + ?rn/2) ?rnCmn ] (91) + ?tkmT imn?1(rm ??rm/2) ?rmCmn + ?tknT imn+1(rm + ?rn/2) ?rnCmn Note, allowance has been made in the above relationships for differing radial dis- cretisation lengths, ?r. This allows for the number of points within any material to be selected by the user without regard for element length. For convenience the following symbols are defined: ? = k ?c , h = hc ?c (92) An area which hasn?t been discussed thus far is the stability of the model. As the model is an explicit and marching solution without iteration, the solution is not inherently stable. Instead the stability condition dictates that the coefficient of T im in Eqn. 85 never be negative. Thus: 1? ?t Cm ? n 1 Rmn ? 0 (93) As there are principally two governing mechanisms for the heat transfer (conduction and convection), there will be two forms of the stability criteria equation. In order to ensure a stable solution, the most severe of these constraints must be met. The relationships governing these conditions are presented in Eqn.?s 95 and 96. Note, as this criteria will be driven by the material layer with the smallest region size, there is no need to consider a material boundary element. Following Eqn. 93, it can be said that the largest allowable time step for a stable solution is given by: ?tmax ? Cm ( ? n 1 Rmn ) ?1 (94) 51 Substituting the relationships for the thermal resistance for conductivity (see Eqn. 84) yields the result: ?tmax ? Cm ( (rm + ?r/2)k ?r + (rm ??r/2)k ?r ) ?1 ?tmax ? 1 ? rm?r ( (rm + ?r/2) ?r + (rm ??r/2) ?r ) ?1 ?tmax ? 1 ? rm?r ( ?r 2rm ) Which yields the final result of: ?tmax ? ?r2 2? (95) A similar process can be followed for the convection boundary with the result being: ?tmax ? ro?r2 2 [(ro ??r/2)? + hro?r)] (96) The interdependence between the solution stability and the region size can be in- ferred from the above results. It is also noted that while the solution accuracy improves with decreasing region size, the stability conditions then require that the maximum time step also decreases. There is, therefore, a trade-off between accuracy and run time. An a priori view of this model yields the conclusion that within any homogeneous layer of material, the temperature distribution will be smooth. Thus it is possible to approximate the temperature distribution by means of a least squares polynomial with a high degree of accuracy. These polynomial distributions provide the basis for incorporating the temperature distributions into the mechanical model. 52 5 Validation of Analysis Techniques In light of the previous research in this area, it was prudent to validate the anal- ysis presented in this work against existing results. A number of test cases were selected and are discussed in this section. In some cases the results were validated against established results/techniques and/or against the results obtained from a FEM analysis. In the case of the FEM, a commercial software package was utilised to generate the results. A brief overview of the FEM models is presented in ?5.1. 5.1 Description of FEM Models For sake of completeness a brief outline of the process of constructing and the con- tents of the FEM models is presented in this section. For all FEM models the MSC Software Patran/Nastran software suite was utilised (at the time of writing the latest version is v2005b). MSC.Patran is the pre- and post-processor software wherein the models are constructed. MSC.Nastran is the FEM code and proces- sor software. Nastran is a well known FEM code in common use worldwide. The contents of the FEM models are described below. Where applicable, the relevant Nastran input cards are supplied. Further information can be obtained from the Nastran documentation38. In the same manner as the Rayleigh-Ritz analysis, the FEM models were constructed on the basis of axisymmetry. For reasons unknown to the author, Nastran does not possess a general quad element formulated for axisymmetryxx. Thus the models were constructed from 6-Noded TRIA elements (CTRIAX6). The presence of the midside node indicates that these elements are of higher order. As such, they should produce results of sufficient accuracyxxi. The material properties were defined as 3-D orthotropic by means of the MAT3 card. Note that although the model is 2-D, the assumptions of axisymmetry imply that the model is still in effect, a 3-D model. Displacement conditions were applied with SPC cards and pressure loading was specified by PLOADX1 cards. For models involving shrinkage, the shrinkages were utilised as the material CTE?s and a unit temperature change was applied to the model. Temperature distributions were applied at the element nodes by means of a field within Patran. This is an automated feature whereby the distribution is entered into Patran as a mathematical expression which Patran then evaluates at the nodes during the preprocessingxxii. xxTo confuse the matter, an axisymmetric QUAD element is available for hyper-elastic materials. xxiThey are also far superior to the constant strain tria elements. xxiiIt must be remembered that the thermal loads within Patran are additive. i.e. if a line of nodes along a material boundary have thermal conditions specified by fields from either side, the final applied nodal temperatures will be twice the desired result. 53 The thermal models utilised the card CONV to define the free convection boundary condition. A fixed temperature boundary was specified at the inner wall. 5.2 Thermal Model Validation - Heat Transfer of an Or- thotropic Tube The thermal problem analysed by Kardomateas 8 is exactly described by the thermal analysis presented here. Therefore in order to validate this analysis it was necessary only to utilise the same conditions and compare the results. The problem describes an orthotropic tube with a fixed, known temperature at the inner wall of 100?C and free convection at the outer wall. The ratio of convection coefficient to conduction coefficient is 0.15m?1. The tube has a inner radius of 20mm and a thickness of 16mm. The thermal conductivity is given as 1.12 ? 10?6. A schematic of the loading is presented in Fig. 5 -1. For the simulation, a refined discretisation was employed of 1000 nodes through the thickness. In accordance with the stability requirements, the time step was calculated at 1.14? 10?4s. In order to compare results, the time was non-dimensionalised according to the relationship: ?t = t? (ro ? ri)2 (97) The results are presented in Fig. 5 -2. As can be seen, excellent correlation was obtained. 54 Figure 5 -1: Sketch of Loading Conditions on Orthotropic Tube 0 0.2 0.4 0.6 0.8 1 20 30 40 50 60 70 80 90 100 Non?dimensional Radius Temperature, [? C ] Kardomateas Numerical Analysis t = 1.0 t = 0.5 t = 0.25 t = 10.0 Figure 5 -2: Temperature Distribution Through the Thickness of an Orthotropic Tube 55 5.3 Mechanical Model Validation - Orthotropic Tube with Temperature Profile The second half of the work of Kardomateas8 calculates the stress states based of the temperature results presented in ?5.2. By increasing the length of the tube, the mechanical analysis presented in this work will begin to conform to the assumption of constant strain at points away from the end. In doing so, it will be possible to compare the results obtained from both analyses (Rayleigh-Ritz and Kardomateas). The ideal position to compare the analyses is on the midplane of the tube, as the shear stress must be zero at this location. The temperature distributions as utilised in ?5.2 were also applied to a FEM model of this tube. The material properties of the tube are presented in Table 5 -1. The solution pos- sessed the following parameters: N 5 M 5 P 8 Q 3 Nl 3 Nm 1 Nm and Nl refer to the number of discretized regions in the radial and axial direc- tions respectively. The other parameters are the indices utilised in the displacement function expansions: u(Z,R) = M? m=0 N? n=0 WmnZmRn w(Z,R) = P? p=0 Q? q=0 WpqRpZq There was some difficulty in reproducing the stress and displacement results of Kar- domateas8. There appeared to be agreement for the steady-state stresses (?t = 10). There was no agreement, however, for the cases with a varying through thickness temperature profile. However, excellent agreement between the FEM analysis and the Rayleigh-Ritz solution was obtained at all positions in time. The results obtained from the Rayleigh-Ritz and FEM simulations are presented in Fig.?s 5 -3 - 5 -5. The transient through-thickness radial displacements obtained from the Rayleigh-Ritz analysis are presented in Fig. 5 -6. 56 20 22 24 26 28 30 32 34 36?80 ?60 ?40 ?20 0 20 40 Radial Position, [mm] Stress Tensor, [MPa ] ? zz ? zz (FEM) ??? ??? (FEM) ? rr ? rr (FEM) Figure 5 -3: Orthotropic Tube Stress State (?t = 0.25) 20 22 24 26 28 30 32 34 36?80 ?60 ?40 ?20 0 20 40 Radial Position, [mm] Stress Tensor, [MPa ] ? zz ? zz (FEM) ??? ??? (FEM) ? rr ? rr (FEM) Figure 5 -4: Orthotropic Tube Stress State (?t = 0.5) 57 20 22 24 26 28 30 32 34 36?60 ?50 ?40 ?30 ?20 ?10 0 10 20 30 40 Radial Position, [mm] Stress Tensor, [MPa ] ? zz ? zz (FEM) ??? ??? (FEM) ? rr ? rr (FEM) Figure 5 -5: Orthotropic Tube Stress State (?t = 10) 20 22 24 26 28 30 32 34 36 ?0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 Radial Position, [mm] Radial Displacement, [mm ] t = 10.0 t = 1.0 t = 0.5 t = 0.25 Figure 5 -6: Transient Through-Thickness Radial Displacements of Orthotropic Tube 58 Table 5 -1: Material Properties Property Magnitude Ez (GPa) 13.7 E? (GPa) 55.9 Er (GPa) 13.7 Grz (GPa) 4.9 ??z 0.277 ?zr 0.4 ?r? 0.068 ?z (?C?1) 40? 10?6 ?? (?C?1) 10? 10?6 ?r (?C?1) 40? 10?6 5.3.1 Remarks The heat transfer analysis developed in this work was validated in ?5.2. The tem- perature results presented in Fig. 5 -2 indicate an excellent agreement with the analytical analysis with less than 1% error. Additionally, the analysis was found to be computationally efficient (Approx. 10 min on a standard desktop PC). As the thermal analysis developed in this work is relatively simple, it is highly possible that more efficient and accurate (or both) techniques exist to predict this type of heat transfer. However, as the problems under consideration are simple, this analysis was found to be adequate to describe them. The added advantage of the simplicity of the analysis is the ease with which modifications can be made. By applying the temperature distributions from Kardomateas 8 (for which there was excellent agreement) to the mechanical analysis it was expected that in similar fashion, the results would agree. This was not the case. It was found that the Rayleigh-Ritz results agreed well with the results from FEM at all points in time which were analysed. However, it appears that the only point in time where these results appear to agree with the published results is at the point of steady state temperature. Noticeably, the behaviour of the radial stress calculated by Kardomateas8 at ?t = 0.25 shows an inflection point. This behaviour cannot be seen in the Rayleigh-Ritz re- sults. Additionally the radial displacement as originally presented appears to de- crease with time at the inner wall. As the tube is subject to an elevated temperature which increases throughout the thickness with time, this behaviour is not intuitive. The behaviour of the radial displacement obtained with the Rayleigh-Ritz analysis, presented in Fig. 5 -6, does follow what would be expected; thus the radial displace- ment begins at its minimum at ?t = 0.25 and increases steadily to its maximum value at ?t = 10. Considering this information in conjunction with the correlation between the results from this work and FEM in this case, as well as the other validation cases presented later in this section, the conclusion is that the mechanical results 59 published by Kardomateas8 are erroneous. Further, as the results of Kardomateas 8 and the present work agree at steady state (with a constant temperature distribu- tion) the error can be attributed to the manner in which the temperature profile has been incorporated within the mechanical model of Kardomateas8. 5.4 Rayleigh-Ritz Analysis with the Constant Axial Strain Assumption Following the problem formulation, as the length of the tube increases, the current mechanical analysis must conform to the assumption of constant axial strain for areas away from the ends. Thus the analysis should be exactly equivalent to the one utilised in Stone et al.7. To assess this, a tube with the following characteristics was modelled. L = 2000mm, ri = 250mm, t = 250mm The tube material properties are listed in Table 5 -2. The tube is subject to shrink- ages associated with a post-cure operation, the material free shrinkages are also listed in Table 5 -2. Table 5 -2: Material Properties Property Magnitude Ez (GPa) 10.80 E? (GPa) 10.80 Er (GPa) 6.03 Grz (GPa) 4.10 ?z? 0.315 ?zr 0.34 ??r 0.34 sz ?1500? 10?6 s? ?1500? 10?6 sr ?5000? 10?6 The input data for the solution was: N 5 M 5 P 8 Q 4 Nl 3 Nm 1 ? 2 60 The parameter ? refers to the scaling parameter in the geometric series distribution in Eqn. 66. In order to calculate the results for the analysis neglecting shear effects (i.e. the constant axial strain assumption), the mechanical analysis presented in this work was reformulated. By removing the dependence on the parameters M and Q in the displacement variations (see Eqn.?s 16 and 17) as well as removing the shear strain from the formulation (i.e. simply neglecting ?rz in Eqn. 2), the shear deformation is prevented and therefore the tube behaves as if it has infinite length. For comparison, the results were compared across the thickness at the tube midplane as the shear stress and strain at this position must be zero. The resulting comparison is presented in Fig. 5 -7. As is expected, excellent agreement was attained. 250 300 350 400 450 500 ?15 ?10 ?5 0 5 10 15 20 Radial Position, [mm] Stress Tensor, [MPa ] ? zz (G rz =?) ? zz (Present) ??? (Grz=?) ??? (Present) ? rr (G rz =?) ? rr (Present) ? rz (Present) Figure 5 -7: Comparison between Constant Strain and Finite Length Models at the Tube Midplane 61 5.5 Results for a Finite Length, Orthotropic Tube with Re- gard to the Equilibrium Conditions The simulation conducted utilising the analysis presented in this work in ?5.4, in- corporated a finite length tube and, although, the results were presented on the midplane, the simulation results are valid for the entire tube. Results obtained from this simulation are presented in this section to investigate the ability of this work to accurately satisfy the required equilibrium conditions. To this end, the displacement along the length at a position halfway through the thickness is presented in Fig. 5 -8. The corresponding stress states at this position are presented in Fig. 5 -9. The stress state through the thickness at the free end is presented in Fig. 5 -10. As can be seen in Fig. 5 -10, the axial stress decays to zero. The shear stress is also equal to zero at this point. Thus, the equilibrium conditions for the tube are satisfied. Fig. 5 -9 shows a substantial variation in the stress along the length. In particular it can be seen that for the majority of the tube the shear stress is zero. It then peaks in a small region near the free end before returning to zero. This behaviour is intuitive and follows what would be expected. The variation in the magnitude of the stress along the length can clearly be seen; in some cases the stress away from the ends being greater than that predicted at the midplane. As this behaviour is impossible to predict with the ?infinite length? assumptions, it is a clear indicator of the merit of this work. 62 0 200 400 600 800 1000 1200 1400 1600 1800 2000 ?3.5 ?3 ?2.5 ?2 ?1.5 ?1 ?0.5 0 Axial Position, [mm] Displacement, [mm ] Axial Displacement Radial Displacement Figure 5 -8: Displacement Variation Along Tube Length (R = 0.5) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 ?4 ?3 ?2 ?1 0 1 2 3 Axial Position, [mm] Stress Tensor, [MPa ] ? zz ??? ? rr ? rz Figure 5 -9: Axial Variation of Stress (R = 0.5) 63 250 300 350 400 450 500 ?15 ?10 ?5 0 5 10 15 Radial Position, [mm] Stress Tensor, [MPa ] ? zz ??? ? rr ? rz Figure 5 -10: Radial Variation of Stress at Free End (z = L) 5.5.1 Remarks The results presented in ?5.4 and ?5.5 demonstrate the ability of the present model to predict the behaviour of a finite, orthotropic tube. The success of this simulation can be determined by comparing the results to the equilibrium conditions. The results presented indicate that all required equilibrium conditions are satisfied with high accuracy. The axial and shear stress at the free end is zero, the radial and shear stress at the inner and outer walls are zero as well. In addition to satisfying the equilibrium conditions, the simulations were conducted for relatively few DOF?s. Thus the final solution was both accurate and efficient. 64 5.6 Thick-Walled Isotropic Pressure Vessel Calculations In order to validate the implementation of the pressure loading, a simulation was completed for an isotropic thick-walled cylinder under internal pressure. In this way, the results from the simulation can be compared to the Lame? equations (which are, for convenience, cited in Roark6) which are reproduced in Eqn.?s 98 - 100. The correlation is presented in Fig. 5 -11. ?zz = 0 (98) ??? = p r2i (r2o + r2) r2(r2o ? r2i ) (99) ?rr = ?p r2i (r2o ? r2) r2(r2o ? r2i ) (100) The tube under consideration possessed the following characteristics: E = 209GPa, ? = 0.3, ri = t = 120mm, p = 10MPa. As can be seen, the tube is fairly thick (t/ri = 1). As the results of Roark6 are based on a long tube, symmetry boundary conditions were enforced at both ends of the tube. This limits the axial displacement, essentially forcing the shear strain to be zero, thereby modelling a long tube. The solution parameters are: N 8 M 8 P 8 Q 8 Nl 1 Nm 1 120 140 160 180 200 220 240 ?10 ?5 0 5 10 15 20 25 30 35 Radial Position, [mm] Stress Tensor, [MPa ] ? zz (Rayleigh?Ritz) ? zz (Roark) ??? (Rayleigh?Ritz) ??? (Roark) ? rr (Rayleigh?Ritz) ? rr (Roark) ? rz (Rayleigh?Ritz) Figure 5 -11: Isotropic Cylinder With Internal Pressure 65 5.6.1 Remarks As seen from Fig. 5 -11, the results obtained from the current analysis coincide with those obtained from the Lame? equations. The accuracy provided by the current analysis must also be evaluated in conjunction with the fact that the pressure vessel analysed in this section possesses a high thickness ratio, i.e. t/r = 1, thus the thick-walled effects will be highly prevalent. This indicates that the current analysis technique can be utilised to adequately predict the effects found in thick-walled pressure vessels. 5.7 Free End Stresses Attained in an Unrestrained Multi- Material Tube under Internal Pressure As discussed earlier, in the presence of a free edge in composites, numerical analyses tend to predict a stress singularity. While it is not the purpose of this work to investigate this phenomenonxxiii it is interesting to investigate the results produced by this analysis for these effects. As these phenomena have been under rigorous analysis, a number of test cases exist. The work of Ye and Sheng39 compares the results obtained from their model (which is specifically formulated to predict edge-effects) to results obtained by Bhaskar et al.40 from a 3-D FEM analysis. Ye and Sheng39 consider a multi-layer tube with free ends under internal pressure. Therefore a stress singularity will be predicted at the material interface at the free end. This analysis becomes a good test for the present method as the results obtained can be compared to both Ye and Sheng39 and the FEM results. To this end, the analysis of Ye and Sheng39 was reproduced utilising this work. These results are conducted for a two layer tube (with a 0?/ 90? configuration) subject to an internal pressure with no end restraint applied. The material properties utilised for the simulation are specified in Table 5 -3. Table 5 -3: 0? Material Layer Properties Property Magnitude E11 (GPa) 137.9 E22 (GPa) 14.5 E33 (GPa) 14.5 G13 (GPa) 5.86 ?12 0.21 ?23 0.21 ?31 0.21 xxiiiThe analysis is involved and deserving of its own work alone. 66 The solution parameters are: N 4 M 4 P 10 Q 3 Nl 2 Nm 12 ? 1.2 In the analysis of Ye and Sheng39 the results are presented as normalised parameters. As the stresses are normalised according to the internal pressure (i.e. ??rr = ?rr/p) this analysis utilised a unit internal pressure, thus the calculated stresses can be compared directly without any need for manipulation. As stated by Ye and Sheng39, the tube internal radius to thickness ratio was taken as 10 and the length is equal to the inner radius. Therefore, for simplicity, a unit thickness (i.e. 1mm) was selected. As all the analyses involved (Rayleigh-Ritz, FEM and Ye and Sheng39) are fully linear, this parameter normalisation can be performed without incurring any error. The results are presented in Fig.?s 5 -12 - 5 -16. As the area of interest is the material interface, results for the stress along the length are presented at this location in Fig. 5 -12. As can be seen for the majority of the length of the tube the stress state is constant. However, in an area near the free end a large variation in stress can be seen, which is indicative of the free-edge effect. At this scale, the correlation between the 3 sets of results is good. For a more detailed investigation of the results near the free end, details of Fig. 5 -12 are presented in Fig.?s 5 -13 and 5 -14. Fig. 5 -13 and Fig. 5 -14 show that the results obtained from the present model correlate well with the results from Ye and Sheng 39 and those obtained from Bhaskar et al.40. The corresponding displacements are presented in Fig. 5 -15 and Fig. 5 -16. It is interesting to note that although a stress singularity is evident at the material interface at the free end, no discernible disruption can be found in the displacement field. 67 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ?1 ?0.5 0 0.5 1 1.5 2 Normalised Axial Position Stress, [MPa ] ? zz (Rayleigh?Ritz) ??? (Rayleigh?Ritz) ??? (Ye et al) ??? (Bhaskar) ? rr (Rayleigh?Ritz) ? rr (Ye et al) ? rr (Bhaskar) ? rz (Rayleigh?Ritz) Detail 2 Detail 1 Figure 5 -12: Axial Stress Variation in Multi-Layer Tube under Internal Pressure 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 1.8 1.82 1.84 1.86 1.88 1.9 1.92 Normalised Axial Position Stress, [MPa ] ??? (Bhaskar) ??? (Ye et al) ??? (Rayleigh?Ritz) Figure 5 -13: Axial Stress Variation in Multi-Layer Tube under Internal Pressure (Detail 1) 68 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1 ?0.9 ?0.85 ?0.8 ?0.75 ?0.7 ?0.65 ?0.6 ?0.55 Normalised Axial Position Stress, [MPa ] ? rr (Ye et al) ? rr (Bhaskar) ? rr (Rayleigh?Ritz) Figure 5 -14: Axial Stress Variation in Multi-Layer Tube under Internal Pressure (Detail 2) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ?2 0 2 4 6 8 10 12 14 16 18 x 10 ?4 Normalised Axial Position Displacement, [mm ] Axial Displacement Radial Displacement Figure 5 -15: Displacement Variation Along Length of Tube at the Material Inter- face, (r = 0.5) 69 10 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 11 ?2 0 2 4 6 8 10 12 14 16 18 x 10 ?4 Radial Position Displacement, [mm ] Axial Displacement Radial Displacement Figure 5 -16: Displacement Variation Through the Tube Thickness at the Free End, (z = L) 5.7.1 Remarks While it was not the intention of this work to study free-edge effects, they can never the less not be avoided. It was shown that for a homogeneous orthotropic tube the equilibrium conditions were satisfied adequately. However, if the tube is non-homogeneous, complexities are introduced via edge effects. Due to the formulation of the numerical analysis, at the location of a material interface and a free-edge, a stress singularity will be predicted. The Rayleigh-Ritz technique as utilised in this work was not formulated to account for extremely large stress gradients (such as those expected at a singularity), as discussed earlier, the more complex the displacements, the more descriptive capability is required for the assumed displacement functions. Thus, singular behaviour does not lend itself to the chosen polynomial expansions. It is therefore expected that the Rayleigh-Ritz analysis presented in this work would not describe this problem particularly well. However, results were presented in this section that indicated that the results obtained from the Rayleigh-Ritz analysis were in good agreement to existing results (including those from a FEM analysis). In fact, the results could be said to be at least as good as both the FEM 40 and the results presented by Ye and Sheng39. In addition, although the analysis utilized to obtain these results was fairly refined (in comparison to some of the other analyses presented in this work), the total DOF?s 70 was not extremely large. Therefore, the final solution possesses a good efficiency. In general, a large portion of the investigations into the free-edge effect have been conducted by use of FEM. As the computer processing power has increased in re- cent times, so the refinement of the FEM models has increased. Thus at present, investigations into free-edge effects make use of FEM models with extremely refined meshes in the area of the free-edge41?43. In cases such as these, the primary purpose of the work (and in turn, the FEM model) is to determine the stress pattern at the free-edge. It is, therefore, suggested that if free-edge effects are of primary importance to the analysis, use is made of FEM as the mesh density in present use is not feasible for the Rayleigh-Ritz model, although the Rayleigh-Ritz analysis would require fewer DOF?s to obtain a similar descriptive ability. However, in the general analysis where free-edge effects are present but not the subject of investigation, the Rayleigh-Ritz model provides an accurate solution away from the edge effects as well as a reasonable approximation in the area of the edge effects. 5.8 Convergence Assessment of the Rayleigh-Ritz and FEM approaches It was stated in ?1.2 that past research1,23 has found the Rayleigh-Ritz solution to be more efficient than the FEM approach. In fact, this was a motivation for use of the Rayleigh-Ritz technique above the use of a FEM approach. This statement must therefore be investigated. To accomplish this, an orthotropic tube under shrinkage was modelled by both approaches using varying Degrees-Of-Freedom (DOF?s). In order to evaluate the validity of the solutions, the stresses on the boundary of the tube were observed. From the equilibrium conditions it is known that in the absence of applied stress, the boundary stresses in question must be zero. Thus it is possible to determine the accuracy of the obtained results. The selected stresses therefore are: ? the axial and shear stress at the tube free end (?zz, ?rz|Z=1) ? the radial and shear stress at the tube inner wall (?rr, ?rz|R=0) ? the radial and shear stress at the tube outer wall (?rr, ?rz|R=1) In order to best determine the analysis convergence, the maximum stress was selected at the desired locations for both the FEM and Rayleigh-Ritz solutions, i.e. the maximum axial stress at the free end was selected. Since in this case all the selected stresses should be zero, this represents the maximum error of the solutions. For the analysis, a similar tube to that utilised in ?5.5 was modelled, except that a tube 71 of length 1000 mm was modelled instead of a length of 2000 mm (as in ?5.5). The results are presented in Fig. 5 -17 to Fig. 5 -22. In order to visualise the data, trend lines have tentatively been fitted to the discrete data points. As can be seen, both approaches obtain more accurate solutions as the DOF?s in- crease. It is shown that for all cases the convergence of the Rayleigh-Ritz model is far quicker than the FEM model. Thus it can be said that the Rayleigh-Ritz model is more efficient than a FEM solution. For the FEM model, each node has 6 DOF?s. However, the assumptions of axisym- metry remove the 3 DOF associated with out-of-plane translation and rotations about angular and radial axes. The Rayleigh-Ritz analysis?s descriptive ability is directly related to the order of the displacement functions. Therefore the DOF?s can be calculated according to the relationship: DOF = [(N + 1)(M + 1) + (P + 1)(Q+ 1)]NmNl Which essentially yields the size of the displacement coefficient vector. 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Degrees of Freedom ? zz FEM Rayleigh?Ritz FEM Trendline Rayleigh?Ritz Trendline Figure 5 -17: Stress Residuals for the Axial Stress at the Tube Free End 72 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Degrees of Freedom ? r z FEM Rayleigh?Rtiz FEM Trendline Rayleigh?Ritz Trendline Figure 5 -18: Stress Residuals for the Shear Stress at the Tube Free End 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Degrees of Freedom ? rr FEM Rayleigh?Ritz FEM Trendline Rayleigh?Ritz Trendline Figure 5 -19: Stress Residuals for the Radial Stress at the Tube Inner Wall 73 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Degrees of Freedom ? r z FEM Rayleigh?Ritz FEM Trendline Rayleigh?Ritz Trendline Figure 5 -20: Stress Residuals for the Shear Stress at the Tube Inner Wall 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.5 1 1.5 2 2.5 Degrees of Freedom ? rr FEM Rayleigh?Ritz FEM Trendline Rayleigh?Ritz Trendline Figure 5 -21: Stress Residuals for the Radial Stress at the Tube Outer Wall 74 0 500 1000 1500 2000 2500 3000 3500 4000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Degrees of Freedom ? r z FEM Rayleigh?Ritz FEM Trendline Rayleigh?Ritz Trendline Figure 5 -22: Stress Residuals for the Shear Stress at the Tube Outer Wall 5.9 On The Use of Rayleigh-Ritz vs. FEM The Rayleigh-Ritz solution and the Finite Element Method possess many similari- ties, thus it is of primary importance to understand the differences between the two techniques to evaluate the merit of this work. The results presented in this work show that the Rayleigh-Ritz analysis is able to predict all the requirements specified in the objectives of the work. Additionally, the model is found to be efficient in solution. The results presented in ?5.8 clearly show that the Rayleigh-Ritz model produces solutions with much fewer DOF?s than a FEM model with comparable accuracy. It must be understood that although a discretisation scheme was implemented in this work, the underlying philosophy is to utilise as few discretised regions as possible to adequately describe the problem. The reason behind this is simple; suppose that a system undergoes a displacement which follows a cubic equation. This problem can be satisfied obviously by a single cubic displacement function. However, two regions of quadratic order are necessary to describe this behaviour, and except in specific cases, the cubic function will still not be described exactly. Without consideration for the additional constraint matrices required by the discretised approach, the result is that an extra DOF is required to accurately capture the behaviour. Thus the model is at its most efficient when utilising as high a series order as possible; which in turn implies as few discretisations as possible. Of course, if the discretisation is extended to the nth degree it can be argued that the Rayleigh-Ritz model would be able to duplicate the analysis of a Finite Element 75 Code. While in theory this is possibly true, in practice this is not the case. Com- mercial FEM software is formulated and implemented with the expectation that it will be utilised to solve problems with very large DOF?s (for example the Nastran code which is particularly adept at handling large problems). Thus the software packages are optimised to this end. The Rayleigh-Ritz model as presented in this work obviously includes none of these optimisations. Consequently the Rayleigh- Ritz analysis will be able to handle considerably fewer DOF?s than the FEM solution due to system limitations. Conversely, it is possible to perform the Rayleigh-Ritz analysis with far fewer dis- cretisations than that required by a FEM approach. Due to the formulation of the FEM, a limitation is placed upon the element aspect ratio (which should be less than about 3:1). This is a rigid constraint and can result in the need to utilise a large number of elements to attain an accurate solution. As can be seen from the formu- lation of the Rayleigh-Ritz analysis, there is no such limit in this approach; rather, the discretisations are simply utilised to better describe the system behaviour. As an illustration, consider a tube with wall thickness of 8 mm and length of 1000mm. For the sake of argument, the tube consists of 4 layers of material of equal thickness. If this tube is to be modelled with FEM it is desired that there is a fairly fine mesh density in the radial direction in order to capture the through-thickness behaviour. Thus, assuming it is possible that 2 quadratic elements are sufficient per layer of material, the radial edge length of these elements is 1mm. In order for the FEM model to give accurate results, the element aspect ratio should be at most 3:1, thus the maximum edge length in the axial direction can be 3mm. This leads to the necessity to model approximately 333 elements in the axial direction. This results in an extremely large total DOF for the problem. Consequently, if the same tube is analysed by means of the Rayleigh-Ritz approach, it is necessary to include only one region per material layer in the radial direction with an appropriate descriptive capability (in order for similar accuracy to the FEM model, a fourth order polynomial would suffice). Additionally the limitation on aspect ratio is far less strict than FEM, thus it may be necessary to include only 6 axial regions in the simulation, as was the case with a similar tube modelled in ?6.7. Thus the total DOF of the simulation will be far less than FEM. Obviously, changes in the geometric parameters can aggravate this problem severely (increases in length, decrease in layer thickness, etc.). This situation illustrates the potential advantages to utilising the Rayleigh-Ritz solution over FEM. At this point it can be argued that it is possible to implement some form of mesh refinement policy whereby a very coarse mesh is present away from the point of interest and is refined to attain a fine mesh in the region of interest. However such claims will then add credence to the arguments initially set out in the introduction. Such forms of mesh refinement require a fair amount of time to implement. Additionally once the mesh in completed, the results will only be able to be studied at selected areas. Thus if other areas are to be studied, the model will require to be re-meshed, committing more time spent in analysis. By utilising the Rayleigh-Ritz model, all areas within the tube can be studied with similar accuracy 76 and in the event that an area requires more descriptive capability it is a trivial task to modify the model and re-perform the analysis. With the above argument in mind, it is felt that there is sufficient merit in the present Rayleigh-Ritz analysis to justify it?s use on modelling multi-layered, finite length tubes. 77 6 Results In order to illustrate the full ability of the analysis presented in this work, some selected results are presented in this section. 6.1 Analysis of Ring Stiffened Cylinders A situation which occurs frequently in reality is that of a long tube with supporting frames along its length, generally referred to as ring-stiffened cylinders. While the tube is long, the presence of the supporting frames means that it is not possible to accurately predict the response of these tubes with the assumption of infinite length. The analysis in this work is, however, capable of describing them. To illustrate this, the response of two such tubes is calculated within this section. The tubes possess the same nominal dimensions but possess different material make- ups. These tubes are envisaged to perform the same function, that of transporting a fluid at elevated temperature and under pressure. Ring stiffeners are present at regular intervals along the length. The results within this section indicate the thermal response of both tubes in ad- dition to the stress states present in the tubes. This solution illustrates a number of points. Firstly, the ability of the analysis to predict such problems (where finite length, pressure loading and thermal loading are combined). Secondly, the effect of the different material make-up on the tube performance can be assessed. This may give an indication as to how to design such vessels. 6.2 Predicting the Response of Ring Stiffeners Ring stiffeners are an important structural element and there are a number of tech- niques for incorporating them into a mechanical analysis. Some of these techniques are discussed below. The frames can essentially be described as an additional tubular section in the model. While this has the advantage of accurately predicting the response and stiffness of these elements, it requires considerable computational effort. Thus it is only truly advantageous when the stress in the frame is of interest. As this is not the case in this analysis, this method is not feasible. An alternative is to specify the frame by enforcing a displacement constraint. How- ever, specifying zero displacement implies that the frame is infinitely stiff which is typically far from reality. For approximate analyses this may be sufficient (if say, the only points of interest were far away from the ring stiffener), but this analysis does not fall into this category. Additionally, displacement constraints impose some 78 discontinuity in the model thereby introducing difficulty in modelling the displaced shape with polynomial functions. An additional alternative is to replace the ring stiffener with an equivalent exter- nal pressure thereby modelling the interface pressure between the elements. Thus, there are no sharp discontinuities in the model and there is no need to model the ring stiffener itself. This approach is therefore computationally efficient whilst still preserving the elasticity of the frame. This technique, in application, is equivalent to the methods by which band loads are applied to cylinders 17,21 and it is the chosen approach for this analysis. 6.3 Modelling Ring Stiffened Cylinders The tubes under consideration are long tubes with supporting frames of the type shown in Fig. 6 -1. It is possible to model this tube by a representative region midway between the stiffeners indicated in Fig. 6 -1. Assuming that the ring stiffen- ers prevent axial motion of the tube (as is often encountered), symmetry boundary conditions can be placed at both ends of the representative regionxxiv. Figure 6 -1: Long Tube with Supporting Frames at Regular Intervals Along the Length By virtue of the modelling techniques described in the previous section, the trans- formation from the tube representative region into the analysis region is presented in Fig. 6 -2. Specifically, the boundary conditions and the interface pressure between the ring stiffener and the cylinder is shown in Fig. 6 -2. xxivWith these boundary conditions, it is also possible to consider a representative region between the ring stiffeners (i.e. from one ring stiffener to the next). Additionally, the representative region could be further split in two along the length by virtue of symmetry if desired. 79 Figure 6 -2: Transformation from Representative Region into Analysis Region 6.4 Material Properties The tubes are constituted from four materials, the properties of which are presented in Table 6 -1. The abbreviations utilised in this table represent: CB - Corrosion Barrier CSM - Chopped Strand Mat WR - Woven Roving FW - Filament Winding The corrosion barrier is present to provide protection for the main structural fibres from chemical corrosion. This layer is typically resin rich, often a resin with glass veil (which is an extremely low areal mass glass fabric). However, even if glass veil is present, the glass content is of sufficiently low volume that the corrosion barrier can be treated as isotropic without any significant loss in accuracy. The chopped strand mat and woven roving are in the form of pre-produced material layers (as opposed to rolls of fibre) which have been laid up around the tube. In reality this is actually a difficult task to accomplish whilst maintaining the cylindrical nature of the tube. In this analysis however, it will be assumed that the final tubes are perfectly cylindrical and symmetrical about the angular axis. The latter assumption is necessary as the random nature of the chopped strand mat fabric is ignored and the properties are assumed to be uniform. 80 Filament winding typically produces the best result for creating cylinders as fibres are wound onto a mandrel. As the process can be mechanised, a large deal of controllability and repeatability exists with this technique. It is commonly the technique utilised to create composite tubing. Usually the fibres are wound at an angle of +/ ? 55? to the axial axis of the tube. The reason for this is that by winding at this angle, the tube possesses a circumferential-to-axial strength ratio of approximately 2:1, which is, of course, the ratio between the hoop and axial stresses in thin-walled pressure vessel theory. The material properties of the corrosion barrier are based on those of the resin, the references for these properties can be found in App. C. The stiffness properties of the chopped strand mat are obtained from test data 44. The thermal properties of the chopped strand mat (i.e. ?, k, c) were based on calculation. The properties of the woven roving and the filament wound layers were also based on calculation. In order to illustrate the techniques utilised to estimate the material properties in Table 6 -1, a sample calculation of the woven roving is presented in App. C. The same techniques are also utilised for the other materials. Table 6 -1: Material Properties for Ring Stiffened Cylinder Simulation Property CB CSM WR FW Ezz (GPa) 3.5 9.67 15.04 7.97 E?? (GPa) 3.5 9.67 15.04 16.23 Err (GPa) 3.5 6.12 4.96 8.16 Grz (GPa) 1.3 1.98 1.84 3.05 ?z? 0.35 0.34 0.101 0.415 ?zr 0.35 0.39 0.355 0.209 ??r 0.35 0.39 0.355 0.054 ?z (?C?1) 65? 10?6 26.79? 10?6 20.4? 10?6 23.9? 10?6 ?? (?C?1) 65? 10?6 26.79? 10?6 20.4? 10?6 4.8? 10?6 ?r (?C?1) 65? 10?6 78.42? 10?6 69.6? 10?6 42.3? 10?6 k (W/m ?K) 0.1780 0.2154 0.2697 0.4561 c (J.kg?1) 1400 1223 1103 945 ? (kg.m?3) 1126 1352 1564 1974 81 6.5 Tube Geometric Properties and Loading Conditions As discussed earlier, both the tubes which are modelled in this section possess the same nominal dimensions. The representative analysis section is sketched in Fig. 6 -3. As can be seen, the tube length is 1000 mm. The inner radius is given as 100 mm, with a total thickness of 8 mm. The individual make-ups of the tubes are discussed in the following sections. With regard to the thermal loading, the tubes are subject to an inner wall temper- ature of 75 ?C above ambient temperature. As the tubes possess no insulation, they are allowed to dissipate heat to the atmosphere by free convection with a coefficient of hc = 10 W/m2K as indicated in Fig. 6 -3. The tubes are subject to an internal pressure of 10MPa. The ring stiffener is located symmetrically about 500mm along the length. The ring stiffener applies a pressure of 5 MPa across 10 mm about this position as shown in Fig. 6 -3. Figure 6 -3: Schematic of Ring Stiffened Tube 82 6.5.1 Filament Wound (FW) Tube Beginning at the inner wall, the filament wound (FW) tube consists of a 2 mm corrosion barrier. The tube then has a 2 mm layer of chopped strand mat. The final layer consists of 4 mm of filament wound fibre. The make-up is illustrated graphically in Fig. 6 -4. Figure 6 -4: Make-Up of Filament Wound Tube 6.5.2 Chopped Strand Mat (CSM) Tube Beginning at the inner wall, the chopped strand mat (CSM) tube consists of a 2mm corrosion barrier. The following layer is a 2mm layer of chopped strand mat. A layer of woven roving of 0.5 mm thickness follows. The woven roving is orientated such that the fibres run along the length of the tube and in the circumferential direction. The final layer of the tube is a 3.5 mm layer of chopped strand mat. The make-up of this tube is illustrated graphically in Fig. 6 -5. Figure 6 -5: Make-Up of Chopped Strand Mat Tube 83 6.6 Results for the FW Tube 6.6.1 Thermal Response of FW Tube The thermal response of the filament wound tube is presented in this section. The simulation possessed the following parameters: nCB 30 nCSM 30 nFW 60 ?T 0.01s Where ni is the number of nodes in material layer i. The result is presented in Fig. 6 -6. 0.1 0.101 0.102 0.103 0.104 0.105 0.106 0.107 0.108 0.109 0.11 0 10 20 30 40 50 60 70 80 Radial Position, [m] Temperature Change, ? T [ ? C ] t = 50s t = 100s t = 150s t = 200s t = 250s t = 500s t = ? Figure 6 -6: Through Thickness Temperature Distribution - FW Tube It is interesting to note that at steady state, the tube temperature profile varies with the radial coordinate. Another significant fact is the clear change in curvature of the temperature profile at the material interface, particularly at the filament winding interface. For purposes of implementing the temperatures profiles in Fig. 6 -6 in the mechanical analysis, the temperature profiles within each material layer were represented by a least squares fit polynomial. Even with low order polynomials (in this case, cubic functions) a good correlation was obtained; in this particular case a correlation coefficient of near unity was attained, (i.e. r2 = 0.99999). Additionally, as each of 84 the material layers is represented in the mechanical analysis by a radial coordinate which ranges from 0 to 1, the polynomial fits can be conducted with ease without having to account for the physical radial dimensions. 6.6.2 Mechanical Results of FW Tube The results from the mechanical model for the filament wound tube are presented in the following sections. The results are presented for the steady state temperature profile. Displacement and stress results are presented at selected positions of interest throughout the tube. This solution was conducted with the following parameters: N 5 M 5 P 9 Q 5 Nl 12 Nm 3 The displacement results are presented in ?6.6.2, followed by the stress states along the length of the tubes length in ?6.6.2 and the stress states through the thickness in ?6.6.2. 6.6.2.1 Displacement Results of FW Tube Fig. 6 -7 shows the displacement variation along the length of the tube, calculated along the interface between the CSM layer and the FW layer (see Fig. 6 -4) which is halfway through the thickness of the tube. Details of Fig. 6 -7 showing the areas of interest (in the vicinity of the ring stiffener) in more detail are presented in Fig.?s 6 -8 and 6 -9. 85 0 100 200 300 400 500 600 700 800 900 1000 ?0.2 0 0.2 0.4 0.6 0.8 1 1.2 Axial Position, [mm] Displacement, [mm ] Axial Displacement Radial Displacement Detail 2 Detail 1 Figure 6 -7: Displacement Variation Along Tube Length - FW Tube (r = 104mm) The first result that is evident from Fig.?s 6 -7 - 6 -9 is that the tube symmetry conditions are indeed satisfied. Clearly in Fig. 6 -8, the axial displacement is zero at the tube midplane. Additionally, in Fig. 6 -7 the axial displacement is zero at both ends of the tube. Fig. 6 -9 indicates that the axial displacement is antisymmetric about the tube midplane as would be expected in these conditions. Further, the radial displacement is also symmetric about the tube midplane as can be seen in Fig. 6 -7. From these results, the presence of the ring stiffener is clearly evident in the radial displacement variation along the length of the tube. 86 400 420 440 460 480 500 520 540 560 580 ?0.06 ?0.04 ?0.02 0 0.02 0.04 Axial Position, [mm] Displacement, [mm ] Axial Displacement Figure 6 -8: Displacement Variation Along Tube Length - FW Tube, Detail 1 (r = 104mm) 400 420 440 460 480 500 520 540 560 580 0.84 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 Axial Position, [mm] Displacement, [mm ] Radial Displacement Figure 6 -9: Displacement Variation Along Tube Length - FW Tube, Detail 2 (r = 104mm) 87 6.6.2.2 Stress Variation with Length - FW Tube The stress variation along the length corresponding to the displacements in Fig. 6 -7 is presented in Fig. 6 -10. Details of Fig. 6 -10 in the region of the ring stiffener are presented in Fig.?s 6 -11 and 6 -12. From Fig.?s 6 -10 - 6 -12 it can be seen that, as expected, all the stress components exhibit symmetrical behaviour about the tube midplane, except the shear stress which is antisymmetrical. The antisymmetrical behaviour of the shear stress will be discussed further in ?7. As is expected, the stress magnitude over the majority of the tube is constant, but the stress varies considerably in the vicinity of the ring stiffener. Of all the stress components, the circumferential stress is the greatest in magnitude. Of course, as the tube is under considerable internal pressure, this behaviour is also expected. As the tube essentially does not possess ?capped? ends, there is no axial component of stress arising from the internal pressure, therefore the axial stress is small in magnitude. 0 100 200 300 400 500 600 700 800 900 1000 ?20 0 20 40 60 80 100 Axial Position, [mm] Stress, [MPa ] ?zz ??? ? rr ? rz Detail 1 Detail 2 Figure 6 -10: Stress Variation Along Tube Length - FW Tube, (r = 104mm) 88 400 420 440 460 480 500 520 540 560 580 600 62 64 66 68 70 72 74 Axial Position, [mm] Stress, [MPa ] ??? Figure 6 -11: Stress Variation Along Tube Length - FW Tube, Detail 1 (r = 104mm) 400 420 440 460 480 500 520 540 560 580 600 ?15 ?10 ?5 0 5 10 15 Axial Position, [mm] Stress, [MPa ] ? zz ? rr ? rz Figure 6 -12: Stress Variation Along Tube Length - FW Tube, Detail 2 (r = 104mm) 89 6.6.2.3 Through-thickness Stress Variation - FW Tube Fig. 6 -13 presents the stress variation through the thickness of the tube at the edge of the ring stiffener (z = 495mm). A detail of Fig. 6 -13 is presented in Fig. 6 -14 which shows the behaviour of the radial and shear stresses in larger scale. The through thickness strains at this position are shown in Fig. 6 -15. As there is a discontinuity in loading at the boundary of the ring stiffener, the stresses are calculated for the region bordering the ring stiffener (i.e. just outside the ring stiffener). From Fig. 6 -14 it can be seen that the radial stress is equal to ?10MPa at the inner wall and returns to a value near zero at the outer wall. The shear stress begins and ends approximately at zero. As expected, the stresses ?zz and ??? are discontinuous at the material boundaries as a consequence of the different material stiffnesses. In general, the value of the radial stress at the tube thickness extremities is a measure of the solution accuracy. As the stresses in Fig.?s 6 -13 - 6 -14 are calculated just before the ring stiffener, the radial stress should be equal to zero at the outer wall. However, in this case there are additional factors which need to be considered. The antisymmetric behaviour of the shear stress in the vicinity of the ring stiffener affects the behaviour of the radial stress in this region, to the point whereby the radial stress will not attain the exact solution. Therefore, the fact that the radial stress is not zero at the outer wall in the vicinity of the ring stiffener can not be used to infer anything regarding the simulation accuracy. These effects are discussed further in ?7. In looking at Fig. 6 -14 it can be seen that the radial and shear stress are essen- tially continuous between the material layers, as is required by equilibrium. Small discontinuities in these stresses can be seen, however, they are well within the limits of accuracy for this simulation. 90 100 101 102 103 104 105 106 107 108 ?50 0 50 100 150 200 Radial Position, [mm] Stress, [MPa ] ? zz ??? ? rr ? rz Figure 6 -13: Stress Variation Through the Tube Thickness - FW Tube (z = 495mm) 100 101 102 103 104 105 106 107 108 ?20 ?15 ?10 ?5 0 5 10 15 20 Radial Position, [mm] Stress, [MPa ] ? rr ? rz Figure 6 -14: Stress Variation Through the Tube Thickness - FW Tube, Showing only Radial and Shear Stress, (z = 495mm) 91 100 101 102 103 104 105 106 107 108 ?6 ?4 ?2 0 2 4 6 8 x 10 ?3 Radial Position, [mm] Strain [mm/mm ] ? zz ??? ? rr ? rz Figure 6 -15: Strain Variation Through the Tube Thickness - FW Tube, (z = 495mm) 6.7 Results for the CSM Tube 6.7.1 Thermal Response of CSM Tube The thermal response of the CSM tube is presented in this section. The simulation possessed the following parameters: nCB 30 nCSM1 30 nWR 8 nCSM2 50 ?T 0.01s The result is presented in Fig. 6 -16. As with the FW tube, the temperature profile at steady state in the CSM tube also varies with the radial coordinate. In this case, the low fibre volume fraction of CSM means that the slope change in the temperature profile at the material interfaces is barely perceptible. 92 0.1 0.102 0.104 0.106 0.108 0.11 0 10 20 30 40 50 60 70 80 Radial Position, [m] Temperature Change, ? T [? C ] t = 50s t = 100s t = 150s t = 200s t = 250s t = 500s t = ? Figure 6 -16: Through Thickness Temperature Distribution - CSM Tube 6.7.2 Mechanical Results of CSM Tube The results from the mechanical analysis for the chopped strand mat tube are pre- sented in the following sections. Again, results are presented for the steady state temperature profile. Displacement and stress results are presented at selected posi- tions of interest throughout the tube. This solution was conducted with the following parameters: N 5 M 5 P 7 Q 5 Nl 12 Nm 4 The displacement results are presented in ?6.7.2, followed by the stress states along the tubes length in ?6.7.2 and the stress states through the thickness in ?6.7.2. 6.7.2.1 Displacement Results of CSM Tube The displacement variation along the tube length is presented in Fig. 6 -17. For these results, the displacements were calculated at the outer material interface, i.e. r = 104.5mm. Details of Fig. 6 -17 in the region of interest around the ring stiffener are presented in Fig.?s 6 -18 and 6 -19. 93 0 100 200 300 400 500 600 700 800 900 1000 ?0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Axial Position, [mm] Displacement, [mm ] Axial Displacement Radial Displacement Detail 2 Detail 1 Figure 6 -17: Displacement Variation Along Tube Length-CSM Tube,(r = 104.5mm) 350 400 450 500 550 600 650 ?0.1 ?0.05 0 0.05 0.1 Axial Position, [mm] Displacement, [mm ] Axial Displacement Figure 6 -18: Displacement Variation Along Tube Length - CSM Tube, Detail 1 94 350 400 450 500 550 600 650 1.45 1.5 1.55 1.6 1.65 1.7 Axial Position, [mm] Displacement, [mm ] Radial Displacement Figure 6 -19: Displacement Variation Along Tube Length - CSM Tube, Detail 2 As can be seen from Fig.?s 6 -17 - 6 -19, the displacements in the CSM tube ex- hibit similar behaviour to that of the FW tube. Again, all the required symmetry conditions are satisfied. 6.7.2.2 Stress Variation with Length - CSM Tube Fig.?s 6 -20 - 6 -22 represent the stress variation along the tube length at different positions through the tube thickness. The stress variation along the inner wall of the tube can be seen in Fig. 6 -20. As expected, the radial stress is constant and reacts the applied internal pressure and the shear stress is zero. It is interesting to note that at the inner wall, even though it is relatively far from the point of application of the external pressure (i.e. the ring stiffener interface), there is still a large variance in stress magnitude in the vicinity of the ring stiffener. The stress variation along the length at the outer material interface (corresponding to the displacements shown in Fig. 6 -17) is presented in Fig. 6 -21. These results can be compared to the results for the FW tube in Fig. 6 -10 in terms of behaviour trends as, in both cases, the results are approximately halfway through the thickness and at a material interface. In terms of behavioural trends, the results in Fig. 6 -21 exhibit very similar behaviour to those presented in Fig. 6 -10. Again, all the stress components, barring the shear stress, are symmetrical about the tube midplane. As mentioned earlier, the behaviour of the shear stress is discussed in ?7. 95 The global stress state at the outer wall is presented in Fig. 6 -22. As can be seen, there are noticeable discontinuities in the axial and circumferential stresses at the extremities of the external pressure band. In the previous section it was mentioned that the antisymmetric behaviour of the shear stress effects the behaviour of the other stress components in the vicinity of the ring stiffener such that the exact solution is not attained. The cause behind the discontinuities in Fig. 6 -22 is exactly the same. A detail of Fig. 6 -22 is presented in Fig. 6 -23 showing, in more detail, the behaviour of the radial and shear stresses in the vicinity of the ring stiffener. As is evident, relatively large discontinuities exist in the shear stress at this scale. Additionally, the radial stress oscillates about the required equilibrium condition (in this case, 5MPa). Due to the manner in which the analysis is formulated, the antisymmetric behaviour of the shear stress is expected; however, while this shear stress behaviour will be attained, it does not follow the correct solution (which predicts a zero shear stress along the outer surface). Therefore, as the stress components are inter-dependent, it is unlikely that if one of the stress components does not follow the exact solution, that the other stress components will be able to. Thus, the shear stress behaviour essentially prevents the other stress components from attaining the continuity which is required. Again, the issue is addressed at length in ?7. It can be seen that this discontinuous behaviour in shear stress only has a local effect on the final results as no significant discontinuities can be noted at the scale of the stress plots in Fig.?s 6 -20 and 6 -21. 0 100 200 300 400 500 600 700 800 900 1000 ?20 ?10 0 10 20 30 40 Axial Position, [mm] Stress, [MPa ] ? zz ??? ? rr ? rz Figure 6 -20: Stress Variation Along Tube Length - CSM Tube, (r = 100mm) 96 0 100 200 300 400 500 600 700 800 900 1000 ?50 0 50 100 150 200 250 Axial Position, [mm] Stress, [MPa ] ? zz ??? ? rr ? rz Figure 6 -21: Stress Variation Along Tube Length at the Outer Material Interface - CSM Tube, (r = 104.5mm) 0 100 200 300 400 500 600 700 800 900 1000 ?20 0 20 40 60 80 100 120 140 160 Axial Position, [mm] Stress, [MPa ] ? zz ??? ? rr ? rz Detail 1 Figure 6 -22: Stress Variation Along Tube Length at the Outer Wall - CSM Tube, (r = 108mm) 97 475 480 485 490 495 500 505 510 515 520 525 ?8 ?6 ?4 ?2 0 2 4 Axial Position, [mm] Stress, [MPa ] ? rr ? rz Figure 6 -23: Stress Variation Along Tube Length - CSM Tube, Detail 1 (r = 108mm) 6.7.2.3 Through-thickness Stress Variation - CSM Tube The stress state through the thickness of the tube at selected locations is presented in Fig. 6 -24 - 6 -27. The strain through the thickness is presented in Fig. 6 -27. Fig. 6 -24 presents the stress through the thickness at the interface of the discrete pressure. As with previous results, Fig. 6 -24 represents the stress state at a point just outside the external pressure band. Fig. 6 -25 presents a detail of Fig. 6 -24 showing the radial and shear stress variation in more detail. Again, equilibrium dictates that the radial stress is equal to ?10MPa at the in- ner wall and zero at the outer wall. From Fig. 6 -25 it can be seen that this is approximately the case. As in the case of the FW tube, the radial stress does not return exactly to zero at the outer wall. This occurs for exactly the same reason as in the FW tube and the discussion relating to this phenomenon can be found in ?7. As with the FW tube, the stresses ?zz and ??? are discontinuous at the material boundaries. This is expected as the material layers exhibit the same strains in these directions at the material interface (as required by compatibility) but the material stiffnesses are different. 98 The stress state at the symmetry plane (i.e. either end of the tube) is presented in Fig. 6 -26. A detail presenting the shear and radial stress in more detail is presented in Fig. 6 -27. At this location it can be seen that the radial and shear stresses do satisfy the required equilibrium conditions; they are not affected by the local effects of the ring stiffener. 100 101 102 103 104 105 106 107 108 ?20 0 20 40 60 80 100 120 140 160 Radial Position, [mm] Stress, [MPa ] ? zz ??? ? rr ? rz Figure 6 -24: Stress Variation Through the Thickness - CSM Tube, (z = 495mm) 99 100 101 102 103 104 105 106 107 108 ?10 ?5 0 5 10 Radial Position, [mm] Stress, [MPa ] ? rr ? rz Figure 6 -25: Stress Variation Through the Thickness - CSM Tube, Showing only Radial and Shear Stress, (z = 495mm) 100 101 102 103 104 105 106 107 108 ?50 0 50 100 150 200 250 Radial Position, [mm] Stress, [MPa ] ? zz ??? ? rr ? rz Figure 6 -26: Stress Variation Through the Thickness - CSM Tube, (z = 0mm) 100 100 101 102 103 104 105 106 107 108 ?12 ?10 ?8 ?6 ?4 ?2 0 2 Radial Position, [mm] Stress, [MPa ] ? rr ? rz Figure 6 -27: Stress Variation Through the Thickness - CSM Tube, Showing only Radial and Shear Stress, (z = 0mm) 6.8 Transient Thermal Stresses in Composite Tubes As discussed in the introduction, transient effects play an important part in the thermal behaviour of composite tubes. With this in mind, the transient thermal behaviour of the two tubes utilised in the previous section is analysed in this section. The exact details of the analysis are discussed in the following sections. 6.8.1 Tube Geometric Parameters The tubes considered in this analysis are essentially the same tubes as those utilised in the previous section. The only modification is that only the thermal loading is considered, i.e. the pressure loadings are not present. In terms of the assumptions of the previous analysis, this implies that the tubes no longer represent ring-stiffened cylinders. In deference to the previous analysis where the tubes were considered to be long with periodic ring stiffeners, this analysis is conducted on a finite length tube of total length 1000 mm with free ends. Due to symmetry, half of this tube was modelled with the appropriate boundary conditions, i.e. a tube of length 500 mm was analysed. 101 6.8.2 Loading Parameters In the previous section the tubes were subjected to an elevated temperature, i.e. thermal expansion. In this analysis the tubes are subject, rather, to thermal con- traction. Rather than re-performing a thermal analysis, the results obtained in the previous section can easily be modified for this purpose. Therefore, consider the case of the tubes transporting fluid which, rather than being at an elevated temper- ature to the ambient conditions (wherein the tubes are assumed to be stress-free) are transporting fluid at a temperature below ambient. In this case, heat will be transferred to the tube from the environment as opposed to the reverse situation. However, as all the thermal properties of the system are exactly the same as those in the previous thermal analysis, the tube behaviour will be similar. In fact, the thermal response of the tubes in this case will be exactly equivalent to that shown in the previous analysis (see Fig.?s 6 -6 and 6 -16) mirrored, as it were, about the line ?T = 0?C. 6.8.3 Transient Stress Variation Results In general, the through-thickness stress state is critical at the material interfaces, particularly when considering the inter-laminar normal stress. With this in mind, the transient behaviour of the radial stress in the tubes is compared in Fig. 6 -28 at the outer material interface for both tubes. The outer material interface was selected as it allows the investigation at a material interface in addition to being at a similar position through the thickness in both tubes. This is therefore a good position to draw a comparison between the tubes. In order to compare the stresses, they are obtained midway along the length of the tubes. From Fig. 6 -28 it is clearly evident that the radial stress in both tubes is small in magnitude. Consequently, the variation exhibited by both tubes is small in mag- nitude. Thus these results might question the importance of transient effects in situations such as this. However, the variation in relative terms is significant, with the FW tube exhibiting a 28.2% change in magnitude over the time period while the CSM tube exhibits a 29.7% change. This, therefore, indicates that although the magnitude of the stress is not large, there is a substantial variation (in relative terms) with time. 102 50 100 150 200 250 300 350 400 450 500 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time, [s] Radial Stress, ? rr , [MPa ] ? rr , FW Tube ? rr , CSM Tube FW Tube Trendline CSM Tube Trendline Figure 6 -28: Transient Radial Stress Variation at the Axial Midplane for the CSM and FW Tubes It should be recognised that the tubes under consideration in this section possess a free end. It can therefore be expected that the free edge effect will be evident in the tubes. As discussed earlier, the result of primary importance in the presence of free edge effects is the normal inter-laminar stress (which exhibits singular behaviour at the free edge). This high stress near the free edge is thought to be the driving factor behind delamination failures. Now, as discussed earlier, the results obtained from this analysis in the region of the free edge will not converge to the exact solution due to the stress singularities. They will, however, be sufficiently accurate so as to allow some behaviour patterns for the tubes to be inferred. The results of radial stress variation along the length at the outer material interface for the filament wound tube is presented in Fig. 6 -29. The corresponding results for the chopped strand mat tube are presented in Fig. 6 -30. Again, Fig.?s 6 -29 and 6 -30 show that the radial stress does not vary significantly in magnitude with time. However, the singular behaviour in the radial stress is clear. Noticeably, the radial stress approaches positive infinity with an almost vertical gradient near the free end of the tubes. The almost vertical behaviour shown in Fig.?s 6 -29 and 6 -30 is particularly interesting as this type of behaviour is usually difficult to obtain with polynomial functions. In fact, it is somewhat surprising that this analysis can predict such effects so well. Therefore, these results indicate the merits in employing the discretisation scheme. 103 0 100 200 300 400 500?2 0 2 4 6 8 10 12 Axial Position, [mm] Radial Stress, ? rr , [MPa ] T = 50 s T = 250 s T = 500 s T = ? Figure 6 -29: Transient Radial Stress Variation Along the Tube Length at the Outer Material Interface, FW Tube 0 100 200 300 400 500?1 0 1 2 3 4 5 6 7 Axial Position, [mm] Radial Stress, ? rr , [MPa ] T = 50 s T = 250 s T = 500 s T = ? Figure 6 -30: Transient Radial Stress Variation Along the Tube Length at the Outer Material Interface, CSM Tube 104 Up to this point, only the transient behaviour of the radial stress has been presented. It was found that the radial stress exhibited a significant variation in relative terms across the time interval in question. Therefore it would be expected that this be- haviour would also be evident in the other stress components. To this end, the transient behaviour of the axial and circumferential stresses of both tubes are com- pared in Fig.?s 6 -31 and 6 -32 respectively. As with Fig. 6 -28, these results are also generated at midpoint along the length of the tubes. 50 100 150 200 250 300 350 400 450 500 ?15 ?10 ?5 0 5 10 15 20 25 Time, [s] Axial Stress, ? zz , [MPa ] FW Tube ? Inner Wall CSM Tube ? Inner Wall FW Tube ? Outer Wall CSM Tube ? Outer Wall FW Tube ? Inner Wall Trendline CSM Tube ? Inner Wall Trendline FW Tube ? Outer Wall Trendline CSM Tube ? Outer Wall Trendline Figure 6 -31: Transient Axial Stress Variation at the Axial Midplane for the CSM and FW Tubes 105 50 100 150 200 250 300 350 400 450 500 ?20 ?15 ?10 ?5 0 5 10 15 20 25 Time, [s] Circumferential Stress, ? ?? , [MPa ] FW Tube ? Inner Wall CSM Tube ? Inner Wall FW Tube ? Outer Wall CSM Tube ? Outer Wall FW Tube ? Inner Wall Trendline CSM Tube ? Inner Wall Trendline FW Tube ? Outer Wall Trendline CSM Tube ? Outer Wall Trendline Figure 6 -32: Transient Circumferential Stress Variation at the Axial Midplane for the CSM and FW Tubes As can be seen in Fig.?s 6 -31 and 6 -32, the axial and circumferential stresses also undergo a substantial variation in magnitude in both the CSM and the FW tubes. In terms of the axial and circumferential stresses at the outer wall, both tubes exhibit very similar behaviour. It is interesting to note that although the magnitude of the radial stresses in both tubes are small (to the point where they could possibly be neglected), the other stress components within the tubes are significant, as indicated in Fig.?s 6 -31 and 6 -32. Additionally, although the tubes under consideration are not particularly thick (r/t > 10), there is a considerable difference between the magnitudes of the stresses shown in Fig.?s 6 -31 and 6 -32 for both tubes. Therefore the variation in stress through the tube thickness is also large. These facts indicate the merit in performing analyses of this type, as, even though the tubes may not be considered thick-walled, they still exhibit large through-thickness stress variations, with stresses of significant magnitude to influence the failure of the tubes. 106 7 Discussion 7.1 Ring Stiffened Tubes - Thermal Results As can be seen from Fig.?s 6 -6 and 6 -16, the thermal behaviour of the two tubes is very similar. This is expected for a number of reasons. Firstly, the through-thickness conduction co-efficient is determined by a rule-of- mixtures approach. As such, it is essentially a ratio between the matrix and fibre conduction coefficients based on the fibre volume fraction. Therefore, as all the materials in both tubes consist of the same glass and resin in conjunction with the fact that the volume ratios between the materials do not differ extremely, the out- of-plane conduction coefficients will not vary considerably. This fact, coupled with the fact that the tubes possess the same geometrical parameters implies that the thermal problems are fairly similar between the tubes. Therefore it is expected that the results exhibit some similarities. From a point of view of further analysis, this is convenient as it allows the stress states in the tubes to be compared directly as there is little difference between the loading of the two tubes. The main result from this section is the fact that at steady state the through- thickness temperature distribution is not constant. Thus the actual thermal loading differs considerably from the simple assumption of constant through-thickness tem- perature change. 7.2 Ring Stiffened Tubes - Mechanical Results In terms of the expected mechanical behaviour of the tubes, the results presented indicate that both tubes behave as expected. As noted, in both tubes the dominant in-plane stress is in the circumferential direction as a result of supporting the in- ternal pressure. In terms of magnitude, the stress in the filament winding is higher than that obtained in the chopped strand mat. As the filament winding possesses a higher circumferential stiffness than the chopped strand mat, this is expected. Consequently, the filament wound tube exhibits a lower radial displacement than that seen in the chopped strand mat. Importantly it must be noted that in both tubes, the radial displacements are within the generally accepted limits of validity for a linear analysis of vessels under pressurexxv. xxvIt is generally accepted that if the out-of-plane displacement of a plate under pressure exceeds a value equal to half the plate thickness, it is necessary to perform a nonlinear analysis33. While this is mainly to include membrane effects, it is safe to say that if the radial displacement of the pressure vessel is below half the thickness, the analysis is well within the limits for linear analyses. 107 Some small discrepancies between the radial and shear stress components can be seen in the through-thickness states of both tubes at the material interfaces. As a result of the C0 compatibility, there is no condition which states that the displacement slopes must be equal at the material interfaces. Therefore, as the stress is directly related to the displacements, a discontinuity in displacement slope will lead to a mismatch in stresses. This therefore becomes a measure of the solution accuracy, i.e., if the displacements are not accurate, the discrepancy in stress becomes more pronounced. However, from the results presented it can be seen that the mismatch shown is well within limits for an accurate solution. In looking at the results presented for the ring stiffened tubes, a large variation in stress can be seen to occur both along the length and through the thickness of the tubes. Obviously, this behaviour cannot be predicted by the ?constant axial strain? assumption. This indicates the increase in solution validity obtained by considering a finite length tube. The antisymmetric behaviour which was noted in the shear stress variation along the length in the tubes is addressed in the following section. 7.3 Stress Discontinuities In the initial assumptions utilised in the development of the mechanical analysis is the fact that the stress tensor is symmetric, i.e. ?rz = ?zr. This assumption is embodied in the manner in which the constitutive relations are written (See Eqn. 2), and consequently leads to a symmetric compliance matrix, C. In elementary solid mechanics, the relations of the stress components on a differential volume element are studied. From this situation it is apparent that, in order to attain moment equilibrium, the complementary shear stress components must be equal. It has, however, been shown12 that this does not have to necessarily be true in all instances; particularly, when a stress component derivative is unbounded (i.e. at a stress discontinuity or stress singularity). In this case the stress tensor may be unsymmetric whilst still satisfying all the equilibrium conditions. This has important repercussions for the present work as it is apparent that if a situation is encountered whereby an unbounded stress derivative is present, this analysis will still attempt to resolve the situation with a symmetric stress tensor. And hence, the exact solution may not be attained. This problem is not unique to only this mechanical analysis, as all standard, approximate numerical techniques utilise the same assumption of a symmetric stress tensor and hence, will possess the same problem. Obviously this also extends to the finite element method, which, it must be remembered, is frequently utilised for the analysis of free-edge effects. With specific regard to stress discontinuities, Raju et al. 12 discusses the problem of a discrete pressure applied to a semi-infinite plate, a case for which the exact elasticity solution is known16. It must be recognised that while this geometry differs 108 to that discussed in this work, and the exact form of the elasticity solutions may differ, both systems must exhibit similar behaviourxxvi. This problem is therefore analogous to the results presented in this work for the discrete external pressure in ?6 and can be thus utilised to explain some of the behaviour noted in these results, particularly the behaviour of the shear stress. In looking at the results for the ring stiffened tubes (See ?6.6.2 and ?6.7), some possibly unexpected behaviour can be noticed in the shear stress variation along the tube length (See for example Fig. 6 -12). It appears as though the shear stress exhibits some form of sawtooth behaviour. In fact, this behaviour is to be expected and will be found in any of the standard numerical techniques 12. The elasticity solution for a discrete pressure over a semi-infinite plate indicates that the shear stress must be zero along the edge of the plate. Consequently, for the case presented here, this would mean that the through-thickness shear, ?rz, must be zero along the outer wall of the tube. However, if the extremities of the pressure bands are approached (in a mathematical sense) from either side, it is found that a discontinuous limit exists at this interface for the shear stress. Raju et al.12 shows that this discontinuity therefore must be interpreted as an unsymmetric stress tensor. As mentioned earlier, numerical techniques cannot describe this discontinuity ex- actly; to this end Raju et al.12 investigates the behaviour of the finite element method in this region. The results show that the shear stress is actually predicted as possessing discontinuous behaviour. This behaviour is not unlike that shown in ?6. It is noted by Raju et al.12 that while the result does not correlate with the exact solution the integral of the shear stress across the pressure interval yields the correct result (in this case, zero). The same can be said for the results presented in this work. In looking at Fig. 6 -12 it can be seen that the shear stress is antisymmetric about the interval of applied pressure. Thus, the integral of this stress across this interval yields zero. It is conjectured by Raju et al.12 that the reason for this boundary being satisfied in an average sense lies with the manner in which the finite element is formulated. The finite element method is a displacement based solution, therefore boundary con- ditions are applied by enforcing nodal displacements. Now, from the relationships presented earlier, the displacements are essentially related to the stresses by means of an integral, i.e. the force/displacement problem is the integral of the stress/strain problem. Therefore it makes sense that the integral of the stress yields the correct boundary condition. Obviously, the Rayleigh-Ritz solution shares a commonality with the finite element method in that it is also a displacement based technique. Therefore this behaviour is expected to be typical of the mechanical analysis pre- sented in this work (as can be seen to be the case). As the numerical technique utilised in this work is unable to predict the exact response in this situation, the merit of including the through-thickness shear in xxviThis statement can, of course, be arrived at by employing the classical assumption of a plate being equivalent to a cylinder of infinite radius. 109 the analysis might be questioned. However, although the exact solution cannot be attained, the only manner by which the mechanical analysis can predict the response for this type of situation requires the through-thickness shear term to be present. It is therefore imperative that provision is made for through-thickness shear effects if cases such as these are to be studied. 7.4 Influence of Stress Discontinuities on other Stress Com- ponents In the previous section the existence of a stress discontinuity in the results in ?6 was shown. Additionally, the behaviour of the shear stress (for which the discontinuity is present) was explained. However, as this predicted response does not agree with the exact solution, it is necessary to fully assess the influence this result has on the other stress components in this area. In terms of the radial stress at the outer wall of the ring stiffened tubes, the result is dictated by the equilibrium requirements. Therefore at all points outside of the band of external pressure, the radial stress must be zero. Within the region of external pressure, the radial stress must be equal in magnitude and opposite in direction to the applied pressure. In Fig. 6 -20 is can be seen that the above situation is not realised fully. While the radial stress is in an average sense equal to the applied external pressure, there is much variation in the magnitude of the stress across this boundary. However, it can be seen that the radial stress exhibits some form of Fourier series type convergence and the discontinuities in radial stress correspond to the points at which the shear stress is discontinuous. Now, since the radial and shear strains are related to one another by the variation in radial displacement (see Eqn. 9), it is unlikely that if the shear stress is not predicted correctly, that the radial stress will be. Therefore the overshoot and oscillation shown for the radial stress at this location is not unexpected. And in fact, similar behaviour will be found in a FEM simulation for similar reasons. It does however introduce problems for data processing. As some of the stress results presented here are calculated at the border of the radial pressure regions, the stress results on the outer walls are not predicted exactly. Thus, as was shown, the radial stress at the outer wall of the tubes may not be zero in some cases. The reason for this is evident from Fig. 6 -21. Similar discontinuities are also evident on the other stress components at the tube outer wall. From the other results presented for stress variation along the tube length, it can be seen that this is, however, a localised effect. 110 7.5 Transient Effects in Composite Tubes From previous research (Kardomateas8) it was expected that the variation in stress for the composite tubes would vary considerably with time. From the results pre- sented in Fig. 6 -28 - 6 -32 indicate that while the stresses do not change by orders of magnitude, there is still a significant variation in the stress magnitude with time. Additionally, it can be seen that the stress states within the CSM tube exhibit higher magnitude at early points in time (50 s) compared to the later times (500 s). This behaviour was also shown in the work of Kardomateas 8 and it indicates that it is important to consider the transient thermal stresses as they may be of higher magnitude than those found at steady-state. This behaviour is most probably due to the fact that the temperature profile through the thicknesses of the tubes varies more severely at early times, thereby creating more of a mismatch between the ex- pansion at the inner and outer walls. However, the radial stress in the FW tube increases steadily with time. From Fig.?s 6 -29 and 6 -30, the singular behaviour in the radial stress is clear. There- fore it can be seen that this analysis is able to provide a reasonable approximation for these free-edge effects with relatively low computational cost. As discussed in the introduction, previous work has shown that the radial stress tends towards infinity at the material interface at a free-edge. This ?infinite? radial stress is a mathematical abstraction, as it is impossible for a real system to attain this position. Clearly, the material will yield and plasticly redistribute the load, or, the material will rupture. However, while the infinite stress cannot be attained, the high radial stresses in the vicinity of the free edge are thought to be the primary cause of delamination failure in laminates. To fully appreciate the importance of this effect it must be recognised that delamina- tion failure is in some regards akin to crack growth. And as is well known, a singular stress state is present at the tip of a crack, and under certain circumstances, the crack growth can be unstable. A similar situation is present in delamination failure. Therefore, if the large radial stress causes delamination in a small region near the free end, this may be a sufficient condition to promote catastrophic failure in the material if the delamination is unstable. From these results it is found that the magnitude of the radial stress away from the ends is small, even in terms of low through-thickness strength. Therefore the merit of this analysis might be questioned. It must be remembered however, that this analysis in no way represents a ?worst-case? scenario for the tubes, and it is possible the other loading conditions can induce through-thickness stresses of higher magnitude. However, this further indicates the importance of the peak in stress at the end of the tube which is predicted to be an order of magnitude higher than the stress at the midplane. 111 Additionally, even though the radial stress component is small, Fig.?s 6 -31 and 6 -32 indicate that the other through-thickness stress are of significant magnitude and also vary considerably through the tube thickness. As the tubes in this analysis are not particularly thick this indicates the importance of considering through-thickness effects. If these stresses are additive with the expected tube loading, the resulting stresses will be higher than originally expected by a significant margin. 7.6 Enforcing Stress Equilibrium Conditions by Means of Constraints The formulation of the Rayleigh-Ritz analysis in this work, allows higher order boundary conditionsxxvii to be satisfied, either by applying constraints in the manner presented in this work, or by modifying the displacement functions such that they satisfy these conditions inherently. The latter approach is frequently utilised in works where orthogonal polynomials are selected as the basis functions. Typically, processes such as these are utilised to satisfy natural boundary conditions however, it is in fact possible to formulate constraints such that they will enforce a user specified stress at a particular location. To some extent this is a heuristic process at it implies some prior knowledge of what the stress tensor should be in order to obtain an accurate solution. However, the underlying philosophy is that if the correct stress distribution is known, then by enforcing it, an exact solution would be obtained. However, from the conclusions of the preceding sections on edge effects and stress singularities, it can be said that this process will not always produce the desired results. Obviously, at the locations where the stress is specified, the analysis will conform to these constraints and the results will reflect the inputs. It was shown though, that at some points within a simulation it is impossible for this analysis to predict the correct results (the particular case in point is at a position of stress discontinuity). Therefore the analysis will not be able to attain what it perceives as the position of minimum potential energyxxviii. Firstly, this will lead to an inaccurate global simulation and secondly, to some point, can even cause the analysis to be numerically unstable. Therefore it must be noted that no attempt has been made in the current formulation to implement a technique to enforce stress states. It is interesting to note that an approach similar to the one described here was followed by Zhu and Lam23 wherein the stress compatibility conditions are satisfied within the analysis (which is also based on the Rayleigh-Ritz technique). It is also interesting that although the reasoning in this work is not explicitly stated, when the xxviiIn this context, higher order refers to the differentiability of the boundary conditions. In this case, C1 and above. xxviiiAs opposed the true position of minimum potential energy, which would require an unsym- metric stress tensor 112 results near a free edge were investigated a relaxed boundary condition was applied to the model. Essentially, rather than enforcing the correct boundary conditions throughout the model (which was the technique followed elsewhere in the work), an average condition was used. This average condition is essentially the boundary conditions written in integral form, similar to that described in ?7.3. 7.7 The Possible Benefits of this Analysis for Isotropic Sys- tems Although this work is aimed primarily at use in the analysis of composite structures, it is important not to lose sight of the fact that this work can also be used for analysis of isotropic structures. It was shown in ?5.6 this model predicted the response of an isotropic pressure vessel with a high degree of accuracy. It must be remembered that there still exists a large demand for analysis of isotropic tubing as it is in frequent use. While the secondary effects in composite vessels are not present, thermal loading can still create stresses within isotropic cylinders 6. And while there are a plethora of analyses for different types of loading on such isotropic vessels, this work is a convenient tool for performing such analyses by possessing the ability to consider a combination of loadings with relative ease. The added advantage is that it is not necessary to implement a large number of analytical equations to analyse the tubes, which requires little work on behalf of the user. 113 8 Conclusions 8.1 Summary It was shown that by allowing 2-dimensional variation for the independent dis- placements and by accounting for through-thickness shear, an analysis method was formulated which allows the analysis of finite length, orthotropic tubes. A simple heat transfer analysis was implemented to describe the transient thermal behaviour of a composite tube. This analysis was based on a explicit, one dimen- sional model. The results presented in this work, indicate that this analysis produces accurate results, whilst it was found that the analysis is also computationally effi- cient. Thus this analysis is able to produce accurate representations of the thermal state of a composite tube at any point in time. The heat transfer analysis was validated against existing results based on an ana- lytical analysis8. The results between these analyses were shown to agree closely. The mechanical analysis was validated in a number of cases against existing results, including those obtained from a commercial finite element software package. It was shown that the results obtained from the mechanical analysis showed excellent agree- ment with the existing results. The only case where the results obtained from this analysis deviated from existing results is for the results presented by Kardomateas8. However, in this case the results from this analysis coincided with results from a FEM analysis. It was therefore concluded that there was some error in the original presentation of the results by Kardomateas 8. The efficiency of the mechanical analysis in this work was compared to a commercial software package and was found to possess a far greater efficiency than the FEM solution. Thus, for similar accuracy, the solution presented in this work would possess far fewer degree?s of freedom compared to the FEM solution. It was shown that while this work was not intended to analyse such effects, the mechanical analysis presented in this work possessed the ability to obtain fairly ac- curate results in the presence of free-edge effects. Whilst it is effectively impossible for a numerical analysis to describe this phenomenon exactly, it was shown that re- sults obtained from this analysis correlated well with existing results (which included results obtained from a FEM analysis, the typical tool of choice in the investigation in free-edge effects). It was however noted that if the edge-effects are of primary im- portance, the FEM solution is more applicable to solving these problems, primarily as a result of the formulation and implementation differences between the solutions. However, if the stresses are not the primary focus of the analysis, the solution offered in this work can provide good results at comparatively low computational cost. Novel results for two ring-stiffened, finite length, composite cylinders were presented. From an analysis viewpoint, this problem demonstrated the ability of the present analysis to model finite length, composite cylinders under a combination of loading. 114 The effect of the ring stiffeners was incorporated into the mechanical analysis by a applying a discrete band on external pressure to the outer surface of the tubes. The heat transfer analysis was utilised to determine the thermal response of the two geometrically similar tubes which possess different material make-up. These thermal loads were combined with internal and external pressure to fully analyse the tubes. From these results, the necessity to incorporate through-thickness shear effects was demonstrated again. It was found that the stresses varied significantly in magnitude in the vicinity of the ring stiffener. In both tubes, stresses of significant magnitude with regard to material failure were obtained. It was also shown that, although the tubes are not particularly thick, the stress variation through the thickness of the tubes is large. This enforces the importance of considering through-thickness effects. Additional results were presented which investigated the transient behaviour of the above two tubes. In this simulation, the tubes were modelled as finite length tubes with free ends under the action only of a thermal contraction. It was shown that the stresses obtained in both tubes varied by significant amounts with time. It was also found that the stresses obtained at early times within the simulation could exhibit higher magnitudes than those obtained a later times. Again it was found that a large variation in stress was present through the thickness of both tubes. These results indicate the importance of considering not only the through-thickness effects, but the transient effects as well. As the tubes are constructed from composite materials, a free-edge effect was ex- pected. This effect was illustrated in the behaviour of the radial stress along the length in both tubes at the outer material interface. The tendency for the radial stress to become infinite at the free edge was clearly illustrated. This behaviour is in line with that presented in existing results and works. This infinite radial stress is thought to be the primary cause of delamination failure in composite laminates. These results again indicate that the analyses developed in this work are able to provide reasonable approximations for these effects with low computational cost. 115 8.2 Conclusions The conclusions which can be drawn from this work are listed below. 1. It is possible to utilise the thermal analysis to predict the temperature response of a composite tube. 2. The mechanical model allows the stress response of a composite tube to be predicted. 3. Provision was made for a combination of loading including pressure effects and thermal expansion. 4. It was shown that the model adequately satisfied the equilibrium conditions for a homogeneous, finite length orthotropic tube with free ends. 5. In the case of non-homogeneous tubes it was verified that the mechanical model predicted singular stress states at the free edge. This is a well known phenomenon, termed the free-edge effect. 6. Although it was not the intent of this work to accurately predict the edge effects, it was shown that the model utilised was able to predict them with at least as good an accuracy as other presented techniques (whose primary purpose was to investigate the edge effects). 7. It was also found that in the case of a discreet pressure load, the mechanical model predicted a stress discontinuity for the shear stress. This behaviour is in line with that of FEM for similar loadings. Due to the underlying assumptions in the formulation of the mechanical model (primarily the assumption of a symmetric stress tensor) this type of behaviour cannot be avoided. 8. In line with similar work it was shown that although there was a stress dis- continuity, the integral over the discrete pressure band yielded the correct analytical result of zero. 9. The effects of this stress discontinuity were shown to affect the other compo- nents of stress. The influence of these discontinuities was shown to only act over a local area though. 10. The models presented allow the modelling of composite tubes of finite length under the action of any combination of the loads: - Pressure Loading - Transient Thermal Loading - Shrinkages 11. The transient behaviour of two composite tubes under thermal loading was investigated. It was shown that transient effects are important as the stress states varied considerably in magnitude with time in both tubes. It was also found that the stresses at early times can possess higher magnitudes than those obtained at later times. 116 12. It was shown that even for small radial stress states arising from thermal loading, the other stress components could be of much higher magnitude, to the point where these stresses could contribute significantly to the failure of a vessel. 8.3 Recommendations At this point in the development of this work it is possible to investigate a number of different avenues of related research. A small number of such areas are proposed in this section. With regard to extending the analytical analysis presented in this work, the imme- diate area of interest is in the dependence of the material properties on the thermal states of the vessel. Throughout this work it has been assumed that the material properties (both mechanical and thermal) are temperature independent. In reality this is not strictly true, particularly in the cure and post-cure processes. There exists therefore an opportunity to investigate the temperature variation of these proper- ties and incorporate them into the model. The resulting model will then become a nonlinear analysis by virtue of this dependence. Another extension which could be performed is to investigate other geometrical shapes than cylinders. This requires that the problem be reformulated (most prob- ably in another coordinate frame which better describes the geometry). In typical pressure vessels, the vessels consist of a roughly cylindrical section which is capped by two ends which are either hemispherical or elliptical. By considering these shapes the analysis would then be able to predict the response of such vessels. This work as well as the proposed research extensions are all analytic approaches; however, it is still necessary to validate these results with the behaviour experienced in reality. This is not a trivial task as the nature of composite tubes presents diffi- culty in extracting strains throughout the thickness of the tube. Hence, it becomes difficult to determine the through-thickness stress state. However, there are a num- ber of techniques which can be employed to perform this analysis. These results can then be utilised to correlate results obtained from the present work with those found in reality. 117 References 1. A. Dasgupta and K.H. Huang. A layer-wise analysis for free vibrations of thick composite spherical panels. Journal of Composite Materials, 31(7):658?671, 1997. 2. L.K. Jain and Y. Mai. Stresses and deformations induced during manufactur- ing. part I: Theoretical analysis of composite cylinders and shells. Journal of Composite Materials, 31(7):672?695, 1997. 3. R. Tucker, P. Compston, and P.-Y.B. Jar. The effect of post-cure duration on the mode I interlaminar fracture toughness of glass-fibre reinforced vinylester. Composites: Part A, 32:129?134, 2001. 4. C.T. Herakovich. Mechanics of Fibrous Composites. John Wiley & Sons, Inc., first edition, 1998. 5. A.E.H. Love. A Treatise on the Mathematical Theory of Elasticity. Dover Pub- lications, fourth edition, 1944. 6. W.C. Young and R.G. Budynas. Roark?s Formulas for Stress and Strain. McGraw-Hill, seventh edition, 2002. 7. M.A. Stone, I.F. Schwartz, and H.D. Chandler. Residual stresses associated with post-cure shrinkage in GRP tubes. Composites Science and Technology, 57:47?54, 1997. 8. G.A. Kardomateas. Transient thermal stresses in cylindrically orthotropic tubes. Journal of Applied Mechanics, 56:411?417, 1989. 9. H.W. Wiersma, L.J.B. Peeters, and R. Akkerman. Prediction of springforward in continuous-fibre/polymer L-shaped parts. Composites Part A, 29A:1333?1342, 1998. 10. C.C. Chamis. Lamination residual stresses in multilayered fiber composites. NASA TN D-6146, 1971. 11. C. Mittelstedt and W. Becker. Interlaminar stress concentrations in layered structures: Part I - a selective literature survey on the free-edge effect since 1967. Journal of Composite Materials, 38(12):1037?1062, 2004. 12. I.S. Raju, J.D. Whitcomb, and J.G. Goree. A new look at numerical analyses of free-edge stresses in composite laminates. NASA TP 1751, 1980. 13. J-H Byun, C-H Lee, S-W Song, S-K Lee, B-S Kim, and C.R. Joe. Stitching effect on flexural and interlaminar properties of mwk textile composites. In ICCM 15 - Fifteeneth Conference on Composite Materials. ICCM. 14. J. Zhang and Y. Wei. A hierachical multiscale model for in-plane mechanical strength of stitched composite laminates. In ICCM 15 - Fifteeneth Conference on Composite Materials. ICCM. 118 15. M.A. Stone, I.F. Schwartz, and H.D. Chandler. Residual stresses arising from partial in-service post-cure of GRP structures. Composites, 25(3):177 ? 181, 1994. 16. S. Timoshenko and J.N. Goodier. Theory of Elasticity. McGraw-Hill Book Company, Inc., second edition, 1951. 17. A.P. Misovec and J. Kempner. Approximate elasticity solution for orthotropic cylinder under hydrostatic pressure and band loads. ASME Papers, 37(1):101? 108, 1970. 18. L.P. Kolla?r, J.M. Patterson, and G.S. Springer. Composite cylinders subjected to hygrothermal and mechanical loads. International Journal of Solids and Structures, 29(12):1519?1534, 1992. 19. Collected papers on instability of shell structures. NASA TN D-1510, 1962. 20. N.F. Knight Jr. and M.P. Nemeth. Stability analysis of plates and shells: A collection of papers in honor of Dr. Manuel Stein. NASA/CP-1998-206280, 1998. 21. K. Chandrashekhara and P. Gopalakrishnan. Elasticity solution for a multilay- ered transversely isotropic circular cylindrical shell. Journal of Applied Mechan- ics, 49:108?114, 1982. 22. K.H. Ip, W.K. Chan, P.C. Tse, and T.C. Lai. Vibration analysis of orthotropic thin cylindrical shells with free ends by the Rayleigh-Ritz method. Journal of Sound and Vibration, 195(1):117?135, 1996. 23. C. Zhu and Y.C. Lam. A Rayleigh-Ritz solution for local stresses in composite laminates. Composites Science and Technology, 58:447?461, 1998. 24. J.N. Reddy and C.M. Wang. An overview of the relationships between solutions of the classical and shear deformation plate theories. Composites Science and Technology, 60:2327?2335, 2000. 25. A.K. Noor and W.S. Burton. Assessment of computational models for multilay- ered composite shells. Applied Mechanics Reviews, 43(4):67?97, 1990. 26. T. Kant and K. Swaminathan. Estimation of transverse/interlaminar stresses in laminated composites - a selective review and survey of current developments. Composite Structures, 49:65?75, 2000. 27. G.P. Dube, S. Joshi, and P.C. Dumir. Nonlinear analysis of thick shallow spheri- cal and conical orthotropic caps using Galerkin?s method. Applied Mathematical Modelling, 25:755?773, 2001. 28. A.J.M. Ferreira, C.M.C. Roque, and P.A.L.S. Martins. Radial basis functions and higher-order shear deformation theories in the analysis of laminated com- posite beams and plates. Composite Structures, 66:287?293, 2004. 29. M. Ganapathi, B.P. Patel, and D.P. Makhecha. Nonlinear dynamic analysis of thick composite/sandwich laminates using an accurate higher-order theory. Composites Part B: Engineering, 35:345?355, 2004. 119 30. L.P. Kolla?r and G.S. Springer. Stress analysis of anisotropic laminated cylinders and cylindrical segments. International Journal of Solids and Structures, 29(12): 1499?1517, 1992. 31. M.A. Stone and H.D. Chandler. Errors in double sine series solutions for simply supported symmetrically laminated plates. International Journal of Mechanical Science, 38(5):517?526, 1996. 32. C.C. Chao, T.P. Tung, and Y.C. Chern. Buckling of thick orthotropic spherical shells. Composite Structures, 9(2):113 ? 137, 1988. 33. R.D. Cook, D.S. Malkus, M.E. Plesha, and R.J. Witt. Concepts and Applications of Finite Element Analysis. John Wiley & Sons, Inc., fourth edition, 2001. 34. B. Baharlou and A.W. Leissa. Vibration and buckling of generally laminated composite plates with arbitrary edge conditions. International Journal of Me- chanical Science, 29:545?555, 1987. 35. J. Schapery. Thermal expansion coefficients of composite materials based on energy principles. Journal of Composite Materials, 2:380?404, 1968. 36. M.R. Spiegel. Schaums Mathematical Handbook of Formulas and Tables. McGraw-Hill, 2nd edition, 1998. 37. A.F. Mills. Heat Transfer. Prentice Hall, 1999. 38. MSC.Nastran 2005 Quick Reference Guide. 2005. 39. J.Q. Ye and H.Y. Sheng. Free-edge effect in cross-ply laminated hollow cylinders subjected to axisymmetric transverse loads. International Journal of Mechanical Sciences, 45:1309?1326, 2003. 40. K. Bhaskar, T.K. Varadan, and C. Jacob. Free-edge stresses in laminated cylin- drical shells due to axisymmetric transverse loads. Journal of the Aeronautical Society of India, 52:26?38, 2000. 41. G. Mishuris and A. ?Ochsner. Edge effects connected with thin interphases in composite materials. Composite Structures, 68:409?417, 2005. 42. C. Mittelstedt and W. Becker. Interlaminar stress concentrations in layered structures: Part II - closed-form analysis of stresses at laminated rectangular wedges with arbitrary non-orthotropic layup. Journal of Composite Materials, 38(12):1063?1090, 2004. 43. S.P. Walker. Thermal Effects of the Bearing Behaviour of Composite Joints. PhD thesis, University of Virginia. 44. B. Mutnuri. Thermal conductivity characterization of composite materials. Mas- ter?s thesis, West Virginia University, 2006. 45. MatWeb, The Online Material Database. http://www.matweb.com, July 2006. 46. DoW Plastics. Derakane Epoxy Vinyl Ester Resins - Chemical Resistance and Engineering Guide. 1995. 120 47. B. Mutnuri. Thermal conductivity characterization of composite materials. Mas- ter?s thesis, West Virginia University, 2006. 48. C.C. Chamis. Simplified composite micromechanics equations for hygral, ther- mal and mechanical properties. NASA TM 83320, 1983. 121 A Detail Derivation of Mechanical Model A detailed derivation of the some of the mathematical relations are provided in this appendix. The philosophy is that this section will serve as a reference to the development of the analysis allowing it to be re-created, rather than presenting the analysis in its entirety. A.1 Formulation of Strain Relationships Following the analysis in ?3 the strain-displacement can be written as: ? = B ?D (101) Performing the differentiation and expanding the above product yields: ? = ? ??? ?z ?? ?r ?rz ? ??? = ? ???????????????? N? n=0 M? m=0 WmnmZm?1Rn 1 L P? p=0 Q? q=0 WpqRpr?1Zq P? p=0 Q? q=0 WpqpRp?1 1 ti Zq N? n=0 M? m=0 WmnnZmRn?1 1 ti + P? p=0 Q? q=0 WpqqZq?1Rp 1 L ? ???????????????? (102) Eqn. 102 can be expanded into: ? = ? ????????? 0 ? ? ? MZM?1RN 1 L 0 ? ? ? 0 0 ? ? ? 0 1 r ? ? ? RP r?1ZQ 0 ? ? ? 0 0 ? ? ? PRP?1 1 ti ZQ 0 ? ? ? NZMRN?1 1 ti 0 ? ? ? QZQ?1RP 1 L ? ????????? ? ???????? Wmn|m=n=0 . . . WMN Wpq|p=q=0 . . . WPQ ? ???????? (103) 122 A.2 Formulation of Energy Potential Relationship The potential energy for the tube system may be represented as: ? = 1 2 ? ? V ? ?TC? dV (104) In this case the strain matrix can be written as the actual strains less the shrinkage values. Hence, ? = ? ??? ?z ? sz ?? ? s? ?r ? sr ?rz ? ??? = ?1 ? ?2 (105) Following the elementary relationships of linear algebra, the following holds : ?T = ?T1 ? ? T 2 (106) The matrix product can then be expanded as: ?TC? = ( ?T1 ? ? T 2 ) C (?1 ? ?2) = ?T1 C?1 ? ?T2 C?2 ? ?T2 C?1 ? ?T1 C?2 (107) Since the material constitutive matrix, C, is symmetric, by virtue of matrix prop- erties, ?T2 C?1 = ?T1 C?2. Therefore Eqn. 107 can be rewritten as: ?TC? = ?T1 C?1 ? ?T2 C?2 ? 2?T1 C?2 (108) The stationary potential energy is arrived at by the following relationship: ?? ?W = 0 = 1 2 ? ? V ? ? ?W ( ?TC? ) dV (109) Since the term, ?2, is only dependent on the material free shrinkages and not the displacement coefficients, W , it will fall away in the above differentiation. Therefore with some manipulation, Eqn. 109 can be written as: ? ? V ? ? ?W ( ?T1 C?1 ) dV = ? ? V ? ? ?W ( ?T2 C?2 ) dV (110) For sake of convenience, Eqn. 110 will be written as: ?U ?W = ?V ?W (111) 123 Considering the left hand side of Eqn. 110, by the chain rule: ?U ?W = ? ? V ? ? ?W ( ?T1 ) C?1 + ?T1 C ? ?W (?1) dV (112) By implementing the strain-displacement relationship, the above equation can be rewritten as: ?1 = B ?D ? ?W ( ?T1 ) C?1 + ?T1 C ? ?W (?1) = ? ?W ( W TBT ) CBW +W TBTC ??W (BW ) = BTCBW +W TBCBT = 2 ( BTCB ) W (113) Therefore it can be said: ?U ?W = 2 ? ? V ? BTCB dV.W = 2KW (114) ?V ?W = 2 ? ? V ? ? ?W ( W TBTC?2 ) dV = 2 ? ? V ? BTC?2 dV = 2L (115) Which yields the familiar matrix form of : KW = L. A.3 Formulation of the Loading matrix The loading matrix is given by: L = ? ? V ? BTC?2 dV (116) The product, C?2 can be written as: C?2 = ? ??? C11 C12 C13 0 C21 C22 C23 0 C31 C32 C33 0 0 0 0 C66 ? ??? ? ??? sz s? sr 0 ? ??? ? ??? C11sz + C12s? + C13sr C21sz + C22s? + C23sr C31sz + C32s? + C33sr 0 ? ??? = ? ??? A0 A1 A2 0 ? ??? (117) 124 Hence, BTC?2 is given as: [ BTz 0 0 BTrz1 0 BT? BTr BTrz2 ] ? ??? A0 A1 A2 0 ? ??? = [ BTz A0 BT? A1 +BTr A2 ] (118) In expanded notation, the above matrix is written as: ? ? [ BTz A0 BT? A1 +BTr A2 ] dr dz = ? 1 0 ? 1 0 ? ????????????? 0 . . . RnmZm?1 1 L A0 ZqA1 1 r + 0 . . . Rp r ZqA1 + pRp?1Zq 1 ti A2 ? ????????????? rtiLdRdZ = ? ????????????? 0 . . . 1 n+ 2 A0t2i + ati n+ 1 A0 A1Lti q + 1 . . . A1Lti (p+ 1)(q + 1) + pLtiA2 (p+ 1)(q + 1) + A2La (q + 1) ? ????????????? (119) 125 A.4 Formulation of Stiffness Matrix The Stiffness Matrix is given as: K = ? ? V ? BTCB dV (120) Considering the product: CB = ? ??? C11 C12 C13 0 C21 C22 C23 0 C31 C32 C33 0 0 0 0 C66 ? ??? ? ??? Bz 0 0 B? 0 Br Brz1 Brz2 ? ??? = ? ??? C11Bz C12B? + C13Br C21Bz C22B? + C23Br C31Bz C32B? + C33Br C66Brz1 C66Brz2 ? ??? (121) The full product can then be written as: BTCB = [ BTz 0 0 BTrz1 0 BT? BTr BTrz2 ] ? ??? C11Bz C12B? + C13Br C21Bz C22B? + C23Br C31Bz C32B? + C33Br C66Brz1 C66Brz2 ? ??? = [ BTz C11Bz +BTrz1C66Brz2 BTz (C12B? + C13Br)+ BT? C21Bz +BTr C31Bz +BTrz2C66Brz1 BT? (C22B? + C23Br)+ BTrz1C66Brz2 BTr (C32B? + C33Br) + BTrz2C66Brz2 ] (122) For elegance the above matrix will be written as: K = [ K1 +K2 K3 +K4 +K5 K3 +K4 +K5 K6 +K7 +K8 +K9 ] (123) As can be seen from Eqn. 123, the symmetry of the stiffness matrix is implied. The integrations of the individual stiffness matrices in Eqn. 123 is performed in the next section. 126 A.4.1 Integrations of Stiffness Matrix As can be seen from the form of the individual matrices required that the indices of the exponents of the non-dimensionalised displacements use separate counters. i.e. the indices m,n become j, k and the indices p, q become u, v. K1 = 1? 0 1? 0 ? ??? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? Rk+njmZm+j?2 1 L2 ? ???C11Ltir dZ dR = ? ???? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? ( t2i k + n+ 2 + tia k + n+ 1 ) mj m+ j ? 1 C11 L ? ???? (124) K2 = 1? 0 1? 0 ? ?? 0 ? ? ? 0... . . . ... 0 ? ? ? knRk+n?2t?2i Zm+j ? ??C66Ltir dZ dR = ? ???? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? knL ti(m+ j + 1) ( ti k + n + a k + n? 1 ) C66 ? ???? (125) K3 = 1? 0 1? 0 C12 ? ??? 0 ? ? ? . . . . . . . . . ? ? ? Rk+pjzq+j?1 1 Lr ? ???Ltir dZ dR = ? ??? 0 ? ? ? . . . . . . . . . ? ? ? C12 jti (q + j)(p+ k + 1) ? ??? (126) K4 = 1? 0 1? 0 C13 ? ??? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? Rk+p?1pjZq+j?1 1 Lti ? ???Ltir dZ dR = ? ???? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C13 ( pjti (q + j)(p+ k + 1) + apj (q + j)(k + p) ) ? ???? (127) 127 K5 = 1? 0 1? 0 C66 ? ??? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? nq tiL Rn+p?1zm+q?1 ? ???Ltir dZ dR = ? ???? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C66nq ( ti (n+ p? 1)(m+ q) + a (n+ p)(m+ q) ) ? ???? (128) K6 = 1? 0 1? 0 ? ???? C22 r2 ? ? ? C22RpZq r2 . . . . . . . . . C22RuZv r2 ? ? ? C22 1 r2 Rp+uZq+v ? ????Ltir dZ dR = ? ????? C22Lti ? 1 r dR ? ? ? C22Lti 1 q + 1 ? Rp r dR . . . . . . . . . C22Lti v + 1 ? Ru r dR ? ? ? C22Lti q + v + 1 ? Rp+u r dR ? ????? (129) K7 = 1? 0 1? 0 ? ??? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C33puZv+q 1 t2i Ru+p?2 ? ??? rLti dZ dR = ? ???? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C33puL ti(v + q + 1) ( ti u+ p + a u+ p? 1 ) ? ???? (130) 128 K8 = 1? 0 1? 0 ? ??? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C66 vq L2 Ru+pZv+q?2 ? ???Ltir dZ dR = ? ???? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C66 vqti L(v + q ? 1) ( t u+ p+ 2 + a u+ p+ 1 ) ? ???? (131) K9 = 1? 0 1? 0 ? ??? 0 ? ? ? 0 . . . . . . . . . 0 ? ? ? C32pRu+p?1Zq+v u+ p tir ? ???Ltir dZ dR = ? ??? 0 ? ? ? . . . . . . . . . ? ? ? C32 L q + v + 1 ? ??? (132) 129 B Estimation of Secondary Poisson?s Ratio, ?23 The procedure for estimating the secondary Poisson?s ratio is presented in this sec- tion. The estimation begins by considering an applied transverse strain ?2. This situation is sketched in Fig. B -1. Figure B -1: Cross Section of a Composite Laminate As in the micromechanics approach, it is assumed that the transverse stress in the fibre and matrix are equalxxix, therefore: ?f = ?m = ?2 (133) By virtue of Hooke?s Law the transverse strain can be written for the matrix and the fibre as: ?f = ?2 Ef2 , ?m = ?2 Em (134) Note that provision is made in the above relationship for anisotropic fibres by virtue of the inclusion of the stiffness term Ef2 . Possibly the most common example of such fibres would be carbon fibres. Considering the global laminate strain the following relationships can be determined: ?2 = ?2 E2 (135) xxixmuch like springs in series 130 ?2 = ?2 ( Vf Ef2 + Vm Em ) (136) Where Vf and Vm are the volume fractions of the fibre and matrix respectively. In the same manner as the transverse strain, the through thickness strain can be written as: ?3 = (?3)fVf + (?3)mVm = (?2)f (?23)fVf + (?2)m?mVm (137) Incorporating the above relationships, the Poisson?s ratio ?23 can be written as: ?23 = ? ?3 ?2 = ? (?2)f (?23)fVf + (?2)m?mVm( Vf Ef2 + VmEm ) (138) In the fibre the following relationships apply: ?2f = ?2 Ef2 (139) ?3f = ?2 Ef2 ?23 (140) Therefore: ?23 = ? Vf?23f Ef2 + Vm?m Em Vf Ef2 + Vm Em (141) As can be seen, the only unknown in Eqn. 141 is the fibre Poisson?s ratio ?23f . As glass fibres are isotropic, the property, ?23f can simply be calculated according to the relationship: ?23f = ? = E 2G ? 1 (142) 131 C Sample Calc for Woven Roving A sample calculation to illustrate the methods utilised for estimating the composite material properties is presented in this section. The calculation is performed for a woven roving laminate, however, the technique is also applicable for other laminates. To begin, the properties of the constituent materials is required. The properties used in this case are presented in Table C -1. The properties for the glass were obtained from MatWeb45 for generic E-glass fibre. The Young?s modulus and Pois- son?s ratio for the resin were obtained from Stone et al. 7. The resin specific heat capacity, density and CTE were obtained from manufacturer?s data46. Data for the conduction coefficient for this resin could not be obtained, therefore the value for a similar vinyl-ester resin was utilised 47. The shear moduli for both materials was calculated from standard relationships for isotropic materials. In all cases, if better data is available (i.e. test data) it should be utilised. Table C -1: Laminate Constituent Material Properties Property Resin Fibre E (GPa) 3.5 72.4 ? 0.35 0.20 G (GPa) 1.3 30.2 ? (?C?1) 65? 10?6 5.0? 10?6 k (W/m ?K) 0.178 1.3 c (J.kg?1) 1400 810 ? (kg.m?3) 1126 2540 In order to predict the final laminate properties it is first necessary to determine the properties of a uni-directional laminate. This is accomplished in the following section. C.1 Uni-directional Properties for the Woven Roving Lam- inate For the woven roving laminate a fibre volume fraction of 31% which corresponds to a mass fraction of approximately 50%, which is a commonly attained value in practise. The in-plane properties can be calculated according to standard approaches4: E1 = Ef Vf + Em (1? Vf ) E1 = 72400? 0.31 + 3500? (1? 0.31) = 24859MPa (143) 132 E2 = Ef Em EmVf + Ef (1? Vf ) E2 = 72400? 3500 3500? 0.31 + 72400? (1? 0.31) = 4965MPa (144) G12 = ( Vf Gf + Vm Gm ) ?1 G12 = ( 0.31 30.2 + 1? 0.31 1.3 ) ?1 = 1843MPa (145) ?12 = ?fVf + ?m (1? Vf ) ?12 = 0.2? 0.31 + 0.35? (1? 0.31) = 0.304 (146) The Poisson?s ratio ?23 can be calculated according to the method presented in App. B: ?23 = Vf?23f Ef2 + Vm?m Em Vf Ef2 + Vm Em ?23 = 0.31? 0.2 72400 + (1? 0.31)? 0.35 3500 0.31 72400 + 1? 0.31 3500 = 0.347 (147) The CTE?s can be calculated according to the method of Schapery35: ?1 = EfVf?f + EmVm?m EfVf + EmVm ?1 = 72400? 0.31? 5? 10?6 + 3500? 0.69? 65? 10?6 72400? 0.31 + 3500? 0.69 = 10.83? 10 ?6 (148) ?2 = ?fVf + ?mVm + (?m ? ?f ) (Ef?m ? Em?f ) VfVmEfVf + EmVm ?2 ? 106 = 5? 0.31 + 65? 0.69+ (65? 5)(72400? 0.35? 3500? 0.2) 0.31? 0.6972400? 0.31 + 3500? 0.69 ?2 = 59.12? 10?6 (149) With the above properties it is possible to determine the laminate properties as shown in the following section. 133 C.2 Woven Roving Laminate Properties As discussed earlier, the laminate through-thickness stiffness is assumed to be equal to the transverse stiffness of the uni-directional laminate utilised in the preceding section. Therefore: Er = E3 = E2 = 4965MPa (150) Similarly, the through-thickness shear modulus, Grz, is equal to the in-plane shear modulus of the uni-directional laminate: Grz = G12 = 1843MPa (151) The in-plane laminate properties are calculated by the transformed stiffness ap- proach described previously. However, for this laminate, as the layers of fabric are orthogonal, the calculations can be greatly simplified. Thus the in-plane properties can be calculated as: E? = Ez = A11 ? A212/A22 ttot (152) ?z? = A12 A22 (153) Where: A = N? k=1 ?Qktk (154) Where ?Q is the stiffness matrix for each layer in the laminate. As stated previously, the orthogonal nature of this laminate allows the expression to be simplified to: A11 = 0.5 (Q11 +Q22) (155) A22 = 0.5 (Q11 +Q22) (156) A12 = 0.5 (Q12 +Q21) (157) 134 and Q = S?1 (158) Where S is the compliance matrix given by: S = ? ?????? 1 E1 ??12 E1 0 ??12 E1 1 E2 0 0 0 1G12 ? ?????? (159) Therefore: A11 = 15191MPa A12 = 1535MPa A22 = 15191MPa Ez = E? = 15036MPa (160) ?z? = 0.101 (161) The CTE?s can be calculated as: ? = A?1 N? k=1 ?Qk?ktk (162) Which yields: ?z = ?? = 20.4? 10?6 /?C (163) The through-thickness CTE can be calculated according to the procedure presented in Herakovich4: ?r = 1 ttot K? k=1 ( ?Sk13?kz + ?Sk23?k? + ?Sk36? kz? + ?kr ) tk (164) 135 Where: ?kz = ?Qk11(?z ? ?kz) + ?Qk12(?? ? ?k?) (165) ?k? = ?Qk12(?z ? ?kz) + ?Qk22(?? ? ?k?) (166) ? kz? = ?Qk16(?z ? ?kz) + ?Qk26(?? ? ?k?) (167) Performing the necessary substitutions: ?r = 69.6? 10?6/?C (168) It is worth noting that the through-thickness CTE predicted above is larger than that of the resin alone. This behaviour is due to the additional in-plane effects which are accounted for in the above technique. The through-thickness Poisson?s ratios can be calculated in a similar manner to the through-thickness CTE4: ?zr = ? ( A?111 F1 + A?112 F2 + A?116 F6 ) ttotA?111 (169) Where the elements A?1ij are the elements in the matrix A?1. The elements Fi can be determined as: Fi = K? k=1 [( ?Sk13 ?Qk1i + ?Sk23 ?Qk2i + ?Sk36 ?Qk6i)tk] (170) ?A?111 = 6.65? 10?5 (171) ?A?112 = ?6.70? 10?6 (172) ?A?116 = 0 (173) performing the substitutions yields: ?zr = 0.355 (174) A similar procedure, also presented in Herakovich 4, can be followed for the Poisson?s ratio, ??r. In this case, as the material is orthotropic, the ratio is exactly equal to ?zr. 136 The laminate thermal properties can be calculated according to the relationships presented in Chamis48: ? = ?f Vf + ?m (1? Vf ) ? = 2540? 0.31 + 1126? (1? 0.31) = 1564 kg.m?3 (175) c = ?f Vf cf + ?m (1? Vf ) cm ? c = 2540? 0.31? 810 + 1126? (1? 0.31)? 1400 1564 = 1103 J/kg.K (176) k = ( 1? ? Vf ) km + km ? Vf 1? ? Vf (1? km/kf ) k = (1??0.31) 0.178 + 0.178 ? 0.31 1? ? 0.31(1? 0.178/1.3) = 0.2697W/m.K (177) 137