SOME COUPLED DEFORMATION, DIFFUSION AND FRACTURE PROBLEMS IN SMALL SCALE STRUCTURES BY TANMAY K. BHANDAKKAR B.E., NAGPUR UNIVERSITY, 2002 M.E., INDIAN INSTITUTE OF SCIENCE, 2004 A DISSERTATION SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN THE SCHOOL OF ENGINEERING AT BROWN UNIVERSITY PROVIDENCE, RHODE ISLAND MAY 2011 © Copyright 2011 by Tanmay K. Bhandakkar This dissertation by Tanmay K. Bhandakkar is accepted in its present form by the School of Engineering as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date________________ ________________________________ Huajian Gao, Advisor Recommended to the Graduate Council Date________________ ________________________________ Eric Chason, Reader Date________________ ________________________________ Pradeep R. Guduru, Reader Approved by the Graduate Council Date________________ ________________________________ Peter Weber, Dean of the Graduate School iii Vita Tanmay K. Bhandakkar was born in Khamgaon, Maharashtra, India on August 19, 1980. He received his B.E. degree from Nagpur University in 2002, and his M.E. degree from Indian Institute of Science, in 2004. He entered the solid mechanics group at Brown University for pursuing Ph.D. degree in January 2006. iv Acknowledgements This thesis is an outcome of my rewarding stay at Brown for my Ph.D. in Mechanics of Solids. During my graduate studies, I have been privileged to work under the guidance of Prof. Huajian Gao. He has been very patient and encouraging with me throughout the course of the work. His approach of formulating and solving problem through physical and mathematical insights is very impressive and I have tried to imbibe it, through out all these years. I will always remember his welcome and enthusiasm, whenever I have knocked his office door for guidance. I am sincerely grateful for the opportunity to work with him and cherish this association. I would like to express my gratitude to the members of the thesis committee, Prof. Chason and Prof. Guduru, for taking time out of their schedule to review and critique my work. The fruitful collaboration with Prof. Chason, Jae Wook Shin, Prof. Curtin and Dr. Audrey Chng is gratefully acknowledged. I would like to thank the faculty at Brown for the courses, which contributed towards my education. I would like to thank Dr. Yujie Wei, who mentored me when I joined Brown. I would also like to thank all the current and former members of the research group of Prof. Gao. During the course of my stay at Brown, I have made many good friends and I would like to thank them for making my stay enjoyable. Especially I would like to convey my thanks to Prateek and Nikhil for their wonderful company. This research was supported by the National Science Foundation through Grant CMMI-0758535, and I would like to thank them for their support. I wish to express my sincere thanks to Patricia Capece, Peggy Mercurio, Stephanie Gesualdi, Jim Scheuerman, Nancy Baker, for their help over the years. v Finally, I would like to thank my family - my parents and my sister who have been a source of constant support and encouragement. vi Contents List of Tables……………………………………………………………………………...xi List of Figures……………………………………………………………………………xii 1. Overview and Introduction……………………………………………………………..1 1.1 Volmer-Weber thin film growth..…………………………………………………2 1.2 Stresses and cracking in Li battery electrodes…………………………………….4 1.3 Dislocation shielding of a cohesive crack…………………………………………6 1.4 Organization of this thesis………………………………………………………...7 2. Basic mathematical techniques and Formulation……………………………………..10 2.1 Integral equation method in fracture mechanics…………………………………10 2.2 Cohesive zone models……………………………………………………………17 2.3 Coupling between stress and diffusion…………………………………………..20 2.3.1 Mass conservation, driving force and kinetic laws of diffusion…………...20 2.3.2 Chemical potential…………………………………………………………22 2.3.3 Grain boundary diffusion…………………………………………………..22 2.3.4 Bulk diffusion……………………………………………………………...23 3. Volmer-Weber thin film growth: Formulation………………………………………. .25 3.1 Compressive stress evolution during Volmer-Weber growth……………………25 3.2 Problem formulation……………………………………………………………..27 3.3 Finite and inhomogeneous GB diffusivity……………………………………….31 3.4 Infinite GB diffusivity……………………………………………………………35 3.4.1 Analytical solution…………………………………………………………35 3.4.2 Numerical solution…………………………………………………………39 vii 4. Volmer-Weber thin film growth: Results……………………………………………...41 4.1 Finite and inhomogeneous GB diffusivity……………………………………….41 4.1.1 Formation of GB diffusion wedges during thin film growth………………42 4.1.2 Effect of inhomogeneous GB diffusivity on compressive stress evolution during thin film growth…………………………………………………….48 4.2 Infinite GB diffusivity……………………………………………………………51 5. Crack nucleation in battery electrodes under diffusion induced stresses……………..56 5.1 Background………………………………………………………………………56 5.2 Diffusion induced stresses in electrodes…………………………………………58 5.2.1 Strip electrode……………………………………………………………...59 5.2.2 Cylindrical electrode……………………………………………………….64 5.3 Cohesive model of crack nucleation in the electrode……………………………68 5.3.1 Strip electrode……………………………………………………………...69 5.3.1.1 Localization spacing p …………………………………………..71 5.3.2 Cylindrical electrode ………………………………………………………75 5.3.2.1 Localization spacing p …………………………………………..77 5.4 Critical conditions for crack nucleation………………………………………….79 6. Dislocation shielding of a cohesive crack…………………………………………….86 6.1 Dugdale cohesive crack………………………………………………………….87 6.1.1 Problem Formulation………………………………………………………87 6.1.2 Limiting solutions for high and low cohesive strengths…………………...93 6.1.3 Numerical results for Dugdale cohesive crack…………………………….94 6.2 Discrete Dislocation Simulations………………………………………………...99 viii 7. Concluding remarks………………………………………………………………….103 A. Numerical algorithm for hypersingular integral equations...……………………….109 B. Green’s function kernel in a strip (Periodic case)…………………………………..113 C. Numerical scheme for integral equations in a strip…………………………………118 D. Green’s function kernel for a single dislocation in a strip………………………….121 E. Green’s function kernel for a circular dislocation loop in a cylinder………………124 F. Numerical scheme for integral equations in a cylinder……………………………..130 G. Asymptotic interpolation solution for dislocation interaction with a cohesive crack.135 Bibliography……………………………………………………………………………139 Parts of this thesis have appeared in the following peer-reviewed publications: ƒ Bhandakkar, T.K., Chason, E., Gao, H., 2009. Formation of crack-like diffusion wedges and compressive stress evolution during thin film growth with inhomogeneous grain boundary diffusivity, International Journal of Applied Mechanics 1, 1-19. ƒ Bhandakkar, T.K., Chason, E., Gao, H., 2010. Analytical model of transient compressive stress evolution during growth of high diffusivity thin films on substrates, Philosophical Magazine 90, 3037-3048. ƒ Bhandakkar, T.K., Chason, E., Gao, H., 2010. Cohesive modeling of crack nucleation under diffusion induced stresses in a thin strip: implications on the critical size for flaw tolerant battery electrodes, International Journal of Solids and Structures 47, 1424-1434. ƒ Bhandakkar, T.K., Chng, A.C., Curtin, W.A., Gao, H., 2010. Dislocation shielding of a cohesive crack, Journal of the Mechanics and Physics of Solids 58, ix 530-541. x List of Tables 4.1 GB diffusivities of Ag, Cu and Sn at room temperature…………………………….51 5.1 Material properties of Si and operating parameters……….…………………….......64 xi List of Figures 1.1: Schematic illustration of stress evolution during high mobility Volmer-Weber growth – interruption……………………………………………………………………...3 1.2: SEM micrographs of germanium films after (a) 1, (b) 10, and (c) 110 discharge/charge cycles and optical micrograph of the sample after 300 cycles (d) [29].. 4 2.1: Crack of length 2a in an infinite plane with a far field stress σ ∞ . The problem can be solved by superposition of a uniform stress and a crack with no far field loading but with crack face pressure equal to σ ∞ …………………………………………………………11 2.2: (a) Ahead of a crack tip, the surface separation δ is related to the cohesive stress σ (δ ) with the maximum cohesive stress σ c referred to as cohesive strength of the material. When the tip opening displacement δ t is equal to the maximum cohesive interaction distance δ c , crack growth is assumed to occur. Γ is the cohesive energy and the cohesive zone size is L. (b) Cohesive zone technique applied to the centre crack problem of Figure 2.1. Cohesive zone acts in the region [− a,−a + c ] and [a − c, a ] with a given cohesive zone law σ = σ (δ ) .............................................................................................17 2.3: A Closed contour lying on a surface………………………………………………...20 3.1: Schematic illustration of a growing polycrystalline thin film with grain size L and increasing film thickness h(t ) on a substrate. During the film growth, atoms diffuse into the GBs due to the tensile coalescence stress in the film associated with GB formation, xii leading to the formation of crack-like GB diffusion wedges. Furthermore, the higher chemical potential at the film surface can overdrive adatoms into GBs and result in a final compressive stress in the film [26-27]………………………………………………….. 28 3.2: Energy landscape of atoms near the surface-GB junction for the process of adatom insertion into a GB, where Δμ d is the chemical potential between surface and the GB, and ΔG ∗ is the activation barrier of diffusion between surface and GB………………...30 3.3: Two auxiliary universal functions in the analytical model of compressive stress evolution during growth of a thin film. The functions (a) λ (h L ) and (b) δ (h L ) along with their fitting function given in Equations (3.35) and (3.36) respectively [85]………37 4.1: Inhomogeneous GB diffusivity profiles represented by the Dirac function g (z ) = 2.0 [1 + exp(βz h0 )] . The higher the value of β , the more inhomogeneous the GB diffusivity becomes………………………………………………………………………41 4.2: Formation of a GB diffusion wedge during film growth ( H = 1 + H& τ ) under uniform GB diffusivity ( β = 0 ). Evolutions of (a) the GB traction and (b) the GB wedge profiles at different time instants during deposition for parameter choices of H& = 1 , A = 1 , h0 L = 0.5 . The stress has been normalized as σˆ = (σ − σ s ) (σ 0 − σ s ) and the GB opening displacement has been normalized as uˆ = uE ∗ [2πh0 (σ 0 − σ s )] , where σ 0 is the initial film stress, σ s is stress corresponding to the excess surface chemical potential during deposition and h0 is the initial film thickness. The time has been normalized as τ = t t GB where t GB is given by Equation (3.22) [86]……………………………………………...43 4.3: Formation of a GB diffusion wedge during film growth ( H = 1 + H& τ ) under strongly inhomogeneous GB diffusivity ( β = 4 ). Evolutions of (a) the GB traction and (b) the GB xiii wedge profiles at different time instants during film deposition for parameter choices of H& = 1 , A = 1 , h0 L = 0.5 . The stress has been normalized as σˆ = (σ − σ s ) (σ 0 − σ s ) and the GB opening displacement has been normalized as uˆ = uE ∗ [2πh0 (σ 0 − σ s )] , where σ 0 is the initial film stress, σ s is stress corresponding to the excess surface chemical potential during deposition and h0 is the initial film thickness. The time has been normalized as τ = t t GB where t GB is given by Equation (3.22) [86]……………………………………45 4.4: Evolution of the normalized stress intensity factor near the tip of a GB diffusion wedge with the increasing film thickness under different deposition rates and uniform GB diffusivity ( β = 0 ). Other parameters used in the calculations are A = 1 and h0 L = 0.5 . ( ) The stress intensity factor has been normalized as Kˆ = K σ 0 πh0 , and the deposition rate has been normalized as H& = h& (h0 t GB ) , where σ 0 is the initial film stress, h0 is the initial film thickness and t GB is given by Equation (3.22) [86]………………………… 46 4.5: Effect of inhomogeneous GB mobility on the evolution of the normalized stress intensity factor near the tip of a diffusion wedge with increasing film thickness at the deposition rate of H& = 1 . Other parameters used in the calculation are A = 1 and ( h0 L = 0.5 . The stress intensity factor has been normalized as Kˆ = K σ 0 πh0 , and the) deposition rate has been normalized as H& = h& (h0 t GB ) , where σ 0 is the initial film stress, h0 is the initial film thickness and t GB is given by Equation (3.22) [86]………………. 47 4.6: A film deposition and growth interruption process. The thin film is deposited at a normalized deposition rate of H& = h& (h0 t GB ) from an initial thickness of 1 up to time τ1 followed by growth interruption till time τ 2 . The time has been normalized as τ = t t GB where t GB is given by Equation (3.22). h0 is the initial film thickness………………….48 xiv 4.7: Evolutions of stress-thickness during one cycle of deposition followed by growth interruption with varying degrees of GB inhomogeneity represented by different values of β in Equation (4.1). The results are shown at (a) a slower deposition rate of H& = 1 and (b) a faster deposition rate of H& = 100 . The deposition rate and time have been normalized as H& = h& (h0 t GB ) and τ = t t GB respectively, where σ 0 is the initial film stress, h0 is the initial film thickness and t GB is given by Equation (3.22). Other numerical parameters used in the calculation are σ 0 = 70 MPa , σ s = −26 MPa , A = 1 , h0 L = 0.5 [86]……………………………………………………………………………………….50 4.8: Variation of the steady state stress with deposition rate during film growth under varying degrees of GB inhomogeneity represented by different values of β in Equation (4.1). The deposition rate has been normalized as H& = h& (h0 t GB ) , where h0 is the initial film thickness and t GB is given by Equation (3.22). Other numerical parameters used in the calculation are σ 0 = 70 MPa , σ s = −26 MPa , A = 1 , h0 L = 0.5 [86]…………………51 4.9: Comparison of the predicted Stress-Thickness curve from analytical and numerical solutions with corresponding experimental data measured during a three cycle depositon- interruption experiment on Sn film deposited on Sn/Si substrate at room temperature [28,85]…………………………………………………………………………………... 53 4.10: Comparison of steady state stress versus deposition rate from analytical model and experiments on growth of Sn film on Sn/Si substrate at room temperature [85]………. 54 5.1: Schematic illustration of crack nucleation in a strip electrode of width 2h during galvanostatic solute (a) intercalation and (b) extraction, modeled as diffusion along the thickness direction (y-axis). The crack nuclei are uniformly spaced with period p and xv modeled as cohesive zones obeying the triangular traction-separation law (Equation (5.30))…………………………………………………………………………………… 60 5.2: Snapshot profiles of solute concentration and diffusion induced stress in a strip electrode. (a) Concentration during insertion, (b) DIS during insertion, (c) concentration during extraction and (d) DIS during extraction. The concentration is normalized as cˆ = cFD (Ih ) during insertion and cˆ = (c − c1 )FD (Ih ) during extraction (Equations (5.11) and (5.15)), while DIS is normalized as σˆ = 3(1 − ν )FDσ (EΩIh ) (Equations (5.12) and (5.16)) [96]……………………………………………………………………………… 63 5.3: Schematic illustration of crack nucleation in a cylindrical electrode of diameter 2rc during galvanostatic (a) intercalation and (b) extraction, modeled as diffusion along the radial direction (r-axis). The axisymmetric crack nuclei are uniformly spaced with period p and modeled as localized cohesive zones obeying the triangular traction-separation law (Equation (5.30))…………………………………………………………………….65 5.4: Snapshot profiles of solute concentration and diffusion induced stress in a cylindrical electrode. (a) concentration during insertion, (b) DIS during insertion, (c) concentration during extraction and (d) DIS during extraction. The concentration is normalized as cˆ = cFD (Irc ) during insertion and cˆ = (c − c1 )FD (Irc ) during extraction (Equations (5.23) and (5.27)), while DIS is normalized as σˆ D = 3(1 − ν )FDσ D (EΩIrc ) (Equations (5.24) and (5.28))…………………………………………………………………………………… 68 5.5: Effect of localization spacing p on the distribution of axial stress σ x in the region between two adjacent cohesive zones of length a h = 0.7 during insertion in a strip electrode (Figure 5.1a). (a) At the spontaneous localization spacing p h = 2.36 , the maximum axial stress between two adjacent localization zones is equal to the cohesive xvi strength σ c . (b) If the localization spacing is taken to be p h = 5 , the maximum axial stress between two adjacent localization zones is seen to be greater than the cohesive strength σ c . (c) If the localization spacing is taken to be p h = 1 , the maximum axial stress between two adjacent localization zones is seen to be lower than the cohesive strength σ c [96]………………………………………………………………………... 72 5.6: Localization spacing p as a function of the normalized cohesive strength for the formation of (a) centre cohesive zones during solute insertion and (b) edge cohesive zones during solute extraction, in a strip electrode. In both (a) and (b), the cohesive strength σ c is normalized by the peak stresses σ peak I = EΩIh [18(1 − ν )FD ] during solute insertion and σ peak E = EΩIh [9(1 − ν )FD ] during solute extraction [96]…………………...74 5.7: Localization spacing p as a function of the normalized cohesive strength for the formation of (a) centre cohesive zones during solute insertion and (b) edge cohesive zones during solute extraction, in a cylindrical electrode. In both (a) and (b), the cohesive strength is normalized by the peak DIS σ peak = EΩIrc [12(1 − ν )FD ] ……………………..78 5.8: The critical conditions for crack nucleation expressed as relationships between the normalized half-width of electrode, the normalized cohesive strength and the normalized critical size of cohesive zone at nucleation. The blue lines plot the critical dimension of electrode while the dashed green lines plot the critical size of cohesive zone at crack nucleation as functions of the normalized cohesive strength. (a, c) plot the critical conditions for nucleation of center cohesive zones of length 2a under solute insertion while (b, d) plot those of symmetric edge cohesive zones of length a under solute extraction in a strip electrode of width 2h . In (a) and (b), the cohesive strength is xvii I normalized by the peak stresses σ peak = EΩIh [18(1 − ν )FD] during solute insertion and E σ peak = EΩIh [9(1 − ν )FD ] during solute extraction. In (c) and (d), the cohesive strength is normalized by the size-independent reference stress σ ref = EΩIl ft [18(1 − ν )FD ] where l ft is the characteristic length scale defined in Equation (5.38) [96]………………………. 80 5.9: The critical conditions for crack nucleation expressed as relationships between the normalized radius of electrode, the normalized cohesive strength and the normalized critical size of cohesive zone at nucleation. The blue lines plot the critical dimension of electrode while the dashed green lines plot the critical size of cohesive zone at crack nucleation as functions of the normalized cohesive strength. (a, c) plot the critical conditions for nucleation of center cohesive zones of radius a under solute insertion while (b, d) plot those of edge annular cohesive zones of size a under solute extraction, in a cylindrical electrode of radius rc . In (a), (b), the cohesive strength is normalized by the peak DIS σ peak = EΩIrc [12(1 − ν )FD ] . In (c) and (d), the cohesive strength is normalized by the size-independent reference stress σ ref = EΩIl ft [12(1 − ν )FD ] where l ft is the characteristic length scale defined in Equation (5.38)……………………….. 81 6.1: A Dugdale cohesive crack interacting with a number of dislocations. (a) The crack lies in an isotropic elastic medium and is subjected to a remote mode-I loading K I∞ . The crack tip cohesive zone size is L . N edge dislocations with Burgers vectors (b cosφ , b sin φ ), j = 1K N , are located at positions (R ,θ ) with respect to the cohesive j j j j zone tip O . The cohesive zone starts at O ′ where the crack-opening displacement is calculated and compared to δ c as the condition for crack initiation. (b) Dugdale’s xviii (rectangular) cohesive law with cohesive strength σ c , work of separation Γ , and critical displacement δ c = Γ σ c ………………………………………………………………………..87 6.2: Dislocation shielding ΔKˆ I of a Dugdale (rectangular) cohesive crack by (a) a single edge dislocation in front of the crack-tip and (b) two edge dislocations symmetrically located on a slip-plane inclined at ± 60o with respect to the cohesive zone tip and crack plane. The results are plotted as a function of the normalized cohesive strength σˆ c for different distances R = 5nm, 10nm, 100nm between the dislocation and the cohesive zone tip O . For high cohesive strengths, the cohesive crack behaves like a singular crack and the selected dislocation structures anti-shield the crack, while for very low cohesive strengths, the same dislocation structures act to shield the cohesive crack with the value evaluated through Equation (6.16) [103]……………………………………………….. 96 6.3: The ratio between the actual cohesive zone length L and the characteristic cohesive zone length LD defined in Equation (6.10) for a Dugdale (rectangular) cohesive crack in the presence of a single edge dislocation in front of the crack-tip. The results are plotted as a function of the normalized cohesive strength σˆ c for different distances R = 5nm, 10nm, 100nm between the dislocation and the cohesive zone tip O . Interestingly, L remains comparable in magnitude to LD for the full range of cohesive strengths and approaches LD in the high strength limit [103]…………………………. 97 6.4: Variation of the net shielding ΔKˆ I of two dislocations in front of a crack-tip with cohesive strength σ c = 1GPa as a function of distance R2 of the second dislocation, while the first dislocation remains fixed at a distance R1 = 10 nm . Here, ΔKˆ exact is the exact solution, ΔKˆ appx is calculated by superposing toughness contributions from each xix individual dislocation treated separately and ΔKˆ LT is the prediction from the Lin- Thomson solution [103]………………………………………………………………… 98 6.5: Simulated dislocation shielding ΔKˆ of a trapezoidal cohesive crack of Tvergaard and Hutchinson [54] at various cohesive strengths σ c . The Lin-Thomson solution predicts a net shielding of ΔKˆ LT = 0.2662 . Although tends to at very high , significant deviations exists even at = 3.6GPa. The DD simulation results [103] are compared with the asymptotic interpolation solution derived in the Appendix G, with the net shielding obtained by adding contribution from individual dislocations given in Equation (G.6), where the cohesive zone size L is taken to be the characteristic cohesive zone size LD (Equation 6.10). The elastic material is assumed to have Young's modulus E = 140GPa, Poisson's ratio ν = 0.33, and plane strain modulus E * = 157GPa. The TH trapezoidal cohesive law has constant energy Γ = 1.125J/m2, δ 1 = 0.25δ c , δ 2 = 0.5δ c and fracture toughness KI,c = E *Γ = 0.42 [103]…………………………………………………. 101 B.1: A periodic array of edge dislocations with Burgers vector (1, 0) located at position (± np, η ) ( n = 0, 1, 2, .... ) in an isotropic elastic strip. (a) and (b) depict two sub-problems used to calculate the stress field for the original problem by the superposition principle. ……………………………………………………………………………….113 E.1: A circular prismatic dislocation loop of unit Burgers vector in the z-direction and radius R with center located at position (0, z 0 ) in an isotropic elastic cylinder. (a) and (b) depict two sub-problems used to calculate the stress field for the original problem by the superposition principle………………………………………………………………….124 xx G.1: Comparison of the numerically calculated shielding ΔKˆ I of a Dugdale (rectangular) cohesive crack (Equation 6.12) by a single edge dislocation in front of the crack-tip and the corresponding asymptotic interpolation solution (Equation G.6). The distance between the dislocation and the cohesive zone tip is taken to be R = 10 nm , while the material properties are E ∗ = 157 GPa , b = 0.25 nm and Γ = 2 J m 2 [103]……………...138 xxi Chapter 1 O1O Overview and Introduction Materials with small dimensions have become ubiquitous in technology and are increasingly driving research towards simultaneous realization of further size reduction, improved performance and better understanding of material behavior in the presence of surfaces, interfaces, and constraints. The application of such small dimensional materials is to be seen everywhere. Material selection in such systems is usually based on electronic, magnetic or optical properties of material and the design is often decided by the functionality. Ultrathin, submicron Copper/Aluminum films on Silicon substrate is one such system, now widely used in microelectronic circuits. Although electronically suitable, this system often suffers from mechanical stresses through a variety of mechanisms such as intrinsic stress during deposition, thermal mismatch, lattice mismatch and interface stresses [1]. Continued size reduction is accompanied with film stresses which exceed the material limiting values causing inelastic deformation or permanent failure through cracking. The inelastic deformation manifests itself in the form of dislocation activity within grains and diffusional flow between free surface and grain boundaries (GBs) of the film [2-5]. These can have an undesirable effect on the 1 2 performance and structural integrity during service and warrants research and improved understanding. In this thesis-work, we focus on three systems: (i) thin film deposition with mass transport between free surface and GBs as the dominant mechanism, (ii) stress evolution and subsequent pulverization in Lithium (Li) alloy battery electrodes during charging- discharging operations, based on the interstitial diffusion mechanism, and (iii) dislocation shielding of a cohesive crack. 1.1 Volmer-Weber thin film growth Our interest here is in studying the mechanism of diffusional flow between free surface and grain boundaries (GBs) in the context of stress development during thin film deposition. The broader aim behind the study is to enhance our understanding of the stress generating mechanisms during film growth process. Gao and co-workers [6-10] have shown that constrained diffusional creep in thin films causes material to be transported from the film surface into the GBs to form GB diffusion wedges, and that a crack-like singular stress field develops near the roots of the GB diffusion wedges [6-7]. The constrained GB diffusional creep mechanism captured the stress-temperature behavior of unpassivated thin copper films on substrate [11] and provided so far the only feasible explanation for experimentally observed glide dislocations on planes parallel to the film-substrate interface [12]. Compared to thermal cycling of as-grown films, thin metal films grown by the Volmer-Weber (V-W) mode can exhibit even more complex stress histories, which can be detected by measuring the substrate curvature and using Stoney’s formula [13]. The nature of stress in the film is related to three stages of V-W growth: (i) nucleation and growth of discrete islands; (ii) coalescence of existing islands 3 and formation of GBs leading to a continuous film; (iii) thickening of the continuous film with complete coverage of the substrate surface. For high mobility materials, the three stages of V-W growth have been associated with experimentally observed initial compressive stress, tensile stress and final compressive stress regimes as shown in Figure 1.1 [14-17]. deposition on deposition off tensile maximum II III Stress*thickness I time thickness interruption post-coalescence incremental tensile increase compressive stress Figure 1.1: Schematic illustration of stress evolution during high mobility Volmer-Weber growth – interruption. Surface tension effects and adatom-substrate interactions are considered as possible mechanisms for initial compressive stress during the first stage of V-W growth [18-21]. The intermediate tension during the second stage of V-W growth has been linked to island coalescence [22-24]. Origins of the final compressive stress during the third stage of V-W growth are still being debated, with various proposed mechanisms attributed to surface stress effect due to change in adatom populations during growth [20-21], trapping of adatoms in interstitial sites near step edges/ledges [25] and excess chemical potential at the film surface due to high adatom concentrations during growth [26-28]. In Chapters 3-4, using the mechanism of excess surface chemical potential during growth [26-28] in 4 conjunction with constrained GB diffusional creep model [6], we will study the effects of GB mobility on stress evolution during the third stage of V-W growth. The predictions are compared with experiments, with results in support of the proposed mechanism of “mass transfer between surface and GBs” for stress evolution during stage three of the V- W thin film growth. 1.2 Stresses and cracking in Li battery electrodes Figure 1.2: SEM micrographs of germanium films after (a) 1, (b) 10, and (c) 110 discharge/charge cycles and optical micrograph of the sample after 300 cycles (d) [29]. The ever dwindling fossil-fuel supplies, related pollution effects and projected future energy needs of humanity have combined to make the development of radically improved energy storage systems a worldwide imperative. Rechargeable Lithium-ion (Li-ion) battery cells is one such device used widely as a secondary battery systems for portable electronic devices due to their high energy density, high operating voltages, low self discharge, and low maintenance requirements compared to conventional aqueous rechargeable cells such as nickel–cadmium and nickel metal hydride [30]. Li-ion batteries are based on the classical intercalation reaction during which lithium is inserted into or 5 extracted from both cathode and anode. Significant research has been directed towards finding materials with ever greater capacity for accommodating Li atoms to increase battery efficiency [31]. Graphitized carbon, the most common anode for Li-ion batteries, exhibits relatively small volumetric change, stable working voltage and good cycle performance. However, the chemical compound LiC6 limits the theoretical capacity of Li-C anodes to 372 mAhg-1 [32], which is deemed insufficient for high power applications [33]. Compared to the Li-C anodes, a variety of Li-alloy (LixM, M:= Sn, Si, Ge, Al) anodes show substantially higher theoretical capacity, high Li packing density and safe thermodynamic potential [33]. For example, as a potential candidate for anode in Li-ion batteries, silicon has a theoretical capacity of 4200 mAh g–1 with the formation of Li4.4Si alloy [34], which is substantially higher than that of carbon [29]. However, Si electrodes suffer from serious irreversible capacity and poor cyclability due to huge volume changes associated with the lithium ion insertion/extraction processes. This volume change compounded with high Li packing density often results in fast disintegration (cracking or ‘crumbling’) of the material [29,35-37] (Figure 1.2). This phenomenon, commonly referred to as “decrepitation”, has become a major obstacle for practical applications of Si and other high capacity materials in rechargeable Li-ion batteries. Existing strategies to prevent decrepitation of Si have mainly focused on using composite materials and reducing the alloy particle size. In the former approach, an electrochemically active phase is homogeneously dispersed within an electrochemically inactive matrix [33,38-39], with the inactive phase designed to accommodate the large strains generated by the active phase while maintaining the structural integrity of the composite electrode during the alloying/de-alloying processes. In the latter approach, the 6 emphasis is on reducing the size and experimenting with the geometry of the electrode. It has been shown that cracking of LixSn can be avoided in multiphase anodes with small particle sizes [40]. Nanostructured Si and Ge thin films exhibited superior performance during charge/discharge cycling compared to their bulk counterparts [41-42]. Anodes made of Cu nanorods plated with Fe3O4 showed an improvement by a factor of six in power density over planar electrodes while maintaining the same discharge time [43]. Recently, silicon nanowire electrodes of diameter less than 100nm have been found to accommodate volumetric strains as large as 400% without pulverization, while maintaining a discharge capacity close to 75% of its theoretical capacity with little fading during cycling [44]. These experiments are strongly suggesting that size reduction is an effective strategy in creating fracture resistant, high capacity battery electrodes. In Chapter 5, we will model the insertion/extraction of Li during galvanostatic (constant current) operation in an electrode as diffusion of interstitial atoms at the continuum level [45-50], while employing the cohesive zone technique [51-57] to investigate crack nucleation in an initially crack-free electrode due to diffusion induced stresses (DIS). 1.3 Dislocation shielding of a cohesive crack Permanent failure of thin films through cracking or void formation as a stress relieving mechanism often occurs at the interfaces and has motivated much work on crack nucleation/propagation and understanding interfacial toughness of known critical interfaces in these systems [1]. On the modeling front, cohesive models of fracture originally developed to remove the crack tip stress singularity by accounting for cohesive interactions [51] or plastic deformation [52] have gradually evolved as an attractive tool to model crack initiation and propagation without an ad hoc failure criterion [53-57]. 7 More recently, cohesive models of fracture have been increasingly used to study crack propagation in cases where plasticity is accounted for by modeling individual discrete dislocations. This capability is especially important for modeling crack nucleation and growth in confined small scale structures like thin films on substrates and/or under cyclic fatigue conditions caused by thermal cycling. Discrete dislocation (DD) plasticity has a natural length scale (the Burgers vector) which is absent in conventional length-scale independent continuum plasticity models. Many experimentally observed trends have been captured by DD-cohesive fracture models [58-65]. The effect of dislocations to reduce the stress intensity factor of a crack tip is referred to as shielding and analytical expressions are derived for a mathematically singular crack by Lin and Thomson [66] and Weertman [67]. A natural question arises as to whether the concept of dislocation shielding of a singular crack is directly applicable to that of a cohesive crack and what the corresponding conditions are. If different, how does it impact the use and interpretation of cohesive models in modeling crack initiation and propagation? In Chapter 6, we will make an attempt to clarify the above questions through consideration of a Dugdale cohesive crack interacting with dislocations under mode I dominant conditions. 1.4 Organization of this thesis The rest of the thesis is organized as follows. In Chapter 2, we present the basic formulation used throughout the thesis. In Sections 2.1 and 2.2 we discuss integral equation methods and cohesive zone models. In Section 2.3, diffusion at interfaces viz. GBs, and bulk diffusion is illustrated in relation to the stress in the system. In Chapter 3, we formulate the V-W growth problem under constrained GB diffusion model. In Section 3.1, we summarize the experimental observations, and set up the model through 8 appropriate governing equations, initial conditions and boundary conditions in Section 3.2. In Sections 3.3 and 3.4, the model is specialized to finite-inhomogeneous or infinite GB mobility with numerical algorithm and/or analytical solutions. In Chapter 4, the results are presented based on the numerical algorithm and the analytical solutions of the previous Chapter. Section 4.1 deals with a finite-inhomogeneous GB mobility, with Sub- Sections 4.1.1 and 4.1.2 devoted to the formation of GB diffusion wedges and compressive stress evolution during thin film growth, respectively. In Section 4.2, the infinite GB mobility solution is compared with the Sn film deposition and growth interruptions experiments [28]. In Chapter 5, crack nucleation studies are performed in strip and cylindrical electrodes due to diffusion induced stresses (DIS) during galvanostatic opearation. In Section 5.1, we discuss the background for the study and in Section 5.2, DIS is evaluated for both strip (SubSection 5.2.1) and cylindrical geometries (SubSection 5.2.2) under galvanostatic boundary conditions. Section 5.3 sets up the cohesive model for crack nucleation and condition for localization spacing in the electrode. In Section 5.4, the critical conditions for crack nucleation are discussed for thin strip and cylindrical electrodes. In Chapter 6, the problem of dislocation shielding of a cohesive crack is studied. In Section 6.1, the dislocation shielding of a mode I Dugdale is analyzed. In Sub-Sections 6.1.2, the limiting solutions under very high and low cohesive strengths are derived and numerical analysis covers the intermediate cohesive strength regime in Sub-Section 6.2.2 through a couple of representative examples. In Section 6.2, the asymptotic formula (Appendix G) based on the analytical solutions of the previous Section are compared with the prediction of DD simulations 1 for a large number of 1 DD simulations are carried by Dr. Audrey Chng, of Institute of High Performance Computing, Singapore. 9 dislocations interacting with a cohesive crack governed by the Trapezoidal traction- separation law [54]. Appendices A-G, provide supplementary information dealing with details of numerical solutions encountered in Chapters 3-6. Finally in Chapter 7, the most important results of this thesis are summarized and possible avenues for future research are discussed. Chapter 2 Basic 2. mathematical techniques and Formulation Stress-diffusion coupling and cohesive zone models in fracture mechanics are two central concepts underlying all the problems attempted in the thesis. The models developed with these concepts are formulated using the well-established integral equation methods in fracture mechanics. The first three sections provide an introduction to these concepts, methods, and solution techniques, which serve as a background for the particular problems covered in subsequent chapters. 2.1 Integral equation method in fracture mechanics Integral equation method is one of the most widely used techniques to solve crack problems in fracture mechanics [68]. Its attractiveness over alternative numerical schemes such as finite element method stems from the ability to deal with mixed- boundary value problems, efficiency (it reduces a partial differential equation in two- dimension to a one-dimensional integral equation) and accuracy. 10 11 We first demonstrate the general procedure of formulating and solving crack problems using integral equation method through a finite straight crack (Griffith crack) of length 2a in an infinite elastic isotropic plane and subjected to a far-field tension σ ∞ as shown in Figure 2.1. The problem undertaken here is one of the simplest in its class with a well-known analytical solution [69], but is comprehensive enough to demonstrate the ideas such as dislocation density, Green’s function, Cauchy kernel, and singular integral equation. σ∞ σ∞ σ∞ y y x = + x 2a 2a σ∞ 2a σ∞ σ∞ (a) (b) Figure 2.1: Crack of length 2a in an infinite plane with a far field stress σ ∞ . The problem can be solved by superposition of a uniform stress and a crack with no far field loading but with crack face pressure equal to σ ∞ . Due to the completely linear nature of the problem, we can apply superposition and split the original problem into a plane with no crack (Figure 2.1a) and a crack in the plane subjected to pressure (Figure 2.1b). The key idea behind the split is to obtain a problem with a known traction distribution on crack faces. Next we make use of dislocations, another well-known defect [70] to account for the crack opening due to the applied tractions on its face. This idea to represent crack opening using dislocations was 12 pioneered by Bilby and Eshelby [71] and has since then being developed and used widely in the fracture mechanics community [67-68]. The dislocations considered here are fictitious i.e. they don’t actually represent plastic deformation and the Burgers vector magnitude varies continuously with position as opposed to the real Burgers vector whose magnitude and orientation are dictated by the crystal structure in the metals. The gradient of Burgers vector is called dislocation density and is denoted by B(ξ ) in the formulation below. Thus the crack problem in Figure 2.1b is transformed to calculating an appropriate dislocation density that accounts for the traction on its faces. Finally in order to relate the dislocation density to the traction on the crack faces, we make use of the Green’s function solution: the stress field produced at (x, y ) by a dislocation with unit Burgers vector located along the crack line in the medium of interest. In Figure 2.1, the Green’s function along the x-axis due to a dislocation with Burgers vector (0,1) located at (ξ ,0 ) along the crack line in an infinite isotropic elastic plane is given as [70] E 1 G ( x, ξ ) = . (2.1) ( ) 4π 1 − v ( x − ξ ) 2 Using linearity of the problem, the pressure along the crack faces in Figure 2.1b can be expressed as a ∫ G(x, ξ )B(ξ )dξ = −σ −a ∞; −a≤x≤a, (2.2) or a 1 B(ξ )dξ = − ( 4π 1 − v 2 ) ∫ −a (x − ξ ) E σ∞; −a≤x≤a. (2.3) The term (x − ξ )−1 in the Green’s function is called Cauchy kernel and the point x = ξ is referred to as the singular point (i.e. as ξ → x , the Cauchy kernel tends to infinity). This 13 singular nature of the Green’s function is common to all crack problems and is an outcome of the singular nature of the dislocation stress field. Equation (2.3) is hence called as the singular integral equation with unknown dislocation density B(ξ ) . In general the singular integral equation has the following form, a a 1 ∫ B(ξ )dξ + k (x, ξ )B(ξ )dξ = g (x ); ∫ −a≤x≤a. (2.4) −a (x − ξ ) −a The function k (x, ξ ) is the smooth-bounded part of the Green’s function and g (x ) is a known function that depends on the tractions on the crack surface. Before moving further, we have to take into account the singularity of B(ξ ) at the crack-tips [68]. For Figure 2.1a, B(ξ ) has singularities at both crack ends and hence can be expressed as R(ξ ) 1 − ξ 2 a 2 where R(ξ ) is a smooth function in ξ . Thus Equation (2.3) can be rewritten as, a 1 R(ξ ) ( 4π 1 − v 2 ) ∫ −a (x − ξ ) 1 − ξ 2 a 2 dξ = − E σ∞; −a≤x≤a (2.5) Even though the Cauchy term becomes infinite at ξ = x and the dislocation density itself becomes infinite as ξ → ± a , the integral in Equation (2.5) is well-defined in terms of the Cauchy Principal Value. This is a mathematical reflection of the fact that stress field in the interior of the crack face is finite. The crack-opening displacement can be expressed in terms of displacement density as, R(ξ ) a a Δu (x ) = ∫ B(ξ )dξ = ∫ dξ −a≤x≤a, (2.6) x x 1− ξ 2 a2 where we have made use of the fact that crack-opening displacement is zero at the end x = a . This places a constraint on the dislocation density, 14 R(ξ ) a a ∫ B(ξ )dξ = ∫ −a −a 1− ξ 2 a2 dξ = 0 . (2.7) Equations (2.5-2.7) are normalized over the interval [-1,1] to facilitate numerical () [ ] integration. Thus with xˆ → x a , ξˆ → ξ a and Rˆ ξˆ → ER(ξ ) σ ∞ (1 − ν 2 ) , Equations (2.5- 2.7) become, 1 1 Rˆ ξˆ () ∫ −1 ( ˆ x − ξˆ) 1−ξ ˆ 2 dξˆ = −4π ; − 1 ≤ xˆ ≤ 1 , (2.8) EΔu ( xˆ ) 1 () R ξˆ ( σ ∞ 1 −ν 2 a ) ∫ = xˆ 1−ξ ˆ 2 dξˆ − 1 ≤ xˆ ≤ 1 , (2.9) 1 () Rˆ ξˆ ∫ −1 1−ξ ˆ 2 dξˆ = 0 . (2.10) Equation (2.8) can be inverted in closed form, but this is usually difficult for generally encountered crack problems and we often have to resort to numerical methods. The numerical treatment of Equations (2.8-2.10) is done through discretization of the interval [-1,1] into finite number of points N . Equation (2.8) is imposed at N − 1 points xˆ i , i = 1,..., N − 1 called as the collocation points and the integrals are evaluated using the points ξˆ j , j = 1,..., N called as the integration points. All the existing numerical integration methods (quadrature) are based on such discretization and differ only through the respective location of the points [72]. Amongst these, the Gaussian quadrature methods are found to be most efficient in terms of accuracy for the same number of points N [72]. The presence of square-root singularity doesn’t present any numerical difficulty as it is treated as a weight in the Gauss-Chebyshev quadrature [68,72] and can be very accurately integrated depending on the number of points N and behavior of 15 () Rˆ ξˆ . For the integrals in Equations (2.8-2.10), the collocation and integration points are xˆ i = cos(iπ N ), i = 1,..., N − 1 and ξˆ j = cos(( j − 0.5)π N ), j = 1,.., N , respectively, based on the Gauss-Chebyshev quadrature, leading to the integration rule f (ξˆ) 1 π N ∫ dξˆ ≅ N ∑ f (ξˆ ) , j (2.11) −1 1 − ξˆ 2 j =1 Thus the discretized version of Equations (2.8) and (2.10) is π N ( ) Rˆ ξˆ j ∑ N j =1 ( xˆ i − ξˆ j ) = −4π ; xˆ i , i = 1,..., N − 1 (2.12) ∑ Rˆ (ξˆ ) = 0 , π N j (2.13) N j =1 which together comprises a system of N × N linear equations in the unknowns Rˆ ξˆ j . ( ) ( ) Once the values of Rˆ ξˆ j , j = 1,..., N are known, the crack-opening displacement is calculated based on Equation (2.9) as, EΔu (xˆ ) ∑ k (xˆ, ξˆ )R(ξˆ ) , π N = (2.14) ( σ ∞ 1 −ν a N 2 ) j =1 j j where ⎧⎪0 − 1 ≤ ξˆ < xˆ ( ) k xˆ , ξˆ = ⎨ . (2.15) ⎪⎩1 xˆ ≤ ξˆ ≤ 1 The stress-field of the original problem in Figure 2.1 is given as a σ ij (x, y ) = σ ija + ∫ Gij (x, y, ξ )B(ξ )dξ , (2.16) −a where σ ija = σ ∞ δ i 2δ j 2 is the solution of part (a) and Gij (x, y, ξ ) is the ij th stress component at (x, y ) due to a dislocation at (ξ ,0 ) with Burgers vector (0,1) in an infinite isotropic 16 elastic plane [70]. Another very important quantity in Linear Elastic Fracture mechanics is the stress-intensity factor K , which can be expressed in terms of the crack-opening displacement gradient using Irwin’s asymptotic relationship [69] ⎛ ∂Δu KI = E Lim 2π (a − x )⎜ − (x,0)⎞⎟ = E 2 Lim 2π (a − x )B(x,0 ) , (2.17) x=a ( 4 1 −ν 2 x →a ) ⎝ ∂x ( ⎠ 4 1 − ν x→a ) or in terms of normalized variables, KI Rˆ (1) = . (2.18) σ∞ π a x=a 4 Similarly, the stress-intensity factor at the other end can be expressed as, KI Rˆ (− 1) =− . (2.19) σ∞ π a x=−a 4 ( ) The values of Rˆ (± 1) can be calculated using the known values Rˆ ξˆ j , j = 1,.., N by Krenk’s interpolation formula [73] sin[(N − 0.5)( j − 0.5)π N ] ˆ ˆ ( ) N 1 Rˆ (1) = N ∑ j =1 sin [( j − 0.5)π (2 N )] Rξj , (2.20) sin[(N − 0.5)( j − 0.5)π N ] ˆ ˆ ( ) N 1 Rˆ (− 1) = N ∑ j =1 sin [( j − 0.5)π (2 N )] R ξ N +1− j . (2.21) The nature of singularity for the dislocation density depends on many factors like the physical situation or the type of crack i.e. centre/edge under investigation. Based on the singularity, the location of integration-collocation points ξˆ j - xˆ i and the expressions in Equations (2.11, 2.20-2.21) will change and are noted accordingly in the Appendices. In Chapter 3, we will encounter a Green’s function kernel with singularity 1 ( x − η ) , which 3 is stronger than Cauchy singularity due to the exponent being greater than 1. The methods of this section will require modification and are discussed in Appendix A. 17 In summary, an integral equation applied to crack problems requires finding the traction distribution on crack faces, appropriate Green’s function (which depends on the geometry and material properties of the medium) and identifying the singularity behavior of the dislocation density distribution. Once these things are in place, numerical methods [68] described here can be applied for the solution of a given problem. δ = Separation distance σ (δ ) σc σ (δ ) σ∞ σ (δ ) Γ δ δ L δ c c δt = Tip opening displacement δc 2a (a) (b) Figure 2.2: (a) Ahead of a crack tip, the surface separation δ is related to the cohesive stress σ (δ ) with the maximum cohesive stress σ c referred to as the cohesive strength of material. When the tip opening displacement δ t is equal to the maximum cohesive interaction distance δ c , crack growth is assumed to occur. Γ is the cohesive energy and the cohesive zone size is L. (b) Cohesive zone technique applied to the centre crack problem of Figure 2.1. Cohesive zone acts in the region [− a,−a + c ] and [a − c, a ] with a given cohesive zone law σ = σ (δ ) . 2.2 Cohesive zone models The original Griffith theory of fracture ignores the unphysical prediction of infinite stress at a crack tip and invokes thermodynamics through energy balance to obtain a fracture criterion [74]. Barenblatt [51] proposed an alternative, where the surface just ahead of the crack interacts through atomic or molecular cohesive forces opposing separation, as can be seen from Figure 2.2a. The relationship between the cohesive forces σ and the δc separation δ is called the cohesive law with cohesive energy Γ = ∫ σ (δ )dδ . Crack 0 extension occurs when the crack tip opening δ t equals the maximum cohesive interaction 18 range δ c . The maximum cohesive stress is called the cohesive strength σ c and bounds the stress in the vicinity of the crack-tip. First principle calculations can be used to infer the exact form of the cohesive law σ = σ (δ ) . Nevertheless, it is often found that, once the cohesive energy Γ and the cohesive strength σ c are given, the exact profile of the cohesive law is immaterial. In literature, the popular cohesive laws are Dugdale/rectangular [52], trapezoidal [54], polynomial/exponential [53], triangular [75]. Historically, Dugdale [52] proposed the cohesive zone model around the same time as Barenblatt, but in an entirely different context. In an effort to model plastic deformation ahead of a crack-tip in a metal sheet, Dugdale assumed that yielding is confined in a narrow zone directly ahead of the crack tip and analyzed the model by viewing the effect of yielding as making the crack longer by an amount equal to the plastic zone size L , with cohesive stresses in the plastic zone acting on the extended crack surface so as to restrain the opening. In his picture, the cohesive strength σ c corresponds to the yield strength σ y of material. Both the applied load and the cohesive stresses create inverse square root singularities at the outer tip of the plastic zone, but these singularities are of opposite sign and the zone size L is chosen so that they cancel each other and bounded stresses result at the outer tip. Cohesive models of fracture have also been extensively used to model crack nucleation in an initially crack-free material [53-57]. In such models, a traction-separation law is assumed to act between separating surfaces and the criterion of crack nucleation is based on the maximum surface separation within the cohesive zone reaching a critical value. In the following, we will modify the integral equations developed in the previous section to deal with cohesive cracks through the example of centre crack in an infinite 19 plane discussed previously (Figure 2.2b). The crack surfaces now interact through cohesive zones at the crack ends with size c and cohesive zone law σ = σ (δ ) between the surfaces. Making use of the superposition principle to account for the cohesive stress and noting that the surface separation δ = Δu (x ) (Equation (2.6)), the traction condition of Equation (2.3) becomes ( ⎧ 4π 1 − v 2 ) ⎡⎢σ ⎛a ⎞⎤ + σ ⎜ B(ξ )dξ ⎟⎥ ⎪− ⎪ E ⎢⎣ ∞ ∫ ⎜ ⎝x ⎟⎥ ⎠⎦ − a < x < −a + c ⎪ a 1 ( ⎪ 4π 1 − v 2 B(ξ )dξ = ⎨− )σ ∫ −a (x − ξ ) ⎪ E ∞ −a+c< x 1 ) are called hypersingular and special numerical techniques are needed n to solve such problems [84]. Details of the numerical method for solving Equations (3.24) and (3.25) are given in Appendix A. 3.4 Infinite GB diffusivity In the limit of infinitely fast GB diffusion, atomic flux in a GB is assumed to be fast enough to eliminate any gradient in the normal stress σ ( z , t ) along the GB so that we can treat σ ( z , t ) as a variable which only depends on time, i.e. σ ( z, t ) = σ (t ) . Since the GB diffusivity δ GB DGB tends to infinity and the GB traction gradient tends to zero, we cannot directly evaluate the flux into GB from Equation (3.9). In this case, it is more convenient to equate the net increase of atomic mass in the GB to the net flux into the GB, 2Ω 5 / 3C s Γs [σ (t ) − σ sd ] h(t ) d 2 dt ∫ u(ς , t )dς = 0 kT . (3.28) Equations (3.6) and (3.28) are to be solved with the boundary and initial conditions given by Equations (3.16b) and (3.18) respectively. 3.4.1 Analytical solution The solution to equation (3.6) has the general form σ 0 − σ (t ) ⎛z h⎞ u( z, t ) = ∗ L f⎜ , ⎟ (3.29) E ⎝h L⎠ where f (z h , h L ) defines the GB displacement profile with f (1, h L ) = 0 in accordance with Equation (3.16b). Combining Equations (3.28) and (3.29), we obtain the governing 36 equation for GB stress evolution as ⎡ Ω1 / 3 ⎤ & + Ω σ sd , 1/ 3 σ& (t )h(t ) + ⎢(1 + δ )h& + σ ( t ) = σ (1 + δ )h (3.30) 2λt D ⎥⎦ 2λ t D 0 ⎣ where LkT tD = ∗ (3.31) 2 E Ω 4 / 3Cs Γs represents the time for the insertion of an atomic layer into the GB during deposition and 1 λ (h L ) = ∫ f ( zˆ, h L )dzˆ , (3.32) 0 h λ ' (h L ) δ (h L ) = (3.33) L λ (h L ) are two auxiliary universal functions that only depend on the ratio between the film thickness and grain size. To determine the universal function λ (h L ) , we insert Equation (3.29) into Equation (3.6) and obtain an integral equation 1 1 P( z L , hςˆ L) f ' (ςˆ, h / L )dςˆ = 1 2π ∫0 (3.34) over 0 ≤ z ≤ h . For a given value of h L , this integral equation can be solved to determine the crack-like GB displacement profile f (z h , h L ) and then to obtain λ (h L ) by integrating the result according to Equation (3.32). Figure 3.3a plots the calculated universal function λ (h L ) along with a fitting function ⎡ ⎛h⎞ 5 ⎛h⎞ 4 ⎛h⎞ 3 ⎛h⎞ 2 ⎛h⎞ ⎤ λ = 0.5 − exp⎢− 0.0344⎜ ⎟ + 0.408⎜ ⎟ − 1.85⎜ ⎟ + 4.1⎜ ⎟ − 5.04⎜ ⎟ − 0.745⎥ . (3.35) ⎢⎣ ⎝L⎠ ⎝L⎠ ⎝L⎠ ⎝L⎠ ⎝L⎠ ⎥⎦ It is seen that λ (h L ) rises quickly within the range 0 ≤ h L ≤ 1 2 and asymptotically approaches its limiting value λ (∞ ) = 1 / 2 . For example, λ (1 2) = 0.416 , which is already 37 within 20% of the limiting value λ (∞ ) = 1 2 . In contrast, the asymptotic limit of δ (h L ) is δ (∞ ) = 0 . Figure 3.3b shows that δ (h L ) can be fitted with the function ⎡ ⎛h⎞ 5 ⎛h⎞ 4 ⎛h⎞ 3 ⎛h⎞ 2 ⎛h⎞ ⎤ δ = exp ⎢− 0.048⎜ ⎟ + 0.56⎜ ⎟ − 2.5⎜ ⎟ + 5.4⎜ ⎟ − 6.27⎜ ⎟ + 0.43⎥ . (3.36) ⎢⎣ ⎝L⎠ ⎝L⎠ ⎝L⎠ ⎝L⎠ ⎝L⎠ ⎥⎦ Once the universal functions λ (h L ) and δ (h L ) are calculated and the film growth history h(t ) is given, Equation (3.30) can be integrated to obtain the GB stress evolution under the initial condition σ (0) = σ 0 . Numerical Equation (3.36) Numerical Equation (3.35) Figure 3.3: Two auxiliary universal functions in the analytical model of compressive stress evolution during growth of a thin film. The functions (a) λ (h L ) and (b) δ (h L ) along with their fitting function given in Equations (3.35) and (3.36) respectively [85]. In the experiments by Shin and Chason [28], the aspect ratio between the film thickness and grain size is larger than 1 2 . For h L ≥ 1 2 , we have 0.5 ≥ λ ≥ 0.416 and 0.178 ≥ δ ≥ 0 . Our numerical solutions (next chapter) indicate that the variations of λ and δ can be neglected in this range so that λ (h L ) and δ (h L ) can be approximated by their corresponding limiting values λ (∞ ) = 1 2 and δ (∞ ) = 0 , so that Equation (3.30) is simplified as 38 ⎡ Ω1 / 3 ⎤ & + Ω σ sd . 1/ 3 σ& (t )h(t ) + ⎢h& + σ ( t ) = σ h (3.37) t D ⎥⎦ 0 ⎣ tD Under a constant deposition rate h& (Equation (3.12)), the solution to Equation (3.37) is ⎡ Ω1 / 3 ⎤ ⎢1+ & ⎥ ⎛ h0 ⎞ ⎣⎢ ht D ⎦⎥ σ (t ) = (σ 0 − σ ss )⎜ ⎟ + σ ss , (3.38) ⎝ h ( t ) ⎠ where σ 0 + Ω1/ 3σ sd t D h& σ ss = (3.39) 1 + Ω1/ 3 t D h& is the steady state stress during deposition. This expression is actually identical to that derived by Chason et al. [26] based on a different set of assumptions and shows that the steady-state stress depends on the deposition rate, with value ranging from σ 0 to σ sd . In the model of Chason et al. [26], GBs are assumed to slide freely along the film-substrate interface whereas here the GB displacement is constrained at the film-substrate interface in accordance with the constrained GB diffusion model [6]. At high deposition rates, the film stress is equal to σ ss = σ 0 as there is insufficient time for adatom diffusion into the GB. At low deposition rates, we have σ ss = σ sd as the adatom diffusion has plenty of time to fully equilibrate between the surface and GB. When growth is interrupted at a specified film thickness h , the stress evolution equation becomes Ω1/ 3 Ω1/ 3 σ& (t )h + σ (t ) = σ si , (3.40) tI tI and the solution is 39 ⎛ Ω1 / 3 t ⎞ σ (t ) = (σ 0i − σ si ) exp⎜⎜ − ⎟⎟ + σ si , (3.41) ⎝ tI h ⎠ where σ 0i is the stress state at the moment deposition is stopped and LkT tI = ∗ (3.42) 2 E Ω 4 / 3C g Γg denotes the time required for the removal of an atom layer out of the GB when growth is interrupted. Thus the analytical model predicts exponential stress recovery during growth interruption with time constant tI h tR = , (3.43) Ω1 / 3 which evidently depends on the film thickness. Equation (3.38) and (3.41) together describe stress evolution during a cycle of deposition and growth interruption for thin film growth in the limit of infinite surface and GB diffusivities. 3.4.2 Numerical solution  Numerical methods can also be used to solve Equations (3.6) and (3.28) for a given history of film growth. In numerical calculations, it will be convenient to normalize the variables in the present problem as z ς t h(t ) & & t D u ( z, t ) x= ,η = , τ = , H (τ ) = , H = h 1 / 3 , uˆ ( x,τ ) = . (3.44) h h tD h0 Ω 2π h0 (σ 0 − σ sd ) E ∗ For a constant growth rate, Equation (3.12) is normalized as ⎛ H& Ω1 / 3 ⎞ H (τ ) = 1 + ⎜⎜ ⎟⎟τ . (3.45) ⎝ h0 ⎠ 40 Combining Equations (3.6) and (3.28) in terms of the normalized variables leads to 1 1 h H (τ ) ∂ 2 uˆ (η ,τ ) L ∂uˆ (η ,τ ) L Ω ∫ H& uˆ (0,τ ) − 0 1 / 3 η 0 ∂η∂τ dη + ∫ 4πh0 H (τ ) 0 P ( x,η ) ∂η dη = 4πh0 , (3.46) which is an integro-differential equation in ∂uˆ (η ,τ ) ∂η . The time derivatives are resolved through finite difference and the integrals are evaluated using Gauss-Chebyshev quadrature. Further details regarding the numerical method for solving such type of integral equation can be found in Appendix A. At each time step, the displacement solution from Equation (3.46) is substituted back into Equation (3.6) to obtain the GB stress. Chapter 4 3. Volmer-Weber thin film growth: Results 4.1 Finite and inhomogeneous GB diffusivity Figure 4.1: Inhomogeneous GB diffusivity profiles represented by the Dirac function g ( z ) = 2.0 [1 + exp(β z h0 )] . The higher the value of β , the more inhomogeneous the GB diffusivity becomes. In this section, we examine the effects of deposition and inhomogeneous GB diffusivity on the formation of GB diffusion wedges and stresses in the film using the algorithm discussed in the Section 3.3. Although the formulation is valid for arbitrarily inhomogeneous grain boundary mobility, here we focus our attention on a class of inhomogeneous GB mobility profiles described by the Dirac function, 41 42 2 .0 g ( z) = (4.1) [1 + exp(β z h0 )] which is shown in Figure 4.1 for different values of β . It can be seen that g (z ) decays with depth z into the grain boundary depending on the value of β , with β = 0 corresponding to the case of homogeneous GB diffusivity. Thus β represents the degree of GB inhomogeneity: the higher the value of β , the stronger the GB inhomogeneity. 4.1.1 Formation of GB diffusion wedges during thin film growth To investigate the phenomenon of GB diffusion wedge during thin film growth, we consider a film growing at a deposition rate H& according to Equation (3.23) and monitor the evolution of GB normal traction and GB width at different time instants during the growth. In addition to the characteristic time scale t GB (Equation (3.22)), the growth process introduces an additional time scale represented by 1 H& (Equation (3.21)). The deposition rate H& is found to control the effect of GB diffusion on the GB traction. If H& is smaller or comparable to 1 , GB diffusion is found to strongly influence the GB traction, otherwise the GB stress is largely unaffected and remains tensile. This deposition rate allows sufficient time for atoms to diffuse into the GBs during the growth process. Figure 4.2 shows the GB normal traction and width profile at different time instants of growth for a film of uniform GB diffusivity at a deposition rate of H& = 1 . This deposition rate allows sufficient time for atoms to diffuse into the GBs during the growth process. 43 0.7 0.6 H = 1.16 0.5 0.4 σˆ 0.3 H = 1.32 0.2 H = 1.48 0.1 H = 1.96 H=3 0 0 0.2 0.4 0.6 0.8 1 x (a) 0.16 H=3 0.14 H = 1.96 0.12 0.1 H = 1.48 uˆ 0.08 H = 1.32 0.06 H = 1.16 0.04 0.02 0 0 0.2 0.4 0.6 0.8 1 x (b) Figure 4.2: Formation of a GB diffusion wedge during film growth ( H = 1 + H& τ ) under uniform GB diffusivity ( β = 0 ). Evolutions of (a) the GB traction and (b) the GB wedge profiles at different time instants during deposition for parameter choices of H& = 1 , A = 1 , h L = 0.5 . The 0 stress has been normalized as σˆ = (σ − σ s ) (σ 0 − σ s ) and the GB opening displacement has been normalized as uˆ = uE ∗ [2πh0 (σ 0 − σ s )] , where σ 0 is the initial film stress, σ s is stress corresponding to the excess surface chemical potential during deposition and h0 is the initial film thickness. The time has been normalized as τ = t t GB where t GB is given by Equation (3.22) [86]. Similar to the formation of GB diffusion wedges under annealing in [6], we found 44 that the normal stress along the GB (Figure 4.2(a)) decays with time (and with increasing film thickness in the present case) and eventually equilibrates with σ s , the stress corresponding to excess surface chemical potential, at steady state. The GB diffusion wedge during deposition (Figure 4.2(b)) also resembles a mode I crack profile similar to GB diffusion wedges under annealing conditions [6]. In faster growing films, there is less diffusion into the GBs and less evidence of a diffusion wedge. 0.8 0.6 H = 1.16 σˆ 0.4 H = 1.32 0.2 H = 1.48 H = 1.96 0 H=3 0 0.2 0.4 0.6 0.8 1 x (a) 0.16 0.14 H=3 0.12 H = 1.96 0.1 H = 1.48 uˆ 0.08 H = 1.32 0.06 0.04 H = 1.16 0.02 0 0 0.2 0.4 0.6 0.8 1 x 45 (b) Figure 4.3: Formation of a GB diffusion wedge during film growth ( H = 1 + H& τ ) under strongly inhomogeneous GB diffusivity ( β = 4 ). Evolutions of (a) the GB traction and (b) the GB wedge profiles at different time instants during film deposition for parameter choices of H& = 1 , A = 1 , h0 L = 0.5 . The stress has been normalized as σˆ = (σ − σ s ) (σ 0 − σ s ) and the GB opening displacement has been normalized as uˆ = uE ∗ [2πh0 (σ 0 − σ s )] , where σ 0 is the initial film stress, σ s is stress corresponding to the excess surface chemical potential during deposition and h0 is the initial film thickness. The time has been normalized as τ = t t GB where t GB is given by Equation (3.22) [86]. In Figure 4.3, we consider the effect of inhomogeneous GB diffusivity on the evolution of GB diffusion wedges during growth. Figures 4.3(a) and 4.3(b) shows the GB traction and width profile for a thin film growing with a deposition rate of H& = 1 and an inhomogeneous GB diffusivity parameter β = 4 . We observe that, similar to the homogeneous GB diffusivity case, the GB stress (Figure 4.3(a)) decays with increasing film thickness, but the decay rate is slower due to the lower average GB mobility. After the film thickness is doubled via deposition, significant portion of the GB has been equilibrated with σ s , but the portion of GBs closer to the substrate still retains a significant fraction of the initial stress. The GB diffusion wedge (Figure 4.3(b)) also adopts a crack like opening profile during deposition, similar to the case of uniform GB diffusivity. However, compared to the uniform GB diffusivity case, the GB diffusion wedge exhibits a smaller slope, indicating a lack of diffusion activities along the GBs once atoms are incorporated. Buehler et al. [10] have shown that the stress intensity factor near the tip of GB diffusion wedges can be used as a criterion for dislocation nucleation and demonstrated that the theoretical predictions are in good agreement with the molecular dynamics 46 simulations. Therefore, the stress intensity factor associated with the GB diffusion wedge is a very important indicator of dislocation plasticity in thin films. 0.6 0.5 dH/dτ = 1 0.4 Kˆ 0.3 dH/dτ = 5 0.2 dH/dτ = 10 0.1 0 dH/dτ = 100 -0.1 1 1.5 2 2.5 3 H(τ) Figure 4.4: Evolution of the normalized stress intensity factor near the tip of a GB diffusion wedge with the increasing film thickness under different deposition rates and uniform GB diffusivity ( β = 0 ). Other parameters used in the calculations are A = 1 and h0 L = 0.5 . The ( ) stress intensity factor has been normalized as Kˆ = K σ πh , and the deposition rate has been 0 0 normalized as H& = h& (h0 t GB ) , where σ 0 is the initial film stress, h0 is the initial film thickness and t GB is given by Equation (3.22) [86]. In Figures 4.4 and 4.5, we consider the effects of film growth and inhomogeneous GB diffusivity on the evolution of stress intensity factor at the tip of the diffusion wedge. Figure 4.4 shows the evolution of stress intensity factor at different deposition rates during film growth under uniform GB diffusivity. For the same film thickness, higher deposition rate means less diffusion into the GB and lower stress intensity factor. For H& = 100 , the GB is virtually closed and the stress intensity factor is close to 0 . The stress intensity factor also decreases with increasing inhomogeneity parameter β of GB diffusion (Figure 4.5). This can again be understood from the fact that lower mobility 47 leads to less GB diffusion which in turn leads to reduced stress intensity. We thus conclude that the GB stress intensity factor in a growing film depends on the deposition rate and the GB mobility. For slowly growing high mobility films, the singular stress concentration at the roots of GB wedges is comparable to that of a GB diffusion wedge under annealing conditions and can lead to nucleation of dislocations, as previously shown by molecular dynamics simulations [10] and experiments [12]. The dislocation nucleation can cause additional stress relaxation and affect the measured stress evolution in the film during growth. 0.7 β=0 0.6 β=1 β=4 0.5 0.4 Kˆ 0.3 0.2 0.1 0 1 1.5 2 2.5 3 H(τ) Figure 4.5: Effect of inhomogeneous GB mobility on the evolution of the normalized stress intensity factor near the tip of a diffusion wedge with increasing film thickness at the deposition rate of H& = 1 . Other parameters used in the calculation are A = 1 and h0 L = 0.5 . The stress ( ) intensity factor has been normalized as Kˆ = K σ πh , and the deposition rate has been 0 0 normalized as H& = h& (h0 t GB ) , where σ 0 is the initial film stress, h0 is the initial film thickness and t GB is given by Equation (3.22) [86]. 48 4.1.2 Effect of inhomogeneous GB diffusivity on compressive stress evolution during thin film growth H (τ ) H& 1 1 τ1 τ2 τ Figure 4.6: A film deposition and growth interruption process. The thin film is deposited at a normalized deposition rate of H& = h& (h0 t GB ) from an initial thickness of 1 up to time τ 1 followed by growth interruption till time τ 2 . The time has been normalized as τ = t t GB where t GB is given by Equation (3.22). h0 is the initial film thickness. In order to study the effect of inhomogeneous GB diffusivity on compressive stress evolution during thin film growth, we have conducted numerical simulations, based on the theoretical framework laid down in Section 3.3, for a cycle of film deposition followed by growth interruption as shown in Figure 4.6. We vary the GB inhomogeneity parameter β in Equation (4.1) and study its effect over a range of deposition rates. The results are shown in Figures 4.7(a) and 4.7(b) for slower and faster deposition rates. At a slower deposition rate of H& = 1 in Figure 4.7(a), atoms have sufficient time to diffuse into the GB and the Stress-Thickness plot for β = 1, 2 during deposition is not much different from the uniform GB diffusivity case, while for β = 4 , the strongly inhomogeneous GB diffusivity leads to slower stress relaxation and higher Stress- Thickness values. The effect of GB inhomogeneity is more apparent in the subsequent stress recovery during growth interruption shown in the right panel of Figure 4.7(a). The 49 greater the value of β , the smaller the stress recovery is, which is consistent with the effect of GB mobility on experimentally observed stress relaxation behaviors during deposition-growth interruption [14-16]. For the case of β = 4 , there is virtually no stress recovery, indicating the limit of very low mobility material [14-17]. At the higher deposition rate of H& = 100 , atoms have less time to diffuse into the GB and the Stress- Thickness curves become insensitive to β with a tensile steady state stress, as shown in the left panel of Figure 4.7(b). During growth interruption, the tensile stress in the film is relaxed to various degrees, with smaller values of β leading to faster stress relaxation. (a) 50 (b) Figure 4.7: Evolutions of stress-thickness during one cycle of deposition followed by growth interruption with varying degrees of GB inhomogeneity represented by different values of β in Equation (4.1). The results are shown at (a) a slower deposition rate of H& = 1 and (b) a faster deposition rate of H& = 100 . The deposition rate and time have been normalized as H& = h& (h t ) 0 GB and τ = t t GB respectively, where σ 0 is the initial film stress, h0 is the initial film thickness and t GB is given by Equation (3.22). Other numerical parameters used in the calculation are σ 0 = 70 MPa , σ sd = −26 MPa , σ si = 0 MPa , A = 1 , h0 L = 0.5 [86]. Figure 4.8 shows the steady state stress as a function of deposition rate for β = 0 (uniform GB diffusivity) and β = 1, 4 (inhomogeneous GB diffusivity). In all three cases, the steady state stress varies from compressive to tensile value as we increase the deposition rate. The tensile value at very high deposition rates is independent of the GB diffusivity, while the compressive stress at lower deposition rates is strongly influenced by it. These trends are qualitatively similar to experimental observations for Ag [26] and Ni [79]. 50 β=0 40 β=1 Steady state Stress (MPa) 30 β=4 20 10 0 -10 -20 -30 0 1 2 3 10 10 10 10 Normalized Deposition rate (dH/dτ) 51 Figure 4.8: Variation of the steady state stress with deposition rate during film growth under varying degrees of GB inhomogeneity represented by different values of β in Equation (4.1). The deposition rate has been normalized as H& = h& (h t ) , where h is the initial film thickness 0 GB 0 and t GB is given by Equation (3.22). Other numerical parameters used in the calculation are σ 0 = 70 MPa , σ s = −26 MPa , A = 1 , h0 L = 0.5 [86]. Table 4.1. GB diffusivities of Ag, Cu and Sn at room temperature. Film Ag Cu Sn δ GB DGB (m 3 s) 9.8 × 10 −28 [87] 3.6 × 10 −30 [88] 3.6 × 10 −22 [89] 4.2 Infinite GB diffusivity In this section, the analytical solution (Equations (3.38) and (3.41)) and the numerical solution (Equation (3.46)) in the limit of infinite GB diffusivity are compared with the experimentally measured compressive stress evolution during deposition and growth interruption of Sn films on substrate [28] Sn is a high mobility material. The room temperature for Sn corresponds to 60% of its melting temperature. Table 4.1 shows that the GB diffusivity of Sn is several orders of magnitude higher than other high mobility materials such as Ag and Cu at room temperature. The experiments by Shin and Chason [28] consist of three cycles of electrodeposition of 2 μm of Sn on Sn/Si substrate at a rate of 2.7 nm / s , each followed by 10 minutes of growth interruption. In the experiments, a 1μm thickness Sn was initially vapor deposited on Si substrate (i.e. h0 = 1 μm ). Electrodeposition on top of this seeding layer ensures that the starting material always had the same initial stress and grain structure. This eliminated the nucleation and island coalescence stage and enabled 52 us to directly model the compressive stress evolution during the third stage of growth [28]. The assumption of infinite GB and infinite surface diffusivity in our model can be checked by comparison of the relevant time scales during Sn film deposition [28]. In our kTh 2 L model, the relevant time scales are τ GBD = for GB diffusion [6-7,27] and ΩE ∗δ GB DGB τ D = h h& for surface deposition. The ratio τ GBD kThLh& = (4.2) τD ΩE ∗δ GB DGB is estimated to be around 0.1 for Sn under the experimental conditions in [28]. This confirms that the GB diffusion time scales are indeed small compared to that associated with deposition. To the best of our knowledge, there is no reported measurement for surface diffusivity (DS ) of Sn at room temperature. Philibert [90] gave an empirical relation for surface diffusivity of fcc metals as: ⎧⎪740 exp(− 15TM T ) cm 2 s T > 0.75TM DS = ⎨ , (4.3) ⎪⎩0.14 exp(− 7TM T ) cm 2 s T < 0.75TM where TM is the melting temperature. Although Sn (TM = 505K) has a tetragonal structure, as a rough estimate we use Equation (4.3) to obtain an approximate estimate for the surface diffusivity of Sn at room temperature to be DS = 1.07 ×10 −10 m 2 s . We assume the characteristic time for surface diffusion as τ S = L2 DS . For the Sn film in [28], the ratio τ S τ D = L2 h& (DS h ) is estimated to be around 5x10-5, which supports our assumption of infinite surface diffusivity. Figure 4.9 shows excellent agreements among Stress-Thickness curves predicted from the analytical solution, numerical analysis and corresponding experimental data from 53 three cycles of electrodeposition and growth interruption of Sn film [28]. The parameters used in the theoretical models are listed below: • Plane strain modulus E ∗ = 58 GPa , • Atomic volume Ω = 1.28 × 10 −29 m 3 , • Deposition time in a cycle TD = 2 μm / 2.7 nm / s = 740 s [28], • Stoppage time in a cycle TS = 600 s [28], • Initial stress σ 0 = 70 MPa , • σ sd = −15 MPa , • σ si = 1.9 MPa , • Grain size L = 2 μm , • Scaled deposition rate H& = 0.003 . Figure 4.9: Comparison of the predicted Stress-Thickness curve from analytical and numerical solutions with corresponding experimental data measured during a three cycle depositon- interruption experiment on Sn film deposited on Sn/Si substrate at room temperature [28,85]. 54 Based on H& = 0.003 and Equation (3.44), t D during deposition, is estimated to be 2.6 × 10 −4 s . The time t I during growth interruption is found to be 4.27 × 10 −3 s . From these values and Equations (3.31) and (3.42), it is estimated that C s Γs = 9.29 × 10 4 s −1 and C g Γg = 5.66 × 10 3 s −1 , with (Cs Γs ) (C g Γg ) = 16 . This can be attributed to different concentrations and attempt frequencies of atoms at the surface and within GBs. Figure 4.10: Comparison of steady state stress versus deposition rate from analytical model and experiments on growth of Sn film on Sn/Si substrate at room temperature [85]. During the cyclic electrodeposition and growth interruption, both theoretical and experimental results indicate that, as soon as the deposition is initiated, the Stress- Thickness drops quickly and approaches a steady-state linear variation with film thickness. Comparison between the analytical predictions and experimental measurements supports the proposed mechanism of atomic diffusion into GBs. During 55 deposition, the film always attains the same steady state stress so that the Stress- Thickness is linearly proportional to the film thickness (Figure 4.9). Figure 4.10 compares the relationship between the steady state stress and deposition rate of Sn film predicted from Equation (3.39) and experimental measurement [28]. The analytical model correctly predicts compressive steady state stresses at slow deposition rates and tensile stresses at higher deposition rates. Therefore, for high mobility materials like Sn, the steady state stress is only a function of deposition rate and independent of the film thickness. Also, the influence of film thickness on transient stress behaviors during deposition and growth interruption is predicted by the analytical model. In particular, the characteristic time required to reach the steady state is predicted to increase with film thickness, in agreement with experiment (Figure 4.9). Chapter 5 5. Crack nucleation in battery electrodes under diffusion induced stresses 5.1 Background In the literature, the insertion/extraction of Li in an electrode has often been modeled as diffusion of interstitial atoms at the continuum level [45-50]. Most of the existing studies have focused on analysis of stresses induced by the diffusion of Li in a host particle, a subclass of problems more broadly referred to as the diffusion induced stresses (DIS) [48-50,77,91-93]. Relatively few theoretical studies have explicitly considered the mechanisms of crack nucleation and propagation. Huggins and Nix [94] considered a bilayer plate structure in which the top layer is subjected to a swelling transformation strain while the bottom layer contains a pre-existing crack. They showed that the swelling in the top layer causes a biaxial tensile stress in the bottom layer and used the Griffith criterion to predict a terminal thickness of the plate below which the pre-existing crack will not propagate. 56 57 Experiments have shown that crack nucleation usually occurs during the first intercalation-deintercalation cycle for the high capacity electrodes. However, the condition for crack nucleation under diffusion induced stresses in an initially crack-free electrode has not been addressed previously. Motivated by the increasing importance of this problem, here we develop a cohesive model of crack nucleation under diffusion induced stresses in battery electrodes under galvanostatic charge and discharge. Compared to the Huggins-Nix model [94], the electrode in our analysis is assumed to be initially crack-free, and the dynamic evolution of DIS and crack nucleation under galvanostatic (constant current) charging/discharging conditions will be explicitly modeled. We adopt the triangular traction-separation law [75] in modeling crack nucleation under DIS. In this model, if the maximum stress in the electrode is less than the cohesive strength of the material, there exists no cohesive zone at all. When the maximum stress exceeds the cohesive strength, the deformation in the electrode starts to localize into an array of cohesive zones. Such localized deformation is thought to be initially reversible, and crack nucleation is assumed to occur only when the maximum surface separation within the cohesive zone reaches a critical value. We will focus on crack nucleation under the most severely loaded state, which will be shown to correspond to a steady state phase after the initial transient has passed but before the maximum stoichiometric solute concentration of the host material is reached. Thin strip and cylindrical geometries are considered for the electrode. The latter geometry is important due to the development of various forms of cylindrical electrodes such as nanorods, nanopillars and nanowires [43-44] and allows us to understand three dimensional (3D) 58 geometry effects on the crack nucleation phenomenon when compared with strip geometry (2D) results. 5.2 Diffusion induced stresses in electrodes The transport of solute in the electrode is modeled as a concentration driven bulk diffusion process (Section 2.3.3), ∂c = D∇ 2 c , (5.1) ∂t where D is the diffusivity and c is the molar concentration of solute and ∇ is the gradient operator. The initial condition is solute free electrode (c = 0) and the boundary condition along the surface of the electrode depends on the mode of battery operation. Under potentiostatic (constant voltage) control c =c s , (5.2) where c s is the known concentration on the surface. Under galvanostatic (constant current) control (Equations (2.28), (2.29), (2.31)), I J n = J ⋅ n = − D∇ n c = − , (5.3) F where I is the known current density on the surface and F = 96486.7 C mol is the Faraday’s constant. The electrode material is considered an isotropic linear elastic solid and the deformation is assumed quasi-static, as atomic diffusion in solids is a slower process compared to elastic deformation. Following an analogy between Diffusion induced stress (DIS) and thermal stresses [48-50,77,91-93], the stress-strain relations can be expressed in index notations as [76], 59 ε ij = 1 E [ ] (1 + ν )σ ij − νσ kk δ ij + Ωc δ ij 3 (5.4) where ε ij are strain components, σ ij are stress components, E is the Young’s modulus, ν is the Poisson’s ratio, Ωc 3 is the swelling transformation strain caused by the insertion of solute atoms into the host. Equation (5.4) can be expressed in terms of stress components, ⎛ Ω(3λ + 2 μ ) ⎞ σ ij = 2 με ij + ⎜ λε kk − c ⎟δ ij (5.5) ⎝ 3 ⎠ where μ = E 2(1 + ν ) and λ = Eν [(1 + ν )(1 − 2ν )] are Lame’s constant. Strain-displacement relations for small strain theory are given as [76] ε ij = 1 2 ( u i, j + u j , i ) (5.6) and the static equilibrium equations in the absence of body forces, are [76] σ ij , j = 0 . (5.7) In these expressions, index after comma denoting partial differentiation with respect to the appropriate coordinate, i.e., u i , j = ∂u i ∂x j . Thus given a molar concentration field based on the solution of Equation (5.1) under appropriate boundary conditions (Equations (5.2-5.3)), stress in an elastic body can be calculated using Equations (5.4-5.7). 5.2.1 Strip electrode Figure 5.1 shows a two-dimensional (2D) strip electrode with width 2h , subject to insertion and extraction of an interstitial species such as Li along the thickness (y) direction. The diffusion equation (5.1) for strip geometry simplifies as, 60 ∂c ∂ 2c =D 2 . (5.8) ∂t ∂y a y y 2h 2a 2h x x p p a (a) (b) Figure 5.1: Schematic illustration of crack nucleation in a strip electrode of width 2h during galvanostatic solute (a) intercalation and (b) extraction, modeled as diffusion along the thickness direction (y-axis). The crack nuclei are uniformly spaced with period p and modeled as cohesive zones obeying the triangular traction-separation law (Equation (5.30)). Solving Equations (5.5.-5.7), leads to following axial stress in the electrode [95], EΩc ( y , t ) h h EΩ EΩy σ D ( y, t ) = − + ∫ c( y ′, t )dy ′ + ∫ y ′c( y ′, t )dy ′ . (5.9) 3(1 − ν ) 6(1 − ν )h −h 2(1 − ν )h 3 −h We consider the variations of solute concentration and the corresponding DIS during charging and discharging. Under galvanostatic boundary conditions (Equation (5.3)) and initially solute free electrode, as shown in Figure 5.1, ∂c ∂c I −D =D =− , (5.10) ∂y h ∂y −h F the solute concentration during insertion can be found as [78] c( y, t ) Dt 3 y 2 − h 2 2 ∞ (− 1)n cos⎛ nπ y ⎞ exp⎧− Dn 2π 2t ⎫ Ih FD h 2 = + 6h 2 − 2 π ∑ n2 ⎜ ⎟ ⎝ h ⎠ ⎨ h2 ⎬, (5.11) n =1 ⎩ ⎭ and the associated DIS is 61 σ D ( y, t ) = EΩ Ih ⎡ h 2 − 3 y 2 2 ∞ (− 1)n cos⎛ nπ y ⎞ exp ⎧− Dn 2π 2t ⎫⎤ . (5.12) ⎢ 3(1 − ν ) FD ⎣ 6h 2 + 2 π ∑ n2 ⎜ ⎟ ⎝ h ⎠ ⎨ h2 ⎬⎥ n =1 ⎩ ⎭⎦ At the end of charging, the stress approaches a steady state while the solute concentration rises steadily with time. This situation persists until the saturation limit of material is reached. The steady state solution then acts as the initial condition for the extraction process, Ih ⎡ 3 y 2 − h 2 ⎤ c( y ,0) = c1 + , (5.13) FD ⎢⎣ 6h 2 ⎥⎦ where Itc c1 = , (5.14) Fh t c denoting the charging time. During extraction, the solute concentration evolves as c( y , t ) − c1 Dt 3 y 2 − h 2 4 ∞ (− 1)n cos⎛ nπ y ⎞ exp⎧− Dn 2π 2t ⎫ Ih FD =− 2 − h 6h 2 + 2 π ∑ n2 ⎜ ⎟ ⎝ h ⎠ ⎨ h2 ⎬ (5.15) n =1 ⎩ ⎭ with the associated DIS σ D ( y, t ) = EΩ Ih ⎡ 3 y 2 − h 2 4 ∞ (− 1)n cos⎛ nπ y ⎞ exp⎧⎪− Dn 2π 2 t ⎫⎪⎤ ⎢ 3(1 − ν ) FD ⎣⎢ 6h 2 − 2 ∑ ⎜ ⎟ ⎨ ⎬⎥ . (5.16) π n =1 n2 ⎝ h ⎠ ⎪⎩ h2 ⎪⎭⎦⎥ Figure 5.2 plots the variations of solute concentration and the associated DIS during the first charging and discharging cycle. Due to the symmetry of the problem, the results are plotted only over half of the strip width. During insertion, the solute concentration continuously rises with time (Figure 5.2a) while the stress approaches a steady state with tension near the center and compression near the free surface of the electrode (Figure 5.2b). The peak tensile stress occurs at the centre with magnitude equal to (Figure 5.2b) 62 I σ peak = EΩIh 18(1 − ν )FD (5.17) when reaching the steady state. During extraction, the surface current is reversed, and the solute concentration continuously decreases with time (Figure 5.2c) while the stress approaches a steady state with compression near the center and tension at the surface of the electrode. In this case, the peak tensile stress at the surface reaches (Figure 5.2d) σ peak E = EΩIh 9(1 − ν )FD (5.18) at steady state. (a) (b) (c) (d) Figure 5.2: Snapshot profiles of solute concentration and diffusion induced stress in a strip electrode. (a) Concentration during insertion, (b) DIS during insertion, (c) concentration during 63 extraction and (d) DIS during extraction. The concentration is normalized as cˆ = cFD (Ih ) during insertion and cˆ = (c − c1 )FD (Ih ) during extraction (Equations (5.11) and (5.15)), while DIS is normalized as σˆ = 3(1 − ν )FDσ (EΩIh ) (Equations (5.12) and (5.16)) [96]. To appreciate the level of diffusion induced stress in high capacity electrodes, we consider silicon nanowire electrodes with an average diameter of 89nm. Recent experiments have shown such electrodes can be charged to the near theoretical capacity of 4227 mAh g–1 at a charge and discharge rate of 20 hours per half cycle [44], corresponding to a surface current density of I = 0.011 A m 2 . With concentration change in LixSi from x = 0 to 5.4, reaching a volume change as large as 59 %, the partial molar volume of Li can be estimated as [48] 4 Ω= = 2 × 10 −5 m 3 mol (5.19) 4.4cmax where cmax = 2.0152 × 104 mol/m3 is the stoichiometric maximum concentration of Li [34]. Other material properties are listed in Table 5.1. If the lower limit of Young’s modulus of fully lithiated Si is used, the peak tensile stress is estimated to be 0.1 GPa during Li insertion and 0.2 GPa during Li extraction. For faster charge and discharge rates, the nanowire electrodes begin to show irreversible capacity losses even during the first cycle [44]. For example, at the charge and discharge rate of 5 hours per half cycle, there is an irreversible capacity loss at the first cycle, but the capacity is then stabilized at 3500 mAh g–1 for the subsequent 20 cycles [44]. In this case, the surface current density is I = 0.036 A m 2 , corresponding to a peak stress of 0.35 GPa during Li insertion and 0.7 GPa during Li extraction. Note that the above solutions to DIS neglect a number of nonlinear coupling effects and may be oversimplified in a number of ways. The reader is referred to Christensen and 64 Newman [46-47] for some detailed discussions. By introducing a coupling between internal stresses and activation energy for diffusion, Haftbaradaran et al. [97] also discovered a class of nonconventional solutions to DIS with a surface choking instability once the product between electrode dimension and charging rate exceeds a critical value. Table 5.1. Material properties of Si and operating parameters. Parameter Symbol Value Source (dimension) Young’s modulus (lithiated Si) E (GPa) 30-80 Poisson’s ratio ν 0.22 [1] 2 Diffusion coefficient D (m /s) 2 × 10 −18 [29] 3 Stoichiometric maximum concentration c max (mol/m ) 2.0152 × 10 4 [34] Fracture energy Γ (J/ m2) 2 5.2.2 Cylindrical electrode † Figure 5.3 shows a cylindrical electrode with diameter 2 rc subject to insertion and extraction of an interstitial species such as Li. Under the assumption of diffusion along radial (r) direction of the electrode, Equation (5.1) reduces to [78], ∂c 1 ∂ ⎛ ∂c ⎞ =D ⎜r ⎟, (5.20) ∂t r ∂r ⎝ ∂r ⎠ and the corresponding axial stress in the electrode based on Equations (5.4-5.7) is [95], † Results for cylindrical electrodes are submitted for publication and currently under review 65 EΩ ⎡ 2 rc ⎤ σ D (r , t ) = ⎢ 2 r ′c(r ′, t ) dr ′ − c(r , t )⎥ . ∫ (5.21) 3(1 − ν ) ⎢ rc 0 ⎥⎦ ⎣ z z r r rc 2a a p p rc (a) (b) Figure 5.3: Schematic illustration of crack nucleation in a cylindrical electrode of diameter 2rc during galvanostatic (a) intercalation and (b) extraction, modeled as diffusion along the radial direction (r-axis). The axisymmetric crack nuclei are uniformly spaced with period p and modeled as localized cohesive zones obeying the triangular traction-separation law (Equation (5.30)). Similar to strip electrode, we consider the variations of solute concentration and the corresponding DIS during charging and discharging of a cylindrical electrode. The initial solute concentration in the electrode is assumed to be zero. Under galvanostatic boundary conditions (Equation (5.3)) and initially solute free electrode, as shown in Figure 5.3, ∂c I ∂c −D =− ; −D =0, (5.22) ∂r rc F ∂r 0 the solute concentration during insertion can be found as [78] 66 c(r , t ) 2 Dt r 2 1 ∞ J 0 (α n r rc ) ⎧ Dtα n2 ⎫ Irc FD rc2 = + 2 − −2 2rc 4 ∑ n =1 α n J 0 (α n ) 2 exp⎨− 2 ⎬ , ⎩ rc ⎭ (5.23) and the associated DIS is EΩ Irc ⎡ 1 r 2 ∞ J 0 (α n r rc ) ⎧ Dtα n2 ⎫⎤ σ D ( y, t ) = ⎢ − 2 +2 3(1 −ν ) FD ⎣⎢ 4 2rc ∑ n =1 α n J 0 (α n ) 2 exp ⎨ − 2 ⎬⎥ , (5.24) ⎩ rc ⎭⎦⎥ where J 0 (r ) is the Bessel function of the first kind and α n are the roots of J 1 (α ) . At the end of charging, the stress approaches a steady state while the solute concentration rises steadily with time. This situation persists until the saturation limit of material is reached. The steady state solution then acts as the initial condition for the extraction process, Irc ⎡ r 2 1 ⎤ c ( r ,0) = c1 + ⎢ − ⎥, (5.25) FD ⎣ 2rc2 4 ⎦ where 2 It c c1 = , (5.26) Frc t c denoting the charging time. During extraction, the solute concentration evolves as c(r , t ) − c1 2 Dt r 2 1 ∞ J 0 (α n r rc ) ⎧ Dtα n2 ⎫ Irc FD =− 2 − 2 + +4 rc 2rc 4 2 ∑ n=1 α n J 0 (α n ) exp ⎨− 2 ⎬ , ⎩ rc ⎭ (5.27) with the associated DIS EΩ Irc ⎡ r 2 1 ∞ J 0 (α n r rc ) ⎧ Dtα n2 ⎫⎤ σ D (r , t ) = ⎢ − 3(1 −ν ) FD ⎢⎣ 2rc2 4 − 4 ∑ n =1 α n J 0 (α n ) 2 exp ⎨− 2 ⎬ ⎥ . (5.28) ⎩ rc ⎭⎥⎦ Figure 5.4 plots the variations of solute concentration and the associated DIS during the first charging and discharging cycle. During insertion, the solute concentration continuously rises with time (Figure 5.4a) while the stress approaches a steady state with tension near the center and compression near the free surface of the electrode (Figure 67 5.4b). The peak tensile stress occurs at the centre when reaching the steady state, with magnitude equal to (Figure 5.4b) σ peak = EΩIrc 12(1 − ν )FD . (5.29) During extraction, the surface current is reversed, and the solute concentration continuously decreases with time (Figure 5.4c) while the stress approaches a steady state with compression near the center and tension at the surface of the electrode. The peak tensile stress occurs at the surface with the same magnitude as Equation (5.29) when reaching the steady state (Figure 5.4d). (a) (b) (c) (d) Figure 5.4: Snapshot profiles of solute concentration and diffusion induced stress in a cylindrical electrode. (a) concentration during insertion, (b) DIS during insertion, (c) concentration during extraction and (d) DIS during extraction. The concentration is normalized as cˆ = cFD (Irc ) 68 during insertion and cˆ = (c − c1 )FD (Irc ) during extraction (Equations (5.23) and (5.27)), while DIS is normalized as σˆ D = 3(1 − ν )FDσ D (EΩIrc ) (Equations (5.24) and (5.28)). Analogous to strip electrode, we can estimate the peak stress for the cylindrical electrode based on the Equation (5.29) with reference to the recent experiments on silicon nanowire electrodes [44]. For a charge and discharge rate of 20 hours per half cycle [44], corresponding to a surface current density of I = 0.011 A m 2 , peak tensile stress during insertion and extraction is estimated to be 0.16 GPa. For a faster charge and discharge 2 rates of 5 hours per half cycle, the surface current density is I = 0.036 A m , corresponding to a peak stress of 0.53 GPa during Li insertion and extraction. The estimated peak tensile stresses in strip and cylindrical electrodes are close to the theoretical strength of pure silicon and have most likely exceeded the cohesive strength of Lithiated silicon. Why can the Si nanowires sustain such extreme mechanical stresses without fracture? In the following section, we consider crack nucleation under diffusion induced stresses in an initially crack-free electrode. 5.3 Cohesive model of crack nucleation in the electrode Figure 5.1 and 5.3 show an emergent array of cohesive zones in the initially crack-free electrode. The cohesive zones are uniformly spaced at a period of p near the center of the electrode during solute insertion and at the edge of the electrode during solute extraction. The cohesive zone is assumed to obey the triangular traction-separation ( σ − δ ) law [75], ⎧σ (1 − δ δ c ) δ ≤ δc σ =⎨ c , (5.30) ⎩0 δ > δc 69 where σ c is the cohesive strength of the electrode material and δ c is the maximum range of cohesive interaction. The fracture energy of the material, Γ = σ c δ c 2 , is assumed to be a material constant typically on the order of surface energy in the absence of significant plastic deformation. The emergent cohesive zones are modelled as continuous distributions of dislocations [71]. Similar to the Dugdale model [52], the cohesive zone size is determined based on the condition that there exists no singularity at the tip of the cohesive zone and crack nucleation is assumed to occur when the maximum surface separation reaches δ c = 2Γ / σ c . Finally, the spacing p between the cohesive zones will be determined by the condition that the stress everywhere in the electrode must not exceed the cohesive strength. 5.3.1 Strip electrode During solute insertion, the cohesive zones would develop at the centre of the electrode as soon as the axial stress exceeds the cohesive strength. Within the cohesive zone, the traction and the surface separation follows Equation (5.30), i.e. a ⎛ σ a ⎞ σ D ( y, t ) + ∫ P ( y,η ) B(η , t )dη = σ c ⎜1 − c ∫ B(η , t )dη ⎟, − a ≤ y ≤ a (5.31) ⎜ 2Γ ⎟ −a ⎝ y ⎠ where the first term σ D ( y , t ) is the diffusion induced stress and the second term is the stress associated with localized deformation within the cohesive zones modeled as continuous distributions of infinitesimal dislocations with density B (η , t ) satisfying a ∫ B(η , t )dη = 0 . (5.32) −a 70 The kernel function P ( y,η ) in Equation (5.31) corresponds to the axial stress at location (0, y ) induced by an array of edge dislocations of unit Burgers vector in the x − direction, located at position η along the y-axis and periodically distributed along the x-axis with a period equal to p . The expression for P ( y,η ) is given in Appendix B. The requirement of no singularity at the tip of cohesive zone and the crack nucleation condition based on the maximum surface separation leads to lim B(η , t ) a − η = 0 , (5.33) η →a a ∫ B(η , t )dη = 2Γ / σ 0 c . (5.34) During solute extraction, the tensile stress region is shifted to the surface of the electrode while the centre of the electrode is under compression. Therefore, in this case the emergent cohesive zones are placed periodically along the edge of the electrode with governing equation h ⎛ σc y ⎞ σ D ( y, t ) + ∫ [P( y,η ) − P( y,−η )]B(η , t )dη = σ c ⎜⎜1 + 2Γ ∫ B(η , t )dη ⎟⎟, h−a≤ y ≤ h (5.35) h−a ⎝ h−a ⎠ In this case, the size of the cohesive zones is determined based on lim B (η , t ) η − h + a = 0 , (5.36) η →h − a and the corresponding crack nucleation condition is h−a ∫ B(η , t )dη = 2Γ / σ h c . (5.37) Normalizing all stress variables by σ c and all length variables by h in Equations (5.12), (5.16) and (5.31)-(5.37), we can identify a characteristic length scale as 71 ⎧ Γ(1 − ν )F 2 D 2 ⎫ 13 l ft = ⎨ 2 2 ⎬ . (5.38) ⎩ E (1 + ν )Ω I ⎭ 5.3.1.1 Localization spacing p   In the present cohesive model, the deformation in the electrode is spontaneously localized into a periodic array of cohesive zones along the length of the strip if the peak DIS in the electrode exceeds the cohesive strength of material. The strain localization leads to stress relaxation in the vicinity of each cohesive zone. Therefore, the post-localization stress distribution depends on the localization spacing p . If the localization spacing is too large, the stress between two adjacent cohesive zones will not be brought down to the level of cohesive strength. Oppositely, if the localization spacing is too small, the region between two adjacent cohesive zones will be over-shielded such that the maximum stress falls below the cohesive strength. (a) (b) 72 (c) Figure 5.5: Effect of localization spacing p on the distribution of axial stress σ x in the region between two adjacent cohesive zones of length a h = 0.7 during insertion in a strip electrode (Figure 5.1a). (a) At the spontaneous localization spacing p h = 2.36 , the maximum axial stress between two adjacent localization zones is equal to the cohesive strength σ c . (b) If the localization spacing is taken to be p h = 5 , the maximum axial stress between two adjacent localization zones is seen to be greater than the cohesive strength σ c . (c) If the localization spacing is taken to be p h = 1 , the maximum axial stress between two adjacent localization zones is seen to be lower than the cohesive strength σ c [96]. This is illustrated in Figure 5.5 with contour plots of axial stress in the electrode for the cohesive zone length a h = 0.7 during insertion. Due to the symmetry of the problem, the results are plotted only over a quarter of the region between two adjacent cohesive zones. When the cohesive zones are too widely spaced, as shown in Figure 5.5b for p / h = 5 , the maximum stress in the intermediate region between two adjacent cohesive zones is greater than the cohesive strength. When the cohesive zones are too narrowly spaced, as shown in Figure 5.5c for p / h = 1 , the stress in the intermediate region is over-relaxed to below the cohesive strength. We assume that the localization process would naturally select the cohesive zone spacing such that the maximum stress in the region between two adjacent cohesive zones is exactly equal to the cohesive strength σ c . In the case shown in 73 Figure 5.5, this critical spacing is found to be p / h = 2.36 (method to determine p / h will be discussed shortly), as can be seen in Figure 5.5a. For the formation of centre cohesive zones during solute insertion (Figure 5.1a), the maximum axial stress in the region between two adjacent cohesive zones occurs along the axis of the strip right in the middle of the two localization zones. Hence, in this case the cohesive zone spacing p is determined by solving a σ x ( p 2 ,0 ) = σ D (0, t ) + ∫ H ( p 2 ,0,η )B(η , t )dη = σ c (5.39) −a together with Equations (5.31-5.34). The kernel function H ( x, y,η ) corresponds to the axial stress at (x, y ) induced by an array of edge dislocations of unit Burgers vector in the x − direction located at position η along the y-axis and periodically distributed along the x-axis with a period equal to p . The expression for kernel function H ( x, y,η ) is given in Appendix B. For a given cohesive strength σ c , we solve Equations (5.31-5.34), followed by checking Equation (5.39) and employing the method of bisection to determine p . For the formation of edge cohesive zones during solute extraction (Figure 5.1b), the maximum axial stress in the region between two adjacent cohesive zones occurs along the free surface at the midpoint between the two localization zones. The cohesive zone spacing p is then determined by solving h σ x ( p 2 , h ) = σ D (h, t ) + ∫ [H ( p / 2, h,η ) − H ( p / 2, h,−η )]B(η , t )dη = σ h−a c (5.40) together with Equations (5.35-5.37). Figures 5.6a and 5.6b plot the cohesive zone spacing as a function of the cohesive strength during solute insertion and extraction, respectively. 74 (a) (b) Figure 5.6: Localization spacing p as a function of the normalized cohesive strength for the formation of (a) centre cohesive zones during solute insertion and (b) edge cohesive zones during solute extraction, in a strip electrode. In both (a) and (b), the cohesive strength σ c is normalized by the peak stresses σ peak I = EΩIh [18(1 − ν )FD ] during solute insertion and σ peak E = EΩIh [9(1 − ν )FD ] during solute extraction [96]. Once the localization spacing p is determined, the critical conditions for crack nucleation are obtained from Equations (5.31-5.34) for center cracks during solute insertion and from Equations (5.35-5.37) for edge cracks during solute extraction by a numerical scheme detailed in Appendix C. When the peak DIS exceeds the cohesive strength, corresponding to the normalized cohesive strength smaller than unity, strain localization occurs spontaneously as a periodic array of cohesive zones with spacing given by Figure 5.6. When the peak DIS is smaller than the cohesive strength, corresponding to the normalized cohesive strength exceeding unity, strain localization does not occur spontaneously. However, in this case there may exist metastable solutions with isolated localization zones. This type of failure may occur at locations with pre-existing defects/weaknesses. The governing equations and numerical algorithm for such metastable, isolated localizations are similar to those 75 for spontaneous, periodic localizations, except the Green’s function kernel is replaced with K ( y ,η ) which corresponds to the axial stress at location (0, y ) induced by a single edge dislocation at (0,η ) with a unit burgers vector in the x − direction [98]. The expression for K ( y ,η ) is given in Appendix D. Figure 5.8 indicates that the critical strip width and cohesive zone size transition smoothly between multiple localization to isolated localization regimes. 5.3.2 Cylindrical electrode The emergent cohesive zones (Figure 5.3) are axi-symmetric and modeled using prismatic dislocation loops. During solute insertion, the cohesive zones would develop at the centre of the electrode (Figure 5.3a) as soon as the stress exceeds the cohesive strength. Within the cohesive zone, the traction and the surface separation obey Equation (5.30), i.e. a ⎛ σc a ⎞ σ D (r , t ) + ∫ P(r , R) B( R, t )dR = σ c ⎜1 − B ( R, t )dR ⎟, ∫ 0≤ R≤a (5.41) ⎜ 2Γ r ⎟ 0 ⎝ ⎠ where the first term σ D (r , t ) is the diffusion induced stress and the second term is the stress associated with localized deformation within the cohesive zones modeled as continuous distributions of prismatic dislocation loops of radius R with density B( R, t ) . The kernel function P(r , R ) in Equation (5.41) corresponds to the axial stress at location (r , 0) induced by an infinite array of coaxial prismatic dislocation loops of unit Burgers vector in the z-direction and radius R, located at positions z = np; n = −∞, ..,∞ along the axis. In practice, since the stress field associated with a prismatic loop in a cylinder 76 decays cubically [99] with distance along the axis of the cylinder, we determine P (r , R ) from 5 prismatic loops with spacing p along the cylindrical axis. Similar approximation has been adopted previously in studying interaction among periodic array of cracks in a layer [100]. The expression for P (r , R ) is given in Appendix E. The non-singulairty condition at the cohesive zone tip and the crack nucleation condition based on the maximum surface separation lead to, lim B(r , t ) a − r = 0 , (5.42) r →a a ∫ B( R, t )dR = 2Γ / σ 0 c . (5.43) During solute extraction, the tensile stress region is shifted to the surface of the electrode while the centre of the electrode is under compression. Therefore, the emergent cohesive zones are placed periodically along the edge of the electrode with governing equation rc ⎛ σ r ⎞ σ D (r , t ) + P(r , R )B(R, t )dη = σ c ⎜1 + c B(R, t )dR ⎟, ∫ ∫ rc − a ≤ R ≤ rc (5.44) ⎜ 2Γ ⎟ rc − a ⎝ rc − a ⎠ In this case, the size of the cohesive zones is determined based on lim B(r , t ) r − rc + a = 0 , (5.45) r →rc −a and the corresponding crack nucleation condition is rc − a ∫ B(R, t )dR = 2Γ / σ rc c . (5.46) 77 5.3.2.1 Localization spacing p Consistent with our assumption for the cohesive zone spacing in strip electrode, localization process would naturally select the spacing in cylindrical electrode such that the maximum axial stress in the region between two adjacent cohesive zones is exactly equal to the cohesive strength σ c . For the formation of centre cohesive zones during solute insertion (Figure 5.3a), the maximum axial stress in the region between two adjacent cohesive zones occurs along the axis of the cylinder right in the middle of the two localization zones. Hence, the cohesive zone spacing p is determined by solving a σ z (0, p 2) = σ D (0, t ) + ∫ H (0, p 2 , R )B(R, t )dR = σ c (5.47) 0 together with Equations (5.41-5.43). The kernel function H (r , z , R ) corresponds to the axial stress at (r, z ) induced by an array of five coaxial circular prismatic dislocation loops of unit Burgers vector in the z − direction and radius R , located at z = np; n = −2, .., 2 along the axis. The expression for kernel function H ( r , z , R) is given in Appendix E. For a given cohesive strength σ c , we solve Equations (5.41-5.43), followed by checking Equation (5.47) and employing the method of bisection to determine p . For the formation of edge cohesive zones during solute extraction (Figure 5.3b), the maximum axial stress in the region between two adjacent cohesive zones occurs along the free surface at the midpoint between the two localization zones. The cohesive zone spacing p is then determined by solving rc σ z ( rc , p 2 ) = σ D (rc , t ) + ∫ H (r , p / 2, R )B(R, t )dR = σ c c (5.48) rc − a 78 together with Equations (5.44-5.46). Figures 5.7a and 5.7b plot the cohesive zone spacing as a function of the cohesive strength during solute insertion and extraction, respectively. Once the localization spacing p is determined, the critical conditions for crack nucleation are obtained from Equations (5.41-5.43) for center cracks during solute insertion and from Equations (5.44-5.46) for edge cracks during solute extraction by a numerical scheme detailed in Appendix F. (a) (b) Figure 5.7: Localization spacing p as a function of the normalized cohesive strength for the formation of (a) centre cohesive zones during solute insertion and (b) edge cohesive zones during solute extraction, in a cylindrical electrode. In both (a) and (b), the cohesive strength σ c is normalized by the peak DIS σ peak = EΩIrc [12(1 − ν )FD ] . Metastable solution in the form of nucleation of isolated localization zone also exists for cylindrical electrodes. The governing equations and numerical algorithm for such metastable, isolated localizations are similar to those for spontaneous localizations, except the Green’s function kernel is replaced with S (r , R ) which corresponds to the axial stress at location (r ,0) induced by a single coaxial circular prismatic dislocation loop of a unit burgers vector in the z − direction and radius R with center located at z = 0 . The 79 expression for S (r , R ) is given in Appendix E. Figure 5.9 indicates that the critical radius and cohesive zone size undergo smooth transition between multiple localization to isolated localization regimes. 5.4 Critical conditions for crack nucleation The results are shown as blue solid lines relating the normalized electrode dimension and the cohesive strength scaled by the corresponding peak stress in the electrode in Figures 5.8a-b, 5.9a-b for strip and cylinder geometry respectively. A comparison of the results between strip and cylindrical electrodes reveals that the geometry does not alter the relationships between cohesive strength, characteristic electrode size and cohesive zone size qualitatively. For given values of electrode dimension and cohesive strength, crack nucleation is predicted to occur along and above the blue lines in the sense that there exist a solution with maximum surface separation within the localization zones exceeding the cohesive interaction range δ c . If the normalized values of electrode dimension and cohesive strength are below these lines, crack nucleation is predicted not to occur, not because of the absence of strain localization but because the maximum surface separation in the cohesive zone cannot reach δ c . In this case, the localized deformation within the cohesive zone is fully reversible and material would recover as soon as the diffusion induced stress in the electrode is reduced. For given values of the normalized cohesive strength and electrode dimension, the corresponding sizes of cohesive zone at crack nucleation are shown as green dashed lines. 80 (a) (b) (c) (d) Figure 5.8: The critical conditions for crack nucleation expressed as relationships between the normalized half-width of electrode, the normalized cohesive strength and the normalized critical size of cohesive zone at nucleation. The blue lines plot the critical dimension of electrode while the dashed green lines plot the critical size of cohesive zone at crack nucleation as functions of the normalized cohesive strength. (a, c) plot the critical conditions for nucleation of center cohesive zones of length 2a under solute insertion while (b, d) plot those of symmetric edge cohesive zones of length a under solute extraction, in a strip electrode of width 2h . In (a) and I (b), the cohesive strength is normalized by the peak stresses σ peak = EΩIh [18(1 − ν )FD ] during E solute insertion and σ peak = EΩIh [9(1 − ν )FD ] during solute extraction. In (c) and (d), the cohesive strength is normalized by the size-independent reference stress σ ref = EΩIl ft [18(1 − ν )FD ] where l ft is the characteristic length scale defined in Equation (5.38) [96]. 81 (a) (b) (c) (d) Figure 5.9: The critical conditions for crack nucleation expressed as relationships between the normalized radius of electrode, the normalized cohesive strength and the normalized critical size of cohesive zone at nucleation. The blue lines plot the critical dimension of electrode while the dashed green lines plot the critical size of cohesive zone at crack nucleation as functions of the normalized cohesive strength. (a, c) plot the critical conditions for nucleation of center cohesive zones of radius a under solute insertion while (b, d) plot those of edge annular cohesive zones of size a under solute extraction, in a cylindrical electrode of radius rc . In (a), (b), the cohesive strength is normalized by the peak DIS σ peak = EΩIrc [12(1 − ν )FD ] . In (c) and (d), the cohesive strength is normalized by the size-independent reference stress σ ref = EΩIl ft [12(1 − ν )FD ] where l ft is the characteristic length scale defined in Equation (5.38). Most interestingly, Figures 5.8 and 5.9 shows that, during both solute insertion and extraction, there exists a critical electrode dimension below which crack nucleation becomes impossible irrespective of the cohesive strength of the material. This critical 82 dimension is found to be hcrI = 7.3l ft , rcrI = 7.3l ft during solute insertion and hcrE = 6.5l ft , rcrE = 8.2l ft during solute extraction; see Figures 5.8a-b and 5.9a-b. Therefore in strip geometry, crack nucleation is more likely to occur at the surface of electrode during solute extraction, as opposed to nucleation at the center of electrode during solute insertion. In cylindrical geometry, we find that crack nucleation is more likely to occur at the center of electrode during solute insertion, as opposed to nucleation at the edge of electrode during solute extraction. In reality, stress concentration induced by surface roughness may further facilitate crack nucleation at the surface of electrodes during solute extraction. Combining these results, a critical dimension for flaw tolerant electrodes is identified as 13 ⎧⎪ Γ(1 − ν )(FD )2 ⎫⎪ L ft = 13l ft = 13⎨ ⎬ . (5.49) ⎪⎩ E (1 + ν )(ΩI )2 ⎪⎭ The significance of this equation is that it predicts an initially crack-free electrode would remain crack free below the critical dimension. Once the electrode dimension exceeds this critical dimension, nucleation of cracks during solute transport would become possible. Our results significantly extends and generalizes the previous analysis by Huggins and Nix [94] who considered a bilayer plate structure in which the top layer is subjected to a swelling transformation strain eT while the bottom layer contains a pre-existing crack. They showed that the swelling in the top layer causes a biaxial tensile stress in the bottom layer and used the Griffith criterion to predicted a critical thickness 23 ⎛ 3K (1 −ν ) ⎞ 2 H h−n = ⎜⎜ IC ⎟⎟ (5.50) π ⎝ EeT ⎠ 83 of the plate below which the pre-existing crack will not propagate. Here K IC is the fracture toughness of the plate. Using the relation K IC = EΓ (1 − ν 2 ) , we can rewrite the Huggins-Nix critical length as 207 Γ(1 −ν ) H h −n = (5.51) π E (1 + ν )eT 2 Comparing Equation (5.51) and our result in Equation (5.49), we can see some qualitative similarities. Both models predict that the critical electrode dimension for fracture resistant electrode should scale with the ratio between the fracture energy and Young’s modulus of the material. However, the scaling is linear in the Huggins-Nix model while it is nonlinear with a power index of 1/3 in our model. Figures 5.8a-b and 5.9a-b show that the critical dimension for flaw tolerant electrodes corresponds to the minimum dimension of the electrode required to prevent crack nucleation irrespective of the cohesive strength of the material. This is a “fail safe” or “design for robustness” concept. If the cohesive strength of the material is known, the critical dimension of the electrode for crack nucleation may be higher than that predicted by Equation (5.49). The plots in Figures 5.8a-b and 5.9a-b are based on cohesive strength normalized by the peak tensile stresses during solute insertion and extraction. Since these peak stresses also depend on the electrode size, as shown in Equations (5.17-5.18), (5.29), it is difficult to see from Figures 5.8a-b and 5.9a-b the actual critical dimension for crack nucleation as a function of the cohesive strength. In order to decouple the effect of electrode size and cohesive strength, we introduce size-independent reference stress σ ref = {EΩIl ft 18(1 − ν )FD , EΩIl ft 12(1 − ν )FD} for the strip and cylinder geometry respectively. The results based on a new normalization of the cohesive strength with 84 respect to σ ref are shown in Figures 5.8c and 5.8d for strip geometry and Figures 5.9c and 5.9d for cylindrical geometry. Figures 5.8c-d and 5.9c-d indicate that the critical dimension for crack nucleation increases almost linearly with cohesive strength at large values of σ c / σ ref . It is interesting that the critical dimension also increases at small values of σ c / σ ref . This is because, under the assumption of constant fracture energy Γ = σ cδ c 2 , the cohesive interaction range δ c increases as σ c decreases. The “fail safe” electrode size defined in Equation (5.49) corresponds a specific combination of δ c and σ c = σ c∗ that is most susceptible to crack nucleation. When σ c > σ c∗ , the electrode becomes more resistant to cracking due to higher strength of the material. In contrast, when σ c < σ c∗ , the increasing resistance to cracking is due to higher interaction range δ c , which makes it difficult to separate cohesive surfaces. In this sense, crack nucleation can be referred to as “strength-controlled” in the range σ c > σ c∗ and “separation-controlled” in the range σ c < σ c∗ . Take the material parameters listed in Table 5.1. For silicon nanowires electrodes to be charged to near theoretical capacity at a charge-discharge rate of 20h per half cycle without fracture, the critical dimension for flaw tolerant electrodes, L ft , is estimated to be 413 nm. The nanowire electrodes adopted in experiments by Chan et al. [44] indeed fall in this range. At faster charge-discharge rates, this critical dimension scales with increasing surface current density according to I −2 / 3 . For the nanowires of 89nm in diameter adopted in experiments [44], we estimate that crack nucleation would occur at the surface of the nanowire at the charging-discharging rate of 2h per half-cycle. Indeed, 85 the voltage profiles at different power rates show that significant capacity loss begins at the charging/discharging rate between 5h and 10h per half cycle [44]. While the reason for such capacity loss is not completely clear, the observation would be more or less consistent with our analysis if the formation of surface cracks during discharging is assumed to be a significant cause. Earlier experiments have also shown that decrepitation was suppressed in amorphous Si films about 100 nm in thickness [41] and in amorphous Ge films around 60-250 nm in thickness [42], leading to superior charge/discharge cycling performance compared to the bulk materials. The broad agreement between experiments and our analysis for the critical electrode dimension for crack nucleation points to the fact that crack nucleation can indeed be suppressed in nanostructured electrodes and design of fracture resistant electrodes can greatly improve the cycling performance of Li batteries. Chapter 6 Dislocation shielding of a cohesive crack In the classical problem of dislocation shielding of a singular crack, the effective stress intensity factor K eff at the crack tip is expressed as K eff = K ∞ + K ⊥ , (6.1) where K ∞ and K ⊥ are stress intensity factors induced by the externally imposed loading and dislocations, respectively. Since crack growth occurs when K eff reaches a critical value K c , the crack is said to be shielded by dislocations when K ⊥ < 0 and anti-shielded ⊥ when K > 0 . Lin and Thomson [66] and Weertman [67] derived analytical solutions of K ⊥ for a singular crack to quantify the shielding effects of dislocations. However, a recent trend is to use cohesive crack models to study crack initiation and propagation using discrete dislocation plasticity [58-65, 101-102]. In spite of an extensive and growing literature, fundamental questions remain unanswered: in the presence of dislocations, does a cohesive crack behave differently from a mathematically sharp singular crack? If different, how does it impact the use and interpretation of cohesive models in modeling crack initiation and propagation? 86 87 The present chapter is an attempt aimed to clarify the above questions. We focus on a Dugdale cohesive crack interacting with dislocations under mode I dominant conditions and derive closed form analytical solutions. The practical significance of the theoretical study for fracture in metals is demonstrated by DD simulations of a large number ( > 10 3 ) of edge dislocations interacting with a cohesive crack described by the trapezoidal traction-separation law of Tvergaard-Hutchinson [54]. 6.1 Dugdale cohesive crack 6.1.1 Problem Formulation y K I∞ (b cos φ j , b sin φ j ) σ −σc σc Rj Γ θj O′ O x O δc δ L (a) (b) Figure 6.1: A Dugdale cohesive crack interacting with a number of dislocations. (a) The crack lies in an isotropic elastic medium and is subjected to a remote mode-I loading K I∞ . The crack tip ( ) cohesive zone size is L . N edge dislocations with Burgers vectors b cos φ j , b sin φ j , j = 1K N , ( are located at positions R j ,θ j ) with respect to the cohesive zone tip O . The cohesive zone starts at O ′ where the crack-opening displacement is calculated and compared to δ c as the condition for crack initiation. (b) Dugdale’s (rectangular) cohesive law with cohesive strength σ c , work of separation Γ , and critical displacement δ c = Γ σ c . 88 Consider a semi-infinite crack located along the negative x axis in an isotropic elastic ( ) medium (Figure 6.1a) with a plane strain modulus E ∗ = E 1 − ν 2 , E and ν being the Young’s modulus and Poisson’s ratio, respectively. The crack is subjected to a far-field loading characterized by a mode-I stress-intensity factor K I∞ . The rectangular cohesive law of Dugdale [52], with work of separation per unit area Γ , cohesive strength σ c and critical opening displacement δ c = Γ / σ c (Figure 6.1b), is assumed to govern the traction within the cohesive zone OO’ of length L (to be determined). The origin of the x − y coordinate system (point O) is located at the tip of the cohesive zone. We consider N edge dislocations ( j = 1,K, N ) with Burgers vectors (b cosφ , b sin φ ) j j located at distances R j from the cohesive zone tip and inclined at angles θ j with respect to the crack plane. Following Dugdale [52], crack initiation is assumed to occur when the crack opening displacement δ T at the beginning of the cohesive zone O ′ exceeds the range of cohesive interaction δ c , i.e. when δ T ≥ δ c , and the length L of the cohesive zone is determined based on the removal of singularity at O . The stress intensity factor at O and the crack opening displacement at O ′ are calculated as K I∞ + K Ic + K I⊥ = 0 , (6.2a) Γ ΔuO∞′ + ΔuOc ′ + ΔuO⊥′ = δ c = , (6.2b) σc where superscripts " ∞ , c , ⊥ ” denote, respectively, contributions from (i) the far-field loading characterized by stress-intensity factor K I∞ , (ii) the cohesive tractions acting inside the cohesive zone OO ′ (Fig. 1a), and (iii) the dislocation-crack interaction. 89 We now evaluate the various terms in Equation (6.2). The crack opening displacement at O′ due to the far field stress-intensity factor K I∞ is (e.g. [74]) 8K I∞ L ΔuO∞′ (− L) = . (6.3) E ∗ 2π The mode I stress intensity factor and the crack-opening displacement for a semi-infinite crack subjected to a distribution of normal traction p(x) [69] along the crack surfaces 0 2 p(t ) KI = π ∫ −t dt , (6.4a) −∞ 4 ⎛ 0 − x + − t ⎞⎟ ΔuO′ ( x) = p (t ) log⎜ ∫ dt . (6.4b) πE −∞ ∗ ⎜ − x − − t ⎟⎠ ⎝ In the Dugdale model, the cohesive traction is p c (t ) = −σ c for − L ≤ t ≤ 0 and vanishes over the rest of the crack faces. The corresponding stress intensity factor at O and crack opening displacement at O ′ are: 8L K Ic = −σ c , (6.5a) π 8σ c L ΔuOc ′ (− L) = − . (6.5b) πE ∗ The normal traction p ⊥ (t ) induced by N dislocations is the sum of tractions generated by the self stress field of each dislocation along the crack faces. For the j th dislocation at position (Rj, θj) (Figure 6.1a), the normal traction p⊥j (t ) is given by [70] p ⊥j = E ∗b 1 {− R [( cos φ j sin θ j t − R j cos θ j )2 − R j 2 sin 2 θ j ] [( 4π t − R cos θ j j 2 ) + R j 2 sin 2 θ j 2 ] j (6.6) ( + sin φ j t − R j cos θ j )([ t − R j cos θ j )2 + 3R j 2 sin 2 θ j ]} where b is the magnitude of the Burgers vector and φ j its orientation with respect to the 90 x-axis. We assume that the net interaction between the N dislocations and the cohesive crack is mode I dominant and neglect the effect of shear tractions. From Equations (6.4) and (6.6), the stress intensity factor and the crack opening displacement due to p ⊥ (t ) are found by the method of residues to be E ∗b ( ) N K I⊥ = ∑ f1 φ j , θ j (6.7a) j =1 4 2πR j where ( ) ( f1 φ j ,θ j = − sin θ j cos φ j − 3θ j 2 − 2 sin φ j cos θ j 2 , ) ( ) (6.7b) and ∑ f (φ ,θ , L R ) N b ΔuO⊥′ (−L) = (6.7c) π 2 j j j j =1 where ⎧⎪ ⎛ cos(θ j 2) ⎞ ⎛ cos(θ j 2) ⎞ ⎫ f2 (φj ,θ j , L Rj ) = sinφj ⎨tan−1⎜ ⎟ + tan−1⎜ ⎜ L R + sin(θ j 2) ⎟ ⎜ L R − sin(θ j 2) ⎟ j j( ⎟ −πH L R − sin θ 2 ⎪ ⎬ ( )) ⎪⎩ ⎝ j ⎠ ⎝ j ⎠ ⎪⎭ [ L Rj (L Rj )cos(φj −θ j 2) + cos(φj − 3θ j 2) sinθ j ] − (L R + cosθ ) + sin θ j j 2 2 j (6.7d) with ⎧0 L R j − sin θ j 2 < 0 H ( ⎪ L R j − sin θ j 2 = ⎨ ) . (6.8) ⎪⎩1 L R j − sin θ j 2 > 0 Generally speaking, the displacement field of a dislocation is not uniquely defined, because the displacement jump associated with a dislocation is defined uniquely only when the source position of that dislocation is known. In writing Equation (6.7c), we 91 have implicitly assumed there is no dislocation-induced displacement jump within the cohesive zone. If a dislocation-induced displacement jump occurs within the cohesive zone, our model will need to be modified; in particular, Equation (6.7c) will need to be modified by adding an additional contribution associated with the displacement jump. However, we argue that this case is physically unstable for most practical cases. For example, considering dislocation nucleation from a crack tip, it would be energetically more favorable to have the slip plane of the emitted dislocation intersect the end of, rather than within or near the tip of, a cohesive zone in order to avoid leaving behind a residual dislocation core which will be a high energy defect. This case is similar to that of dislocation nucleation from the tip of a diffusion wedge considered by Buehler et al. [10] where it was shown that dislocation nucleation from a diffusion wedge, which has the same stress singularity as a crack, is almost twice more difficult as that from a crack tip. We thus expect dislocation nucleation at the end of the cohesive zone, or, as occurs in atomistic simulations, that the crack and cohesive zone advance ahead of the source after the dislocation is nucleated. A similar situation occurs when a dislocation glides to a cohesive zone. The resulting high energy configuration can be easily released by having the cohesive zone advance ahead of the dislocation, so that the core energy is completely released. For the case of an edge dislocation directly ahead of the crack tip and created by some sort of climb process, the displacement jump corresponding to the extra plane of atoms would lie along the crack line and would physically correspond to a blunted crack. The two planes of atoms at the now-blunted crack tip would become potential cohesive planes, and the shielding effects due to stresses would be nearly identical to those 92 considered here – as if the dislocation had been moved just above or below the crack line by one plane spacing. This latter change would be quite negligible except when the dislocation is very, very close to the crack tip. Using Equations (6.3), (6.5) and (6.7), Equation (6.2) can be rewritten as E ∗b ( ) N 8L K I∞ − σ c +∑ f1 φ j , θ j = 0 , (6.9a) π j =1 4 2πR j 8K I∞ L 8σ c L b N Γ 2π − + ∑ f 2 (φ j ,θ j , L R j) = . (6.9b) E ∗ πE ∗ π j =1 σc We make two observations about Equation (6.9). First, for the cohesive crack in the absence of dislocations, we recover the classical solution of Dugdale [52] π ΓE ∗ LD = , K I∞ = E ∗Γ . (6.10) 8 σ c2 Second, we observe that all length scales in Equation (6.9) can be normalized by (E b ) Γ , ∗ 2 the cohesive strength σ c by Γ b , and the toughness K I∞ by E ∗Γ . Introducing the following dimensionless parameters, ΓR j ΓL σ b K I∞ Rˆ j = ∗ 2 , Lˆ = ∗ 2 , σˆ c = c , Kˆ I = , (6.11) Eb E b Γ E ∗Γ Equation (6.9) can be recast in the following form 8 Lˆ ( ) N 1 1 Kˆ I = π σˆ c − ∑ ˆ f1 φ j , θ j , (6.12a) 4 2π j =1 Rj 8Kˆ I σˆ c 8Lˆ ∑ f (φ ,θ , Lˆ / Rˆ ), σˆ N 1= Lˆ − σˆ c2 + c (6.12b) 2π π π j =1 2 j j j 93 which can be solved simultaneously for Lˆ and Kˆ I , given Rˆ j and σˆ c . Therefore, as the first main result of our study, we identify that the solution to Equation (6.9) can be completely characterized in terms of 3 dimensionless parameters Rˆ j , σˆ c , Kˆ I (after Lˆ is eliminated). 6.1.2 Limiting solutions for high and low cohesive strengths In the limit of high cohesive strength σˆ c → ∞ , we have vanishing cohesive zone length Lˆ → 0 and the Dugdale cohesive crack solution in Equation (6.12) reduces to the Lin- Thomson solution for a singular crack in Equation (6.13) ( ) N 1 1 Kˆ ILT = 1 − ∑ f1 φ j , θ j . (6.13) 4 2π j =1 Rˆ j In order to investigate how a finite σˆ c changes the behaviour of dislocation-crack interaction, let us first consider the limit of low cohesive strength σˆ c → 0 . In this limit, the cohesive zone size becomes very large, as shown in Equation (6.10), while the product Lˆ σˆ c2 remains finite, in which case Equation (6.12) is reduced to 8Lˆ N Kˆ I = σˆ c − 1 ∑ 1 ( ) f1 φ j , θ j , (6.14a) π 4 2π j =1 Rˆ j 8 Kˆ I σˆ c 8 Lˆ 1= Lˆ − σˆ c2 . (6.14b) 2π π Solving Equations (6.14a) and (6.14b) simultaneously for Kˆ I leads to 2 1 ⎡N ⎤ Kˆ I = 1+ ⎢∑ 1 ( ) f1 φ j , θ j ⎥ . (6.15) σˆ c →0 32π ⎢ j =1 Rˆ j ⎥ ⎣ ⎦ 94 An important observation can be immediately drawn from Equation (6.15): in the low strength limit, dislocations always shield the crack irrespective of the sign of the Burgers vector, i.e. Kˆ I > 1 . Therefore, we have arrived at the second contribution and the main theme of the present work: in the limit of low cohesive strength, the interaction between dislocations and a cohesive crack is fundamentally different from the classical singular crack solution. 6.1.3 Numerical results for Dugdale cohesive crack To illustrate the transition between high cohesive strength (Lin-Thomson) and low cohesive strength shielding, we consider a few simple representative examples resulting in pure mode I loading on the crack. Figure 6.2a shows the geometry and the numerically calculated net shielding effect ΔKˆ I = Kˆ I − 1 of a single edge dislocation in front of a Dugdale crack and Figure 6.2b shows similar results for two edge dislocations located symmetrically on slip-planes inclined at ± 60 o with respect to the crack plane. The distance between the cohesive zone tip O and the dislocation is taken to be R = 100 nm, 10nm, 5nm, and the normalized fracture toughness Kˆ I is numerically solved as a function of the normalized cohesive strength σˆ c from Equation (6.12). For E ∗ = 157 GPa , b = 0.25 nm and Γ = 2 J m 2 , the corresponding values of Rˆ are {20.382, 2.0382, 1.0191}. Based on Equations (6.13) and (6.15), the high and low strength limiting solutions for the geometry in Figure 6.2a are 1 1 Kˆ I = Kˆ ILT = 1 − ; Kˆ I = 1+ . (6.16a) σˆ c →∞ 8πRˆ σˆ c →0 8πRˆ and for the geometry in Figure 6.2(b) are 95 3 9 Kˆ I = Kˆ ILT = 1 − sin θ cos(θ 2) ; Kˆ I = 1+ sin 2 θ cos 2 (θ 2) , (6.16b) σˆ c →∞ 8πRˆ σˆ c →0 8πRˆ (a) 0.05 100nm 0.02 10nm 5nm 0 ΔKˆ I = 0.001 σˆ c →0 0.01 ΔKˆ ILT = −0.044 -0.05 ΔKˆ I -0.1 y b − 0.14 -0.15 O R x − 0.198 -0.2 -5 -4 -3 -2 -1 0 1 10 10 10 10 10 10 10 (b) σˆ c = π 8 Lˆ D ( ) 0.2 100nm 10nm 0.1 0.094 5nm 0.05 0 Δ Kˆ I = 0.005 σˆ c →0 ΔKˆ ILT = −0.1 -0.1 ΔKˆ I -0.2 y b -0.3 R − 0.314 60o O 60o x -0.4 R − 0.44 b -0.5 -5 -4 -3 -2 -1 0 1 10 10 10 10 10 10 10 σˆ c = π 8 Lˆ D ( ) 96 Figure 6.2: Dislocation shielding ΔKˆ I of a Dugdale (rectangular) cohesive crack by (a) a single edge dislocation in front of the crack-tip and (b) two edge dislocations symmetrically located on a slip-plane inclined at ± 60o with respect to the cohesive zone tip and crack plane. The results are plotted as a function of the normalized cohesive strength σˆ c for different distances R = 5nm, 10nm, 100nm between the dislocation and the cohesive zone tip O . For high cohesive strengths, the cohesive crack behaves like a singular crack and the selected dislocation structures anti-shield the crack, while for very low cohesive strengths, the same dislocation structures act to shield the cohesive crack with the value evaluated through Equation (6.16) [103]. Figures 6.2a and 6.2b clearly show a transition in the net shielding between the Lin- Thomson solution, which is antishielding in all of these cases, to the low strength solution in which a dislocation always shields the crack. For sufficiently high strengths σˆ c ≥ 1.0 , dislocations affect the toughness of a cohesive crack in a similar way as a singular crack. A normalized cohesive strength of σˆ c ≈ 0.1 is selected for characterizing the transition from high strength (singular crack) to low strength regimes. For b = 0.25 nm and Γ = 2 J m 2 , σˆ c = 0.1 corresponds to σ c = 8 00 MPa . This indicates that a fairly high cohesive strength will be required for a cohesive crack to behave as a singular crack in the sense of Lin-Thomson [66] analysis. The generality of this conclusion will be demonstrated by some practical DD simulations shortly. Also, as R increases the cohesive crack behaviour remains similar to the Lin-Thomson solution over a wider range of cohesive strengths. This trend is consistent with our expectations that, as R increases, the character of the crack tip (singular or cohesive) become less important. We may also interpret Figures 6.2a and 6.2b in terms of the characteristic cohesive zone length LD defined in Equation (6.10) if we express the normalised cohesive strength as 97 π σˆ c = (6.17) 8LˆD where L Γ LˆD = D∗ 2 , (6.18) Eb is similarly defined as Lˆ in Equation (6.11). Therefore, the normalized cohesive strength σˆ c is directly related to the normalized cohesive zone length Lˆ D . Hence, Figures 6.2a and 6.2b also shows that the cohesive crack behaviour tends towards that of a singular crack as Lˆ D decreases with respect to the characteristic length (E ∗b 2 ) Γ . 1.5 100nm 1.45 10nm 5nm 1.4 1.35 L 1.3 L D 1.25 y 1.2 b 1.15 O R x 1.1 1.05 1 -5 -4 -3 -2 -1 0 1 10 10 10 10 10 10 10 ( σˆ c = π 8 Lˆ D ) Figure 6.3: The ratio between the actual cohesive zone length L and the characteristic cohesive zone length LD defined in Equation (6.10) for a Dugdale (rectangular) cohesive crack in the presence of a single edge dislocation in front of the crack-tip. The results are plotted as a function of the normalized cohesive strength σˆ c for different distances R = 5nm, 10nm, 100nm between the dislocation and the cohesive zone tip O . Interestingly, L remains comparable in magnitude to LD for the full range of cohesive strengths and approaches LD in the high strength limit [103]. 98 Figure 6.3 shows the variation of the normalized cohesive zone length L LD of the Dugdale crack by a single edge dislocation in front of the crack-tip with respect to the normalized cohesive strength σˆ c . Interestingly, the actual cohesive zone length L remains comparable in magnitude to the characteristic cohesive length LD defined in Equation (6.10) over a wide range of cohesive strength, and approaches LD in the high strength limit. -0.05 -0.1 -0.15 ΔKˆ I -0.2 y ΔKˆ exact -0.25 b b ΔKˆ appx O R1 R2 x ΔKˆ LT -0.3 0 100 200 300 400 500 R2(nm) Figure 6.4: Variation of the net shielding ΔKˆ I of two dislocations in front of a crack-tip with cohesive strength σ c = 1GPa as a function of distance R2 of the second dislocation, while the first dislocation remains fixed at a distance R = 10 nm . Here, ΔKˆ 1 exact is the exact solution, ΔKˆ appx is calculated by superposing toughness contributions from each individual dislocation treated separately and ΔKˆ LT is the prediction from the Lin-Thomson solution [103]. Equations (6.12) are evidently non-linear, so superposition of the shielding from multiple dislocations does not hold, in contrast from the singular crack limit where 99 elasticity and superposition do hold. To examine superposition, we consider two collinear dislocations in front of the crack, as shown in Figure 6.4. We vary the position of the second dislocation while keeping the cohesive strength and the distance between the first dislocation and the origin O fixed at σ c = 1GPa, R1 = 10 nm . Figure 6.4 shows ΔKˆ exact , ΔKˆ appx and ΔKˆ LT , corresponding to the exact solution to the normalized fracture toughness in the presence of two dislocations, an approximate solution obtained by superposing the toughness due to each dislocation treated separately, and the classical Lin-Thomson solution, respectively. When the second dislocation is far away from the crack, its effect is small relative to the first dislocation. As the second dislocation is brought closer to the crack, the reasonable match between ΔKˆ exact and ΔKˆ appx suggests that the superposition principle can still be used to characterize multiple dislocations interaction with a cohesive crack with reasonable accuracy. For dislocations close to the crack, i.e. in Figure 6.4 when the second dislocation is at R2 < 50nm, superposition clearly fails. 6.2 Discrete Dislocation Simulations ‡ Our analysis based on the Dugdale cohesive crack suggests that the cohesive strength σ c plays a critical role in the transition between the Lin-Thomson solution and the low strength solution given in Equation (6.15). Although our analytical solutions have been derived based on the Dugdale cohesive crack model, we expect that the basic conclusions should be generally valid for dislocations interaction with cohesive cracks. The generality ‡ The DD simulation results reported in this chapter are carried out by Dr. Audrey Chng from Institute of High Performance Computing, Singapore and included here for the sake of completeness. 100 of our analysis and the practical significance of our findings is demonstrated through a comparison with Discrete Dislocation (DD) simulation [103], where a cohesive crack described by the trapezoidal traction-separation law of Tvergaard and Hutchinson [54] (TH) interacts with a large number (1200) of immobile dislocations following the DD framework developed by van der Giessen and Needleman [101]. In the simulation, a dislocation structure formed after some crack growth under a particular cohesive strength is frozen and further crack growth is examined as a function of variation in cohesive strength under constant cohesive energy [103]. Figure 6.5 shows the effect on the net dislocation shielding ΔKˆ of lowering the cohesive strength σ c of the TH trapezoidal cohesive law from 7.2GPa, where the LT ( ) result ΔKˆ LT = 0.2662 is accurate, down to 0.9GPa under constant cohesive energy. Raising cohesive strength has the effect of decreasing the difference between the cohesive ( ΔKˆ ) and singular ( ΔKˆ LT ) crack predictions, but ΔKˆ − ΔKˆ LT remains significant even at σ c = 3.6GPa. The dislocation closest to the origin is located at a distance Rmin = 62.5nm from the cohesive zone tip. Estimating the length of the cohesive zone using L D = π Γ E 8σˆ 2 , Figure 6.5 suggests Rmin / LD > 50 before ΔKˆ → ΔKˆ LT in the presence of a large number of dislocations. These conclusions involving a large number of dislocations interacting with a TH cohesive zone are consistent with those for a single dislocation or a pair of dislocations interacting with the Dugdale cohesive zone. In general, regardless of the number of dislocations or the exact form of traction- separation law, the cohesive crack tends to the singular behaviour as cohesive strength increases or as dislocations move away from the crack tip. 101 Based on the asymptotic solutions of dislocations interaction with a cohesive crack in the high and low strength limits, a simple asymptotic interpolation scheme is adopted in the Appendix G to derive a general approximate solution of dislocations interaction with a cohesive crack. Predictions of the approximate solution are also shown in Figure 6.5, and the agreement with the DD simulations is quite good. Our asymptotic solution significantly improves upon the Lin-Thomson solution in providing a more reasonable estimate of dislocation shielding of a cohesive crack for all cohesive strengths. 0.7 DD simulation 0.6 Asymptotic Lin-Thomson 0.5 0.4 ΔKˆ I 0.3 0.2 0.1 0 0 2 4 6 8 σc (GPa) Figure 6.5: Simulated dislocation shielding ΔKˆ of a trapezoidal cohesive crack of Tvergaard and Hutchinson [54] at various cohesive strengths σ c . The Lin-Thomson solution predicts a net shielding of ΔKˆ LT = 0.2662 . Although ΔKˆ tends to ΔKˆ LT at very high σ c , significant ( ) deviations ΔKˆ − ΔKˆ LT exists even at σ c = 3.6GPa. The DD simulation results [103] are compared with the asymptotic interpolation solution derived in the Appendix G, with the net shielding obtained by adding contribution from individual dislocations given in Equation (G.6), where the cohesive zone size L is taken to be the characteristic cohesive zone size LD (Equation 6.10). The elastic material is assumed to have Young's modulus E = 140GPa, Poisson's ratio ν = 0.33, and plane strain modulus E * = 157GPa. The TH trapezoidal cohesive law has constant energy Γ = 1.125J/m2, δ 1 = 0.25δ c , δ 2 = 0.5δ c , and fracture toughness KI,c= E *Γ = 0.42 MPa m [103]. 102 It is important to note that our comparison above, is based on a frozen dislocation structure. The toughening that emerges from a DD simulation involves the evolution of a dislocation structure in the presence of a particular cohesive zone. Thus, in addition to differences in shielding for a particular dislocation structure, changing the cohesive zone strength may also change the dislocation structure itself, and thus modify the overall toughness. In previous DD modelling of fracture in thin films, Chng et al. [64] found that the overall toughness was nearly independent of cohesive zone strengths in the range 0.9GPa to 1.8 GPa. Chapter 7 8. Concluding remarks A fundamental understanding of inelastic deformation mechanisms and stresses in small structures is essential in view of their importance in nearly all modern technologies. The subject is not only important from the standpoint of reliability but also from the point of understanding material behavior as length dimension gets smaller. In the work presented in this thesis, we have tried to analyze material response at small length scales for a number of systems. The main conclusions of the work are summarized as follows: Stress evolution in Volmer-Weber thin film growth: The thin film growth model developed here is based on a rigorous integro-differential equation formulation of constrained GB diffusion [6-7] and follows the proposed mechanism of excess surface chemical potential during deposition which drives adatoms into grain boundaries [26-28]. 1. A numerical algorithm has been developed to treat the hypersingular integro- differential equations that arise as a result of inhomogeneous GB diffusivity in the thin film growth model. 103 104 2. Crack-like GB diffusion wedges are shown to form during the thin film growth. The stress-intensity factor near the tip of a wedge in a slowly growing high mobility film is comparable to that of a GB wedge under annealing conditions, and can affect the stress evolution in the film during growth through additional stress relaxation. 3. In the limit of infinite GB diffusivity, an analytical model has been developed to describe compressive stress evolution which shows excellent agreement with numerical solution of the problem as well as experimentally measured Stress- Thickness evolution during periodic deposition and growth interruption of Sn films on Si/Sn substrates [28]. 4. The finite and infinite GB mobility simulation results have successfully reproduced the experimentally observed effects of GB mobility and film thickness on transient Stress-Thickness behaviors during both deposition and interruption stages, as well as the steady state stress as a function of deposition rate and film thickness. 5. Apart from the mechanism proposed by Chason et al. [26-28], other important competing theories for the post coalescence compressive stress evolution during high mobility thin film deposition include (a) surface stress effect due to change in adatom population during growth [20-21]; (b) trapping of adatoms in interstitial sites near step edges and ledges [25]. Although these mechanisms qualitatively explain the post coalescence compressive stress and subsequent relaxation during growth interruption, they lack the predictive features and have so far failed to quantitatively explain experimental observations on the effects of GB mobility, 105 deposition rate and film thickness on the stress evolution and recovery behaviors. In contrast, the success of constrained GB diffusion model in providing a unified explanation of stress-temperature behaviors [11] and parallel glide dislocations observed in thin films under thermal cycles [12], and Stress-Thickness behaviors during film growth and stoppage [17,26,79] provides a strong evidence that atomic transport between film surface and GBs may play a crucially important role in stress evolution during thin film deposition. 6. Mechanisms like dislocation activity and bulk diffusion in the interior of crystal [104] and surface diffusion [7,82] not treated in our model can also influence the material response. Atomistic simulations to study surface chemical potential in the presence of deposition may help to provide a better physical representation for surface kinetic parameters A , C s Γs and C d Γd in our continuum model. These topics provides avenue for further research in the topic of stress evolution during deposition. Crack nucleation under diffusion induced stresses in battery electrodes: We have developed a cohesive model of crack nucleation in an initially crack-free electrode under galvanostatic intercalation and deintercalation processes for strip and cylindrical geometry. 1. Our analysis shows that nanoscale size is indeed a key to suppressing crack formation under large diffusion induced stresses in high capacity battery electrodes such as Si. The results for cylindrical electrode are similar to the strip one, implying relatively small influence from electrode geometry. 106 2. An important outcome of our analysis is a critical length scale { [ H ft = 13 Γ(1 − ν )F 2 D 2 E (1 + ν )Ω 2 I 2 ]} 13 for the electrode to remain flaw tolerant. When the characteristic dimension of the electrode is below this length scale, crack nucleation becomes impossible irrespective of the cohesive strength of material. 3. Our model shows that the critical length scale for flaw tolerance should scale inversely with the density of charging current. The faster the charging, the higher the stress and the more likely the crack will be nucleated. This is consistent with experimental observations. Based on the flaw tolerant dimension H ft , a potential design criterion for flaw tolerant electrode is I 2 / 3 H ≤ crit , where I is the operating current and H the dimension of electrode. For a “fail-safe” design of electrode, one might first determine the allowable value of I 2 / 3 H by a standard testing protocol and then design the actual electrode dimension and current density according to the scaling relation. 4. We would like to point out a number of limitations in our current model. First, we have only considered highly idealized electrode geometry. Even for such simple geometry, we have not yet investigated the stability of the crack after it is nucleated, where it might be necessary to consider different boundary conditions on the newly formed crack surfaces (which are reported to form solid electrolyte interphase and degrade the battery performance). Second, the effect of hydrostatic stress can influence the chemical potential [48,77,92] as well as the activation energy for diffusion [105], introducing strongly non-linear terms in the governing equation and boundary conditions. The full coupling between stress and diffusion 107 is expected to be important to model electrode in a strongly confined environment. Once the coupling is introduced, the crack nucleation problem will have to be treated in a more sophisticated nonlinear framework, which would substantially increase the complexity of the problem. Third, the current formulation is strictly valid only in the limit of dilute solutions since otherwise the distribution of Li ions may significantly alter the physical properties of the material such as the elastic modulus. Fourth, the lithiation induced phase changes needs to be further studied. Although not treated in the current study, phase changes during lithiation are another major cause for stresses in addition to DIS. The micrometer length of electrode is important in this context, as a critical thickness of 2 μm is required for the crystallization of amorphous Si to crystalline Li15Si4 [106]. All these issues may significantly affect the electrode failure process in some way. Nevertheless, the simple model adopted here may serve as a first step to bring out some essential features of an appropriate length scale and the associated scaling laws for crack nucleation. Dislocation shielding of cohesive cracks We have conducted analytical and numerical studies on dislocations shielding of a Dugdale cohesive crack. 1. Our theoretical analysis shows that the normalized toughness Kˆ I = K I∞ E ∗ Γ is completely characterized by the normalized distance of dislocations from the cohesive zone tip Rˆ j = ΓR j E ∗b 2 and the normalized cohesive strength σˆ c = σ c b Γ . 108 2. We have derived closed form asymptotic solutions in the limits of high and low cohesive strengths, and constructed an approximate general solution based on an asymptotic interpolation scheme. 3. Based on these solutions, we have shown a surprising result that, while the classical singular crack model predicts that a dislocation shields or anti-shields a crack depending on the sign of the Burgers vector, at low cohesive strengths a dislocation always shields a crack irrespective of the Burgers vector. Thus, there exist qualitative differences between the cohesive and singular crack solutions. Our results suggest that the classical Lin-Thomson solution serves as an upper bound to the amount of anti-shielding for a cohesive crack. 4. The practical significance of our findings is demonstrated through a comparison of the asymptotic solution with DD simulation §, where the net shielding effect of a large number ( > 10 3 ) of edge dislocations on a cohesive crack described by the Tvergaard-Hutchinson trapezoidal traction-separation law is analyzed. The simulated shielding showed increasing deviations from the classical Lin-Thomson predictions as we reduce the cohesive strength σ c from 7GPa to around 1GPa, indicating increasing importance of the crack tip character as the cohesive zone size exceeds about 2% of the distance of the nearest dislocation from the cohesive zone tip. 5. While the analysis deals with Dugdale and trapezoidal cohesive laws, we expect that the basic conclusions to remain true for other types of cohesive laws. § DD simulations are carried by Dr. Audrey Chng, IHPC, Singapore. Appendix A A. Numerical algorithm for hypersingular integral equations We describe below the numerical algorithm for solving the coupled integral equations described by Equations (3.24) and (3.25). As discussed in Section 3.3 of Chapter 3, that the GB wedge profile is expected to adopt a crack-like opening profile near z = h , which implies that ∂uˆ ∂x and ∂σˆ ∂x will have the asymptotic behaviors ∂uˆ ∂x ~ (1 − x) −1 / 2 and ∂σˆ ∂x ~ (1 − x) 3 / 2 as x → 1 . Following [84], we can express the unknown function ∂σˆ ∂x as ∂σˆ ( ) ∞ ( x, τ ) = ∑ d jU j ( x ) 1 − x 2 3/ 2 g ( x) (A.1) ∂x j =0 where U j ( x) = sin (( j + 1) x ) sin( x) are Chebyshev polynomials of the second kind. The series is truncated at N terms in the actual numerical implementation. The unknown function ∂uˆ ∂x can be expressed as ∂uˆ ( x,τ ) Ru ( x,τ ) = (A.2) ∂x 1− x2 109 110 where Ru ( x,τ ) is a smooth function in [0,1] . For a stable algorithm, an implicitly discretized form of ∂σˆ ( x,τ ) ∂τ and ∂uˆ ( x,τ ) ∂τ are approximated by ∂σˆ ( x,τ + Δτ ) σˆ ( x,τ + Δτ ) − σˆ ( x,τ ) ≈ , (A.3) ∂τ Δτ ∂uˆ ( x,τ + Δτ ) uˆ ( x,τ + Δτ ) − uˆ ( x,τ ) ≈ . (A.4) ∂τ Δτ Combining Equations (3.24-3.27) and (A.1-A.4), we have H& Δτ dη N −1 1 ∑d ja j + 2 ∫ H (τ + Δτ ) 0 Q( x,η ) Ru (η ,τ + Δτ ) = σˆ ( x,τ ) , (A.5) j =0 1 −η 2 dη H& Δτ 1 N −1 1− x ∑d b − ∫ R j j u (η ,τ + Δτ ) + H (τ + Δτ ) 1 + x Ru ( x,τ + Δτ ) = uˆ ( x,τ ) , (A.6) j =0 x 1 −η 2 where aj = ∫ x ( U j (η ) 1 − η 2 ) 3/ 2 dη + U j ( 0) + H& Δτ (1 − x) 5 / 2 (1 + x) 3 / 2 U j ( x) 0 g (η ) AH (τ + Δτ ) g (0) H (τ + Δτ ) g ( x) Δτ ∂ 2 P ( x ,η ) Δτ ∂P( x,0) U j (0) 1 H (τ + Δτ ) 3 ∫0 ∂ 2η − U (η )(1 − η 2 3/ 2 ) d η − , (A.7) H (τ + Δτ ) 3 ∂η j g (0) bj = Δτ H (τ + Δτ ) 2 { [ ] [ − 2 x sin ( j + 1) cos −1 ( x) − ( j + 1) 1 − x 2 cos ( j + 1) cos −1 ( x) . ]} (A.8) Note ∂ 2 P( x,η ) ∂η 2 can be decomposed into hypersingular and regular term as ∂ 2 P ( x,η ) 2 ⎡ ∂ 2 P ( x ,η ) 2 ⎤ = + ⎢ − ⎥. (A.9) ∂η 2 (x −η) 3 ⎣ ∂η 2 (x −η )3 ⎦ Hypersingular Regular 111 Integrals involving regular term in Equation (A.7) are computed using the Gauss- Legendre quadrature [72]. The hypersingular integral in Equation (A.7), 1 U j (η ) ∫ (x −η) 0 3 (1 − η 2 ) 3 / 2 dη , j = 0,1,.., N − 1 , is evaluated in closed form using the analytical formulae derived in [84]. Integrals involving Ru ( x,τ ) in Equations (A.5) and (A.6) are evaluated using the Gauss- Chebyshev quadrature [68] by extending the integral range [0,1] to [− 1,1] via even continuation Ru (− x,τ ) = Ru ( x,τ ) x ∈ [0,1] . Using 2 N + 1 points Gauss-Chebyshev quadrature, we have dη dη 1 2 N +1 π 1 1 1 ∫0 K ( x,η ) Ru (η ,τ ) 1 − η 2 = 2 −∫1 K ( x,η ) Ru (η ,τ ) 1 − η 2 = 2 ∑i =1 2 N + 1 K ( x, η i ) Ru ( η i ,τ ) (A.10) where the integration points are given by ⎡ (i − 0.5)π ⎤ η i = cos ⎢ i = 1,2,..., 2 N + 1 . (A.11) ⎣ 2 N + 1 ⎥⎦ In Equation (A.5), K ( x,η ) = Q( x,η ) and Q( x,0) = 0 ; while in Equation (A.6) ⎧0 η