Mathematical Details on Singular Integral Equation Method for Solving Crack Problems

This article provides a detail derivation of a singular Fredholm integral equation for the solution of a mixed mode crack problem in a nonhomogeneous medium. The integral equation derived here has already been addressed by F. Delale and F. Erdogan (Delale & Erdogan 1983), one of the most cited and pioneer papers in fracture mechanics that uses singulalr integral equation method (SIEM) to solve crack problems. However, probably due to its limit of paper length, some mathematical details are not provided to bring this powerful method, SIEM, to its full strength. In this paper we fill in the mathematical gaps, and both analytical and numerical parts are addressed in details. Some discussions from the view point of differential equations are given, and new numerical outcomes under different loading functions are provided.


Motivation and Background
Rigorously speaking, most numerical approximation methods (including finite elements, finite difference, and finite volume) in engineering applications can be considered as "hybrid" methods.That is, there must be a beginning stage of theoretical developments upon which the discretization steps are constructed so that accuracy, convergence, and stability can be ensured.However, due to the crack-tip singularity (Kaya & Erdogan, 1987;Muskhelishvili, 1953;Muskhelishvili 1963), most computational methods are not adequate to cope with the crack problems.One of the better methods to solve crack problems is by singulalr integral equation method (SIEM) due to its nature of "capturing the crack-tip singularity" (Muskhelishvili, 1953;Sneddon, 1972;Sneddon & Lowengrub, 1969).Usually the solution outcomes obtained by SIEM are considered to be benchmarks for other numerical outcomes (Erdogan, 1978;Kaya & Erdogan, 1987;Konda & Erdogan, 1994).In order to incorporate the crack-tip singularity it takes a somewhat lengthy mathematical derivation in SIEM (Muskhelishvili, 1953;Muskhelishvili. 1963).Though SIEM has been broadly used in solving crack problems, perhaps because of its lengthy mathematical derivation, the method itself does not seem to be very well documented.In other word, whenever SIEM is used, the theoretical derivations always seem to be very brief and sometimes even skipped in many cases.It is this motivation that we present a detail account of mathematical derivation on SIEM.
The example used here is a mode-I crack problem in nonhomogeneous materials, also well known as functionally graded materials, FGMs (Chan, Paulino, & Fannjiang, 2001;Erdogan, 1978;Erdogan, 1995;Erdogan & Ozturk, 1992;Jin & Noda, 1994;Konda & Erdogan, 1994).This is another SIEM advantage-in addition to its accuracy and the ability to incorporate the crack-tip singularity, it can be very flexible to model crack problems in nonhomogeneous materials.This mode-I crack problem in nonhomogeneous materials has already been solved by Delale and Erdogan (Delale & Erdogan 1983).It is highly recommended that readers should have an access to the paper (Delale & Erdogan 1983) when studying this paper.We actually conducted an REU (Research Experience for Undergraduates) grant funded by National Science Foundation by studying paper (Delale & Erdogan 1983) for the "Solid Mechanic Team" in summer 2017 and 2018.We have filled in every theoretical step with mathematical details.

Steps to Singulalr Integral Equation Method
Basically the theoretical procedures that eventually lead to a type of Fredholm integral equation can be divided into following steps (Chan, Paulino, & Fannjiang, 2001;Erdogan, 1978;Erdogan & Ozturk, 1992;Muskhelishvili, 1953;Sneddon, 1966;Sneddon & Lowengrub, 1969): 1. Linear elasticity and geometry of the crack problem 2. Governing partial differential equation (PDE) 3. Fourier transform and inverse Fourier transform 4. Solution to the ordinary differential equation (ODE) 5. Boundary conditions imposed 6. Fredholm integral equation These six steps of formulation of the crack problem is pretty uniform in solving problems in linear elasticity fracture mechanics (LEFM), and it can be considered as a standard solution technique to the partial differential equations (PDEs) that arise in LEFM.We shall follow above outlines in giving detail formulation of the crack problem.Once the Fredholm integral equation is obtained, the discretization steps start, and it becomes the numerical part of SIEM.It is another important component of SIEM.However, due to the page limit we briefly address the numerical procedures, and we may pursue the numerical part in a more detail manner at another occasion.Some of the numerical results are given in terms of crack opening displacement profile and the mode-I stress intensity factors (SIFs).

Linear Elasticity and Geometry of the Crack Problem
A plane elasticity problem is considered.The nonhomogeneous material is in terms of the Young's modulus E = E(x, y), and the Poisson's ratio ν is assumed to be constant (Delale & Erdogan, 1983;Konda & Erdogan, 1994).The basic framework of analysis, background and notations is given below (Gdoutos, 1993;Sneddon & Elliott, 1946;Sneddon & Lowengrub, 1969;).

Strain-Displacement Relations
, where i ≡ x, j ≡ y; or in terms of matrix,

Condition of Strain Compatibility
As a plane problem is considered, following strain compatibility is obtained.
Sum of the two equations in (1) leads to equation (2), which is the statement of condition of compatibility: 2.1.3(Generalized) Hooke's Law Hooke's law relates the strains and stresses through the following equations: The Young's modulus is taken as a function of material position (x, y), i.e., E = E(x, y) (Delale & Erdogan, 1983;Erdogan, 1978;Erdogan & Ozturk, 1992;Konda & Erdogan, 1994).Figure 1 shows a crack sitting in a nonhomogeneous material It is worth of noting is that when the crack problem is solved, we will apply a "superposition principle" (Bueckner, 1958;Delale & Erdogan, 1983;Gdoutos, 1993), see Figure 2.This shall clarify why the boundary conditions in equation ( 30) are imposed.(Bueckner, 1958)

Airy Stress Function
The introduction of Airy stress function is to make a "vector" problem to be a "scalar" problem.That is, instead of solving a 2 × 2 linear system of second order of PDEs, we solve a fourth order of PDE.
Using classical elasticity , we let the Airy stress function be Φ = Φ(x, y) such that: Relations between Airy stress function Φ = Φ(x, y) and strains are linked by (generalized) Hooke's law as: From the compatibility equation (3), we have1 Term ♠ undergoes the following process: Term ♣ takes very similar route: Finally, term ♡ can be rewritten as Substitution of equations ( 8), ( 9) and ( 10) in the compatibility equation (3) leads to which is the governing PDE of the elasticity problem for generalized plane stress.The PDE for plane strain can be obtained by replacing E and ν by E 1−ν 2 and ν 1−ν , respectively.

Simplification of the PDE
Clearly the fourth order PDE (11) depends on the function form of the Young's modulus E(x, y).As a mathematical convenience we let E = E(x, y) = E 0 e βx+γy (Delale Erdogan, 1983;Erdogan, 1978;Erdogan & Ozturk, 1992) where β and γ are material parameters, then the PDE ( 11) is simplified to By letting γ = 0, the PDE (12) above reduces to the simpler form which is the form that will be investigated in the remainder of this manuscript.

Fourier Transform and Inverse Fourier Transform
One of the most useful methods for solving PDEs is the Fourier transform, which transforms the PDE into an ordinary differential equation () By solving the ODE then we solve the associated PDE by taking its inverse Fourier transform.Let the Fourier transform and the inverse Fourier transform be defined as that is, Φ(x, y) is the inverse Fourier transform of the function f (y, α).
Considering each of the terms in equation ( 13), we have: Summing equations ( 17) to ( 21) leads to the following ODE:

Solution to the ODE
The corresponding characteristic equation to the ODE ( 22) is and function f (y, α) takes the form of e m i (α)y , where m i (i = 1, 2, 3, 4)2 are the roots of equation ( 23).We can assume that the real parts of m 1 and m 2 are positive and for the polynomial in equation ( 23) has only even powers.Without loss of generality we may find3 that We like to point out some observations: • Roots m i 's may take very complicated algebraic form, but sum or product of m i 's can be expressed in simpler form.
Depending on the Fredholm kernel, we may want to find other forms of m i 's (e.g. . ., etc.) rather than each exact m i .
• Roots m i 's are in terms of material parameters.
As we are only interested in the case that the displacements and stresses vanish as (x 2 + y 2 ) → ∞, the so called far-field condition, solution to ODE (22) takes the form4 Before the boundary conditions are imposed in the next sub-section, we need to express each stress component in terms of the solution obtained from ODE. Equations ( 6), ( 16), and (26) lead to the following expression for the stress components.

Boundary Conditions Imposed
The following boundary conditions corresponding to a Mode I crack problem are imposed: where p(x) is a known function (loading).The material function is simplified to5 and v is the y-component of the displacement which needs to be determined.
The determination of A 1 and A 2 facilitates completing the problem on the upper half plane (y > 0).A 1 and A 2 are determined by (any two) boundary conditions on y = 0 , −∞ < x < ∞.Consider that the half (upper) plane is subjected to traction Equations ( 28) and ( 29) give the following equations ( 33) and (34), respectively.
By taking Fourier transform one obtains Equation ( 35) plus equation ( 36) multiplied by iα m 2 gives that Thus A 1 can be obtained as To get A 2 we consider equation ( 35) plus equation ( 36) multiplied by iα m 1 , which gives so that A 2 is obtained as

Fredholm Integral Equation
Recall the boundary conditions in equation ( 30).Condition (a), together with equation ( 29), leads to In order to get both A 1 and A 2 we need another equation.Thus an unknown function g(x) is introduced and defined as: Note that g(x) = 0 for |x| > a and Hence, using Hooke's law together with equations ( 27) and ( 28), we also have Using equations ( 42) and ( 45), and the fact that Separating terms for e −m 1 y and e −m 2 y , and noting that y > 0, we get Now the function g(x) can be evaluated as Multiplying equation ( 48) by e βx , we obtain The above equation says that e βx (g(x)−c ′ (x)) is the Inverse Fourier transform of function h(α).By taking Fourier transform and using the fact that g(x) = 0 for |x| > a , we have: From equation ( 41), A 2 = −m 1 A 1 m 2 ; using this expression in the equation above, one obtains By simplifying the terms inside the bracket, we get which leads to By substituting −a e (β+iα)t g(t)dt into σ yy in equation ( 28) and using boundary condition (b) σ yy (x, +0) = p(x) in equation ( 30), we obtain which is the expression for the loading on the crack face.The derived integral equation ( 54) is the target that we have been seeking, and asymptotic analysis needs to be conducted here to determine the type of singularity that the kernel K(y, α) may possess (Mikhlin, 1964;Muskhelishvili, 1953;Sneddon & Lowengrub, 1969).
2.6.1 Asymptotic Analysis of Kernel K(y, α) in Equation ( 54) Equations ( 24) and ( 25) show that the asymptotic behavior of m 1 and m 2 is Let K(y, α) be the following part of the integrand in the inner integral of equation ( 54): Accordingly, let the inner integral of equation ( 54) be denoted as and let K ∞ be the asymptotic value of K(y, α) , that is Substituting equation (55) into equation ( 56), which states that m 1 → |α| and m 2 → |α| as |α| → ∞, we obtain Using addition and subtraction, equation ( 57) can be rewritten as The second integral in the above equation may be expressed as where the last equality is obtained from integrating by parts twice.It is worth of pointing out that equation ( 60) is vital in SIEM, for this is the step that the singularity is isolated out, and this "singular integral" is carried out analytically. Let Substituting equations ( 56), (59), and ( 61) into (60), we obtain By letting y → +0 in the first integral of equation ( 63) 7 , one obtains The integral in equation ( 64) can be split into two integrals: So we obtain h(x, y, t) as Considering the limit y → +0 and substituting equation ( 65) into (54), we obtain Replacing E 0 2 by 4µ 0 1+κ , we get The equation can be rewritten in the form where the Fredholm kernel is defined by Therefore, we have reached the singular integral equation ( 68).It has a Cauchy singularity, and the integral is evaluated according to Cauchy principal value.Due to the complicated form of the regular kernel k(x, t) in equation ( 69), one has to undergo a numerical approximation, which will be addressed in the next section.

Numerical Part of SIEM
The numerical part of SIEM actually takes about the same amount effort as the theoretical part.Due to the focus of this work is on the derivation of the Fredholm singular integral equation, we briefly address the numerical part of SIEM here.
Our plan is to deliver another work on the details of the numerical part of SIEM, same as this paper about the theoretical part, in summer 2018 REU program.

Numerical Procedures of SIEM
The numerical solution of equation ( 68) can be divided in the following six steps (same step number as the theoretical part): • Normalization • Representation of the Density Function

• Tchebyshev Polynomial Expansion
• Formation of Linear System of Equations

• Evaluation of Singular and Hypersingular Integrals
• Evaluation of Non-singular Integral We shall briefly address some of six steps below.

Normalization
For numerical solution we normalize the interval (−a, a) by defining so that the normalized version of integral equation ( 68) becomes accompanied with the single-valued condition ∫ 1 −1 ϕ(s)ds = 0 , where a is the crack length, β , κ , and µ 0 are material parameters, ϕ is the slope of the crack profile, n is the Fredholm kernel, q is the loading function, s is the integration variable, and r is the collocation variable.

Representation of the Density Function
Because of the Cauchy singularity, we select our density function to be This step is actually the key step that how SIEM can capture the crack singularity.The representation of the density function by a function 1 √ 1−s 2 is a mathematical outcome from solving a Cauchy singular integral equation (Hadamard, 1952;Martin, 1991;Martin, 1992;Muskhelishvili, 1953;Muskhelishvili, 1963;Sölingen, 1939).

Tchebyshev Polynomial Expansion
The Φ(s) in equation ( 71) is approximated by Tchebyshev polynomials (Hochstrasser, 1972): where T n (s) are the Tchebyshev polynomials of first kind, and A n 's are coefficients that need to be determined numerically.
The upper index N − 1 is related to the number of collocation points we choose.

Formation of Linear System of Equations
By the approximation of Φ(s) in equation ( 72), if we choose N collocation points, then we will need to solve an N × N linear system of the unknown A n 's.

Evaluation of Singular Integrals
The main formula used in evaluation of the singular integrals is () where the T n (r) and U n (r) are the Tchebyshev polynomials of the first and second type, respectively.

Evaluation of Non-singular Integrals
Some quadrature methods are used in evaluation of the non-singular integrals.As the integral domain is from 0 to infinity, convergence test at this step needs to be carefully conducted to make sure the numerical results are stable.

Stress Intensity Factors
Since the propagation of a crack starts around its tips, it is very important to study and determine the stress intensity factors (SIFs) at both crack tips.The mode-I SIFs can be calculated from (Erdogan, 1995;Gdoutos, 1993) K I (a) = lim To compute the SIFs, a useful integral involving Tchebyshev polynomials is given here [?]: Notice that the integral is no longer singular, as |r| > 1.We report the SIFs in Figure 3.

Crack Surface Displacement
We also compute the crack surface displacements at various material parameters βa, and they are shown in Figurer.The numerical outcomes have been compared with the results by Delale and Erdogan (Delale & Erdogan, 1983) and they are consistent.

Conclusion Remarks
Our contribution of this paper falls into the review category, but singulalr integral equation method (SIEM) is a method that a lot of original research can take advantage of it; just like finite element method (FEM) has been so widely used in fracture mechanics.However, compared with other numerical methods, SIEM sometimes is called a "semi-analytical" method due to its nature involving heavy analytical work.As SIEM is so suitable for dealing with crack singularity and its flexibility in the underlying linear elasticity theory, it has been used so frequently in fracture mechanics, and authors believe that SIEM deserves a comprehensive and a well documented treatment.Our work in this paper is focus on theoretical part, and we also briefly address the numerical part.(Authors plan to have another focus work on the numerical part in the near future.)We also provide some solution results including crack surface displacement profiles and stress intensity factors (SIFs).By a comprehensive and well organized presentation of SIEM, we hope that fracture mechanics community can benefit from our work.

Figure 4 .
Figure 4. Crack surface displacements in an infinite nonhomogeneous plane under uniform crack surface loading σ yy (x, 0) = −p 0 at various material parameters βa, and a is the half crack length; ν = 0.3.