# Comparison between Visibility and Diffraction Criteria on SIF andJ-Integral Value for Mode I Crack Using RKPM

Authors: Masood Hajali; Caesar Abishdid; Arash Nejadpak
Journal: Journal of Structures
 (2013)

---

## Abstract

Recently, mesh-free methods are increasingly utilized in solving various types of boundary value problems. Much research has been done on mesh-free methods for solving differential equation problems including crack and also obtained satisfactory results. Among these methods, reproducing kernel particle method (RKPM) has been used increasingly in fracture mechanic problems. TheJ-integral and the stress intensity factor (SIF) are the most important parameters for crack problems. In this study J-integral has been used to calculate the SIF in the crack tip. The mode I SIF at the crack tip in a work-hardening material is obtained for various dilation parameters using RKPM. A comparison between two conventional treatments, visibility and diffraction on SIF and J-integral value, is conducted. Visibility and diffraction methods increase the accuracy of RKPM results and effect on the J-integral results at the crack tip. In comparing between the visibility and diffraction methods to modify the shape functions, the diffraction criterion seems to have better results for the J-integral and SIF value.


---

## Body

## 1. Introduction

In fracture mechanic problems, the finite element formulations have been well developed and several amounts of research has been accomplished. Standard finite element approaches for crack problems are usually ineffective due to the mesh-based view and propagation of the crack during the loading process. Mesh-free methods eliminate some or all of the traditional mesh-based views of the computational domain and rely on a particle view of the field problem. The major difference between finite element methods is that the domain of interest is discretized only with nodes, often called particles. In recent years, much research has been done on mesh-free methods for solving differential equation problems including crack and also obtained satisfactory results. Among these methods reproducing kernel particle method (RKPM) has been used increasingly in fracture mechanic problems. Boundary value problems (BVPs) often have essential boundary conditions (EBCs) that involve derivatives, for example, in beams and plates, where slopes are commonly enforced at the boundaries. Such problems are solved numerically using mesh-free techniques like the RKPM and the EFGM.In fracture mechanic problems, the concept of energy release rate was first introduced by Cherepanov [1] and Eshelby [2], but it was Rice who first used this independent path integral in fracture mechanics problems. In 1968, Rice [3] presented the concept of energy release rate by means of J-integral.  The J-integral represents a way to calculate the strain energy release rate, or work (energy) per unit fracture surface area, in a material. An important feature of the J-integral is that it is path independent and it helps to calculate the J-integral at a far distance from the crack tip. In linear elastic fracture mechanics the J-integral has a direct relationship with the stress intensity factors (SIFs). In this study the J-integral has been used to calculate the SIF at the crack tip.There have been two widely used criteria, visibility and diffraction which are used for dealing with the internal discontinuity in fracture mechanics. A visibility criterion has been developed by Belytschko et al. in 1994 [4] to modify the shape function of the particles near to discontinuities or nonconvex parts such as cracks. Visibility is the easiest way to introduce discontinuities in mesh-free methods. In this method, the boundaries of the domain and any internal discontinuities can be considered as an opaque barrier. Also, discontinuity means considering that the influence domain of a shape function for a particle and the connection line between this particle and other particles can be assumed as a light beam in the influence domain of the shape function. A diffraction criterion has been developed by Belytschko et al. in 1996 [5] to modify the shape function of the particles in the vicinity of the discontinuities such as cracks.


## 2. RKPM Shape Functions and Their Derivatives

The reproduced kernel function ofu(x) can be written as
(1)uR(ξ)=∫Ωu(x)ϕ¯(ξ;ξ-x)dx,
where ϕ-(ξ;ξ-x)is the modified kernel function on domain Ω that is expressed by
(2)ϕ¯(ξ;ξ-x)=C(ξ;ξ-x)ϕ(ξ-x),(3)ϕa(ξ-xi)=1aϕ(ξ-xia),
where ϕa(ξ-x) is window function, C(ξ;ξ-x) is a correction function, and a is the dilation parameter of the kernel function. Dilation parameter is defined in order to make more flexibility for the window function and this parameter will control the expansion of the window function on the domain. The correction function C(ξ;ξ-x) proposed by Liu et al. is shown by a linear combination of polynomial including some unknown coefficients. These unknown coefficients will be computed after imposing the boundary conditions. In order to get the equations for reproducing an arbitrary function, consider the following Taylor series expansion:
(4)u(x)=∑α=0∞(-1)αα!(ξ-x)αu(α)(ξ).
Substituting (4) into (1) leads to
(5)uR(ξ)=∑α=0∞(-1)αα!(∫Ω(ξ-x)αϕ¯(ξ;ξ-x)dx)u(α)(ξ).
In order to simplify (5), the αth degree moment matrix of functionϕ-a(ξ,ξ-x)is defined by
(6)m-α(ξ)=∫Ω(ξ-x)αϕ-a(ξ;ξ-x)dx.
Then (6) will be rewritten in the form of
(7)uR(ξ)=m0¯(ξ)u(ξ)+∑α=1∞(-1)αα!mα¯(ξ)u(α)(ξ).
In order to exactly reproduce the nth order polynomial function, the following conditions need to be satisfied:
(8)m-o(ξ)=1,m-α(ξ)=0,α=1,2,…,n.
Or in summary:
(9)m-α(ξ)=δαo;α=0,1,2,…,n.
If a correction function including n+1unknown coefficient is defined, n+1 equations of (9) can be satisfied simultaneously. The correction function is defined by
(10)C(ξ,ξ-x)=∑α=0nβα(ξ)(ξ-x)α.
It can be also expressed in matrix form:
(11)C(ξ;ξ-x)=PT(ξ-x)β(ξ),
where PT(ξ-x) is a set of basic functions and including n+1 components and β(ξ) is a set of unknown coefficient. Substituting (11) into (9) and considering definition of moment matrix in (6) lead to
(12)∫Ω〈(ξ-x)α(ξ-x)α+1⋯(ξ-x)α+n〉ϕa(ξ-x)dx×{βo(ξ)β1(ξ)⋮βn(ξ)}=δαo,α=1,2,…,n.
From (12) the unknown coefficient sets of βi(ξ) are obtained. Equation (12) can also be rewritten as
(13)〈mα(ξ)mα+1(ξ)⋯mα+n〉{βo(ξ)β1(ξ)⋮βn(ξ)}=δαo.
Or it can be shown in matrix form as (14) and (15):
(14)[mo(ξ)m1(ξ)⋯mn(ξ)m1(ξ)m2(ξ)⋯mn+1(ξ)⋮⋮⋯⋮mn(ξ)mn+1(ξ)⋯m2n(ξ)]{βo(ξ)β1(ξ)⋮βn(ξ)}=[10⋮0],(15)M(ξ)β(ξ)=P(0).
Moment matrix M can be shown as
(16)M(ξ)=∫ΩP(ξ-x)PT(ξ-x)ϕa(ξ-x)dx.
Since the window function is always positive, all the components of moment matrix are linearly independent with respect to ϕa. Therefore, the moment matrix is nonsingular. Hence, simultaneously solving (15), the unknown coefficient sets of βi(ξ) are obtained:
(17)β(ξ)=M-1(ξ)P(0).
After obtaining the unknown coefficient sets βi(ξ) the correction function can be easily calculated from (10). After obtaining the unknown coefficient sets of βi(ξ), the correction function is determined and the function u(x) or its derivatives can be obtained using the reproducing function. Equation (17) can be discretized in order to apply to various problems. Equation (18) is the result of the discretization in the reproducing equation using the trapezoid integration method:
(18)uR(ξ)=∫Ωu(x)ϕ¯(ξ;ξ-x)dx=∑i=1NPu(x)ϕ¯(ξ;ξ-x)Δx,
where  NP is the total number of particles distributed throughout the domain  Ω. Equation (18) can be rewritten as
(19)uR(ξ)=∑i=1NPψi(ξ)ui,
and ψi(ξ) is called shape function:
(20)ψi(ξ)=PT(ξ-xi)βi(ξ)ϕa(ξ-xi)Δxi,
where i is the particle number on the domain, xi is the coordinate of that particle, ϕa(ξ-x) is the kernel function, and ψi(ξ) is defined as the shape function of particle i with coordinate of ξ. Two most regular kernel functions which are used in mesh-free methods are Gaussian and spline functions. All the spline functions are symmetric around x=0 axis. In this study, we employ the cubic spline function as the kernel function, which is
(21)S3(x)={23-4x2-4x30≤x<1243-4x+4x2-43x312≤x<10x≥1.
By considering (19) and deriving from (20) with respect to xk, the definition of the derivative of the Ith shape function becomes
(22)ψi,xk(ξ)={P,xkT(ξ-xi)βi(ξ)ϕa(ξ-xi)+PT(ξ-xi)βi,xk(ξ)ϕa(ξ-xi)+PT(ξ-xi)βi(ξ)ϕa,xk(ξ-xi)}Δxi.


## 3. Calculation of Multidimensional Shape Functions in RKPM

In what proceeded, the relationships needed to reproduce the desired function and its derivatives in a continuous and discrete computational domain. Now, the shape functions will be calculated for the multidimensional problems. The kernel function in one-dimension was(23)ϕa(ξ-xI)=1aϕ(ξ-xIa).
By multiplying the above kernel function in each dimension, the multidimensional kernel function is obtained:
(24)ϕa(ξ-xI)=∏j=1nsd1ajϕ(ξ-xI/jaj),
where nsd is the number of dimensional Euclidean space and ajis dilation parameter for jth dimension. The multi-dimension correction function is defined based on (11), where P is based multidimensional vector function which is included in
(25)PT(ξ-xI)=[1ξ1-xI/1ξ2-xI/2⋯(ξnsd-xI/nsd)N]
which is obtained by applying the condition of completeness:
(26)β(ξ)=M-1(ξ)P(0),
where M is
(27)M(ξ)=∫ΩP(ξ-x)PT(ξ-x)dx1⋯dxnsd.
So the correction function will be
(28)C(ξ;ξ-xI)=PT(0)M-1(ξ)PT(ξ-xI).
And the modified kernel function will be
(29)ϕ-a(ξ;ξ-xI)=C(ξ;ξ-xI)ϕa(ξ-xI).
And finally, the multidimensional reproducing function form will become
(30)uR(ξ)=∫Ωϕ-a(ξ;ξ-x)u(x)dx1⋯dxnsd.
And it will be discretized to determine the shape function:
(31)ψi(ξ)=ϕ-a(ξ;ξI-xI3)ΔVI,ψi(ξ)=PT(0)M-1(ξ)P(ξ-xI)ϕa(ξ-xI)ΔVI,
where  ΔVI  is the corresponding volume of the Ith particle and, in two-dimensional problems, is each particle’s ration area.


## 4. Ration Area of Each Particle in Two-Dimensional RKPM

Ration area in two-dimensional problems for each particle is shown byΔSI which is
(32)ΔVI=ΔSI.
The ration area is considered one in element-free Galerkin method (EFG). EFG method was introduced by Belytschko and Tabbara in 1992 [6] but it is not accurate because for all particles the ration area is one. But in RKPM the ration area of each particle is calculated separately and is different for each particle. Initially, all the particles which are scattered in the domain are connected to each other. Then for each particle perpendicular lines between the particles with other particles are drawn, and the area of the polygons trapped for each particle is considered as ΔSI (Figure 1).Figure 1
Determination of ration area for each node.In the FORTRAN programming that has been developed in this study a different method has been used to calculate the ration area of each particle. First the domain has been divided into small raster rectangular (1000*2000). The distance between coordinates of the center of each rectangular is compared with all particles. The rectangle which has the smallest distance from each particle is allocated to the ration area of that particle.


## 5. Display of Two-Dimensional Shape Functions and Their Derivatives

A5×5 plate including 2500 particles is selected. The particles are uniformly scattered on the plate. The dilation parameter is considered to be equal to one. Figures 2, 3, and 4 show the two-dimensional shape functions and their derivatives for different particles.Shape functions of a vertex particle.
(a)
(b)Shape functions of a central and a middle particle.
(a)
(b)Derivatives of the shape function for central particle along theY and X direction.
(a)
(b)


## 6. Modification of RKPM Shape Functions Using Visibility and Diffraction Criteria

Through engineering problems, the domain of the problem may contain nonconvex boundaries, particularly the fracture ones having discontinuous displacement fields.  In such conditions, the shape functions associated with particles, whose supports intersect the discontinuity, should be modified. One of these criteria is the visibility introduced by Belytschko et al. [4] and  Krysl and Belytschko [9]. In this approach, if the assumed light beam meets the discontinuity line, the shape function after the barrier will be cut. Therefore, a discontinuity is applied to the geometry. For example, if a crack is considered, the influence domain of particles I and J close to the crack tip using visibility criterion can be shown as Figure 5. As can be seen, the particles that at particles I or J cannot be seen by an observer will be removed. In other words, the shape function of the particles which the crack or discontinuity prevent from reaching the light beam will be modified to amount to a zero.Figure 5
Modified influence domain of the particle next to the line of discontinuity using visibility criterion.Figures6(a) and 6(b) show the window function of a node next to the crack using visibility criterion [7].Modified window function of a node next to the crack using visibility criterion [7].
(a)
(b)Diffraction criterion (Belytschko et al. in 1996 [5]) is based on the bending of the light beam which has been described in the visibility criterion around a tip discontinuity. Consider the end of the discontinuity line in Figure 7. If the distance between the crack tip and the end of the arc is called d for particle I, then a circle with the center being the crack tip and radius of d is drawn. Areas outside the circle and behind the discontinuity are removed and the amount of the shape function in these areas will be zero (Figure 7).Figure 7
Modified influence domain of the particles next to the line of discontinuity using diffraction.Figures8(a) and 8(b) show the window function of a node next to the crack using visibility criterion [7].Modified window function of a node next to the crack using diffraction criterion [7].
(a)
(b)


## 7. SIF andJ-Integral

The main purpose of fracture mechanics is to determine the status of cracks in different loading conditions. Stress, strain, displacement, and energy fields are required to obtain a driving force for crack growth. SIF andJ-integral are two important concepts of crack problems. SIF is used to quantify the stress field around the crack tip. Many methods have been developed to determine the stress intensity factor. One of these methods to calculate the stress intensity factor is J-integral. The material can be cracked in three different modes including opening, shearing, and tearing mode.The first mode is related to the tensile stress which is orthogonal to the page and is called opening mode. The second mode is applied in the same way as shear stress on the page and is called shearing or sliding mode. The third mode which includes out-of-plane shear stress of the crack is called tearing mode. Other situations of the loading are combinations of these three modes. If a node is considered with distancer and angle of α with the x-axis in the vicinity of the crack edge (see Figure 9), then the stress field in this node is calculated according to the Irwin method in different crack modes. Therefore, stress field in the crack tip for linear elastic materials is calculated by
(33)6ij=K2πrfij(θ),
where K parameter is the SIF for different modes in the crack tip, and shown KI, KII, and KIII are for the first, second, and third modes. Values of these coefficients are determined according to the dimensions and loading condition of the problem. Therefore, the SIF relationship is calculated from the analysis of the geometrical and loading condition. KI, KII, and KIII are physically the intensity of force transfer at the crack tip due to the creation of the crack in the material. SIF plays an important role as a failure parameter. Rice [3] also showed that this integral has linear elastic attitude with the energy release rate and was independent of the path around a crack. The two-dimensional J-integral was defined as [10]
(34)J=∫Γ(Wdx2-σijnjui,1ds),j=1,2,
where W is strain energy density, σ is stress tensor, n is the normal to the curve Γ, and u is the displacement vector. The strain energy density is given by
(35)W=∫0εijσijdεij.
Also, J-integral can be obtained in terms of SIF of the first, second, and third modes:
(36)J=1E′(KI2-KII2)+12μKIII2,E′={E1-ν2planestrainEplanestress,
where μ is shear modulus, E′ is modulus of elasticity, and ν is Poisson ratio. An important feature of the J-integral is that it is path independent and this helps to calculate the J-integral in a far distance from the crack tip. If Γ is considered as path independent around an inclined crack tip which has angle of α with the x-axis, then J-integral can be shown in matrix form as
(37)J=∫ΓWdη-∫Γ{σnτn}[∂un∂ξ∂un∂ξ]ds,
and σn and τn are the stresses in an arbitrary direction which has angle of αwith x-axis:
(38)σn=σxxcos2α+σyysin2α+τxysinαcosα,τn=(σyy-σxx)sinαcosα+τxy(cos2α-sin2α)un and νn are displacement in the same direction:
(39)un=ucosα+vsinα,Vn=-usinα+vsinα.
Substituting (38) and (39) in (37) the J-integral will be easily calculated. Figure 10 shows that Q1Q2Q3Q4 is considered as a path for the J-integral in a fully elastic domain. First, the shape of the integral path is described and then the value of integral is calculated on each separate path for two plane stress and plane strain conditions.Figure 9
Crack tip coordinate [8].Figure 10
J-integral path at the crack tip.It is obvious that stress isσ=D·εin elastic condition and it can be stated in matrix form:
(40){σxxσyyσxy}=[D11D120D12D22000D11].{εxxεyyεxy}=[D11εxx+D12εyyD12εxx+D22εyyD33εxy].
Then strain energy density is calculated from
(41)W=12σTε.
Substituting (40) in (41), strain energy density will be
(42)W=12εxx(D11εxx+D12εyy)+13εyy(D12εxx+D22εyy)+12D33εxy2,W=12D11εxy2+12D22εyy2+D12εxxεyy+12D33εxy2.
And the J-integral on the closed path is
(43)J=∫Q1Q2+∫Q2Q3+∫Q3Q4+∫Q4Q1.
Then SIF is calculated from (44) for plane-stress and plane-strain conditions
(44)J=K2E⟹K=JEplane-stressJ=K2E(1-ν2)⟹K=JE1-ν2plane-strain.


## 8. Edge Crack Modeling in RKPM

With what was stated previously and using a FORTRAN program that was written for solving the liner elastic on a steel plate with specified dimension using RKPM, the stress, strain, and displacement field inx and y directions in all computational particles and calculation of SIF under plane-stress and plane-strain conditions were obtained. Penalty method is used to apply the boundary conditions. Penalty coefficient, β, is adopted as 106E, in which E is Young’s modulus. A rectangular steel plate is selected with size of 3×10ft2 (0.91×3.05 m2). An edge crack is considered with a length of 7 inch (0.178 m) in the middle of the plate. A tensile stress of 35 ksi (241.3 MPa) is applied at the bottom and the top of the plate. The loading increment is assumed 7 ksi (48.3 MPa). Roller constraint is used for the plane in front of the crack and pin constraint is used for the front face of the plate (Figure 11). Spline 3rd degree is used as a window function.Figure 11
Domain and boundary conditions.The modulus of elasticity of the plate is 30,000 ksi (207,000 MPa), Poisson ratio of 0.3, and hardening parameter is 10. The problem geometry is included 752 particles uniformly scattered on the surface of the plate and 73 particles positioned on the circles with angles of 30 degrees around the crack tip as shown in Figure12.Figure 12
Particle arrangements.Dilation parameter has a great impact on the performance of the RKPM and especially on the shape function. The main duty of this parameter is determining the extent of the window function on the computational domain. Also, the window functions with smaller influence domain will have better performance in a continuous mode. Window functions with a very small range of influence domain will only include the same particle that has been defined on it, in which case the answers will not be sustainable. The dilation parameter needs to be defined such that it covers two particles of the domain. This is due to the fact that the window functions which cover just one particle will lead to numerical instability in the reproducing relations in the discrete case. Generally, it is hard to find any relation between a dilation parameter andΔx, Δy (distance between particles in x and y directions), and the optimized dilation parameter might be different for different problems. Here, we compare the results of SIF for various dilation parameters between RKPM and analytical solution.To analyze a plate with elastic-plastic behavior under static load, a force has to be applied in consecutive phases. At each phase, the amount of strain is obtained, and, with having all the parameters required in the previous stage, the status of the plate is determined in the new phase. By status of the plate, it is meant the stress in the plate and the amount of plastic strain. For this plate, dilation parameters are compared for two visibility and diffraction criteria. The Gaussian number is considered to be equal to 3. Figures13 and 14 show SIF results versus dilation parameter for both RKPM and analytical methods. It can be concluded that diffraction criterion has the better results and closer results to the analytical solution.Figure 13
SIF versus dilation parameter using visibility criterion.Figure 14
SIF versus dilation parameter using diffraction criterion.The FORTRAN program developed for the elastic-plastic material is also able to recognize the yielded particles according to Von-Mises criterion. For the same plate with dimensions3×10 ft2 including 825 particles, J-integral diagram is shown for plane-strain conditions. The crack tip region was refined using more particles in the circle arrangements. In this problem, dilation parameter is 0.1 for the refined particles and is 0.13 for the rest of particles. Tensile stress of 35 ksi (241.3 MPa) is applied at the bottom and the top of the plate. In each 7 ksi (48.3 MPa) loading increment, the J-integral values are calculated. Figures 15 and 16 show the J-integral values versus tensile stress for plane-strain condition. From Figures 15 and 16, it can be concluded that there is an almost parabolic relationship between J-integral value and far-field tensile stress. Diffraction criterion seems to have better results for the J-integral versus tensile stress. Also, the result of J-integral value in diffraction method is almost 10 percent more than visibility method as shown in Figure 17.Figure 15
J-integral values versus tensile stress using visibility criterion.Figure 16
J-integral values versus tensile stress using diffraction criterion.Figure 17
Comparison ofJ-integral values versus tensile stress for diffraction and visibility criterion.


## 9. Conclusions

(1)
When the dilation parameter increases for all particles in plastic analysis, the difference of theJ-integral in fully plastic regions and fully elastic regions will increase. The reason for this is that when the crack is analyzed in the elastic-plastic condition, for the particles in crack tip with less dilation parameter J-integral is calculated in fully plastic and domain of influence domain does not enter to elastic region. Also, the SIF versus dilation parameter graphs show that increasing the density of particles at the crack tip using a circle or a star arrangement will result in more realistic answers for SIF.
(2)
TheJ-integral, which is the energy discharge rate for a fully elastic analysis, is less than the elastic-plastic analysis, using the Ramberg-Osgood model. The reason is that, in the elastic-plastic condition, the material at the crack tip experiences more strain than in the fully elastic condition and more energy discharge per very small crack growth.
(3)
In comparing between the visibility and diffraction methods to modify the shape functions, the diffraction criterion seems to have better results for the SIF in both the elastic and plastic analysis.
(4)
There is an almost parabolic relationship betweenJ-integral and far-field tensile stress for mode I edge crack and plane strain condition.



---

*Source: 978684-2013-08-28.xml*