Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Determination of stress intensity factors for laboratory test specimens using numerical methods Ge, Baolai 1994

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1994-0435.pdf [ 3.33MB ]
Metadata
JSON: 831-1.0080857.json
JSON-LD: 831-1.0080857-ld.json
RDF/XML (Pretty): 831-1.0080857-rdf.xml
RDF/JSON: 831-1.0080857-rdf.json
Turtle: 831-1.0080857-turtle.txt
N-Triples: 831-1.0080857-rdf-ntriples.txt
Original Record: 831-1.0080857-source.json
Full Text
831-1.0080857-fulltext.txt
Citation
831-1.0080857.ris

Full Text

DETERMINATION OF STRESS I N T E N S I T Y FACTORS FOR LABORATORY TEST S P E C I M E N S USING NUMERICAL METHODS By Ge Baolai B. Sc. (Mechanics), Shanghai Jiaotong University, P. R. China, 1984 M. Sc. (Mathematics), Lakehead University, Canada, 1992 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 1994 © Ge Baolai, 1994 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mechanical Engineering The University of British Columbia 2324 Main Mall Vancouver, Canada V6T 1Z4 Date: / / -7 / 0o£ / n «H A b s t r a c t Using the Airy stress function, plane bnear elasticity problems reduce to solving a biharmonic problem of the first type. Analytic solutions to the biharmonic problems of the first type in a domain with a re-entrant sharp corner are discussed in this thesis. General solutions to the biharmonic equation in the vicinity of the sharp corner are given in terms of a series of some special eigenfunctions, which are obtained from two auxiliary eigenvalue problems corresponding to symmetric and anti-symmetric stress states. In the case of a line crack, the leading term in each series contains the singularity in stresses near the crack tip, and its coefficient determines the stress intensity factor of Mode I or Mode II loading. In this thesis, stress intensity factors are evaluated with the specimens of two simple ge-ometries - circular disc and rectangle - that are subjected to three types of loadings: uniformly distributed forces, concentrated forces and those yielding pure bending. Stress fields are found using the series stress functions with the boundary collocation method. The least squares approach is applied to the system of linear equations arising from the boundary collocation procedure. Compared to other sophisticated methods, the present approach employing the technique of boundary collocation is straightforward and relatively efficient in numerical computation. Numerical values of stress intensity factors are compared to the results obtained using other less efficient methods. A sound and overall agreement is found in most cases. The approach presented in this thesis along with the results may be used as an auxiliary tool to determine the fracture toughness of laboratory test specimens. n Table of Contents A b s t r a c t ii Table of Content s iii List o f Tables vi List o f F igures viii A c k n o w l e d g e m e n t xii 1 In troduct ion 1 1.1 Crack and Stress Singularities in Plane Elasticity Problems 1 1.2 Griffith Crack Problems 3 1.3 Literature Review 7 1.4 The Objectives 10 1.5 Scope of The Thesis 11 2 M a t h e m a t i c a l Foundat ions 13 2.1 Series Solution of Biharmonic Equation 13 2.2 Stress Functions for The Problems with V-Shaped Boundary 14 2.3 The Fourier Expansion Method 19 2.4 Discussion on The Fourier Expansion Method 20 3 Numer ica l Approaches and Analys is 24 3.1 The Localness of The Basic Solutions 24 3.2 Numerical Approaches 25 iii 3.3 Error Measurements 27 3.4 Boundary Collocation Method 28 4 N u m e r i c a l E x p e r i m e n t s 31 4.1 Notation and Conventions 32 4.2 Comparison of Numerical Results with Previous Work 35 4.3 Rectangular Plate under Uniform Tension 44 4.4 Rectangular Plate under Concentrated Loading 52 4.5 Rectangular Plate under In-plane Pure Bending 66 4.6 Rectangular Plate under Concentrated Shearing 72 5 Discuss ion and Conclusion 78 Bibl iography 82 A p p e n d i c e s 85 A T h e Special Crack Stress Function and S t r e s s e s 85 A. I Stress Functions - Mode I 85 A.2 Stress Functions - Mode II 86 A.3 Stresses in Cartesian and Polar Coordinates 87 B A n Al ternat ive Form of Boundary Condi t ions 89 B.l The General Form of Biharmonic Problems of The First Type 89 B.2 A Rectangular Plate under General Loadings 91 B.3 Rectangular Plate under Uniform tension 95 B.4 Circular Disc with Radial Crack 96 C T h e Auxi l iary Eigenvalue Prob lems 98 D Coupled Poisson's Equat ions 101 iv E The Least Squares Solution of Linear Systems 105 v List of T a b l e s 4.1 Ten leading coefficients with using 20, 40 and 80 truncated terms (circular disc, uniform tension) 37 4.2 Comparison of leading coefficients with Wu's results [37] (circular disc, uniform tension) 37 4.3 Ten leading coefficients with using 100, 140 and 200 truncated terms (circular disc, diametric loading) 41 4.4 Comparison of leading coefficients with Wu's results [37] (circular disc, diametric loading) 41 4.5 Ten leading coefficients with using 150, 180 and 200 truncated terms (circular disc, crack opening loading) 43 4.6 Comparison of leading coefficients with Wu's results [37] (circular disc, crack opening loading) 43 4.7 Ten leading coefficients with using 30, 40 and 80 truncated terms (rectangle, uniform tension) 49 4.8 Comparison of boundary conditions at selected boundary points (rectangle, uni-form tension) 51 4.9 Ten leading coefficients with using 25, 70 and 90 truncated terms (rectangle, loaded at crack opening) 58 4.10 The agreement of boundary conditions at selected boundary points (rectangle, loaded at crack opening) 59 4.11 Ten leading coefficients with using 60, 80 and 100 truncated terms (rectangle, loaded at left corners) 64 VI 4.12 The agreement of boundary conditions at selected boundary points (rectangle, loaded at left corners) 65 4.13 Ten leading coefficients with using 34, 40 and 80 truncated terms (rectangle, pure bending) 70 4.14 The agreement of boundary conditions at selected boundary points (rectangle, pure bending) 71 4.15 Comparison of boundary conditions at selected boundary points (rectangle, pure shearing) 76 4.16 Ten leading coefficients with using 20, 30 and 36 truncated terms (rectangle, pure shearing) 77 C.17 The eigenvalues corresponding to Mode I and Mode II cases for selected vertex angles 99 E.18 Comparison of the least squares solution of dimensionless K\ and errors in the Mode I compact loading case for c/a = 0.5, h/a — 0.8, l/h — 0.25 using totally 100 terms and 159 collocation points 107 vn List of F i g u r e s 1.1 Elliptic crack in an infinite plate by Inglis 3 1.2 A line crack in an infinite plate 4 1.3 Section of a plate with a Griffith crack: Mode I displacement 5 1.4 Section of a plate with a Griffith crack: Mode II displacement 5 1.5 Section of a plate with a Griffith crack: Mode III displacement 6 2.6 Infinite domain with a F-shaped angular corner 15 2.7 Circular disc with a radial crack less than the radius. The origin of the polar coordinates is located at the crack tip 22 2.8 Rectangle with an edge crack 22 3.9 Motz' finite different mesh near the crack tip 0. with • representing near-points and o representing far-point. The crack opening is towards left. Grid points are not specified along the crack 26 4.10 Circular disc under uniform tension a0 35 4.11 The variation of K\ with the number of collocation points (circular disc, uniform tension) 36 4.12 The variation of K\ with the number of truncated terms (circular disc, uniform tension) 38 4.13 Root-mean-square error ||e|| with the number of collocation points (circular disc, uniform tension) 39 4.14 Maximum error HeHoo with the number of collocation points (circular disc, uni-form tension) 39 4.15 Circular disc under diametric loading 40 viii 4.16 Circular disc loaded at the crack mouth 42 4.17 Rectangular plate under uniform tension 44 4.18 Distribution of collocation points along the edges of cracked rectangle 45 4.19 The variation of dimensionless K\ with the number of collocation points (rect-angle, uniform tension) 46 4.20 The variation of dimensionless K\ with the number of truncated terms (rectangle, uniform tension) 47 4.21 The root-mean-square errors ||e|| with the number of collocation points (rectan-gle, uniform tension) 48 4.22 The maximum errors ||e||oo with the number of collocation points (rectangle, uniform tension) 48 4.23 The variation of dimensionless Ki with the crack length ratio c/a (rectangle, uniform tension) 49 4.24 Dimensionless tensile stresses ay on different sections (rectangle, uniform tension). 50 4.25 The isochromatics of ay (rectangle, uniform tension) 50 4.26 The standard compact specimen for the test of fracture toughness 52 4.27 The specimen and the mathematical model for concentrated loading at crack mouth 53 4.28 The specimen and the mathematical model for concentrated loading at upper side edge 53 4.29 The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.25 (rectangle, loaded at crack opening). . . 55 4.30 The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.5 (rectangle, loaded at crack opening). . . 55 4.31 The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 1 (rectangle, loaded at crack opening). . . . 56 ix 4.32 The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a — 0.5 (rectangle, loaded at crack opening). 57 4.33 The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a — 0.67 (rectangle, loaded at crack opening). 57 4.34 Dimensionless tensile stresses on different sections (rectangle, loaded at crack opening) 58 4.35 The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.25 (rectangle, loaded at left corners). . . . 60 4.36 The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.5 (rectangle, loaded at left corners) 61 4.37 The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 1 (rectangle, loaded at left corners) 61 4.38 The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a — 0.5 (rectangle, loaded at left corners). 62 4.39 The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a = 0.67 (rectangle, loaded at left corners). . 62 4.40 The dimensionless tensile stresses on different sections (rectangle, loaded at left corners) 63 4.41 Rectangular plate with edge crack subjected to pure bending 66 4.42 The variation of dimensionless K\ with the number of collocation points (rect-angle, pure bending) 68 4.43 The variation of dimensionless K\ with the aspect ratio h/a for fixed crack length ratio c/a = 0.5 (rectangle, pure bending) 68 4.44 The variation of dimensionless K\ with crack length ratio c/a for h/a — 0.5 and h/a > 1 (rectangle, pure bending) 69 4.45 Dimensionless tensile stresses ay on different sections (rectangle, pure bending). . 70 4.46 A simple specimen for Mode II test and numerical model 72 4.47 The variation of K\\ with the crack length ratio cja for different aspect ratios (rectangle, pure shearing) 74 4.48 The variation of A'n with the number of collocation points for cja = hja = 0.5 (rectangle, pure shearing) 75 4.49 The isochromatics of shear stresses around the crack tip and the corners of loaded edges (rectangle, pure shearing) 75 B.50 An arbitrary boundary C of certain region Q in x ~- y plane 89 B.51 A rectangular plate with an edge crack OA perpendicular to the edge AB 92 C.52 Variation of minimum real positive eigenvalues with the vertex angle for sym-metric (Mode I) and antisymmetric (Mode II) cases 99 D.53 Finite difference grid for the upper half cracked rectangle OABCDO. The crack is along OA and OD is the symmetry line. The boundary points are labeled with o and interior points are labeled with • 102 XI A c k n o w l e d g e m e n t I would like to thank Professor Henry Vaughan for his guidance and encouragement during the research and the writing of this thesis. I appreciate many interesting and heuristic discussions with him in the past two years. Also I wish to thank the Natural Science and Engineering Research Council of Canada for their financial support. XI1 Chapter 1 Introduct ion The interest in crack problems in the mathematical theory of elasticity arises from the theory of brittle fracture. Although numerical methods such as the finite element and finite difference techniques have been widely used in solving various complicated problems, analytical solutions are still of significance when high accuracy and numerical efficiency are concerned in solving crack problems. 1.1 Crack and Stress Singularit ies in Plane Elasticity P r o b l e m s In two dimensional elasticity theory, the stresses satisfy the equilibrium equations ^ + ^ = 0 . (1.2) ox ay By introducing an Airy stress function <j), these equations are identically satisfied if letting d24> d24> d*4> Imposing also the following compatibility condition V 2 ( a x + a y ) = 0 (1.4) where V 2 is the usual Laplace or harmonic operator, we have the following biharmonic equation in 4> V V = 0 . (1.5) 1 Chapter 1. Introduction 2 Given stresses on the boundary, solving for the stress field is equivalent to solving the following biharmonic problem of the first type V V = 0 in ft, 4> = f on dSl, (1.6) 4>n ~ 9 °n dil where / and g are determined from given stresses conditions and cf>n is the normal derivative of (j> outwards to the boundary dil. This form is sometimes more convenient for the sake of mathematical analysis and numerical computations. Based on harmonic analysis, it is known that the solution of the second order elliptic problem is analytic provided the boundary and boundary values of (j> are both analytic. Singularities can occur only when the boundary or some part of its data is not "smooth". Similar conclusions can be extended to the biharmonic problem. The second derivatives of <f> (even its first derivative) and in turn the stresses may become unbounded in the neighbourhood of crack tip, where geometric singularity occurs. In fracture mechanics, the first mathematical approach to crack problems was proposed by lnglis[9] in early 1900s. By considering an elliptic hole inside an infinite plate under uniform tension a (Fig.1.1), Inglis showed that the maximum stress in the plate is at the extremity of the major axis and is given by (°y)m.x = o(l + 2a/b) (1.7) where a and b are respectively the major and minor axes of the ellipse. As b tends towards zero, Equation (1.7) may be written in terms of the curvature as follows / 2a\1< / 2 K U x ^<r\l + —) (1-8) where p — b2/a is the radius of curvature of the ellipse at the apex of the major axis. As b tends towards zero, p also approaches to zero and thus K)m a x*2a(-J (1.9) Chapter 1. Introduction 3 t t t t t t t t t t Figure 1.1: Elliptic crack in an infinite plate by Inglis This indicates that for a line crack (see Fig. 1.2), the tensile stress at the crack tip becomes infinitely large. 1.2 Griffith Crack P r o b l e m s The Griffith crack is of primary importance in fracture mechanics. Consider an elliptic crack of length 2c in an infinite plate under tensions that are perpendicular to the major axis of the ellipse (Fig. 1.1). The surface of the crack is free from the applied stress and there is a prescribed tensile stress a at infinity so that r x y ( x , 0 ± ) = ay(x,0±) - 0, |x| < c (1.10) and ax(x,y) ~> 0, Txy(x,y) -* 0, ax(x,y) -» a, r -•-> oo (1.11) where r — (x2 -f- y2)1'2. In practice we are interested in the more general case in which a at infinity is not constant. It is interesting to note that making use of the superposition principle Chapter 1. Introduction 4 t t t t t t t t t t U H I TIT I I Figure 1.2: A line crack in an infinite plate of linear elasticity theory, this problem is equivalent to the one below Txy(x,0±) = 0, ay(x, 0± ) = -a, \x\ < c ax(x,y)->0, Txy(x,y) -> 0, ay(x,y) -> 0, r -> oo. (1.12) (1.13) This problem can in turn be generalized to the case in which the internal pressure opening the crack varies along the crack. It has been accepted that two dimensional crack problems can be categorized into three fundamental ones corresponding to three distinct modes of displacement: opening mode (Mode I), shearing mode (Mode II) and tearing mode (Mode III), as shown in Fig. 1.3 Fig. 1.4 and Fig. 1.5. It can be seen that any crack deformation can be represented by the appropriate superposition of these three cases. It also follows that for each of the three modes there is an associated crack tip stress field. Using the analytic method of Westergaard[29] Irwin [10] showed that near the straight crack tip, the primary stress components are of the forms h{6) + 0(1), Txy~^f2(8) + 0(1), oz Km MB)+ 0(1) (1.14) where 9 is the polar angle and r = \Jx2 j- y2. Since the leading terms contain a singular Chapter 1. Introduction Figure 1.3: Section of a plate with a Griffith crack: Mode I displacement Figure 1.4: Section of a plate with a Griffith crack: Mode II displacement 6 X Figure 1.5: Section of a plate with a Griffith crack: Mode III displacement term 1/y/r, the stresses are unbounded as r —> 0. This agrees with the results given by Inglis that at the crack tip, the tensile stress becomes infinitely large in the behavior of 0(ll^/r) as r approaches zero. An important feature of these asymptotic expressions is that the stress distributions around all cracks are similar and depend only on r and 9. The difference between one cracked component and another lies in the magnitudes of the constants K\, K\\ and Km, defined as stress intensity factors by Ki= Urn ay(a:,0)v^7rx (1.15) x—>0+ Kn= lim CTxy(a;,0)V^7r2 (1.16) #111= lim CTX2(z,0)v^7rx. (1.17) x—»0 + These stress intensity factors indicate the magnitudes of the crack tip stress field in the three modes. In each case, as the applied forces increase, the stress intensity factor wiD approach to its critical value, which indicates the material property, termed fracture toughness. Amoag these three modes, the most significant one is Mode I, since in most cases the stress behaviors are dominated by Mode I cases. Consequently methods of estimating K\ have received much Chapter 1. Introduction Chapter 1. Introduction 7 attention. In general stress intensity factors are determined by a number of factors. Previous research, mainly due to Irwin, shows that for a specific material, the stress intensity factor K (for each mode) can be expressed in general as K = f(°,c) (1.18) where o represents certain kind of measure of the applied load and c is the length of crack. The function / depends on the configuration of the cracked component and the way in which the loads are applied. Many functions have been determined for various specimen configurations and are available from fracture mechanics literatures. In this thesis, we will estimate K\ and K\\ of some specimens by using a relatively simple method which can also determine the stress field in the entire region of the cracked component. 1.3 Literature Rev iew In the early evolution of fracture mechanics, most methods used to estimate the stress intensity factors were based on either analytic or experimental approaches. In the analytic methods, a fundamental contribution was made by Westergaard in early 1930s. Sneddon [22] was the first to derive the stress field expansions for the two cases of a Griffith crack and a circular crack using integral transform methods. Williams [33] later extended the work of Sneddon and derived series solutions to the Inharmonic problem of the first type in an infinite domain with a reentered corner. Since 1960s, many important and valuable analytic solutions have been collected and given by Sih [19]. Most of them are based on the complex function method. The application of complex analysis is powerful in solving two dimensional linear elasticity problems. The approach using complex functions and conformal mapping was used in handling various simple geometric singularities with great successes. For the problems with simple geometry, such as circular disc, quasi-analytic solutions are possible(See, for example, [19|). Unfortunately when dealing with the problems with more general and complicated geometry and singularities, this Chapter 1. Introduction 8 approach faces severe challenges. Also, for a given problem, the transformed problem is usually not simple either and special efforts have to be made. Based on these disadvantages, today this approach is not used as widely as before and replaced by other general-purposed and efficient numerical methods such as the finite element and finite difference methods. Another popular approach was that using integral transforms. Some interesting, but not always of practical importance, problems such as the circular disc problems were studied by Rooke and Tweed[18] using Mellin transforms. Other theoretical discussions are seen, for example, in a systematic introduction to the application of integral transforms in elasticity theory and fracture mechanics by Sneddon [23]. The variations of the integral transform method have led to a number of simple approximations for estimating stress intensity factors for different configurations and loadings. In general these approximate methods are efficient sometimes, but they do not give detailed stress distribution. The integral transform method itself often suffers from the disadvantage that the inverse transform is required and is generally very complicated. Except for limited cases, analytic results are rarely available and numerical approximation is required. The J-integral method is frequently used to estimate the stress intensity factors. It is shown by Rice [16], [15] that the stress intensity factor can be related to a path independent integral. But in order to evaluate the J-integral, the displacement field must be determined first. Using the Fourier expansion approach, Vaughan and Wu [28] proposed a few solutions of circular disc problems in Mode I cases and some other problems. They used a specially selected series of Airy functions. On the boundary of the disc, applied loads were resolved as radial traction and tangential shear forces which were further expanded in Fourier cosine and sine series in polar angle respectively. The radial and shear stresses inside the disc were also expanded into Fourier cosine and sine series. An infinite linear system in the unknown coefficients of the Airy function was setup by equating the radial and shear stresses to the boundary conditions and then solved by considering only the leading diagonal subsystem. However, the Fourier expansion method can not be easily applied to more general cases. Chapter 1. Introduction 9 Either analytically or numerically solving the problem with Airy stress function in elasticity problems requires the solution of a biharmonic problem of the first type. Solving biharmonic problems using finite difference methods was studied in detail by Bjorstad [3] in early 1980s. The accuracies and complexities of various difference schemes were discussed in his work. Applying the usual 13-point scheme or other alternatives (see for example [3]) to the crack problem, however, faces the difficulty in obtaining a higher order approximation near the crack tip. Greenspan and Schultz [8] used the approach proposed by Smith [20], [21] (and studied in more detail later by Ehrlich [4], [5], [6], McLaurin [13] and Vajtersic [26] [27], etc.) to a rectangle plate with an edge crack under uniform tension. This approach separates the biharmonic problem into two coupled Poisson's equations and appropriate boundary conditions. Greenspan and Schultz used nine-point stencil near the crack tip, which generated a local truncation error of order 0(h6)) h being the local mesh size, to the Laplacian operator, and used uniform five-point stencil in the rest of the region. By using this technique, the user benefits from reducing the fourth order biharmonic problem to the second order Poisson type problems which are relatively easier to solve. The finite elements method is probably the most common method used today. In engi-neering, the elastostatic problems are usually solved through displacement models. Geometric singularities are usually handled by local refinement and/or the use of special elements. Nowa-days many commercial and research packages such as ANSYS have the facility to evaluate stress intensity factors. Although it has become a routine procedure that displacement field is solved first, it is not necessary that the problem be solved following this pat tern. Theoretical work was done investigating the finite element solution of biharmonic problems by Banerjee and Benson [1]. In the finite element methods that start with stress functions, the stress function is dealt with directly, and hopefully, the stress field and in turn the stress intensity factor can be obtained more accurately. Successful use of the finite element method to solve crack problems requires the user's knowledge and experience in selecting elements and generating a grid in the crack tip area. In the treatment of crack tip singularities, local refinement is usually used. Chapter 1. Introduction 10 However, as pointed out by Strang and Fix [25], in some cases where geometric singularities exist, the sole refinement may not be good enough and results may still be poor. In such cases special elements or basis functions are strongly recommended. Being aware of the deficiency of pure refinement techniques, some people turn to looking for special functions that approximate the solution in the vicinity of crack tip instead of numerical discretizations. This leads to the development of various "subtract-singularity-out" techniques. The achievement was first made in the studies of the second order elliptic problem of the form V • (pVu) -f qu — f in it, u — g on fi by Motz [14]. He proposed an approach in 1940s using a singular function to approximate the solution in the neighbourhood of boundary singular point while using usual finite difference scheme everywhere else. Woods [36] used the same idea and extended it to the problems with more general types of singularities. Fundamental contributions to the analytic solution of the first biharmonic problem with geometric singularities, e.g. a crack, was first made by Williams [33]. He used an infinite series to approximate the solution of the biharmonic function in the neighbourhood of the crack tip. Important work was also done by Whiteman [30], [31], [32]. Following Motz's idea, Whiteman [30] solved the problem in the rectangle plate with an edge crack by using finite difference approach outside a subregion containing the crack tip and using Williams' series in the subregion of the crack tip. 1.4 T h e Object ives In this thesis we will apply a semi-analytic technique that uses special Airy functions in the form of infinite series to the solution of stress fields and corresponding stress intensity factors to some plane crack problems. Particularly we are interested in the stresses and the stress intensity factors for a rectangle plate with an edge crack under different loadings which are to be used as specimens in laboratory tests. Also we will present some theoretical results in solving the fist biharmonic problem in the Airy stress function in the domain with a V-shaped crack and discuss the numerical approaches Chapter 1. Introduction 11 that are to be used. Furthermore, we will discuss some difficulties that arise from the use of our technique, and shed new understanding on why this type of difficulty occurs. 1.5 Scope of T h e Thes i s The main aim of this thesis is to derive a semi-analytical approach to the solution of the fist biharmonic problem in a domain with a straight reentrant corner or edge crack and present an applicable and relatively efficient numerical technique for practical uses. In this chapter fundamental concepts related to fracture mechanics have already been sum-marized. A literature review has also been given. Important work done by individuals will be referred to again throughout the following chapters. Chapter 2 presents the mathematical foundation upon which the primary method we used is based. The general solutions to the biharmonic equation in various forms in an infinite domain containing a "branch cut" or the line crack are obtained through the technique of separation of variables. Special Airy stress functions in the form of infinite series are developed. These expressions are extensively used in later discussions. A brief review of the use of Fourier expansion method employing these basis functions is given and followed by discussions. Chapter 3 discusses the property of those stress solutions and the numerical aspects in the use of these stress functions. Error measurements used in our numerical computations are introduced. Chapter 4 contains numerical results to some problems of particular interests in practice. The values of Mode I and Mode II stress intensity factors for two kinds of geometries - circular disc with a radial crack and rectangular plate with an edge crack, under various loadings, are given and followed by remarks and discussions. Finally Chapter 5 will summerize the achievements we have made in applying our approach to the assessment of stress intensity factors and in the meantime the determination of stress fields. Some open questions will be brought to the readers for further research. In this thesis, we try to keep the mathematical formulae as simple and general as possible. Chapter 1. Introduction 12 Detailed expressions will be put in the appendices for curious readers. Unless specified oth-erwise, all the variables such as dimension, stresses and loads are nondimensionali^ed for the simplicity of mathematical expressions. The stress expressions in terms of the special stress functions for Mode I and Mode II cases in polar coordinates are given in Appendix A. An alternative form of boundary conditions for the biharmonic problem of the first type is de-rived in Appendix B. Corresponding examples are given for the crack problems in rectangular plate and circular disc. Appendix C discusses the auxiliary problem in developing the stress functions that satisfy the crack surface stress free conditions on the boundary of a F-shaped angular corner. Appendix D introduces an approach to the solution of the first biharmonic problem as coupled Poisson's equations with appropriate boundary conditions, since Poisson's equation or Laplace equation in a cracked domain is relatively easy to solve and (local) analytic solutions are available. Appendix E discusses the numerical aspects in solving the linear least squares problem arising from boundary collocation with special stress functions in solving crack problems. C h a p t e r 2 Mathemat ica l Foundat ions In this chapter, we will discuss the solution of the first biharmonic problem arising from two dimensional elastostatic plane stress-strain problems in an infinite domain with an edge crack. Section 2.1 introduces the construction of the solutions of the biharmonic equation using specially selected harmonic functions. Section 2.2 discusses the Airy stress functions in an infinite domain with a F-shaped boundary, on which stresses vanish. Section 2.3 and 2.4 discuss the Fourier expansion method in solving crack problem using the special Airy stress functions. 2.1 Series Solut ion of Biharmonic Equation In polar coordinates, the periodic solutions of Laplace's equation V 2 u — 0 with real separable variables belong to the family H = {log r, r m cos m0, r m cos md,rn sin n0,r n sin n0} (2.19) where m and n are arbitrary real numbers. The generic solutions to the biharmonic equation V 2 V 2 u( r ,0 ) = O (2.20) can be formed using these harmonic functions. We first present the following results in the form of lemmas. L e m m a 2.1 For any harmonic function u, which satisfies V u — 0 and is sufficiently differ-entiate, the functions Vi = r2u, v2 = rcos&u and v3 = rsmOu satisfy the biharmonic equation. 13 Chapter 2. Mathematical Foundations 14 This can be verified directly. Lemma 2.1 points out a way of constructing generic solutions to the biharmonic equation from those to the Laplace equation. With the solutions in the forms given in (2.19), we claim L e m m a 2.2 The functions of the forms u = rmcos(n0) or u = 7r"sin(7i0) are biharmonic if \m\ — \n\ or \m ~ 2| — \n\, for m,n 6 R. Proof. The proof uses the results of (2.19) and the trigonometric identities. From Lemma 2.1 and linearity, we have more generally L e m m a 2.3 If u — u(r,6), v — v(r,0) are the harmonic functions, then the function w — (r2 — c2)u -+ v satisfies the biharmonic equation with c some constant. 2.2 Stress Funct ions for T h e Prob lems wi th F - S h a p e d B o u n d a r y Consider an infinite domain with a V-shaped boundary, that is, an angular corner with a < 2TT (see Figure. (2.6)). We now seek a particular type of stress function so that the normal and shear stresses given by this stress function vanish on the crack surfaces. Let us first look at the following simple conclusions. Propos i t ion 2.4 For any biharmonic function (f> = ^>(r,0), it is an Airy stress function such that on the reentrant boundary dilc on which 0 — ± a / 2 , the derived normal stress og — 0 and shear stress TTQ = 0 if and only if the function cf) and its normal derivative <pn are constants on dtlc. Proof. The proof is straightforward. Following is only for the sufficient condition. Since <f> — constant for all r, dd> d2d> S = « ? s 0 - (221) So dr ffe^-K^ = 0. (2.22) Chapter 2. Mathematical Foundations 15 Also, on dilc, Thus Therefore Figure 2.6: Infinite domain with a F-shaped angular corner. d6 dn d2(f> = drdO = ° ' d_ (\ dcj>\ 1 00 1 d2<p dr rTe \rdO) r drdd 0. • (2.23) (2.24) (2.25) This in fact has pointed out a way to construct special Airy functions that satisfy the stress-free conditions on straight crack surfaces. In Mode I cases, the stress distribution is symmetric, and therefore, by Lemma 2.2, let < h M ) - r A [ C l c o s ( A - 2 ) 0 + C2 cos X0] (2.26) where Cj, C2 are arbitrary constants, and A is chosen such that <f> and its normal derivative vanish at ± a / 2 . In the recent work by Vaughan and Wu [28] where a = 27r, A is assigned Chapter 2. Mathematical Foundations 16 two sets of values: Ai = | + k, k — 1,2,. . . , and A2 = 2 , 3 , . . . , and correspondingly, the stress function are represented by two infinite series <£(r, 0) = J2 AmrTn \m cos(m - 2)0 - (m - 2) cos(mfl)] m + X]^ n r n i c o 8 ( n ~ 2 )^ _ c o B n < ? ] (2-27) n with m = 3/2, 5 / 2 , . . . and n — 2 , 3 , . . . For the Mode I stress function (2.26), we may extend the above results to more general cases when a < 2TC. Propos i t ion 2.5 For an infinite domain with an angular corner or a V-shaped crack with angle a < 27r, the stress function (f>\ — <j>(r, 0; A) for Mode I of the form given by (2.26) will satisfy the crack surface stress free condition if A is chosen as the solution of the following eigen equation sin(A !)<.* - (A l ) s i n a . (2.28) with A / 1. Proof. Let <fi(r, 9) = rXf{0) = rx \cx cos(A - 2)0 + c2 cos X0\ (2.29) By Proposition 2.4, we are determining f(0) such that * ( r , ± f ) - & ( r , ± | ) - 0 or / ( ± | ) ^ / ' ( ± | ) _ 0, (2.30) That is C] cos(A - 2 ) - + c 2 c o s A - = 0, (2.31) Cl(A - 2) sin(A - 2 ) | + c2A sin A | = 0. (2.32) We now want to find A so that there are nontrivial solutions for C] and c2- This requires that the determinant vanish Acos(A- 2 ) - sin A - - (A - 2) cos A - sin (A - 2 ) - - 0 (2.33) Chapter 2. Mathematical Foundations 17 which can be simplified as sin(A - l ) a = - (A - l ) s i n a . D (2.34) When a < 2-K the solution of A (A ^ 1) has to be found numerically and may not always be real. Here we are interested only in the cases when a — 2ir, corresponding to an edge line crack. In this case, A has simple solutions. From (2.28), we have s i n ( A - l ) a = 0. (2.35) and therefore, A has the following values A = 1 + - , k = ±l,±2,... (2.36) In order that dominant terms in the stresses be of order O^r"1'2) as r approaches to zero, we take k positive. Thus, we have x = i2' i3; l> ••• (237) Having had the solutions for A, we now turn to determining the corresponding values for Cj and c2 in (2.26). Substituting A into the homogeneous boundary condition (2.31) for / ( # ) , we have, with a — 2ir, cos(A - 2)TT = COSATT = 0, for A = - + jfe, A: ^ 1 , 2 , . . . (2.38) and correspondingly, (2.32) becomes Cj(A 2) + c2A = 0. (2.39) Without loss of generality, we may define Cl = A, (2.40) c2 = 2 - A (2.41) and write ^AX) = Yl a*rX lA c o s ( A " 2)9 ~ (A - 2 ) c o s(A 60J • ( 2 - 4 2 ) A Chapter 2. Mathematical Foundations 18 where the coefficients a\ are undetermined. For A — 2, 3 , . . . , we have from (2.31) ci + c2 = 0 (2.43) and the corresponding expression for <fi takes the form <A2) -•=• Yl b ^ tcos(A " 2)° - ™s{\0)\ . (2.44) X Corresponding to each A there is a solution for cf>, thus, we obtain the stress function in series expression (2.26) for Mode 1 that satisfies the stress free conditions on the crack surfaces. For the anti-symmetric cases, e.g. Mode II cases, by Lemma 2.2, let the stress function <t>\{r, 0) = rx [cj sin(A - 2)0 + c2 sin \8] (2.45) Similar result holds for Mode II cases. Without proof, we give the following proposition. P r o p o s i t i o n 2.6 For an infinite domain with an angular corner or a V-shape crack with angle a < 2ir, the stress function <fi\ — <j>(r,0; A) for Mode II of (2-45) will satisfy the crack surface stress free condition if A is chosen as the solution of the following eigenvalue problem s in(A- l)a = ( A - l ) s i n a . (2.46) with A ^ 1. The topic on the solution of (2.28) and (2.46) is beyond the scope of this thesis. However, it will be interesting to know for what a, the corresponding stresses become unbounded as r —* 0. It is found that when a > 180°, Re(X) < 2 for Mode I cases, and consequently the stresses are unbounded. For Mode II cases, the stresses become unbounded if a > 254°. Numerical solutions and discussions are given in Appendix C. Similar arguments apply to higher order problems with boundary singularities More gener-ally, the singularity in a 2mth-order problem is determined by the principal part of the operator. The leading term in the solution is of the form rxip(0; A), where tp is a smooth function of 0 such as harmonic functions above and A is an eigenvalue of an associated problem, both of which depend on the boundary conditions. Chapter 2. Mathematical Foundations 19 2.3 T h e F o u r i e r E x p a n s i o n M e t h o d The first biharmonic problem can be stated as the following general form V 4 i i ( r ,0) = O infi Fu - f on dn, (2.47) Gu — g on dQ, where F and G represent, for example, the second order Unear operators that give the normal and shear stresses in terms of u. Obviously, the solution for u is not unique. However, we are interested only in the unique solutions of the second derivatives of u, which determine stresses. Assume the solution has the form °° u u(r> e) = E 9mrm<(>m{0), m - 1 + - , A = 1, 2 , . . . (2.48) m with 4>m(0) defined earlier. The coefficients qrn are to be determined from boundary conditions. In the recent work of Vaughan and Wu [28j, two boundary functions, namely the normal and shearing forces, are expanded into Fourier cosine and sine series / - E ^ c o s ^ (2-49) 11=0 OO 9 - E ^ s i n t ^ (2.50) Apply the operators F and G to (2.48) term by term and then expand (with respect to 0) in Fourier cosine and sine series respectively to obtain OO OO OO / OO where FU ~ E Q™ E fl"^(r) c o s /^ - E I E *n°W(r) J COS fid m fi=0 n=0 \ m ' OO OO OO / OO \ Gu = ] T qm ^ bmu(r) Sin i/^ = E ( E 9mbmv(r) J s i n vQ m u=0 v=0 \ m / amM = - I" F{rm{0)cf>m(0))cosf,0d0) /x = 0 , l , . 5T Jo bmu = - [^ G (rm(0)<t>m(0)) sin vOM, v = 1, 2 , . . 7T JQ Chapter 2. Mathematical Foundations 20 For circular disc with a radial crack whose length equals to radius a, r(0) — a and amfi and bmu are easy to calculate. Applying the boundary conditions leads to two infinite sets of linear equations oo 5^9mamM = //i (2.51) m oo ^ qmbmv = 9u (2.52) m for fi,v — 1,2, . . . , which are to be solved simultaneously. Assume that the series (2.48) converges for all r and 6 (This is not guaranteed in general) and decreases fast enough, we use only the first finite terms. Take fi,L> — J, • • •, n and m = 1, • • •, In. This gives a 2n x 2n linear system Ax = b where AtJ = atJ or b,j (2.53) x = (ax, • • •, aM, 6i, • • •, b2n-M)T, ( 2 - 5 4 ) b=- (h,-,fn,9u-Jn)T- (2.55) The coefficient matrix A is nonsingular and there exists a unique solution for x. 2.4 Discuss ion on T h e Fourier Expans ion M e t h o d One advantage of the Fourier expansion method is that the singularities of the second type uii-smooth da ta in the boundary conditions, for instance, compact loadings, can be handled better through Fourier expansions. Such unsmoothness is essentially "averaged out". This technique has been used with great success by Vaughan and Wu in their recent work [28] in solving the problem of a cracked disc under "concentrated loadings". Generally speaking, the success of the Fourier expansion method strongly depends both on the shape of the boundary and the Chapter 2. Mathematical Foundations 21 boundary values. Recall that when solving the problems in which variables are separable, the solution u(x,y) is expanded in a series of eigenfunctions u(x,y) - ^TcTlNTl(x)Mn(y), x,y G 12 n where Mn(y) are the eigenfunctions (sine or cosine functions), Nn(x) are known functions and cn are to be determined. Given boundary function / (e.g. stress function or its normal derivative or stresses), the problem is reduced to finding c n in the boundary condition Bu = / , x,y 6 dil, where B is a linear operator, e.g. identity or partial differentiate operator. Expand both sides on the boundary into Fourier series to obtain following problem \ n m I m where Hm are sine or cosine functions of the form u / \ • m 7 r u i \ m 7 r Bm(x) = sm- x, or HTn(x) = cos- x b - a b - a and Fm are the Fourier coefficients defined by rb Fm{y)- j v(x)Hrn(x)dx, m = 0 , l , . . . Ja for v(x) defined over [a,6](e.g. [0,7r]). Typical regions on which this method works well are usually limited to rectangle and circular disc when Cartesian and polar coordinates are used. In these cases Nn become constant on the boundary. The evaluation of the Fourier coefficients on the left-hand side, however, is not easy in general for complicated functions. For example, for a disc with a radial crack of length less than radius a(as shown in Fig. (2.7)), for which the origin is off the center of the circle by e, the equation of the edge is given in polar coordinates by E(r, 0) = r 2 2re cos 6 + e2 - a2 = 0 or 1 / 2 , r = R(0) = \e cos 0 ± (a2 - e2 sin (?) Chapter 2. Mathematical Foundations 22 Figure 2.7: Circular disc with a radial crack less than the radius, coordinates is located at the crack tip. The origin of the polar Figure 2.8: Rectangle with an edge crack. Chapter 2. Mathematical Foundations 23 For the rectangle with an edge crack, as shown in Fig. (2.8), for the point at the top edge, r(6) — hj sin#. Both cases make the evaluation of the Fourier coefficients Fm fairly complicated. Therefore in solving the problem with irregular boundaries more general methods such as finite element and finite difference or other alternatives will be more advisable. Chapter 3 Numer ica l Approaches and Analys is This chapter discusses in a general way how and how well the solution of a crack problem using the special stress functions derived in the previous chapter can be obtained numerically. 3.1 T h e Localness of T h e Basic Solut ions In the previous chapter, series solutions (2.27) that satisfy the homogeneous boundary condi-tions on the F-shaped reentered boundary in an infinite domain were derived, in the case of a finite region, these basis functions generally will not satisfy the boundary conditions specified on the boundary that encloses the region. This suggests that for a crack problem in a finite region, these basic solutions are locally defined. It is therefore interesting to know when these solutions may be applied globally and when they may not. Note that in the neighbourhood of crack tip O, the series solution is absolutely convergent when r < TQ for some fixed r0. However, the convergence of the series is not assured globally. A frequently used approach in the numerical solution of crack problems is to use singular solutions such as (2.27) only locally in the neighbourhood of the crack tip and apply conventional t reatments , e.g. the finite element or finite difference methods elsewhere. In the finite element methods, the use of special basis functions and the alternatives such as (2.27) has led to the development of the so-called singular elements. In their theoretical work on the local treatment of boundary singularities in the finite element methods, Strang and Fix [25] suggested that in the vicinity of a singular point, for example, the crack tip, singular basis functions in the form of power series or logarithms are preferred. Furthermore, they showed that there exists an approximation in the spanned finite element 24 Chapter 3. Numerical Approaches and Analysis 25 space consisting of the usual piecewise polynomials together with some specially selected sin-gular basis functions, such that the error in the approximate solution is bounded. In the finite difference methods, there corresponds to the "subtract-singularity-out" schemes, typically, the Motz's and Woods' methods. In these methods, the singular basis functions are applied only at the grid points in a small neighbourhood of the crack tip and classic finite difference schemes are applied at the grid points outside the neighbourhood. It has been well accepted that the finite difference or finite element (without "singular elements") method alone will lead to a poor approximation to the true solution unless very fine local i J inement is performed. 3.2 N u m e r i c a l A p p r o a c h e s When solving crack problems in a domain in which one dimension is much bigger than another, such as a long strip, the boundary collocation method introduces further numerical difficulties. The major problem is that the coefficient matrix may become ill-conditioned due to the use of massive collocation points and large amount of terms (of power series) required by the accuracy concern. This is because, on the one hand, the elements in the A;th column in the matrix, which correspond to the fcth term in the series, become extremely small (or extremely large if r is not scaled so that r < 1) as k increases, so that the coefficient matrix becomes nearly not full rank due to the non appreciable roundoff effects. On the other hand, the radius of two distinct collocation points that are far from the origin are so close to each other(as in the long strip plate case), that the two corresponding rows in the matrix become nearly identical. In such cases one may apply the series only locally in the small neighbourhood dflo °f the singular point or crack tip and use other form of basis functions smoothly connected to the truncated series on the boundary dilo °f fio a n d then collocate at globally selected points. One can also use the finite element approach. For the finite element solution of the biharmonic problems, we refer the reader to Banerjee and Benson's work [lj. Another choice is to use the finite difference approach with local analytic solutions. Most of the methods use the approach of "subtract-singularity-out", originally due to Motz' and Woods' Chapter 3. Numerical Approaches and Analysis & 1 < 1 ( (i »-Figure 3.9: Motz' finite different mesh near the crack tip 0. with • representing near-points and o representing far-point. The crack opening is towards left. Grid points are not specified along the crack. work (see [14] and [36]). In solving the second order problem of the form V • (pS/u) f qu — / in fi, u ~ g on fi in a region with an edge crack, Motz divided the region !f2 into two parts: a neighbourhood fio and the whole region except Q.o, that is, fi\fio, a n d used a set of mesh points. The mesh points in ilo a r e called near points, while an equal number of adjacent boundary points in il\ilo a r e called far points (as shown in Fig. 3.9). Using an infinite series oo oo u = J2qxux(r,6) - ]TgAT-A/(0; A), A A which captures the singular behavior of the solution near the crack tip, in truncated form at both near- and far-points, Motz obtained an expression which maps the approximate values at far-points to those at near-pints, that is, un - Muf, where superscripts n and / refer to near and far points and M is the matrix representation of the mapping. By setting up the finite difference equations at mesh points in fi\fio) Motz had Chapter 3. Numerical Approaches and Analysis 27 a set of linear equations which gives u at mesh points and un via the mapping. The coefficients qx can then be solved simultaneously. In Woods' method, which was originally applied to the BVPs with Laplace (or Poisson's) equation, the harmonic function u is split so that u = (f>(a,(3) + if) where <p is a harmonic function with two parameters having the desired form of the singularity at the singular point 0. Woods assumed that a suitable choice of a and /3 will make </> capture the entire singularity, so that ip is well behaved. There is no difficulty in applying the Motz's or Woods' method to biharmonic problems in a cracked plated. For more details, the reader may look at Bernal and Whiteman's work [2] in which both the Motz' and Woods' methods were applied to a rectangular region with an edge crack. 3.3 E r r o r M e a s u r e m e n t s In the boundary collocation method, for any function w evaluated at a set of boundary points (£m,Vrn.), m — 1 , . . . , M, introduce the discrete Li norm notation / l M \ V 2 H I = ( M E K » I 2 J • (3-56) Define 11*11! = (N|2 + K| |2)1 / 2 (3.57) for u on the boundary, where un is the normal derivative of u outwards to the boundary. Also define as maximum norm Halloo = max( |wm | ) (3.58) m for w on the boundary. Similar notation is introduced for u and u7l as IMU.i = maxdluHoo + Halloo) (3.59) m Now for the boundary conditions u ~ / , un ~ g on the boundary, define compound errors Chapter 3. Numerical Approaches and Analysis 28 measured in II • IIi and « i [jj ( E I"™ - f(U,fhn)\2 + E K«»)™ - $(6»,»/m)|2] j (3.60) ||e||oo = m a x ( | u m •- / (£,„, r)m)\, \(un)m - s(£m,77m)|) (3.61) To measure the error in stresses in the collocation procedure, define IN* = ( ^ ( E K» - a(^> '/-)i2 + E iT- - r(£-> ^) i 2 ) ) (3-62) ||e||oo = niax(|CTTO - ^(Cm,Vm)\,\TTn - T{(m,rjm)\) (3.63) where a and r are given on the boundary and am and r m are the normal and shear stresses evaluated at the mth point ( ^ m , ^ m ) , respectively. If write (3.62) and (3.63) in terms of stress function (j>, then these two are in fact semi-norms for </>. Equation (3.60) through (3.63) are used to measure errors in our numerical experiments. Equation (3.60) and (3.61) are convenient for mathematical analysis while (3.62) and (3.63) are of more practical meaning. By the definition (3.60), it is obvious that the error as measured in the £2 norm vanishes if and only if the approximate solution u and un are identical with the given function / and g respectively at each boundary point. However, in most cases, we required that ||e||i be minimized rather than ||e||x = 0. 3.4 B o u n d a r y Col locat ion M e t h o d Suppose that dimensions of the region under concern are comparable in two directions (x and y), that is, the region is not strip like, we may assume that the solution is admissible in the whole region including the boundary (See the discussions in the next section). An ideal case is that the region is a circular disc with a radial crack. In particular when the radial crack length is equal to the radius of the disc, the stress function depends only on the polar angle 6 on the edge of the disc. The problem can then be solved in a nice way using the Fourier expansion method (See Vaughan and Wu[28j). Chapter 3. Numerical Approaches and Analysis 29 In more general cases, where the region has an irregular shape, numerical methods are preferred. Numerically, we have to use the series in the truncated forms, containing finite terms so that the effect of round-off(which is machine dependent) will not be too significant. Then we choose a reasonable number of collocation points, say, M / 2 points, on the boundary, at which two boundary conditions are forced to be satisfied. This leads to an M X N linear system from which A\ can be solved. Since the collocation points are chosen only on the boundary, it is often referred to as the boundary collocation method. When M = N, the linear system is square and solved exactly. At each collocation point, the boundary conditions, namely the normal stress and shear stress or their alternative forms, are satisfied exactly. However this may not be always good in practice. First of all, although the boundary conditions are satisfied exactly at selected points, the boundary conditions may be seriously violated everywhere else due to unsufficient number of collocation points. This is especially the case when data changes significantly along the boundary. Secondly, if one includes more terms with a corresponding increase in the collocation points, the coefficient matrix, which is formed by the entries defined by the values of r f(0; A) and/or its derivatives, is found to contain many very small numbers and is effectively ill-conditioned. An inversion of the matrix is then not possible. Instead of setting equal numbers of terms and boundary conditions, we choose more collocation points and in turn more equations than the number of coefficients and solve the linear system using the least squares method(See discussions in Appendix E). We are now turning to the procedure of boundary collocation. With the truncated series expression of the stress function N 0(r,0)-£>ArA/(*;A) A we choose a certain number of collocation points from which we obtain M equations. This leads to an M X N linear system. If M — N, then ||e||i = 0. By adding more collocation points, one can expect that the error ||e||i may increase. We are looking for such an N, if it exists, tha t ||e||i is bounded as the number of collocation equations increases. The problem can then be phrased as follows: Chapter 3. Numerical Approaches and Analysis 30 In the stress function of the form, 4> = ct>{r,d;P)^Y.1,rMr,0) = Y.^Tmfn{9), r,0eSl, (3.64) n n find p = (pi, • • • ,PN) with N to be determined, such that for the given conditions (constraints) defined on dil < M / M ) , ^=g(r,0), r,9 e dil (3.65) on and given constant e, the error defined at the total M collocation points / , / M Up,, \\X^2 IMIi= I j i j lE 1^ " / ( '» . M 2 + EK^)»-«( '» . M'JJ J <e. (3-66) We perform the following procedures. Initially, truncate the series into Ni terms and select boundary points to get M\ collocation equations. Gradually adding more points and in turn setting up a sequence of sets of collocation equations (sequentially of M x , . . ., MR equations), we obtain a sequence of solutions to the vector of coefficients {p}k — {-^Xi," •' > A\N }k}> k = 1 , . . . , K. Note that only the leading term A:i/2 is of significance since it is directly related to the stress intensity factors. With this in mind, we perform the procedure until a stable A3/2 is obtained, denoted by ^3/2,^1 • Then we increase the number of terms to N2 and perform the same procedure to get A 3 / 2 I JV 2 . Repeat the process until a steady state in {A3/2n}, n = l,...,Ni is achieved. The last number Ni, of truncated terms can then be considered the optimal one. As mentioned earlier, in order to avoid the round-off effect, the number of terms Nk, k — 1 ,2 , . . . should not be taken to large. By solving the problem in the least squares sense the error defined by (3.60) is measured in global sense. The approximate (or the least squares) solution on the boundary apparently does not match the given conditions exactly but only compromises the constraints in the least squares sense. Since we are primarily interested in the stress field near the crack tip, if small changes in boundary conditions causes only small changes in the leading coefficients, in particular A3/2, even though the boundary conditions are not satisfied perfectly in some local portion of the boundary, the approximate solution shall be considered as the best fit to the problem. One may also judge the least squares solution by examining the resultant force on the boundary. Chapter 4 N umerical Exper iments This chapter presents numerical experiments on several model problems, which represent the physical problems of using laboratory specimens to test for the the fracture toughness by ap-plying different loadings. The major work of this chpater, and also the main objective of this thesis, is to calculate stress intensity factors for the test specimens with the application of different loadings. As mentioned in Chapter 1, the stress intensity factor in general is a function of load and geometry, that is K :• f(0,C) where a represents certain kind of applied force and c is the crack length. The function / depends on the configuration of the cracked component and the way in which the loads are applied. For a rectangular plate with an edge crack subjected to certain (distributed) loadings, the stress intensity factor is written as K = Kuf(c/a,h/a,l/a) (4.67) where KQ is some reference quantity with the dimension of stress intensity factor and / is a dimensionless function. Here variables a is the width of the rectangle, c is the crack length and I is the length over which the load is applied. The emphasis of the numerical work in this chapter is to find function / numerically for various cases. The primary geometry of the testing specimens concerned in this chapter will be a rectangle with an edge crack dividing the rectangle into two symmetric parts. The loadings include uniform tension, concentrated loads and those yielding pure bending. In our model problems, 31 Chapter 4. Numerical Experiments 32 we are primarily concerned with the calculation of Mode I stress intensity factors with the variation of geometric parameters and loading manners. Numerical experiments also include the a t tempt to Mode II cases. In addition to stress intensity factors, stress distribution and error estimations are also presented. In all of the model problems we solve for the stress function using the boundary collocation method. The linear systems arising from collocation were solved in the least squares sense with the approach of orthogonal transformations. The numerical experiments were carried out on Sun SPARC workstations with SunOS re-lease 4.1.2. All computations are in double precision. All of the graphics representing numerical results were generated with TecPlot. This chapter is organized as follows: Section 1 introduces the errors measurements for later use, and describes the conventions of nondimensional representation of quantities involved in all of the model problems. Section 2 presents the comparisons of our numerical results with those obtained in [37] using Fourier expansion method. Section 3 through Section 5 present our new Mode I results for the test specimens under different loadings. Finally, Section 6 contains a few experiments for Mode 11 loading with discussion. 4.1 N o t a t i o n and Convent ions Most of the numerical results presented in this chapter are dimensionless. Normally, we write boundary conditions in the following form a(s) = a0f(s) (4.68) where OQ is the reference quantity, called reference stress or stress scaling factor in this thesis, with unit of force/area, and f(s) is a dimensionless function. As introduced in Chapter 2, the special Airy stress function(with dimension) is written in general * = l > m V > m ( r , * ) (4.69) m Chapter 4. Numerical Experiments 33 where ipm(r,9) are dimensionless and <4*t have the dimension of force. The general form of boundary condition is B<i> - aQf. (4.70) where B is a linear operator, for example, B — d2/dx2. Define dimensionless stress function (f> by 0 - * M ) , (4.71) then, the series expression of our special stress function is written m - Y,Ami>(r,6) (4.72) m where Am are dimensionless coefficients to be determined. The dimensionless boundary condi-tion is then B4> = YiAmB*m = f- (4.73) m Correspondingly, the dimensionless stresses are denned by a - a/au (4.74) where a in general denotes normal and shear stresses. For instance, the boundary conditions for a circular disc under uniform tension are expressed in computation as aT — 1, TTI) = 0, at r -- a for 0 < 9 < 2ir. Without causing confusion, we also use the notation without bar to refer to dimensionless stresses. With dimensionless coefficients Am determined, as mentioned in Chapter 2, the stress in-tensity factor is then given by K = cA\ = cA*au (4.75) 2 2 Chapter 4, Numerical Experiments 34 where c = 3/2 for Mode I cases and c — 1 for Mode II cases. Define also dimensionless stress intensity factor K by K = K/K0 (4.76) where K0 is the reference stress intensity factor in terms of o0 and geometric parameters, for example, KQ — ooy/c, where c is the length of crack. For different problems, KQ will be defined differently. Without causing confusion, in some cases, we simply use A'[ or K\\ without bar to refer to dimensionless K\ and K\\, assuming the reader understands such simplification. To measure the errors in numerical discretizations, in the previous chapter we have intro-duced errors measured in the discrete L2 norm and L^ norm. Again define respectively for a total number of m computed values of stresses, the root-mean-squares error = [ i ( E K - *(-*)ia + E to - T('*)ia)) > (4-77) and maximum error ||£||oo = max|£;fc| m = mzx(\ak-a(sk)\}\Tk-T(sk)\) (4.78) where o(sk), T(sk) are the given normal and shear stresses at boundary point sk, respectively, and ak and rk are the calculated normal and shear stresses at that point. With dimensionless stresses, define corresponding dimensionless errors by 1/2 [i ( E I** ~ *('*)!' + E to - *('*)l2)) . (4-79) and max \ek\ m m*x{\ak-a(sk)\,\fk-f(Sk)\) (4.80) Chapter 4. Numerical Experiments 35 By the definitions, these are relative errors with respect to the reference stress or stress scaling factor OQ mentioned previously. 4 .2 C o m p a r i s o n of Numer ica l Resu l t s wi th Prev ious Work For comparison, we applied the boundary collocation approach to the circular disc problems studied by Vaughan and Wu[28]. Case 4.2 .1 Circular disc under uniform tension. As studied in [28], consider a circular disc of radius a with a radial crack of length c = a subjected to uniformly distributed stress er0 on the circumference, as shown in Figure 4.10. Using the stress function introduced in Chapter 2, the asymptotic stresses in front of the crack Figure 4.10: Circular disc under uniform tension cr0. 3 . fa\ll2 tip are 0Y = oe = - A i aQ I -2 2 \r The stress intensity factor is given by TTS — 0, at 0 = 0, as r - > 0. (4.81) 3 (a\ll2 , Ai — -AiOQ\ - \ V27rr ( - : ) ' — - \Z2TT A 3 ao\/a 2 5 3 r— --=• -y/2*AsKa 2 2 (4.82) Chapter 1 Numerical Experiments 36 3.185 3 155 20 40 60 80 100 120 Number of collocation points 140 160 180 Figure 4.11: The variation of Ki with the number of collocation points (circular disc, uniform tension). where KQ is the reference stress intensity factor. Because of the symmetry, the problem is solved in the upper half plane only. The dimen-sionless leading coefficient ^ 3 is found to be 1.495, which agrees closely with the value given in 2 [28] multiplied by \f^K. The first ten coefficients Am and Bn in the first and second series (in the notations used in Wu[37] and Vaughan and Wu[28]) are tabulated in Table 4.1 with the use of different number of truncated terms. The values of the first ten coefficients in forty truncated terms are compared with the results in [37] in Table 4.2 The variation of stress intensity factor with the number of collocation point using different number of truncated terms is investigated and shown in Figure 4.11. It indicates that , with enough terms, six or more in each series, the dimensionless stress intensity factor approaches to an asymptotic value as the number of collocation points increases. It is seen in Figure 4.11 that dimensionless K\ achieve the accu-racy of 0.0025 if more than forty points are used. Thus, the number 40 can be considered as a truncation value for the selection of collocation points. Figure 4.12 shows a quick convergence of K\ with respect to the number of truncated terms using fixed number of collocation points. Chapter 4. Numerical Experiments 37 M = Am+l 1.49576 -0.490013 -0.0067706 0.0015948 -0.0005349 0.0002234 -0.0001157 8.8686-05 -2.2202-05 -2.0686-05 = 20 % + i 0.474267 0.0890304 -0.0101775 0.0025807 -0.0008839 0.0003730 -0.0001732 -0.0003152 0.0003092 4.1862-05 M = A . i m-|- 5 1.49548 -0.489853 -0.0067869 0.0016016 -0.0005361 0.0002190 -0.0001020 5.2187-05 -2.8723-05 1.6776-05 = 40 Bm+l 0.474096 0.0890454 -0.0101814 0.0025658 -0.0008480 0.0003176 -0.0001237 4.5450-05 -1.1973-05 -2.6481-06 M = Am+k 1.49537 -0.489794 0.0067921 0.0016041 -0.0005375 0.0002198 -0.0001025 5.2468-05 -2.8832-05 1.6742-05 = 80 Bm + l 0.474036 0.0890502 -0.0101876 0.0025701 -0.0008510 0.0003196 -0.0001251 4.6356-05 -1.2776-05 -1.7590-06 Table 4 .1: Ten leading coefficients with using 20, 40 and 80 truncated terms (circular disc, uniform tension). Ge ^4T O + i Bm+1 1.49548 0.474096 -0.489853 0.0890454 -0.0067869 -0.0101814 0.0016016 0.0025658 -0.0005361 -0.0008480 0.0002190 0.0003176 -0.0001020 -0.0001237 5.2187-05 4.5450-05 -2.8723-05 -1.1973-05 1.6776-05 -2.6481-06 Wu Am+i_ Bm+l 1.49540 0.474066 -0.4897740 0.0889071 -0.0067858 -0.0101788 0.0016022 0.0002567 0.0005372 -0.0008482 0.0002199 0.0003173 -0.0001005 -0.0001257 5.0266-05 4.3982-05 -2.8274-05 -9.4248-06 1.570805 -3.1416-07 Table 4.2: Comparison of leading coefficients with Wu's results [37| (circular disc, uniform tension). Chapter 4. Numerical Experiments 38 318 6 8 10 12 Number of terms in the 1 st and 2nd series Figure 4.12: The variation of K\ with the number of truncated terms (circular disc, uniform tension). The errors in the boundary collocation are plotted in Figure 4.13 and 4.14. As one may expect, as the number of terms involved in the series increases, the errors, both measured in the root-mean-square and maximum senses, drop significantly. Figure 4.14 also indicates tha t , in order that the maximum error at certain particular point be bounded, 30 or more terms must be used. C a s e 4 .2 .2 Circular disc under diametric loading As shown in Figure 4.15, the concentrated load P is replaced by a sinusoidal distribution of radial tensile stress over an arc on the circumference of angle a — 7r/K, K — 4 / c f l , k 6 N^(positive integers). The reference stress or stress scaling factor in this case is defined as P 1 P oo aat a at (4.83) The stress intensity factor is given by Ki -sjl-KAzOnsfa 2 2 Chapter 4. Numerical Experiments 39 n <D (6 ? c (0 1 •5 0 IT 0010 0.009 0.008 0007 0 006 0.005 0.004 0.003 0.002 0.001 0 000 n r 10 terms 14 terms 20 terms 34 terms 40 terms _l_ _l_ 60 80 100 120 Number of collocation points 140 160 -< ~o 180 Figure 4.13: Root-mean-square error ||e|| with the number of collocation points (circular disc, uniform tension). 0 05 10 terms 14 terms 20 terms 30 terms 40 terms < 60 80 100 120 Number of collocation points 140 160 180 Figure 4.14: Maximum error HeU^ with the number of collocation points (circular disc, uniform tension). Chapter 4. Numerical Experiments 40 Figure 4.15: Circular disc under diametric loading. 3 / = - , 1 P /-= -V2itAz y/a 2 * aat O 1 -V^AS-KQ 2 a (4.84) where t is the thickness of the disc and KQ is the reference stress intensity factor. Unlike Case 4.2.1, we have found in our numerical experiments that the convergence rate of the series seems fairly slow. More than fifty terms in each series have to be used to achieve the convergent value of K\ at certain accuracy. The slow convergent behavior was also observed by Wu[37], in which Wu used Fourier expansion method to solve for those coefficients. The first ten coefficients of Am and Bn using different number of truncated terms for a — 7r/21 are tabulated in Table 4.3. Table 4.4 compares the results with those obtained by Wu[37]. In [37], the results are given in four decimal digits only. It can be seen that Am agree with Wu's results very well while Bn differ from Wu's in higher order terms. C a s e 4 .2 .3 Circular disc under crack opening loading As shown in Figure 4.16, the concentrated load is replaced by a sinusoidal distribution of vertical tensile stress over an arc of circumference of angle a = -K/K, K — Ak \- 1, k £ iV+ near the Chapter 4. Numerical Experiments 41 M = A , i 0.677417 -0.0712769 -0.0940771 -0.0626092 0.0489328 0.0402994 -0.033773 -0.0289136 0.0252785 0.0223942 = 100 Bm + l -0.17343 0.0510315 0.308686 0.002054 -0.313463 0.000165363 0.308099 -7.41027-05 -0.302041 -9.90934-05 M = Am + k 0.674891 -0.0697791 -0.0942235 -0.0625382 0.0488905 0.0403274 -0.0337928 -0.0288989 0.0252673 0.022403 = 140 Bm+i -0.175245 0.0512957 0.308401 0.00224908 -0.313636 0.000290014 0.307983 1.16497-05 -0.302124 -3.82565-05 M --A . • 0.673737 -0.069248 -0.0942277 -0.0624855 0.0488546 0.0403197 -0.033786 -0.0288833 0.0252549 0.0223993 = 200 Bm+i -0.175788 0.051464 0.308152 0.00236614 -0.313565 0.000368668 0.30782 7.01571-05 -0.302064 7.93027-06 Table 4.3: Ten leading coefficients with using 100, 140 and 200 truncated terms (circular disc, diametric loading). Ge Am+i_ Bm+X 0.67374 -0.17579 -0.069248 0.051464 -0.094228 0.30815 -0.062485 0.0023661 0.048855 -0.31357 0.04032 0.00036867 -0.033786 0.30782 -0.028883 7.0157e-05 0.025255 -0.30206 0.022399 7.9303e-06 Wu Am+i_ Bm M 0.672 -0.176 -0.069 0.051 -0.094 0.308 -0.063 0.002 0.049 -0.314 0.040 0.358 -0.034 0.309 -0.029 0.630 0.025 -0.304 0.023 0.265 Table 4.4: Comparison of leading coefficients with Wu's results [37| (circular disc, diametric loading). Chapter 4. Numerical Experiments 42 Figure 4.16: Circular disc loaded at the crack mouth. crack opening. While in [37] and [28], P is replaced with a sinusoidal distribution of shear stress plus a radial tensile stress, which constitute an equilibrium force system. The stress intensity factor is given as in case 4.2.2 by K\ =• -\/2TrAso0y/a 2 2 3 r— A 1 P n = -\2ltAs v a 2 2 a at 3 1 = -V2nA3-K0 (4.85) 2 2 a where t is the thickness of the disc and KQ is the reference stress intensity factor. The first ten coefficients of Am and Bn using different number of truncated terms for a — 7r/21 are tabulated in Table 4.3. Table 4.4 compares the results with those obtained by Wu[37]. In [37], the results are given in four decimal digits only. These results differ from Wu's slightly, but the magnitudes are in the same order. The requirement of a large number of terms in computations in Case 4.2.2 and Case 4.2.3 implies that the boundary collocation method becomes less efficient in these cases and other Chapter 4. Numerical Experiments 43 M = 4»+i 2.05663 -1.11571 0.0941759 -0.0459668 0.0274806 -0.0183698 0.0131969 -0.00997126 0.00781862 -0.006305 150 Bm+i 1.06033 -0.051336 0.0994753 -0.0809569 0.065869 -0.0548056 0.0463156 -0.0395754 0.0339305 -0.0291029 M = Anfl 2.05248 -1.11361 0.0940231 -0.0458911 0.0274353 -0.0183396 0.0131755 -0.00995536 0.00780643 -0.00629543 180 Bm + I 1.05884 -0.0514795 0.0994843 -0.0809178 0.0658648 -0.054781 0.0463167 -0.0395619 0.0339365 -0.029097 M = A , , 2.05301 -1.11384 0.0940315 -0.0458957 0.0274381 -0.0183415 0.0131768 -0.00995638 0.00780721 -0.00629604 200 Bm+i 1.05884 -0.0513412 0.099421 -0.0808637 0.0658285 -0.0547451 0.0462904 -0.0395344 0.0339155 -0.0290744 Table 4.5: Ten leading coefficients with using 150, 180 and 200 truncated terms (circular disc, crack opening loading). Ge Am+k 2.05301 -1.11384 0.0940315 -0.0458957 0.0274381 -0.0183415 0.0131768 -0.00995638 0.00780721 -0.00629604 #m+l 1.05884 -0.0513412 0.099421 -0.0808637 0.0658285 -0.0547451 0.0462904 -0.0395344 0.0339155 -0.0290744 Wu A^k 2.047 -1.111 0.094 -0.046 0.027 -0.018 0.013 -0.010 0.008 0.006 Bm^\ 1.057 -0.052 0.100 -0.082 0.067 -0.056 0.048 -0.042 0.038 -0.034 Table 4.6: Comparison of leading coefficients with Wti's results (37| (circular disc, crack opening loading). Chapter 4. Numerical Experiments 44 efficient approaches are demanded. 4.3 Rectangular P la te under Uniform Tension Consider now an ideal case, a rectangular plate with an edge crack subjected to uniform tensions <7o loaded on its two opposite sides, as shown in Figure 4.17. In terms of stress function, the <?o B t t t t M t t G O rrrrrrn D E Figure 4.17: Rectangular plate under uniform tension. boundary conditions are given by 32$ d2$ dy2 dxdy 0 on OA, AB and CD, d2* d2* OQ, -7.—— — 0 on BC . (4.86) (4.87) dx2 ~u ' dxdy Because of the symmetry, the problem is solved in the half region 0 ABC DO. Figure 4.18 shows the distribution of collocation points along the edges. Note that the selection of collocation points must exclude the points A and D because the stress function and its derivatives vanish at 6 — 0 and IT. Special treatments are needed at corners B and C because the normal derivative is not denned. In terms of stresses, one more condition in normal stress is added at each of the two corners. The boundary conditions specified by (p and its (outward) normal derivative (frn on the physical boundary 0 ABC DO are given in the following standard statement of the Chapter 4. Numerical FJxperiments 45 biharmonic problem of the first type (See Appendix B) V 4 # = 0 in O ABC DO $ = $ n = 0 on OA $ = 4>n •=- 0 on AB * = ^ i x +c)\ *n = 0 on BC 2 * = ~Y~(X •+ c)2> *»» = ^oa o" CD Use the truncated series introduced in Chapter 2 M cj)(r, 0) = ] T Amim [m cos(m - 2)6* - (m - 2) cos m0] m JV (4.88) + £ ] - B n f [cos(n - I ) * - c o s n0] (4.89) for rn. = 3/2, 5 /2 , . . . , M , n = 2 , . . . , N, 0 < 0 < TT, where <f> — $/CTQ and £ = r/tf, d is some scaling metric, for example, d = (a2 -f 62)1 '2 . Solving for (dimensionless) unknown coefficients Am and i?n using boundary collocation, the asymptotic stresses in front of crack tip are given by 3 / A 1 / 2 <7X = ay ~ -A3CT0 - , Txy = 0, at 0 = 0 as r -> 0. (4.90) 2 2 \ r / The stress intensity factor K\ is given by 3 / d \ 1 / 2 / K\ = - A 3 cr0 - V 27rr 2 2 \ r / 3. (dY/2 r~ (4.91) Figure 4.18: Distribution of collocation points along the edges of cracked rectangle Chapter 4. Numerical Experiments 46 15 30 45 60 75 90 105 120 135 Number of collocation points Figure 4.19: The variation of dimension less K\ with the number of collocation points (rectangle, uniform tension). Define the reference stress intensity factor K0 - aov^rc, (4.92) then Ki = - A j ( j j K0 (4.93) The dimensionless K\, defined as K\jKQ is calculated and presented in the following figures. The variation of dimensionless K\ with the number of collocation points for the case cja — h/a — 0.5 is given in Figure 4.19 for the use of different number of truncated terms. For each case of truncated terms, K\ is bounded as the number of collocation points increases. As the number of truncated terms increase, K\ converges to the value of 3.0. This behavior is shown in Figure 4.20 for fixed number of collocation points. The corresponding errors defined earlier for the above case are plotted with respect to the number of collocation points in Figure 4.21 and 4.22 respectively. Both errors are bounded as the number of collocation points increases. As we expect, the more truncated terms, the Chapter 4. Numerical Experiments 47 3.100 2.800 20 30 Number of terms 40 Figure 4.20: The variation of dimensionless K\ with the number of truncated terms (rectangle, uniform tension). smaller the errors. The first ten coefficients of totally 30, 40 and 50 truncated terms for c/a — h/a — 0.5 are tabulated in Table 4.7. Unlike what we have seen in the circular disc cases, these coefficients for each case seem not to have the tendency to decay. Figure 4.23 presents the variation of dimensionless K\ with the crack length ratio c/a for different aspect ration h/a. These results agree excellently with those given in [17]. Figure 4.24 gives the dimensionless tensile stress distributions on four horizontal sections. Because of the crack tip singularity, tensile stress increase dramatically near the crack tip. The isochromatics of stresses ax are plotted in Figure 4.25. Table 4.8 shows the dimensionless stresses calculated at selected boundary points. An excellent agreement can be found between the given and calculated values. The crack length ratio c/a and aspect ratio h/a are taken as the same value 0.5 in this case. A total of forty truncated terms are used in the series. Chapter 4. Numerical Experiments 48 0.008 n 0.007 0.006 i1 i r— ~i—r_ 40 terms 50 terms 80 terms 60 80 100 Number of collocation points Figure 4.21: The root-mean-square errors ||e|| with the number of collocation points (rectangle, uniform tension). 0.050 0.040 oooo n—f- - 1 — i — r - —i 1 1 r-40 terms 50 terms 80 terms 60 80 100 Number of collocation points Figure 4.22: The maximum errors HeH,*, with the number of collocation points (rectangle, uniform tension). Chapter 4. Numerical Experiments 49 M = m + 5 0.990142 -0.573283 0.0266888 0.146361 -0.248051 -0.374282 -0.0213386 0.0789041 -0.12184 -0.313915 :30 Bm + l 0.157647 0.125654 -0.520773 -0.0783619 1.64643 1.01943 -0.481084 -0.125902 2.08096 2.18395 M = m + 5 1.00079 -0.584494 0.0296892 0.147633 -0.25178 -0.380675 -0.0146433 0.0725379 -0.121763 -0.353751 :40 Bm+i 0.165562 0.124866 -0.528473 -0.0798253 1.68301 1.00166 -0.469402 -0.136316 2.3264 2.07408 M = 1.00312 -0.586718 0.0301527 0.148008 -0.252607 -0.3827 -0.0135544 0.0729004 -0.112773 -0.381019 :80 Bm+l 0.167125 0.12463 -0.529909 -0.0818353 1.69484 0.995543 -0.454491 -0.220696 2.43747 2.17493 Table 4.7: Ten leading coefficients with using 30, 40 and 80 truncated terms (rectangle, uniform tension). 03 04 05 Crack length ratio c/a 08 Figure 4.23: The variation of dimensionless Ki with the crack length ratio c/a (rectangle, uniform tension). Chapter 4. Numerical Experiments 50 Figure 4.24: Dimensionless tensile stresses ay on different sections (rectangle, uniform tension). 0.50 0 2 5 ooo -0.25 -050 -/ S—\ 1 / 1 / X \ / \ 1 \ 1 f \ \ / / V i /^\ \ I ! -' ' N/ | •! l ! / \ V—y / x \ / \ / ^ s ( \ J ) K s \ V / / s \ / N 0 5 0 0 2 5 0 0 0 0 2 5 0 5 0 Figure 4.25: The isochromatics of ay (rectangle, uniform tension). Chapter 4. Numerical Experiments 51 X 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 -0.3 -0.35 -0.4 -0.45 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 y 0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 Oca] 0.00056723 -0.00158 -0.00047336 0.0015484 0.00053905 -0.00039135 -0.00031257 1.0045 0.99792 1.0004 0.99976 1.0011 1.0016 0.99946 0.99755 0.99885 1.002 1.003 1.0002 0.99674 0.99711 1.0016 1.0045 1.0007 0.99437 0.99702 1.0066 0.99368 0.0046376 0.0025018 -0.0062063 -0.0072887 0.0064975 0.0089161 -0.007139 Teal -0.00072779 0.00067197 0.0011563 -0.0016913 -0.0017849 0.00437 -0.0011737 0.0085196 -0.0062623 0.0078947 0.0026789 -0.0050803 -0.0035647 0.0023071 0.0041482 0.00071913 -0.0029337 -0.0025513 0.00079758 0.0027082 0.00097008 -0.0017581 -0.0015538 0.00077884 0.00021712 -0.0021272 0.0042745 -0.0095338 -0.0012212 0.0047737 0.0037534 -0.0033303 -0.0074398 -0.0050355 0.016378 ^given 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ^gi ven 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Table 4.8: Comparison of boundary conditions at selected boundary points (rectangle, uniform tension). Chapter 4. Numerical Experiments 52 4.4 R e c t a n g u l a r P l a t e u n d e r C o n c e n t r a t e d Load ing Consider the compact specimen shown in Figure 4.26. It has been used as a standard specimen h h i r K c K \ k P 1 w 1 Vr I \ \ ) d C "^ -M i 1 ° 1 p Figure 4.26: The standard compact specimen for the test of fracture toughness. for the test of fracture toughness in the past. Numerical methods, mostly finite element methods have been applied to the problems with this specimen. Generally speaking the approach using stress function to the problems involving concentrated forces faces mathematical difficulties. In 1960s and 1970s Wilson[34], [35] published his results of stress intensity factors with the variation of crack length using a simplified mathematical model. He used a solid rectangle and replaced the concentrated load P acting on the holes with a uniform shear stress distribution on the left sides (Figure 4.26, right). Using a similar stress function(Williams[33] stress function), Wilson solved the problem using boundary collocation method. The specimen we proposed here, as shown in Figure 4.27 and 4.28 is slightly different from the standard one. Following the idea of Vaughan and Wu in [28] and Wu[37], we replace the concentrated load P by a distributed shear stresses over certain length / on the left edges. We use two set-ups, one with the lugs welded at the top left edges (Figure 4.27) and the other with the lugs welded at the crack mouth (Figure 4.28). Chapter 4. Numerical Experiments 53 Figure 4.27: The specimen and the mathematical model for concentrated loading at crack mouth. , 1 1 1 t (o 1 • 1 1 \h o B i T • h ' A h , T ' ' • | : / • i • i t * i 0 D Figure 4.28: The specimen and the mathematical model for concentrated loading at upper side edge. Chapter 4. Numerical Experiments 54 Define the reference stress CT° "' r° " l i (4M^ where t is the thickness of the plate and / is the length of loaded edge. The stress intensity factor K\ is given by 3 . (dVl2 3 A h ftoid ( P r~\ , = 2 A I T V — U ^ J - (4-95) Define thereafter the reference stress intensity factor KQ C a s e 4 .4 .1 Loading at crack opening The mathematic statement of the problem can also be expressed as follows V 4 $ = 0 in O ABC DO $ = $ n = 0 on OA $ = 0, $ n = r0y if 0 < y < I, $ „ = T01 if / < y on AB (4.97) $ = r 0 (x + c), * „ = 0 on BC $ = r0la, $ n = r0 / on CD It must be pointed out tha t , strictly speaking, one cannot have constant shear at corners. To avoid stress discontinuity at corners, one may define the above shear stress distribution by "rounding" it off to zero at corners. In our numerical experiments, zero shear stresses are specified at upper corners B and C, and the crack opening corner /I is never used since it introduces the trivial stress conditions (free stress conditions are satisfied automatically along OA). The variations of dimensionless K\ with the crack length ratio c/a for different aspect ratio h/a are shown in Figure 4.29, 4.30 and 4.31 for different loading length ratios, l/h = 0.25, 0.5 and 1, respectively. These results agree with those in [17] very well. Chapter 4. Numerical Experiments 55 0 00 0 10 0 20 0.30 0.40 0 50 Crack length ratio c/a 0.60 0 7 0 0.80 Figure 4.29: The variation of dimensionless Ki with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.25 (rectangle, loaded at crack opening). 0 3 0.4 0 5 Crack length ratio c/a Figure 4.30: The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.5 (rectangle, loaded at crack opening). Chapter 4. Numerical Experiments 56 25 20 l/h = 0.75 15 </> c ifi c CO 0) fi V) -o <D N S 1 0 S 5 Z M/a=l h/a = 0.8 h/a = 0 fr h/a = 0.5 hla = 0,2 0.3 0.4 0.5 Crack length ratio c/a 0.8 Figure 4.31: The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 1 (rectangle, loaded at crack opening). To see the effect of the approximation of loadings, we varied the length of loading for fixed aspect ratios and calculated the values of dimensionless K\ by varying the crack length ratios. The results are presented in Figure 4.32 (with h/a = 0.5) and Figure 4.33 (with h/a = 0.67). It is found that K\ is nearly independent of /, the length of loaded edge. To see the efficiency of the method applied to this case, we tabulate the first ten coefficients of Am and Bn in 50, 140 and 180 terms(half of each used in the first and the second series) in Table 4.9. Once again, we do not see the decay of the coefficients in the higher order terms. Figure 4.34 shows the dimensionless tensile stresses on different horizontal sections. The stresses on the boundary are evaluated at the collocation and non-collocation points and are given in Table 4.10. It can be seen that a pretty good accuracy has been achieved over the evaluation points. C a s e 4 .4 .2 Loading at ends of left edges Chapter 4. Numerical Experiments 57 03 04 05 Crack length ratio c/a Figure 4.32: The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a — 0.5 (rectangle, loaded at crack opening). 0.3 0.4 0 5 Crack length ratio c/a Figure 4.33: The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a = 0.67 (rectangle, loaded at crack opening). Chapter 4. Numerical Experiments 58 ft I o (0 i i • — r y/h = 0.031 y/h = 0.125 y/h =0.281 y/h = 0.500 -_ —s-^**" s r-1 1 1 1 / i i i i i i i i *, i •"* -^ -—--x. vs > i i i i i i 2 -' -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 x/a Figure 4.34: Dimensionless tensile stresses on different sections (rectangle, loaded at crack opening). M -. Am+l 0.55607 -0.80596 0.26682 0.14199 -0.67346 -1.8707 1.2183 -1.4636 0.32418 -6.6031 = 50 5 m f l 0.44864 -0.006442 -0.98448 0.078128 6.5485 2.5376 -5.5559 15.698 -7.8198 121.81 M -An + J 0.56372 -0.81774 0.27052 0.14202 -0.6884 -1.8987 1.2447 -1.4661 0.29393 6.729 - 70 Bm+\ 0.45767 -0.011774 -0.98913 0.070653 6.7262 2.3858 -5.2248 14.304 -2.2075 103.36 M Am+k 0.56959 -0.82569 0.2729 0.14622 -0.69031 -1.9227 1.258 -1.4914 0.24987 -6.804 - 90 Bm + l 0.46131 -0.0093913 -1.0011 0.050937 6.8031 2.3715 -5.0762 13.898 1.1617 92.229 Table 4.9: Ten leading coefficients with using 25, 70 and 90 truncated terms (rectangle, loaded at crack opening). Chapter 4. Numerical Experiments 59 X 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 y 0.11538 0.15385 0.19231 0.23077 0.26923 0.30769 0.34615 0.38462 0.42308 0.46154 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.46154 0.42308 0.38462 0.34615 0.30769 0.26923 0.23077 0.19231 0.15385 0.11538 Ocal -0.0013638 0.001585 -0.00035975 -0.00032831 0.0011835 -0.00080093 0.00075082 -4.4336-05 0.00041876 -0.0021308 0.00011652 0.0018447 -0.0012833 -0.0015457 -0.001853 -0.0024848 -0.0030067 -0.0031431 -0.0025512 -0.0011713 -9.2777-05 -0.005495 0.023133 -0.023059 0.028598 -0.040713 0.049902 -0.042522 0.0055925 0.048641 -0.081676 TCal 0.00045077 -0.00089723 0.0007453 -3.4544-05 -0.00085599 0.0014721 -0.0017312 0.0017437 -0.0017707 0.0016557 -8.224-05 0.0033247 0.0012584 -5.2589-05 -0.00020663 -0.00028449 -0.00077199 -0.0018414 -0.00298 0.0035965 -0.0011059 -0.046193 0.0054228 0.012838 -0.031612 0.08281 -0.34007 -1.033 -1.0433 -0.94589 -1.0015 ^given 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 'given 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -1 -1 -1 Table 4.10: The agreement of boundary conditions at selected boundary points (rectangle, loaded at crack opening). Chapter 4. Numerical Experiments 0.00 0 10 0.20 0.30 0.40 0 50 Crack length ratio c/a 0.60 070 0 8 0 Figure 4.35: The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h — 0.25 (rectangle, loaded at left corners). The mathematic statement of the problem can also be expressed as follows V 4 * = 0 in O ABC DO $ = $ n = 0 on OA $ = 0, $ n =. 0 if 0 < y < h - I, $ n = T0y if y > h - I on AB (4.98) $ - T0(X + c), $ n - 0 on BC $ = rQla, # „ = T0l on CD The variations of dimensionless K\ with the crack length ratio c/a for different aspect ratio h/a are shown in Figure 4.35, 4.36 and 4.37 for different loading length ratios, l/h = 0.25, 0.5 and 1, respectively. As in Case 4.4.1, we varied the length of loading for fixed aspect ratios and calculated the values of dimensionless K\ by varying the crack length ratios. The results are presented in Figure 4.38 (h/a — 0.5) and Figure 4.39 (h/a -— 0.67). One can see that K\ is nearly independent of the length of loading. Chapter 4. Numerical Experiments 61 0.3 04 05 Crack length ratio c/a Figure 4.36: The variation of dimensionless K\ with crack length ratio c/a for different aspect ratio h/a with loading length l/h = 0.5 (rectangle, loaded at left corners). 25 P t520 2-v> a> CO -o CD gio o l/h = l Va = 0.2 i/a = 0.5 Va = 0.67 Va = 0.8 Va = 1 03 04 05 Crack length ratio c/a 0.8 Figure 4.37: The variation of dimensionless Ki with crack length ratio c/a for different aspect ratio h/a with loading length l/h = 1 (rectangle, loaded at left corners). Chapter 4. Numerical Experiments 62 15.0 112.5 ai10 0 w 7.5 x> N O z 5.0 2.5 0.0 h/a = 0.5 ^ ^ > — a — l/h = l l/h = 0.5 l/h = 0.25 0.0 0 1 0 2 0.3 0 4 0.5 Crack length ratio c/a 0.6 0 7 0.8 Figure 4.38: The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a — 0.5 (rectangle, loaded at left corners). 15.0 I1" V) S 1 0 0 I" 5 0 2.5 0 0 h/a - 0.6f & 67 l /h-1 l/h - 0.5 l/h - 0.25 l /h- 0.1 0.0 0.1 0 2 0 3 0.4 0.5 Crack length ratio c/a 0 6 0 7 0 8 Figure 4.39: The variation of dimensionless K\ with crack length ratio c/a for different loading length l/h with fixed aspect ratio h/a — 0.67 (rectangle, loaded at left corners). Chapter 4. Numerical Experiments 63 14 12 10 8 6 4 2 0 2 4 y/h = 0.031 y/h = 0.125 y/h = 0.281 y/h = 0.500 / • \ Figure 4.40: The dimensionless tensile stresses on different sections (rectangle, loaded at left corners). To see the efficiency of the method applied to this case, we tabulate the first ten coefficients of Am and Bn in 50, 140 and 180 terms(half of each in the first the and second series) in Table 4.11. Once again, we do not see the decay of the coefficients in the higher order terms. Figure 4.40 shows the dimensionless tensile stresses on different horizontal sections. The stresses on the boundary are evaluated at the collocation and non-collocation points and are given in Table 4.12. It can be seen that a pretty good accuracy has also been achieved over the evaluation points. Chapter 4. Numerical Experiments M = Am+\ 1.08094 -1.57502 0.523242 0.385584 -1.81628 -2.33845 0.173632 -0.399948 -0.254305 -16.5446 = 60 Bm + l 0.863543 0.0880142 -2.44303 1.71426 10.2407 4.98047 0.705137 -8.19895 70.5434 77.8475 M A n + i 1.09383 -1.5937 0.528956 0.391143 -1.83861 -2.36599 0.180076 -0.43536 -0.18451 -16.9962 = 80 Bm+i 0.873796 0.0898527 -2.47075 1.72615 10.3998 4.91442 1.15786 -9.43853 75.3924 66.0496 M = A , ! m + 5 1.09956 -1.60104 0.530592 0.394343 -1.84657 -2.38341 0.1883 -0.446795 -0.188036 -17.015 100 BmH 0.877041 0.0919236 -2.48403 1.72519 10.4699 4.91963 1.19005 -9.57849 76.2055 62.9429 Table 4.11: Ten leading coefficients with using 60, 80 and 100 truncated terms ( loaded at left corners). Chapter 4. Numerical Experiments X 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 y 0.11538 0.15385 0.19231 0.23077 0.26923 0.30769 0.34615 0.38462 0.42308 0.46154 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.46154 0.42308 0.38462 0.34615 0.30769 0.26923 0.23077 0.19231 0.15385 0.11538 "cal 0.0052404 -0.00066613 0.0013158 0.0048854 -0.0026731 0.005876 -0.0014475 0.0043205 -0.00079751 0.0011952 -3.6139-05 0.0035125 0.0032861 0.0020485 0.0030147 0.0046834 0.0056906 0.0060288 0.0071065 0.0034517 1.7858-05 0.0051159 -0.021281 0.01865 -0.021694 0.013164 -0.0044496 -0.010735 0.0011471 0.026304 -0.060878 7"cal 0.0011362 0.001114 -0.0021653 0.0014494 0.00021394 -0.0017253 0.0025468 -0.0031646 0.0050793 -0.010343 -0.00045182 -0.0029401 -0.0012881 -0.00063921 -0.00075404 -0.00024142 0.00096593 0.0019508 0.0018615 0.00094688 0.0019654 -1.0046 -1.0001 -1.0045 -0.98735 -1.0203 -0.97869 -1.0039 -1.0241 -0.96004 -1.0087 ^given 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 "^given 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 Table 4.12: The agreement of boundary conditions at selected boundary points ( loaded at left corners). Chapter 4. Numerical Experiments 66 4.5 Rectangular P la te under In-plane Pure Bending Consider the specimen of a cracked thin plate simply supported at two ends, as shown in Figure 4.41. Apply concentrated forces P parallel to the crack. As shown Figure 4.41, assume a linearly M T r ^ ^W o W^ ^afi D M Figure 4.41: Rectangular plate with edge crack subjected to pure bending. variant tensile stress cr(x) on the section, on which P is acting, antisymmetrically about the centerline of the plate, that is (with the origin located at the crack tip) (* + <-f) (4.99) (4.100) o(x) = -2CT 0 ( x + c - - ) / a where a 0 is the reference stress. From equilibrium condition, we find _ 6 M _ 3(b0 - b)P where t is the thickness of the plate and M is the bending moment across the section on which P is applied. Here assume P is far enough from the crack so that a nearly pure bending state can be assumed in the cracked plate under concern. For simplicity of expression and comparison, we use M as the reference applied load in the following discussion. Chapter 4. Numerical Experiments 67 With the assumption of normal stress distribution mentioned above, the mathematical form of the problem can also be stated V 4 $ = 0 in O ABC DO $ = $ n = 0 on OA $ = $ n = 0 on AB (4.101) $ = ^ 2 ( x + cf (2(z f c) 3 a ) , $ n = 0 on AC 6a $ = - — , $ „ = 0 on CD 6a The stress intensity factor is given by 3 , /oV / 2 ,s— 3 , f2Trd\1/2 / 6 M ,-\ . . = iM—) (PI^) ( 4 1 0 2 » where d is the geometric scaling metric. Also, define the reference stress intensity factor 6M A ' o - ^ (4-103) Dimensionless K\ is then expressed by K\jK0 as usual. The variations of dimensionless K\ with the number of collocation points using different number of truncated terms in the stress function are plotted in Figure 4.42. For each case corresponding to certain number of truncated terms, a steady state in K\ is observed as enough collocation points, forty, say, are used. The variation of dimensionless Ki with the aspect ratio h/a for a fixed crack length ratio c/a — 0.5 using 40 truncated terms of the stress function is shown in Figure 4.43. It is interesting to see that the dimensionless A'I decreases a bit first from h/a — 0.5 and when h/a > 0.75, achieves a constant value of around 1.5. This implies that when "pure bending" status is obtained by using a relatively long beam plate, the stress intensity factor will be independent of the aspect ratio. In Figure 4.44, the variations of dimensionless K\ with the cn.rk length ratio c/a — 0.5 for different aspect ratio h/a are presented. The results shown in the figure agree with the Chapter 4. Numerical Experiments 68 2.0 45 0) c CD 1.5 <j> S 1.0 N C g 03 C CD E C o 0.5 -0.0 - $ = = $ = ^ O - o _l_ - o - o 20 terms 30 terms 40 terms 50 terms 80 terms _i_ _1_ 1 r -i_ 20 25 30 35 40 45 Number ot collocation points 50 55 <) Figure 4.42: The variation of dimensionless K\ with the number of collocation points (rectangle, pure bending). 200 i i • o 175 I £ 1.50 CO m 1.25 N C o 2 1 00 (D E o 2 0 75 -o- - O -0.50 '—'—'—'—'—'—'—'—L- ' • ' 0 50 0.75 1.00 -o-1 25 Aspect ratio h/a 1.50 1.75 -o 200 Figure 4.43: The variation of dimensionless K\ with the aspect ratio hja for fixed crack length ratio c/a — 0.5 (rectangle, pure bending). Chapter 4. Numerical Experiments 69 3.0 o2.5 f — 20 I 1 5 g to c Q) E 0.5 -i ^ H =3= h/a h/a; 0.5 : 1 A - " Jk J ^ k 000 010 0.20 0.30 0.40 0.50 Crack length ratio c/a 0.60 0.70 080 Figure 4.44: The variation of dimensionless K\ with crack length ratio c/a for h/a = 0.5 and h/a > 1 (rectangle, pure bending). results given in [17] very precisely for c/a > 0.15, while for the extreme case when c/a < 0.15, significant difference is observed. The first ten coefficients Am and Bn from totally 34, 40 and 80 truncated terms used in computations are given in Table 4.13. We see in this case, the coefficients vary smoothly. The dimensionless tensile stresses on four horizontal sections in the plate are shown in Figure 4.45 for c/a — 0.5 and h/a = 1. In Table 4.14 the stresses on the boundary are evaluated at a set of selected points. It can be seen that the relative errors in nondiniensional stresses are less than 0 .1%. Chapter 4. Namerical Experiments 70 M = Am+k 0.483914 -0.393837 -0.181152 -0.0616686 -0.0591681 -0.0653153 -0.0366069 -0.0501515 -0.0283842 -0.0320773 34 Bm+i 0.0973237 0.289767 0.214795 0.11066 0.302639 0.213497 0.282267 0.285204 0.235194 0.263675 M = 0.490375 -0.398203 -0.182004 -0.0638293 -0.0577243 -0.0692128 -0.0351316 -0.0572441 -0.0321942 -0.041434 40 Bm+l 0.0989372 0.28814 0.220046 0.103604 0.3077 0.212379 0.284011 0.334869 0.249896 0.395619 M = 0.500044 -0.405125 -0.184082 -0.0661999 -0.0573942 -0.0740365 -0.0292523 -0.0688286 -0.0194089 -0.0622632 80 Bm+i 0.101833 0.288532 0.228107 0.0970704 0.330756 0.187677 0.316849 0.321139 0.241421 0.525205 Table 4.13: Ten leading coefficients with using 34, 40 and 80 truncated terms (rectangle, pure bending). (0 a) 4.0 3.0 2.0 S 1.0 <D N "5 c o E T3 C o Z 0.0 -1.0 2.0 -3.0 40 -0 i 1 1 1 1 . _ . y/h = 0.04 ___ y/h = 0.16 -_ - y/h = 0.36 -_ - y/h = 0.64 — y/h = 1 i i i — i — ( * 50 -0.25 000 025 x-coordinates of each section 050 Figure 4.45: Dimensionless tensile stresses ay on different sections (rectangle, pure bending). Chapter 4. Numerical Experiments 71 X 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 -0.3 -0.35 -0.4 -0.45 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 y 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0.875 0.75 0.625 0.5 0.375 0.25 0.125 Ocal -0.00014638 0.00038623 -1.0077-05 0.00014909 0.00028852 -4.9881-06 0.00020009 -1 -0.89999 -0.80005 -0.69981 -0.60033 -0.4997 -0.39997 -0.3004 -0.19963 -0.09988 -0.00051902 0.10029 0.20033 0.29938 0.40021 0.50043 0.59922 0.70064 0.79964 0.90009 0.99998 -0.0011523 -0.00065734 0.00072381 -0.00034816 0.00032243 -0.0022849 0.0029532 Teal 5.5477-06 1.9035-06 -3.8212-05 0.00014028 -0.00026641 0.00012597 0.00053261 -2.5821-05 0.00013961 0.00023647 -0.00013912 0.00054924 -0.00022479 0.00022167 0.00032311 -9.1689-05 0.00018296 0.00023689 2.2281-05 0.00019658 9.3674-05 9.2333-05 0.00034738 -0.00025412 0.00054307 -0.00018203 0.00030451 -8.7867-05 -0.00053037 0.00010821 -0.00073184 0.0003573 0.00074465 -0.0009319 -0.0017412 ^given 0 0 0 0 0 0 0 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 -0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0 0 0 0 0 0 Tgi ven 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Table 4.14: The agreement of boundary conditions at selected boundary points (rectangle, pure bending). Chapter 4. Numerical Experiments 72 4.6 Rectangular P la te under Concentrated Shearing Consider the specimen shown in Figure 4.46 subjected to two self-equihbrium couples. The P Figure 4.46: A simple specimen for Mode II test and numerical model. problem is solved in the upper half region cut along the crack and symmetry axis(indicated with dashed lines). The concentrated force P is replaced in the numerical model by a uniform shearing stress distribution at the right top edge. The concentrated force Q yielded at the pinned point to resist the sliding is replaced with a normal stress distribution on the left edge. Using the (dimensionless) antisymmetric series introduced in Chapter 2 in the truncated form M 4>(r, 0) = Y. AmTm [sin(m - 2)0 - sin(m0)] m N + 5Z 5"r" in s i n ( n - 2)e - (n -2)sin n0} n the stress intensity factor K\\ is given by (4.104) Ku = Aapo ( - j \f^K 2 Lp \ a J \at (4.105) Chapter 4. Numerical Experiments 73 where p0 — P/at is the reference stress or scaling factor and t is the thickness of the plate. As in the previous sections, define the reference stress intensity factor by P r-at If we approximate Q by a uniform normal stress distribution along AB and approximate P by a uniform shear stress distribution at the right end of BC over a length of Lp, the corresponding mathematical statement of the first biharmonic problem is stated as follows V 4 * = 0 in O ABC DO $ = $ n = 0 on OA 2 $ = M - > $ „ = 0 on AB (4.106) * Qoh2 x , . , ,. , j <s — ——, $ n = q0h it — c ^. x < a - c - Lp; $n = Po(x + c -\- Lp — a) + ah otherwise, on BC ah2 $ = — + (y - h)(ah - p0Lp) * n ^ 0 on CD where q0 — 2p0Lp/h. Note again here that the shear stress distribution is expected to be "smoothed" out to zero at sharp corner B to avoid the discontinuity. The variation of dimensionless K\\ with the crack length ratio c/a for different aspect ratios h/a is given in Figure 4.47 with using 30 truncated terms. For large aspect ratios h/a > 1.0 A'n increases smoothly. However, when h/a < 0.5, A'n oscillate with c/a. As encountered in the Mode I cases, for small crack length ratios, the convergence of A'n is fairly slow and the collocation procedure has the difficulty to satisfy the specified boundary conditions uniformly, and relatively large errors were observed. It is also interesting to see that at for h/a < 0.5, there is a valley or turning point in the values of A'n for c/a ~ 0.5. Similar phenomenon was reported in Jones and Chisholm[ll], in which the Mode II stress intensity factors were obtained using a compound mixed mode specimen. The convergent behavior of A'n in the collocation procedure with the least squares approach were observed in most cases in which the crack length ratio c/a is not too small or too large. Chapter 4. Numerical Experiments 74 3.5 p o rti >. v> c <u c (0 in 1" 0) <» T3 o ^ in (0 0) c o in r (U f-Q 3.0 2.5 2.0 1.5 1 n O b 0.0 I - h/a = 0.4 - h/a = 0.5 • h/a = 1.0 - h/a = 1.5 O -— O --o— --o-~ , -—-O" o— - O ^> "O () _l_ 0.0 0.1 02 0.3 0.4 0.5 Crack length ratio c/a 0 6 0 7 0.8 Figure 4.47: The variation of K\\ with the crack length ratio c/a for different aspect ratios (rectangle, pure shearing). As an example, the variation of dimensionless Ku with the number of collocation points is presented in Figure 4.48 for c/a = 0.5 and h/a — 0.5 using 30 truncated terms in the series of Mode II stress function. The corresponding stresses evaluated at a set of selected boundary points are shown in Table 4.15. A sound agreement can be found between the given values and calculated ones. The isochromatics of shear stress TT0 are shown in Figure 4.49 The coefficients of the series of Mode II stress functions do not converge as rapidly as those for the Mode I cases. Table 4.16 lists the first ten coefficients in each of the truncated series for c/a — 0.5, h/a — 0.5 using totally 20, 30 and 36 truncated terms in the collocation procedure. One can see the large difference among the corresponding coefficients in the higher order terms. Comparing the experiments introduced in the previous sections, we had more difficulties with Mode II problems than Mode I. Particularly, we had fairly slow convergent rate of K\\ with respect to the total number of truncated terms in the series of Mode II stress function and with the number of collocation points. In order to obtain better agreement with the prescribed boundary conditions, we tried to include more terms and more collocation points. However, we Chapter 4. Numerical Experiments 75 2.20 I m 2 1 5 2.10 a> x> o 2 c o in c 0) E Q 2.05 2 00 ' <• -o-—1 I l _ 20 40 60 80 Number of collocation points 100 120 Figure 4.48: The variation of Ku with the number of collocation points for c/a — h/a = 0.5 (rectangle, pure shearing). 0.50 0.25 ooo -0.25 0 50 • — ' — ' — > - _1 1 u_ 0 50 0 2 5 0 0 0 0.25 0 5 0 Figure 4.49: The isochromatics of shear stresses around the crack tip and the corners of loaded edges (rectangle, pure shearing). Chapter 4. Numerical Experiments 76 0.5 0.5 0.5 0.5 0.5 0.5 0.5 y 0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 °cal 0.010676 0.035494 0.054743 0.044253 -0.0083181 -0.057426 0.0032644 Teal -0.083626 -0.033371 0.02975 0.074948 0.064027 -0.040785 0.0015905 ^given 0 0 0 0 0 0 0 ''"given 0 0 0 0 0 0 0 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 0 -0.0625 -0.125 -0.1875 -0.25 -0.3125 -0.375 -0.4375 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.011421 -0.023588 -0.064611 -0.023397 -0.060321 -0.1507 -0.14438 -0.06928 -0.048317 -0.058015 -0.013991 0.045094 0.030096 -0.023094 0.031004 0.90716 0.94875 1.0211 0.69926 0.16741 -0.11329 -0.062871 0.052908 0.073765 0.04358 0.019205 0.0093409 0.007047 -0.0079383 -0.030721 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 0.9871 1.0092 0.99337 0.99058 1.0052 1.0211 1.0118 -0.0021814 0.011406 0.0032638 0.0018163 0.010723 0.013807 0.018767 Table 4.15: Comparison of boundary conditions at selected boundary points (rectangle, pure shearing). Chapter 4. Numerical Experiments 77 M = m+5 -0.17309 0.10335 0.0050966 -0.025773 -0.10191 -0.086894 -0.084073 -0.12319 0.18286 -0.68854 20 Bm + l 0.012901 -0.025471 -0.012997 0.010935 0.02716 0.059426 -0.032214 -0.12937 -0.14979 -0.11705 M ^ Am+L -0.17313 0.094361 0.010496 -0.058125 -0.0014494 -0.0040579 0.046879 0.39739 0.27809 0.027594 30 Bmi\ 0.015173 -0.031377 -0.0043019 -0.036009 0.049767 -0.021251 0.051771 -0.19505 0.021758 -0.026131 M --Am+k -0.17293 0.097938 0.010584 -0.037888 -0.042662 -0.03267 -0.052534 0.11644 0.24402 0.10508 = 36 Bm+i 0.012178 -0.032495 -0.0072063 -0.018262 0.045916 0.010074 0.053517 -0.1907 -0.014547 -0.19054 Table 4.16: Ten leading coefficients with using 20, 30 and 36 truncated terms (rectangle, pure shearing). then had to diminish the series to avoid the roundoff effects. Thus, during the course of the experiments, we were always trying to oscillate between increased terms and increased round off. Nevertheless, the Mode II results obtained in this sections have enriched the resources that can be used as numerical references for future potential laboratory work for fracture toughness tests. The results presented here are new and not reported before. Due to the lack of reference data either numerically and experimentally, further studies are expected in the future research. C h a p t e r 5 Discuss ion a n d Conc lus ion Numerical results in the previous chapter have shown the success in calculating stress intensity factors in most Mode I(symmetric) cases with the use of series stress function. The work presented in this thesis is an extension of the previous work of Vaughan and Wu[28] on a circular disc to more general cases. New results of stresses and stress intensity factors have been obtained for a rectangular plate subjected to different loadings. In particular, stresses and stress intensity factors in a compact specimen, that is, a rectangular plate under concentrated tensions, have been investigated. This example has provided us with a new and convenient experiment for determining the fracture toughness of a given material in the laboratory. Also available are the useful results for pure bending, in which an in-plane simply supported cracked rectangular plate is loaded with concentrated forces parallel to the crack. With concentrated loads, which are easy to produce in a laboratory test, this case has provided another possible way of determining fracture toughness of materials. As one has noticed, boundary conditions in all of the model problems are specified by appropriate stress distributions. In particular, in the case of compact specimen the concentrated tensions are replaced by the distributions of shearing stresses in two parts of the loading edges. This approximation is based on the assumption that small changes in the data on the boundary would cause only correspondingly small changes in the stresses at the crack tip. In other words, this assumption is made in the St. Venant sense. It is interesting to see that , in compact tension cases, the stress intensity factors at the crack tip are independent of the way of distributed loads assumed along the left side loading edges. This means, from a practical point of view, that the stresses and stress intensity factors in the vicinity of crack tip has been well captured by 78 Chapter 5. Discussion and Conclusion 79 our numerical models. This observation encourages one to choose relatively simple numerical model, for instance, a distribution of shearing stress over the whole left side edge, in calculating the stress intensity factors. Numerical experiments have also shown that the smoother the distributed stresses on the boundary are, the better the agreement of the (least squares) solution and given boundary conditions is. Convergent and stable values of the stress intensity factors can only be obtained when the boundary conditions are reasonably assured. Here by convergent and stable values we mean the values of stress intensity factors independent of the number of truncated terms used in the series of stress function and the number of collocation points employed respectively. In the cases where concentrated forces are approximated by a small range stress distribution, the collocation of stresses may not be able to obtain a good approximation of the solution. With this in mind, we also used alternative forms of boundary conditions in terms of the stress function and its outward normal derivative(See Appendix C), which are more general for abstract and theoretical studies in mathematics and other fields. By integrating the second derivatives of the biharmonic function specified on the boundary, which are stresses in our cases, we may "round off" the "corners" or sudden step changes in the boundary data in terms of stresses, and therefore may expect to improve our numerical results. By doing so, however, we must remind ourselves that we achieve the accuracy up to the first, derivative rather than the second derivatives on the boundary, and a good approximation to the solution of the biharmonic function on the boundary does not always guarantee a good approximation of the second derivatives or stresses. To our particular problems, the boundary collocation is better in general over the Fourier expansion method[28] and the finite element methods based on the following concerns. First of all, the boundary collocation method itself is very straightforward and easy to implement in numerical computations. Secondly, unlike the Fourier expansion and finite element methods, no integrations are required either analytically or numerically. The coefficient matrix is obtained directly from boundary collocation and this makes the boundary collocation approach very Chapter 5. Discussion and Conclusion 80 efficient in practical uses. Of more importance is that the collocation approach we used employs a series of basis functions that are biharmonic everywhere and simulate a stress-free re-entrant boundary around the crack tip. The singularities of stresses in the neighbourhood of the crack tip are then captured by the basis functions. However, the boundary collocation method faces numerical difficulties when dealing with the problems in which there exist significant changes in the boundary data, for example, concentrated forces in our cases. Concentrated forces must be approximated by a distribution of stresses. Accurate solutions to the stress field is not available near the loaded area. Also the series stress function used in our collocation approach is derived from an infinite domain with a reentered boundary. With bounded finite domains as we have in our cases, these stress functions should be considered as local solutions in the neighbourhood of the crack tip. This may explain why the numerical experiments deteriorate in the cases when the length of crack is too small or too large. The results for the Mode II model in Section 4.6 are new and not reported before. Compared to the compound specimen used by Jones and Chisholm[ll], the specimen and the corresponding numerical model given in this thesis is simpler. Pure shearing at the crack tip is emulated and antisymmetric stress state is assured. Due to the lack of reference data , further study in this model problem is needed both numerically and experimentally. Having succeeded in many cases, particularly for Mode I loading, we have also found that the boundary collocation approach with the least squares method is not successful universely. It does not take care of Mode II problems as nicely as the Mode I cases. We have also found different convergent rates of the series evaluated at different collocation points. Generally speaking, at the points far from the crack tip at which the origin of the polar coordinates is located the series have slower convergence rates than those at the points close to the crack tip. In some cases, we have observed strong oscillations in the series evaluated at the points near the two corners of the rectangle. Another valuable experimental aspect is the special effort that has been made in solving the system of Unear equations obtained from boundary collocation. In order to obtain better Chapter 5. Discussion and Conclusion 81 approximation to the given boundary values, large numbers of collocation points are used in each case and this leads to an overdetermined linear system. Because of this large number of collocation points densely located and the use of many terms in the power series of the basis functions, the coefficient matrices in our experiments almost all suffer from the problem of ill-conditioning. To obtain the least squares solution to an overdetermined linear system with the above problem, special efforts must be made(See for example, [7|,|12] and [24]). Rather than using the associated normal equation, we applied the QR factorization to the system and solved it directly. An example that shows the failure of using normal equation in solving our problems and the success of using QR method is given in Appendix E. Finally we are left with some problems for future study. In order to have a guidance in selecting collocation points, we need an efficient priori error estimation theory. Unfortunately, unlike finite element methods few error estimations for boundary collocation method are avail-able. Such a work may be somewhat theoretical and challenging but will be definitely helpful in numerical calculations. Special efforts may also be helpful to make in the future to smooth out the relatively large errors near the corners of the rectangle. A possible way is perhaps to use the technique of spectral analysis combined with multigrid or multi-level approaches. Further study in a theoretical ground is needed. Since the main purpose of our work is to provide an efficient way of testing for fracture toughness of materials in laboratory tests with numerical support, the major work in the future will be the laboratory experiments. Having had the stress intensity factors given in terms of apphed forces and given geometries, we are ready to obtain the fracture toughness by increasing the specified forces to critical values. However, since the numerical results we have obtained are all based on the linear elastic theory, plastic deformation must be taken into account in practical experiments. B i b l i o g r a p h y R. N. Banerjee and M. W. Benson. "An approximate Inverse Based Multigrid Approach to Biharmonic Problem". Inter. J. Computer Math., 40:201-210, 1991. M. J. M. Bernal and J. R. Whiteman. "Numerical Treatment of Biharmonic Boundary Value Problems with Re-entered Boundaries". Cornp. J., 13(1), 1970. P. E. BJ0rstad. "Numerical Solution of the Biharmonic Equation". PhD thesis, Stanford University, 1980. L. W. Ehrlich. "Solving the Biharmonic Equation as Coupled Finite Difference Equations". SIAM J. Numer. Anal., 8(2):278-287, 1971. L. W. Ehrlich. "Coupled Harmonic Equations, SOR and Chebyshev Acceleration". Math. Comp., 26(18):335-343, 1972. L. W. Ehrlich and M. M. Gupta. "Some Difference Schemes for the Biharmonic Equation". SIAM J. Numer. Anal., 12(5):773-790, 1975. G. H. Golub and C. F. Van Loan. "Matrix Computations". Johns Hopkins University Press, 1989. D. Greenspan and D. Schultz. "Fast Finite Difference Solution of Biharmonic Problems". Comm. ACM, 5(5):347-350, 1972. C. E. Inglis. "Stress in A Plate Due to The Presence of Cracks and Sharp Corners". Proc. Inst. Naval Architects, 60, 1913. G. R. Irwin. "Analysis of Stresses and Strains near The End of A Crack Traversing A Plate". Trans. ASME, Appl. Mech., 1957. D. L. Jones and D. B. Chisholin. "An Investigation of The Edge-Sliding Mode in Fracture Mechanics". Engng. Frac. Mech., 7:261 270, 1975. C. L. Lawson. "Solving Least Squares Problems". Prentice Hall, Englewood Cliffs, 1974. J. W. McLaurin. "A General Coupled Equation Approach for Solving the Biharmonic Boundary Value Problem". SIAM J. Numer. Anal., 11(1):14 33, 1974. H. Motz. "Treatment of Singularities of Partial Differential Equations by Relaxation Meth-ods". Q. Appl. Math., 4:371 377, 1946. J. R. Rice. "A Path Independent Integral and the Approximation Analysis of Strain Concentration by Notches and Cracks". J. Appl. Mech., 35:379 386, 1968. 82 Bibliography 83 [16] J. ft. Rice. "Mathematical Analysis in the Mechanics of Fracture". In H. Liebowitz, editor, Fracture Volume II. Academic Press, 1968. [17] D. P. Rooke and D. J. Cartwright. Compendium of Stress Intensity Factors. Procurement Executive, Ministry of Defense, London, Her Majesty's Stationary Office, 1972. [18] D. P. Rooke and J. Tweed. "The Stress Intensity Factors of A Radial Crack in A Point Loaded Disc". Int. J. Engng. Sci., 11:285-290, 1973. [19] G. C. Sih. "Handbook of Stress Intensity Factors". Lehigh University, 1973. [20] J. Smith. "The Coupled Equation Approach to the Numerical Solution of the Biharmonic Equation by Finite Differences, I". SIAM J. Numer. Anal, 5(2):323-339, 1968. [21] J. Smith. "The Coupled Equation Approach to the Numerical Solution of the Biharmonic Equation by Finite Differences, I". SIAM J. Numer. Anal, 7(1): 104-111, 1970. [22] I. N. Sneddon. "The Distribution of Stress in The Neighbourhood of A Crack in An elastic Solid". Proc. Roy. Soc, pages 187 229, 1946. London. [23] I. N. Sneddon. "The Application of Integral Transforms in The Theory of Elasticity". Springer-Verlag, New York, 1975. [24] G. W. Stewart. Introduction to Matrix Computations. Academic Press, New York, 1973. [25] G. Strang and G. J. Fix. An Analysis of The Finite Element Method. Prentice-Hall, Inc., 1973. [26] M. Vajtersic. "A Fast Algorithm for Solving the First Biharmonic Boundary Value Prob-lem". Computing, 23:171 178, 1985. [27] M. Vajtersic. "Some Iterative Poisson Solvers Applied to Numerical Solution of The Model Fourth-Order Elliptic Problem". Aplikace Matematiky, 30(3):176 185, 1985. [28] H. Vaughan and Q. Wu. "Stresses in A Circular Disc Containing A Radial Crack". Q. J. Mech. Appl. Math., 44(Pt. 3):413-428, 1991. [29] H. M. Westergaard. "Bearing Pressures and Cracks". Trans. ASME, J. Appl Mech., 61:49-53,1939. [30] J. R. Whiteman. "Numerical Treatment of a Problem from Linear Fracture Mechanics". In A. R. Laxmoore and D. R. J. Owen, editors, Proceeding of the 1st International Conference on Numerical Methods in Fracture Mechanics, Swansea, 1978. [31] J. R. Whitenian. "Finite Difference and Singularities in Elliptic Problems". In J. Gradwell and R. Wait, editors, A Survey of Numerical Methods for Partial Differential Equations. Clarendon Press, 1989. Bibliography 84 [32] J. R. Whiteman. Two dimensional biharmonic problems. In J. Gradwell and R. Wait, ed-itors, A Survey of Numerical Methods for Partial Differential Equations. Clarendon Press, 1989. [33] M. L. William. "Stress singularities Resulting from Various Boundary Conditions in An-gular Corners of Plates in Extension". J. Appl. Mech., 74:526 528, 1952. [34] W. K. Wilson. "Analytic Determination of Stress Intensity Factors for the Manjoine Brittle Fracture Test Specimen". Technical Report WERL 0029-3, Westinghouse Research Report for AEC, 1965. [35] W. K. Wilson. "Stress Intensity Factors for Deep Cracks in Bending and Compact Tension Specimens". Eng. Frac. Mechs, 2(2):169-171, 1970. [36] L. C. Woods. "The Relaxation Treatment of Singular Points in Poisson's Equation". Q. J. Mech. Appl. Math, 6(Pt . 2):163 185, 1953. [37] Q. Wu. "Stresses in a Circular Disc Containing a Radial Crack". Master's thesis, Univer-sity of British Columbia, 1990. A p p e n d i x A T h e Special Crack S t r e s s Function and S t r e s s e s This appendix gives the details of stress expressions derived from the crack-surface-stress-free stress functions for Mode I and Mode II. A . l Stress Funct ions — M o d e I For Mode I cases, the general solution takes the form in polar coordinates with the origin located at crack tip </>(r, 0) = >T Amrm [m cos(m - 2)0 - (m - 2) cos(m0)] I Y, BnTn |cos(n - 2)0 - cosn0] (A.l07) n with m = | , | , . . . , and n — 2 , 3 , . . . . Here we adopt the index notation used by Vaughan and Wu in [28] although alternative forms may be more convenient. By the general relations d24> d f\a<f>\ the general expressions for stresses are given by _ Iffy j _ d ^ > _ d2(j  r dr r2 d62 ' oT = Y Amrm"2m(m 1) [(m - 2) cos m0 - (m 4) cos(m - 2)0] m + ^ 5 n r " ~ 2 ( n - l ) [ n c o s n 0 - ( n - 4 ) c o s ( n - 2 ) 0 | , (A.109) n ag — Y^ Amrm 2m(rn-- l ) [ m c o s ( m - 2)0 (m 2)cosm0j m + ] T f i n r " 2n(/i - l ) [ c o s ( n - 2 ) 0 - c o s n0 ] , (A.110) n TTQ — _ y ^ ^ m ^ m 2fn(rn - l ) (m - 2) [sin rri9 - sin(m 2)0] 771 85 Appendix A. The Special Crack Stress Function and Stresses 86 - ^ f l n r n " 2 ( n - l ) [ 7 i s i i i r i 0 - ( n - 2 ) s i n ( n - 2 ) 0 ] , ( A . I l l ) n Note that in the vicinity of crack tip, as r —* 0, (£~ A 3 J~5 ( - c o s - + c o s - 0 | (A.112) Y 3 \2 2 2 y V ; and correspondingly, the asymptotic forms of stresses ar ~ g A | - r - s ( 5 c o s - - c o s y J + 0 ( 1 ) , (A. 113) <je ~ %Alr ~2 ( 3 c o S 2 + C O S y ) + ° ( 1 ) ' (A.114) r r f l ~ ^ | r - 5 f sin y + s i n - ) + 0 ( 1 ) . (A. 115) In particular, in front of the crack tip at 0 — 0, as r --> 0 Or = ae ~ 7.4?r""2 , Trg = 0. (A .116) 2 2 A . 2 S t r e s s Funct ions - M o d e II In mode II cases, the stress distributions are ant isymmetric and the corresponding general solution of the stress function is of the form 4>(r, 0) = YJ Amrm lsin(m - 2)0 - sin(mfl)] m + Y^BnTn[nsm(n - 2)0 - (n - 2)sinn0] (A. 117) n with 771 = | , | , . . . , and n — 3 , 4 , . . . . The general expressions of stresses for Mode II are <rr = ^2Amrm'2(Tn - l)[m sin m0 - (m 4)s in(m 2)0} m + ^ 5 n r n _ 2 7 i ( n - l ) [ ( n - 2 ) s i n n 6 l ( ri — 4) sin(n — 2)0], (A.118) n ag = V ] ylmrm~277i(7n - 1) [sin(m - 2)0 - sin m.0] m + ^ f l n r n 2 n(n l)[7isin(n - 2 ) 0 ( n 2) sin n0] , (A.119) n Tre = - *y\Arnrm~'2(m, -• l)[(7n - 2)cos(m 2)0 mcosmfl] m - ^ Bnrn 2 71(71 - l)(ri - 2) [cos(n - 2)0 - cos7i0j , (A.120) Appendix A. The Special Crack Stress Function and Stresses 87 Similar to Mode I, in the vicinity of crack tip, as r - > 0, I r l ( ' 0 . 3 S ~ - i s r J { sin - -\- sin — 6 Y 2 \ 2 2 and correspondingly, 1 i / ft ^/J*1 <rr ~ - - ^ 3 7 - ^ 2 f 5 s i n - - sin —- ) f 0 ( 1 ) , a e ~ - - A J T - 2 I s m - + sin — J -I-0(1), 1 , - 1 / e 36»\ ^ v r r e ~ - A | r 2 I cos - + 3 cos — I + 0 ( 1 ) . In front of the crack tip at 0 — 0, OT ~ Oy — 0 , TTg ~ A 3 T- 2 . 2 A . 3 Stresses in Cartes ian and Polar Coordinates In two dimensional problems, denote the matrix rep of stress tensor in Cartesian by S and the matrix rep in polar coordinates by ' x ' xy CT, TT TT8 Of) Let T be the change of basis matrix (from polar to Cartesian coordinates), that is T -cos 0 sin 6 sin 9 cos 6 Then, the transform if of the following form (A.121) (A.122) (A.123) (A.124) (A.125) (A.126) (A.127) (A.128) S ^ T1 ST. (A.129) Appendix A. The Special Crack Stress Function and Stresses 88 Stress components can be written explicitly as follows ox — oT cos2 9 - 2TrflCos0sin 0 + og sin2 6, (A.130) Oy — aT sin2 6 + 2T r^cos^sin^ + oe cos2 d, (A.131) rxy — (aT - ag) cosOshiO i~ TTO(COS2 8 ~ sin2 0). (A.132) Given the stress expressions from the stress function defined earlier in polar coordinates, the stresses in Cartesian coordinates can then be obtain via (A.129) or (A.130) through (A.132) for the applications in which Cartesian coordinates are more convenient. A p p e n d i x B A n Alternat ive Form of B o u n d a r y C o n d i t i o n s Instead of being specified by the second derivatives of the biharmonic function, namely stresses, an alternative form of boundary conditions, which is more convenient for mathematical analysis, can also be specified by the biharmonic function and its (outward) normal derivative. B . l T h e General Form of Biharmonic P r o b l e m s of T h e F i r s t T y p e Consider the general boundary conditions on an arbitrary curved boundary as shown in Fig. (B.50), in terms of stress function (j> Figure B.50: An arbitrary boundary C of certain region ft in x — y plane 4>yynx - 4>xyny - a(s)nx - T(s)ny = ti, (B.133) ~4>xynx + <t>xxny — cr(s)ny -\ T(S)TIX - t2> (B.134) 89 Appendix B. An Alternative Form of Boundary Conditions 90 nx -ny -dy ds' dx ds' d A ds -- h, - h-where the normal and shearing forces &(s) and T(S) are given and the direction cosines are given by (B.135) (B.136) The boundary conditions can be written in the form d (B.137) (B.138) Integrating these two equations along the boundary, we obtain the expressions of the first derivatives of the stress function <f> at an arbitrary point s in the following form 4>x{s) = - / h<ls f dx f dv = / o{s) — ds - / T{s)-f-ds^ /!{*) +cu (B.139) Jsu ds Jao ds 4>y(s) = / tids r" dv f dx = o{s)-f-ds I- T(s)~ds ..- f2{s) \ c2, (B.140) Jsu ds J3u ds Notice that if the region under concern is multiply connected, the integration has to be done over each contour that forms the boundary of the region, and consequently the integrals will not be single-valued. If this is the case, some arbitrary choices of the stress function and its derivatives can be taken, as far as the stresses are uniquely determined. We now restate the boundary conditions by specifying the values of biharmonic function </> and its normal derivative <f>n on the boundary. Note that 44 — = (f>xnx + 4>yny dn dy dx — Yx ~ <Py ds ds = g(s). ( B . H i ) Appendix B. An Alternative Form of Boundary Conditions 91 Since d(f> - (j)xdx + (pydy, = f1(s)dx + f2(s)dy, (B.142) integration along the boundary gives dx [" dy #5) - j'flis)^dS + J°f2(s)fsds 'SQ "• •* J So = f(s) + const. (B.143) Now the plane elasticity problem can be stated as the following equivalent form in a region Q V 4 u( r ,0 ) = O in 0 < £ = / +cons t on <9fi, (B.144) 7 = 9 o i dil, dn which is more general in mathematics and has been investigated by many people in several other areas as well. Note tha t , the stress function and its normal derivative are not determined uniquely. How-ever, we need consider only the biharmonic functions with single-valued second partial deriva-tives. One of the advantages of using the alternative form of boundary conditions rather than the those in terms of stresses is that the new boundary conditions are in general "smoother" than the original ones and this may benefit the numerical computations. B.2 A Rectangular P la te under General Loadings Let us carry on to consider a rectangular plate with an edge crack under certain kind of loadings (Fig. B.51). The boundary conditions in terms of stresses on the edges are specifies as follows (for symmetry, only the upper half region () ABC DO is concerned) On OA: ay - 4>*x - 0, (B.145) Appendix B. An Alternative Form of Boundary Conditions 92 Figure B.51: A rectangular plate with an edge crack OA perpendicular to the edge AB. On AB: On BC: On CD: Oxy •- -<t>xy - 0; Ox - 4>yy ~ <T\(y), °xy --- '4>xy ~ n(y); °y - 4>xx - a-ii*), &xy - -<t>xy "; T2(x)\ 0<x - 4>yy = °3(y), °xy = -4>xy = T-i(y). (B.I46) (B.147) (B.148) (B.149) (B.150) (B.151) (B.152) We now integrate these functions to obtain the expressions of <f> and its normal derivatives on specific edges. 1) On OA(-c < x < 0,1/ = 0), since <pxx — 0 and 4>xy — 0, the first derivatives <fix and 4>y have Appendix B. An Alternative Form of Boundary Conditions 93 to be constants. Let <px = c01 and <fiy = CQ2, then the normal derivative </>0,n - (t>y - C02- (B .153) Here we use numeric subscripts 0, 1, 2 and 3 to denote edges OA, AB, BC and CD respectively. Integrate (j) along OA to get (with dy/dx — 0 on CM) <t>o -z J 4>xdx - cQlx + c03. (B.154) where C03 is some constant. For free stress condition required on the crack surface, without loss of generality, we set c0i = co2 = c03 = 0. Thus we have </> = 0„ = O. (B.155) 2) On AB(x = - c , 0 < y < h), integrate (pyy and (f>xy to get 4>i,x = </>o,xU - / ridy Jo = - P'ndy (B.156) Jo 01,1/ = <£o,j/U + / £ T l^2 / JO - / " f f i ^ (B.157) ./o 01 ^ 00 U + / 01,^1/ Jo = 00U "t / (0O,yU + / CTl*/j = 0oU I £ (£tridri\dy, (B.158) 01 ,n ~ 01,x - - / " n i t y (B.159) Jo Thus, Appendix B. An Alternative Form of Boundary Conditions 94 3) On BC{ - c < x < a — c,y — h), integrate cf>xx and (pxy to get 02,x = 01,x|.B + / O ^ Z = - I TXdy+ I a2dx (B.160) = / a1dy+ I r2dx (B.161) Thus, 4) On CL>, Thus, 02 = 4>I\B + (p2,xdx /x / yx > = <t>i\B \ I (f a2dA dx i 4>liX\B(x ^ c) (B.L62) 02,n = 02,y = / cridt/ + / r2dx. (B.J63) 03,x = 02,x|<7 - / T3rfj/ /•/» /-(a c) yy — - I Tidy \ / a2( ix / Tidy (B .164) JO J c Jk 03,y = 02,y |c + / O^dy Jh j-h yx ry - I aLdy \ / T2dx f / a3dy (B.165) J0 J-c V/i 03 = 02 |C + / 03,y< J^/ = 02|C + / f 03,y|c + / ^3^2/] - <h\c + f (JV <73dri\ dy-r <l>2,y\c(y - h) (B.166) 03,n = 03,x y/i /-(a c) /-y - - I ndy \ / <7 2 ^ / T3<fy- (B.167) Jo J-c Jh Appendix B. An Alternative For in of Boundary Conditions 95 B . 3 Rectangular P la te under Uniform tens ion Let us look at one of our model problems - a rectangular plate with an edge crack under uniform tensions o loaded on two opposite edges. In this case, 0\ — T\ = 0, o2 — o, T2 — 0, 0$ — r 3 — 0. 1) On OA: 4>=0, (B.168) ^„ = 0; (B.169) 2) On AB: 3) On BC: 4) On CD: 4> = I" [y(a1(ri)dr})dy = 0, (B.170) Jo Jo 4>n = ^0,X\A = 0; (B.171) 4> - J'J* {o2(£)d£)dx / {od£)dx a, - | ( x + c)2 , (B.172) 4>n = <t>l,v\B = 0; (B.173) Jo Jo = r r {ad£)dx a 2a\ (B.174) <Pn - 4>2,y\c /a c adx = oa. (B.175) Appendix B. An Alternative Form of Boundary Conditions 96 B.4 Circular Disc wi th Radial Crack Finally let us look at the boundary conditions for a circular disc with radial crack under general loadings. For simplicity, consider only the case when the crack length is equal to the radius of the circular disc. From the stress-stress function relations r dr r2 dO2' dr2' T dr on the edge of the disc, where r — a, the radius of the disc, \rdej a2 do a oroo Integrate TT$ along the edge of the disc to obtain the following equation * d<f> Jn 1 fe d2cf> r i r dd> i r Jo a* Jij 60 a Jo drdO do <fi{a,6) + c i - - b r -\ c2 a2 - ?«•••>-; (&L + '- (B'177) where c = c\ -f c2, ci , c2 are some constants. Also from the stress expressions Substitute (B.178) into (B.177) to obtain the following equation involving 4> and its second derivative with respect to 0 f Jo TT8(a,0)de=-^<p(a,6) + \ ( ^ ) ~ar(a,0) + c, (B.179) a2 \ 802 } Denote ifs(0) = (p(ay0)i then above equation becomes a second ordinary differential equation in v-^ + i, = a2 laT(a,0) - J%re(a,0)d9 - c) . (B.180) where ar, rrg are given. The boundary conditions for this differential equation are determined so that at crack opening, where 0 — -K, ip(%) ^ d>(a,ir) = 0, (B.181) V>fl(ir) = fo(a,0) = O. (B.182) Appendix B. An Alternative Form of Boundary Conditions 97 Here homogeneous conditions are assumed for the stress function and its normal derivative are assumed (see Section 2.2). Note that constant c is independent of the boundary conditions on the crack, it can then be set to zero. When 'ip is solved, the normal derivative of (f> along the disc edge can then be obtained from (B.177). However solving this differential equation is not easy in general, unless for pure mathematical discussions, in which the boundary conditions of this form may be convenient, this is seldom done in practice. A p p e n d i x C T h e Auxil iary Eigenvalue P r o b l e m s In Section 2.2, two auxiliary problems were introduced, in which the values of A in the equation sin(A 1)« - i (A - l ) s i n a (C.183) are to be found for a given inner angle a of the crack. For simplification we use notation A for A — 1 below. This equation can be written in a more symmetric form sin At* sin a — - - ± — - (C.184) or equivalently the eigenequation where &mz = Az (C.185) A . z sin a a = Xa. (C.186) (C.187) One may look at the more general problem of finding solutions of both A and a, where a £ R. Define function / ( x ) = ^ ( f ) . (C.188) x Assuming A is real, which is of more significance in general, the problem is reduced to looking for such pairs (x,y) tha t f(x) — f(y). Numerical solutions of A, denoted by Aj and A2 for Mode 1 and Mode II cases respectively, are tabulated in Table C and shown in Figure C.52. For z in (C.185) complex, solving for Re 2 is somewhat complicated. Muller's method is generally suggested for searching complex roots of (C.185). Although initial guess is not necessary in Muller's method, a guideline in searching for the solution of z with smallest real 98 Appendix C. The Auxiliary Eigenvalue Problems 99 a 360 350 340 330 320 310 300 290 280 Ai 0.5 0.500053 0.500426 0.501453 0.50349 0.506933 0.512221 0.519854 0.530396 A2 0.5 0.529355 0.562007 0.598192 0.638182 0.682295 0.730901 0.784441 0.84344 a 270 260 250 240 230 220 210 200 190 Ai 0.544484 0.562839 0.586279 0.615731 0.65227 0.697165 0.751975 0.818696 0.900044 A2 0.908529 N/A 1.06021 1.14891 1.24804 1.35949 1.48581 1.63053 1.79893 a 180 170 160 150 Ai 1 1.1251 1.28841 1.53386 A2 2 2.25313 2.65888 N/A Table C.17: The eigenvalues corresponding to Mode I and Mode II cases for selected vertex angles. 3.0 2.5 £ 2.0 CO 0 1 E ZJ E 1.0 0.5 0.0 V ^\ \ M o d e I M o d e II "^- ^ ~"~~"~-— ... 150 180 210 240 270 300 Inner angle of the vertex 330 360 Figure C.52: Variation of minimum real positive eigenvalues with the vertex angle for symmetric (Mode I) and antisymmetric (Mode II) cases Appendix C. The Auxiliary Eigenvalue Problems 100 part is somehow helpful. Let z — x -f iy, x,y £ A!, then sin z — (elz - e"lz) /2i. Define complex function f(z) = sin az ± zsin a (C.189) where a 6 R is the vertex angle mentioned earlier, then Re f(z) - - (e u y | eay) sin ax ± x sin a, (C.190) Im f(z) ^ — (e " y - e a y) cos ax ± 1/ sin a (C.191) To find the complex root of f[z\ a) — 0 with smallest real part , one may first draw contours of Re f(z) and Im f(z) in the z-plane within the range of interest for each given a, and find an initial guess by locating the intersection of contours Re / ( z ) = 0 and Im f(z) — 0. Then one can use efficient numerical solvers to obtain accurate solutions. A p p e n d i x D Coupled Poisson's Equations It has been shown(See, for example, [20], [2.11), that the solution of the biharmonic problem of the first type can be replaced by the solutions of the following coupled equations V 2u - v in Q V2v - 0 in ft (D.192) u = </> on dil un — 'ip on dQ where v is some smoothly defined harmonic function. Note that in order to solve the second Poisson's equation separately, the boundary condition for v must be specified. An iterative scheme, introduced, for example, in [13] for solving the above coupled problems takes the form y2u(fc|l/2) = v(h) j n JJ u(k+l/2) _ ^ o n d n u{k + i) _. au(k) _j_ n a)u(k+1/'2) (smoothing process) (D.193) V 2 ^ - i - i / 2 ) = 0 i n n v(ku/2) = &u(kn) 7(-JiJ*+1) - V ) on dU v(k+i) = pvW + (} _ p)v(k+i/2) ( s m o o t h i n g process ] where w0' — v(°> — 0, 7 / 0 is the coupling constant and a, /3 6 [0,1] are smoothing parameters . Iteration (D.193) converges and is consistent with the biharmonic equation for any 7 ^ 0[13]. For the finite difference method, the standard five-point scheme can be easily applied to each of the above Poisson equation and various Poisson solvers can be employed. Note that the discretization of (D.192) or (D.193) is separatable. it is trivial to discretize the first Poisson equation, while the second one requires that v be approximated along the 101 Appendix D. Coupled Poisson's Equations 102 boundary, which involves the approximation of Au and un along the boundary. For each grid point on the boundary, a fictitious point must be added outside the boundary. As treated in Vajtersic[26], we may use linear extrapolation at these fictitious points. Denoting the values of u and v by subscript 0 in a bi-index notation, wherever they lie on the left and lower sides of the rectangle, we demonstrate in the following how to approximate v along the crack boundary AO, shown in Figure D.53. Without loss of generality, let us look at the two points at the left B \ : A (0,0) -(<s) o (> -4- o -*- J c D Figure D.53: Finite difference grid for the upper half cracked rectangle O ABC DO. The crack is along OA and OD is the symmetry line. The boundary points are labeled with o and interior points are labeled with • . lower corner of the rectangle. With bi-index system, i = 0,. . . , Nx, j — 0,. . . , Ny, and uniform mesh size h, we have, at point (1,0) A/,ti|i>0 = 72 (0o." • '0 i .o + 02,o + u i , i - 4</>i,u) - -rj- l ^2(0°.° "" 301.° + 02-°)' I l / \ where A/, denotes the discrete Laplacian. Here linear extrapolation is assumed at the fictitious point (0, - 1 ) (not shown in the figure). From (D.193) with 7 = 2//t([l3], [27]) Vi,o = A/,u | i | 0 r{un\ 1,0 i>i,o) 2uh 1 2 Appendix D. Coupled Poisson's Equations 103 2xti , 2 ~ - + faM 4- ^ i , o (D.194) where (/>lfi = (</>0,o - 4^ 1 | 0 + 4>2lo)/h2. Similarly, *>o,i - ~^r I V-o.i f ^ o , i (D.195) where ^ 0 | i = (</>o,o - 4^0,i + $0,2) I h2. At point (1 , 1), setting A/Hi , i = 72 (^ ,1 t -"i,o I '"2,2 I '"A^H-I - 4 v M ) = 0, (D.196) gives 4v1(i - -y2,i - i>i,2 - v0,i I v1 |0 2 2 , - ^ " 1 , 1 + </>o,i + </>i,o + -^('0o,i I- V-i,o). (D.197) Similarly, we have -•i>i,i ) 4w2,o - v-*,i v2-i - v-ijo - ^ 2 , 2 1^ ,1 + 1^ ,0 - (D.J 98) Write the discrete form of (D.193) in matrix form Lu — h v ••[- b 2 Lv - —Mu + c where L denotes the N X N (N = Nx X Nx) matrix operator of Laplacian with five-point scheme and M = diag(D \ I, D, ...,D,D + I)NXN, with D — diag(l,0,...,(), l ^ x A / * and / is the identity matrix of order Nx. Here u and v are the unknown vectors containing the approximate values of u(x,y) and v(x,y) at grid points, and the vectors b and c arise respectively from the discretization of the first and second Poisson equations with special treatment (D.194) through (D.198). Appendix D. Coupled Poisson's Equations 104 From (D.193) an iterative scheme for solving the coupled systems takes the following form Zti( f c t l / 2) ... / iV f c ) - | b (D.199) / i l ) , a / ' ] (1 a)U<* f ,/2> (D.200) Lv{kil/2)- ^Mu(hil) -\ c (D.201) v(*+i) = £.„(*) t. ( i .. /3)u(fc+1/2). (D.202) The numerical aspects in solving u and v with the iteration (D.199) though (D.202) have been studied by many people, for example, Ehrlich [4|, [5J, [6], Greenspan and Schultz [8], McLaurin [13], Smith [20], [21] and Vajtersic [27], etc. Greenspan and Schultz [8] applied the coupled iterative approach to the solution of Airy stress function in a rectangle with an edge crack(Mode I case). They used nine point scheme in the vicinity of the crack tip and solved the coupled equations using the iterations of (D.199) through (D.202) with SOU method in each iteration. For the efficiency and complexities of various methods for solving the above iterative equations, the reader is referred to the cited references. A p p e n d i x E T h e Leas t S q u a r e s So lu t ion of L inea r S y s t e m s Except for few cases in which closed form of analytic solution are rarely available, uniform satisfaction of boundary conditions in most problems is impossible, and some approximation and compromise must be made. For general discussion, write the stress function for our crack problems of the form N « = E**X»M) (E.203) n where Xn a r e basis functions (fundamental solutions) and coefficients qn are to be determined so that on the boundary, with / and g given, boundary conditions N Fu(s) ^ £ g „ f ' x » • f(s), n N Gu(s) - ^ g „ G x „ - g(s) n are satisfied. Here F and G are some differential operators. We choose reasonably distributed points and collocate the boundary condition to obtain N $ > » # * ) = / ( * ) , (E.204) n N J]g„V'W = 5K) (E-205) n for s, on the boundary. Here "reasonable distributed" means the collocation points are chosen so that the significant changes in the boundary values are well captured. The collocation approach mentioned above leads to an M X N linear system Ax = b, where £ is a column vector of unknown coefficients, and alj — $j(s ,) or atJ — ipj(s{) for st on the boundary. Define, with || • || denoting the L2 vector norm, E{x)=f±i\\b-Ax\\j , (E.206) 105 Appendix E. The Least Squares Solution of Linear Systems 106 where M is the number of components of vector b or the number of equations from the collo-cation procedure, then minimizing E{x) is equivalent to solving the following normal equation ATAx = ATb. If we split the linear system so that N £ 9 n # * ) = / (* ) , i=l,...,M/2, N n Examining the coefficient of AT A, we find M / 2 + 1 . . . . . M . M / 2 Af/2 (E.207) (E.208) (E.209) (E.2I0) If the collocation points are condensely distributed, then for s, 6 [a,/?], for example, 0 6 [0,7r], /3 - a T f13 fP M */tt ./at It can be seen that if functions <j> and V a r e orthogonal, then A A will be approximately diagonal dominant, and well conditioned. Unfortunately, the stress functions that we used do not have this property in general. Considering the effect of rounding errors, the orthogonal transformation based method such as QR method is generally recommended rather than the normal equation approach for finding the least squares solution(See [7], [12] and [24]). In the QR method, the M X N coefficient matrix A is decomposed into A^QR (E.212) where Q is of order M X M and orthogonal, R is a upper triangle of order M x N. Then the linear system Ax — b can be written Rx - QTb (E.213) Decompose R and the right hand side into r Ri 0 R (E.214) Appendix E. The Least Squares Solution of Linear Systems 107 Method normal equation QR based method Kj -1.38798 1.86502 Ikll 0.93459 0.0299331 lleIU 3.35803 0.330086 Table E.18: Comparison of the least squares solution of dimensionless Ki and errors in the Mode I compact loading case for c/a = 0.5, h/a —• 0.8, l/h — 0.25 using totally 100 terms and 159 collocation points. (E.215) where Rx is the upper triangular of order TV x N, C\ E RN, c2 fc RM N • Then ||6 - Ax\\2 = \\QTb - QTAx\\2 (E.216) = \\QTb-Rx\\2 (E.217) = ( H d - ^ x l ^ + Hcall2)2. (E.218) So the minimizer is given by x^R^Cl. (E.219) The minimum error is therefore ||c2||. In our numerical experiments, we have found a few cases in which the QR factorization based method is superior to the method using normal equation in solving the least squares problems arising from the boundary collocation procedure. As an example, we present the comparison of the values of dimensionless A'I and corresponding errors in the boundary conditions in the case of a rectangular plate under concentrated loads of tension with c/a = 0.1, h/a = 0.8 and l/h — 0.25 using 100 truncated terms(See Section 4.3). Table E.18 shows the failure of using normal equation approach. Without more knowledge of practical problems, it will be hard to tell tha t the dimensionless K\ obtained from the normal equation is absolutely wrong in this particular example. QTb = c2 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080857/manifest

Comment

Related Items