Simulation on the process of fatigue crack initiation in a martensitic stainless steel Vom Fachbereich Maschinenbau der Universität Kassel zur Erlangung des akademischen Grades eines Doktor-Ingenieurs (Dr.-Ing.) genehmigte DISSERTATION von Master of Engineering Xinyue Huang Hauptreferent: Prof. Dr. rer. nat. Angelika Brückner-Foit Koreferent: Dr.-Ing. Igor Altenberger Prüfer: Prof. Dr.-Ing. Berthold Scholtes Prüfer: Prof. Dr. Xueren Wu Tag der mündlichen Prüfung: 18.4.2007 Tag der Einreichung: 26.4.2007 Acknowledgement The present work has been carried out at the Division of Quality and Reliability, Institute of Material Engineering, Department of Mechanical Engineering, University of Kassel. I would like to express my deep attitude and appreciate towards all those who have helped to make the work finish with success. I am grateful to Prof. Angelika Brückner-Foit for giving me the opportunity to work in the Institute under her guidance and supervision. Her advices, comments and discussions have always been very helpful and stimulating. Financial support by the German Research Foundation (Deutsche Forschungsgemeinschaft) is gratefully acknowledged. I sincerely thank my colleague, Dr. Stefanie Anteboth, for her suggestions and for her help in finite element software, Dr. Yasuko Motoyashiki and Mr. Micheal Besel, for their comments, discussions and suggestions. I would like to extend my thanks to Mrs. Heike Hammann for her continuous support in administration matters and Mr. Ralf Herbold for his maintenance of the computer system. To all the colleagues in the Institute, I express many thanks for their help, suggestions, as well as the pleasant and friendly atmosphere. I thank my family, for their support through these years. i ABSTRACT The present research deals with the computer simulation for the microcrack initiation process for a martensitic steel subjected to low-cycle fatigue. As observed on the specimen surface, the initiation and early propagation of these microcracks are highly microstructure dependent. This fact is taken into account in the mesoscopic damage accumulation models in which the grains are modelled as single crystals with anisotropic material behaviour. The representative volume element generated by a Voronoi tessellation process is used to simulate the microstructure of the polycrystalline material. Stress distributions are analyzed by a finite element method with elastic and elasto-plastic material properties. The simulation is first carried out on two-dimensional models and then on a simplified three-dimensional model where the three-dimensional slip system and stress state are taken into account. Continuous crack initiation is simulated by defining potential crack path within each grain and the number of cycles to crack initiation is estimated on the basis of the Tanaka-Mura and the Chan equations. The simulation model yields the relation of crack densities versus the number of cycles and the results are compared with experimental data. For all of the strain ranges considered the simulation results coincide well with the experiment data. ii KURZFASSUNG Die vorliegende Arbeit beschäftigt sich mit der Computersimulation des Rissinitiierungsprozesses für einen martensitischen Stahl, der der niederzyklischen Ermüdung unterworfen wurde. Wie auf der Probenoberfläche beobachtet wurde, sind die Initiierung und das frühe Wachstum dieser Mikrorisse in hohem Grade von der Mikrostruktur abhängig. Diese Tatsache wurde in mesoskopischen Beschädigungsmodellen beschrieben, wobei die Körner als einzelne Kristalle mit anisotropem Materialverhalten modelliert wurden. Das repräsentative Volumenelement, das durch einen Voronoi-Zerlegung erzeugt wurde, wurde benutzt, um die Mikrostruktur des polykristallinen Materials zu simulieren. Spannungsverteilungen wurden mit Hilfe der Finiten-Elemente-Methode mit elastischen und elastoplastischen Materialeigenschaften analysiert. Dazu wurde die Simulation zunächst an zweidimensionalen Modellen durchgeführt. Ferner wurde ein vereinfachtes dreidimensionales RVE hinsichtlich des sowohl dreidimensionalen Gleitsystems als auch Spannungszustandes verwendet. Die kontinuierliche Rissinitiierung wurde simuliert, indem der Risspfad innerhalb jedes Kornes definiert wurde. Die Zyklenanzahl bis zur Rissinitiierung wurde auf Grundlage der Tanaka-Mura- und Chan- Gleichungen ermittelt. Die Simulation lässt auf die Flächendichten der einsegmentige Risse in Relation zur Zyklenanzahl schließen. Die Resultate wurden mit experimentellen Daten verglichen. Für alle Belastungsdehnungen sind die Simulationsergebnisse mit denen der experimentellen Daten vergleichbar. iii CONTENT PREFACE .…………………………………………………………………………………… 1 CHAPTER 1 INTRODUCTION …………………………………………………………… 3 1.1 Fatigue Behavior and Fatigue Tests ……………………………………………………… 3 1.1.1 Three Stages of Fatigue .……………………………………………………………… 3 1.1.2 Fatigue Test - Strain Cycling and Stress Cycling ...……..…………………………… 6 1.1.3 Damage Accumulation during Multiple Crack Initiation ………………………… 7 1.2 Mechanism of Crack Initiation …………………………………………………………… 9 1.2.1 Mechanism of PSB Formation ……………………………………………………… 9 1.2.2 Mechanism of Crack Initiation from PSB ………………………………………… 11 1.2.3 Mechanism of Crack Initiation from Inclusions …………………………………… 12 1.3 Models of Crack Initiation ……………………………………………………………… 13 1.3.1 Conventional Models ……………………………………………………………… 13 1.3.2 Microstructure-Based Models ……………………………………………………… 14 1.3.3 Models Based on Probability ……………………………………………………… 17 1.4 Modeling of Polycrystal Materials ……………………………………………………… 18 1.4.1 Representative Volume Element …………………………………………………… 18 1.4.2 Mesoscopic Mosaic Models ……………………………………………………… 19 CHAPTER 2 EXPERIMENTAL DATA AND STATISTICAL ANALYSIS …………… 21 2.1 Material ………………………………………………………………………………… 21 2.1.1 Microstructure ………………………………………………………………… 22 iv 2.1.2 Mechanical Properties ……………………………………………………………… 23 2.2 Low Cycle Fatigue Tests ………………………………………………………………… 23 2.3 Experiment Results ……………………………………………………………………… 25 2.3.1 Fatigue Life ……………………………………………………………………… 25 2.3.2 Elasto-Plastic Behavior Obtained from Experiment Data ………………………… 25 2.3.3 Cyclic Deformation Behavior of F82H …………………………………………… 27 2.4 Observation on the Surface of Fatigue Specimens ……………………………………… 28 2.4.1 Morphology of Microcracks on Specimen Surface ……………………………… 28 2.4.2 Statistics for Characteristics of Microcracks ……………………………………… 31 2.4.3 Characteristics of One-Segment Cracks …………………………………………… 33 2.4.3.1 Crack Length ………………………………………………………………… 33 2.4.3.2 Crack Orientation …………………………………………………………… 33 2.4.3.3 Crack Density as Function of Cycles ……………………………………… 33 2.5 Characteristics of Crack Initiation ……………………………………………………… 34 2.6 Scatter of Experimental Data ………………………………………………………………35 CHAPTER 3 IDEAS AND HYPOTHESES OF MODELING ………………………… 37 3.1 Material Model ………………………………………………………………………… 37 3.2 Fatigue Model …………………………………………………………………………… 38 3.3 Parameter Studies …………………………………………………………………… 39 3.3.1 Critical Shear Stress Study ..………………………………………………………… 39 3.3.2 Microstructure Parameter Study …..………………………………………………… 40 CHAPTER 4 CONSTRUCTION OF SIMULATION MODELS ……………………… 41 v 4.1 Model Outline …………………………………………………………………………… 41 4.2 Representative Volume Element Model ………………………………………………… 43 4.2.1 Determination of Slip System ……………………………………………………… 43 4.2.2 2D-RVE Model ……………………………………………………………………… 45 4.2.3 3D-RVE Model ……………………………………………………………………… 47 4.2.4 RVE Size and Voronoi Boundary Effects ...………………………………………49 4.3 Model for Finite Element Analysis ……………………………………………………… 51 4.3.1 Coordinate Systems ………………………………………………………………… 51 4.3.1.1 Coordinate Systems of 2D-RVE Model …………………………………… 51 4.3.1.2 Coordinate Systems of 3D-RVE Model …………………………………… 52 4.3.2 Boundary Conditions ……………………………………………………………… 54 4.3.3 Element and Mesh ………………………………………………………………… 55 4.3.3.1 Element and Mesh of Two-Dimensional Model .…………………………… 56 4.3.3.2 Element and Mesh of Three-Dimensional Model …………………………… 57 4.4 Material Properties ……………………………………………………………………… 58 4.4.1 Stress-Strain Response of Elastic Material ..………………………………………… 58 4.4.2 Stress-Strain Response of Elasto-Plastic Material ..………………………………… 59 4.5 Modeling of Crack Initiation Process ...………………………………………………… 62 4.5.1 Fatigue Model ……………………………………………………………………… 62 4.5.2 Average Resolved Shear Stress …………………………………………………… 62 4.5.2.1 Transformation of Stress Tensors …………………………………………… 62 4.5.2.2 Average Resolved Shear Stress ...…………………………………………… 63 vi 4.5.3 Crack Initiation Process …………………………………………………………… 66 4.5.4 Summary of Simulation Procedures and Applied Criteria .………………………… 67 4.6 Verification of simulation model ………………………………………………………… 68 4.6.1 Similarity of Mosaic Model to Studied Material …………………………………… 71 4.6.1.1 Structure of 2D-RVE ……………………………………………………… 71 4.6.1.2 Structure of 3D-RVE ……………………………………………………… 73 4.6.2 Stress-Strain Response of Models ………………………………………………… 74 CHAPTER 5 RESULTS OF 2D SIMULATIONS ..……………………………………… 76 5.1 Stress Distribution in Uncracked RVE …………………………………………………… 77 5.1.1 Stress Distribution in Elastic Models ………………………………………………… 77 5.1.1.1 Von Mises Stress Distribution …………………………………………… 77 5.1.1.2 Shear Stress Distribution ………………………………………………… 80 5.1.2 Stress Distribution in Elasto-Plastic Model …………………………………………81 5.2 Relations of Crack Density versus Number of Cycles …………………………………… 82 5.2.1 Tentative Parameters ……………………………………………………………… 82 5.2.2 Parameters Study of Elasto-Plastic Model …………………………………………… 83 5.2.2.1 Critical Shear Stress ...……………………………………………………… 83 5.2.2.2 Shear Modulus ...…………………………………………………………… 86 5.2.3 Relation of Crack Density to Cycles ………………………………………………… 87 5.2.4 Effect of Microstructures …………………………………………………………… 89 5.3 Effect of Stress Redistribution on Crack Initiation Sequence …………………………… 91 5.3.1 First Example ………………………………………………………………………… 91 vii 5.3.2 Second Example …………………………………………………………………… 93 5.3.3 Third Example …………………………………………………………………… 94 5.4 Crack Patterns of Elasto-Plastic Model …………………………………………………… 96 CHAPTER 6 RESULTS OF 3D SIMULATIONS ………………………………………… 99 6.1 Stress Distribution in Uncracked RVE …………………………………………………… 99 6.1.1 Stress Distribution in Elastic Model ..………..……………………………………… 99 6.1.2 Stress Distribution in Elasto-Plastic Model ……………………………………… 100 6.2 Crack Patterns ………………………………………………………………………… 103 6.2.1 Results of Elastic Model …………………………………………………………… 103 6.2.2 Results of Elasto-Plastic Model ……………………………………………… 105 6.3 Relations of Crack Density with Number of Cycles ……………………………… 107 6.3.1 Results from Tanaka-Mura Equation ……………………………………………… 107 6.3.2 Results from Chan Equation ……………………………………………………… 111 6.4 Risk of Crack Initiation ………………………………………………………………… 113 6.5 Discussion ……………………………………………………………………………… 114 CHAPTER 7 CONCLUSION ……………………………………………………………… 117 APPENDIX A …………………………………………………………………………… 121 APPENDIX B …………………………………………………………………………… 126 APPENDIX C …………………………………………………………………………… 136 APPENDIX D …………………………………………………………………………… 141 REFERENCES ……………………………………………………………………………… 143 viii PREFACE It is found that most of the failure of engineering components is related to fatigue damage. To satisfy structural functions, the components have inevitably notches and/or holes, where the local stress level is higher than the average stress because of the stress concentration. As observed, some macro-cracks may form around these areas on or near the component surfaces after some loading cycles, even if the loading amplitude is much lower than the estimated safe load based on the static fracture analysis. Fatigue fracture may happen when the macro-crack has grown to a critical length and the remaining ligament cannot sustain the loading of the next cycle. In some cases, the length of the macro-crack was not long enough to be detected by common detecting devices when failure happened. To ensure the safety of engineering components the fatigue behavior of materials has received great attention. The fatigue behavior, however, is so complicated that greater efforts are still needed, especially in the regime of small cracks. It is found that the chemical composition, metallurgical phases, microstructure dimensions, processing and surface treatment can alter the fatigue behavior of small cracks significantly, not to mention the combined influence of temperature and environment media. For the most important fatigue stage, crack initiation, there is still no general law which can take these important factors into account. The present study is an attempt to find a quantitative method which is able to predict the fatigue life to crack initiation. This work is based on a mesoscopic model and focuses on the simulation of the multiple crack initiation behavior of a particular material. 1 This work is organized in seven chapters. In Chapter 1 some aspects about the fatigue behavior of metals, especially the crack initiation mechanisms and models related to the present study, will be reviewed and discussed. The simulated material and fatigue experiments carried out in a previous project will be introduced in Chapter 2. The ideas and hypotheses about the simulation work are explained in Chapter 3. The details about the construction of the two- dimensional and the three-dimensional models are described in Chapter 4. Chapter 5 is dedicated to the simulation results obtained from the two-dimensional models and Chapter 6 to the three- dimensional results. The conclusion is in the last chapter, Chapter 7. 2 CHAPTER 1 INTRODUCTION In this chapter, some important aspects of and developments in the crack initiation behavior of metals are reviewed and discussed. The background of fatigue research is briefly introduced in Section 1.1. In Section 1.2, microcrack initiation mechanisms are described and the influence factors related to the material microstructure are discussed. The existent prediction models of crack initiation are reviewed in Section 1.3. The available mesoscopic models and the model structures are summarized in Section 1.4. 1.1 Fatigue Behavior and Fatigue Tests 1.1.1 Three Stages of Fatigue In general, the fatigue process is considered to be composed of three stages: crack initiation, stable crack propagation and unstable propagation, which is followed by final fracture. The influencing factors on fatigue life Nf, the number of cycles to fracture, comprise the applied loading levels and types, loading history, material property, material processing history, surface treatment and also the service environment, such as temperature and the surrounding media. Fatigue life is spent mostly in the first two stages, i.e. crack initiation and propagation. The distinction of crack initiation and propagation is strongly linked to the size scale of the cracks concerned [1]. Technically the stage of crack initiation is originally referred to the period from uncracked material to the occurrence of detectable macro-cracks. It is possible, in practice, to distinguish the two stages quantitatively by the cracks measurable in experiments and during in- 3 service inspection. The size of the detectable crack is, normally, in the scale of millimeters. In this case, the period after crack initiation and before the final fracture is the propagation stage, which is, nowadays, called long crack propagation. The damage accumulation process under fatigue loading can be roughly divided into two different scenarios: Scenario A: A high number of microcracks initiate on the surface continuously during almost the whole fatigue life. Before the formation of macro-cracks, several different damage mechanisms exist simultaneously, such as microcrack nucleation, propagation and coalescence, or the combination and competition of these modes [2-4]. If a macro-crack is formed the fracture is imminent. This phenomenon is called multiple crack initiation behavior in literature. In this case the crack initiation life is comparatively long, sometimes up to 80% of the failure life. Therefore, this fatigue behavior is referred to as crack initiation dominated. The development of quantitative relations between models of the physical process of crack initiation and macroscale measurements of fatigue life is still at an early stage. Scenario B: Under certain conditions (for example when the loading level is low [2, 4]) only a few cracks initiate and then one primary crack propagates to a critical length. The final failure is caused by the primary crack propagation and the crack propagation life is relatively long. In this case the fatigue behavior is crack propagation dominated. For the problems that belong to Scenario B, the crack propagation model is a proper solution. Crack propagation is much better understood than crack initiation. According to linear elastic fracture mechanics (LEFM) the long crack propagation behavior can be described in a power law proposed by Paris mKCdNda )(/ ∆= (1- 1) 4 where a is the crack length and da/dN is the crack propagation rate per cycle. ∆K is the stress intensity factor range and its value depends on the applied loading, crack geometry and crack length. C and m are the material constants and can be obtained from crack propagation tests. In general, the complete da/dN curve is presented in a log-log diagram, as shown schematically in Fig. 1- 1 (a). It consists of three regions, I, II and III. Region I is the so-called near threshold region. When ∆K is lower than the stress intensity threshold, ∆Kth, the crack is supposed to stop growing. In region II the crack grows following the Paris law, which is a straight line in log-log coordinates. The crack propagates rapidly in the region III, leading to the final fracture. When the Paris law is used to deal with small crack propagation, it is found that the small crack behavior is quite different (if the size of a small crack is down to the scale of microstructure it is referred to as a microstructurally small crack). The crack propagation rate varies within a wider scatter band, as shown in Fig. 1- 1 (a) and the threshold of small crack is lower than that predicted by long crack experiment. In the last decades the small crack growth behavior has received intensive attention. From abundant experimental investigations and microscopy observations, it has been found that the abnormal behavior of small cracks is caused by the applied analysis method and by the nature of microstructures, as reviewed by Miller [5]. A small crack grows fast inside a grain but when it reaches the grain boundary the growth rate is slowed down and the crack possibly stops growing, as illustrated in Fig. 1- 1 (b) [6]. From the investigations of β Titanium alloys (bcc) [7-8], it was found that high angle grain boundaries stop microcrack propagating into the next grain, but low angle grain boundaries do not. The high angle grain boundaries become the barriers of microcrack growth. The small cracks show intermittent growth behavior. The grain misorientation plays an important role in this behavior. 5 (a) (b) Fig. 1- 1 Schemes of (a) crack propagation curve and (b) small crack growth behaviour Small crack Long crack ∆Kth II III I da /d N da /d N a ∆K 1.1.2 Fatigue Test-Strain Cycling and Stress Cycling The first study on fatigue test was made by Albert in 1829 with a device which applied cyclic loadings to a chain made of iron in order to find the number of cycles until fracture [1]. These kinds of experiments can be referred to as total-life fatigue tests and are still carried out nowadays for the study of the fatigue behavior of engineering components. With the development of theories about fatigue and fracture, more advanced test machines have been invented for more comprehensive fatigue tests. A wide spectrum of materials has been tested with different loading and environment conditions. The experimental methods most commonly used for investigation of the essential fatigue behavior of materials are rotary bending and tension-compression (push-pull) fatigue tests. With respect to loading conditions, the push-pull test can be divided into two types, stress cycling and strain cycling fatigue tests [9]. The data obtained from both of these tests are used for the fatigue resistant design of engineering structures. The so called stress cycling experiment, which is also referred to as high-cycle fatigue test, is analogous to the situation where the stress 6 level in components is much lower than yield stress. The strain cycling fatigue experiment, which is referred to as the low-cycle fatigue test, is more interesting for the purpose of fatigue life evaluation because the stress state in the specimen is more similar to that near the root of notches in the component, where the local stress level can be close to or even higher than the yield strength and the plastic deformation may occur. For both stress and strain cycling tests, the influence factors on fatigue life are: loading amplitude, loading ratio R (minimum to maximum amplitude), loading frequency f (or loading rate), temperature and environment media. The symmetrical push-pull loading, i.e. R is -1, is often applied. The fatigue behavior of material can differ with loading rate but when the loading rate is lower than a critical value (depending on the material), the fatigue behavior of the material is almost rate-independent. The fatigue behavior tested by stress cycling is different from that by strain cycling. In strain cycling fatigue, the strain amplitude is constant during the experiment. For most aluminum alloys and some types of steels, cyclic hardening behavior, i.e. the stress level increases with fatigue cycles, is often observed. If the stress level varies in the other way around, i.e. the stress level decreases with cycles, as observed in some hardened or strengthened materials, e.g. martensitic steel, cyclic softening happens [9-11]. 1.1.3 Damage Accumulation during Multiple Crack Initiation Multiple crack initiation has received significant attention recently. The availability of multiple sites for crack initiation makes it a common feature in many kinds of material failures, such as in the conventional fatigue for steels [12-14], Ni-based superalloy [15], α-irons [16-17] and cast aluminum alloy [18], as well as in thermal fatigue [19] and in the fatigue of welding [20]. From these studies, the following common characteristics are found: - The most dominant cracks were observed in larger grains. 7 - The crack initiation mechanism varies with temperature and chemical composition. For example in [16], at low temperature the initiation mechanism was intergranular initiation and at room temperature most cracks initiated were transgranular cracks. In [14], the crack growth in the ferrite phase was stopped by the pearlite phase and ferrite grain boundaries. - There are several types of cracks, one-segment cracks (i.e. the microcracks with no kinks) and multi-segment kinked cracks, observed on the surface. The process of damage accumulation is the combination of crack initiation, growth and coalescence. - In order to evaluate fatigue damage accumulation quantitatively, the crack density, i.e. the number of microcracks per unit surface area of the specimen, is introduced. The crack density varies with the normalized cycle N/Nf in a typical way, as shown in Fig. 1- 2[12]. This crack density increases at the beginning, reaches the maximum value and then starts to decrease, which indicates that crack coalescence happens. - In some combined conditions of loads and materials [21-23], the slip bands in early fatigue life become deeper and wider with the increasing number of cycles but no crack is detected. Then a microcrack appears after a comparatively short fatigue interval and grows up to the size of a whole grain. This implies that the damage is accumulated by some smaller defects, which are induced by dislocation motion, as will be discussed in the next section. In some materials, it is observed that the number of cycles for a microcrack growing along the slip band and reaching the first barrier is much smaller than the fatigue life, like a ‘sudden’ crack initiation [21]. 8 C ra ck d en si ty N/Nf Fig. 1- 2 A schematic drawing of the crack density varying with the normalized number of fatigue cycles 1.2 Mechanism of Crack Initiation Fatigue cracks are found initiating not only at the sites of inclusions, scratches or some other defects, but also from the well-polished surface of fine and uniform materials under fatigue loading, according to laboratory investigations. An early research about fatigue damage on an apparently defect-free surface was performed by Ewing and Humfrey [24]. In the experiments with Swedish iron (ferrite) subjected to rotary bending fatigue they found some slip bands on the surface. The slip was particularly intense along the slip bands. These slip bands were named ‘persistent slip bands’ (PSB) and crack initiated from these PSBs. In this section the mechanisms of PSB formation and crack initiation along PSBs will be described. 1.2.1 Mechanism of PSB Formation Since the phenomenon of persistent slip band (PSB) was discovered, many researches have been devoted to the investigation of how the PSB forms and the relation of PSB formation with fatigue loading and fatigue life. From the study on the behavior of single crystal fatigue (mostly fcc metals, such as copper and nickel), it has been found that the PSB formation is the result of the cyclic deformation of crystals. Some fundamental results were given by Mughrabi [25], who 9 studied the cyclic deformation behavior of a pure copper single crystal under fully reversed cyclic loading. It was found that the PSB is related to the amplitude of resolved plastic shear strain γpl. Fig. 1- 3 shows the peak value curve of cyclic resolved shear stress τs versus γpl for Cu single crystal oriented for single slip. The τs-γpl curve shows different characteristics in three regions, A, B and C. When the applied plastic shear strain γpl is lower than γpl,AB (in region A), the shear stress τs increases with γpl and approaches to a critical shear stress τ*s, saturation stress, from which a plateau of the curve starts. In region A, fine slip markings can be observed but there is no progressively accumulated damage. The persistent slip bands form at the beginning of region B, from where the plastic shear strain amplitude is larger than γpl,AB. The volume portion of PSB increases with the amplitude of γpl proportionally [26] in region B and PSBs will cover the whole grain when γpl approaches to γpl,BC. A similar mode of behavior has been found later in some other fcc and bcc single crystals and also in some polycrystalline materials [1]. In fcc single crystal oriented for multiple-slip, however, the hardening behavior is somewhat different. From the research of Gong et al. [27] it was found that the plateau of shear stress- plastic shear strain curve (range B in Fig. 1- 3) disappeared for Cu single crystal oriented for multiple-slip and the well-defined PSBs are not commonly found. This implies that the PSBs form dominantly in the single slip plane. Fig. 1- 3 Schematic of hear stress-plastic shear strain curve of fcc single crystal sτ * sτ ABpl ,γ plγBCpl ,γ A B C 10 There are fewer studies on the mechanism of PSB formation in bcc materials, although the PSB was first discovered in a low carbon α−iron. The fatigue behavior of pure bcc materials is very different from that of pure fcc materials. The strain hardening curves of a pure bcc single crystal show very strong temperature, strain rate and impurity dependence [9]. If the temperature is higher than a transition temperature Tk or the strain rate is low enough, or when impurity elements are added, even if the amount is very small, the cyclic deformation behavior of bcc changes significantly. Under these conditions, the cyclic deformation in bcc material can be quite similar to that of fcc material. In the research of Sommer [16] on the low cycle fatigue behavior of α-iron, with the added carbon content of only 74 wt ppm, the PSBs were observed, where the experimental temperature was above 343K and the strain rate was slower than 1×10-4 /s. The persistent slip bands were observed on the surfaces of a low carbon steel [28] and in the ferrite phase of a steel [29] as well. The PSB is a group of slip planes usually spreading in the whole grain size. The cyclic loading induced surface roughness, the extrusions and intrusions which were identified by scanning electron microscope, are located at the sites of PSBs. By advanced microscopic techniques, such as the atomic force microscope (AFM), the microstructure of persistent slip bands can be well-observed [30-31]. The profile of extrusions is approximately triangular and they grow during fatigue life in the direction of the active slip. Microhardness measurement on the PSB revealed that the PSB is softer than the matrix, therefore the material deformation is mostly carried by PSBs. 1.2.2 Mechanism of Crack Initiation from PSB The crack initiation in defect-free materials is found mostly taking place at PSBs. The locations of crack nucleation are reported at the PSB-matrix interface [32], at the root of intrusions [33] and extrusions [34]. The surface roughness is the result of irreversible dislocation 11 motion instigated by fatigue load cycles. From the observation of transmission electron microscopy (reviewed by Suresh[1]) it is revealed that a dipole consisting of edge dislocations of opposite signs will annihilate to form a vacancy if the space between them is smaller than a critical value. The annihilation of dipoles is responsible for the extrusions and the interstitial-type dipoles are responsible for intrusions. The crack nucleation mechanism proposed by Essmann et al. [35] gives a detailed microscopic description of the irreversible glide in PSBs based on the analysis of dislocation movement. The hypothesis is that the annihilation of vacancy-type dipoles is the dominant point- defect generation process and that the annihilation of dislocations within slip bands is the origin of irreversibility. This irreversible slip can happen when a screw dislocation glides along different paths forwards and backwards and consequently the PSBs are formed by the irreversible slip. The extrusion is the PSBs emerging on the surface. The cracks nucleate at the intersection of the PSBs and surface. The crack initiation mechanism proposed by Lin and Ito [36] and Tanaka and Mura [37] is also based on the formation of intrusions and extrusions and the irreversibility of dislocation motion. But here, the dislocation motion is described on the two adjacent parallel slip planes. The proposed mechanism has the background of experiment observations where it was found that the slip plane during the tensile loading and the one during compressive loading were closely spaced but, in fact, distinct from each other [1]. 1.2.3 Mechanism of Crack Initiation at Inclusions There are two typical damage modes regarding the crack initiation at inclusions: (i) the debonding of the inclusion from the matrix when the adhesion between inclusion and matrix is weak [38] and (ii) the breaking of the inclusion when the inclusion strength is lower than the matrix [39]. Matrix microcracks nucleate at the sites of the interfaces between the inclusion and 12 matrix. Crack initiating near inclusions can also be of the slip band type [38]. The size of inclusion is found to be a critical factor. For example, the work of Laz and Hillberry [39] on 2024-T3 aluminum alloy indicates that the size of cracked inclusions is larger than 5µm. Another phenomenon often observed in the material containing inclusions is the subsurface crack initiation when the inclusion size is large. Very small inclusions (<1 µm, for example) can suppress crack initiation through slip homogenization [9]. 1.3 Models of Crack Initiation 1.3.1 Conventional Models The early prediction models for crack initiation are based on the low-cycle fatigue experiment and the Manson-Coffin equation. From the fatigue experiment, the relation of loading level (stress range ∆σ in stress cycling and strain range ∆ε or plastic strain range ∆εp in strain or plastic strain cycling) and the number of cycles to specimen failure Nf can be obtained. The general form of the relation, found by Coffin and Manson e.g. for low-cycle fatigue, is in the form of a power function as the following equation, c ff p N )2( 2 'εε =∆ (1-2) where 2Nf is the number of load reversals to failure, the fatigue ductility coefficient and c the fatigue ductility exponent. and c are material parameters. Equation (1-2) is still widely used nowadays. However, the microstructural influence cannot be described by Eq. (1-2). One model proposed by Cheng and Laird [40] has a similar form: ' fε ' fε ' 2 CN f p =∆ αγ (1-3) 13 where ∆γp is the plastic shear strain range, C’ and α are material constants. Eq. (1-3) was developed on the basis of PSB formation but it does not provide any explicit microstructural parameter. 1.3.2 Microstructure-Based Models The life prediction model proposed by Tanaka and Mura [37] for the crack initiation from slip bands yields the relations between the number of cycles to crack initiation and material parameters. This model is based on the assumption that irreversible dislocation pile-ups cause crack initiation. In order to incorporate slip irreversibility, the deformation within slip bands is modeled by two closely adjacent layers of dislocation pile-ups. The dislocations in each layer have different signs, as shown in Fig. 1- 4. It is assumed that the motion of dislocations formed by previous forward loading in layer I are irreversible and that the reverse plastic flow is taken up by the motion of dislocations with the opposite sign on layer II. The positive back stress induced by the positive dislocation pile-ups in layer I facilitates the pile-up of the negative dislocations in layer II. This process continues with loading cycles. The progress of dislocation accumulation is calculated by using the theory of continuously distributed dislocations. The strain energy of dislocation is accumulated to the same extent in each forward and reverse loading. When the accumulated energy reaches the amount of fracture energy, a crack initiates. According to the Tanaka-Mura model, the quantitative equation to estimate the crack initiation life Nc for the stage I crack is derived as: 2)2()1( 8 kd GW N sc −∆−= τνπ (1-4) where G is the shear modulus and ν the Poisson’s ratio, Ws is the specific fracture energy for a unit area consisting of the surface energy and the plastic fracture work. ∆τ is the resolved shear stress range, which is the stress range from the minimum shear stress to the maximum shear 14 stress. k is the frictional stress and d is the length of the slip band, d = 2a. For the crack initiation life prediction of a material, d is the grain size. As parameters d, ∆τ, k and Ws in Eq. (1-4) are microstructure-related factors, the microstructure effects on crack initiation are introduced into the model. Fig. 1- 4 Tanaka-Mura model Since the predicted tendency of fatigue life varying with grain size coincides with what has been observed in experiments, there has recently been an increase in the application of the Tanaka-Mura equation. Hoshide and colleagues [4, 41] applied Eq. (1-4) to simulate the fatigue behavior for copper, steel and titanium alloys under multi-axial loading. Zimmermann and Rie [42] used this model for the simulation for aluminum alloy, iron and carbon steel under strain control fatigue loading. Alexandre et al. [43] applied the model to the analysis of the Inconel 718 alloy. Tryon and Cruse [21] applied the Tanaka-Mura model to develop a probabilistic evaluation of crack nucleation life. In order to predict the microcrack length at initiation and incorporate other microstructural sizes, a modified Tanaka-Mura equation is proposed by Chan [44], -a Slip band Interstitial dipole Vacancy dipole I II a y x Grain boundary 15 )()( )2)(1( 8 2 2 2 d c d h k GN ⋅−∆−= τνλπ (1- 5) where three more additional variables are introduced: c is half of the length or the depth of an incipient crack (the size can be a part of the slip band), h is the width of the slip band and parameter λ is a constant, λ = 0.005. This model is developed on the assumption that the dislocation dipoles contribute to the crack formation and the criterion of crack formation is seq dW γ2= where γs is the surface energy and Weq is the strain energy stored in the dislocation dipoles of a single slip band. For the convenience of description, Eq. (1-5) is rewritten as following: G d hc kd GNi ⋅⋅⋅−∆−= 2 2 )( )2()1( 8 λτνπ (1-6) The left term in the right-hand side of Eq. (1-6) is similar to the Tanaka-Mura equation (see Eq. (1-4)) but the specific fracture energy Ws is replaced by G d hcW ss 2)(⋅== λγ (1-7) In Eq.(1-7), the fracture energy Ws is considered to be only composed of surface energy γs. The above two models include the most influencing microstructural parameters for the life prediction of microcrack initiation. Most of the parameters can be determined by standard experiments and only a few are needed to be defined by additional investigation. The observation of crack nucleation on the surface by means of an atomic force microscope (AFM) can give more detailed information because of its high resolution. Based on this new technology, Harvey et al. [45] proposed a model to predict fatigue crack initiation life by means of these microscopic parameters: epys th i hEf KN εσ ∆ ∆= 4 2 (1- 8) 16 where ∆Kth is the long crack propagation threshold, σys is the yield strength of material, E is the Young’s modulus and ∆εp is the plastic strain range. The value of these parameters can be obtained with standard tests. he is the slip spacing and f is the fraction of plastic strain range of applied loading, which is related to the slip height δ. δ and he can be measured from the records of AFM photographs on the surface. It is supposed that a crack will initiate when the cumulative slip height reaches the threshold of crack-opening displacement. 1.3.3 Models Based on Probability Due to the pronounced influence of the microstructure, crack initiation is a stochastic process and is dealt with by initiation probability in some models [21, 47-51]. Some of them are derived from empirical equations based on the investigation of specimen surfaces [46-47] and the number of cracks is the random variable. Some other models combine the microstructure-based model with the stochastic model, such as the model proposed by Tryon and Cruse [21]. In their model the Tanaka-Mura model is employed for the evaluation of the number of cycles to crack initiation to a grain size. The model applied by Morris et al. [48] is to simulate crack initiation from inclusions. In this case the crack initiation life is a function of microstructure parameters [49], such as the inclusion size, the distance of the inclusion to the grain boundary and the effective stress. In the statistical model of Ihara et al. [50] the energy stored in the material of each cycle was the basic random variable. The criterion for crack initiation is the formation of a PSB when the accumulated energy is higher than a critical value. A stochastic model, recently developed by Meyer and Brückner-Foit [51], is focused on the influence of microstructure parameters on low-cycle fatigue life using a planar random cell structure. In this model the crack initiation probability depends not only on the strain amplitude, but also on the microstructure variables, such as individual grain size and orientation. 17 1.4 Modeling of Polycrystal Materials 1.4.1 Representative Volume Element In order to establish the macroscopic relation of mechanical and physical properties to the real material microconstituents and microstructures, the concept of representative volume element (RVE) is introduced [52]. An RVE is a material volume which is statistically representative of the infinitesimal material point and its neighborhood. The RVE can have many kinds of micro-elements, such as grains separated by grain boundaries, inclusions, voids, microcracks and other similar defects. It provides a mesoscopic analysis tool which links the macroscopic homogenous material and its inhomogeneous microstructure. The size of an RVE is macroscopically infinitesimal compared to the scale of the bulk material and its boundary conditions, so that the local stress state can be accurately represented. On the other hand, the RVE is microscopically large enough to represent the real material microconstituents and microstructures, and the micro-damage evolving process. As shown in Fig. 1- 5, L is the scale of bulk material, H is the scale of boundary conditions on bulk, D is the scale of the RVE and d the microstructure scale inside the RVE. The magnitudes of L, H, D and d satisfy the relations of D<. Under the given loading, the slip system with maximum Schmid factor will become active. But, as observed, the microcracks are only initiated along martensitic laths in the material F82H. In a study of a Fe-Ni-Co-Cr-Mo-C (carbon content 0.23%) alloy [69], it was shown that an aligned group of martensitic laths shares {110} slip planes, which lie along the axis of laths. This indicates that {110} planes are the only candidates where microcrack can initiate. 4.2.2 2D-RVE Model The global coordinate system is XYZ. The lattice coordinate system is fixed on each grain. The three crystallographic axes of bcc crystal are assigned to be the three axes 123 of lattice coordinate system, as shown in Fig. 4- 2. Grain orientation is represented by angle φ, which is the angle between 1-axis of lattice coordinate system and X-axis of the global coordinate system, as shown in Fig. 4- 3. Supposing )011( slip plane to be selected from the {110} family as the primary slip plane, the trace of )011( slip plane in a 2D-RVE is the straight line which is 45˚ to the grain 1-axis in 12- plane from the top view of the bcc lattice, as shown in Fig. 4- 2 and Fig. 4- 3. In this case the orientation angle of PCP α is equal to φ+45°. In cubic lattice, the slip plane (110) is perpendicular to )011( . Therefore the shear stress in this plane is the same as in (110). 45 Fig. 4- 2 Lattice coordinate system and crystallographic axes of bcc Fig. 4- 3 Orientation of )011( slip plane in two-dimensional model In the global coordinate system XYZ the 2D-RVE is represented by a square with unit thickness. The 12-plane of the lattice coordinate system is in the XY-plane of the global coordinate system. One example of a 2D-RVE with grain aggregate is shown in Fig. 4- 4. The 3-axis of all the lattice coordinate systems is in the direction of the Z-axis of the global coordinate system, which is perpendicular to the paper. The polygons in solid lines are theVoronoi cells and represent prior austenite grains. The grain orientations are determined by angles produced by a uniform random number generating process. The dash line in each grain represents the potential crack path (PCP) of the grain and indicates, also, the orientation of martensitic lath. The slip direction is in the a = c c a a 3 245° Viewing direction O 1 )011( plane X Y 1 2 O PCP O’ α 45° φ 46 Orientation of martensitic laths and PCPs Prior austenitic grains Grain boundaries O orientation of martensitic lath. It is known that a two-dimensional Voronoi tessellation contains more cells with a small aspect ratio than the real grain structure. The longest line may coincide with one of the edges of the random cell leading to unrealistic grain paths. Therefore, it was decided to define PCP by drawing a straight line through the centre of gravity of each individual grain. The PCP divides one real grain into two virtual grains, which are the substructures for finite element analyses. Fig. 4- 4 Illustration of a 2D-RVE with grain aggregate and PCPs Y X 4.2.3 3D-RVE Model The simplified 3D-RVE is the extension of the planar cell structure of a 2D-RVE in the third dimension. The cells are generated by the two-dimensional Voronoi process, as described in Section 4.2.2, and lie in XY-plane, as illustrated by bold solid lines in Fig. 4- 5. The 3D-RVE represents a thin layer of material on the specimen surface. It allows three-dimensional orientations 47 to be assigned to grains and constraint conditions to be applied in the third dimension. The plane which lies in z = 0 is called the base plane and is connected to the bulk material. The plane with z = thickness represents the free surface of the specimen, as shown in Fig. 4- 5. The fine lines in Fig. 4- 5 represent the potential crack paths in the grains. Fig. 4- 5 3D-RVE model Unlike the two-dimensional model, the “crack path” of the three-dimensional model is a plane which goes through the thickness of RVE. The PCP can be found by extending the (110) plane, which is the slip plane selected from the {110} family for the 3D-RVE, until it goes through the base plane and free surface. The PCP plane is determined by the intersection area of the (110) plane with XY-plane, as shown in Fig. 4- 6. In principle, the slip plane is tilted with respect to the global Z-axis direction because of the random grain orientation. However, the potential crack path is always supposed to be perpendicular to the XY-plane from the view of geometric dimensions. This simplification is reasonable since the thickness of the 3D-RVE is so small compared to the other two dimensions that the values of x and y coordinates for the traces of the same (110) plane Width Length Thickness Free surface O Y X Crack path Base plane Z 48 on the base plane and on free surface are almost the same. The approximation simplifies the model creation and the finite element analyses. This simplification, however, is not applied to the stress analyses on the plane of crack path. In other words, the stress components on PCP planes are presented in three-dimensional space and are associated with the grain orientations. Fig. 4- 6 Scheme of the potential crack paths of 3D-RVE 1 3 2 X Z Yo PCP plane (110) plane 4.2.4 RVE Size and Voronoi Boundary Effects Generally speaking, the number of grains in an RVE is not limited. The minimum number of grains required in a model depends on the following conditions. (i) The effective number of grains Ne in a Voronoi model is smaller than the given number of grains N. The phrase of ‘effective number’ means the number of grains which are not on model boundaries and ‘given number’ means the total number of grains in a model. The shapes of the grains on model boundaries are different from the grains inside the models. From the statistics of the numbers of Ne and N, it is found that the fraction of Ne /N decreases with N increasing, as shown in Fig. 4- 7. Because the area of grains on 49 model boundaries should smaller than that of inner grains, a model should consist of at least more than 40 grains. 20 30 40 50 20 40 60 80 100 120 Number of grains, N N e/N [% ] Fig. 4- 7 Fraction of the number of grains on RVE boundary to the number of grains in a model (ii) The size of RVE can be evaluated by the effective property, for example, the Young’s modulus, as reviewed in Chapter 1, Section 1.4. Because of the anisotropic property of an individual grain, the overall property of an RVE is much related to the number of grains and the grain anisotropic factor A’, 1211 442' CC CA −= (4-2) where C11, C12 and C44 are constants of material constitutive matrix, which will be described in Section 4.4. The factor A’ varies from 1 to 4 and A’=1 represents isotropic material. The larger the A’, the more grains are required. From the study on the Young’s modulus of an Al2O3 polycrystal with Voronoi models of 5~1000 grains [71] the standard deviation drops to 1.1% for the model with 40 grains. From this point of view, a model with 40 grains is large enough for the Al2O3 polycrystal. (iii) For the simulation of crack initiation, the number of grains should be large enough to yield the required crack density. The crack density Cd depends on 50 AN kCd = (4-3) where k is the number of cracks, N the number of grains in the model and A the average grain area. For example, if the required crack density is 4 mm-2 and the average grain area is 0.0025 mm2, then from Eq. (4-3) we get k/N = 0.01. Since k and N are integers N must be larger than 100 in this case. In order to simulate the continuous crack initiation process which can be compared to experimental data, a RVE with many grains is highly recommended. On the other hand, however, the model size should not be too large considering the FE meshing and the nonlinear material properties applied in the simulations, which will increase computation cost dramatically for large models. The final decision of the model size is a balance between the above mentioned influence factors. 4.3 Model for Finite Element Analyses 4.3.1 Coordinate Systems 4.3.1.1 Coordinate Systems of 2D-RVE Model The local rectangular coordinate system xiyizi, (the subscript i means the number of the corresponding virtual grain, i = 1,2, …, n, where n is the number of virtual grains) is located at the origin point of the global coordinate system XYZ, as shown in Fig. 4- 8. If the shear stress on the slip band )011( is to be studied, the local coordinate system is supposed to be defined as such: the xi-axis is in the direction of the ith PCP and the yi-axis is perpendicular to the xi-axis. In the two- dimensional model the zi-axis is always in the same direction of the Z-axis of the global coordinate system. In this case the orientation angle of PCP α is equal to φ+45°. 51 Fig. 4- 8 Schematic illustration of global, local and crystallographic axes coordinate systems on 2D-RVE X Global coordinate system Grain i Lattice coordinate system 123 Y Z zi xi yi 45° 1 2 Local coordinate system xiyizi Potential Crack path O αi 4.3.1.2 Coordinate Systems of 3D-RVE Model For simulating the randomly distributed grain orientations, a local coordinate system xiyizi is created in each grain. The axes of this local coordinate system coincide with the cubic lattice axes 123 of the grain as shown in Fig. 4- 9. The stress components from the FE analysis are given in the local coordinate systems. The relations between the global coordinate system XYZ and local coordinate systems xiyizi depend on the orientation of grains. It is well known that the slip direction of bcc crystal is <111>, which consists of the four slip directions in the (110) plane, namely [ ]111 , [ ]111 , [ ]111 and [ ]111 as shown in Fig. 4- 9 (from A to D), along the space diagonals AC or BD in the cubic lattice. Dislocation motion will occur along any of these two possible gliding lines when the magnitude of resolved shear stress is high enough to overcome the critical friction stress. 52 Fig. 4- 9 Slip plane (110) and gliding directions From the finite element analysis, the calculated stress tensors are given in the grain axes 123. The resolved shear stresses along these lines are obtained by transforming the stress components in xiyizi to the required plane and directions. For this purpose, two additional coordinate systems, (in ''' iii zyx [ ]111 and [ ]111 ) and (in '''''' iii zyx [ ]111 and [ ]111 ), are introduced. One of the a coordinate systems ''' iii zyx is shown in Fig. 4- 10. The ' ix -axis poi s in the gliding directions A The 'iy -axis is perpendicular to the ' ix -axis. Both axes 'ix and 'iy are on the (110) slip plane and the 'iz -axis points out to the normal direction of the slip plane and follows the right-hand rule. Similarly, the other coordinate system '''''' iii zyx can be defined but the '' ix points e direction of line BD. dditional nt C. in th Y Z xi 1 zi 3 [ ]111 [ ]111 [ ]111 D A C B O [ ]111 (110) yi 2 X 53 Fig. 4- 10 Coordinate system on slip plane along gliding direction''' iii zyx D C O z’i y’i x’i zi xi B yiA 4.3.2 Boundary Conditions The boundary conditions are defined as following: in order to simulate the strain control uniaxial tension in the specimen, nodal displacements are assigned to stretch the RVE to the maximum strain amplitude, 'H/2, in X direction. Since all the RVEs in the study include random size and orientation grain structure, the microstructure is not a periodic repetition of unit cells and not even a symmetric plane can be found in these models. Thus, neither periodic boundary condition nor symmetric boundary condition is used. The applied boundary conditions consist of a restrained boundary and a loaded boundary for the two-dimensional models, as shown in Fig. 4- 11. The displacements in X direction of nodes with coordinate x = 0 are restrained to zero. The displacements of nodes with coordinate x = xmax (xmax is equal to the model width) have the maximum value which produces equivalent strain amplitude on RVE in X direction. The other two boundaries y = 0 and y = ymax (ymax is equal to model length) are traction-free. There are two additional constraints for the three-dimensional models. The displacement of nodes on the base plane z = 0 in Z direction are restrained to zero because this plane is connected to the bulk material, 54 while the plane of z = zmax (zmax is equal to model thickness) is traction-free. The displacement in the Y direction of nodes at the lower-right corner x = xmax , y = 0 and the lower-left corner x = 0, y = 0, are constrained, as illustrated in Fig. 4- 12. Fig. 4- 11 Boundary conditions on 2D-RVE Fig. 4- 12 Boundary conditions on 3D-RVE Y O X X Y O Free surface Z 4.3.3 Element and Mesh The finite element model is created by means of the commercial software PATRAN [72]. Each virtual grain in the RVE is defined as an element set. The elements and nodes are generated by an automatic meshing process. The through-thickness crack will be introduced into the model by node 55 releasing for both 2D-RVE and 3D-RVE models. The modified RVE with the new crack is the model for the next step of simulation. 4.3.3.1 Element and Mesh of Two-Dimensional Model Whether the FE output value is sufficiently accurate or not depends greatly on the meshing policy, i.e. element type and meshing density. Which type of element should be chosen for stress analysis depends on the model geometry and boundary conditions. For Voronoi polygons the 4-node element was recommended according to Watanabe [73]. Because of the nature of polygon shapes, the 4-node rectangular element cannot be used. The ‘Paver’ method is chosen for its more consistent meshing in the surrounding of intersection points of curves. The elements generated with the ‘Paver’ method are mostly 4-node quadrilateral solid elements, but in the corner of a polygon the 3-node triangle solid element is often used instead. The ‘Paver’ meshing policy is designed in such way that the element dimensions on grain boundaries, PCPs and model boundaries are smaller in size [ ]111 than those of elements in the inner areas of grains, as shown in Fig. 4- 13. Because the stress gradient on the grain boundaries, PCPs and model boundaries are higher than in the inner areas of grains, a finer meshing is definitely necessary here. How fine the mesh density must be depends on the required accuracy of the analysis. In general, a finer meshing gives better convergence than coarser meshing. But if the mesh is too fine, the shape of ‘Paver’ elements may be distorted and consequently might cause the analysis accuracy to decrease. Furthermore, the number of degrees of freedom increases dramatically in a fine mesh and thus the running cost increases. Efforts should be made to getting a good combination of adequate accuracy and practical running time. To find the most suitable balance between these factors, a small model with 20 grains was tested with different meshes. The best meshing density was found after a few trials. With this meshing policy one virtual grain can be modeled by 200 (for 56 small grains) to 700 (for large grains) 4-node (and a few 3-node) linear elements. In the two- dimensional analysis, both plane-strain and plane-stress states are investigated. Fig. 4- 13 Mesh on RVE 4.3.3.2 Element and Mesh of Three-Dimensional Model The elements in the three-dimensional model are generated by extruding the two-dimensional planar elements in the third dimension to the required thickness. In the 3D-RVE 6-node solid prismatic elements and 8-node solid hexahedral elements are produced. The number of nodes of a three-dimensional model is twice that of a two-dimensional model if they have the same number of elements. When the size of the two-dimensional model approaches the performance limit of a computer system, the FE analysis for the corresponding three-dimensional model can not be carried out successfully. Under this circumstance the meshing density of the three-dimensional model has to be reduced. Therefore, the three-dimensional model in the simulation has fewer grains and a coarser mesh than those of two-dimensional models. 57 4.4 Material Properties 4.4.1 Stress-Strain Response of Elastic Material The stress state at a certain point is presented in a vector form in the present study as shown in Eq. (4-4a). The strain state at a point in a deformed body is presented in a vector form in Eq. (4- 4b). { } { }Txzyzxyzyx τττσσσσ = (4-4a) { } { }Txzyzxyzyx γγγεεεε = (4-4b) The generalized Hooke’s law gives the elastic stress-strain relation as [74]: { } { }εσ ijC= where Cij is the constitutive matrix or stiffness matrix. It comes from the fourth order tensor Cijkl, which is called constitutive tensor and has 81 material constants. Because of the symmetry of stress and strain tensors, the number of material constants in Cijkl is reduced to 36 and they are represented in the form of 6×6 stiffness matrix Cij. From the strain energy density theory it can be proved that the stiffness matrix Cij is symmetric to its diagonal line, i.e. Cij = Cji. Hence there exist 21 independent elastic constants in Cij for an anisotropic material. If the material has three planes of elastic symmetry the independent constants are reduced to 9 and the matrix Cij becomes: (4-5) ⎥⎢ 000C ⎥⎥ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎣ ⎡ = 66 55 44 33 2322 131211 0 00 000 000 C Csym C CC CCC Cij 58 The material with three symmetric planes is called orthotropic material, which displays different values of stiffness in mutually perpendicular directions. If the crystal has more symmetric planes, there are fewer constants in matrix Cij. There are three independent constants in a cubic crystal. In this study the linear elastic stress-strain response of individual grains of the investigated material is assumed to be orthotropic. The components of Eq. (4-5) are chosen from the material stiffness matrix of a single crystal pure iron [75]: C11 = C22 = C33 = 233 GPa, C12  = C13 = C23 = 135 GPa, C44 = C55 = C66 = 118 GPa. In the FE analysis code ABAQUS, the local coordinate systems are designed to define the material orientation and the output stress (strain) components. When the local coordinate system is defined at the slip direction (i.e. not at the lattice coordinate system), for example, a transformation of the stiffness matrix is imposed. The transformation matrix between two coordinate systems is described in Appendix A1. 4.4.2 Stress-Strain Response of Elasto-Plastic Material Experiments showed non-linear stress-strain curves with pronounced plastic deformation for most strain ranges for F82H steel. Therefore, a purely elastic simulation does not seem to be appropriate. An elasto-plastic response is selected for the present simulations. In order to distinguish between the contributions of elastic and plastic deformation, the elasto-plastic material behaviour is approximately described by a piecewise linear stress-strain curve. In many cases, the elasto-plastic response of ductile metals can be simplified into three parts, one elastic part, one hardening plasticity part and a perfect plasticity part. One typical example is shown schematically in Fig. 4-14. The bold solid line represents the stress-strain curve of a material obtained from a tension test. E is the Young’s modulus, E’ is the tangential modulus of the linear hardening plasticity, σy0 is the stress where plastic deformation appears and σu is the stress where perfect plasticity starts. 59 In the present study the cyclic stress-plastic strain relation, i.e. the relation of stress amplitude ∆σ/2 and plastic-strain amplitude ∆εp/2, for the plastic deformation part is derived from experimental data by E p 222 σεε ∆−∆=∆ (4-6) The experimental data can be found in Table 2-4 (Chapter 2, Subsection 2.3.2). Fig. 4-14 A scheme shows a typical tri-linear stress-strain curve The hardening plasticity region (part II in Fig. 4-14) is fitted approximately by a four-piece curve, as shown in Fig. 4-15, where the solid lines are the fitting curves and the experimental data are shown as solid triangles. The initial yield point is obtained from the plasticity part of the stress- strain curve at vanishing plastic strain. The perfect plasticity starts from the ultimate tensile strength Rm (540MPa) where the plastic strain is estimated as 10.0%. The applied stress-plastic strain data are listed in Table 4-1. In order to deal with the combined stress states, the elastic limit is presented by a yielding criterion. The yielding criterion used in the present elasto-plastic model is the von Mises criterion, 1 E E´ 1 I II III ε σ εeO εp σu σy0 60 which states that yielding will happen when the maximum shear strain energy at a point in the material reaches a critical value. Since the shear strain energy is proportional to the second invariant of the deviatoric stress tensor J2, the criterion can be expressed as: 0)( 22 =−= κσ Jf (4-7) whereκis a critical value of yielding. J2 can be expressed in term of stress components in the following equation, ( ) ( ) ( )[ ] 2222222 61 zxyzxyxzzyyxJ τττσσσσσσ +++−+−+−= Fig. 4-15 Experimental data of stress amplitude vs. plastic-strain amplitude and fittings Table 4-1 Stress-plastic strain data for simulation 0 100 200 300 400 500 600 0 0.005 0.01 0.015 ∆εp/2 [-] ∆σ /2 [M P a] 2 43 1 Stress [MPa] Plastic strain [-] 240 0.0 392 0.001 450 0.002 490 0.01 545 0.1 61 4.5 Modeling of Crack Initiation Process 4.5.1 Fatigue Model One of the applied models in the simulation is the Tanaka-Mura equation Eq. (1-4) (see Chapter 1 [37]). On the basis of this model, the number of cycles to crack initiation in the ith individual grain ∆Ni can be estimated by Eq. (4-8), which is the rewritten form of Eq. (1-4).: 2)2()1( 8 cresi s i i d GWN ττνπ −∆−=∆ i = 1, 2,… n (4-8) where, n is the number of prior austenite grains, di is the length of the ith slip band and ires τ∆ is the resolved shear stress range on the ith slip band. The crack initiation life Nk is the sum, ∑ = ∆= k j j ik NN 1 k = 1, 2, … l (4-9) where l is the number of simulation loops. The corresponding crack density , i.e. the number of cracks in the unit area, is defined as, kd C A kC kd = k = 1, 2, … l (4-10) where A is the area of the model. Similarly, the number of cycles to crack initiation can be estimated by the Chan equation )()( )2)(1( 8 2 2 2 i i i i cres i d c d hGN i ⋅−∆−=∆ ττνλπ (4-11) where h is the width of the slip band, c is half of the length of a nucleated crack and λ is a constant λ = 0.005. In the present study c = d/2. 4.5.2 Average Resolved Shear Stress 4.5.2.1 Transformation of Stress Tensors 62 The stress components at nodes on RVE are calculated by ABAQUS. If the local coordinate system does not coincide with the crack path, the stress components from FE output need to be transformed. The stress components in the two coordinate systems iii zyx iii zyx and can be presented as, ''' iii zyx { } { } { } { }Txzyzxyzzyyxx T xzyzxyzzyyxx ''''''' τττσσσσ τττσσσσ = = If the stress components in the coordinate system are known, then the stress in the coordinate system can be calculated by a stress transformation from iii zyx ''' iii zyx { } [ ]{ }σσ T=' (4-12) where [T] is the transform matrix. More details about tensor transformation can be found in [76] and in Appendix A2. In 2D-RVE, the shear stress xy'τ is the required component. In 3D-RVE the stress transformation Eq. (4-12) is applied to the two slip systems and the two sets of shear stress ''' zxτ and ''' ′′′′ zxτ along the two different slip lines are computed for the subsequent analysis. 4.5.2.2 Average Resolved Shear Stress Since the amplitude of the resolved shear stress on slip bands is inhomogeneous in the simulation models, the average shear stress iτ on the ith PCP is taken as the resolved stress, ∫= iA i i i dAzyxA ),,(1 ττ i = 1, 2, … n (4- 13) where Ai is the area of the ith slip band and τi(x,y,z) is the shear stress distribution function on the slip band. 63 A. Average Resolved Shear Stress of the Two-Dimensional Model In a two-dimensional model, because of the unit thickness of the model, Eq. (4-13) becomes ∫= iL i i i dLyxL ),(1 ττ (4-14) where Li is the length of the slip band on the ith grain. The average stress iτ can be derived from a discrete form of Eq. (4-14) as follows ∑− = + ∆+= 1 1 1 ) 2 (1 im j ij ijij i i LL τττ ∑− = ∆= 1 1 im j iji LL where Li is the length of the crack path in the model, and ∆Lij the distance between two adjacent nodes along the path on the ith PCP. The notations of τij and τij+1 represent the shear stresses at two adjacent nodes, where j is an index and mi is the number of nodes in the ith path. The shear stress component τij obtained from the output of FE analysis corresponds to the applied strain amplitude in the tension part. The Tanaka-Mura equation is based on the assumption that the compression amplitude results in the same extent of dislocation motion but in the reverse direction. Thus the ires τ∆ in Eq. (4-8) and Eq. (4-11) refers to the whole stress range of fatigue. Therefore, the total resolved shear stress range ires τ∆ on a PCP is equal to 2 iτ . B. Average Resolved Shear Stress of the Three-Dimensional Model The discrete form of Eq. (4-13) for the simplified uniform-thickness three-dimensional model is ∑− = ∆= 1 1 1 im j ijij i i AA ττ (4-15) 64 where ijτ is the average stress on the slip plane of the jth three-dimensional element on the ith PCP. This plane is denoted as the P-plane of the element and it is located on the PCP plane, as shown in Fig. 4- 16. Since there is only one layer of elements, the stress distribution on this plane is linear and the average stress ijτ can be derived, ) 22 ( 2 1 11 base ij free ij base ij free ij ij ++ +++= τττττ where free ijτ and freeij 1+τ are the shear stresses at two adjacent nodes on the free surface of the ith P and base ij CP τ and baseij 1+τ on the base plane, as illustrated in the scheme of Fig. 4- 16. The area of a PCP plane Ai is the sum of the areas of element P-planes ∆Aij, ∑− = ∆= 1 1 m j iji AA ∆Aij is surrounded by four element edges, freeijL∆ , baseijL∆ , hij and hij+1, as shown in Fig. 4- 16. In the simplified model, = = free ijL∆ baseijL∆ ijL∆ , hij = hij+1 = h, where h is the thickness of the model. Therefore, hLA ijij ⋅∆=∆ Fig. 4- 16 Elements on PCP plane Free surface PCP Plane free ijτ base ij 1+τ base ijτ free ij 1+τ hij+1 hij free ijL∆ base ijL∆ P 65 It is assumed that the deformation of an element is negligible compared to the element nominal size. Then Eq. (4-15) has the following form: (4-16) There are two slip directions on one PCP plane in the three-dimensional RVE. The average shear stresses on both of the two directions are derived by Eq. (4-16). The one with the higher value among these two average shear stresses leads to the lowest value of the number of cycles to crack initiation and hence determines the onset of the damage accumulation process in this particular direction. ij m j free ijijijij⎨ ⎟⎟⎜⎜+⎟⎟= ∑ freebasebase i i LL i ∆ ⎭⎬ ⎫ ⎩ ⎧ ⎠ ⎞ ⎝ ⎛ + ⎠ ⎞ ⎜⎜⎝ ⎛ +− = ++1 1 11 222 11 τττττ 4.5.3 Crack Initiation Process The number of cycles to crack initiation ∆Ni is calculated for all potential crack paths in the grain structure. The first micro-crack is initiated along the crack path with the shortest life, and is introduced into the RVE model. Cracks too close to the RVE boundary have to be excluded for stability reasons. In the next step a stress analysis is performed for the RVE with one crack. Compared to the undamaged RVE, the average stress level drops in a displacement-controlled situation, and the stress is redistributed in the vicinity of the cracked grain. Based on this new stress field, the number of cycles to crack initiation is again calculated for all remaining crack paths. As in the case of the undamaged structure, the crack path with the minimum value of the number of cycles can be identified. Then the next crack is introduced into the RVE. The simulation is stopped if the crack density reaches a critical value or if the FE model becomes unstable. 66 It should be noted that a crack can initiate only when the resolved shear stress range ires τ∆ is higher than 2τc. If ires τ∆ is lower than 2τc no dislocation pile-up will occur. However the same value of ∆Ni is obtained from Eq. (4-8) or Eq. (4-11) for negative or positive values of ires τ∆ - 2τc. Negative values may occur if the local shear stress is low. On the other hand, if ires τ∆ is equal to 2τc by chance, the number of cycles ∆Ni may be an infinitively large number. In order to eliminate the possibility of incorrect results the number of cycles is set to a large value, 107, i.e. the fatigue limit, once either of these two cases occurs. 4.5.4 Summary of Simulation Procedures and Applied Criteria All procedures of simulation for the two- and the three-dimensional models are summarized in the following steps: Step 1 : Create an uncracked RVE with grain structure, PCP and material properties; Step 2 : Create FE models with mesh and boundary conditions; Step 3 : Run FE analysis; Step 4 : Get output of stress components from FE data file and calculate the average shear stress on all crack paths; Step 5 : The number of cycles on all PCPs are calculated from Eq. (4-8) or Eq. (4-11); Step 6 : The PCP with the minimum value of ∆Ni becomes a crack; Step 7 : This crack is introduced into RVE and the RVE structure is modified; Step 8 : Go to Step 3. The above steps will be repeated until one of the following events occurs: • The number of load cycles for the given strain amplitude reaches the preset number, which is estimated to be more than 20% of the failure life of experimental data; • Half of the PCPs are cracked (except the ones near to RVE edges); 67 • The number of cycles to crack initiation of all remaining PCPs are larger than 107; • The model becomes unstable. The next two steps deal with the data processing of simulation results: Step 9 : The crack density versus the number of cycles is derived from Eq. (4-9) and Eq. (4- 10), respectively; Step 10 : The stress distribution is obtained by using the PATRAN post-process from the output file of the finite element program. 4.6 Verification of Simulation Model In order to clarify whether the models used in the present simulation can suitably represent the studied material or not, a validity check is carried out. The following methods are applied: o Statistics of the geometrical parameters of the model microstructure, such as the average grain size, the distribution of grain size and the distribution of martensitic lath orientation; o FE analyses for the stress-strain response of the RVE. The size of an RVE is represented by the number of grains contained in the model. For a two- dimensional model it represents the area corresponding to the real material. For example, a 100 grain 2D-RVE represents a surface area of 0.27 mm2, which is considered to be large enough to grasp the main features of the crack initiation process. The simplified three-dimensional model is created by extruding the two-dimensional model in the third dimension. The area of 2D-RVE represents the same area of three-dimensional model when they have the same number of grains. The applied loading, i.e. the strain amplitude ∆ε/2, varies from 0.25% to 0.38%, which represents the strain range ∆ε from 0.50 to 0.76%. 68 The effects of the model structure are studied using several different 2D-RVEs. A short name is assigned to each model, as listed in Table 4-2. The model 3D200 in Table 4-2 is created for the analysis of overall stress-strain response of the simulation model. The material parameters are the same for all these models. Table 4-2 Characteristics of RVE models Short name of model Number of grains Represented length [µm] Represented area [mm2] 2D100_1 100 520 0.27 2D100_2 100 520 0.27 2D100_8 100 520 0.27 2D80 80 465 0.216 3D80 80 465 0.216 3D200 200 735 0.54 The uncracked 2D100_1 model, along with the FE mesh, is shown in Fig. 4- 17. The grain structures are represented by dark gray lines, the mesh presented by light gray lines and the potential crack paths by bold black lines. The angle between crack path and the lattice coordinate 1-axis is 45º. From the automatic meshing of PATRAN pre-processing about 85000 elements and 87000 nodes are generated. The chosen meshing density is considered to be fine enough to get accurate stress analyses. 69 Fig. 4- 17 RVE 2D100_1 with FE mesh Since the number of degrees of freedom of a three-dimensional model is higher than that of a two-dimensional model, the model size of the 3D-RVE is reduced to 80 grains. In order to make a direct comparison between the results of two-dimensional and three-dimensional models, an isostructural 2D-RVE is created, which has the same grain structure and similar grain orientation as the 3D-RVE. The 3D-RVE is generated with an alternative Voronoi parameter. Therefore, the grain configuration is different from the two-dimensional models described above. The finite element mesh of the 3D80 model is shown in Fig. 4- 18 where the grain structures are displayed with solid lines and the crack paths with dash lines. About 39,000 three-dimensional solid elements and 80,000 nodes are generated in 3D80 by the PATRAN automatic meshing procedures. The number of elements in RVE 2D80 is the same as 3D80 but only 40,000 nodes are generated. 70 Fig. 4- 18 3D80 RVE with grain structure and FE mesh 4.6.1 Similarity of Mosaic Model to Studied Material The geometrical similarity of the Voronoi tessellation to the real grain structure of studied material is checked by a comparison based on the statistics in terms of grain size and orientation. The grain size and orientation of the martensitic laths of the studied steel are available in the database. The statistical distribution of the grain sizes of Voronoi models is obtained from the measurement of a certain amount of Voronoi cells with the same counting methods as for the material F82H. 4.6.1.1 Structure of 2D-RVE The grain size distribution of 2D-RVE models is obtained from the statistics of 326 Voronoi cells. The relative frequency and the probability distribution, from the statistical data of the F82H grains and the Voronoi cells, are shown in Fig. 4- 19 (a) and (b) respectively. The distribution of martensitic lath orientation in models is obtained by the statistics on the angles generated from the 71 random number generator. The histogram of relative frequency of orientation of 2D-RVEs and material F82H are shown in Fig. 4- 20. It can be found that, as shown in Fig. 4- 19 (a), the Voronoi cells have a larger portion of grain size from 0.04 to 0.08 mm than that of the real material, F82H steel. Therefore, the probability distribution of cells increases faster in the interval of 0.04~0.08 mm, see Fig. 4- 19 (b). The portion of cell size of the Voronoi models, in the range of either larger than 0.08 mm or smaller than 0.04 mm, is smaller than that of F82H. The average size of the Voronoi model is 51.3 mm for 2D-RVEs. 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 Grain size [mm] R el at iv e fre qu en cy 2D model F82H (a) 0.00 0.20 0.40 0.60 0.80 1.00 0 0.03 0.06 0.09 0.12 Grain size [mm] Pr ob ab ili ty d is tri bu tio n 2D model F82H (b) Fig. 4- 19 Distribution of grain size, from the F82H steel and simulation models 72 The orientation distribution of the martensitic lath from the material F82H and from the two- dimensional models is, in general, very similar, as shown in Fig. 4- 20, with the only exception of the interval 50º ~ 60º, where a few more laths are orientated in real material than those in models. 0 0.05 0.1 0.15 0.2 10 20 30 40 50 60 70 80 90 Orientation of martensite lath [°] R el at iv e fre qu en cy Model F82H Fig. 4- 20 Orientation distribution of the martensitic lath from F82H steel and models 4.6.1.2 Structure of 3D-RVE From the statistics in term of cell size of the three-dimensional Voronoi models (on 199 cells), the relative frequency of cell size can be obtained. The histogram in Fig. 4- 21 shows the cell size distribution of 3D-RVE models together with the distribution of grain size of the real material. Similar to the phenomenon observed in 2D-RVEs, there are more grains with intermediate grain sizes between 0.04 mm to 0.07 mm in 3D-RVEs than that in F82H material. The average cell size of Voronoi tessellation is 53.6 mm. 73 0.00 0.10 0.20 0.30 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 Grain size [mm] R el at iv e fre qu en cy 3D model F82H Fig. 4- 21 Grain size distribution in 80 grains RVE 3D80 4.6.2 Stress-Strain Response of Models The idea is to apply a uniaxial tension strain on the uncracked RVE model with the same boundary conditions as described in Chapter 4, Section 4.3.2. The stress distribution σ(x,y,z) is obtained from the output of elasto-plastic FE analysis and is represented by the element stress σι over the model, i = 1, 2, … m, where m is the number of elements of the model. The average stress σ of the model is derived from the stress distribution output according to the following equation, ∫ ∑ = == V m i iiVV dVzyx V 1 1),,(1 σσσ where V is the model volume and Vi is the element volume. In a stochastic mosaic model the average σ will be approaching the σmatl, which is the simulated macroscopic stress of the material, when the model is sufficiently large. The difference between the macroscopic stress σmatl and the stress from experimental data, σtest, indicates the similarity of the mechanical response of the model to the simulated material. 74 Based on this theory, a three-dimensional, elasto-plastic model with 200 grains, 3D200, is constructed and the obtained results are shown in Fig. 4- 22. Comparing the results of the model with the experimental data of the F82H material, one concludes that the mechanical behaviour of the 200-grain model is quite similar to the real material. However the average stress σ from the two-dimensional models with lesser number of grains, such as the 100-grain models, may scatter more or less around the average stress of the 3D200 model. The fluctuated magnitude may depend on the microstructure of each particular model. In general, the designed model suitably analogues the material behaviour. Fig. 4- 22 Stress-strain response from the simulation model and experimental data 200 300 400 500 0.1 0.2 0.3 0.4 0.5 ∆ε/2 [%] ∆σ /2 [M Pa ] Average stress Test data 75 CHAPTER 5 RESULTS OF 2D SIMULATIONS The simulation results presented in this chapter are obtained from the following four kinds of two-dimensional models: (1) plane-strain elastic model; (2) plane-stress elastic model; (3) plane- strain elasto-plastic model and (4) plane-stress elasto-plastic model. The simulation is carried out in two steps: at first tentatively assuming that the parameters in [41] are suitable to the studied material and secondly correcting the simulation by using the optimum parameter(s) which will be determined in the parameter study. The stress distributions obtained from finite element (FE) analyses for the uncracked models are shown in Section 5.1. The crack density, i.e. the crack number per unit area, is presented as the function of the number of cycles. The results of crack density varying with the number of cycles obtained from simulation are presented in Section 5.2. The method and results of parameter studies are also given in this section. In Section 5.3 the stress redistribution caused by the introduced crack in the RVE and its influence on crack initiation sequence are exhibited. In Section 5.4 crack patterns from elasto-plastic simulation are presented. To achieve a fast and efficient simulation, a program package is developed in the present research, based on programs compiled at the University of Karlsruhe in a previous project. The program package consists of some functional subroutines and even the ABAQUS code. It can 76 create FE models, perform FE analyses, process data obtained from ABAQUS and predict the number of cycles to crack initiation, continuously and automatically. By running these programs a great number of simulations can be carried out in an acceptable time period. The description of these in-house programmes can be found in Appendix B. All the images of stress distribution in RVE models, with or without crack, are obtained by PATRAN post-processing from the output of finite element analysis. The unit of the stress scale for all the images is pascal (Pa). Some of these images display stress distribution and nodal displacement simultaneously. In order to make the introduced cracks more distinguishable, the magnitudes of nodal displacement in some images are magnified by a factor. The loading direction of simulation models is assigned along the horizontal axis for the convenience of PATRAN processing. It should be noted that in the experimental scans the loading axis is in the vertical direction. 5.1 Stress Distribution in Uncracked RVE The von Mises stress and shear stress for the four models are calculated. According to the images processed by PATRAN, the stress distribution contours for the four kind models are quite similar although the stress magnitudes are different. Thus, only results from two plane-strain models, model (1) and (3) will be presented in this section. 5.1.1 Stress Distribution in Elastic Models 5.1.1.1 Von Mises Stress Distribution Fig. 5- 1 (a) shows the image of von Mises stress distribution from results of FE analysis (plane-strain) with the orthotropic material response in the 2D100_1 subjected to a strain range ∆ε = 0.60%. The yellow lines represent the grain boundaries and the white lines the crack paths. The inhomogeneous stress distribution caused by grain misorientation is noticeable. The von 77 Mises stress varies within each grain and from grain to grain. The difference between the highest and the lowest value is about a factor of three. A very impressive feature is the high stress concentration at triple-points and along the neighbouring grain boundaries, as pointed out by arrows in Fig. 5- 1 (a). The stress concentration at triple-points is a phenomenon often observed in experiments. One example is chosen in the area inside frame A of Fig. 5- 1 (a), which is magnified and shown in Fig. 5- 1 (b). The highest stress gradient occurs on the high-angle boundary of two adjacent grains, of which the orientation difference is 42°. In general, the high stress gradient areas are found mostly near grain boundaries. The stress level and distribution within one grain is related to the orientation of the grain and the misorientation of its neighbour grains. Some typical examples can be found in Fig. 5- 1 (c), which is the magnified area within frame B in Fig. 5- 1 (a). The numbers in the blue boxes are grain indexes. The grains numbered 12, 45 and 70 are oriented in the similar directions, which can be recognized from their crack path orientations. The stress levels of these three grains are similar, in the range from 400 MPa to 470 MPa. The grains numbered 61, 68 and 93, however, undergo quite different stress levels, which vary from 340 MPa to 600 MPa in the elastic model, although these grains are also similarly oriented. The highest stress levels occur in those grains, where they are surrounded by the grains with the misorientation angle of about 45º. Two groups of grains are selected for illustrating, for example, in Fig. 5- 1 (c) the grains numbered 58, 87 and 10 and the grains numbered 12, 70, 6, 45 and 90. The crack paths of grains numbered 58, 87 and 10 are oriented at about 0º or 90º. (Because of the symmetry of crystal bcc, the two grains with an angle difference of 90º have the same orientations in the two-dimensional models.) The orientations of crack path of grains numbered 12, 70, 6, 45 and 90 are all around ±45º to the loading axis. Grains 10, 58 and 87 are surrounded by grains 12, 70, 6, 45 and 90. The von Mises 78 stress levels of the grains numbered 10, 58 and 87 (oriented ±45º to loading axis) are the highest in the model. The stress distribution contours under differently applied strain ranges are very similar in the elastic analyses, but the stress amplitude is correspondingly higher if the applied strain range is higher. B A (a) Von Mises stress (Pa) Loading direction (b) Local area in square A of (a) (c) Local area in square B of (a) Fig. 5- 1 Von Mises Stress distribution, 2D100_1, ∆ε = 0.60%, elastic model, PE 93 61 19 44 68 6 45 70 1058 90 12 87 Grain boundary 79 5.1.1.2 Shear Stress Distribution Fig. 5- 2 shows the local shear stress distribution in 2D100_1 when applied strain ∆ε is  0.60%. The numbers in the blue boxes are grain indexes. In the two-dimensional FE analysis, the output stress components are in the local coordinate system xiyizi of each grain. As has been defined in Chapter 4, the xi-axis of local coordinate system is in the direction of the crack path, which is the slip plane of the grain in the two-dimensional model. Hence Fig. 5- 2 displays the distribution of local shear stress τxy in grain slip plane. It is found that the local shear stresses are fairly constant within one grain, but typical values differ by a factor of ten or even more from grain to grain. The shear stress magnitude strongly depends on the path orientation. If the orientation angle of a crack path is in the maximum shear stress direction, i.e. about ±45º to the loading direction, the shear stress in it is comparatively high, such as in grains 23, 86, 3 and 94 in Fig. 5- 2. The shear stress magnitude is also related to the misorientation angles between neighbour grains. For instance, the two grains, numbered 85 and 97, have nearly the same orientations but different shear stress levels because they are differently influenced by neighbour grains. 85 97 3 23 86 94 Fig. 5- 2 Shear stress distribution, 2D100_1, ∆ε = 0.60%, elastic model, PE 80 5.1.2 Stress Distribution in Elasto-Plastic Model The obtained von Mises stress distribution in RVE 2D100_1 from elasto-plastic, plane-strain FE analysis is shown in Fig. 5- 3 (a), where ∆ε = 0.60%. Fig. 5- 3 (b) shows the local shear stress distribution along crack paths. Comparing Fig. 5- 1 (a) and Fig. 5- 3 (a), Fig. 5- 2 and Fig. 5- 3 (b), one can find that the von Mises stress contours and shear stress contours in the elastic model and the elasto-plastic model are quite similar. The stress concentration at grain triple-points in the elasto-plastic model caused by grain misorientation is very similar to that in elastic model. But the stress magnitude in the elasto-plastic model, either von Mises stress or shear stress, is much lower than that in the elastic model. The difference between the maximum and the minimum von Mises stress is about a factor of two in the elasto-plastic model. It indicates that the stress fluctuation in the elastic-plastic model is not as pronounced as that in the elastic model. The plastic strain magnitudes ∆εp/2 of the specimens are in the range from 0.0684% to 0.1775% corresponding to the applied strain amplitudes ∆ε/2 of 0.25% ~ 0.38%, as shown in Table 2-4 (Chapter 2, Section 2.3). The tangential modulus E’ in the plastic deformation regime drops below 150 GPa, which is much lower than Young’s modulus of F82H, 202 GPa, as shown in Fig. 4-15 (Chapter 4, Section 4.4). Under a given displacement boundary condition, the stress magnitude determined by elasto-plastic model, therefore, is much lower than that by the elastic model, especially in the area where significant stress concentration occurs. The simulation with elasto-plastic material properties is considered to be a better way since the stress state is more similar to that in a real material under the applied strain amplitudes. 81 (a) Von Mises stress (b) Local plastic shear strain Fig. 5- 3 Stress distribution, ∆ε = 0.60%, 2D100_1, elasto-plastic model, PE 5.2 Relations of Crack Density Versus Number of Cycles The simulation yields the crack density Cd, i.e. the number of cracks per unit area, as the function of the number of load cycles N, when the number of load cycles to crack initiation ∆N is determined. The elastic model is not a suitable one for the stress analysis when a crack is introduced into the model because it is not able to simulate the stress singularity at crack tips. The stress magnitude determined by the elastic model is unreasonable high at crack tips, as will be shown in the next section. Therefore only results from elasto-plastic models, under plane-strain and plane-stress, are used for the quantitative analysis in this section. 5.2.1 Tentative Parameters In elastic models, parameters G, Ws and τc in Eq. (4-8) are material constants. In the first step of research, the tentative parameters from [41] for plain carbon steel, as listed in Table 5-1, were applied. In the elasto-plastic models, the shear modulus G in Eq. (4-8) is no longer a constant in the plastic deformation region. When the applied loading exceeds yield strength the shear 82 modulus G’ drops continuously in a fashion as the decreasing of tangential modulus E’. Therefore, in the first step of research, a simple two-linear curve of G’ was used. The resolved stress level ∆τ in the slip band and the length of slip band d are variables. The latter only depends on the model structure but the former depends on the model structure and the global stress state of the model. Table 5-1 Material constants in Eq. (4-8) [41] Shear modulus of elasticity G [GPa] Poisson ratio ν Specific fracture energy Ws [kJ/m2] Critical shear stress τc [MPa] 81 0.3 2.0 108 The simulation results with the tentative parameters on a two-dimensional, elasto-plastic, plane-strain model can be found in [77]. The agreement between simulation and experimental data is quite good for the intermediate strain range, for example, ∆ε = 0.60%. This is rather remarkable as all the parameters were estimated only on the basis of literature data. Some discrepancy, however, exists for low strain ranges (e.g. 0.50% and 0.55%), where the simulation tends to overestimate the crack density. For high strain ranges, however, it underestimates crack densities. Efforts in the second step are devoted to find the proper parameters for better simulations. 5.2.2 Parameters Study of Elasto-Plastic Models 5.2.2.1 Critical Shear Stress In Eq. (4-8) critical shear stress τc is a decisive parameter. In order to find the optimum value of τc which may lead to the best simulation results, a method based on variance estimation is applied. 83 The idea is: A group of trial values, τc1, τc2,…τcM, are chosen and the two-dimensional simulations are repeated with these trial values. For the ith strain range ∆ει and the jth value τcj, the crack densities from simulation i jkdC ) ~( versus the numbers of cycles are calculated, i = 1, 2, … n, j = 1, 2, … M, where n is the number of given strain ranges and M is the number of trial values of τ i jkN cj. The subscript k corresponds to the crack density and the number of cycles at the kth observation time, k = 1, 2, … ni, where ni is the number of observations. In this way the crack density and the number of cycles N for all trial values of τ dC ~ cj under the applied strain ranges ∆ει can be obtained. The relation of from simulation results can be expressed by a power function, NCd −~ bi j i jd NAC )() ~( = (5-1) where A and b are fitting coefficients. From the curves, the variance between the simulation results and the experimental data under the i NCd −~ ijQ th strain range can be calculated by ∑ = −= in k i jkd i kdi i j CCn Q 1 2])~()[(1 (5-2) where is the crack density of experimental data with respect to the kikdC )( th observation time. One illustration is given in Fig. 5- 4. The average variance jQ over all the strain ranges represents the degree of the approximation of the simulation results to the experimental data of a trial τcj, ∑ = = n i i jj Qn Q 1 1 (5-3) By comparing the variances of all trial value of τcj, the minimum average value of variance minQ can be found and the corresponding value of τcj is considered to be the optimum value of τc. 84 ),...,( 21min MQQQMinimumQ = (5-4) Fig. 5- 4 Scheme for illustrating the estimation of τc The selected values of τcj are from 130 MPa to 170 MPa for the plane-strain model and 80 MPa to 120 MPa for the plane-stress model. Obtained results are shown in Table 5-2 for plane- strain and Table 5-3 for plane-stress models respectively. The sign minus in Table 5-2 and Table 5-3 indicates: 0)()~( <− idijd CC It is found that the average variance Q reaches the minimum when τc = 160 MPa for the plane-strain model and τc = 100 MPa for plane-stress. Both of the values are comparable to the data in literature (108 MPa for the 0.37% carbon steel and 146 MPa for SAE 1045 normalized steel in [4, 41]). Because of the existence of a large amount of dislocations in the martensite phase and the dislocations are possible obstacles to the movement of mobile dislocations, a higher τc is expected. Considering that the results shown in this chapter are from the two-dimensional 0 10 20 30 10 100 1000 N Cr ac k de ns ity [m m -2 ] Test data 2D Elasto-plastic, 50 2D Elasto-plastic, 128 Fit - 5 0 Fit - 128 0)()~( 1 >−+ idijd CC 0)()~( 1 <−+ idijd CC τcj+1 τcj τcj+1 τcj 85 analyses and the stress state and grain misorientation effect in the three-dimensional model may change the tendency of the Cd-N curves, the further discussion of τc, therefore, will be resumed in the next chapter. Table 5-2 Variances between results of simulation model (PE) and experiment with various τcj (MPa) τc = 130 τc = 140 τc = 150 τc = 160 τc = 170 ∆ε = 0.50% 144.8122 62.0416 15.3668 2.6834 -3.9772 ∆ε = 0.55% 81.2944 28.8193 9.0484 0.4145 -4.0399 ∆ε = 0.60% 45.3783 19.5605 6.7906 -7.7841 -16.8199 ∆ε = 0.64% 61.8912 25.7441 9.4831 5.5603 -12.3919 ∆ε = 0.76% 11.1013 -20.2289 -34.0101 -46.5799 -60.9787 Average Q 68.8955 23.1873 1.3358 -9.1412 -19.6415 Table 5-3 Variances between results of simulation model (PS) and experiment with various τcj (MPa) τc = 80 τc = 90 τc = 100 τc = 110 τc = 120 ∆ε = 0.50% 315.2422 149.3560 3.9772 7.2564 0.5992 ∆ε = 0.60% 156.4593 94.3441 25.9893 10.5894 -11.9391 ∆ε = 0.76% 10.7224 -20.3063 -33.0215 -47.0384 -62.4658 average Q 160.8080 74.4646 -1.0183 -9.7309 -24.6019 5.2.2.2 Shear Modulus In a second step, the shear modulus G’ for plastic deformation regime is presented by a fitting curve of quadrilinear form obtained from FE stress analysis with the 200 grains 3D elasto- plastic model. The data obtained are listed in Table 5-4. The details of the determination of G’ can be found in Appendix D. 86 Table 5-4 Shear moduli in plastic deformation Shear strain γxy [-] Shear stress xyτ [MPa] Shear moduli G’ [GPa] 0.003 172 58.3 0.004 211 51.3 0.012 268 22.3 0.1 300 2.9 5.2.3 Relation of Crack Density to Cycles Fig. 5- 5 (a) – (c) show the crack density curves obtained from the two-dimensional elasto- plastic model under plane-stress with different critical shear stresses (the integers in the legend) and from experiments (Test data in the legend). These curves show that the number of one- segment cracks increases with the number of cycles and the strain amplitude. The values of the number of cycles to crack initiation vary with the values of τc. The larger the value of τc, the longer the initiation life is, as shown in Fig. 5- 5. By comparing the crack density curves obtained with the same value of τc for different strain ranges, in Fig. 5- 5 (a), (b) and (c), it is found that a single value of τc does not fit all strain ranges. A large τc fits to the low strain ranges but not to the high strain ranges, and vice versa. For example with τc = 110 MPa, the predicted crack initiation life N is close to the experimental data for the strain range of 0.50% as shown in Fig. 5- 5 (a), but not to the strain ranges of 0.60% and 0.76%, as shown in Fig. 5- 5 (b) and (c). The results of the parameter study in Tables 5-2 and 5-3 show the same tendency. 87 (a) ∆ε=0.50% 0 10 20 30 10 100 1000 10000 N C ra ck d en si ty [m m -2 ] Test data 2D Elasto-plastic, PS, 70 2D Elasto-plastic, PS, 90 2D Elasto-plastic, PS 110 Fit, 70 Fit, 90 Fit, 110 (b) ∆ε=0.60% 0 10 20 30 10 100 1000 10000 N C ra ck d en si ty [m m -2 ] Test data 2D Elasto-plastic, PS, 70 2D Elasto-plastic, PS, 90 2D Elasto-plastic, PS,110 Fit, 70 Fit, 90 Fit, 110 (c) ∆ε=0.76% 0 10 20 30 10 100 1000 10000 N Cr ac k de ns ity [m m -2 ] Test data 2D Elasto-plastic, PS, 70 2D Elasto-plastic, PS, 90 2D Elasto-plastic, PS, 110 Fit, 70 Fit, 90 Fit, 110 Fig. 5- 5 Crack density curves from elasto-plastic simulation with different τc 88 5.2.4 Effect of Microstructures The simulations with anisotropic elastic and elasto-plastic material properties are carried out using three RVE models, 2D100_1, 2D100_2 and 2D100_8. The three models have different grain morphologies and orientations, as shown in Fig. 5- 6 (a) – (c). (a) 2D100_1 (b) 2D100_2 (c) 2D100_8 Fig. 5- 6 Three different RVE models The derived relations between crack density Cd and the number of cycles to the initiation N from the three RVE structures show a similar tendency. The crack density curves, Cd-N, from the elastic model and the elasto-plastic plane-strain model do not show a strong influence of the microstructure on the crack initiation life. The results from the elastic plane-strain model for ∆ε = 0.50% are shown in Fig. 5- 7 (a). The scatter from elasto-plastic plane-stress model, however, is rather significant, as presented in Fig. 5- 7 (b). The most remarkable scatter, which is comparable to the scatter of the test data, happens in the simulation with the model 2D100_8 (in Fig. 5- 6 (c)). But in general, the scatter between the simulation results is much lower than that of experimental data. Besides microstructure parameters, it seems that there are additional reasons responsible for the scatter of experimental data. Inhomogeneous material property is the most likely one, which 89 is not taken into account in the present model. Inhomogeneous material property may result from tiny microscopic defects e.g. composition segregations, or submicrostructure such as second phrase precipitates, or from internal stress caused by the martensitic phase transformation. Due to the microscopic inhomogeneity the local τc of the material can also vary from point to point. Another possible reason is that scatter arises from the procedure of fatigue data processing. Elastic model, ∆ε = 0.50% 0 10 20 30 100 1000 10000 N C ra ck d en si ty [m m -2 ] 2D100_1, PE, 108 2D100_2, PE, 108 2D100_8, PE, 108 (a) Simulation results from elastic plane-strain model Elasto-plastic model, ∆ε=0.50% 0 10 20 30 100 1000 10000 N C ra ck d en si ty [m m -2 ] 2D100_1, PS, 100 2D100_2, PS, 100 2D100_8, PS, 100 Test data (b) Simulation results from elasto-plastic plane-stress model Fig. 5- 7 Crack density obtained from simulation of three RVEs, elastic-plastic models, PS 90 5.3 Effect of Stress Redistribution on Crack Initiation Sequence In the simulation the crack are initiated along the crack paths according to defined failure criteria. The initiated crack is introduced into the RVE by node separation, whose length is identical to the crack path, and simulates a one-segment crack observed in the experiments. As shown in Eq. (4-8), the number of cycles to crack initiation ∆Ni is a function of local resolved shear stress range ∆τi and the length of the crack path di ∆Ni = f(∆τi , di ) (5-5) This implies that a good candidate for an initiated crack within the model is a large crack path oriented at about ±45 º to the loading axis and with a high stress level. This is consistent with the crack patterns from simulation. Most of the initiated cracks are from large grains and are oriented about ±45° to the loading axis. After a crack has been introduced into the model, the stress is redistributed. The crack initiation sequence might be influenced by this stress redistribution. Three typical kinds of phenomena which violate the above prediction are found from both of elastic and elasto-plastic models. They are described in more detail using the following three examples. The crack initiation sequences described in this section are obtained from the elastic model. The crack pattern of the elasto-plastic model will be shown in the next section. 5.3.1 First Example This example shows how the stress redistribution induced by the initiated crack changes the crack initiation sequence. It is taken from the simulation with the elastic model 2D100_1 subjected to the stain range ∆ε = 0.76%. As shown in Fig. 5- 8 (a), the first crack, crack 1, is initiated in grain 23, in which the local shear stress is at the highest level. After crack 1 has been initiated (see Fig. 5- 8 (b)) the shear stress distribution, compared to the uncracked RVE, is changed. According to the simulation results from the uncracked RVE, the crack initiation 91 process should follow the sequence in column “Uncracked RVE” of Table 5-5, where the numbers are crack path indexes. These crack paths are all long ones and oriented in the preferential directions. The first crack, crack 1, is path 208 and the next crack would be path 20, along which the shear stress is the highest after path 208, if the stress would not redistribute. But the next initiated crack, crack 2, is path 131 instead of path 20. Path 131 is located in the adjacent grain of crack 1 and oriented along the band of high stress caused by the crack tip of crack 1. Although the stress redistribution area is limited to only a few adjacent grains, the shear stress field along the crack path 131 is strongly influenced. It leads to the initiation of crack 2 and changes the crack initiation sequence. If crack paths in adjacent grains are not within the high stress zone of the initiated crack, however, the sequence in Table 5-5 will be followed, as the crack 3 shown in Fig. 5- 8 (b). The two cracks Fig. 5- 8 (b), crack 1 and crack 2, are likely to coalesce and to form a kinked crack. Table 5-5 Sequences of crack initiation from simulation, RVE 2D100_1 Crack path index Crack No. Grain No. Uncracked RVE With crack 1 1 23 208 - 2 86 20 131 3 3 123 20 92 (a) First crack (b) Second and third cracks Fig. 5- 8 Crack initiation process and shear stress redistribution, 2D100_1, ∆ε = 0.76%, elastic model, PE 1 2 1 3 5.3.2 Second Example The simulated cracks are mostly initiated from long crack paths with orientations around ±45º, but that is not the case in all the simulations. It is found that some small cracks do initiate. This is true for the simulations of the two-dimensional elastic and elasto-plastic models when the applied strain amplitude is high. One example is taken from the simulation with model 2D100_1 subjected to strain range ∆ε = 0.76%, with 10 cracks initiated. As can be seen in Fig. 5- 9, crack 10 is an exception. Obviously the local stress field of the crack path has been enhanced by the high stress field around the tips of the two previously initiated cracks, crack 4 and crack 7 in this example. Crack 10 is initiated in the slip band, which is very short and oriented in a very low angle to the loading direction. It tends to link the two adjacent cracks to form a three-segment crack. 93 10 9 8 7 6 5 4 3 2 1 Fig. 5- 9 Crack initiation process in 2D elastic model, PE, ∆ε = 0.76% 5.3.3 Third Example This example will illustrate an unexpected change of the crack initiation sequence found in the simulation. In the elastic model the crack initiation sequence is considered to be the same for different strain ranges. It is true for most of the results but not for all results. Sometimes the initiation sequences are different. One example is found in 2D100_1, plane-strain. The crack path indices of two strain ranges are given in Table 5-6. The crack initiation sequence of 4, 5 and 6 is path 321, 77 and 123 for ∆ε = 0.76%. For ∆ε = 0.50%, however, the sequence is changed to path 123, 321 and 77. This phenomenon is depicted with the help of a figure, Fig. 5- 10. The two N curves, N1 and N2, are associated to two crack paths, path 1 and path 2, the lengths of which are d1 and d2 respectively, and d1 < d2. N1 and N2 curves monotonously decrease with the increase of stress amplitude, but the decreasing rates of the two N curves are different. This depends on two variables (as indicated by Eq. (5-5)), the average shear stress in the crack path and the length of the crack path. N’1 and N’2 can be determined from the two N curves. For a certain strain range, N1 of path 1 is lower than N2 of path 2. For another certain strain range, which is α times higher 94 than the former strain range, N’1 is higher than N’2. The changes of crack initiation sequence occur only under some particular conditions. The details along with a mathematical proof can be found in Appendix C. Table 5-6 Crack initiation sequence Crack initiation sequence Path indices for ∆ε = 0.50% Path indices for ∆ε = 0.76% 1 208 208 2 131 131 3 20 20 4 321 123 5 77 321 6 123 77 Fig. 5- 10 Illustration of the number of cycles N varying with shear stresses in the crack path 0 5000 10000 15000 20000 100 120 140 160 180 200 Shear stress [MPa] N N1(d1) N2(d2) d2 > d1 α α N1’>N2’ N10, α>0, β>0, γ>0, A>0, d>0, i.e. all the constants are positive. Additionally, the following requirements should also be satisfied: τ >τc, γτ>τc, i.e. 1<= λτ τ c Table C-1 Variables and factors of two paths under two strain ranges Path 1 Path 2 Stress when ε1 applied 1τ =τ 2τ =γ 1τ =γτ Stress when ε2 applied 1τ =ατ 2τ =αγ 1τ =αγτ Path size d1=d d2=βd1= βd Suppose the numbers of cycles to crack initiation for the two paths are ∆Ν1 and ∆Ν2 when the model is subjected to ε1 and ∆Ν’1 and ∆Ν’2 when it is subjected to ε2. Αccording to Eq. (C-1), the above conditions can be written as: 22 11 1 )()( cc d A d AN ττττ −=−=∆ , 22222 )()( cc d A d AN τγτβττ −=−=∆ 22 11 1 )()( ' cc d A d AN τατττ −=−=∆ , 22222 )()(' cc d A d AN ταγτβττ −=−=∆ If ∆ N1<∆ N2 when ε1 is applied and ∆ N’1>∆ N’2 when ε2 is applied, the crack sequences are not the same. The following derivation is aimed to find the valid conditions under which the above situation occurs. 1. Applied strain ε1 If ∆ N1< ∆ N2 when the applied strain is ε1, the following inequality should be valid: 22 )()( cc d A d A τγτβττ −<− Since A>0, d>0, the above inequality is written as: 22 )( 1 )( 1 cc τγτβττ −<− 137 ∴ 22 )()( cc τγτβττ −>− Since τ > τc >0, β>0, γτ > τc therefore the positive square root is the solution: )( cc τγτβττ −>− That is: cτβτγβ )1()1( −<− (C-2) Discussion: (1) If 01 >−β , we get 1>β and if 01 >−γβ , we get 1>γβ , i.e. 1),max( << γλβ . The following inequality, derived from Eq. (C-2), is valid: cτγβ βτ 1 1 − −< (2) If 01<−β and 01 >−γβ , the inequality Eq. (C-2) cannot be satisfied. (3) If 01<−β and 01 <−γβ , we get 1<β and 1<γβ , i.e. βγ 11 << . The following inequality is valid: cτγβ βτ 1 1 − −> (4) If 01 >−β and 01<−γβ , 1>β , 1<γβ , we get ),max( λβ βγ 1<< . The following inequality is valid: cτγβ βτ 1 1 − −> 2. Applied strain ε2 If ∆ N’1> ∆ N’2 when the model is subjected to ε2 and the following inequality should be valid: 22 )()( cc d A d A ταγτβτατ −>− Since A>0, d>0, the above inequality becomes: 22 )( 1 )( 1 cc ταγτβτατ −>− ∴ 22 )()( cc ταγτβτατ −<− If α>1, then ατ >τc, αγτ >τc, the positive square root is the solution: 138 )( cc ταγτβτατ −<− That is: ατγβτβ )1()1( −<− c (C-3) Discussion: (1) If 01 >−β and 01 >−γβ , we get: 1>β , 1>γβ , i.e. 1),max( << γλβ , and The following inequality is valid: α τ γβ βτ c⋅− −> 1 1 (2) If 01<−β and 01 >−γβ 1<β , 1>γβ , we get βγ 1> . The following inequality is valid: α τ γβ βτ c⋅− −> 1 1 (3) If 01 <−β and 01<−γβ , we get 1<β , 1<γβ , i.e. βγ 11 << . The following inequality is valid: α τ γβ βτ c⋅− −< 1 1 (4) If 01 >−β and 01<−γβ , the inequality Eq. (C-3) cannot be satisfied. 3. Conclusions Both of the inequalities for ε1, Eq. (C-2) and ε2, Eq. (C-3) should be satisfied simultaneously. From discussions 1.(1) and 2.(1), we find, when α>1, 1>β and 1),max( << γλβ , a τ that satisfies the following inequality can be found, c c τγβ βτα τ γβ β 1 1 1 1 − −<<⋅− − (C-4) Therefore inequality Eq. (C-4) is valid. 5. Example 139 An example is given as follows for further illustration. Consider two crack paths, d=d1=8.812, d2=14.769, i.e. β=1.676, when α=1.529, τc =108 and τ =τ1=254, τ2= 220, i.e. γ =0.866, λ=0.491, 7724.01 =β . These parameters satisfy the conditions: α>1, 1>β , 1),max( << γλβ and )262( 1 1)171( 1 1 c c τγβ βτα τ γβ β − −<<⋅− − The ratio of the number of cycles of the two paths under the two strain ranges can be calculated by 986.0 )( )(/ 2 2 21 =− −=∆∆ c cNN ττ τγτβ 112.1 )( )('/' 2 2 21 =− −=∆∆ c cNN τατ ταγτβ respectively. Therefore, for strain range ε1, the number of cycles of path 1 ∆N1 is smaller than that of path 2 ∆N2 , i.e. ∆N1<∆N2, then path 1 is the first initiated crack. But for strain range ε2, the number of cycles of path 2 ∆N’2 is smaller than that of path 1 ∆N’1, i.e. ∆N’2 <∆N’1, thus path 2 is the first initiated crack. 140 Appendix D Determination of Shear Moduli In a RVE the average τ will approach to the τm in an isotropic material under the given strain when the RVE is sufficiently large. Therefore, a 200 grains three-dimensional model is used. The RVE model is subjected to the shear strain γxy in XY-plane. The shear modulus G’ is determined by the average shear stress xyτ , which is calculated by FE analysis, divided by the given shear strain γxy. xyxyG γτ /= (D-1) The average shear stress xyτ is derived from Eq. (D-2), i m i i V xyxy VV dVzyx V xy∑∫ === 1 1),,(1 τττ (D-2) where ),,( zyxxyτ is the stress distribution function on the three-dimensional model, ixyτ is the shear stress component at the center of the ith element, i=1, 2, m. The shear moduli G’, derived by Eq. (D-1) and (D-2) at several strains from the obtained shear stresses, are shown in Table D-1. The value in the first line of Table D-1 is in the elastic regime, which is very close to the value of macroscopic elastic shear modulus, 81 GPa. The shear stress-strain curve is shown in Fig. D-1. 141 Table D-1 Shear moduli in plastic deformation Shear strain γxy [-] Shear stress xyτ [MPa] Shear moduli G’ [GPa] 0.001 78.9 78.9* 0.003 172 58.3 0.004 211 51.3 0.012 268 22.3 0.1 300 2.9 * This value is the elastic shear modulus 0 100 200 300 400 0.00 0.04 0.08 0.12 Shear strain [-] Sh ea r st re ss [M P a] Fig. D-1 Shear stress-strain curve derived from the simulation model, 3D200, elasto-plastic model 142 REFERENCES [1] Suresh, S., Fatigue of materials, 2nd version, Cambridge, 1998 [2] Brückner –Foit, A., Huang, X., and Motoyashiki, Y., Mesoscopic simulation of damage accumulation under fatigue loading, ECF15 (The 15th European Conference of Fracture), Stockholm, 2003 [3] Weiss, J. and Pineau, A, Continuous and sequential multiaxial low-cycle fatigue damage in 316L stainless steel, ASTM STP1191, Advances in Multiaxial Fatigue, Eds: McDowell, L. and Ellis, R., pp. 183–203, 1993 [4] Hoshide, T. and Socie, D.F., Crack nucleation and growth modelling in biaxial fatigue, Engng. Fract. Mech., Vol. 29, pp. 287–99, 1988 [5] Miller, K.J., The behavior of short fatigue cracks and their initiation, Fatigue & Fracture of Engng Mater. & Struct., Vol. 10, No. 2, pp. 93–113, 1987 [6] Blochwitz, C. and Richter, J., Plastic strain amplitude dependent surface path of microstructurally short fatigue crack in face-centred cubic metals, Mater. Science and Engng, Vol. A367, pp. 120–29, 1999 [7] Hu, Y.M., Floer, W., Krupp, U., Christ, H. –J., Microstructurally short fatigue crack initiation and growth in Ti –6.8Mo –4.5Fe –1.5Al, Mater. Sci. Engng, Vol. A278, pp. 170–80, 1999 [8] Long, M., Crooks, R. and. J Rack, H, High –cycle fatigue performance of solution –treated metastable β Titanium alloys, Acta mater. Vol. 47, No. 2, pp. 661–69, 1999 [9] Polak, J., Cyclic plasticity and low cycle fatigue life of metals, Elsevier, pp. 28–32, 1991 [10] Dougherty, J.D., Srivatsan, T.S. and Padovan, J., Cyclic stress response, strain resistance and fracture behavior of modified 1070 steel, Engng Fract. Mech., Vol. 53, No. 6, pp. 829 –847, 1996 [11] Farsetti, P. and Blarasin, A., Int. J. Fatigue, Vol. 10, pp. 153–161, 1988 [12] Carstensen, J.V., and Magnin, T., Characterization and qualification of multiple crack growth during low cycle fatigue, Int. J. Fatigue, Vol. 23, pp. 195–200, 2003 143 [13] Wu, Z. and Sun, X., Multiple fatigue crack initiation coalescence and growth in blunt notched specimens, Engng Fract. Mech., Vol. 59, No.3, pp. 353–59, 1998 [14] Hoshide T., Yamada, T., Fujimura, S. and Hayashi, T., Short crack growth and life prediction in low cycle fatigue on smooth specimens, Eng Fract. Mech., Vol. 21, No.1, pp. 85–101, 1985 [15] Boyd –Lee, A.D., Fatigue crack growth resistant microstructures in polycrystalline Ni –base superalloys for aeroengines, Int. J. of Fatigue, Vol. 21, pp. 393–405, 1999 [16] Sommer, C., Mughrabi, H., and Lochner, D., Influence of temperature and carbon content on the cyclic deformation and fatigue behavior of α –iron, Parts I + II, Acta Materiallia, Vol. 45, No. 5, pp. 1527–46, 1998 [17] Vasek, A. and Polak, J., Low cycle fatigue damage accumulation in ARMCO –iron, Fatigue Fract. Engng Mater. Struct., Vol. 14, No. 2/3, pp. 193–204, 1991 [18] Han, S. –W. Kumai, S. and Sato, A., Effects of solidification structure on short fatigue crack growth in Al–7%Si–0.4%Mg alloy castings, Materials Science and Engineering, Vol.A332, pp. 56–63, 2002 [19] Johansson, T., Kullig, E., Brückner-Foit, A. and Riesch-Oppermann, H., A fracture mechanics model for interacting cracks in thermal fatigue, Proc. ECF 11, Vol. I, EMAS, pp. 257–62, 1986 [20] Frise PR, Bell R., Modelling fatigue crack growth and coalescence in notches, Int. J. Pressure Vessel & Piping, Vol. 51, pp. 107–26, 1992. [21] Tryon R.G. and Cruse T.A., A reliability –based model to predict scatter in fatigue crack nucleation life, Fatigue & Fracture of Engng Mater. & Struct.,Vol. 21, pp. 257–67, 1998 [22] Brückner-Foit, A., Meyer, S., Möslang, A. and Diegele, E., Stochastic simulation of fatigue damage evolution in a martensitic steel, FATIGUE 2002, Stockholm, Sweden, May, 2002 [23] Huang, X., Wu, X.R., Ding, Ch. and Hu, B.R., Small Crack propagation behavior of a forged titanium alloy Ti –6Al –4V, FATIGUE 2002, Stockholm, Sweden, May, 2002 [24] Ewing, J.A. and Humfrey, J.C.W., The fracture of metals under repeated alternations of stress, Philosophical transactions of the Royal Society, London, Vol. A200, pp. 241–50, 1903 [25] Mughrabi, H., The cyclic hardening and saturation behavior of copper single crystals, Mater. Sci. Engng, Vol. 33, pp. 207–23, 1978 144 [26] Winter, A.T., A model for the fatigue of cooper at low plastic strain amplitude, Philosophical Magazine, Vol. 29, pp. 719 –34, 1974 [27] Gong, B., Wang, Z. and Wang, Z., Cyclic deformation behavior and dislocation structures of [001] copper single crystals –I cyclic stress –strain response and surface feature, Acta Mater., Vol. 45, No. 4, pp. 1365–77, 1997 [28] Pohl, K., Mayr, P., and Macherauch, E., Persistent slip bands in the interior of a fatigued low carbon steel, Scripta Metallurgica, Vol. 14, pp. 1167–69, 1980 [29] Sasaki, S. and Ochi, Y., Some experimental studies of fatigue slip bands and persistent slip bands during fatigue process of low carbon steel, Eng Fract. Mech. Vol. 12, pp. 531–40, 1979 [30] Man, J., Obrtlı´k, K., Blochwitz, C., Pola´k, J., Atomic force microscopy of surface relief in individual grains of fatigued 316L austenitic stainless steel, Acta Materialia Vol. 50, pp. 3767– 80, 2002 [31] Nakai, Y., Fukuhara, S. and Ohnishi, K., Observation of fatigue damage in structural steel by scanning atomic force microscopy, Int. J. Fatigue, Vol. 19, No. 1, pp. S223–36, 1997 [32] Ma, B. –T. and Laird, C., Overview of fatigue behaviour in copper single crystals –II. Population size, distribution and growth kinetics of stage I cracks for tests at constant strain amplitude, Acta Mater., Vol. 37, pp. 337–48, 1989 [33] Katagiri, K., Omura, A., Koyanagi, K., Awatani, J., Shiraishi, T., and Kaneshiro, H., Early stage crack tip morphology in fatigued copper, Metallurgical Transactions, Vol. 8A, pp. 1769 – 73, 1977 [34] Cheng, A.S., and Laird, C., Fatigue life behavior of copper single crystals, Part I: Obser – vations of crack nucleation, Fatigue of Engng. Mater. & Struct., Vol. 4, pp. 331–41, 1981 [35] Essmann, U., Gösele, U. and Mughrabi, H., A model of extrusions and intrusions in fatigued metals, I. Point –defect production and the growth of extrusions, Philosophical Magazine, Vol. A 44, pp. 405–26, 1981 [36] Lin, T.H. and Ito, Y.M., Mechanics of fatigue crack nucleation mechanism, J. Mech. Phy. Solid, Vol. 17, pp. 511–23, 1969 [37] Tanaka, K. and Mura, T., A dislocation model for fatigue crack initiation, J. of Applied Mechanics, Vol. 48, pp. 97–103, 1981 145 [38] Kunio, T. Shimizu, M., Yamada, K., Sakara, K. and Yamamoto, T., The early stage of fatigue crack growth in martensitic steel, Int. J. Fract., Vol. 17, No. 2, pp. 111–19, 1981 [39] Laz, P.J., Hillberry, B.M, Fatigue life prediction from inclusion initiated cracks, Int. J. Fatigue, Vol.30, No. 4, pp. 263–70, 1997 [40] Cheng, A.S. and Laird, C., fatigue life behavior of copper single crystals, part I and Part II. Fat. Fract. Eng. Mater. Struct., Vol. 4, pp. 331–41, pp. 343–53, 1981 [41] Hoshide T. and Kusuura K., Life Prediction by Simulation of Crack Growth in Notched Components with Different Microstructures under Multiaxial Fatigue, Fatigue & Fracture of Engng Mater. & Struct., Vol. 21, pp. 201–13, 1998 [42] Zimmermann, A. and Rie, K. –T., Simulation of microcrack formation and growth during cyclic loading considering microstructure, Fatigue 2002 –Proceedins of the 8th international fatigue conference, Vol. 4/5, 2002 [43] Alexandre F., Deyber S. and Pineau A., Modelling the optimum grain size on the low cycle fatigue life of a Ni based superalloy in the presence of two possible crack initiation sites, Scripta Materialia, Vol. 50, pp. 25–30, 2004 [44] Chan, K., A microstructure –based fatigue –crack –initiation model, Metallurgical and Material transactions A, Vol. 34A, pp. 43–58, 2003 [45] Harvey, S.E., Marsh, P.G. and Gerberich, W.W., Acta Metall. Mater., Vol. 42, pp. 3493–502, 1994 [46] Goto, M., Statistical investigations of the behavior of microcracks I carbon steel, Fatigue Fract. Engng Mater. Struct., Vol. 14, No.8, pp. 833–45, 1991 [47] Ochi, Y., Ishii, A, Sasaki, K., An experimental and statistical investigation of surface atigue crack initiation and growth, Fatigue Fract. Eng Mater. Struct., Vol. 8, pp. 32–39, 1985 [48] Morris, W.L., James, M.R. and Buck, O., Computer simulation of fatigue crack initiation, Eng Fract. Mech., Vol. 13, pp. 213–21, 1980 [49] Chang, R., Morris, W.L. and Buck, O., Fatigue crack initiation at intermetallic particles in aluminium alloy – A dislocation pile –up model, Scripta Metallurgica, Vol. 15, pp. 1941–94, 1979 146 [50] Ihara, C. and Tanaka, T., A stochastic damage accumulation model for crack initiation in high- cycle fatigue, Fatigue Fract. Engng Mater. Struct., Vol. 23, pp. 375–80, 2000 [51] Meyer, S., Brückner –Foit, A., Möslang, A., A stochastic simulation model for micro –crack initiation in a martensitic Steel, Computational Materials Science, Vol. 26, pp. 102–10, 2003 [52] Nemat –Nasser, S., Modelling Small Deformations of Polycrystals, edited by John Gittus and Joseph Zarka, Elsevier Applied Science Publishers, 1986 [53] Lemaitre, J., Formulation and identi:cation of damage kinetic constitutive equations. In: Krajcinovic, D., Lemaitre, J. (Eds.), Continuum Damage Mechanics Theory and Applications. Springer, Wien–New York, 1987 [54] Ren, Z. –Y., Zheng, Q. –S., A Quantitative study of minimum sizes of representative volume elements of cubic polycrystals—numerical experiments, J. Mech. & Phy. Solids, Vol. 50, pp. 881 – 93, 2002 [55] Ahmadi, A.and Zenner, H., Simulation of microcrack propagation for different load sequences, FATIGUE 2002, Stockholm, Sweden, May, 2002 [56] Bomas, H., Hünecke, J., Laue, S. and Schöne, D., Anrisslebensdauervorhersage am Beispiel glatter und gekerbter Proben des Werkstoffes Cm15, Mechanistenorientierte Lebensdauer – vorhersage für zyklisch beanspruchte metallische Werkstoffe, Abschlusskolloquium des Schwerpunktprogramms der DFG, pp. 53–68, 2004 [57] Anteboth, S., Brückner –Foit, A., Hoffmann, M.J. and Sutter, U., Simulation des elektro – mechanischen Verhaltens von PZT mit realer Domänenstruktur, Ceramic Forum Int., DKG, Vol. 13, pp. 179–81, 2005 [58] Kozaczek, K.J, Petrovic, S.G. and Ruud, C.O., Mcllree, A.R., Modelling of stress distribution on the microstructural level in alloy 600, Fatigue and crack growth: Environment effects, modeling studies and design considerations, PVP, ASME1995, Vol. 306, pp. 223–32, 1995 [59] Weyer, S., Brückner –Foit, A. and Fröhlich, A., Overall properties of ceramics subjected to compressive loading, Int. Conf. on Engng Ceramics and Struct., Ceramic Engng and Science Proc. Vol. 21, pp. 101–107, 2000 147 [60] Weyer S., Fr€ohlich A., Riesch-Oppermann H., Cizelj L., Kovac M., Automatic finite element meshing of planar Voronoi tessellations, Engineering Fracture Mechanics Vol. 69, pp. 945–58, 2002 [61] Andersson, J., The influence of grain size variation on metal fatigue, Int. J. Fatigue, Vol. 27, pp. 847–52, 2005 [62] Kumar S., Kurtz S.K., Banavar J.R. and Sharma M.G. Properties of a three –dimensional Poisson –Voronoi tesselation: a Monte Carlo study, J. Statistics Phy., Vol. 67, pp. 523–51, 1992 [63] Bertsch, J., Mikroskopische Untersuchung der Bildung von Ermüdungsrissen an zwei ferritisch –martensitischen Stählen im unbestrahlten und vorbestrahlten Zustand, Technik und Umwelt Wissenschaftliche Berichte FZKA 5984, Institut für Materialforschung, Forschungs –zentrum Karlsruhe, 1997 [64] Meyer, S., Entwicklung eines stochastischen Schädigungsmodells für einen martensitischen Stahl unter Ermüdungsbelastung, PhD Dissertation, Universität Kassel, 2002 [65] Armas, A.F., Petersen, C., Schmitt, R., Avalos, M., Alvarez-Armas, I., Mechanical and microstructural behavior of isothermally and thermally fatigued ferritic/martensitic steels, Journal of Nuclear Materials, Vol. 307-311, pp. 509–13, 2002 [66] Maday, M. –F. Phenomenological aspects of fatigue cracking in as –received and hardened F82H modifed steel exposed to lithiated water with dissolved hydrogen at 240°C, Journal of Nuclear Materials, Vol. 283–287, pp. 68–93, 2000 [67] Li, J.C.M., Microstructure and properties of materials, Ed. J.C.M. Li, Univ. of Rochester, World Scientific, 2000 [68] Kim, H.J., Kim, Y.H. and Morris Jr., J.M., Thermal mechanisms on grain and packet refinement in a lath martensitic Steel, ISI J International, Vol. 38, No. 11, 1998. [69] Guo, Z., Sato, K., Lee, T. –K. and Morris, Jr., J.W., Ultrafine grain size trough thermal treatment of lath martensitic steels, Ultrsfine Grained Materials, TMS 2000, eds., Rajiv S. Mishra, S.L. Semiatin etc., pp. 51–62, 2000 [70] Kelly, P.M., Jostsons, A. and Blake, RG., Acta Metall. Mater., Vol. 38, pp. 1075, 1990 148 [71] Weyer, S., Mikromechanisches Modell zur Berechnung effektiver Materialeigenschaften von geschädigten Polykristallen, Dissertation, Univ. Karlsruhe, Germany, 2001 [72] MSC Software, PCL PATRAN 2001, Reference Manual [73] Watanabe, O., Zbin, H.M. and Takenouchi, E., Crystal plasticity: micro –shear banding in polycrystals using Voronoi tessellation, Int. J. Plasticity, Vol. 14, No.8, pp. 771 –789, 1998 [74] Chen, W.-F. and Saleeb, A.F., Constitutive equations for engineering materials, Vol. 1, Elasticity and Modeling, Elsevier, 1982 [75] Kraska, M., Textursimulation bei großen inelastischen Verformungen mit der Technik des repräsentativen Volumenelements(RVE), Pro Business, pp. 35–37, 1996 [76] Ameen, M., Computational Elasticity, Theory of elasticity and finite and boundary element methods, Alpha Science Int. Ltd., pp. 68–77, 2005 [77] Brückner–Foit, A., and Huang, X., Numerical simulation of micro –crack initiation of martensitic steel under fatigue loading, Int. J. Fatigue. Vol. 28, pp. 963–71, 2006 [78] Hunsche, A. and Neumann, P., Quantitative measurement of persistent slip band profiles and crack initiation, Acta Metallurgica, Vol. 34, pp. 207–17, 1986 [79] Venkataraman, G, Chung, Y-W., Nakasone, Y. and Mura, T., Free energy formulation of fatigue crack initiation along persistent slip bands: calculation of S-N curves and crack depths, Acta Metall. Mater., Vol. 38, No. 1, pp. 31-40, 1990 [80] Goldstein, H., Classical Mechanics, Addison –Wesley, 1974 [81] Cook, R.D., Malkus, D. S., Plesha, M. E., Concepts and applications of finite element analysis, 3rd. edition, John Wiley & Sons, 1989 149 Erklärung Hiermit versichere ich, dass ich die vorliegende Dissertation selbständig und ohne unerlaubte Hilfe angefertigt und andere als die in der Dissertation angegeben Hilfsmittel nicht benutzt habe. Alle Stellen, die wörtlich oder sinngemäß aus veröffentlichten oder unveröffentlichten Schriften entnommen sind, habe ich als solche kenntlich gemacht. Kein Teil dieser Arbeit ist in einem anderen Promotions- oder Habilitationsverfahren verwendet worden. (M.E. Xinyue Huang)