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. The -integral and the stress intensity factor (SIF) are the most important parameters for crack problems. In this study -integral has been used to calculate the SIF in the crack tip. The mode 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 -integral value, is conducted. Visibility and diffraction methods increase the accuracy of RKPM results and effect on the -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 -integral and SIF value.
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 -integral. The -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 -integral is that it is path independent and it helps to calculate the -integral at a far distance from the crack tip. In linear elastic fracture mechanics the -integral has a direct relationship with the stress intensity factors (SIFs). In this study the -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 of can be written as where is the modified kernel function on domain that is expressed by where is window function, is a correction function, and 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 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: Substituting (4) into (1) leads to In order to simplify (5), the th degree moment matrix of functionis defined by Then (6) will be rewritten in the form of In order to exactly reproduce the th order polynomial function, the following conditions need to be satisfied: Or in summary: If a correction function including unknown coefficient is defined, equations of (9) can be satisfied simultaneously. The correction function is defined by It can be also expressed in matrix form: where is a set of basic functions and including components and is a set of unknown coefficient. Substituting (11) into (9) and considering definition of moment matrix in (6) lead to From (12) the unknown coefficient sets of are obtained. Equation (12) can also be rewritten as Or it can be shown in matrix form as (14) and (15): Moment matrix can be shown as Since the window function is always positive, all the components of moment matrix are linearly independent with respect to . Therefore, the moment matrix is nonsingular. Hence, simultaneously solving (15), the unknown coefficient sets of are obtained: After obtaining the unknown coefficient sets the correction function can be easily calculated from (10). After obtaining the unknown coefficient sets of , the correction function is determined and the function 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: where NP is the total number of particles distributed throughout the domain . Equation (18) can be rewritten as and is called shape function: where is the particle number on the domain, is the coordinate of that particle, is the kernel function, and is defined as the shape function of particle 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 axis. In this study, we employ the cubic spline function as the kernel function, which is By considering (19) and deriving from (20) with respect to , the definition of the derivative of the th shape function becomes
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 By multiplying the above kernel function in each dimension, the multidimensional kernel function is obtained: where nsd is the number of dimensional Euclidean space and is dilation parameter for th dimension. The multi-dimension correction function is defined based on (11), where is based multidimensional vector function which is included in which is obtained by applying the condition of completeness: where is So the correction function will be And the modified kernel function will be And finally, the multidimensional reproducing function form will become And it will be discretized to determine the shape function: where is the corresponding volume of the th 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 which is 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 (Figure 1).
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
A 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.
(a)
(b)
(a)
(b)
(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 and close to the crack tip using visibility criterion can be shown as Figure 5. As can be seen, the particles that at particles or 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.
Figures 6(a) and 6(b) show the 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 for particle , then a circle with the center being the crack tip and radius of 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).
Figures 8(a) and 8(b) show the window function of a node next to the crack using visibility criterion [7].
(a)
(b)
7. SIF and -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 and -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 -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 distance and angle of with the -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 where parameter is the SIF for different modes in the crack tip, and shown , , and 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. , , and 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 -integral was defined as [10] where is strain energy density, is stress tensor, is the normal to the curve , and is the displacement vector. The strain energy density is given by Also, -integral can be obtained in terms of SIF of the first, second, and third modes: where is shear modulus, is modulus of elasticity, and is Poisson ratio. An important feature of the -integral is that it is path independent and this helps to calculate the -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 -axis, then -integral can be shown in matrix form as and and are the stresses in an arbitrary direction which has angle of with -axis: and are displacement in the same direction: Substituting (38) and (39) in (37) the -integral will be easily calculated. Figure 10 shows that is considered as a path for the -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.
It is obvious that stress is in elastic condition and it can be stated in matrix form: Then strain energy density is calculated from Substituting (40) in (41), strain energy density will be And the -integral on the closed path is Then SIF is calculated from (44) for plane-stress and plane-strain conditions
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 in and 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 106, in which is Young’s modulus. A rectangular steel plate is selected with size of ft2 ( 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.
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 Figure 12.
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 , (distance between particles in and 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. Figures 13 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.
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 dimensions ft2 including 825 particles, -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 -integral values are calculated. Figures 15 and 16 show the -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 -integral value and far-field tensile stress. Diffraction criterion seems to have better results for the -integral versus tensile stress. Also, the result of -integral value in diffraction method is almost 10 percent more than visibility method as shown in Figure 17.
9. Conclusions
(1)When the dilation parameter increases for all particles in plastic analysis, the difference of the -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 -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)The -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 between -integral and far-field tensile stress for mode edge crack and plane strain condition.