LIBRARY Michigan State University This is to certify that the dissertation entitled CONTACT MECHANICS OF LAYERED COMPOSITES UNDER AXISYMMETRIC INDENTATION presented by ZHENWEN WANG has been accepted towards fulfillment of the requirements for the PhD. degree in ENGINEERING MECHANICS «9%»; 0%; Major Professor’s Signature [Olin 771 22,05.— U . Date MSU is an Affirmative Action/Equal Opportunity Institution fi~—-.— —-—. - WW- .r fi" PLACE IN RETURN Box to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 2/05 m/CIRWDI—nhdd-DJS CONTACT MECHANICS OF LAYERED COMPOSITES UNDER AXISYMMETRIC INDENTATION By Zhenwen Wang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSIPHY Department of Mechanical Engineering 2005 ABSTRACT CONTACT MECHANICS OF LAYERED COMPOSITES UNDER AXISYMMETRIC INDENTATION By Zhenwen Wang The analytical studies on contact mechanics have been limited to problems of either half space or single layer due to mathematical difficulty, though layered composites subjected to indentation are commonly encountered in industrial applications. The studies of indentation of layered composites, however, are based on numerical approaches. This dissertation provides a theoretical method for the contact mechanics of layered composites subjected to axially symmetric indentation. A new function is introduced in this dissertation to reduce the complexity of the mathematical process, and mathematical solutions are provided for all the problems investigated in this research. However, the mathematical solution for the final integration could only be obtained for the point loading condition. A numerical method was used to evaluate the final results for other loading conditions, such as uniform stress, flat indentation and spherical indentation. In this dissertation, the effects of material property, layer thickness, boundary condition, loading condition, and lamination on contact mechanics were investigated for the cases of a half space, a single layer bonded to a rigid base, a single layer bonded to an elastic half space and two-layered composites bonded to a rigid base. The dissertation also investigated the frictional effect at the contact interface. Both shear slip and normal separation theories were incorporated into the mathematical formulation, allowing the study of debonding at the interface between layers. New transformed shear slip and normal separation coefficients are proposed to study the imperfect bonding interfaces with a finite length. Contact mechanics models have been proposed based on numerical results. These models provide insight into the relationships among total load, maximum displacement, contact radius, layer thickness and material properties, and guidelines for engineering applications. ACKNOWLEDGEMENTS The author is indebted to many people during the research and development of this dissertation. First of all, I own many thanks to my adviser, Dr. Dahsin Liu, for his encouragement through the years to keep from giving up. Without Dr. Liu's tireless work with me over the weekends to accommodate my regular working schedule, I would never have been able to turn in this dissertation. Furthermore, the suggestions and comments from the professors in my guidance committee are invaluable and much appreciated. I also own a debt to my wife, Xinlan, who encouraged me to come back to the program after a year-long interruption, for her patience, understanding and continuous support through the program to achieve my life time goal. Finally, thanks goes to my two lovely children, six year old Claire and three old Collin, who were understanding when father had to go to work on the weekends in the past few years. TABLE OF CONTENTS List of Tables ....................................................................................................... vii List of Figures ..................................................................................................... viii List of Figures ..................................................................................................... viii Chapter One Introduction ................................................................................... 1 1.1 Introduction ..................................................................................................... 1 1.1.1 History ......................................................................................................................................... l 1.1.2 Summary of contact mechanics research .................................................................. 1 1.1.3 Contact mechanics in composite material ................................................................. 6 1.1.4 Approach to layered composites .................................................................................... 8 1.2 Object of this study ......................................................................................... 9 1.3 State of the problem ..................................................................................... 11 1.4 Contents of this dissertation ......................................................................... 11 Chapter Two Development of Theoretical Analysis ............................................ 15 2.1 Elasticity Analysis Based on Cylindrical Coordinate System ........................ 15 2.2 Hankel Transformation and Bessel Functions .............................................. 20 2.3 Displacement Solutions ................................................................................ 22 2.4 Stress Solutions ............................................................................................ 23 2.5 Exact Solution for Half Space under Point Load ........................................... 29 2.6 Traction Boundary Conditions for Contact Problems .................................... 35 Chapter Three Numerical Solution and Verifications ....................................... 38 3.1 Elastic Layer on Rigid Base .......................................................................... 38 3.2 Elastic Layer on Elastic Base ....................................................................... 42 3.3 Two-Layer Composite on Rigid Base ........................................................... 46 3.4 Validations .................................................................................................... 50 3.4.1 Numerical integration vs. theoretical solution ........................................................ 51 3.4.2 Validations among cases ................................................................................................. 53 (i) Comparisons between case (1) and case (2) .......................................... 54 (ii) Comparisons between case (1) and case (3) ......................................... 57 (iii) Comparisons between case (2) and case (4) ........................................ 59 Chapter Four Basic Applications ..................................................................... 63 . 4.1 Effect of Material Property ............................................................................ 63 4.2 Effect of Thickness ....................................................................................... 67 4.3 Effect of Base Condition ............................................................................... 76 4.4 Effect of Loading Condition ........................................................................... 84 4.5 Effect of Lamination ...................................................................................... 93 4.6 Cylinder Shape lndenter ............................................................................... 98 Chapter Five Contact Theories ....................................................................... 102 5.1 Half space ................................................................................................... 102 5.2 Single layer on rigid base ........................................................................... 103 5.3 Single layer on elastic base ........................................................................ 106 5.4 Two layer composite on rigid base ............................................................. 110 Chapter Six Advanced Study ............................................................................ 114 6.1 Shear Loading ............................................................................................ 114 6.2 Contact Problems with Frictions ................................................................. 118 6.3 Bonding Theory .......................................................................................... 122 6.3.1 Entire Imperfect Interface .............................................................................................. 124 6.3.1.1 Verifications ................................................................................... 126 6.3.1.2 Variation of shear slip coefficient Kr .............................................. 128 6.3.1.3 Variation of normal separation coefficient Kz ................................ 131 6.3.2 Imperfect bonding with finite length .......................................................................... 134 6.3.2.1 Verifications ................................................................................... 138 6.3.2.2 Variation of shear slip coefficient kr ............................................... 142 6.3.2.3 Variation of normal separation coefficient kz ................................. 145 Chapter Seven Summary and Conclusions ................................................... 148 7.1 Summary of the Research .......................................................................... 148 7.2 Conclusions ................................................................................................ 148 7.3 Suggestions and Future Study ................................................................... 154 Appendix A Mathmatica ® C for General Derivation ....................................... 157 Appendix B Mathematica® Code for Half Space ............................................ 159 Appendix C Mathematica® Code for Point Loading ....................................... 161 Appendix D Mathematica® Code for Two-layer Half Space ........................... 163 Appendix E Mathematica® Code for Two-layer on Rigid ................................ 166 Appendix F Mathematica® Code for Shear Loading ...................................... 169 Appendix G C++ Source Codes for Numerical Calculations ........................... 172 Appendix G1 C++ source code for half space ..................................................... 173 Appendix G2 C++ source code for single layer bonded to rigid base ............. 181 Appendix G3 C++ source code for single bonded to elastic half space .......... 189 Appendix G4 C++ source code for two-layer bonded to rigid base ................. 201 Appendix G5 C++ code for single layer bonded to an elastic half space subjected to shear loading ...................................................................................... 204 Appendix G6 C++ source code for single layer bonded to an elastic half space with shear slip and normal separation theories ................................................... 218 Bibliography ...................................................................................................... 232 vi List of Tables Table 1 Summary of axisymmetric indentation interface stress distribution and its Hankel transforms ....................................................................................... 37 Table 2 Material properties for steel, aluminum and nylon. ................................ 64 Table 3 Material properties used for thickness study ............. I ............................. 67 Table 4 Material properties for friction study ..................................................... 119 Table 5 Material properties for friction study ..................................................... 126 Table 6 Material properties for friction study ..................................................... 128 Table 7 Material properties for friction study ..................................................... 131 Table 8 Some values for different imperfect bonding length ............................. 137 Table 9 Material properties for friction study ..................................................... 139 vii LIST OF FIGURES Figure 1 Boundary condition of the indentation .................................................... 9 Figure 2 Stresses acting on an element of a solid revolution. ............................ 16 Figure 3 Half space with axisymmetric loading p(r). ........................................... 29 Figure 4 Single layer situated on rigid base: ...................................................... 38 Figure 5 Elastic layer situated on elastic base. ................................................... 42 Figure 6 Two-layer composite on rigid base. ...................................................... 46 Figure 7 Gaussian point integration algorithm. ................................................... 52 Figure 8 Numerical solution vs. exact solution for elastic half space. ............. 53 Figure 9 Comparison of normal stresses between case (1) and case (2) ........... 55 Figure 10 Comparison of shear stress between case 1 and case 2. .................. 55 Figure 11 Comparison of radial displacements between case (1) and case (2). 56 Figure 12 Comparison of normal displacements between case (1) and case (2). .................................................................................................................... 56 Figure 13 Comparison of normal stresses between case (1 ) and case (3) ......... 57 Figure 14 Comparison of shear stresses between case (1) and case (3) ........... 58 Figure 15 Comparison of radial displacements between case (1) and case (3). 58 Figure 16 Comparison of normal displacements between case (1) and case (3). oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo Figure 17 Comparison of normal stresses between case (2) and case (4) ......... 60 Figure 18 Comparison of shear stresses between case (2) and case (4) ........... 60 Figure 19 Comparison of radial displacements between case (2) and case (4). 61 Figure 20 Comparison of radial displacements between case (2) and case (4). 61 Figure 21 Normal stress distributions based on different materials. ................... 65 viii Figure 22 Shear stress distributions based on different materials. ..................... 65 Figure 23 Radial displacement distributions based on different materials. ......... 66 Figure 24 Normal displacement distributions based on different materials. ........ 66 Figure 25 Normal stresses at the contact surface. ............................................. 68 Figure 26 Normal stresses at 0.25 mm below the contact surface. .................... 68 Figure 27 Normal stresses at 0.5 mm below the contact surface. ...................... 69 Figure 28 Normal stresses at 1.0 mm below the contact surface. ...................... 69 Figure 29 Shear stress distribution at the contact surface. ................................. 71 Figure 30 Shear stress distribution at 0.25 mm below the contact surface ......... 72 Figure 31 Shear stress distribution at 0.5 mm below the contact surface ........... 72 Figure 32 Shear stress distribution at 1.0 mm below the contact surface ........... 73 Figure 33 Radial displacement distributions at the contact surface. ................... 73 Figure 34 Radial displacement distributions at 0.5 mm below contact surface... 74 Figure 35 Radial displacement distributions at 1.0 mm below contact surface... 74 Figure 36 Normal displacement distributions at the contact surface ................... 75 Figure 37 Normal displacement distributions at 0.5 mm below contact surface. 75 Figure 38 Normal displacement distributions at 1.0 mm below contact surface. 76 Figure 39 Normal stress distributions at the contact surface. ............................. 77 Figure 40 Normal stress distributions at z=0.5 mm from contact surface. .......... 78 Figure 41 Normal stress distributions at z=1.0 mm from contact surface. .......... 78 Figure 42 Shear stress distributions at contact surface. ..................................... 79 Figure 43 Shear stress distributions at z=0.5 mm from contact surface. ............ 80 Figure 44 Shear stress distributions at z=1.0 mm from contact surface. ............ 80 ix Figure 45 Radial displacement distributions at the contact surface. ................... 81 Figure 46 Radial displacement distributions at z=0.5 mm for contact surface. 81 Figure 47 Radial displacement distributions at z=1.0 mm from contact surface. 82 Figure 48 Normal displacement distributions at the contact surface ................... 82 Figure 49 Normal displacement distributions at z=0.5 mm from contact surface. .................................................................................................................... 83 Figure 50 Normal displacement distributions at z=1.0 from contact surface ....... 83 Figure 51 Normal stress distributions at contact surface. ................................... 85 Figure 52 Normal stress distributions at z=0.5 mm from contact surface. .......... 85 Figure 53 Normal stress distributions at z=1.0 mm from contact surface. .......... 86 Figure 54 Normal stress distributions at z=2.0 mm from contact surface. .......... 86 Figure 55 Shear stress distributions at contact surface. ..................................... 87 Figure 56 Shear stress distributions at z=0.5 mm from contact surface. ............ 87 Figure 57 Shear stress distributions at z=1.0 mm from contact surface. ............ 88 Figure 58 Shear stress distributions at z=2.0 mm from contact surface. ............ 88 Figure 59 Radial displacement distributions at the contact interface. ................. 89 Figure 60 Normal displacement distributions at z=0.5 mm from contact surface. .................................................................................................................... 89 Figure 61 Radial displacement distributions at z=1.0 mm from contact surface. 90 Figure 62 Radial displacement distributions at z=2.0 mm from contact surface. 90 Figure 63 Normal displacement distributions at the contact surface ................... 91 Figure 64 Normal displacement distributions at z=0.5 mm from contact surface. .................................................................................................................... 91 Figure 65 Normal diplacement distributions at z=1.0 mm from contact surface .92 Figure 66 Normal displacement distributions at z=2.0 mm from contact surface. .................................................................................................................... 92 Figure 67 Normal stresses at the bonding interface for nylon/aluminum. ........... 94 Figure 68 Shear stresses at the bonding interface for nylon/aluminum .............. 94 Figure 69 Radial displacements at the bonding interface for nylon/aluminum. 95 Figure 70 Normal displacements at the bonding interface for nylon/aluminum. .95 Figure 71 Normal stresses at the bonding interface for aluminum/nylon. ........... 96 Figure 72 Shear stresses at the bonding interface for aluminum/nylon. ............. 97 Figure 73 Radial displacements at the bonding interface for aluminum/nylon. 97 Figure 74 Normal displacements at the bonding interface for aluminum/nylon. .98 Figure 75 Force vs. deflection for single layer on rigid base. ........................... 104 Figure 76 Stiffness vs. Young’s modulus for single layer on rigid base. ........... 105 Figure 77 Relationship between stiffness and h/a. ........................................... 106 Figure 78 Force vs. deflection for single layer on elastic base. ....................... 107 Figure 79 Stiffness vs. the ratio of Young’s moduli. .......................................... 108 Figure 80 Normalized stiffness vs, the ratio of Young’s moduli. ....................... 108 Figure 81 Stiffness vs thickness changes ......................................................... 109 Figure 82 Force vs. deflection for two-layer composite on rigid base. .............. 111 Figure 83 Relationship between stiffness and the ratio of Young’s moduli. ...... 112 Figure 84 Stiffness vs. thickness ratio of the two layers. .................................. 113 Figure 85 Two-layer half space with shear loading ........................................... 115 Figure 86 Normal stresses at bonding interface between layer 1 and 2 ........... 119 Figure 87 Shear stresses at bonding interface between layer 1 and 2 ............. 120 Figure 88 Radial stresses at the bonding interface between layer 1 and layer 2 .................................................................................................................. 120 Figure 89 Displacements at bonding interface of layer 1 and 2 ........................ 121 xi Figure 90 Displacements at bonding interface of layer 1 and 2 ........................ 121 Figure 91 Shear slip theory of two layer half space .......................................... 123 Figure 92 Normal stress verification ................................................................. 126 Figure 93 Shear stress verification .............. . .................................................... 127 Figure 94 Shear displacement verification ........................................................ 127 Figure 95 Normal displacement verification ...................................................... 128 Figure 96 Normal stresses vs. shear slip coefficient kr ..................................... 129 Figure 97 Shear stresses vs. shear slip coefficient kr ....................................... 130 Figure 98 Radial displacement Ur vs. shear slip coefficient kr .......................... 130 Figure 99 Normal displacement Uz vs. shear slip coefficient kr ........................ 131 Figure 100 Normal stresses vs. normal separation coefficient kz ..................... 132 Figure 101 Shear stress vs. normal separation coefficient kz .......................... 133 Figure 102 Shear displacement Ur vs. normal separation coefficient kz .......... 133 Figure 103 Normal displacement Uz vs. normal separation coefficient kz ........ 134 Figure 104 Two-layer composite of imperfect bonding with finite length .......... 134 Figure 105 Shear slip coefficient kr (é) distribution. .......................................... 136 Figure 106 Shear slip coefficient kz(f,) distribution ............................................ 136 Figure 107 Shear slip coefficient vs. imperfect bonding length ......................... 138 Figure 108 Normal separation coefficient vs. imperfect bonding length ........... 138 Figure 109 Normal stress for all three cases. ................................................... 140 Figure 110 Shear stress for all three cases. ..................................................... 140 Figure 111 Displacement Ur for all three cases. ............................................... 141 Figure 112 Displacement Uz for all three cases. .............................................. 141 xii Figure 113 Normal stress vs. shear slip coefficient kr. ..................................... 143 Figure 114 Shear stress vs. shear slip coefficient kr. ....................................... 143 Figure 115 Radial displacements Ur vs shear slip coefficient kr ....................... 144 Figure 116 Normal displacement Uz vs shear slip coefficient kr ....................... 144 Figure 117 Normal stresses vs. normal separation coefficient kz. .................... 145 Figure 118 Shear stresses vs. normal separation coefficient kz. ...................... 146 Figure 119 Shear displacements vs normal separation coefficient kz. ............. 146 Figure 120 Normal displacements vs. normal separation coefficient kz. .......... 147 xiii CHAPTER ONE INTRODUCTION 1.1 Introduction 1.1.1 History Contact mechanics is a subject that has been discussed for more than a century. Early in 1882, the first academic paper, “On the Contact of Elastic Solids”, was published by Heinrich Hertz, With the increasing demand in the engineering development in railway and marine reduction gears, rolling contact of bearing industry and construction engineering, extensive efforts were made to the contact mechanics study, along with other engineering mechanics theory. Hertz’s theory was restricted to frictionless surfaces and perfect elastic bodies. Over the last half century, research in contact mechanics has focused on the removal of these restrictions for the pure contact problem. At the same time, development of the theories of plasticity and linear visco-elasticity paved the road to investigate the stresses and deformations at the contact of inelastic bodies. With the advances of composite manufacturing technology and Its extensive usage in the aerospace and automotive industry, contact mechanics was applied to composite study in the past few decades. 1.1.2 Summary of contact mechanics research The Boussinesq-Flamant problem was the very first contact problem that was approached with theoretical mechanics. It was defined as a half space loaded by a point load. Flamant (1892) obtained a solution by way of a three dimensional solution for normal loading on a straight boundary. Boussinesq (1892) extended the solution to the case of an inclined force, which can be visualized as a combination of normal and tangential loadings. Dean, Parsons and Sneddon (1944) presented a theoretical approach to semi- infinite elastic solid with defined stress applied at the interior of the solid. The direction of the loading conditions is perpendicular to the surface of the semi- infinite solid. In order to simplify the mathematic work, the authors had to assume the material is incompressible and the loading is uniform stress. They also investigated the compressible elastic solid with numerical method. They found that by using a single numerical factor, a fair accuracy can be achieved for the general case from the equivalent incompressible medium. During the course of investigation of the distribution of stress and displacement in elastic solids, it used to be carried out through guessing appropriate combination of solutions, which satisfy the prescribed boundary conditions in any special case. Love (1939) used this cumbersome method in the case of a conical punch. With some experienced engineers’ pioneer work on rigid indention, Harding and Sneddon (1945) noticed that a systematic application of integral transforms into the rigid indentation of semi-infinite elastic solid will reduce the problem into a pair of integral equations, which fall into the classic mathematical cases studied by Titchmarsh (1937) and Busbridge (1938). In the Harding and Sneddon study, they established the foundation for the mathematical solution for axial symmetric indentation problem. Though the development is limited to the semi-infinite elastic solid, the authors pointed out the possibility of extending the method to the solution for the stress distribution in a plate with finite thickness. Sneddon (1951) presented a solution for a finite thickness plate perfectly bonded to a rigid support. He used a symmetrical approach, mirrored the material and the loading along the bonding interface with a rigid support. Therefore, by solving a finite plate with doubled thickness and symmetrical loading at both top and bottom surfaces, it leaded to the solution equivalent to the single thickness plate perfectly bonded to a rigid base. Later, Sneddon used the theory developed by Harding and Sneddon (1945) to successfully obtain solutions for the Boussinesq problem with flat-ended cylindrical (1946) and rigid cone indenters (1948). With more and more theoretical development driven from the industry, more and more research was focused on sliding contact. After the half-space point solution by Boussinesq and Cerutti, it seemed that integration of the point loading condition would lead to all other loading conditions. However, this direct approach normally leads to a series of intractable integrals with mathematical difficulties. Green (1949) introduced a method for the stress analysis of normally loaded half space, which can be extended to the shear loading of the half space. Based on the Green study, Hamilton and Goodman (1966) developed the equations for the complete stress field due to a circular contact region carrying Hertzian normal pressure and a proportional distributed shearing traction. The shearing traction is represented by the normal pressure multiplied by the frictional coefficient. A complex harmonic function was used in their investigation. The solution for the stresses inside the indented material was developed to illustrate the potential failure point. Stress boundary conditions were specified for the half space solid in their research. Based on some Russian pioneer work, Meijers (1968) developed a symbolic solution for the contact problem of a rigid cylinder pressed on an elastic layer connected perfectly to a rigid base. Hertz contact boundary condition was assumed in his development with no friction between the cylinder and the layer for the generalized plain stress problem. Meijers focused on developing a symbolic solution for both thin and thick layer conditions. He used a basic integral equation for the pressure distribution at the contact interface presented by Aleksandrove (1962) and Koiter (1963), by introducing a complex kernel function for the contact stress and displacement and expanding the kernel function into a power series. The first three coefficients of the power serials were evaluated by using numerical integration for different Poisson ratios. With this approximating work, Meijers was able to develop a general solution for both thin and thick single layer indentations. The details of the numerical processing are quite cumbersome and it is very difficult to extend the approach to layered composite due to the complexity of mathematics. Pao et al (1971) studied the upper and lower bounds of the maximum stress of an elastic layer bonded to a half space elastic foundation by a rigid cylinder. The symmetrical loading technique from Sneddon’s (1951) was used for his calculation. A combination of half space and single layer solutions from Sneddon’s development was used in this investigation. Pao used a segment approach to divide pressure distribution into many small intervals. By assembling a matrix from each interval, he was able to come up with a series of functions to gain a numerical solution. Different boundary conditions, frictionless and rigid bonding, between the layer and rigid base were investigated. The stress and displacement component was evaluated to demonstrate the effects from position ratio, friction between layer and the rigid support, and layer thickness. Ting (1966) presented symbolic and numerical solution of the contact stresses between a rigid indenter and a viscoelastic half-space, which extended the research from elastic media to viscoelastic media. In his approach, he introduced a method which can solve an arbitrary axisymmetric indention on viscoelastic half space. A dynamic rebound of a rigid sphere on viscoelastic half space was studied by Hunter (1960). A numerical method was used to evaluate the result by both authors. Dhaliwal (1970a) developed a solution for the single layer pressed against an elastic foundation. A different numerical method was used to develop the solution, which achieves good agreement. He also developed a symbolic solution for an arbitrary profile (1970b). However, no numerical result was presented in his further investigation. Gupta et al (1973 and 1974) studied stress distributions in plane strain of a layered elastic solid subjected to arbitrary normal and shear boundary loading. Fourier transform of Airy stress function was used in his investigation. His analysis result shows that internal stress distributions as obtained by elliptical contact pressure are valid for most cases. He also pointed out that the half width and maximum pressure on the surface may be quite different from Hertzian solutions in order to solve the problem with his approach. Matthewson (1981) presented a theory for the indentation of a soft coating by a rigid body. The theory is based on the assumption that the coating is perfectly bonded to the substrate, and behaves linearly elastically. The displacement of the contact boundary condition is approximately approached by a finite power series. The approximation used quantities averaged through the coating thickness. This assumption differs from many other publications, where an exact displacement boundary condition was met, and stress on the boundary problem was approached approximately. Jaffar (1993) used Chebyshev series method to study the surface deformation of a bonded elastic layer indented by a rigid cylinder. His result shows close agreement with the exact solution for the half plane problem with the derived asymptotic formula for a thick layer. Giannakopoulos (1997a, 1997b) investigated the indentation of solids with gradients in elastic properties. A power law material model was used in his approach. Due to the difficulties, commercial finite element code was used to study different indentation case. 1.1.3 Contact mechanics in composite material There are a lot of studies and publications in laminate composites study. Pagano ' (1970) presented a theoretical solution for laminate composites. Swanson (2004) studied the Hertzian contact loading of orthotropic materials. He used both procedures, developed by Willis (1966) and Pagano (1970) separately, and further developed a detailed numerical method to achieve the results for laminate composites. This research is to focus on the layered composite. Each layer has homogeneous material property. Literature on laminated composites is not reviewed in details. Layered composite indentation was very challenging in contact mechanics. In the application, a soft coating on a hard substrate is very efficient for protection. In automotive safety design, a soft layer on a hard supporting structure provides the necessary stiffness to meet the design requirements and also provide cushion to minimize the injury of vehicle occupant. To the author’s knowledge, few publications are available in this subject. The major theoretical study was limited to indentation of single layer bonded to an elastic half space. The majority of the industry currently uses finite element approach with commercial code, such as LS-DYNA®, Abaqus®, PAMCRASH®, MADAMO® for dynamic research, Abaqus®, Nastran®, ANSYS®, etc. for static study. Shield and Bogy (1989) studied a multiple region contact solution for a flat indenter on a single layer bonded to a half space. The flat indenter has sharp corners, which has singularity. Since it is a multi-region contact problem, each region has its own solutions. Jergensen et al (1998) studied spherical indentation of composite laminates with controlled gradients in elastic anisotropy. His study provides an insight of indentation of the indentation problem on orthotropic media. There is no known analytical method for the problem he was trying to solve. Finite element approach was used in his investigation. Lu and Liu (1992) presented an lnterlayer Shear Slip Theory for delaminating analysis, which was extended from a layer-wise laminate theory, the previous development by Lee and Liu et al (1991 ). The continuity condition of the interlaminar shear stresses on the composite interfaces was satisfied, which is only valid for analysis of laminates with shearing mode (mode It and mode III). In order to be able to handle the normal mode (mode I), Liu et al (1994) extended the lnterlayer Shear Slip Theory by introducing linear normal separation theory. Combining both linear shear ship theory and linear normal separation theory, the theory became robust to solve complex loading conditions. In this research, both shear slip theory and normal separation theory was implemented into the development as part of the advanced study. There is tremendous literature regarding the classic topic of contact mechanics. The books authored by Gladwell (1980) and Johnson (1985) provided rich reference literature on contact mechanics. Only some critical publications are reviewed here, which represent major progresses of the research and theory development. 1.1.4 Approach to layered composites In order to solve the contact problem, simplifying boundary conditions was critical for mathematical success. For a rigid indenter, the contour of the indenter shape defines the deforming shape contact interface. Specifically, the displacement at the contact interface can be defined, and the normal tractions outside the contact area disappear. It seems that the proper boundary conditions have to be mixed displacement and traction. In order to solve this boundary value problem which involves biharmonic equations, certain transformation is needed. With the mixed boundary condition, mathematical difficulty for the theoretical approach to the contact problem rises. V Figure 1 Boundary condition of the indentation In this research, the focus is on contact mechanics of the axially symmetric indentation. The general setup of the problems is illustrated in Figure 1. Hankel transformation is used in order to solve the biharmonic equation. The details of Hankel transformation are discussed in Chapter Two. 1.2 Object of this study In automotive safety development, head injury is based on the acceleration level. In crash accident, the human head contacts the vehicle interior trim. The human head is a complex structure, which can be simplified as a rigid skull and a deformable skin. In recent development, the dummy head was optimized as a vinyl skin and polyurethane skull structure, which represent the human structure more closely. In these developments, due to the complexity of the mathematical modeling, there is no thorough theoretical study of the structure. All the research work done in the last three decades was based on experiments and experience, though a lot of papers were published. Most of the research in the industry was done on experiments and finite element analysis. The vehicle interior trim normally is constructed by sandwich composite with a metal supporting base. The sandwiched composite is made of foam and plastic materials, such as Ensolite®, rubber foam, urethane, polyurethane, thermoplastics, etc. There were a lot of studies on the structure design, but few literatures from contact mechanics study can be found on indentation study of layered composite. The author’s early experimental study (Liu, Dang and Wang, 1996 and1997) of foam composite under low velocity impact was the initial drive that leads into the research. It was realized that the theoretical studies for the contact mechanics of layered composites are needed to better understand the performance of the composite materials. This dissertation develops a theoretical approach to this problem to provide some design guidance for the composite structures. From the literature review, most of the research is limited to two dimensional, semi-infinite and single layer materials. To the author’s knowledge, the results of indentation to the layered composites have not been seen. The objective of this study is to investigate the layered composites subjected to indentation by an axisymmetric indenter. A new technique is proposed and significantly reduces 10 the complexity of the mathematics, which provides a unique way to look into stress and displacement distribution. The final goal is to develop a model for the composite indentation that can be used as guidance for the design. Also shear slip and normal separation theory are incorporated into this technique to study the bonding interface. An imperfect bonding with finite length was studied with a step function for the shear slip coefficient, which provides insight into the local debonding of the composite material. 1.3 State of the problem This dissertation will study the layered composite indented by an axisymmetric indenter. An elastic theoretical study of the indentation was carried out in this study. Studied cases include a half space, one layer bonded to rigid base, one layer bonded to an elastic half space, and two layer bonded to rigid base. The shear slip theory was incorporated into the above development and investigated in this dissertation as well. 1.4 Contents of this dissertation This dissertation includes seven chapters. Chapter One reviews the studies done in the past, and states the subjects that are to be investigated in this study. Chapter Two discusses the basic theory developed for this study. In this chapter, we develop the basic theoretical derivation, discuss the Bessel function used in the study, and the Hankel transform that is the core part of the mathematical base for this study. We also discuss the exact solution of some simple cases. 11 Lastly, we discuss the details of the approach with a stress boundary condition for the contact problem. Chapter Three discusses the detail of a numerical approach to the final solution for complicated cases. Mathematica® was used in this approach to reduce the mathematical work significantly. Microsoft® C++ programs were developed for the numerical calculation in the final result for this investigation. Verifications were carried out among these derivations. For example, the numerical result is verified against the theoretical solution for point loading. By increasing the thickness of the first layer, a two-layer case would be simplified to a one-layer case, etc. Chapter Four discusses some basic applications. Six different cases are investigated in this chapter: Case 1, material study: with half space and point loading, we investigated aluminum, steel, nylon material to understand the stress and displacement for different material under the same loading condition. Case 2, thickness study: two layer half spaces are used to investigate the thickness effect. The composite consists of an aluminum layer bonded perfectly with a very thick nylon base. Spherical indentation is investigated during the study. 12 Case 3, boundary condition study: in this case, a two layer half space with a combination of aluminum/steel, aluminum/nylon, aluminum/rigid are investigated, using spherical loading. Case 4, loading conditions: an aluminum half space with different loading conditions, point loading, uniform pressure, flat indentation, and spherical shape indenter are studied. Case 5, lamination study: two different cases are investigated — material order and thickness ratio. For material order study, the thickness ratio is set to constant, and the composites of Nylon/Aluminum and Aluminum/Nylon were investigated with spherical indentation. For thickness study, with a fixed material order, both Nylon/Aluminum and Aluminum/Nylon, the thickness ratios are presented in this chapter. Case 6, Shear slip theory: This was studied for two-layer composted indented against a rigid base. Chapter Five presents models with the techniques developed in the previous chapters. The relationship of force, deflection, Young’s modulus, thickness, and contact radius is incorporated into simple models. These models provide simple and valuable information as guidance to design and applications for layered composite. Chapter Six is presents some insight into the advanced study. The pure shear traction at the indentation interface is presented. With friction coefficient to relate 13 normal traction and shear traction, the solution to a more general contact problem can be obtained. Moreover, the shear slip theory and normal separation theory are introduced into the technique developed in the previous chapters. With the shear slip coefficient being set to a step function with Hankel transform, we investigate the limited imperfect bonding for the composite under indentation. Chapter Seven summarizes the results and concluded the study of this dissertation. Some insights into the theory are discussed and suggestions to future study were presented. Lastly, bibliographies and appendices provide references and the information that can’t be covered in detail in the main section of this dissertation due to limited space. 14 CHAPTER TWO DEVELOPMENT OF THEORETICAL ANALYSIS In this chapter, we develop a new technique that will be used subsequently to solve contact problems consisting of the layered composites. Validations of the technique are presented as they are necessary for the justification of later applications. For the simple case of point loading, the new technique leads to the same solution as otherwise presented by Sneddon (1953). For other loading cases, such as uniform pressure, flat contact and spherical punch, the stress and displacement solutions are in integral forms and explicit expressions are not found. Therefore, numerical methods are sought to investigate these cases. 2.1 Elasticity Analysis Based on Cylindrical Coordinate System Based on the free-body diagrams shown in Figure 2, we can establish the following equilibrium equations (Timoshenko and Goodier, 1970), which are also the governing equations for cylindrical bodies under static loading. Note that both equations are independent of i9 and the body force is not considered in the formulation. 60’ +61” +07 —06 =0 (2-1a) 0r 6r r Qflfigfizzo (2-1b) 6r 82 r 15 Figure 2 Stresses acting on an element of a solid revolution. Since linear elastic analysis for isotropic bodies is of interest in this research, the following Hooke’s law is imposed a, = $0. -vJv(xa)da <2-10c) 0 where Jv(xa) denotes a Bessel function of order v of the first kind, i.e. x)=n20_( 1) n(x/2)2n+v n!F(n+v+1) where F(n +1) = n! is the Gamma function when n is an integer. The function f(a) defined by Equation (2-10b) is called the v’h -order Hankel transform of function f(x). Equation (2-10c) is called the Hankel inversion theorem. If we let x=r and 01:6, then Equation (2-10b) can be rewritten as follows 7(5) = jrf(r)Jo(r§)dr (2-10d) 0 The Hankel transformation has the following property w jai— j—{—)Jo(r:)dr - -4: flé) (2-11) 0 V4¢ can be expressed as V2(V2¢). If we take the Hankel transform of Equation (2-8) and use the relationship of Equation (2-11), then we have 21 er W0(§r)dr = (—z——: ) er0(r§)dr (2-12a) 0 dz 0 oo 2 oo IrV2(V2¢)Jo(§r)dr = (5—2—42) er2¢/o(r<§)dr (2-12b) 0 Z 0 Substituting Equation (2-12a) into (2-12b), we have m 2 00 JVV4¢10(§’)dr = (if - :212 IrWo(r§)dr = 0 (2'13) 0 dz 0 00 Let G(§,z) = Ir¢/0(r§)dr , i.e. the Hankel transform of ¢(r,z), then 0 d2 22 (—2--6 ) G(§.ZI=0 (2-14) dz The integration of Equation (2-14) is elementary and its solution can be expressed as G(éj,z) = (A + Ede—Z: + (C + Dz)ez§ (2-15) where A, B, C, and D are functions of 5. Now we can develop stress and displacement expressions based on Equation (2-15). 2.3 Displacement Solutions If we multiply both sides of Equation (2-9a) by rJ0'(r§) and integrate over r from 0 to co, we have the following equation after integration by parts 22 J‘ru,.11(r§)dr = 0 1+2: fig dz lnverting this result by using the Hankel inversion theorem, i.e. Equation (2-10c), we obtain the following expression (X) 1+ dG ur= “ jéz—Jlrred: (2-16a) it 0 dz Similarly, from Equation (2-9b) we have 00 2 d G 2+2 jruzJo(r:)dr-= 2 — ”.520 0 dz lnverting this result by using the Hankel inversion theorem, i.e. Equation (2—10c), we have u. = lei-1" 26— ’1’” ” €20)J0(r~f)d§ (2-16b) 0 dz 1” 2.4 Stress Solutions In a similar manner, we can find the expressions for stress components. Multiplying both sides of Equation (2-6a) by rJ0(r§), Integrating them over r and inverting the equation by using the Hankel inversion theorem, i.e. Equation (2- 10c), we find 23 °° 3 0. =j51u+md G d 3 —(3l+4#)§2%g]Jo(r€)dé (MM 0 Z In a similar fashion, multiplying both sides of Equation (2-6d) by U 1 (r5), integrating them over r and inverting the equation by using the Hankel inversion theorem, Equation (2-100), we find 00 2 rrz = 1521/1 d f +(2 + 2/1)§ZG]J1(r§)d§ mm 0 dz Obtaining the expressions for a, and 0'6 is not as simple as the components we just obtained. If we add Equation (2-63) to Equation (2-60), then we have 2 6 or +09 =[2(/I+#)-a—2—-2#V2]¢z Z Multiplying both sides of this equation by rJ0(§r), and integrating them with respect to r from 0 to 00, we obtain °° 3 2 d G 0’ d0 jrrar +ae)Jo(ré)dr=2(/1+m 3 —2#(—2 —:2)— 0 dz dz dz lnverting the above equation by using the Hankel inversion theorem, i.e. Equation (2-100), we obtain °° d3G d or +09 = 16(1—7+#§2—Q)J0(r§)d¢f (2-16e) 0 d dz Z 24 Now let us rewrite the equilibrium equation, Equation (2-1a), in the following form 62-020,) = r(0',. +09)—r2 E r Substituting Equation (2-16d) and (2-16e) into the above equation, it yields 2dG —(r 0r)=2r1§(1——+#§d7)Jo(r€)d§-r2j§[1;§;2+é(1+2#)—]Jo(r§)dcf Integrating the foregoing equation with respect to0 r, interchanging the order of the integrations and using the integrals eroogur = 111 (r6) and 0 5 jr2J11r6idr = 2%J1(r5)‘ —’-Jo. 6 0 5 we obtain =j61/1— + 6 (J1 + 33,t1)—I~/o(ré‘)d§-2(li + —"—) j62 i’g—Ji (r6>d6 (2-16f) dz3 r 0 dz From Equations (2-16e) and (2-16f), we have ‘62 d6 °° 3G 09 = ,1 j§(-:Z—3G- )—Jo(r r6>d6+ 2” + —”—) 0j62; dG J (r 61d6 (2469) 0 25 As discussed earlier, A, B, C and D of Equation (2-14) are functions of 5. Without losing any mathematical generality, we can let (1+v)(l—2v) E62 [(A+Bz)e‘25 +(C+Dz)e25];7(6) (2-17) 6*(692) : where 5(5) = er(r)J0(r§)dr is a function of 5, and J0 and J, are Bessel 0 functions of its first kind. Substituting Equation (2-17) into Equations (2-16), we have 0; = j6{1A6 + 8(1— 2v + z-6IIe‘i‘f +1—C6 + 00 — 2v — 2611er ifi<6>Jod6 O (2-18a) rrz = j{1A6 — B<2v — 2611625 +166 + D<2v + 26>]er 16(6)J1(r6)d6 0 (2-18b) or = j6{1—A6 + B<1+ 2v — 26m”: + 1C6 + D<1+ 2v + 26%:5 }fi(<§)Jo(r6)d§ 0 — j {ti—A6 + 8(1— 26)]e'25 + [C6 + 00+ 2611er }%17(6)J1 (r6166 0 (2-18c) 26 09 = j2v6(Be‘25 + De“ )17(§)Jo(r§)d§ 0 §(§)J1(r§)d§ ‘r|v-— + j{[(—A6 + 3(1- 26)]e'24‘ + [C6 + D(1+ 26)]e2‘5} O (2-18d) u, = jg—i'gEu—M + 3(1— 26)]e‘25 + [C6 + D<1+ z6>1ezfm6>J1d6 0 (2-18e) u2 = ojl%fl’—){1A6 + 8(2 — 4v + 26m”: + 1C6 + 0(1 + z6)1ezfn3(6)Jo(r6)d6 0 (2-18f) The details can be found from appendix A. For simplicity, let us define a; = [A6 + 30- 2v + 2.6)]e‘24‘ + [—C6 + D(1— 2v — 26)]ez‘5 (2-19a) r; = [A6 — B(2v — 26H?“ + [C6 + D(2v + 2.6)]ez‘f (2-191)) u;=1::"{[—A6 + 30- 26)]63‘24r + [C6 + 0(1 + 26)]e24} (2-190) u’ = —1—2—"{[A6 + 3(2 — 41/ + 26)]?25 +[C6 + D(1+ 2.6)]e24’ } (2-19d) Z Since ,u = , then Equations (2-190) and (2-19d) become E 2(1 +v) 27 2w: = [-Aé‘ + 30- 26)]e‘zf + [C6 + D(1+ 26)]ezs‘ (2—19e) — 2w: = [A6 + 3(2 -4v + Zak—24‘ + [06 + D0 + 26ml? (2-19f) Accordingly, the stress and displacement components can be expressed by the following relations. 07. = Edméflwfidé (2-20a) rm = ETLfléNoUfldcf (2-201,) u,- = °(:jui625(6)J1(r6)d6 (2-20c) u. = Lju2m61Jo (r6)d6 (2-20c) The use of Hankel transformation greatly simplifies the expressions of stress and displacement components and has significant benefit of reducing mathematical complexity in imposing boundary conditions. In the following sections, we will apply these results to a few special contact problems of the layered composites consisting of perfect bonding interfaces. 28 J W I. “. —LLII 2.5 Exact Solution for Half Space under Point Load The solution for a half space subjected to a point load was developed by Sneddon (1953). The same problem is solved below to validate the new technique. ”IV Figure 3 Half space with axisymmetric loading p(r). Since both stress and displacement components have to be zero when 2 approaches infinity, C and D must vanish in Equation (2.17), i.e. (l+v)(1—2v) 2 (A + lane—25 5(6) (221) 1‘35 G*(~§,Z)= Equations (2-18) can then be simplified as 29 a. = j6{1A6 + Bo - 2v + z6>1}e‘“‘ 17(6)Jo(r6‘)d€ (2-22a) 0 6.2 = 11A6 — B<2v — z6)1e‘256(6>Jo(r6)d6 <2-22b) 0 a, = 161-A6 + B<1+ 2v - z6)1e'29‘;7(6)Jo(r6)d6 0 (2-22c) — Zia—:46 + 8(1— z6)1e‘2‘=‘ £6<6>J1 (r61d6 0, + 019 = 361-146 + B<1+ 4v — 261163f 17(§)Jo(ré)d€ (2-22d) u, = GEE—£33146 + 8(1- z6)1e‘zffi(6)J1(r6>d6 (2-22e) uz = gig-Vlmwa—m 26>1e‘2513(6)Jo(r6>d6 (2-220 Similarly, Equations (2-19a), (2-19b), (2-19e) and (2-19f) can be simplified as a; =[A6+ B(1-2v+z§)]e‘z9‘ (2-23a) r; = [A6 — B(2v — 26)]e‘24’ (2-23b) 2m: =[—A6+B(1—26)]e‘z§ (2-23c) 30 — 2w: = [A6 + 3(2 ~4v + Zak-29‘ (2-23d) With the use of Equation (2-17), the boundary conditions can also be simplified. This provides us with a great advantage of solving the contact problems of layered composite, i.e. oz = :1602 17(6)Jo(r¢f)d§ (2-25a> rrz = :IT;17(6)Jo(r€)dé (2-251,) ur = Es‘u:§(€)J1(ré)d6 (2-25c) uz =§u;l_7(§)~10(r§)d§ (2-25d) The boundary conditions for the contact problem defined in Figure 3 are az|z=0 = p(r) when r < a (2-26a) az|z=0 = 0 when r > a (2-26b) 1,212: 0 = 0 (2-26c) 31 If the indenter has a smooth edge, the normal stress component vanishes at the edge. However, if the indenter has a sharp edge, there is singularity around the edge. From Equations (2-26) and (2-25), with Hankel inversion theorem, Equation (2- 10c), we have a: 0, =1 (2-27a) z=O a: TI’Z : 0 (2-27b) z=0 With the above boundary condition and Equations (2-23a) and (2-23b), we have A; + 80— 2v) =1 (2-28a) A if — 2vB = 0 (2-28b) Solving the above equations, we have Therefore, function G defined in Equation (2-21) becomes C(62) = (H ”5):; 2”) (2v + z6)e‘zf 13(6) 32 Substituting the solutions A and B into Equations (2-22), we can find all stress and displacement components in the half space. Details of the analysis are given in Appendix B. 02 = 1(1 +26)€e‘z‘5fi(6)Jo(rcf)d6 (2-29a) 0 rm = 1262625 17(§)J1(rcf)d§ (2.291)) 0 a) = 1(1- z6)6e‘2517(6)Jo (r6)d6 +§ j(—1 + 2v + Z§)e“25fi(é)~/1(ré)d6 0 0 (2-29c) 019 = 121/662‘ .5(6)Jo(r6)d6 —} j (—1 + 2v + z6)e‘zi 13(6)J1 (r6)d6 0 0 (2-29d) u) = 10+”1 gzv'ziie-zé p(6)J1(r6)d6 (2-29e) 0 u. = I (I + v)”; 2" — “5) e‘ziz—2(6)Jo(r6)d6 (2-29f) 0 For near point loading, the normal stress distribution on the contact surface found by Sneddon (1951) is 33 f.)__3_ (2 30a) 27! _3_ (a2 + r2)2 p(r) = (- The Hankel transform of p(r) is (I) 6(6): jrp(r)Jo(r6)dr=(—f—)e““5 (2-301)) 0 71' Substituting Equation (2-30b) into Equations (2-29), we have all stress and displacement components. Details of analysis can be found in Appendix C. 2 2 az:_(_21:)(3z+a)(z+a) +ar (2-31a) 7r __ ["2 +(a+z)2]2 Trz = 4%) "3rz(a+z) (2-32b) [r2 +(a+z)2]§ 0'6? =—(—2£—)[-l-(l-2v)+ (—I+2v)(a+z) + 2vr(a+z)—rz3] (2-330) 70' r r[r2 +(a+z)2]5 [r2 +(a+z)2]5 a',=—P—[l(1—2v)+ (-l+2v)(a+z) _ (a+22)r +rz[2(a+z)2—r2]] 2717' r 5 _ _3_ _ r[r2+(a+z)2]2 [r2+(a+z)2]2 [r2+(a+z)2]2 (2-34d) z (l+v)P 1 a” 22:5- [ (l—2v)(a+z) + rz 3] (2-356) r (—l+2v)+ r[r2 +(a+z)2]§ [r2 +(a+z)2]§ 34 u =(1+V)P (5-4v)az+(3—2v)22+2(1_V))(02+r2) z 27$ 1 (2-36f) [ 3 [r2 +(a+z)2]2 These solutions are identical to the solutions presented by Sneddon (1951 ). 2.6 Traction Boundary Conditions for Contact Problems The most challenging part of using the new technique to solve general contact problems is the integral involved in Equations (2-29). A transformation technique will be necessary to reduce the complexity of integration. However, the type of boundary condition is also of a major concern as the inversion of the transform will become very difficult, if not impossible, if mixed boundary conditions are imposed. In other words, either displacement or traction should be used in each contact problem to reduce the complexity of the dual integrations. Moreover, in reality, all contacts occur with limited contact surfaces and the displacements are unknown in the non-contact areas adjacent to the contact surfaces, even though we can assume there is no deformation at areas far away from the contact areas. Therefore, specifying displacements is not feasible in solving contact problems. On the contrary, since there is no traction on the surfaces of non-contact areas, it automatically leads to specifying traction boundary conditions for studying contact problems. The focus perhaps should be placed on how to define the tractions on the contact surfaces. As the tractions on the contact surfaces depend on geometry and surface condition of the contact problem, there is no easy rule to guess the traction on the contact surface. However, many pioneers have defined tractions with their experience and knowledge in mathematics and 35 physics. The traction boundary conditions have been refined to be close to reality for some commonly used indenters, such as those with point, flat, and spherical noses. Johnson (1985) summarized the pressure distribution for contact problems with a circular region. For a circular contact region with a radius a, it is required to find displacements at a surface point and stresses at an internal point due to the pressure distributed over the circular region. Solutions in closed form can be found for axisymmetric pressure distributions in the form of p = p0(l — r2 / a2 )" . When n = 0, it is a uniform pressure. When n = -%, it is for a flat punch. When n = %, it represents a spherical indention condition. Combining with a point load, 5(r) , we have defined the 2n*r which can be represented by a function of p(r) = — pressure distributions for indenters with regular shapes. We will use Hankel transformation in the following chapters. For point loading, Sneddon (1951) used p(r) = -i 0 instead of p(r) = — am for the ointloadin a roach, 22(a2+,2)3/2 27m 9 9 pp because it lead to singularity solution at the contact point with infinite stress when using function p(r) = — 250;) . With revised stress distribution If r function p(r) = - —}—)- a , it eliminates the singularity with more realistic resun. 36 The contact traction distributions and associated Hankel transforms for several contact problems are summarized in Table 1. Table 1 Summary of axisymmetric indentation interface stress distribution and its Hankel transforms Hankel Transform lndenter Contact Interface Stress 5(5) = Profile Distribution Total Load 00 jrp(r)Jo(r6)dr 0 Near p a P Point P(") = "— -P 5(5) = __e-0~§ Loading 2” (02 + ”2):”2 2n = — h 5 ”mm“ M m w en 0 rsa 2 ‘(6) — 3301 (a6) Pressure ‘ ”a P0 p ‘ 6 1 p(r) = 0 when r>a r2 "I" p(r) = -Po(1- —) 2 Flat a 2 2 _ poasin(a§) Punch " 27661 PO 12(5): ’— when OSrSa é: Ip(r) = 0 when r>a r2 l p(r)=-Ivo(l-—2)2 _ p0 .. Spherical a __ 3,112 mg) ="—3— lndenter 3 P0 05 When OSrSa [05 (308(06) _ sin(a§)] P(r) = 0 when r>a Shape p(r) p0 W e f' I? _ ”(r22 _ r12)p0 P(€) 5 Uniform Pressure p(r) = 0 when rr2 [r2J1(r26) - nJ1(r16)] 37 CHAPTER THREE NUMERICAL SOLUTION AND VERIFICATIONS In this chapter, we use the theory developed in Chapter Two to obtain solutions for several contact problems, eg. a single-layer material situated on a rigid base, a single-layer material situated on an elastic base and a two-layer composite situated on a rigid base. A numerical scheme is developed to calculate the integrations which cannot be executed analytically. Validations of the numerical scheme are carried out through the investigations. 3.1 Elastic Layer on Rigid Base A single elastic layer situated on a rigid half base is one of the most common contact problems in the real world, e.g. soft thin films and protection coatings on metals. In this section, we will develop a generic solution for this contact problem. For any loading condition as shown in Figure 4, we assume there is no friction on the contact surface and the bonding between the elastic layer and the rigid base is perfect. p(r) l I 1) I l 4, Elastic Layer I4—2a-0—J r /////9 /// / Rigid Base v 2 Figure 4 Single layer situated on rigid base. 38 Based on the above assumptions, we have the following boundary conditions: at z = 0 02 = p(r) (3-18) 6,, = 0 (3-1b) at z = h 142 = 0 (3-1c) u, = () (3-1d) Combining Equation (2-18a) with the boundary condition Equation (3—1a), we have ja?6j(6)Jo(r6)d6 = p(r) 0 If we apply Hankel’s inverse transformation to the above equation, we obtain the following relationship a; :1 (Ma) where a; is defined in Equation (2-17a). Similarly, combining Equations (2-18b) and (3-1b), we have 39 erzr(6)Jo(r6)d6 = 0 0 After applying Hankel’s inverse transformation to it, we have 1,2 = 0 (3-2b) where 2'; is defined in Equation (2-17b). Moreover, combining Equations (2-18d) and (3-1c), we have ju;fi(§)Jo(r6)d§ = 0 0 Using Hankel’s inverse transformation, we have “2 = 0 (3-20) where u; is defined in Equation (2-17d). As the fourth exercise, combining Equations (2-18c) and (3-1d), we have (I) jui6fi(6)J1(r6)d6 = 0 0 After using Hankel's inverse transformation, we have u, = o (3-2d) 4O where u: is defined in Equation (2-17c). As given in Chapter Two, we have defined (1 +v)(1—2v) G( , )= 62 E52 [(A + Bz)e_z‘5 + (C + one“ 115(6) If we use Equations (2-17), Equations (3-2) can be rewritten as follows. 1+B(2—v)+D(1—2v)+A6—C6=0 (3-4a) (A + C); - 2(3 — D)v = 0 (3-4b) Age—h?” -C.6e”5 (1 + v)— Be‘h‘f (1 — kg) — Dehé (1 +116) = 0 (3-40) Age—’75 + Céehé: + 36"“ (2 -4v+ h6)+ Dehé‘ (—2 +4v+ kg) = o (3-4d) The variables A, B, C and D can be found from solving the foregoing equations, i.e. = 28”; (—2+4(1 -e2”5 )v2 41262 +v(5—3e2h5 + 2115)) A (3-5a) 6(3+e"”f (3—4v)—4v+2e2hi(s—12v+sv2 +2h24‘2)) B : e2“ (—1 + e2” (—3 + 4v) + 2126) (3%) 3 +e4”9‘(3-4v)-4v+ 2e2h5(5 -12v+8v2 + 211262) 2((3 — 4v)v+ e2” (2 + 4v2 + 12262 + v(-5 + 216») (3-5c) — 6(3+e4”5 (3 —4v)—4v+2e2"5 (5-12v+8v2 + 212262)) 41 3 —4v+e2”5(1 + 2kg) (3-5d) 3+e4”5(3—4v)—4v+ 2e2”5(5 -12v+8v2 + 2h262) B:— 3.2 Elastic Layer on Elastic Base Shown in Figure 5 is an elastic layer situated on a semi-infinite elastic base. It can be considered as an extension from the previous case which is an elastic layer situated on a rigid base. p(r) ll . ll _ L 0 I T r Layer1 28“—’ h Layer2 V Z Figure 5 Elastic layer situated on elastic base. To simplify the boundary conditions, there is no friction at the contact surface and layer 1 is perfectly bonded to layer 2. The thickness of layer 1 is h and the contact radius is a. If we use superscript to represent the layer, the boundary conditions for this contact problem can be summarized as follows atz=0 42 09’ =p(r) (3-6a) :2.) = o (3-6b) at z = h 09’ =02?) (3-6c) A? = r)? (3—6d) 11(1) : “(2) (3-6e) .9) = .9) (s-er) The function of G(<§,z)for layer 1 and layer 2 can be written, respectively, as 61(532) = (1 + ”2):; 2'“ )1(A1 + Bide—“5 + (Cl + Dlz)e251i(6) we I 02(62) = (”’2’“ g2’2)1(A2 + 322)e-z€ +(C2 + 022)” ]17(5) (3-7b) 524‘ For layer 2, the stress at a point far from the origin vanishes. Therefore, 02 and 02 shall be zero and function GZ(§,z) becomes (1+v2)(l—2v2) 2 (A2 + 1322):”; 17(6) (3-7c) 02(622) = 529‘ 43 With the use of Equations (2-18) and Hankel’s inverse transformation, the boundary conditions, i.e. Equations (36), become at z = 0 0:“) =1 (3-8a) 6;“) = o (3-8b) at z = h 0:0) = (5(2) (MC) 6332‘” = 5,52) (3—8d) 14:“) = 1):”) (3-8e) 1):“) z u:(2) (Hf) Substituting Equations (2-17) into the foregoing equations and defining 02 =0 and 02:0, we have 1+Bl(1-2v1)+01(1—2v1)+Alg—Clg=0 (3-9a) —281v1+2D1v1+ A1§+C1§ = 0 (3-9b) 44 Ale-”15 6- Aze‘hl'56 — Clehléé- 3167’": (—1 + 2v] 4116) (3-9c) —32e‘h15(1—2v2 +h16)— Dleh15(—1+2v1+h16)= 0 Ale‘hl4‘6- Aye—”156 + Clehléé + Bze—hlg (2v1 — r216) (3-9d) — ale-”1592122 + 111;) + ole-”15am + 1116) = 0 _ Ale—”150 + v1) + Clehl‘fa + v1) + Age—h1§(1+v2)+ Ble‘h1‘5(1+v1)(1—h16) El 51 E2 E19" + Bze’h15(1+ v2)(-1 +h16) + Dleh1‘5(1+ v1)(-l + 1116) = 0 E25 515 (3-9d) Ale‘h1§(1+ v1) + Cleh1§(l+ v1) _ Aze‘hl‘5(1 + v2) E1 El E2 + 325’": (1 + v2)(-2+4v2 -h16) + 8184115 (1 + v1)(2—4V1 + 1116‘) (3-9e) E16 E15 + Dieh1§(1+V1)(—2+4V1 + 1116) = 0 E16 Solving the six equations simultaneously, we can identify A1, 8,, C1, D1, A2 and 82. With the use of Mathematica ®, the computational work can be significantly reduced since the formulae can be manipulated for the convenience of programming and the potential errors can be reduced to minimum. Details of the computational scheme can be found in Appendix D. 45 3.3 Two-Layer Composite on Rigid Base The formulation of a two-layer composite situated on a rigid half base will allow us to evaluate the effect of interfacial bonding condition on composite performance. In this investigation, we assume there is no friction on the contact surface and there are perfect bonding conditions between layer 1 and layer 2, and between layer 2 and the rigid half base. p(r) H. r Layer 1 26"—> m Layer2 IE2 ///////////////// Rigid Base 172 Figure 6 Two-layer composite on rigid base. The boundary conditions for this problem are as follows at z = O 0 £1) = P0) (3-10a) IS) = 0 (3-10b) 46 atz=h1 0'9) = of,” (3-100) 6);) = z)? ‘ (3-10d) u,(.1) = u,(.2) (3-10e) 119) = u?) (3-10f) at z = h] + h2 14,12) = 0 (3-109) 115,2) = o (3-10h) The function of G(.§,z) for layer 1 and 2 can be written, respectively, as G1 (6.2) = (1+ ”2):; 2’“ ’ [(A1 + Ewe—if + (C1 + D1z)ez‘v‘ 117(6) (Ma) 1 02(62) = (1 + ”2)" 7121qu + 3220f“; + (C2 + Dzz)e2"='117(6)(3-11b) 526 Combining Equations (2-18) and (3-10), and using Hankel’s inverse transformation, we have atz=0 47 0:“) :1 (3-123) 7:)” = o (3-12b) at z = h] 0;“) = e29) (3-12c) 1;“) = 632‘” (3—12d) 11:“) = 15(2) (3-12e) u:(l) = 15(2) (3-12f) at z = h] + h2 11:”) = o (3129) 11;”) = o (3-12h) Substituting Equations (2-17) into the above boundary conditions, we Obtain the following eight equations l+Bl(l-2V1)+D1(1—2v1)+Alé—Clg=0 (3-133) 231v1+201v1+ A1§+C1§ = 0 (3-13D) 48 Ale‘h156 - Aze“”1§6 -— Cle“"1‘5 6 + Czehl‘56 + Bze‘h15(—1 + 2v2 — I216) + 31947150 - 2v1+h1§)- Dleh1§(—1 + 2v] + hlg‘) (3-13c) +Dzeh1‘5(-1+2v2 +h14‘) = 0 Ale‘hlé‘g- Aze'hl‘v‘g + Cle'h15 6 —C2e”156 + Bze_h1§(2v2 —h16) (3-13d) + ale-”15(-2vl + 1116) + Dlehlg (21)] + i116) + Dzeh15(—2v2 — r1116): 0 _ Ale—hl‘ffl +111) + Clehl‘fa + v1) + Aye—”150 + v2) 51 El 52 __ Czeh15(l+ v2) + Ble‘h1§(1+ v1)(l 4116) + Bze'h15(1+ v2)(-1 + 1216) 52 515 526 + Dleh1§(l+v1)(l + I216) _ Dzeh1§(1+ v2)(1 +1116) z 1514‘ 526 (3-13e) 0 Ale—”150 +121) + Clehléa +v1) _ Aze’h15(1+v2) __ Czehléfl +122) El El 52 52 + Bze’hl‘f (1 + v1 )(—2 + 4121 — 1116) + ale-”15 (1 + v1)(2 — 4v] + 1116) 525 526 + Dlehlgfl +v1)(—2 +4121 + 1216) _ Dzeh15(1+ v2)(-2 +4v2 + I116) = 515 526 (3-13f) O 49 _ Aze-(hl ”'2 ”(1 + v2) _ 3267“" “’2 )5 (I + V2 )(2 - 4V2 + (116 + I726) E E 2 25 (3-139) _ C2e(hl ”2); (1 + v2) _ 0262“” ”2’5 (1 + v2)(—2+4v2 + h]: + I126) _ 0 E2 526 _ Aze‘IhI +h2)‘5(1 + v2) + 326“" “’2 “5(1 +122 )(1 4116 - (226) E 52 2‘5 (3-13h) + C2e(h1+h2)5(1+ v2) + Dze(”1+"2)5(1 + v2)(1 +1116 H126) _ 0 E2 525 Solving the eight equations simultaneously, we can obtain solutions for A1, B1, C1, D1, A2, B2, C2 and D2. In order to simplify the mathematical work, it is an advantage to use Mathematica®. Details of the procedures are given in Appendix E. 3.4 Validations To the author’s best knowledge, there is no analytical solution for the stress and displacement components of the contact problems except the case investigated in Chapter Two, i.e. a half elastic space loaded by a point force. Hence, numerical solutions are sought to evaluate the stress and displacement components for the three contact problems mentioned in this chapter. In this research, Gaussian integration method is used for the stress and displacement integrations. With the use of Mathematica to simplify the formulations, C++ programs are developed to evaluate the integrations. The following validations are then carried out: 50 1. Validate the numerical solution for the point-loaded half space problem with comparison to the analytical solution. 2. Validate the numerical solution for the problem of an elastic layer on a rigid half base with comparison to the analytical solution for the point- loaded half space problem. 3. Validate the numerical solution for the problem of an elastic layer on a rigid half base with comparison to the analytical solution for the point- loaded half space problem by largely increasing the thickness of the elastic layer. 4. Validate the numerical solution for the problem of the two-layer composite on a rigid half base with comparison to that of an elastic layer on a rigid half base by largely increasing the thickness of layer 2. 5. Validate the numerical solution for the problem of the two-layer composite on a rigid half base with comparison to that of an elastic layer on a rigid half base by using identical material properties for the two layers. 3.4.1 Numerical integration vs. theoretical solution In Chapter Two, we derive the theoretical solution for the elastic half space under point load. The range of the integration is from 0 to 0° for the parameter 5. However, this kind of integration is not possible in the numerical process. Gaussian integration is thus used to evaluate the integration and to determine the integration upper limit of parameter 9‘ of Equations (2-25). 51 f(§) 64 points 64 points 64 points A A [_A_\ Y \ ...... 1 1 I 1 g 0 a1 32 ...... a19 a20 6 Figure 7 Gaussian point integration algorithm. The numerical integration uses 64 Gaussian points. The integration interval is divided into many small sections, typically 20 sections as shown in Figure 7, to achieve better convergence and accuracy. The parameter I; is normally chosen to be 40 for a unit contact radius. 52 ’’’’’’’’’ ' __o- -(').04 AW 1 -0.06 ——z=0 exact ii“ -0.08 —----z=1 exact g -------- z=2 exact i3 '0'10 --—-z=3exact 1: z=0 numerical o z=1 numerical) A z=2 numerical i _o z=3 numerical i 0.0 1.0 2.0 3.0 4.0 5.0 6.0 rla Figure 8 Numerical solution vs. exact solution for elastic half space. As shown in Figure 8, the numerical integration and exact solution matches very well. The Gaussian integration points and the upper limit of g for the integrations are therefore used subsequently to evaluate more complex contact problems in this research. 3.4.2 Validations among cases Solutions for the following contact cases are obtained: (1) a half space, (2) an elastic layer on a rigid half base, (3) an elastic layer on an elastic half base and (4) a two-layer composite on a rigid half base. Among these solutions, some cases should share the same solutions if the thicknesses and material properties of layers are well chosen. For example, case (2) should become case (1) if the 53 thickness of the elastic layer in case (2) approaches infinity. Similarly, case (3) should reduce to case (1 ), if the elastic layer and the base of case (3) has the same material properties as the half space of case (1 ). Similar comparisons can be applied to cases (3) and (4) and others. Only the aforementioned comparisons will be presented in this research. (i) Comparisons between case (1) and case (2) If we increase the thickness of the elastic layer of case (2) to infinity, we should arrive at the same result of case (1 ). For a spherical indentation, if the layer thickness is 30 mm, Young’s modulus is 7GPa, Poisson’s ratio is 0.33 and the contact radius is 1 mm, the comparisons of stress and displacement components between case (1) and case (2) are shown in Figure 9 to Figure 12. Case (2) reduce to case (1) very well. 54 r—B- 2:0 ha—Ifispac; 40") ' +z=1 half space 1 i E 40.0 +z=2 half space 1 E. —e—z=3 half space . l 6 430.0 ‘ —-)I(—z=4 half space 1 i +z=0 one layer on rigid 40-0 ,. —°—z=1 one layer on rigid I +z=2 one layer on rigid j '50‘0 ‘ —0—z=3 one layer on rigid 1 . . -60.0 , f—z=4 one layer on rrgid i y 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r (mm) 1.0 ~ ~ ] 0.0 '1‘0 ‘ -B—z=0 half space 71 1 A -2.0 +z=1 half space ‘ g +z=2 half space ' g -3 0 - -e—z=3 half space ’ g +z=4 half space i '4-0 +z=0 one layer on rigid 1 —0—z=1 one layer on rigid I -5.0 1 . . ‘ +z=2 one layer on ngrd I -6.0 -0-z=3 one layer on rigid i ___:z=4 one layer on rigid I -7.0 . . 1 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r (mm) Figure 10 Comparison of shear stress between case 1 and case 2. 55 —B— z=0 one layer —e—z=1 one layer ( +z=2 one layer 1- —e—z=3 one layer " +z=4 one layer ....... i +z=0 one layer on rigid I .7 . - 1.2E-02 1.0E-02 8.0E-03 Uz (mm) 4.0E-O3 - 2.0E-03 1 0.0E+OO - 1 0.0 1.0 2.0 Figure 12 Comparison of normal displacements between case (1) and case (2). 6.0E-03 A ////// nnnnnnn nnnnn r‘ . n"- "r‘ a--- "n fihHHS-Q--.'.c "u -- :::::::-“’c..... ~-~ ....,,,.: .... -.. .. _ -a—z=0 one layer 1 1 [ —e—z=1 one layer I i ‘ +z=2 one layer 1 1 I -e—z=3 one layer 1 1 1 +z=4 one layer . i —-—z=0 one layer on rigid ; l —o—z=1 one layer on rigid { ' -—z=2 one layer on rigid { I —o—z=3 one layer on rigid 1| w :-z=4-91)e layer on rieig -‘. 56 -, ‘ -o—z=1 one layer on rigid . 1' —t—z=2 one layer on rigid i " 1 —0—z=3 one layer on rigid 3 :1 A_,_E=4£D? ILL/er QMIQE .' i (ii) Comparisons between case (1) and case (3) The case consisting of an elastic layer situated on an elastic half base, i.e. case (3), can be reduced to the half space case, i.e. case (1 ). The numerical results are given in Figure 13 to Figure 16 based on the foregoing material properties and thickness. 1 —13—z=Oonelayer ‘1 —e—z=1 one layer I I +z=20nelayer I 1 —e—z=3onelayer ‘ _ +z=40nelayer —o—z=0twolayer I I —o—z=1twolayer I; —-—z=2twolayer I —-z=3twolayer I —z=4twolayer 1 I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I 1 r (mm) Figure 13 Comparison of normal stresses between case (1) and case (3). 57 ' —e—z=1onelayer +z=20nelayer 1 —e—z=3onelayer +z=40nelayer ! —-—z=01wolayer 1 1 —+—z=1twolayer I 1 —~—z=2twolayer I, ' +z=3Molayer I“ ( I —z=4twolayer J1 l 0.0 1.0 2.0 3.0 4.0 5.0 6.0 Figure 14 Comparison of shear stresses between case (1) and case (3). 8.0E-04 , 1 mm 1.12.11 1 —e—z=1 one layer ‘I 4.05414 1 I + z=2 one layer 1 —e—z=3 one layer - "-"l' -, ,_ 1 +z=4 one layer 1’ ’5‘ 0.0E+00 I —o—z=0 two layer I g 1" 2:, -o—z=1twolayer|I a. ;. (151‘? if. h '—.—z=21wo|ayer 'I 3 H 3'5 " ' 4015414 N,» 1 _._z=3 two layer I, ,_ ' 1.-,-z=,4;!1veleyer I -8.0E-04 « U. 1 ‘.‘ 2'1 I -1.2E-03 1 1 1 1 I 0.0 1.0 2.0 3.0 4.0 5.0 Um) I _ I Figure 15 Comparison of radial displacements between case (1) and case (3). 58 —e:z_=06ne‘laye? 1 1 —e—z=1 one layer ’1 I1 +z=20nelayer1| l I —e—z=30nelayer I +z=40nelayer I —-—z=0twolayer I —o—z=1twolayer 1 1 : —-—z=2twolayer > I —o—z=3twolayer i 4.015433 . I If—z=4twolayer 1 , ‘ 1 I 2.0E-03 1 1 1 1 1 1 0.0E+00 1 1 1 A 1' 0.0 1.0 2.0 4.0 5 o ' Figure 16 Comparison of normal displacements between case (1) and case (3). (iii) Comparisons between case (2) and case (4) If the two layers of case (4) have identical material properties, it will hold the same results as case (2). Details of the results are shown in Figure 17 to Figure 20. 59 —a—z=0 one layer I, - '10'0 +z=1 one layer 1 1 I 13400 ‘ +z=2 one layer ‘ 1 5 -e—z=3 one layer 1 ‘ 1 :300 +z=4 one layer I ‘ ‘ I b ' —-—z=0two layer I I -40.0 -o—z=1 two layer 1 1 —-—FZIwolayer ) '50'0 * —o—z=3twolayer 1 l ‘ -—z=4twola r I «60.0 . 1 2-- W - ye” , 0.0 1.0 2.0 3.0 4.0 5.0 6.0 1 r (mm) 1.0 1 1 0 0 _ 1'3: :2: :13: 3:: 2:: 3:: 2:: :2: 22:: 22:77:.” ‘" " 1 —1.0 1 7; -2.0 - ‘ —a— z=0 one layer I I g . —e— z=1 one layer l I 4,-3.0 1 —a—z=2 one layer I E —e— z=3 one layer 1 b .4_o + z=4 one layer 1 1 1 --— z=0 two layer ‘ -5.0 —o—z=1 two layer 1 1 + z=2 two layer ' —6.0 —¢— z=3 two layer I 1 -—z=4 two layer 1 -7.0 . T i , 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r (mm) Figure 18 Comparison of shear stresses between case (2) and case (4). 60 -a— z=0 one layer I —o—z=1 one layer I —A— z=2 one layer I f —e— z=3 one layer ‘ +z=4 one layer I : —-— z=0 two layer I - —o—z=1 two layer —~— z=2 two layer _._ z=3 two layer I I —z=4 two layer I I .0 3.0 4.0 5.0 . r (mm) I, I 0.0 1.0 2 I 1.0E-02 , ,, a V, a,” tea-:- I I —e—z=0 one layer I 7 ' —o—z=1onelayer II I +z=20nelayer II I —e—z=3 one layer II I +z=4 one layer I I I I 6.0E-03 I —o—Z=0 two layer I I I? I —~—z=1twolayer II I I —-—z=2two layer I IE. 40E-031 I —o—z=3twolayer I I s l -— z=4 two @ng .I I ZOE-03 l I I W _. _. A_ I I I -2.0E-03 I I I 0.0 1.0 2.0 3.0 4.0 5.0 ‘ I r(mm) I Figure 20 Comparison of radial displacements between case (2) and case (4). 61 In the foregoing comparisons, we have simplified the complex cases to simple cases. We have validated the numerical solutions with the analytical solution for the half space problem, the problem of an elastic layer on a rigid half base with the problem of half space, the problem of an elastic layer on an elastic half base with the problem of half space, and the problem of two-layer composite on a rigid half base with the problem of an elastic layer on a rigid half base. These calculations serve as cross-validations for the derivation and programming of the numerical scheme. With all these exercises, we have successfully developed and validated the numerical scheme. In the next chapter, we will use the numerical scheme to study more applications. 62 CHAPTER FOUR BASIC APPLICATIONS We have validated the solution formulae and numerical scheme for studying contact problems in the previous chapter. In this chapter, we use the numerical scheme to study some basic contact problems, such as effects due to material property, thickness, boundary condition, loading condition and lamination of layers. Material property has strong effect on both stress and displacement distributions. We study it first. Finite thickness, especially thin layer, can alter the stress and displacement distributions in a half space. It is investigated secondly. Base condition and loading condition are another two fundamental elements for mechanics analysis. They are also critical to the response of materials subjected to contact forces. They are examined subsequently. The last part of this chapter is focused on the lamination effect of composite layers. We confine our studies to composites with perfect bonding interfaces. Only the lamination with different material properties and thicknesses are of interest. 4.1 Effect of Material Property Three different materials, such as steel, aluminum and nylon, are used in the following investigations. They represent hard, intermediate hard and soft materials. For simplicity, we choose to investigate the half space problem loaded by a point force. The material, loading condition and boundary condition are summarized in Table 2. We use relatively small force and contact radius in the studies to avoid the possible numerical overflow in the C++ program. These values can be scaled by 100 times to 1,000,000 times without changing the 63 numerical values during calculations. For the rest of the chapters, we continue to use Newton (N) for force and millimeter (mm) for length. Table 2 Material properties for steel, aluminum and nylon. Steel Aluminum Nylon Total Load P 100 N 100 N 100 N lndentor Shape near-point load near-point load near-point load Young’s Modulus 200,000 MPa 7,000 MPa 2,500 MPa Poisson Ratio 0.29 0.33 0.40 Contact Radius 1.0 mm 1.0 mm 1.0 mm Shown in Figures 21 and 22 are the elastic stress distributions. As expected, they are identical for all different materials. These results can also be attributed to the equilibrium equations, Equations (2-2), which are not dependent on any material property. The displacement distributions, however, are dependent on the material property. They are higher for softer materials, which have lower Young’s modulus, as illustrated in Figures 23 and 24. 64 I 0.0 I -2.0 I.L—9_:z=—0STIT _'I I I +z=1STL II I '40 +Z=2 STL I I -e-z=3$TL - I I -6.0 - I +z=4 STI. I A +z=0 AL . I I I g -8.0 I —o—z=1 AL I I I ‘u' 100 I +z=2 AL I I I b - . I —o—z=3 AL . +z=4 AL I 42.0 I _....z=0 Nylon I I {14.0 I I —o—Z=1 Nylon I | I ‘ +Z=2 Nylon I I I -16.0 .I I _._z=3Nylon I I ' —l—Z=4 Nylon I ; -18.0 i flPM ' I I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 -B— z=0 STL I +z=1 STL + z=2 STL I —9— z=3 STL I + z=4 STL —a— z=0 AL I —o— z=1 AL I —-¢— z=2 AL I —o— z=3 AL I + z=4 AL I I —0— z=0 Nylon I —+—Z=1 NYIOD I + z=2 NYIOl'l —o—z=3 Nylon I I :fl'fiflj I 5.0 6.0 I I Figure 22 Shear stress distributions based on different materials. 65 -5.0E-04 - Ur (mm) I 10504 a -9.0E-04 ' I -1.1E-03 It #Ifii I fly I I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I 125-02 -- . I l I 1.05-02 a 8.0E-03 6.0E-03 I Uz (mm) 4.0E-03 , I 2.05-03 I Figure 24 Normal displacement distributions based on different materials. 66 4.2 Effect of Thickness In this section, we study the thickness effect on the stress and displacement distributions in contact problems. An elastic layer situated on an elastic half base and subjected to a spherical indenter is selected for investigation. The elastic layer is made of aluminum while the elastic half base is made of nylon. The thickness of the top layer, aluminum, varies from 1 mm to 8 mm. Table 3 summarizes the contact problems. Table 3 Material properties used for thickness study Layer 1: Aluminum Layer 2: Nylon Total Load 100 N 100 N lndenter Shape spherical spherical Young’s Modulus 7.000 MPa 2,500 MPa Poisson Ratio 0.33 0.40 Thickness 1, 2, 4 and 8 mm infinity Contact Radius 1.0 mm 1.0 mm Using the analytical formulae and numerical scheme developed in Chapter Three, we calculated the stress and displacement distributions and present them in the following diagrams. 67 0'2 (MPaI +z=0(h=2)I I I i +Z=0(h=4)II I -5o.o I , .+§££h:8_)II I -6o.oI L -2- A, H». _ #I__fi _____ I I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I I +z=0.25(h=1)I +z=o.25(h=2)I I—o—z=o.25(h=4)I ; I L+F£§ P=o.7936*u*52*[—‘ a Combining equations 5.3 and 5.4, we have E —0.4529 h P=1.0788*U*E * —2 —‘- (5-5) 2 E —0.1329 . J 5.4 Two layer composite on rigid base We use a similar approach to study the case of two-layer composite situated on a rigid base. With the numerical scheme developed earlier, we can identify the effect of the layer materials, if no simple explicit analytical relationship can be idenfified. 110 1200 * 1000 y = 133105x + 0.0002 R2=1 800 2 § 600 O I.L 400 200 0 fig 0 0.002 0.004 0.006 0.008 Deflection (mm) Figure 87 Force vs. deflection for two-layer composite on rigid base. 111 y = 1 .6383x'0'854 R2 = 0.9951 PIUIEZ 0.0 5.0 10.0 15.0 20.0 Figure 88 Relationship between stiffness and the ratio of Young’s moduli. From the calculated data points, a contact model is presented as follows: —0.854 51] (5-6) P=l.6383*U*E2 *[ 51 Similarly, the stiffness and the thickness ratio are calculated and plotted in Figure 89. The relationship is P=0.9423*U*Ez *[ifl— —0.3962 2 ] Combining Equations 5.6 and 5.6, we have 112 E —0.427 It —0.l83l P=l.2425*U*E2*[—2—) [A] (5-8) Stiffness (PIUIEZ) 51 112 2.5 2.0 , 1.5 y = 0.9423X-08962 R2 = 0.9705 1.0 0.5 0.0 1 1 0.0 1.0 2.0 3.0 4.0 5.0 1 h1/h2 : Figure 89 Stiffness vs. thickness ratio of the two layers. 113 CHAPTER SIX ADVANCED STUDY In this chapter we will apply the techniques developed in previous chapters to study cases involving contact interface and bonding interface. First of all, we will derive a general formulation for cases with shear loading condition. Combining normal and shear loading conditions, we will be able to study frictional effect at indentation surfaces. We will also incorporate a shear slip condition into interfaces of layered materials to simulate cracks in the layers and to investigate the effect of the interfacial bonding on the performance of the layered materials. 6.1 Shear Loading All the equations derived for normal loading in previous chapters can be extended to shear loading. Here, the two-layer half space case is presented to demonstrate the extension. it then is used to study the friction effect at the contact surface. 114 9(r) Layer2 Figure 90 Two-layer half space with shear loading. For simplicity, there is no normal loading at the contact surface and the two layers are perfectly bonded. The thickness of layer 1 is h and the contact radius is a. If we use superscripts (1) and (2) to represent layer 1 and layer 2, respectively, the boundary conditions for this contact problem can be summarized as follows: 2 = 0 020) = 0 (6-1a) ti? = q <6-1b) z = h 0;” = a?) (649 115 r2) = r2” (6—1d) ug) = 11:2) (6-1e) ugh =ug2> (6-1f) The function of G(§,z) for layer 1 and layer 2 can be written, respectively, as 0105.2) = (1+ ”EX; 2‘ 2‘“ ) [(A1 + 812)!“ + (C1 + 012x229” 1cm) (6-2a) l Gztm = (1* ”2)“ ' ”Nu/12 + 3206—25 + (62 + 022W 151:) (6-2b) 5252 For layer 2, the stress at a point far away from the origin vanishes. Therefore, C2 and Dz shall be zero and function 02022) becomes 021:,z>= (”VZXIQZVZ’W + hoe-Zing) , (6-2c) 526 With the use of Equations (2-18) and Hankel’s inverse transformation, the boundary conditions, i.e. Equations (6-1), become 0:“) = 0 (6-3a) .116 7:9) :1 (6-3b) z = h 0:0) : 0:;(2) (MC) 4,“) = 62‘” (6-3d) 21:“) = 15(2) (6-3e) u:(1) : u:(2) (5-30 Similarly, if we use Equation (2-19) and replace 13(5) with (7(6) in the stress function of 605,2), we then have the following equations: 31(1—2v,)+01(1—2v1)+Alg—Clg=0 (6-4a) -]-231vl+201v1 +A1§+Cl§=0 (64D) A1§*A2§*C162h1§§-31(-1+2V1‘h1§)-32(1-2V2 +’714‘) (6-4c) —Die2”15(—1+2vl +1115) = 0 _ 73/116 _ _ _ A16 Azé-Cle 3: 32(2vz h1§>+31(2V1+/71§) (6-4d) - 01e2h15(2v1 + 121:) = 0 117 _ A1(l+v1) + C1(1+vl)e2”lcf + A2(l+v2) + Bl(l+v1)(l-h1§) El 51 52 E15 (6-5e) + 32(1 + V2)(‘1 + 1114:) + 010 + V] )(1 +hl‘kahlé: ___ 0 E26 E15 Al(1+vl) + C1(1+V1)€2h1§ _ A2(1+V2) + 32(1+Vl)(‘2+4"2 -h11§) E1 E1 E2 E1; (6'50 + Bl(l+v1)(2-4v1 44114") + 01(1 +V1)("‘2+4’V1 +hlg)82hlé : 0 E16 £15 From Equations (6-5), we can determine A1, B1, C1, D], A2, and 32. Details of the procedures are given in Appendix F. Based on the above equations, we can evaluate the internal stress in contact conditions involving shear loading. 6.2 Contact Problems with Frictions In the previous section, we developed a technique to investigate contact problems involving shear loading. There is a relation between normal and shear loading in Coulomb friction or sliding friction. It can be defined as 61(5) = f * p(é) (6-6) where p(g) and q(§) are normal and shear stresses, respectively, and f is the friction coefficient on the contact surface. In most cases, we consider f a constant between the two materials contacting each other. We can always separate the contact force into normal and shear components. The two force components can be identified independently as shown in Chapter Four and section 1 of this chapter. Based on linear elasticity, the results from the normal 118 and the shear loadings can be superimposed to represent the total solution. A C++ program is written to evaluate the case mentioned above. When f = 0, it represents a smooth contact interface. The cases of f equals 0, 0.05, 0.10 and 0.20 and 0.40 are examined. The results contributed by the shear force only are presented in the following diagrams. The material properties are listed in table 4. Table 4 Material properties for friction study lndenter Shape Young’s Poisson . and Total Force Modulus Ratio Thickness Matef'a' 1 Spherical 7.0E+4 MP3 0.33 1 mm (Aluminum) P=100 N Material 2 - 43,68.) (compreSSIon) 2.0E+5 MPa 0.29 Half space 0.5 0.0 1 .......................... -O.5 E -1.0 - E4 5 g .. 6 -2.0 I -2.5 l -3.0 ‘ +'zé17E0.05)l -3,5 ‘-&—Z=1 (f=0.1) I +z=1 (f=0.2L': -4.0 ’i‘ h 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r(mm) Figure 91 Normal stresses at bonding interface between layer 1 and 2 119 * +£— z=1 (f=0.05)l *—a—z=1 (f=0.1) I l 7*.Z_=1_ “70:21-1 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I r (mm) .-rn” _________________________ 'pi ,. v'. " lie—3:1 Riki? fat—z=1 (f=0.1) 1 7.6633368); 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r (mm) Figure 93 Radial stresses at the bonding interface between layer 1 and layer 2 120 I... . . 0.0E+OO '7, .szf=:::‘ L -1.0E-06 ’ .A. ‘cacc' . ‘:::.... I -2.0E-06 - x , 114‘“ I h 3-3.0E—06 n‘.issé;;;§3g.}éé5*‘fiiéiafl I l .5. I 5 54.0506 « I } l -5.0E-06 I M . I +z=1(f=0.05)l I -6.0E-06 . +z=1(f=0.1) I I ‘ I ‘—x—z=1(f=0.2) l l JOE—06 , . , . . 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I r(mm) l I I—s—2¥1(1=0.05)l l ' '—-e—z=1(f=0.1)I I Ian—z=1 (1:02) 1 1.0E-05 1 ‘ E:: | l . 5.0E-06 - . I ‘ “::=::; ..... 1' ~-- . I I I 0.0E+00 . ; i 0.0 1.0 2.0 3.0 4.0 5.0 6.0 1 r(mm) 1 Figure 95 Displacements at bonding interface of layer 1 and 2 121 From the results shown in the foregoing diagrams, the stresses and displacements both increase as the friction coefficient increases. 6.3 Bonding Theory Liu et al (1994) used a bonding theory to analyze imperfect composite laminates. The theory was based on a relationship between the mismatch of displacements on the opposite sides of an interface and the shear stress on the laminate interface. The theory is defined as follows: ul" w?) = krarz (6-7a) 149) — u 12> = kzaz (6-7b) where u?) represent the displacement along j direction in layer i, and k, and k2 represent the shear slip coefficient along rdirection and the normal separation coefficient alongz direction, respectively, at the interface. When k, and k2 approach zero, it presents a perfect bonding condition. The interface can consist of an imperfect bonding zone and a perfect bonding zone. A distance d from r=0 is defined to represent the imperfect bonding zone. 122 p(r) H I011 l . Layer1L—28—* h r l Layer2 Figure 96 Shear slip theory of two layer half space In order the solve this problem, the boundary condition in Equations (3-6) are modified as follows 2 = 0 09) = p(r) (6'83) .19 = 0 (6-8b) z==h 09) = a?) (6-8c) 4;) = ,g.) (6-8d) 14$” —u§2) =krorz <6-8e) 123 u?) — u?) = kz oz (6-8f) 6.3.1 Entire Imperfect Interface When d is infinite, k, and k: are finite, the boundary conditions, i.e. Equations (6- 8), can be expressed as follows with the use of Equations (2-18) and Hankel’s inverse transformation: 2 = 0 6;“) =1 (6-9a) 1,25” = o (6-9b) z = h 6:“) = 62(2) (6-9c) r22") = r13” (6-9d) fl” -1132) =60; <6-9e) 11;“) w?” =kzcr§ (6-90 When we substitute Equations (2-19) into the above functions, we have the following six equations: 1+Bl(l-2v])+D](l-2v])+Alé-C1§=O (6-108) 124 _281VI+2D1V1+ A1§+C1§ = 0 Ala: - A2; — C1e2h15§ — BI(—1 + 2v] -/11é)+ 32(1— 2122 + 111.5) — DleZhl‘fG-l + 2v] + I216): 0 A16 - A25 + Clezmé + 32(2vz 7 h1§)+ B1(-ZV1+ h1éf) + 01e2h15(2v1+h1§)= o _ A1(1+v1) + CIe2h15(1 +121) + 31(1 +v1)(l 411;) 151 51 E16 + Dle2h1‘5(1+v1)(1+h15)+A2(1+v2 4%,) E15 E2 (1+v2)(—1+h1:)) = 0 529 + 32 (21(er - hlkré: + A](1+v1) + CIe2h15(1+vI) + B](1+v1)(2—4v1 +1115) El El E15 + 01e2h1‘5(1+v1)(—2+4v1 +1115) +A2(-l+v2 4:26) 1:315 E2 + 320-162 + 2km — 121ng + (1+ V2)” + 4V2 415)) = 0 525 (6-10b) (6-10c) (6-10d) (6-10e) (6-10f) From Equations (6-10), the problem is solved with help of Mathematica®. A C++ program is also written to obtain the numerical results. 125 6.3.1.1 Verifications First of all, if both k, and k2 equal zero, we have the same result as perfect bonding condition as shown in the following diagrams. The material properties are listed in table 5. Table 5 Material properties for friction stud)! lndenter Shape Young’s Poisson Th' k and Total Force Modulus Ratio '6 ness Material 1 . (aluminum) gamete]! 7.0E+4 MPa 0.30 1 mm Material 2 . (nylon) (compreSSIon) 2.5E+3 MPa 0.29 Half space 0.0 ' — — —— _ __ _ w .z _ -1.0 . I -2.0 I E -3.0 I E . 6 -4.0 -5.0 ,-__- ___ ,_ 7 7 +z=1 (kr=0, kz=0) -6.0 . . " - —9—z=1 (perfect bonding) .I -7.0 1 - - 7. - ~ ’7 " ;. f 7 :z ‘i;:"i;t. 0.0 2.0 4.0 6.0 r(mm) Figure 97 Normal stress verification 126 Ur (M Pa) -0.2 ~ -0.4 1' -0.6 j ‘ E 1 1: -1.0 I b I -1.2 1 I '1-4 "1 —o—z=1 (kr=0, kz=0) I I l l '1-5 ‘ —e—z=1 (perfect bonding) 3 . -1.8 as I , ---- if” - ’ ‘f—VI 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r(mm) Figure 98 Shear stress verification 0.0 - __ 0.0 'I — -— —~ I I _ +2=1(kr=0,kz=0) I ‘ 0.0 1 I ~ —e—z=1 (perfect bonding) I 0.0 VI--- . -. __ 7. : 77:77:- I 7:1 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r (mm) Figure 99 Shear displacement verification 127 Uz (M Pa) 0.0 I 0.0 1 0.0 i 0.0 1.0 2.0 + F1 (kl=0, kz=0) —e—z=1 (perfect bonding) i .7, f’ 3.0 4.0 r (mm) 5.0 1‘. I. 6.0 Figure 100 Normal displacement verification 6.3.1.2 Variation of shear slip coefficient Kr If we let kz equal zero, we can study the stress and displacement as a function of kr. Table 6 lists the material properties used in the calculations. Table 6 Material properties for friction study lndenter Shape Young’s Poisson . and Total Force Modulus Ratio Th'c'mess Material 1 . - Spherlcal 7.0E+4 MPa 0.30 1 mm (11:01.1? (nylon) (compressron) 2.5E+3 MPa 0.29 Half space From the following results, we can see that when kr is very small (less than 0.00001), the stress and displacement is close to the perfect bonding. Figure 101 and figure 102 show that kr can affect the shear stress more significantly than the 128 normal stress. Similarly, Figure 103 and 104 show that kr has more effect on the shear displacement than the normal displacement. For kr between 0.0001 and 0.01, some sinusoidal oscillations can be found in the numerical results. These values of kr represent changes of interface from perfect bonding to complete debonding. The sinusoidal oscillations is due to the Bessel function in the numerical evaluation. When the integral range for 5 increases, the magnitude of the vibration decrease and approaches stable as shown in Figure 102. Since the Bessel function converges very slowly, to achieve a perfect numerical result would not be realistic. I 0.0 I -1.0 ~ I A-2.0 1 ’ 7 M II I . m —a—z=1(kr=0) I 1 $73“ +z=1 (kr=0.00001) ‘. I ~40 1 +z=1(kr=0.00005) I ‘ j °_50 , +z=1 (k=0.0001) . I I ' —e-z=1 (kr=0.0002) I 76-0 . —o—z=1 (kr=0.001) -10 ,, +F1 (kr=0.01) I I +z=1 (kr=0.1) I ' ‘8.0 l 1 1 I 0.0 1.0 2.0 3.0 4.0 5.0 6.0? r(mm) Figure 101 Normal stresses vs. shear slip coefficient kr 129 0.0 I .10 I320 l E30 " —a—‘ z=1 kFO I l E +z=1 kr=0.00001 I , I 04.0 ~ +z=1 kr=0.00005 . l —e—z=1 kr=0.0001 . -5-07 , _ —e—z=1 kr=0.0002 I l -- +z=1 kr=0.001) ‘ w I 45-01 ~ —o—z=1 kr=0.01) I l ‘ 70 —)|(-Z=1 kr=9.1) . . l 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I I r (mm) I —B— z=1 (kr=0) 1-05-03 7 +z=1(kr=0.00001) I i +z=1 (kr=0.00005) I 8.0E-04 1 ...... —e— z=1 (kr=0.0002) I ' ~ I“ + z=1 (kr=0.001) 6.0E-04 ,7 +z=1 (kr=0.01) I 4'5"???" .....,,_II‘_II-II‘IIIIII ‘¢""~‘.::.:"-'.i 2:1 (”301) l I I7 . l : E 4.0E'04 fl 1:”, “Hi, I E, 4 w 5 205-04I I I l I 0.0E+00 I I -2.05-04 -I I l ' 4.05-04 5— - -7- -_ _ .7 7. . __ I I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I l I r (mm) Figure 103 Radial displacement Ur vs. shear slip coefficient kr 130 1.0E-02 * FF ~ * - —~- WA ~#—~ I —e—z=1(kr=0) 905-03 5;?“ +2: :1(kr=0.00001I I I 8.0E-O3 +z=1(kr-000005 I A 7.0E-03 .......... ‘ :§:1I$888%2I I I E 6.0E-03 '- .- :z_1(g_-8.o1) I I}: 5.05-03 n.“ 2- ( - .) I I I3 4.0E-03 " . w . I;¢\::"--- I 3.0E-03I ‘ f , 2.05-03 I I 1.0E-03 I I 0.0E+00 I. , — . . W . we a I I 0.0 1.0 2.0 3.0 4.0 5.0 6.0' r (mm) , _ I Figure 104 Normal displacement Uz vs. shear slip coefficient kr 6.3.1.3 Variation of normal separation coefficient Kz We have discussed the stress and displacement variation when kz equals zero. Similarly, we will study the changes of stress and displacement with respect to kz when kr equals zero. The results are presented in the following diagrams and the material properties are listed in table 7. Table 7 Material properties for friction study lndenter Shape Young’s Poisson , and Total Force Modulus Ratio Thickness Material 1 . ialuminum) imam 7.0E+4 MPa 0.33 1 mm Material 2 . (nonn) (tenSIon) 2.553 MPa 0.29 Half space Figures 105 and 106 show that the bonding condition is near perfect when kz is very small, i.e. less than 0.00001. When kz increases, the normal stress at the 131 bonding interface decreases significantly; however, the shear stress is not as much affected by the change of kz as the normal stress. Similarly, the normal displacement at the bonding interface also decreases more significantly than the shear displacement as shown in Figure 108. , +z=1(kz=0) +z=1(kz=0.00001) I I —e—z=1(kz=0.00005) I , + z=1 (kz=0.0001) I I —e—z=1(kz=0.0002) II I I —o—Z=1 (kz=0.001) I I . .I +z=1(kz=0.01) I I +z=1(kz=0.1) IN I - ._ I I I 0.0 I I I 1' , , I 0.0 1.0 2.0 3.0 4.0 5.0 6.0 r (mm) Figure 105 Normal stresses vs. normal separation coefficient kz 132 II II +Z=1ékz=m 160 I +z=1 kz=0.00001) ' +z=1 Ikz=0.00005) I 1 40 I I —e—z=1 kz=0.0001) - r —e—z=1Ikz=0.0002) , +z=1 kz=0.001) 1 1.20 I I +2§1Iszo.01) I I @100 kz—0.1) -0 I 0' . l .. 3 £0.80 ; , u I 00.60 - I 0.40 I 7- I I 0.20 I I ~ I I 0.00 W W -—W 0.0 1.0 2.0 3.0 4.0 5.0 6.0 I r(mm) . Figure 106 Shear stress vs. normal-separation coefficient kz 0.0E+00 -4.0E-04 , , ‘. —9—z=1 kz=0) I"-8 0E-04 +Z=1 kz=0.00001 1 E - kz=0.00005 Iv ¢ o :Fl rig-33?? , h _ _ I “ """"""" z= - . '3 ”E 03 +2=1Ikz=0.01) I —e—z=1 kz=0.1) ' -1.6E-03 I I I “.'\‘.““-‘:-‘;:::3:32;ch 4.05-03 , - W I I I 0.0 1.0 2.0 3.0 4. 5.0 6.0 Figure 107 Shear displacement Urvs. normal separation coefficient kz 133 —a—z=1 Ikz=0) I | IA _ _ I +z=1 kz=0.00001) : E 2'0E 02 +z=1 kz=0.00005) I: I +FI F8883?) - I +F Z= . ‘ =3 '3-05'02 . Ikz=0.01) I i —e—z=1 kz=0.1) I 405-02 I IIIIIII I I -5.0E-02 - IIIIIIII I - — ----- 0.0 2.0 4.0 6.0 r(mm) I I_.-‘__.+ ___ Figure 108 Normal displacement Uz vs. normal separation coefficient kz 6.3.2 Imperfect bonding with finite length II IN = Layer1 L'Za —’ T r Figure 109 Two-layer composite of imperfect bonding with finite length For a finite imperfect interface with a radius d, we will study the following problem to demonstrate some fundamental phenomena. The problem involves a step 134 function. When r S d, kI.(r) = k,. and when r > d, k,(r) = 0 where d is the half crack length as shown in Figure 104 and i represents either r or z. The nonlinear shear slip coefficient and normal separation coefficient are defined as follows: 15(5) = Irk,1.0 15 Figure 110 Shear slip coefficient kr (g) distribution. 0.50 0.40 0.60 1 0.30 ~ 0.20 0.10 0.00 -0.10 Figure 111 Shear slip coefficient kz(§) distribution. kz=1.0 when r510 kz=0.0 when r>1.0 20 136 10 f. 15 We can’t see the direct relationship between the coefficients and the length of the imperfect bonding, though significant efforts were made to find proper functions for the shear slip and normal separation coefficients. We plot the coefficients against 5 and also summarize some critical values from the curves in table 8. Finding better functions for the shear slip and normal separation coefficients still remains for future study. Table 8 Some values for different imperfect bonding length d max kr* «f @743?” { @k1rs; time max kz* f @klfigme 1.0 0.2543 2.4516 5.8843 0.5 3.8317 2.0 1.0174 1.2254 2.9422 2.0 1.9159 3.0 1.9614 0.8172 2.2891 4.5 1.2772 4.0 4.0695 0.6129 1.4711 8.0 0.9579 5.0 6.3585 0.4903 1.1769 12.5 0.7663 6.0 9.1563 0.4086 1.3473 18.0 0.6386 8.0 16.2778 0.3065 0.7355 32.0 0.4790 10.0 25.4341 0.2452 0.5884 50.0 0.3832 137 . i l: l l l l . l 20 . ‘ d-8 F i - l / i i 15 i ii . l . 5 i‘- 10 " .' \ l ' ‘i l. . . ‘ "' x .‘ it u x f l . i". . / _q 9". Mr. ‘ S ‘4. ‘ '. ; A ,4. . g‘ Y . ‘, xv irr‘frx ._ _ ”u __ H . , 2,, H ' ‘ l" ’ ‘1"‘4w" 2'.‘ */J 35%.)7773‘11 x.’:‘.;..,"r;13-"'nr"\ ‘rrur ' .P‘.;.',}:Jf;\'.' x l.’ Li~ ‘ 1 ' .iflr' 7.. . l 2 3 4 5 6 _5 L <5 Kz* 50 (1:10 Figure 113 Normal separation coefficient vs. imperfect bonding length 6.3.2.1 Verifications Verifications for the technique outlined above include the following special 1) perfect bonding on the entire interface as presented earlier 138 2) any d value with both kr and kz equal to zero, i.e. perfect bonding on the entire interface 3) d equals zero with random kr and kz. i.e. perfect bonding on the entire interface The results are shown together in the following diagrams. The material properties are listed in table 9. Table 9 Material properties for friction study lndenter Shape Young’s Poisson and Total Force Modulus Ratio Thickness Matef'a' 1 Spherical 7.0E+4 MPa 0.33 1 mm (aluminum) _ M t _ | 2 P-100 N a ena (compression) 2.5E+3 MPa 0.29 Half space (nylon) 139 1 0.0 ,,,,,, i -1.0 l -2.0 i $43.0 , .2, 1 64.0 ,v. # __ 7__ 2. __ l _50 - —~—z=1.0(kr=0,kz=0) I l —o-Z=1.0(d=0) l I .60 l —e—z=1.0(perfectbonding) l l -7.0 ' , T5”. ‘ “ ; , 0.0 2.0 3.0 4.0 5.0 i r(mm) Figure 114 Normal stress for all three cases. 0.0 1 ’ -0.2 : -o.4 i i a -o.e ; . o. i E -0.8 l E l b ‘1-0* “2 .__ __ __ __ ; -12 . l —a—z=1.0(kr=0,kz=0) l 14 i +z=1.0(d=0) l . _ _ . i ‘ i . —°—z=1.0(perfectbonding) ! -1.6 5 fl .._ f -_ ~ 2222. 0.0 2.0 4.0 6.0 2 ' (mm) _ 6.0 Figure 115 Shear stress for all three cases. 140 U r (mm) 8.0E-03 U2(mm) 7.0E-04 6.0E-04 5.0E-04 4.0E-04 3.0E-04 ~ 2.0E-04 1.0E-04 * 0.0E+00 Figure 116 Displacement Ur for all three cases. ,#_L_ ,___1. 'A 7‘] —+— z=1.0 (kr=0,kz=0) + z=1.0 (d=0) 0.0 7.0E-03 6.0E-03 5.0E-03 4.0E-03 3.0E-03 2.0E-03 * 1.0E-03 - 0.0E+00 Figure 117 Displacement Uz for all three cases. .0. z=1.0 (perfect bonding) j 1.0 2.0 3.0 4.0 5.0 r (mm) 6.0 ; +z=1.0 (d=0) —o— z=1.0 (perfect bonding 7+F1.bagr;o,kpor 1 ) l 0.0 1.0 2.0 3.0 4.0 5.0 _ r (mm) 141 6.0 6.3.2.2 Variation of shear slip coefficient kr Similarly, we will investigate how kr affect the stresses and displacements when kz equals zero. Using the C++ program established earlier and incorporating kr, we have the following numerical results. The material properties are listed in‘ table 8 before. We can see that kr doesn’t poses significant effect on the normal stress and displacement, as shown in Figure 118 and 121. However, kr affects the shear stress and shear displacement significantly. Similar to section 3.1.2, there are some sinusoidal oscillations when kr is between 0.0005 and 0.01. Increasing the integral range doesn’t not eliminate the oscillations due to the slow convergence of the Bessel function in the integration. Beyond this range, the results are close to perfect bonding when kr is very small and the shear stress is very close to zero when the interface is completely debonded. 142 ,+ z=1 (kr=0) i , l ‘ a -2.0 ‘—a—z=1(kr=0.0001) ; ~ g -30 ,+z=1(kr=0.0005) 7; 4.0 —o—z=1(kr=0.0008)! j ‘ D _5.0 o~z=1(kF0.0009) , 45.0 r i+z=1(kr=0.001) l _ . +z=1(kr=0.01) , '7'0 F —o—z=1(kr=0.1) l _ so. * ' 2.224 . 0.0 5.0 10.0 15.0 1 r(mm) Figure 118 Normal stress vs. shear slip coefficient kr. l —9—z=1 (RFD) i . +z=1 (kr=0.0001) ‘ 45— z=1 (kn-10.0005) ‘ +z=1 (kr=0.0008) l e z=1 (kr=0.0009) * +z=1 060.001) + z=1 (kr=0.01) —+—z=1 (kr=0.1) l: -5.0 . F . . 0.0 2.0 4.0 6.0 8.0 10.0 12.0 __ T (m). Figure 119 Shear stress vs. shear slip coefficient kr. 143 Ur (Miss)— 8.0E-04 ——. - - ______ ._ l +¥1 (MEET *7 ' 7.05-04 - I+z=1 (kr=0.0001) . 6.0504 g ' .a—z=1(kr=0.0005) f i +z=1 (kr=0.0008) 5.05-04 - ' 49— z=1 (kr=0.0009) , - w +z=1 (kr=0.001) ? 1 4'0E 04 i +z=1 (kr=0.01) ‘ i 30504 - _, l_._z=1 (kr=0.1) J ' 2.05-04 l ' i 1.0E-04 l 0.0E+00 . r 5‘ - 0.0 2.0 4.0 6.0 8.0 10.0 12.0 i r(mm) i Figure 120 Radial displacements Ur vs shear slip coefficient kr. 9.0E-03 8.0E-03 7.0E-03 Ke— z=1 (kr=0) + z=1 (kr=0.0001) ‘ ~9— z=1 (kr=0.0005) . + z=1 (kr=0.0008) 6.0E-03 - e z=1 (kr=0.0009) 5 0503 . +z=1 (kr=0.001) ' . +z=1 (kr=0.01) 4.05-03 - g—l—z=1 (kr=0.1) ,' 3.0E-03 - l 2.05-03 4 ! 1.05-03 - """"""""" ‘ 0.0E+00 ~ . . , . 0.0 2.0 4.0 6.0 8.0 10.0 12.0 Figure 121 Normal displacement Uz vs shear slip coefficient kr. 144 6.3.2.3 Variation of normal separation coefficient kz If we set kr to zero, we can investigate the variations of stress and displacement against kz. The results are presented in the following diagrams. It seems kz has a significant effect on both normal stress and normal displacement as shown in Figures 122 and 125. The influence of kz to the shear stress and shear displacement is, however, not significant. I+z=1 (kz=0) i rla—z=1 (kz=0.0001) j A~—I—z=1(kz=0.0005) l ‘ a; 30 - i.+z=1(kz=0.0009)'! l N ,+Z-_-1 (kz=0.001) l i b 2.0 ‘ ' j+z=-1(kz=0.01) ‘ l L:o—z=1 (kz=0.1) ' .,. ‘\ l . - “5"“? 1 Figure 122 Normal stresses vs. normal separation coefficient kz. 145 i 1.4 , —9—z=1(kz=0) l +z=1(kz=0.0001) - 12 - ~~B—z=1(kz=0.0005) i ’ —o—z=1(kz=0.0008) . €10 - -—& z=1(kz=0.0009) l a +z=1(kz=0.001) ,1: i €08 l +z=1(kz=0.01) i . —+——z=1 (kz=0.1) g : D 0.6 - L‘ ----- ~— ' ‘; l 0.4 . 7 ‘ 0.2 f‘ .. . 0.0 l r l T i ‘ 0.0 2.0 4.0 6.0 8.0 10.0 12.0 ‘ r (mm) i Figure 123 Shear stresses vs. normal separation coefficient kz. , 0.0E+00 ( -2.0E-04 '5 '4-05'04 E? +z=1 (RFD) ‘ 43.05-04 g, - -' +z=1 (kz=0.0001) i 0.“? -8 0504 # K _ " ' 4F z=1 (kz=0.0005) ‘ l g ' * +F1 (kz=0.0008) 1 5 -1.0E-03 —O z=1 (kz=0.0009) I . _1. 2503 . +z=1 (kz=0.001) l . i +z=1 (kz=0.01) l ‘1-4E‘O3i --+—z=1 (kz=0.1) -1.6E—03 - ‘ i -1.8E-03 l ~+~ ~ -~ ~ -- e 2 - l 0.0 2.0 4.0 6.0 8.0 10.0 Figure 124 Shear displacements vs normal separation coefficient ikz. 146 0.0E+00 ' -5.0E-03 405-02 '1-554’2‘ *:2=1'(‘szO)W , ,+z=1 (kz=0.0001)l, TE 11 i E, -2.0E-02 , . , 5. -0 z=1 (kz=0.0005) -2.5E-02 .. l—o—z=1(kz=0.0008)l ho z=1(kz=0.0009)' '3-054’2 :+z=1(kz=0.001) i l 3.55-02 * +z=1 (kz=0.01) . 40E 02 --*-Z=1 09:01) ." 0.0 2.0 4.0 6.0 8.0 10.0 12.0 r (mm) Figure 125 Normaldisplacements vs. normal separation coefficient kz. 147 CHAPTER SEVEN SUMMARY AND CONCLUSIONS 7.1 Summary of the Research Indentation of layered composites were studied in this research. A new technique is presented to solve the layered composites under different indentation conditions. No theoretical solution can be derived from the final solutions of the stress and displacement except for a very simple case, since the solution carries integration from zero to infinity. Therefore a numerical method is used to evaluate these integrations. Cases of half space, single layer bonded to the rigid base, single layer bonded to an elastic half space, and two-layer bonded to rigid base are investigated in this research. Based on the research, mathematical models are presented for each different case. The mathematical models provide simple guidelines for the layered composite design. In the advanced study section, the shear loading at the contact interface was studied. With co-relationship between shear stress and normal stress by the friction coefficient, friction effect can be studied at the contact interface. Finally, shear slip and normal separation theories were incorporated into the technique developed previously. This study provides some insight into the debonding for the layered composite. 7.2 Conclusions With the theoretical and mathematical development, we are able to study a few parameters that have an effect on the indentation. 148 A half space was used to study the material effect. We validated that material properties have no influence to the stress distribution, while the softer material intends to have higher displacement. An elastic layer perfectly bonded to the elastic half base and indented with spherical indenter was used to study the thickness effect to the stress and displacement. We found that the normal stress distributions at the contact surface were identical for layers with different thicknesses. Apparently, the contact stress distribution is a local phenomenon and is not affected by the thickness of the layer underneath it. As the top elastic layer becomes thicker, the normal stress distribution at the same depth below the contact surface becomes higher. This is because the elastic half base, which is made of nylon, has lower stiffness than the top elastic layer, which is made of aluminum, and tends to release the normal stress distribution. The shear stress distributions vanish identically on the contact surface for all cases due to the free shear traction boundary condition. When the top elastic layer becomes thicker, the shear stress distribution at the same depth becomes lower. This phenomenon is opposite to that of normal stress distribution. The radial displacement distributions change from negative to positive at the contact surface. Whether it is positive or negative, the thinnest top layer always has the largest deformation. The normal displacement distributions are related to the shape of the spherical indenter. As the top layer becomes thicker, the vertical displacement becomes less. The same case, an elastic layer perfectly bonded to elastic half space indented by a spherical indenter, was used to investigate the boundary condition. We 149 found that the normal stress distributions at the contact surface for all three combinations are exactly the same, indicating again the local phenomenon at the contact surface. As the stiffness of the half base increases, the normal stress increases. The shear stress distributions vanish at the contact surface due to the free shear traction boundary condition. When the stiffness of the half base increases, the shear stress distribution decreases for the depth z= 0.5 mm, as discussed in Chapter Four. However, it increases for the depth 2= 1 mm. This is because a harder base can sustain a higher stress than a softer base. The displacement distributions, both radial and normal components have the similar relationship. We studied the effect of loading condition on stress and displacement distributions. The half space problem (made of aluminum) is chosen along with various load conditions, eg. near-point load, uniform stress, flat punch and spherical loading, for the investigations. All cases have an identical total force although the geometry of the indenters is different. At the contact surface, we see numerical oscillations for all cases except the one subjected to a near-point load, which has an exact solution. The spherical indenter has the highest maximum normal stress followed by flat punch and then by uniform stress. The near-point load gives the highest maximum normal stress. It should be noted that the maximum normal stress is located away from the origin except for the near-point load case. The effect of lamination on the stress and displacement distributions was studied. The two-layer composite situated on a rigid half base and subjected to a 150 spherical indenter was chosen for investigations. Two cases were studied: the ratio of Young’s moduli and the ratio of thicknesses between the two layers. For the ratio of Young’s moduli, we used nylon/aluminum (the ratio E1/E2= 0.0356) and aluminum/nylon combinations. For the ratio of thickness, we based it on nylon/aluminum. More specifically, we fixed the thickness of the top nylon layer and changed the thickness of the bottom aluminum layer so the relative location of interest to the loading surface can be maintained. From the calculations, we can see the stress distributions do change with the thickness ratio. The displacement distributions, however, change with the thickness ratio. As the thickness ratio between the first layer and the second layer decreases, the radial displacement distributions increase slightly while the normal displacement decreases. Similarly, we investigate the stacking order of aluminum/nylon, E1/E2= 28.1. As the thickness ratio decreases, the normal stress decreases. However the shear stress increases, and both radial and normal displacement distributions increase. Furthermore, we developed some simple indentation models for the indentation condition developed previously. For half space, we came up with an exact theoretical solution and the following relationships among loading force, material properties, indenter shape, contact area and maximum deflection: a: Near point loading: P = 7! C; E * Uz (l-v ) . (Jr * a)2 Uniform pressure: P = —TE * Uz (l-v ) 151 a Flat indentation: P = 2 E * Uz (l-v ) Spherical indentation: P = 4a 2 E *Uz 3(l-v ) From the numerical solution, we also developed simple models for the layered composites for a spherical indenter shape. The following models are presented: Single layer indented against a rigid base: P=3.8149"‘E"‘Uz"‘(fl —O.6389 a) Single layer indented against an elastic base: E —0.4529 ,1 —0,]329 P =1.0788*U * 152 *[ij {—1) E] a Two layered composites indented against a rigid base: E —0.427 It —0.183l P=1.2425*U*E2 {—2) [A] El 112 As advanced study, we investigated the shear loading at the contact interface. This allowed the author to investigate the friction effect on the internal stress and displacement. We found that both stress and displacement increase with the increase of the friction. The shear slip and normal separation theories were incorporated into the early development of this research. The results show that when kr is very small (less 152 than 0.00001), the stress and displacement is close to the perfect bounding. The shear slip coefficient kr can affect the shear stress more significantly than the normal stress and has more effect on the shear displacement than the normal displacement. When the shear slip coefficient kr was between 0.0001 and 0.01, some sinusoidal oscillations were observed in the numerical results. These values of kr represent changes of interface from perfect bonding to complete debonding. The sinusoidal oscillations were caused by the slow convergence of the Bessel function in the numerical integrations. For the effect of normal separation coefficient, we found that the bonding condition is near perfect when k2 is very small, i.e. less than 0.00001. When kz increases, the normal stress at the bonding interface decreases significantly; however, the shear stress is not as much affected by the change of k2 as the normal stress. Similarly, the normal displacement at the bonding interface also decreases more significantly than the shear displacement. Moreover, we extended the shear slip and normal separation theories to the finite imperfect bonding interface with the introduction of a step function for the shear slip and normal separation coefficients. We investigated the effect from both coefficients. We found that shear slip coefficient kr doesn’t pose a significant effect on the normal stress and displacement. However, kr affects the shear stress and shear displacement significantly. Some sinusoidal oscillations were observed when kr was between 0.0005 and 0.01. Beyond this range, the results were close to perfect bonding when kr was very small and the shear stress was very close to zero when the interface was completely debonded. Similarly, we 153 found that the normal separation coefficient kz has a significant effect on both normal stress and normal displacement, while the influence of kz on the shear stress and shear displacement is not significant. 7.3 Suggestions and Future Study Contact mechanics involves mathematical analysis under some mechanics laws. The development in this dissertation was aimed at filling some voids that were not addressed in the literature. Sneddon’s (1951) work was used as the base in the research development, i.e. using Hankel transformation to solve the biharmonic equations. With the use of computer software packages such as Mathematica® and Microsoft Visual Studio ®, the goal to tackle the sophisticated contact mechanics problems became possible, though tremendous efforts in mathematical manipulation were still necessary. However, Hankel transformation did not allow the technique to be used for other composites, such as anisotropic and/or inhomogeneous fiber composites. The technique was limited to axisymmetric problems. This research in this dissertation touched some fundamental issues. Many related issues are still up for grasp. In this research, we only developed contact models involving a spherical indenter. Other types of indenters should also be investigated. Furthermore, in order to mathematically solve the imperfect bonding with finite length, specially transformed shear slip and normal separation relationships were used. It would be worthwhile to explore other relationships that allow theoretical solutions to the imperfect bonding problems. 154 Finally, a theoretical study for specially mixed boundary conditions would prove to be very useful in advancing the contact mechanics. It would allow the true boundary condition to be exactly matched instead of the semi-guessing involved in the current practice. 155 APPENDICES 156 Appendix A Mathmatica ® C for General Derivation H=EO/(2*(l+v)) A=EO*v/((l+v)*(1-2*v)) G0=(l+v)(1-2*v)/(E0*§‘2)((A+B*z)Exp[- §*z]+(C+D*z)Exp[§*z])*ph Gl=Co11ect[Simplify[D[GO,z]],{Exp[-§*z],Exp[§*z]}] G2:Co11ect[Simplify[o[cl,z]],{Exp[-§*z],Exp[§*z]}] G3=Collect[Simplify[D[G2,z]],{Exp[-§*z],Exp[§*z]}] oz=Collect[Simplify[§*((A+2*u)G3- (3*A+4*u)*§‘2*Gl)*BesselJ[0,§*r]],{Exp[-§*z],Exp[§*z]}] trz=Collect[Simplify[€‘2*(A*G2+(A+2*u)*§‘2*G0)*Besse1J[1,§*r ]],{Exp[-§*z],Exp[§*z]}] or1=Collect[Simplify[§*(A*G3+(A+2*u)*§*2*Gl)*Besse1J[o,§*r]] .{Exp[-§*z],Exp[§*z]}] or2=Collect[Simplify[2(1+0)/r*§‘2*G1*BesselJ[l,§*r]],{Exp[- §*z],Exp[§*z]}] oel=cO11ect[Simplify[A*§*(G3—§‘2*G1)*Besse1J[0,§*r]],{Exp[- §*z] ,Exp [5*2] }] 092=Collect[Simplify[2(A+u)/r*§‘2*Gl*BesselJ[l,§*r]],{Exp[- §*ZJ,EXpl€*ZJ}] ur=Collect[Simplify[(1+0)/u*§‘2*Gl*BesselJ[l,§*r]],{Exp[- 6*21 .EXp[§*21 }] 157 uz=Collect[Simplify[§*(G2(A+2*u)/u*§‘2*G0)*BesselJ[0,§*r]], {EXp[-€*z],EXp[§*z]}l 158 Appendix B Mathematica® Code for Half Space Half space - General Derivation 138: (1+ v) (1-2*v) / (HHE‘Z) Lwrwp (r) BasaelJ[0, Er] dlr C0=0 D0=0 u=E0/(2*(l+v)) A=E0*v/((l+v)*(l-2*v)) G0=(A+B*z)Exp[-§*z]+(C+D*z)Exp[§*z] G1=Collect[D[G0,z], {Exp[-§*z],Exp[§*z]}] G2=Collect[D[Gl,z], {Exp[-§*z],Exp[§*z]}] GB=Collect[D[G2,z], {Exp[-§*z],Exp[§*z]}] ozit=Simplify[§((A+2*u)G3-(3*A+4*u)*§*2*G1)*BesselJ[0,§*r]] trzit=Simplify[§‘2(A*G2+(A+2*u)*§‘2*G0)*BesselJ[l,§*r]] orit=Simplify[§(A*G3+(A+2*u)*5‘2*G1)*BesselJ[0,§*r]- 2(A+u)/r*§‘2*Gl*BesselJ[l,§*r]] oeit=Simplify[A*§(G3- §‘2*Gl)*BesselJ[0,§*r]+2(A+u)/r*§‘2*G1*BesselJ[l,§*r]] urit=Simplify[(A+u)/u*§‘2*Gl*BesselJ[l,§*r]] uzit=Simplify[§(G2—(A+2*u)/u*§‘2*G0)*BesselJ[0,§*r]] oz = foazita: pI-Id§ 0 :12 = farmitar pl-IdE O or: Jonit*wd§ 0 159 oe= fmwit*plid§ 0 ur=Jmurit1~pHd§ 0 112 = Inuit. pHd§ 0 160 Appendix C Mathematica® Code for Point Loading Half space - point loading C0=0 D0=0 A0=2*v/§ BO=1 Ph=(-P/2/N)EXpl-a*€] u=E0/(2*(l+v)) A=E0*v/((l+v)*(1—2*v)) G0=(1+v)(1-2*v)/(Eo*§*2)((A0+Bo*z)Exp[— §*z]+(C0+D0*z)Exp[§*z])*ph G1=Simplify[D[G0,z]] G2=Simplify[D[G1,z]] G3=Simplify[D[G2,z]] ozit=Simplify[§((A+2*u)G3-(3*A+4*u)*§‘2*Gl)*BesselJ[0,§*r]] trzit=Simplify[§‘2(A*G2+(A+2*u)*é‘2*G0)*BesselJ[l,§*r]] orit=SimplifY[€(A*G3+(A+2*u)*§‘2*Gl)*BesselJ[0,§*r]- 2(A+u)/r*§‘2*Gl*BesselJ[l,§*r]] oeit=Simplify[A*§(G3— §*2*G1)*Besse1J[0,§*r]+2(A+p)/r*§‘2*G1*Besse1J[1,§*r]] urit=Simplify[(A+u)/u*§‘2*Gl*BesselJ[1,§*r]] uzit=Simplify[§(Gz-(A+2*u)/u*§‘2*G0)*Besse1J[0,§*r]] #include #include #include #include #include #include #include #include //time calculation headfile #include using namespace std; main(void) { ”Running time calculation start time_t BeginSec, EndSec, Total; // calculate running time time(&BeginSec); //**** Setup input and output file char lnFile1[80], OutFile[80]; printf("\n Please enter the input file name: "); gets(lnFile1); printf("\n Please enter the output file name: "); gets(OutFile); //****Error message of opening files constchar ErrorMsg[]="\n***main: unable to open file: "; fstream lnMatl(lnFile1,ios::in), OutComp(OutFile,ios::out); if (lnMatl.fail()){ //detect if the file can be open cerr<>i; lnGauss>>x[i]; lnGauss>>w[i]; lnGauss.close(); cout<31;i—){ x[i]=x[i-32]; w[i]=w[i-32]; } l/2)establish 0-31 (-1,0) backward from 32-64 for(i=0;i<32;i++){ x[i]=-x[63-i]; w[i]=w[63-i]; } l/****Read imput file double E=0,Nu=0,hc=0.0,a=0.0,p0=0.0,P=0.0; int indenter=0; double r___max=0; // max value of z and r double a0=0, b0=0; l/integration range int Nz=0,Nr=0,Nab=0; ”number of points for z and r while(!lnMatl.eof()){ lnMatl>>indenter; // lndenter Shape lnMatl.ignore(80,'\n'); lnMatl>>E; // Young’s modulus 174 lnMatl.ignore(80,'\n'); lnMatl>>Nu; // Possion ratio lnMatl.ignore(80,’\n'); lnMatl>>hc; // layer thickness lnMatl.ignore(80,'\n'); lnMatl>>a; ”contact radius lnMatl.ignore(80,'\n'); lnMatl>>a0; // integration lower limit lnMatl.ignore(80,'\n'); InMatl>>b0; // integration upper limit lnMatl.ignore(80,'\n'); lnMatl>>Nab; // number of integration zones in [a0,b0] lnMatl.ignore(80,'\n'); lnMatl>>r_max; // ----9 lnMatl.ignore(80,'\n'); lnMatl>>Nz; //number of points for z -------10 lnMatl.ignore(80,'\n'); lnMatl>>Nr; // number of points for r --—----11 lnMatl.ignore(80,'\n'); lnMatl>>P; lnMatl.ignore(80,'\n'); l/calculate p0 if (indenter==1){ p0=P/(3.14159*pow(a,2));//***uniform pressure else if (indente ==2){ //***flat punch - uniform displacement p0=P/(2*3.14159*pow(a,2)); else if (indenter==3){ p0=3*P/(2*3.14159*pow(a,2)); //***Spherical lndenter } lnMatl.close(); OutComp<<"Half Space Indentation" <<"\n\nlndenter shape="< #include #include #include #include #include #include #include #include #include using namespace std; double *f_integrant(int indenter,double E,double Nu,double h,double a,double z,double r,double xi,double p0,double P,double fc[]); main(void) { time_t BeginSec, EndSec, Total; // calculate running time time(&BeginSec); //**** Setup input and output file char lnFile1[80], OutFile[80]; printf(”\n Please enter the input file name: "); gets(lnFile1); printf("\n Please enter the output file name: "); gets(OutFile); I/****Error message of opening files constchar ErrorMsg[]="\n***main: unable to open file: "; fstream lnMatl(lnFile1,ios::in), OutComp(OutFile,ios::out); if (lnMatl.fail()){ //detect if the file can be open cerr<>i; . lnGauss>>x[i]; lnGauss>>w[i]; lnGauss.close(); cout<31;i-){ x[i]=x[i-32]; w[i]=w[i-32]; } for(i=0;i<32;i++){ x[i]=-x[63-i]; w[i]=w[63-i]; } //****Read imput file double E=0,Nu=0,h=0.0,hc=0,a=0.0,p0=0.0,P=0.0; int indenter=0; double r_max=0; // max value of z and r double a0=0, b0=0; ”integration range int Nz=0,Nr=0,Nab=0; ”number of points for z and r while(!lnMatl.eof()){ lnMatl>>indenter; // lndenter Shape lnMatl.ignore(80,'\n'); lnMatl>>E; ”Young's modulus lnMatl.ignore(80,'\n'); lnMatl>>Nu; //Possion ratio lnMatl.ignore(80,'\n'); 182 } lnMatl>>h; // layer thickness lnMatl.ignore(80,'\n'); lnMatl>>hc; // calculated layer thickness lnMatl.ignore(80,'\n'); lnMatl>>a; lnMatl.ignore(80,'\n'); //contact radius lnMatl>>a0; lnMatl.ignore(80,'\n'); // integration lower limit lnMatl>>b0; // intergation upper limit lnMatl.ignore(80,'\n'); lnMatl>>Nab; /l number of integration zones in [a0,b0] lnMatl.ignore(80,'\n'); lnMatl>>r_max; lnMatl.ignore(80,'\n’); ”calculation zone of max r lnMatl>>Nz; //number of points for z lnMatl.ignore(80,'\n'); lnMatl>>Nr; // number of points for r lnMatl.ignore(80,'\n'); lnMatl>>P; lnMatl.ignore(80,'\n'); lnMatl.close(); //*****calculate p0 if (indenter==1){ p0=P/(3.14159*pow(a,2)); else if (indenter==2){ p0=F>/(2*3.14159*pow(a.2)): else if (indenter==3){ } p0=3*P/(2*3-14159*P°W(3:2))§ // total load P l/***uniform pressure //***flat punch - uniform displacement //***Spherical lndenter OutComp<<"Single Layer Indentation" <<"\n\nlndenter shape="< #include #include #include #include #include #include #include I/time calculation headfile #include double *f_integrant(int indenter,double E1 ,double E2,double Nu1,double Nu2,double h1,double h2, double a, double z,double r,double xi,double p0,double P,double fcfl); main(void) { time_t BeginSec, EndSec, Total; // calculate running time time(&BeginSec); //**** Setup input and output file char lnFile1[80], OutFile[80]; printf("\n Please enter the input file name: "); gets(lnFile1); printf("\n Please enter the output file name: "); gets(OutFile); I/****Error message of opening files constchar ErrorMsg[]="\n***main: unable to open file: "; fstream lnMatl(lnFile1,ios::in), OutComp(OutFile,ios::out); if (lnMatl.fail()){ //detect if the file can be open cerr<>i; lnGauss>>x[i]; lnGauss>>w[i]; lnGauss.close(); cout<31;i-){ x[i]=x[i-32]; w[i]=w[i-32]; } //2)establish 0-31 (-1,0) backward from 32-64 for(i=0;i<32;i++){ x[i]=-x[63-i]; w[i]=w[63-i]; } //****Read imput file double E1=0,E2=0,Nu1=0,Nu2=0,h1=0.0,h2=0.0,a=0.0,p0=0.0,P=0.0; int indenter=0; double r_max=0; // max value of r double a0=0, b0=0; ”integration range int Nz1=0,N22=0,Nr=0,Nab=0; //number of points for z and r while(!lnMatl.eof()){ lnMatl>>indenter; // layer thickness lnMatl.ignore(80,'\n'); lnMatl>>E1; //Young's modulus lnMatl.ignore(80,'\n'); lnMatl>>E2; ”Young's modulus 190 lnMatl.ignore(80,'\n'); lnMatl>>Nu1; lnMatl.ignore(80,'\n'); lnMatl>>Nu2; lnMatl.ignore(80,’\n'); lnMatl>>h1; lnMatl.ignore(80,'\n'); lnMatl>>h2; lnMatl.ignore(80,'\n'); lnMatl>>a; lnMatl.ignore(80,'\n'); lnMatl>>a0; lnMatl.ignore(80,'\n'); lnMatl>>b0; lnMatl.ignore(80,'\n'); lnMatl>>Nab; lnMatl.ignore(80,'\n'); lnMatl>>r_max; lnMatl.ignore(80,'\n'); lnMatl>>Nz1; lnMatl.ignore(80,'\n'); lnMatl>>N22; lnMatl.ignore(80,'\n'); lnMatl>>Nr; lnMatl.ignore(80,'\n'); lnMatl>>p0; lnMatl.ignore(80,'\n'); lnMatl>>P; lnMatl.ignore(80,'\n'); //Possion ratio //Possion ratio // layer thickness // layer thickness l/contact radius // integration lower limit // integration upper limit // number of integration zones in [a0,b0] ll max r value //number of points for 2 //number of points for z // number of points for r // stress constant // total force lnMatl.close(); //*****calculate p0 if (indenter==1){ p0=P/(3.14159*pow(a,2));/l***uniform pressure else if (indenter==2){ //***flat punch - uniform displacement p0=P/(2.0*3.14159*pow(a,2)); else if (indenter==3){ p0=3.0*P/(2*3.14159*pow(a,2)); //***Spherical lndenter } OutComp<<"Two Half Space Perfectly Bonded at z=0"; OutComp<<"\n\nlndenter shape="< #include #include #include #include #include #include #include #include lltime calculation headfile #include using namespace std; double *f_integrant(int indenter,double E1 ,double E2,double Nu1, 204 double Nu2,double h1, double a, double z,double r,double xi, double p0,double P,double f,double fc[]); main(void) { time_t BeginSec, EndSec, Total; // calculate running time time(&BeginSec); //**** Setup input and output file char lnFilel [80], OutFile[80]; printf("\n Please enter the input file name: "); gets(lnFile1); printf("\n Please enter the output file name: "); gets(OutFile); //****Error message of opening files constchar ErrorMsg[]="\n***main: unable to open file: "; fstream lnMatl(lnFile1,ios::in), OutComp(OutFile,ios::out); if (lnMatl.fail()){ //detect if the file can be open cerr<>i; lnGauss>>x[i]; lnGauss>>w[i]; lnGauss.close(); cout<31;i-){ x[i]=x[i-32]; w[i]=w[i-32]; } //2)establish 0-31 (-1,0) backward from 32-64 for(i=0;i<32;i++){ x[i]=-x[63-i]; w[i]=w[63-i]; } //****Read imput file double E1=0,E2=0,Nu1=0,Nu2=0,h1=0.0,hc1=0.0,h02=0.0,a=0.0, p0=0.0,P=0.0,f=0; int indenter=0; double r_max=0; // max value of r double a0=0, b0=0; Ilintegration range int Nz1=0,N22=0,Nr=0,Nab=0; I/number of points for z and r while(!lnMatl.eof()){ lnMatl>>indenter; // layer thickness lnMatl.ignore(80,'\n'); lnMatl>>E1; ”Young's modulus lnMatl.ignore(80,'\n'); lnMatl>>E2; //Young's modulus lnMatl.ignore(80,'\n'); lnMatl>>Nu1; l/Possion ratio lnMatl.ignore(80,'\n'); lnMatl>>Nu2; //Possion ratio lnMatl.ignore(80,'\n'); lnMatl>>h1; // layer thickness lnMatl.ignore(80,'\n'); lnMatl>>hc1; II calculated thickness in layer one lnMatl.ignore(80,'\n'); lnMatl>>hc2; // calculated thickness in layer two lnMatl.ignore(80,'\n'); lnMatl>>a; l/contact radius lnMatl.ignore(80,'\n'); lnMatl>>a0; // integration lower limit 206 lnMatl.ignore(80,'\n'); lnMatl>>b0; // integration upper limit lnMatl.ignore(80,'\n'); lnMatl>>Nab; // number of integration zones in [a0,b0] lnMatl.ignore(80,'\n'); lnMatl>>r_max; // max r value lnMatl.ignore(80,'\n'); lnMatl>>Nz1; ”number of points for z lnMatl.ignore(80,'\n'); lnMatl>>N22; //number of points for z lnMatl.ignore(80,'\n'); lnMatl>>Nr; // number of points for r lnMatl.ignore(80,'\n'); lnMatl>>P; // total force lnMatl.ignore(80,'\n'); lnMatl>>f; // friction coefficient lnMatl.ignore(80,'\n'); lnMatl.close(); //*****calculate p0 if (indenter==1){ p0=P/(3.14159*pow(a,2)); //***uniform pressure else if (indenter==2){ //***flat punch - uniform displacement p0=P/(2.0*3.14159*pow(a,2)); } else if (indenter==3){ p0=3.0*P/(2*3.14159*pow(a,2)); //***Spherical lndenter } OutComp<<"Two Layer Composites Half Space with Shear Force at Boundary Condition"; OutComp<<"\n\nlndenter shape="< #include #include #include #include #include #include #include #include lltime calculation headfile #include //using namespace std; double *f_integrant(int indenter,double E1 ,double E2,double Nu1, double Nu2,double h1,double h2,double a, double z,double r, double xi,double p0,double P,double kr,double kz,double d,double fc[]); main(void) { time_t BeginSec, EndSec, Total; // calculate running time time(&BeginSec); //**** Setup input and output file char lnFile1[80], OutFile[80]; printf("\n Please enter the input file name: "); gets(lnFile1); printf("\n Please enter the output file name: "); gets(OutFile); //****Error message of opening files constchar ErrorMsg[]="\n***main: unable to open file: "; fstream lnMatl(lnFile1,ios::in), OutComp(OutFile,ios::out); if (lnMatl.fail()){ //detect if the file can be open cerr<>i; lnGauss>>x[i]; lnGauss>>w[i]; lnGauss.close(); cout<31;i-—){ x[i]=x[i-32]; w[i]=w[i-32]; } II2)establish 0-31 (-1,0) backward from 32-64 for(i=0;i<32;i++){ x[i]=-x[63-i]; w[i]=w[63-i]; } //****Read imput file double E1=0,E2=0,Nu1=0,Nu2=0,h1=0.0,h2=0.0,a=0.0,p0=0.0,P=0.0; int indenter=0; double r_max=0; // max value of r double a0=0, b0=0; ”integration range int N21 =0,N22=0,Nr=0,Nab=0; //number of points for z and r 219 double kr=0.0, kz=0.0, d=0.0; //kr (k1)-coefficient Ur, kz(k2) - coefficient for U2 while(!lnMatl.eof()){ lnMatl>>indenter; lnMatl.ignore(80,'\n'); lnMatl>>E1; lnMatl.ignore(80,'\n'); lnMatl>>E2; lnMatl.ignore(80,'\n'); lnMatl>>Nu1; lnMatl.ignore(80,'\n'); lnMatl>>Nu2; lnMatl.ignore(80,'\n'); InMatl>>h1 ; lnMatl.ignore(80,'\n'); lnMatl>>h2; lnMatl.ignore(80,'\n'); lnMatl>>a; lnMatl.ignore(80,'\n'); lnMatl>>a0; lnMatl.ignore(80,'\n'); lnMatl>>b0; lnMatl.ignore(80,'\n'); lnMatl>>Nab; lnMatl.ignore(80,'\n'); lnMatl>>r_max; lnMatl.ignore(80,'\n'); lnMatl>>Nz1; lnMatl.ignore(80,'\n'); lnMatl>>N22; lnMatl.ignore(80,'\n'); lnMatl>>Nr; lnMatl.ignore(80,'\n'); lnMatl>>P; lnMatl.ignore(80,'\n'); lnMatl>>kr; lnMatl.ignore(80,'\n'); lnMatl>>kz; lnMatl.ignore(80,'\n'); lnMatl>>d; lnMatl.ignore(80,'\n'); // layer thickness ”Young's modulus ”Young's modulus llPossion ratio //Possion ratio // layer thickness // layer thickness llcontact radius // intergation upper limit // number of integration zones in [a0,b0] // max r value //number of points for 2 //number of points for z // number of points for r // total force llcoefficient along r axis llcoefficient along 2 axis //crack radius lnMatl.close(); //*****calculate p0 if (indenter==1){ 220 p0=P/(3.14159*pow(a,2)); //***uniform pressure else if (indenter==2){ //***flat punch - uniform displacement p0=P/(2.0*3.14159*pow(a,2)); else if (indenter==3){ p0=3.0*P/(2*3.14159*pow(a,2)); //***Spherical lndenter OutComp<<"Two Layer Composites Half Space with Shear Slip Theory within aa at bonding interface"; OutComp<<"\n\nlndenter shape="<